Bootstrap w estymacji punktowej i przedziałowej na przykładzie badania średniej.

Zapoznaj się ze składnią i działaniem następujących funkcji w pakiecie R: sample(), rnorm(), replicate(), boot (pakiet boot), boot.ci (pakiet boot).

Zadania

Zadanie 1. Estymacja błędu standardowego średniej próbkowej metodą bootstrap.

library(foreign)
dane <- read.spss("~/Downloads/dane.sav", to.data.frame = TRUE)
attach(dane)
dochod<-Dochod_na_osobe

#punktowo klasycznie
mean(dochod)
## [1] 808.5322
sd(dochod)
## [1] 1064.67
se<-sd(dochod)/sqrt(length(dochod))
se
## [1] 6.146569
sep<-se/mean(dochod)
sep
## [1] 0.007602133
#przedziałowo klasycznie
mean(dochod)-1.96*se   #dolna krawędź
## [1] 796.4849
mean(dochod)+1.96*se   #górna krawędź
## [1] 820.5794
# 1.96 to kwantyl rozkładu normalnego dla 95% ufności

Teraz zrobimy to samą metodą bootstrap:

B=999
mean.dochod=rep(0,B)
nobs=50 #liczba próbek bootstrapowych - zwiększaj ich liczbę: 50, 250, 500, 1000, 10000

 for (i in 1:B) 
{
 boot.data=sample(dochod,nobs,replace=TRUE)
 mean.dochod[i]=mean(boot.data)
 }

mean_boot <- mean(boot.data)

sd_boot <- sd(boot.data)

se2<-sd(boot.data)/sqrt(length(boot.data))

sep2<-se/mean(boot.data)

hist(boot.data)

plot(density(boot.data))

#przedziałowo bootstrapowo  95% ufności:
dolny<-mean(boot.data)-1.96*se2
gorny<-mean(boot.data)+1.96*se2

Wyniki dla różnych wielkości próby:

Dla wielkości próby ‘r nobs=50’ otrzymujemy: - średnią próbkową: ‘r mean_boot’ - odchylenie standardowe: ‘r sd_boot’ - błąd standardowy: ‘r se2’ - procentowy błąd standardowy: ‘r sep2’ - dolną krawędź przedziału ufności: ‘r dolny’ - górną krawędź przedziału ufności: ‘r gorny’

Porównaj wyniki z obu metod.

Jakie są wnioski? Czy różnice są istotne? Jak wielkość resamplingu wpływa na wyniki?

Zadanie 2. Estymacja błędu standardowego średniej próbkowej metodą bootstrap dla różnych wielkości próby.

A teraz z bootstrapem i gotową funkcją boot:

?boot

mean.boot=function(dochod,idx) {
ans=mean(dochod[idx])
ans
}

DOCHOD.mean.boot = boot(dochod,statistic=mean.boot, R=999)
class(DOCHOD.mean.boot)
## [1] "boot"
names(DOCHOD.mean.boot)
##  [1] "t0"        "t"         "R"         "data"      "seed"      "statistic"
##  [7] "sim"       "call"      "stype"     "strata"    "weights"
DOCHOD.mean.boot
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = dochod, statistic = mean.boot, R = 999)
## 
## 
## Bootstrap Statistics :
##     original     bias    std. error
## t1* 808.5322 -0.2324327    6.201993
plot(DOCHOD.mean.boot)

boot.ci(DOCHOD.mean.boot,conf=0.95,type=c("norm","perc"))
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = DOCHOD.mean.boot, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (796.6, 820.9 )   (796.4, 820.8 )  
## Calculations and Intervals on Original Scale

Testy t studenta

Czy dochody na osobę różnią się istotnie w woj. pomorskim i podkarpackim? Porównaj wyniki testu t-studenta z wynikami testu bootstrapowego.

library(MKinfer) 
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
?boot.t.test

dane2<- dane %>%
  filter(Wojewodztwo %in% c("Pomorskie", "Podkarpackie")) 

  boot.t.test(Dochod_na_osobe~Wojewodztwo, R=999, dane2)
## 
##  Bootstrap Welch Two Sample t-test
## 
## data:  Dochod_na_osobe by Wojewodztwo
## number of bootstrap samples:  999
## bootstrap p-value < 0.001001 
## bootstrap difference of means (SE) = -197.0347 (30.77981) 
## 95 percent bootstrap percentile confidence interval:
##  -258.0512 -139.2055
## 
## Results without bootstrap:
## t = -6.3601, df = 2932.3, p-value = 2.331e-10
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -256.3683 -135.5439
## sample estimates:
## mean in group Podkarpackie    mean in group Pomorskie 
##                   648.2102                   844.1663

Test proporcji

Dla danych “Credit” wykonaj test proporcji dla zmiennych “Student” i “Gender” - czyli sprawdź, czy status studenta różni się istotnie w zależności od płci.

Czy wyniki z włączonym bootstrapem różnią się od wyników pojedynczego testu chi2?

