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("C:/Users/weron/Downloads/dane.sav")
attach(dane)
dochod<-Dochod_na_osobe

#punktowo klasycznie
mean(dochod)
## [1] 808.5322
sd(dochod)
## [1] 1064.67
se<-sd(dochod)/sqrt(length(dochod))
se
## [1] 6.146569
sep<-se/mean(dochod)
sep
## [1] 0.007602133
#przedziałowo klasycznie
mean(dochod)-1.96*se   #dolna krawędź
## [1] 796.4849
mean(dochod)+1.96*se   #górna krawędź
## [1] 820.5794
# 1.96 to kwantyl rozkładu normalnego dla 95% ufności

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)

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
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_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
B=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
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)

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 różnych wielkości próby:

Dla wielkości próby ‘r nobs=50’ otrzymujemy:

- średnią próbkową: ‘r mean_boot’

- odchylenie standardowe: ‘r sd_boot’

- błąd standardowy: ‘r se2’

- procentowy błąd standardowy: ‘r sep2’

- dolną krawędź przedziału ufności: ‘r dolny’

- górną krawędź przedziału ufności: ‘r gorny’

Porównaj wyniki z obu metod.

Jakie są wnioski? Czy różnice są istotne? Jak wielkość resamplingu wpływa na wyniki?

Wnioski:

Klasyczna średnia jest stabilna i nie zależy od dodatkowych losowych próbek. Średnia bootstrapowa zmienia się wraz ze wzrostem liczby próbek (nobs). Przy małych próbkach (nobs = 50) średnia jest zawyżona, ale dla większych prób (nobs = 10 000) konwerguje do wartości klasycznej. Klasyczna metoda daje niski błąd standardowy, ponieważ zakłada określoną strukturę danych (np. normalność). W bootstrapie SE jest większy, szczególnie dla małych próbek, ale maleje wraz ze wzrostem liczby próbek. Klasyczna metoda daje węższy przedział ufności, ponieważ zakłada mniejszą zmienność. Bootstrap przy małej liczbie próbek daje bardzo szerokie przedziały, które stają się coraz węższe i bardziej zbliżone do klasycznych wraz ze wzrostem liczby próbek.

Przy małej liczbie próbek bootstrapowej (nobs = 50), różnice między wynikami obu metod są znaczne. Bootstrap daje zawyżone średnie i szerokie przedziały ufności, co wynika z dużej zmienności losowego resamplingu. Dla większych prób bootstrapowych (nobs ≥ 1000), wyniki są coraz bardziej zbieżne z metodą klasyczną, a różnice stają się marginalne.

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 ... wykonano
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.0196568    5.990044
plot(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.3 )   (797.5, 820.5 )  
## Calculations and Intervals on Original Scale

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(MKinfer) 
## Warning: pakiet 'MKinfer' został zbudowany w wersji R 4.3.3
library(tidyverse)
## Warning: pakiet 'tidyverse' został zbudowany w wersji R 4.3.3
## Warning: pakiet 'ggplot2' został zbudowany w wersji R 4.3.3
## Warning: pakiet 'dplyr' został zbudowany w wersji R 4.3.2
## Warning: pakiet 'stringr' został zbudowany w wersji R 4.3.2
## Warning: pakiet 'lubridate' został zbudowany w wersji R 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ 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

# trzeba zmienić Wojewodztwo na tekst
dane$Wojewodztwo <- as.character(as_factor(dane$Wojewodztwo))

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) = -196.2239 (30.74559) 
## 95 percent bootstrap percentile confidence interval:
##  -253.3041 -139.8224
## 
## 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
library(dplyr)
# Filtrowanie danych
dane2 <- dane %>%
  filter(Wojewodztwo %in% c("Pomorskie", "Podkarpackie"))

# Test bootstrap
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) = -195.1811 (30.75426) 
## 95 percent bootstrap percentile confidence interval:
##  -254.5761 -139.9989
## 
## 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

Obie metody wskazują na istotne statystycznie różnice między średnimi dochodami w województwach pomorskim i podkarpackim. Klasyczny test t-Studenta i test bootstrapowy dają spójne wyniki, potwierdzając istotność różnic.

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.2
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.3103
chisq.test(tabelka) # różnice???
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tabelka
## X-squared = 0.87218, df = 1, p-value = 0.3504

Zarówno test klasyczny, jak i test z symulacją wskazują na brak istotnej różnicy w proporcjach statusu studenta (YES/NO) między płciami (Male/Female). P-value są większe niż 0.05, co oznacza brak podstaw do odrzucenia hipotezy zerowej.

Wartości statystyki X2 są nieco różne.

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 ' ' 1

A teraz z włączonym bootstrapem:

library(lmboot)
## Warning: pakiet 'lmboot' został zbudowany w wersji R 4.3.3
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.9369369 0.9229229 0.0000000 0.8658659

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.03167
ggbetweenstats(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?

Czynnik Student:

* Istotny wpływ na średni bilans karty kredytowej.

* Studenci mają wyraźnie różny bilans w porównaniu do osób, które nie są studentami.

Pozostałe czynniki (Ethnicity, Married, Gender):

* Brak istotnego wpływu na bilans karty kredytowej.

Klasyczny test ANOVA i bootstrapowy ANOVA dają zgodne wyniki. Bootstrap zapewnia większą elastyczność w analizie danych z odstępstwami od założeń klasycznego testu. Duża liczba replikacji (B=999) stabilizuje wyniki testu bootstrapowego, czyniąc go wiarygodnym narzędziem analizy.