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.

library(foreign)
dane <- read.spss("C:/Users/Mateusz/Desktop/Semestr III/Nieklasyczne metody statystyki/Zajęcia 6/dane.sav", to.data.frame = TRUE)
attach(dane)
dochod<-Dochod_na_osobe

#punktowo klasycznie
print(paste("średnia próbkowa =", mean(dochod)))

## [1] "średnia próbkowa = 808.532162306147"

print(paste("odchylenie standardowe =", sd(dochod)))

## [1] "odchylenie standardowe = 1064.67018546966"

se<-sd(dochod)/sqrt(length(dochod))
print(paste("błąd standardowy =", se))

## [1] "błąd standardowy = 6.14656886103053"

sep<-se/mean(dochod)
print(paste("procentowy błąd standardowy =", sep))

## [1] "procentowy błąd standardowy = 0.00760213278776553"

#przedziałowo klasycznie
dolna <- mean(dochod)-1.96*se   #dolna krawędź
print(paste("dolna krawędź przedziału ufności =", dolna))

## [1] "dolna krawędź przedziału ufności = 796.484887338527"

gorna <- mean(dochod)+1.96*se   #górna krawędź
print(paste("górna krawędź przedziału ufności =", gorna))

## [1] "górna krawędź przedziału ufności = 820.579437273767"

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

Teraz zrobimy to samą metodą bootstrap:

set.seed(123) 

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 <- round(mean(boot.data),digits=2)
print(paste("średnia próbkowa =", mean_boot))

## [1] "średnia próbkowa = 777.26"

sd_boot <- sd(boot.data)
print(paste("odchylenie standardowe =", sd_boot))

## [1] "odchylenie standardowe = 807.986343605106"

se2<-sd(boot.data)/sqrt(length(boot.data))
print(paste("błąd standardowy =", se2))

## [1] "błąd standardowy = 114.266524533859"

sep2<-se2/mean(boot.data)
print(paste("procentowy błąd standardowy =", sep2))

## [1] "procentowy błąd standardowy = 0.147011112514492"

#przedziałowo bootstrapowo  95% ufności:
dolny<-mean(boot.data)-1.96*se2
print(paste("dolna krawędź przedziału ufności =", dolny))

## [1] "dolna krawędź przedziału ufności = 553.302150487856"

gorny<-mean(boot.data)+1.96*se2
print(paste("górna krawędź przedziału ufności =", gorny))

## [1] "górna krawędź przedziału ufności = 1001.22692666058"

hist(boot.data)

plot(density(boot.data))

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

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

Estymacja dla różnych liczebności próby

set.seed(123) 
B=999
nobs_list <- c(50, 250, 500, 1000, 10000)

bootstrap_statystyki <- function(nobs) {
  mean.dochod <- rep(0, B)
  for (i in 1:B) {

boot.data <- sample(dochod, nobs, replace = TRUE)
  }
  mean_boot <- round(mean(boot.data), digits = 2)
  sd_boot <- sd(boot.data)
  se2 <- sd(boot.data) / sqrt(length(boot.data))
  sep2 <- se2 / mean(boot.data)
  dolny <- mean(boot.data) - 1.96 * se2
  gorny <- mean(boot.data) + 1.96 * se2
  return(c(mean_boot, sd_boot, se2, sep2, dolny, gorny))
}

wyniki <- sapply(nobs_list, bootstrap_statystyki)

rownames(wyniki) <- c(
  "Średnia próbkowa",
  "Odchylenie standardowe",
  "Błąd standardowy",
  "Procentowy błąd standardowy",
  "Dolna krawędź przedziału ufności",
  "Górna krawędź przedziału ufności"
)
colnames(wyniki) <- paste("n =", nobs_list)

wyniki <- round(wyniki, 4)  # Zaokrąglenie dla czytelności
print(wyniki)

##                                     n = 50   n = 250  n = 500  n = 1000
## Średnia próbkowa                  777.2600  852.5500 787.7900  779.0100
## Odchylenie standardowe            807.9863 1425.8677 889.8804 1022.7930
## Błąd standardowy                  114.2665   90.1798  39.7967   32.3436
## Procentowy błąd standardowy         0.1470    0.1058   0.0505    0.0415
## Dolna krawędź przedziału ufności  553.3022  675.7981 709.7849  715.6171
## Górna krawędź przedziału ufności 1001.2269 1029.3029 865.7878  842.4038
##                                  n = 10000
## Średnia próbkowa                  792.1400
## Odchylenie standardowe           1076.0322
## Błąd standardowy                   10.7603
## Procentowy błąd standardowy         0.0136
## Dolna krawędź przedziału ufności  771.0489
## Górna krawędź przedziału ufności  813.2293

Wykresy gęstości metody bootstrap dla pozostałych wielkości próby