# Przykład 2. Test Chi2 dla dwóch zmiennych jakościowych
library(ISLR)
data("Credit")
?Credit
attach(Credit)
# Czy status studenta (YES, NO) różni się istotnie wg płci (Male, Female)?
tabelka<-table(Student,Gender)
chisq.test(tabelka,simulate.p.value = TRUE, B = 2000)
## 
##  Pearson's Chi-squared test with simulated p-value (based on 2000
##  replicates)
## 
## data:  tabelka
## X-squared = 1.2115, df = NA, p-value = 0.3103
chisq.test(tabelka) # różnice???
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tabelka
## X-squared = 0.87218, df = 1, p-value = 0.3504

Testy ANOVA

Wykonaj i zwizualizuj test ANOVA dla danych “Credit” (z pakietu ISLR) test Anova: czy średni bilans na karcie kredytowej różni się istotnie w zależności od pochodzenia, stanu cywilnego, statusu studenta i płci?

## Analysis of Variance Table
## 
## Response: Balance
##            Df   Sum Sq Mean Sq F value    Pr(>F)    
## Ethnicity   2    18454    9227  0.0463    0.9548    
## Married     1     1332    1332  0.0067    0.9349    
## Student     1  5713181 5713181 28.6378 1.484e-07 ***
## Gender      1     4828    4828  0.0242    0.8765    
## Residuals 394 78602117  199498                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A teraz z włączonym bootstrapem:

library(lmboot)
anova_boot<-ANOVA.boot(Balance ~ Ethnicity + Married + Student + Gender,data=Credit,B=999)
## Warning in ANOVA.boot(Balance ~ Ethnicity + Married + Student + Gender, : This function has only been fully tested for one-way and two-way ANOVA.
anova_boot$`p-values`
## [1] 0.9579580 0.9179179 0.0000000 0.8888889

W przypadku Anovy 1-czynnikowej, możemy wykorzystać pakiet wizualizująco - obliczeniowy “ggstatsplot”. Pakiet ten ma w sobie opcję bootstrappingu, która pozwala na obliczenie wartości p-wartości dla testu ANOVA.

library(ggstatsplot)
## You can cite this package as:
##      Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
##      Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167
ggbetweenstats(data=Credit,
  y=Balance,
  x=Ethnicity,
  nboot=999  #liczba prób bootstrapowych
)

Jakie są wnioski? Czy różnice są istotne? Jak wielkość resamplingu wpływa na wyniki? Jakie są różnice między testem ANOVA a testem ANOVA z bootstrapem?

Wnioski i potencjalne różnice

Metoda klasyczna oszacowała średni dochód na 808.53, z błędem standardowym 6.14 i przedziałem ufności [796.48, 820.58], zakładając normalny rozkład danych.

Metoda bootstrap dała średnią 808.53, z błędem standardowym 6.38 oraz przedziałem ufności [796.1, 821.2] (normalny) i [795.9, 820.5] (percentylowy). Wyniki obu metod są bardzo zbliżone, co oznacza, że dane prawdopodobnie spełniają założenia normalności i nie zawierają znacznych wartości odstających.

W teście t-Studenta różnica dochodów między województwami była istotna zarówno w metodzie klasycznej (p < 0.05), jak i bootstrapowej (p < 0.001), ale bootstrap uwzględnił większą zmienność, dając szerszy przedział ufności.

Kluczowe wnioski:

• Bootstrap potwierdził wyniki klasyczne, co wskazuje, że w tym przypadku metoda parametryczna była wystarczająca.

• Brak istotnej różnicy między metodami sugeruje, że dane są dobrze dopasowane do normalnego rozkładu.

• Bootstrap jest bardziej adekwatny w przypadku danych z silną asymetrią lub dużą zmiennością, ale w tej analizie nie wniósł dodatkowych korzyści poza potwierdzeniem stabilności oszacowań.

Resampling ~ Wynik

Resampling ~ Wynik Przy R = 100 przedziały ufności [794.6, 821.4] dla normalnej i [792.6, 820.2] dla percentylowej są mniej stabilne i bardziej podatne na losowe wahania, co potwierdza ostrzeżenie o niestabilności wyników. Przy R = 999 przedziały są dokładniejsze i lepiej odwzorowują zmienność danych, co wynika z większej liczby próbek. W teście proporcji (Student vs. Gender) bootstrap nie wpłynął znacząco na wynik (p= 0.3504 vs. p=0.3113 ), co sugeruje, że w przypadku danych spełniających klasyczne założenia bootstrap nie zawsze przynosi dodatkowe korzyści.

Anova a Anova Bootstraped

Anova a Anova Bootstraped W klasycznym teście ANOVA p-wartości opierają się na założeniu normalności i homoscedastyczności reszt, co w tym przypadku daje np. p=0.9548 dla Ethnicity i p <0.05 dla Student. W teście ANOVA z bootstrapem wartości p są podobne dla większości zmiennych (p=0.9599 dla Ethnicity, p=0.8629 dla Gender), ale różnica jest bardziej zauważalna dla zmiennej Student (p<0.05), co wskazuje na większą istotność przy resamplingu. Bootstrap jest bardziej odporny na naruszenia założeń klasycznego ANOVA, lepiej odzwierciedlając rozkład danych, ale może prowadzić do różnic w wynikach dla zmiennych o niestandardowym rozkładzie. Warto stosować go, gdy założenia klasycznego ANOVA mogą być naruszone.