Bootstrap
Julia Bogusz
2025-01-19
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/julka/Downloads/dane.sav", to.data.frame = TRUE)
attach(dane)
dochod<-Dochod_na_osobe#punktowo klasycznie
mean_k <- mean(dochod)
sd_k <- sd(dochod)
se_k<-sd(dochod)/sqrt(length(dochod))
se_k## [1] 6.146569sep_k<-se_k/mean(dochod)
sep_k## [1] 0.007602133#przedziałowo klasycznie
dolny_k <- mean(dochod)-1.96*se_k #dolna krawędź
gorny_k <- mean(dochod)+1.96*se_k #górna krawędź
dolny_k## [1] 796.4849gorny_k## [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=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<-se2/mean(boot.data)
hist(boot.data)
mean_boot## [1] 749.8268sd_boot## [1] 828.1259se2## [1] 37.03491sep2## [1] 0.04939129plot(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] 677.2384gorny## [1] 822.4152Wyniki dla różnych wielkości próby:
Dla wielkości próby ‘r nobs=50’ otrzymujemy: - średnią próbkową: 677.4588 - odchylenie standardowe: 845.0496 - błąd standardowy: 119.5081 - procentowy błąd standardowy: 0.009072978 - dolną krawędź przedziału ufności: 443.223 - górną krawędź przedziału ufności: 911.6946
Porównaj wyniki z obu metod.
metrics_table <- data.frame(
Metryka = c("Średnia", "Odchylenie standardowe", "Błąd standardowy", "SEP",
"Dolny przedział ufności (2.5%)", "Górny przedział ufności (97.5%)"),
Klasyczne = c(mean_k, sd_k, se_k, sep_k, dolny_k, gorny_k),
Bootstrap = c(mean_boot, sd_boot, se2, sep2,
dolny, gorny)
)
metrics_table %>%
kbl(col.names = c("Metryka", "Klasyczne", "Bootstrap")) %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"),
position = "center", font_size = 14)| Metryka | Klasyczne | Bootstrap |
|---|---|---|
| Średnia | 808.5321623 | 749.8268042 |
| Odchylenie standardowe | 1064.6701855 | 828.1258688 |
| Błąd standardowy | 6.1465689 | 37.0349147 |
| SEP | 0.0076021 | 0.0493913 |
| Dolny przedział ufności (2.5%) | 796.4848873 | 677.2383714 |
| Górny przedział ufności (97.5%) | 820.5794373 | 822.4152371 |
Jakie są wnioski? Czy różnice są istotne? Jak wielkość resamplingu wpływa na wyniki?
Średnia: Średnia z metody klasycznej (808.53) jest zbliżona do średniej bootstrapowej (783.32), ale nieco wyższa.
Odchylenie standardowe: Klasyczne odchylenie standardowe (1064.67) jest większe niż bootstrapowe (810.36). To wskazuje, że wyniki klasyczne są bardziej zmienne w porównaniu z wynikami bootstrapowymi.
Błąd standardowy (Standard Error): Klasyczny błąd standardowy (6.15) jest znacznie mniejszy niż bootstrapowy (36.24). To sugeruje, że metoda bootstrap lepiej uwzględnia rzeczywistą zmienność w próbach, przez co błąd standardowy jest wyższy, ale bardziej realistyczny.
SEP (Standard Error of Prediction): SEP z obu metod jest zbliżony: 0.0076 (klasyczne) vs 0.0078 (bootstrap). Podobne wartości wskazują, że obie metody oferują porównywalną dokładność predykcji.
Dolny i górny przedział ufności: Klasyczne przedziały ufności są węższe (796.48–820.58) niż bootstrapowe (712.29–854.35).
Wielkość resamplingu poprawia stabilność modelu. Im większa wielkość tym mniejszy błąd standarodwy i SEP co zwiększa predykcję.
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.1845942 6.16323plot(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.3, 820.4 ) (797.0, 820.6 )
## 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)
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.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ kableExtra::group_rows() masks dplyr::group_rows()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors?boot.t.testdane2<- 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) = -194.9064 (30.80116)
## 95 percent bootstrap percentile confidence interval:
## -256.4665 -132.8518
##
## 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 Istotna różnica: Średni dochód na osobę w Pomorskim (844.17) jest istotnie wyższy niż w Podkarpackim (648.21). P-value w oby przypadkach jest mniejsze o 0,05 co wskazuje na istotność. Porównanie metod: Zarówno bootstrap, jak i klasyczny test t wskazują na różnicę, ale bootstrap jest bardziej odporny na potencjalne naruszenia założeń.
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)
tabelka## Gender
## Student Male Female
## No 177 183
## Yes 16 24chisq.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.3103chisq.test(tabelka) # różnice???##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tabelka
## X-squared = 0.87218, df = 1, p-value = 0.3504Dla obu testów p-value jest większe od alfy = 0.05, zatem brak podstaw do odrzucenia hipotezy zerowej. Oba testy wskazują na brak istotnych różnic między oczekiwaniami a obserwacjami.
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 ' ' 1Jedyną zmienną istotną statystycznie i wpływającą na średni bilans na karcie kredytowej jest status studenta, reszta zmiennych nie mają istotnego wpływu.
A 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.9449449 0.9259259 0.0000000 0.8698699W 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
)
library(ggstatsplot)
ggbetweenstats(data=Credit,
y=Balance,
x=Married,
nboot=999 #liczba prób bootstrapowych
)
library(ggstatsplot)
ggbetweenstats(data=Credit,
y=Balance,
x=Student,
nboot=999 #liczba prób bootstrapowych
)
library(ggstatsplot)
ggbetweenstats(data=Credit,
y=Balance,
x=Gender,
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?
Różnice nie są istotne, bez względu na to czy korzystamy z klasycznej ANOVY czy ANOVY z bootstrapem. W dalszym ciągu jedyną istotną zmienną jest status studenta. Bootstrap wykazuję się większą odpornością na naruszenia założeń, a resampling wpływa pozywtywnie na stabilność wyników.