Bootstrap
Jakub Busłowski
2025-01-20
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)
library(haven)
dane <- read_sav("dane.sav")
attach(dane)
dochod<-Dochod_na_osobe
#punktowo klasycznie
mean(dochod)## [1] 808.5322sd(dochod)## [1] 1064.67se<-sd(dochod)/sqrt(length(dochod))
se## [1] 6.146569sep<-se/mean(dochod)
sep## [1] 0.007602133#przedziałowo klasycznie
mean(dochod)-1.96*se #dolna krawędź## [1] 796.4849mean(dochod)+1.96*se #górna krawędź## [1] 820.5794# 1.96 to kwantyl rozkładu normalnego dla 95% ufnościWyniki klasycznych obliczeń dla estymacji punktowej średniej i względnego błędu standardowego:
Średnia dla zmiennej dochód: 808,53
Odchylenie standardowe: 1064,67
Błąd standardowy: 6,15
Względny błąd standardowy: 0,76%
Wyniki przedziałowo klasycznych:
Dolna krawędź: 796,48
Górna krawędź: 820,58
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)
mean_boot## [1] 1069.232sd_boot <- sd(boot.data)
sd_boot## [1] 1381.043se2<-sd(boot.data)/sqrt(length(boot.data))
se2## [1] 195.309sep2<-se/mean(boot.data)
sep2## [1] 0.005748584hist(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
dolny## [1] 686.4261gorny## [1] 1452.037Wyniki dla metody bootstrap dla 999 próbek i wielkości próby 50:
Średnia bootstrapowa: 783,18
Odchylenie bootstrapowe: 608,3
Błąd standardowy bootstrapowy: 86,03
Względny błąd standardowy bootstrapowy: 0,78%
Dolna krawędź: 614,56
Górna krawędź: 951,79
Wnioski:
Metoda klasyczna okazała się bardziej precyzyjna, co widać po węższym przedziale ufności [796,48;820,58] w porównaniu do metody bootstrapowej [614,56;951,79]. Średnia klasyczna wyniosła 808,53, co jest wyższe niż średnia bootstrapowa 783,18. Względny błąd standardowy w metodzie klasycznej to 0,76%, natomiast w metodzie bootstrapowej wyniósł 0,78%, co wskazuje na większą zmienność w wynikach bootstrapu. Chociaż bootstrap uwzględnia większą losowość i lepiej oddaje niestandardowe rozkłady danych, wiąże się z większym błędem standardowym (86,03 w porównaniu do 6,15 w metodzie klasycznej). Metoda klasyczna jest bardziej precyzyjna, pod warunkiem spełnienia założeń o normalnym rozkładzie danych, natomiast bootstrap zapewnia większą elastyczność kosztem szerszych przedziałów ufności i większej zmienności wyników.
B3=999
mean.dochod3=rep(0,B3)
nobs3=250 #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)
mean_boot3## [1] 724.2899sd_boot3 <- sd(boot.data3)
sd_boot3## [1] 733.4789se3<-sd(boot.data3)/sqrt(length(boot.data3))
se3## [1] 46.38928sep3<-se3/mean(boot.data3)
sep3## [1] 0.06404794hist(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
dolny3## [1] 633.3669gorny3## [1] 815.2129Wyniki dla metody bootstrap dla 999 próbek i wielkości próby 250:
Średnia bootstrapowa: 798,12
Odchylenie bootstrapowe: 1108,5
Błąd standardowy bootstrapowy: 70,11
Względny błąd standardowy bootstrapowy: 8,78%
Dolna krawędź: 660,71
Górna krawędź: 935,54
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 250), wartości: średnia bootstrapowa staje się bardziej stabilna (798,12 vs. 757,13), odchylenie standardowe wzrasta (1108,5 vs. 608,3), błąd standardowy bootstrapowy również rośnie (70,11 vs. 86,03). Przedziały ufności stają się węższe [660,71;935,54] vs. [614,56;951,79].
Wniosek: Zwiększenie liczby prób resamplingowych poprawia stabilność wyników, ale powoduje wzrost błędów standardowych, co może prowadzić do większej zmienności w oszacowaniach.
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.3195662 6.145183plot(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.9 ) (795.8, 820.4 )
## Calculations and Intervals on Original ScaleBootstrapowy błąd standardowy (6,30) jest zbliżony do klasycznego (6,15), a przedziały ufności (normalny: 796,3–821,0, percentylowy: 796,5–821,0) są niemal identyczne. Minimalny bias (−0,16) potwierdza stabilność i zgodność wyników bootstrap z klasycznymi oszacowaniami.
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(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniondane$Wojewodztwo <- as.character(dane$Wojewodztwo)
dane$Wojewodztwo_Nazwa <- recode(dane$Wojewodztwo,
"2" = "Dolnośląskie",
"4" = "Kujawsko-pomorskie",
"6" = "Lubelskie",
"8" = "Lubuskie",
"10" = "Łódzkie",
"12" = "Małopolskie",
"14" = "Mazowieckie",
"16" = "Opolskie",
"18" = "Podkarpackie",
"20" = "Podlaskie",
"22" = "Pomorskie",
"24" = "Śląskie",
"26" = "Świętokrzyskie",
"28" = "Warmińsko-mazurskie",
"30" = "Wielkopolskie",
"32" = "Zachodniopomorskie")library(MKinfer)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.5
## ✔ ggplot2 3.5.1 ✔ stringr 1.5.1
## ✔ lubridate 1.9.4 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── 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 errorsdane2 <- dane %>%
filter(Wojewodztwo_Nazwa %in% c("Pomorskie", "Podkarpackie"))
boot.t.test(Dochod_na_osobe ~ Wojewodztwo, R = 999, data = 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.5374 (30.78041)
## 95 percent bootstrap percentile confidence interval:
## -252.6584 -132.9560
##
## 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 18 mean in group 22
## 648.2102 844.1663Wartość p wskazuje na istotne statystycznie różnice w dochodach między województwami. Średni dochód w województwie pomorskim jest wyższy o około 196,82 zł na osobę. Test bootstrapowy potwierdza te różnice, a jego wyniki są zgodne z wynikami klasycznego testu t. Zarówno test t-Studenta, jak i bootstrapowy wykazują bardzo niskie wartości p, co wskazuje na istotność różnic. Przedziały ufności obu metod są zbliżone, co potwierdza ich spójność.
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.3258chisq.test(tabelka) # różnice???##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tabelka
## X-squared = 0.87218, df = 1, p-value = 0.3504Zarówno klasyczny test chi-kwadrat, jak i test bootstrapowy pokazują, że nie ma istotnych statystycznie różnic w statusie studenta w zależności od płci. Wysokie wartości p (chi-kwadrat: 0,3504, bootstrap: 0,3298) wskazują na brak podstaw do odrzucenia hipotezy zerowej. Wyniki obu testów są bardzo zbliżone, co potwierdza ich spójność.
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 ' ' 1A 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.9539540 0.9419419 0.0000000 0.8798799W 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.03167ggbetweenstats(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?
Klasyczny test ANOVA i jego wersja bootstrapowa pokazują, że nie ma istotnych różnic w średnim bilansie na karcie kredytowej w zależności od pochodzenia, stanu cywilnego, statusu studenta i płci. Wysokie wartości p dla większości zmiennych potwierdzają brak podstaw do odrzucenia hipotezy zerowej. Jedynym wyjątkiem jest zmienna „Student”, dla której oba testy wykazały istotne różnice (p < 0,001). Zwiększenie liczby prób bootstrapowych nie zmienia wyników – wartości p i wnioski są zgodne z klasycznym testem ANOVA. Bootstrap potwierdza wyniki klasycznego testu i daje dodatkową pewność, szczególnie w przypadku mniejszych prób lub niestandardowych danych.