Bootstrap
Veronika Zhdamarova
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("C:/Users/Veronika Zhdamarova/Desktop/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=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*se2B=999
mean.dochod=rep(0,B)
nobs=250 #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*se2B=999
mean.dochod=rep(0,B)
nobs=500 #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*se2B=999
mean.dochod=rep(0,B)
nobs=1000 #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 średniej próbkowej metodą bootstrap: - średnią próbkową: ‘r mean_boot’ = 808.5322 - odchylenie standardowe: ‘r sd_boot’ = 1064.67 - błąd standardowy: ‘r se2’ = 6.146569 - procentowy błąd standardowy: ‘r sep2’ = 0.007602133 - dolną krawędź przedziału ufności: ‘r dolny’ = 796.4849 - górną krawędź przedziału ufności: ‘r gorny’ = 820.5794
#Wyniki dla różnych wielkości próby:
Dla wielkości próby ‘r nobs=50’ otrzymujemy:
- średnią próbkową: ‘r mean_boot’ = 896.8864
- odchylenie standardowe: ‘r sd_boot’ = 1386.738
- błąd standardowy: ‘r se2’ = 196.1143
- procentowy błąd standardowy: ‘r sep2’ = 0.00685323
- dolną krawędź przedziału ufności: ‘r dolny’ = 512.5023
- górną krawędź przedziału ufności: ‘r gorny’ = 1281.27
Dla wielkości próby ‘r nobs=250’ otrzymujemy:
- średnią próbkową: ‘r mean_boot’ = 762.7767
- odchylenie standardowe: ‘r sd_boot’ = 861.4567
- błąd standardowy: ‘r se2’ = 54.4833
- procentowy błąd standardowy: ‘r sep2’ = 0.00805815
- dolną krawędź przedziału ufności: ‘r dolny’ = 655.9894
- górną krawędź przedziału ufności: ‘r gorny’ = 869.564
Dla wielkości próby ‘r nobs=500’ otrzymujemy:
- średnią próbkową: ‘r mean_boot’ = 822.8599
- odchylenie standardowe: ‘r sd_boot’ = 1301.7
- błąd standardowy: ‘r se2’ = 58.21381
- procentowy błąd standardowy: ‘r sep2’ = 0.007469764
- dolną krawędź przedziału ufności: ‘r dolny’ = 708.7608
- górną krawędź przedziału ufności: ‘r gorny’ = 936.9589
Dla wielkości próby ‘r nobs=1000’ otrzymujemy:
- średnią próbkową: ‘r mean_boot’ = 821.4826
- odchylenie standardowe: ‘r sd_boot’ = 1107.904
- błąd standardowy: ‘r se2’ = 35.035
- procentowy błąd standardowy: ‘r sep2’ = 0.007482288
- dolną krawędź przedziału ufności: ‘r dolny’ = 752.814
- górną krawędź przedziału ufności: ‘r gorny’ = 890.1512
Zapiszę wyniki w postaci tabeli dla łątwiejszego porówania:
# Wyniki bootstrap w postaci tabeli z dodatkowymi wynikami
library(knitr)## Warning: pakiet 'knitr' został zbudowany w wersji R 4.3.2# Tworzenie ramki danych z wynikami
bootstrap_table <- data.frame(
Wielkość_Proby = c("50", "250", "500", "1000", "Średnia próbkowa"),
Średnia_Próbkowa = c(896.8864, 762.7767, 822.8599, 821.4826, 808.5322),
Odchylenie_Standardowe = c(1386.738, 861.4567, 1301.7, 1107.904, 1064.67),
Błąd_Standardowy = c(196.1143, 54.4833, 58.21381, 35.035, 6.146569),
Procentowy_Błąd = c(0.00685323, 0.00805815, 0.007469764, 0.007482288, 0.007602133),
Dolna_Granica = c(512.5023, 655.9894, 708.7608, 752.814, 796.4849),
Górna_Granica = c(1281.27, 869.564, 936.9589, 890.1512, 820.5794)
)
# Tworzenie tabeli w formacie Markdown
kable(bootstrap_table,
col.names = c("Wielkość próby",
"Średnia prób.",
"Odch. stand.",
"Błąd stand.",
"Proc. błąd",
"Dolna granica",
"Górna granica"),
caption = "Wyniki bootstrap dla różnych wielkości próby i średniej próbkowej")| Wielkość próby | Średnia prób. | Odch. stand. | Błąd stand. | Proc. błąd | Dolna granica | Górna granica |
|---|---|---|---|---|---|---|
| 50 | 896.8864 | 1386.7380 | 196.114300 | 0.0068532 | 512.5023 | 1281.2700 |
| 250 | 762.7767 | 861.4567 | 54.483300 | 0.0080582 | 655.9894 | 869.5640 |
| 500 | 822.8599 | 1301.7000 | 58.213810 | 0.0074698 | 708.7608 | 936.9589 |
| 1000 | 821.4826 | 1107.9040 | 35.035000 | 0.0074823 | 752.8140 | 890.1512 |
| Średnia próbkowa | 808.5322 | 1064.6700 | 6.146569 | 0.0076021 | 796.4849 | 820.5794 |
###Wnioski: Średnia próbkowa z obu metod jest zbliżona, co wskazuje na poprawność zastosowania metody bootstrap. Błąd standardowy w metodzie bootstrap może być trochę większy przy małej liczbie próbek.
Wraz ze wzrostem liczby próbek bootstrap (np. z 50 do 1000), przedział ufności staje się węższy, a oszacowania bardziej stabilne. Dla bardzo małej liczby próbek wyniki mogą być niestabilne. Różnice w wynikach są niewielkie, czyli metoda bootstrap jest dobrą alternatywą do klasycznych metod dla takich zadań.
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.2157733 6.117376plot(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.7 ) (796.9, 820.8 )
## 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.3.3library(tidyverse)## Warning: pakiet 'ggplot2' został zbudowany w wersji R 4.3.2## Warning: pakiet 'readr' został zbudowany w wersji R 4.3.2## Warning: pakiet 'dplyr' został zbudowany w wersji R 4.3.2## Warning: pakiet 'stringr' został zbudowany w wersji R 4.3.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.4.4 ✔ 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) = -195.5965 (30.78106)
## 95 percent bootstrap percentile confidence interval:
## -253.6218 -137.7195
##
## 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.1663Odpowiedzi na pytania: Wartość p wskazuje na bardzo istotne statystycznie różnice w dochodach między województwami. Średnia dochodów w województwie pomorskim jest wyższa o około 195.95 zł na osobę. Test bootstrapowy również wskazuje na bardzo istotne różnice w dochodach, z wartościami zbliżonymi do tych uzyskanych w klasycznym teście t. Test t-Studenta, jak i test bootstrapowy wykazują bardzo niskie wartości p, co oznacza, że różnice między województwami są istotne statystycznie. Przedziały ufności są zbliżone, co potwierdza zgodność obu metod.
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.3.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.3208chisq.test(tabelka) # różnice???##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tabelka
## X-squared = 0.87218, df = 1, p-value = 0.3504Odpowiedź na pytanie: W obu testach p-value jest większe od 0,05, co oznacza brak istotnych statystycznie różnic w proporcjach między statusami studenta w zależności od płci.
Testy ANOVA
Wykonaj i zwizualizuj test ANOVA dla danych “Credit” (z pakietu ISLR) test Anova.
## 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 ' ' 1Czy średni bilans na karcie kredytowej różni się istotnie w zależności od pochodzenia, stanu cywilnego, statusu studenta i płci?
Odpowiedź na pytanie: Pochodzenie:p=0.9548 – brak istotnych różnic w średnim bilansie w zależności od pochodzenia. Stan cywilny: p=0.9349 – brak istotnych różnic w średnim bilansie w zależności od stanu cywilnego. Status studenta: p=1.484e-07 – istotne różnice w średnim bilansie między studentami a niestudentami. Płeć: p=0.8765 – brak istotnych różnic w średnim bilansie w zależności od płci. Można wywnioskować, że tylko status studenta wpływa istotnie na średni bilans na karcie kredytowej.
A teraz z włączonym bootstrapem:
library(lmboot)## Warning: pakiet 'lmboot' został zbudowany w wersji R 4.3.3anova_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.9159159 0.0000000 0.8858859Wyniki bootstrapu: Jedyna zmienna, która jest istotna statystycznie i wpływa na średni bilans na karcie kredytowej, to posiadanie statusu studenta
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)## Warning: pakiet 'ggstatsplot' został zbudowany w wersji R 4.3.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?
Odpowiedzi: Większa liczba resamplingu zwiększa dokładność oszacowania wartości p w teście bootstrapowym.
Klasyczny ANOVA zakłada normalność rozkładu i jednorodność wariancji. Ten test jest szybszy w obliczeniach, ale mniej dokładny w przypadku odstępstw od założeń. ANOVA z bootstrapem nie wymaga spełnienia założeń normalności. Jest bardziej dokładny.