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("/Users/Agata/Downloads/dane.sav", to.data.frame = TRUE)
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=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
mean_boot
## [1] 806.0708
sd_boot
## [1] 1039.627
se2
## [1] 10.39627
sep2
## [1] 0.007625346
dolny
## [1] 785.6941
gorny
## [1] 826.4475

Wyniki dla różnych wielkości próby:

Dla wielkości próby ‘r nobs=50’ otrzymujemy: - średnią próbkową: ‘r mean_boot’- 629.6218 - odchylenie standardowe: ‘r sd_boot’ - 644.0206 - błąd standardowy: ‘r se2’ - 91.07826 - procentowy błąd standardowy: ‘r sep2’ - 0.00976232 - dolną krawędź przedziału ufności: ‘r dolny’ - 451.1084 - górną krawędź przedziału ufności: ‘r gorny’ - 808.1352

Dla wielkości próby ‘r nobs=250’ otrzymujemy: - średnią próbkową: ‘r mean_boot’- 692.0217 - odchylenie standardowe: ‘r sd_boot’ - 733.7338 - błąd standardowy: ‘r se2’ - 46.4054 - procentowy błąd standardowy: ‘r sep2’ - 0.008882046 - dolną krawędź przedziału ufności: ‘r dolny’ - 601.0671 - górną krawędź przedziału ufności: ‘r gorny’ - 782.9763

Dla wielkości próby ‘r nobs=500’ otrzymujemy: - średnią próbkową: ‘r mean_boot’- 815.6211 - odchylenie standardowe: ‘r sd_boot’ - 948.8071 - błąd standardowy: ‘r se2’ - 42.43194 - procentowy błąd standardowy: ‘r sep2’ - 0.007536059 - dolną krawędź przedziału ufności: ‘r dolny’ - 732.4545 - górną krawędź przedziału ufności: ‘r gorny’ - 898.7877

Dla wielkości próby ‘r nobs=1000’ otrzymujemy: - średnią próbkową: ‘r mean_boot’- 823.6502 - odchylenie standardowe: ‘r sd_boot’ - 1144.474 - błąd standardowy: ‘r se2’ - 36.19145 - procentowy błąd standardowy: ‘r sep2’ - 0.007462596 - dolną krawędź przedziału ufności: ‘r dolny’ - 752.715 - górną krawędź przedziału ufności: ‘r gorny’ - 894.5854

Dla wielkości próby ‘r nobs=10000’ otrzymujemy: - średnią próbkową: ‘r mean_boot’- 799.2741 - odchylenie standardowe: ‘r sd_boot’ - 1088.278 - błąd standardowy: ‘r se2’ - 10.88278 - procentowy błąd standardowy: ‘r sep2’ - 0.007690189 - dolną krawędź przedziału ufności: ‘r dolny’ - 777.9439 - górną krawędź przedziału ufności: ‘r gorny’ - 820.6044

Punktowo klasycznie: 808.5322 1064.67 6.146569 0.007602133 796.4849 820.5794

Porównaj wyniki z obu metod.

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

Małe próby (np. 50) charakteryzują się dużymi błędami standardowymi i szerokimi przedziałami ufności. Wyniki są mniej precyzyjne i mogą być bardziej podatne na szum w danych. Duże próby (np. 10 000) prowadzą do bardziej stabilnych i precyzyjnych wyników, co przejawia się w mniejszym błędzie standardowym i węższych przedziałach ufności.

Różnice między wynikami dla małych prób a metodą punktową klasyczną są zauważalne, ale stają się coraz mniej istotne w miarę zwiększania rozmiaru próby. Dla dużych prób (np. 1000, 10 000) wyniki z resamplingu są zgodne z wynikami klasycznymi.

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.02432293     6.08367
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.6, 820.4 )   (797.5, 821.2 )  
## 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) 
library(tidyverse)
## ── 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.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.8053 (30.79507) 
## 95 percent bootstrap percentile confidence interval:
##  -257.9747 -132.3578
## 
## 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

Dochody na osobę istotnie różnią się między województwami pomorskim a podkarpackim. Zarówno test t-Studenta, jak i test bootstrapowy wykazują, że średni dochód w woj. podkarpackim jest niższy. Obie metody prowadzą do podobnych wniosków, co świadczy o wiarygodności analizy. Bootstrap dodatkowo potwierdza wyniki, uwzględniając potencjalne naruszenia założeń testu t.

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.3373
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

Nie ma istotnych różnic w rozkładzie statusu studenta między płciami. Wartości p w obu testach (0.3293 i 0.3504) są większe od poziomu istotności 0.05. Jeśli próbka jest mała lub liczności w tabeli są niskie, bardziej wiarygodne są wyniki testu z symulacją Monte Carlo.

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

Jedynie status studenta wpływa istotnie na średni bilans na karcie kredytowej. Studenci mają wyraźnie różne średnie bilanse w porównaniu do osób niebędących studentami. Wykres potwierdza, że status studenta jest kluczowym czynnikiem różnicującym średni bilans na karcie kredytowej. Średni bilans jest wyraźnie niższy dla studentów w porównaniu z osobami niebędącymi studentami. Pochodzenie, stan cywilny i płeć nie mają istotnego wpływu na tę zmienną.

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.9659660 0.9499499 0.0000000 0.8768769

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
)

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?

Brak istotnych różnic w średnim saldzie karty kredytowej (Balance) między grupami etnicznymi. Bootstrap daje dodatkową stabilność wyników, ale w tym przypadku jego wpływ jest niewielki. Klasyczny test ANOVA i ANOVA z bootstrapem są zgodne, co oznacza, że klasyczne założenia zostały spełnione.