set.seed(123) 
par(mfrow = c(2, 2))
nobs_list <- c(250, 500, 1000, 10000)

for (nobs in nobs_list) {
  mean.dochod <- rep(0, B)
  for (i in 1:B) {
    boot.data <- sample(dochod, nobs, replace = TRUE)
    mean.dochod[i] <- mean(boot.data)
  }
  
  # Wykres gęstości dla aktualnej próby
  plot(density(mean.dochod),
       main = paste("n =", nobs),
       xlab = "Średnia próbkowa",
       ylab = "Gęstość")
}

Porównaj wyniki z obu metod.

Bootstrapowa średnia próbkowa (777.26) jest niższa od klasycznej (808.53).
Może to wynikać z losowości próbek w resamplingu, (szczególnie przy n=50), ponieważ metoda klasyczna operuje na pełnym zbiorze danych, podczas gdy bootstrap bazuje na losowych podpróbkach.
Bootstrapowy błąd standardowy (114.27) jest znacznie większy od klasycznego (6.15). To oznacza, że metoda bootstrap przyn=50 charakteryzuje się większą zmiennością wyników w porównaniu do metody klasycznej.
Procentowy błąd standardowy potwierdza te różnice (bootstrap: 0.1470, klasyczna: 0.0076).
Przedział ufności bootstrap róWnież jest znacznie szerszy, co wskazuje na większą niepewność oszacowań przy użyciu bootstrapu.

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

Różnice w oszacowaniach (średnia, błąd standardowy, przedział ufności) są istotne, szczególnie przy małej próbie (n=50).
Metoda bootstrap opiera się na losowym próbkowaniu, co zwiększa zmienność wyników.
Przy większych próbach (n=10000) błąd standardowy bootstrap znacząco maleje (10.76), a przedział ufności staje się węższy ([771.05,813.23]).
Średnia bootstrap staje się również bliższa klasycznej średniej próbkowej. Większe próby redukują wpływ losowości w resamplingu.
W tym przypadku klasyczna metoda wypada lepiej niż bootstrap przy dużej próbie, natomiast należy zauważyć, że wraz ze wzrostem liczebności prób, wyniki bootstrap zbliżają się do wyników klasycznych.

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:

set.seed(123) 
# ?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.08768945    6.035462

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.8, 820.4 )   (796.4, 819.9 )  
## Calculations and Intervals on Original Scale

Testy t studenta

set.seed(123) 
library(MKinfer) 
library(tidyverse)
# ?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) = -195.2333 (30.75198) 
## 95 percent bootstrap percentile confidence interval:
##  -258.3369 -136.9888
## 
## 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

Czy dochody na osobę różnią się istotnie w woj. pomorskim i podkarpackim?

W klasycznym teście t-studenta p-value = 2,331*10^−10, co jest znacznie mniejsze niż typowy poziom istotności alfa=0,05.
Dla bootstrapowej werstji testu p-value<0.001001, co również wskazuje na istotne różnice.
W obu przypadkach wyniki jednoznacznie wskazują, że dochody na osobę w województwach pomorskim i podkarpackim różnią się istotnie.

Porównaj wyniki testu t-studenta z wynikami testu bootstrapowego.

Obie metody wskazują na bardzo małą p-wartość, co oznacza istotną różnicę między grupami.
Należy jednak zauważyć, że w przypadku klasycznego testu wartość ta jest znacznie niższa.
Klasyczny test t-Studenta i bootstrapowy test różnicy średnich uzyskują również bardzo podobne przedziały ufności:
- Klasyczny: (-256.3683, -135.5439)
- Bootstrapowy: (-258.3369, -136.9888)

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.

# Test Chi2 dla dwóch zmiennych jakościowych
library(ISLR)
set.seed(123) 

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

chisq.test(tabelka) 

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tabelka
## X-squared = 0.87218, df = 1, p-value = 0.3504

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

W przypadku klasycznego testu chi-kwadrat p-value=0.3504, co wskazuje, że nie ma istotnych różnic w statusie studenta w zależności od płci (przy poziomie istotności alfa=0,05).
W przypadku testu chi-kwadrat z symulacją bootstrap p-value= 0.3138, co róWnież wskazuje na brak istotnych różnic między grupami.
Zarówno klasyczny test chi-kwadrat, jak i test z symulowanym p-value (bootstrap), wskazują brak istotnych różnic.

Dodatkowy wykres pokazujący tą zależność

library(ggplot2)

tabelka_df <- as.data.frame(tabelka)
colnames(tabelka_df) <- c("Student", "Gender", "Count")

ggplot(tabelka_df, aes(x = Gender, y = Count, fill = Student)) +
  geom_bar(stat = "identity", position = "fill") +
  labs(title = "Proporcje statusu studenta według płci",
       y = "Proporcja",
       x = "Płeć")

