Bootstrap
Zuzia i Maja
14-01-2025
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.
Metoda klasyczna
#punktowo klasycznie
m1 <- mean(dochod)
m1## [1] 808.5322sd1 <- sd(dochod)
sd1## [1] 1064.67se<-sd(dochod)/sqrt(length(dochod))
se## [1] 6.146569sep<-se/mean(dochod)
sep## [1] 0.007602133#przedziałowo klasycznie
d1 <- mean(dochod)-1.96*se #dolna krawędź
d1## [1] 796.4849g1<- mean(dochod)+1.96*se #górna krawędź
g1## [1] 820.5794# 1.96 to kwantyl rozkładu normalnego dla 95% ufnościBootstrap
50 próbek
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*se2250 próbek
B=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_boot2 <- mean(boot.data)
sd_boot2 <- sd(boot.data)
se3<-sd(boot.data)/sqrt(length(boot.data))
sep3<-se/mean(boot.data)
hist(boot.data)
plot(density(boot.data))
#przedziałowo bootstrapowo 95% ufności:
dolny2<-mean(boot.data)-1.96*se2
gorny2<-mean(boot.data)+1.96*se2Porównanie
Porównaj wyniki z obu metod.
Jakie są wnioski? Czy różnice są istotne? Jak wielkość resamplingu wpływa na wyniki?
| Metoda klasyczna | Bootstrap 50 | Bootstrap 250 | |
|---|---|---|---|
| Średnia | 808.5321623 | 800.3602655 | 847.7913357 |
| Odchylenie | 1064.6701855 | 908.4800986 | 989.7019318 |
| Błąd standardowy | 6.1465689 | 128.4784877 | 62.5942462 |
| Błąd procentowy | 0.0076021 | 0.0076798 | 0.0072501 |
| Dolna krawędź | 796.4848873 | 548.5424297 | 595.9734999 |
| Górna krawędź | 820.5794373 | 1052.1781013 | 1099.6091715 |
Porównując wyniki metody klasycznej i bootstrap, zauważamy, że średnia dla metody klasycznej jest wyższa niż dla bootstrapu, niezależnie od liczby resamplingu. Odchylenie standardowe w metodzie bootstrap są mniejsze niż w metodzie klasycznej, jednak wraz ze wzrostem liczby próbek odchylenie standardowe wyraźnie wzrasta. Błędy standardowe i błędy procentowe w metodzie klasycznej są znacznie niższe, co może wynikać z dużej zmienności w próbach bootstrapowych. Liczba resamplingu ma wyraźny wpływ na wyniki – przy większej liczbie prób średnie wartości stabilizują się, co wskazuje na większą dokładność, ale odchylenie standardowe rośnie, co pokazuje większą rozpiętość wyników. Dolna i górna krawędź przedziału ufności również się różnią, co świadczy o zmienności wyników przy różnych próbach. Różnice w średnich między metodą klasyczną a bootstrapem mogą być statystycznie istotne, szczególnie jeśli bootstrap dokładniej reprezentuje populację, dzięki mniejszym błędom standardowym. Wnioski wskazują, że bootstrap oferuje bardziej elastyczny sposób oszacowania statystyk, szczególnie w małych próbkach, podczas gdy metoda klasyczna jest mniej podatna na zmienność. Liczba prób resamplingu odgrywa kluczową rolę – większa liczba prób zapewni bardziej stabilne oszacowania.
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.1314013 6.177947plot(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.5 ) (795.8, 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)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ 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()
## ✖ dplyr::group_rows() masks kableExtra::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.test
dane2<- dane %>%
filter(Wojewodztwo %in% c("Pomorskie", "Podkarpackie"))
# Bootstrapowy test t-Studenta
bootttest <- boot.t.test(Dochod_na_osobe~Wojewodztwo, R=999, data=dane2)| bootttest.estimate | bootttest.boot.p.value | |
|---|---|---|
| mean in group Podkarpackie | 648.2102 | 0.001001 |
| mean in group Pomorskie | 844.1663 | 0.001001 |
Dla obu testów wartość p-value jest mniejsza od ustalonego poziomu istotności - alfa = 0.05 - zatem mamy podstawy do odrzucenia hipotezy zerowej na rzecz alternatywnej. Możemy powiedzieć, że średnie dochody w województwie pomorskim i podkarpackim różnia się istotnie.
ggplot(dane2, aes(x = Wojewodztwo, y = Dochod_na_osobe)) +
geom_boxplot() +
labs(title = "Porównanie dochodów na osobę w województwach",
x = "Województwo",
y = "Dochód na osobę")
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.3123chisq.test(tabelka) ##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tabelka
## X-squared = 0.87218, df = 1, p-value = 0.3504Klasyczny test chi-kwadrat porównuje rozkłady dwóch zmiennych jakościowych (dla niego próbki muszą być wystarczająco duże). Test chi-kwadrat z bootsrapingiem włącza symulację i generuje rozkład statytskiki testowej za pomocą losowych permutacji danych (bardziej elastyczne podejście).
Porównaniu obu metod podlega wartość p-value dla testu proporcji. W obu przypadkach p-value okazało się większe od założonego poziomu istotności (alfa = 0.05), zatem nie możemy odrzucić hipotezy zerowej o braku wystąpienia istotnych zależności między stutusem studenta a płcią. Wyniki dla obu metod są podobne, więc zakładamy, że założenia klasycznegogo testu zostały spęłnione i dobrze przedstawia zalezności pomiędzy analizowanymi zmiennymi.
Testy ANOVAy
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:
## Warning in ANOVA.boot(Balance ~ Ethnicity + Married + Student + Gender, : This function has only been fully tested for one-way and two-way ANOVA.| zm.objasniajace | p.value |
|---|---|
| Ethnicity | 0.9529530 |
| Married | 0.9389389 |
| Student | 0.0000000 |
| Gender | 0.8618619 |
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)## 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=Student,
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?
Dla klasycznej ANOVY wartość p.value jest statystycznie istotna dla zmiennej student, zatem powiemy że posiadanie statusu studenta wpływa na średni bilans na karcie kredytowej. Takie same wnioski można wyciagnąć z analizy ANOVY z boostrapem co zapewnia stabliność wyników. Wielkość resamplingu znacząco wpływa na wyniki analizy statystycznej. Większa liczba obserwacji prowadzi do większej precyzji estymacji średnich, węższych przedziałów ufności i większej mocy testów, co ułatwia wykrywanie istotnych różnic między grupami. W przypadku mniejszych prób oszacowania są bardziej zmienne, przedziały ufności szersze, a wyniki mniej precyzyjne, co może utrudniać interpretację.Optymalizacja wielkości próby jest kluczowa, aby uzyskać wiarygodne i stabilne wnioski.