Bootstrap
Krzysztof Kowalski
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)
dane <- read.spss("dane.sav", to.data.frame = TRUE)
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ściTeraz zrobimy to samą metodą bootstrap:
B=999
mean.dochod=rep(0,B)
nobs=10000 #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] 823.5821sd_boot <- sd(boot.data)
sd_boot## [1] 1109.98se2<-sd(boot.data)/sqrt(length(boot.data))
se2## [1] 11.0998sep2<-se/mean(boot.data)
sep2## [1] 0.007463213hist(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*se2Wyniki dla różnych wielkości próby:
Dla wielkości próby ‘r nobs=50’ otrzymujemy: - średnią próbkową: ‘745.8542’ - odchylenie standardowe: ‘644.4991’ - błąd standardowy: ‘91.14594’ - procentowy błąd standardowy: ‘0.008240979’ - dolną krawędź przedziału ufności: ‘563.3264’ - górną krawędź przedziału ufności: ‘928.3820’
Dla wielkości próby ‘r nobs=250’ otrzymujemy: - średnią próbkową: ‘753.49’ - odchylenie standardowe: ‘809.33’ - błąd standardowy: ‘51.19’ - procentowy błąd standardowy: ‘0.00816’ - dolną krawędź przedziału ufności: ‘650.33’ - górną krawędź przedziału ufności: ‘856.65’
Dla wielkości próby ‘r nobs=500’ otrzymujemy: - średnią próbkową: ‘829.16’ - odchylenie standardowe: ‘970.36’ - błąd standardowy: ‘43.40’ - procentowy błąd standardowy: ‘0.00741’ - dolną krawędź przedziału ufności: ‘743.51’ - górną krawędź przedziału ufności: ‘914.81’
Dla wielkości próby ‘r nobs=1000’ otrzymujemy: - średnią próbkową: ‘833.93’ - odchylenie standardowe: ‘1014.08’ - błąd standardowy: ‘32.07’ - procentowy błąd standardowy: ‘0.00737’ - dolną krawędź przedziału ufności: ‘771.56’ - górną krawędź przedziału ufności: ‘896.31’
Dla wielkości próby ‘r nobs=10000’ otrzymujemy: - średnią próbkową: ‘820.74’ - odchylenie standardowe: ‘1123.25’ - błąd standardowy: ‘11.23’ - procentowy błąd standardowy: ‘0.00749’ - dolną krawędź przedziału ufności: ‘798.33’ - górną krawędź przedziału ufności: ‘843.14’
Porównaj wyniki z obu metod.
Jakie są wnioski? Czy różnice są istotne? Jak wielkość resamplingu wpływa na wyniki?
Różnice między bootstrapem a metodą klasyczną:
Przy nobs ≤ 250, różnice są wyraźne. Bootstrapowe średnie i błędy standardowe mogą być odchylone od wartości klasycznych. Dla nobs ≥ 1000, bootstrap daje wyniki bardzo zbliżone do klasycznych. Większe nobs zwiększa precyzję estymacji (mniejszy błąd standardowy i węższe przedziały ufności). Małe nobs prowadzi do większej zmienności wyników i szerokich przedziałów ufności.
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## uruchamianie serwera httpd dla pomocy ... wykonanomean.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.1186876 6.18934plot(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.5, 820.8 ) (796.8, 820.5 )
## Calculations and Intervals on Original ScaleTesty 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) ## Warning: pakiet 'MKinfer' został zbudowany w wersji R 4.4.2library(tidyverse)## Warning: pakiet 'tidyverse' został zbudowany w wersji R 4.4.2## Warning: pakiet 'forcats' został zbudowany w wersji R 4.4.2## Warning: pakiet 'lubridate' został zbudowany w wersji R 4.4.2## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ 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.3115 (30.7871)
## 95 percent bootstrap percentile confidence interval:
## -258.5122 -137.4419
##
## 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.1663Zarówno klasyczna metoda t-Studenta, jak i test bootstrapowy wskazują na statystycznie istotną różnicę dochodów między województwami. Różnica średnich i przedziały ufności są bardzo zbliżone w obu metodach. Wyniki bootstrapowe są bardziej odporne na odchylenia od założeń normalności rozkładu dochodów, co czyni je bardziej wiarygodnymi w sytuacji nienormalnych rozkładów. Klasyczna metoda daje bardzo podobne wyniki, co wskazuje, że rozkład danych prawdopodobnie jest zbliżony do normalnego (lub próbka jest wystarczająco duża, aby zastosować centralne twierdzenie graniczne).
Średnie dochody w województwie pomorskim są wyraźnie wyższe niż w podkarpackim. Różnica wynosi około 196 zł na osobę, co jest istotne statystycznie (p<0.001). Przedziały ufności z obu metod (klasycznej i bootstrapowej) obejmują bardzo podobne wartości, co świadczy o zgodności wyników. Test bootstrapowy potwierdza wyniki klasyczne, jednocześnie uwzględniając potencjalne nienormalności w danych.
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)## Warning: pakiet 'ISLR' został zbudowany w wersji R 4.4.2data("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.3098chisq.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, jak i test z symulacją (p>0.05), wskazują brak istotnych różnic między płciami pod względem statusu studenta. Wpływ bootstrappingu:
Wyniki są zbliżone, ale test z symulacją generuje nieco inne statystyki (np. X^2 ) oraz p-value. Różnice wynikają z tego, że symulacja nie stosuje poprawki Yatesa i opiera się na generowaniu losowych próbek.
Wrażliwość na liczebności:
Jeśli w tabeli kontyngencji są małe liczebności, metoda bootstrapowa jest bardziej wiarygodna, ponieważ nie polega na założeniach asymptotycznych (przybliżeniu rozkładu statystyki do rozkładu chi-kwadrat).
Wnioski:
Status studenta (YES/NO) nie różni się istotnie w zależności od płci (p>0.05) na podstawie obu testów. Wyniki testu bootstrapowego są zbliżone do klasycznego testu chi-kwadrat, co sugeruje, że liczebności w tabeli kontyngencji są wystarczająco duże, aby klasyczna metoda była wiarygodna.
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)## Warning: pakiet 'lmboot' został zbudowany w wersji R 4.4.2anova_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.9629630 0.9349349 0.0000000 0.8958959W 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)## Warning: pakiet 'ggstatsplot' został zbudowany w wersji R 4.4.2## 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?
Wyniki ANOVA:
Status studenta (p ≈ 1.484×10^−7) jest jedyną istotną statystycznie zmienną, która wpływa na średni bilans na karcie kredytowej. Pochodzenie (Ethnicity), stan cywilny (Married) i płeć (Gender) nie mają istotnego wpływu.
Wyniki ANOVA z bootstrapem:
Wyniki bootstrapowe są zgodne z klasycznym testem ANOVA: status studenta (Student) jest jedyną zmienną istotną statystycznie (p=0.000). Pochodzenie, stan cywilny i płeć pozostają nieistotne.
Porównanie klasycznego ANOVA i ANOVA z bootstrapem:
Oba podejścia wskazują, że tylko status studenta ma istotny wpływ na bilans. lasyczny ANOVA zakłada normalność reszt i jednorodność wariancji między grupami. Bootstrapowy ANOVA jest bardziej odporny na niezgodność z założeniami i uwzględnia ich naruszenia, co może być przydatne w przypadku nietypowych danych. Zwiększenie liczby prób bootstrapowych (np. 999 do 10,000) może poprawić stabilność wyników, szczególnie gdy dane są niestandardow, ale w tym przypadku wyniki są spójne, co sugeruje, że dane nie naruszają istotnie założeń klasycznego ANOVA.