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)
## Warning: pakiet 'foreign' został zbudowany w wersji R 4.3.3
dane <- read.spss("C:/Users/mskox/Downloads/dane.sav", to.data.frame = TRUE)
attach(dane)
dochod<-Dochod_na_osobe
#punktowo klasycznie
mean_1<-mean(dochod)
sd_1<-sd(dochod)
se_1<-se<-sd(dochod)/sqrt(length(dochod))
se_1
## [1] 6.146569
sep_1<-se/mean(dochod)
sep_1
## [1] 0.007602133
#przedziałowo klasycznie
dolny_1 <- mean(dochod)-1.96*se   #dolna krawędź
gorny_1 <- mean(dochod)+1.96*se   #górna krawędź

# 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
#przedziałowo Efron
dolny<-quantile(boot.data,0.975)
gorny<-quantile(boot.data,0.025)
cbind(dolny,gorny)
##          dolny gorny
## 97.5% 1669.535     0

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.

# Tabela porównanie

tabela <- 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_1, sd_1, se_1, sep_1, dolny_1, gorny_1),
  Bootstrap = c(mean_boot, sd_boot, se2, sep2, 
                dolny, gorny)
)

print(tabela)
##                           Metryka    Klasyczne    Bootstrap
## 1                         Średnia 8.085322e+02 6.605794e+02
## 2          Odchylenie standardowe 1.064670e+03 5.445747e+02
## 3                Błąd standardowy 6.146569e+00 7.701450e+01
## 4                             SEP 7.602133e-03 9.304815e-03
## 5  Dolny przedział ufności (2.5%) 7.964849e+02 1.669535e+03
## 6 Górny przedział ufności (97.5%) 8.205794e+02 0.000000e+00

Jakie są wnioski?

Wartości średnich są znacząco różne między metodą klasyczną (8.085) a bootstrapową (9.246), co może sugerować, że bootstrap lepiej odwzorowuje niektóre cechy rozkładu w danych.

Odchylenia standardowe oraz błędy standardowe są znacznie mniejsze w metodzie bootstrapowej niż w klasycznej, co może sugerować większą precyzję estymacji w tej metodzie, choć wymaga to dodatkowej analizy.

Jak wielkość resamplingu wpływa na wyniki?

Większa liczba resamplowań zwiększa stabilność wyników, zmniejsza wariancję oszacowań, ale może prowadzić do nadmiernego dopasowania (overfitting) w przypadku małych próbek. Z kolei mniejsza liczba resamplowań zwiększa losowość w wynikach, ale przyspiesza proces obliczeniowy.

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:

library(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.1472602     6.05451
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.5 )   (797.5, 820.6 )  
## 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)
## ── 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
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) = -194.417 (30.7996) 
## 95 percent bootstrap percentile confidence interval:
##  -250.0798 -135.7153
## 
## 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

Zarówno test t-studenta (p-value znacznie poniżej 0.05), jak i test bootstrapowy wskazują na istotną różnicę w dochodach na osobę pomiędzy województwami pomorskim a podkarpackim. Obie metody oszacowały różnicę na około -197, co pokazuje, że dochody w województwie podkarpackim są znacznie niższe niż w pomorskim.

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.3
data("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.3088
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

W obu testach wartość p-value >0.05, co oznacza brak podstaw do odrzucenia hipotezy zerowej. Test z symulowanym p-value (bootstrapowym) i test z korelacją Yatesa dają podobne wartości p-value (odpowiednio 0.3253 i 0.3504), więc oba te testy prowadzą do tych samych wniosków.

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

Tylko status studenta istotnie wpływa na średni bilans na karcie kredytowej. Pozostałe zmienne (pochodzenie, stan cywilny, płeć) nie mają statystycznie istotnego wpływu.

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.9499499 0.9269269 0.0000000 0.8578579

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.03167
ggbetweenstats(data=Credit,
  y=Balance,
  x=Ethnicity,
  nboot=999  #liczba prób bootstrapowych
)

Czy różnice są istotne?

Różnice nie są istotne zarówni dla klasycznej ANOVA jak i ANOVA z bootstrapem.

Jak wielkość resamplingu wpływa na wyniki?

Większa liczba resamplowań zwiększa stabilność wyników i dokładność przedziałów ufności. W tym przypadku wyniki są stabilne.

Jakie są różnice między testem ANOVA a testem ANOVA z bootstrapem?