Jak widać, faktycznie istotne różnice pomiędzy grupami nie występują.

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?

attach(Credit)
model1 <- lm(Balance ~ Ethnicity + Married + Student + Gender, data=Credit)
anova(model1)

## 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

plot(model1)

Wyniki testu Annova:

Zmienna Ethnicity (Pochodzenie) p=0.9548 (>0,05)
Brak istotnych różnic w średnim bilansie na karcie kredytowej w zależności od pochodzenia.

Zmienna Married (Stan cywilny) p=0.9349 (>0,05)
Brak istotnych różnic w średnim bilansie na karcie kredytowej w zależności od stanu cywilnego.

Zmienna Student (Status studenta) p=1.484e-07 (<0,05)
Występują istotne różnice w średnim bilansie na karcie kredytowej między studentami a osobami niebędącymi studentami.

Zmienna Gender (Płeć) p=0.8765 (>0,05)
Brak istotnych różnic w średnim bilansie na karcie kredytowej w zależności od płci.

Podsumowując jedynie status studenta wpływa istotnie na średni bilans na karcie kredytowej.
Pozostałe zmienne (pochodzenie, stan cywilny, płeć) nie mają istotnego wpływu.

Dodatkowa wizualizacja

library(ggpubr)
library(patchwork)

par(mfrow = c(2, 2))

# Wizualizacja modelu ANOVA
p1<-ggboxplot(Credit, x = "Ethnicity", y = "Balance", color = "Ethnicity", 
          add = "jitter", palette = "jco") +
  stat_compare_means(method = "anova", label = "p.format") +
   labs(x = "Pochodzenie",y = "Bilans")

p2<-ggboxplot(Credit, x = "Married", y = "Balance", color = "Married", 
          add = "jitter", palette = "jco") +
  stat_compare_means(method = "anova", label = "p.format") +
   labs(x = "Stan cywilny",y = "Bilans")

p3<-ggboxplot(Credit, x = "Student", y = "Balance", color = "Student", 
          add = "jitter", palette = "jco") +
  stat_compare_means(method = "anova", label = "p.format") +
   labs(x = "Status studenta",y = "Bilans")

p4<-ggboxplot(Credit, x = "Gender", y = "Balance", color = "Gender", 
          add = "jitter", palette = "jco") +
  stat_compare_means(method = "anova", label = "p.format") +
  labs(x = "Płeć",y = "Bilans")



(p1 | p2) / (p3 | p4)

Faktycznie największe zróżnicowanie rozkładów widać w przypadku zmiennej Student.

A teraz z włączonym bootstrapem:

set.seed(123) 

library(lmboot)
anova_boot<-ANOVA.boot(Balance ~ Ethnicity + Married + Student + Gender,data=Credit,B=999)
anova_boot$`p-values`

## [1] 0.9489489 0.9269269 0.0000000 0.8718719

W przypadku testu Annova z bootstrapem ponownie można wyciągnąć te same wnioski na podstawie uzyskanych wartości p-value:
- Brak istotnych różnic w Balance w zależności od Ethnicity. (p-value=0.9489489 >0,05)
- Brak istotnych różnic w Balance w zależności od Married. (p-value=0.9269269 >0,05)
- Występują istotnych różnice w wartościach zmiennej Balance w zależności od Student. (p-value~0 <0,05)
- Brak istotnych różnic w Balance w zależności od Gender. (p-value=0.8718719 >0,05)

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)
attach(Credit)

ggbetweenstats(data=Credit,
  x=Ethnicity,y=Balance,
  nboot=999  #liczba prób bootstrapowych
)

ggbetweenstats(data=Credit,
  x=Student,y=Balance,
  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

Wyniki testów, zarówno klasycznego ANOVA, jak i z bootstrapem, oraz wizualizacja z pakietu ggbetweenstats jednoznacznie wskazują, że różnice w średnich wartości zmiennej “Balance” pomiędzy grupami etnicznymi nie są istotne statystycznie. Resampling metodą bootstrap nie wpłynął istotnie na wyniki, co sugeruje wiarygodność analizy.

Dla zmiennej “Student” różnice w średniej wartości zmiennej “Balance” są istotne statystycznie i mają duże znaczenie. Studenci mają znacznie wyższe salda kredytowe niż osoby, które nie są studentami.

Pomimo braku wyraźnnych różnic poimędzy testem ANNOVA, a ANNOVA z bootstrapem, należy zwrócić uwagę, że zwiększenie liczby prób bootstrapowych (np. nboot = 999) może istotnie wpłynąć na większą stabilność wyników. Może to być szczególnie przydatne, gdy dane są nieregularne, próby małe, a założenia klasycznego ANOVA są wątpliwe.