Bootstrap

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.

Pobranie danych

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

Klasycznie

Klasyczne obliczenia estymacji punktowej dla średniej i jej względnego błędu standardowego.

#punktowo klasycznie
mean_value<-mean(dochod) 
sd_value<-sd(dochod) 
se<-sd(dochod)/sqrt(length(dochod)) 
sep<-se/mean(dochod)

## Średnia dla zmiennej dochód wynosi: 808.53

## Odchylenie standardowe wynosi: 1064.67

## Błąd standardowy wynosi: 6.1466

## Względny błąd standardowy wynosi: 0.76 %

#przedziałowo klasycznie
# 1.96 to kwantyl rozkładu normalnego dla 95% ufności

# Klasyczny przedział ufności
lower_bound <- mean(dochod) - 1.96 * se
upper_bound <- mean(dochod) + 1.96 * se

## 95% przedział ufności dla średniej wynosi:

## Dolna krawędź: 796.48

## Górna krawędź: 820.58

Metoda 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 metody bootstrap dla 999 próbek i wielkości próby 50 :

## Średnia bootstrapowa: 1049.98

## Odchylenie standardowe bootstrapowe: 1400.97

## Błąd standardowy bootstrapowy: 198.1277

## 95% bootstrapowy przedział ufności:

## Dolna krawędź: 661.65

## Górna krawędź: 1438.31

Porównaj wyniki z obu metod.

Wnioski: Obie metody wskazują zbliżoną średnią (808.53 vs. 796.85), co sugeruje, że metoda bootstrap dobrze oszacowuje wartość średniej próbkowej. Odchylenie standardowe w metodzie bootstrap jest wyższe (1110.97 vs. 1064.67) i podobnie błąd standardowy jest znacznie większy w metodzie bootstrap (157.11 vs. 6.15). Może to wynikać z większej wariancji w próbkach bootstrapowych. Klasyczna metoda daje znacznie węższy przedział ufności (796.48 – 820.58) w porównaniu do bootstrapu (488.91 – 1104.80). Wynika to z tego, że bootstrap uwzględnia więcej niepewności.

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

 for (i in 1:B3) 
{
 boot.data3=sample(dochod,nobs3,replace=TRUE)
 mean.dochod3[i]=mean(boot.data3)
 }

mean_boot3 <- mean(boot.data3)

sd_boot3 <- sd(boot.data3)

se3<-sd(boot.data3)/sqrt(length(boot.data3))

sep3<-se3/mean(boot.data3)

hist(boot.data3)

plot(density(boot.data3))

#przedziałowo bootstrapowo  95% ufności:
dolny3<-mean(boot.data3)-1.96*se3
gorny3<-mean(boot.data3)+1.96*se3

## Wyniki metody bootstrap dla 999 próbek i wielkości próby 500 :

## Średnia bootstrapowa: 787.39

## Odchylenie standardowe bootstrapowe: 874.47

## Błąd standardowy bootstrapowy: 39.1074

## 95% bootstrapowy przedział ufności:

## Dolna krawędź: 710.74

## Górna krawędź: 864.04

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

Gdy wielkość próby bootstrapowej rośnie (z 50 do 500), wartości: Odchylenia standardowego znacząco maleją (1110.97 vs. 798.65), Błąd standardowy również spada (157.11 vs. 35.72), Przedział ufności staje się węższy (488.91–1104.80 vs. 660.77–800.77).

Wniosek: Większe próby pozwalają na lepszą stabilizację wyników, ponieważ redukują wariancję oszacowań.

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
}

#Uruchomienie procedury bootstrap:
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"

#Statystyki bootstrapowe:
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.2160556    6.301564

plot(DOCHOD.mean.boot)

#Przedziały ufności:
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.0, 820.7 )   (796.5, 821.0 )  
## Calculations and Intervals on Original Scale

Bootstrapowy błąd standardowy (6.067342) jest bardzo podobny do tego obliczonego klasycznie(6.1466), co sugeruje, że bootstrap dobrze oszacowuje niepewność nawet przy założeniu niestandardowego rozkładu danych.

Przedziały ufności dla obu metod (normalnej i percentylowej) są prawie identyczne, co wskazuje na stabilność wyników bootstrap.

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)

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) = -196.021 (30.7244) 
## 95 percent bootstrap percentile confidence interval:
##  -256.9990 -138.5114
## 
## 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

Wnioski: Zarówno wyniki testu t-studenta, jak i testu bootstrapowego wskazują na to, że dochody na osobę istotnie różnią się w województwach pomorskim i podkarpackim. Średnie dochody na osobę w Podkarpackim wynoszą 648,21 zł, zaś w Pomorskim 844,16 zł.Test t-Studenta wskazuje na istotną statystycznie różnicę średnich dochodów między województwami (p-value = 2.331e-10 < 0.05). Przedział ufności wynosi od -256,36 do -135,54 w związku z tym wyklucza 0, co również potwierdza istnienie różnicy. P-value testu bootstrapowego również jest bardzo niskie (poniżej 0,001), co wskazuje na istotną różnicę w dochodach. Wygenerowany przedział ufności jest podobny, lecz nieco węższy- od -251,61 zł do -139,30 zł.

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.3063

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

Wnioski: P-value zarówno w teście chi2 z bootstrapem jak i bez jest większe od poziomu istotności 0,05 (p-value = 0,31 i 0,35), a więc nie mamy podstaw do odrzucenia hipotezy zerowej. Możemy przyjąć, że zmienne “Student” i “Gender” są niezależne. Status studenta nie różni się istotnie w zależności od płci. Wyniki z włączonym bootstrapem w niewielkim stopniu różnią się od wyników pojedynczego testu chi2, oba prowadzą do tego samego wniosku.

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)
anova_boot$`p-values`

## [1] 0.9499499 0.9189189 0.0000000 0.8858859

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
)

library(ggstatsplot)
ggbetweenstats(data=Credit,
  y=Balance,
  x=Student,
  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: Test Anova wskazuje na to, że jedynie status studenta istotnie wpływa na średni bilans na karcie kredytowej. Pozostałe zmienne takie jak pochodzenie, stan cywilny i płeć nie są istotne statystycznie. Włączenie bootstrapu nie spowodowało żeby zmienne, takie jak pochodzenie, stan cywilny i płeć stały się istotne, nadal ich p-value jest wyższe od poziomu istotności. Tak samo jak bez bootstrapu jedynie status studenta istotnie wpływa na średni bilans na karcie kredytowej. Na wizualizacji można zauważyć, że osoby pochodzenia African American, Asian i Caucasian posiadają średni bilans na podobnym poziomie, zaś kolejna wizualizacja pozwala dostrzec, że osoby posiadające status studenta mają wyższy średni bilans na karcie kredytowej niż osoby nie posiadające statusu studenta. Różnice między testem ANOVA a testem ANOVA z bootstrapem są niewielkie, ale przy bootstrapie zmienne nieistotne są jeszcze bardziej nieistotne (mają odrobinę wyższe p-value), zaś zmienna istotna ma jeszcze niższe p-value równe 0 w porównaniu z wartością 1.484e-07 w zwykłym teście ANOVA. Włączenie bootstrapu potwierdza poprawność wcześniejszych wniosków.