Bootstrap

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("D:/Users/Mateusz/Desktop/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

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?

set.seed(123) 
dochod <- rnorm(1000, mean = 800, sd = 1000) 

mean_classic <- mean(dochod)
sd_classic <- sd(dochod)
se_classic <- sd_classic / sqrt(length(dochod))
sep_classic <- se_classic / mean_classic
dolny_classic <- mean_classic - 1.96 * se_classic
gorny_classic <- mean_classic + 1.96 * se_classic

results <- data.frame(
  nobs = integer(),
  mean_boot = numeric(),
  sd_boot = numeric(),
  se2 = numeric(),
  sep2 = numeric(),
  dolny = numeric(),
  gorny = numeric()
)

B <- 999
sample_sizes <- c(50, 250, 500, 1000, 10000) 

for (nobs in sample_sizes) {
  mean_dochod <- rep(0, B)
  
  for (i in 1:B) {
    boot.data <- sample(dochod, nobs, replace = TRUE)
    mean_dochod[i] <- mean(boot.data)
  }
  
  mean_boot <- mean(mean_dochod)
  sd_boot <- sd(mean_dochod)
  se2 <- sd_boot / sqrt(nobs)
  sep2 <- se2 / mean_boot
  dolny <- mean_boot - 1.96 * se2
  gorny <- mean_boot + 1.96 * se2
  
  results <- rbind(results, data.frame(
    nobs = nobs,
    mean_boot = mean_boot,
    sd_boot = sd_boot,
    se2 = se2,
    sep2 = sep2,
    dolny = dolny,
    gorny = gorny
  ))
}

print("Wyniki klasyczne:")

## [1] "Wyniki klasyczne:"

print(data.frame(
  mean_classic = mean_classic,
  sd_classic = sd_classic,
  se_classic = se_classic,
  sep_classic = sep_classic,
  dolny_classic = dolny_classic,
  gorny_classic = gorny_classic
))

##   mean_classic sd_classic se_classic sep_classic dolny_classic gorny_classic
## 1     816.1279    991.695   31.36015  0.03842553       754.662      877.5938

print("Wyniki bootstrapowe dla różnych nobs:")

## [1] "Wyniki bootstrapowe dla różnych nobs:"

print(results)

##    nobs mean_boot    sd_boot         se2         sep2    dolny    gorny
## 1    50  817.8213 142.860264 20.20349231 0.0247040425 778.2225 857.4201
## 2   250  815.4398  62.041303  3.92383652 0.0048119267 807.7491 823.1305
## 3   500  817.2049  43.816204  1.95952021 0.0023978322 813.3642 821.0456
## 4  1000  816.9679  30.381089  0.96073439 0.0011759756 815.0849 818.8510
## 5 10000  816.3110   9.473752  0.09473752 0.0001160557 816.1254 816.4967

library(ggplot2)
ggplot(results, aes(x = as.factor(nobs), y = mean_boot)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = dolny, ymax = gorny), width = 0.2) +
  labs(
    title = "Bootstrapowe średnie i przedziały ufności",
    x = "Liczba próbek (nobs)",
    y = "Średnia bootstrapowa"
  ) +
  theme_minimal()

Wnioski: Porównanie wyników klasycznych i bootstrapowych pokazuje, że obie metody dają zbliżone średnie, co wskazuje na zgodność w estymacji wartości centralnej. Różnice pojawiają się w odchyleniu standardowym i błędach standardowych, które w bootstrapie maleją wraz ze wzrostem liczby próbek, poprawiając precyzję estymacji. Przedziały ufności są szersze przy mniejszych próbach bootstrapowych, ale zbliżają się do klasycznych wyników dla większych rozmiarów próby. Bootstrap jest szczególnie przydatny przy małych próbach lub nienormalnych rozkładach, natomiast dla większych próbek różnice między metodami stają się nieistotne, co świadczy o ich zbieżności.

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* 816.1279 0.7655209    30.39423

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%   (755.8, 874.9 )   (756.9, 881.0 )  
## 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.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ── 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.9369 (30.80696) 
## 95 percent bootstrap percentile confidence interval:
##  -260.0526 -136.8628
## 
## 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: Porównanie wyników testu t-Studenta z testem bootstrapowym dla dochodów na osobę w województwach pomorskim i podkarpackim wykazuje dużą zgodność między metodami. Średni dochód na osobę w Pomorskiem (844.17) jest istotnie wyższy niż w Podkarpackiem (648.21).Obie metody wykazują istotną statystycznie różnicę w dochodach na osobę między województwami (p < 0.001). Przedziały ufności z obu metod są bardzo podobne, co świadczy o ich zgodności. Bootstrap dodatkowo uwzględnia potencjalne odstępstwa od normalności rozkładu, co czyni go bardziej elastycznym.

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

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

Wnioski: Porównując wyniki testów proporcji dla zmiennych “Student” i “Gender”, zauważamy, że wyniki z zastosowaniem różnych metod nie wskazują istotnych różnic pomiędzy grupami. Wyniki wszystkich metod, zarówno z testów chi-kwadrat (symulowany i z korektą Yatesa), jak i bootstrapowych, wskazują, że status studenta nie różni się istotnie w zależności od płci. Włączanie bootstrapu w tym przypadku nie zmienia końcowego wniosku, ale zwiększa pewność analizy w przypadku niestandardowych rozkładów danych.

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

Wnioski: Przeprowadzony test ANOVA dla danych “Credit” (z pakietu ISLR) badał, czy średni bilans na karcie kredytowej różni się istotnie w zależności od pochodzenia, stanu cywilnego, statusu studenta i płci.Jedynie status studenta (Student) znacząco wpływa na średni bilans na karcie kredytowej (p < 0.001). Pozostałe zmienne (pochodzenie, stan cywilny, płeć) nie wykazują istotnego wpływu na tę zmienną.

library(ISLR)
data("Credit")

# Wizualizacja różnic za pomocą boxplotów
library(ggplot2)

# Pochodzenie (Ethnicity)
ggplot(Credit, aes(x = Ethnicity, y = Balance, fill = Ethnicity)) +
  geom_boxplot() +
  labs(title = "Średni bilans w zależności od pochodzenia",
       x = "Pochodzenie",
       y = "Bilans") +
  theme_minimal()

# Stan cywilny (Married)
ggplot(Credit, aes(x = factor(Married), y = Balance, fill = factor(Married))) +
  geom_boxplot() +
  labs(title = "Średni bilans w zależności od stanu cywilnego",
       x = "Stan cywilny",
       y = "Bilans") +
  theme_minimal()

# Status studenta (Student)
ggplot(Credit, aes(x = Student, y = Balance, fill = Student)) +
  geom_boxplot() +
  labs(title = "Średni bilans w zależności od statusu studenta",
       x = "Status studenta",
       y = "Bilans") +
  theme_minimal()

# Płeć (Gender)
ggplot(Credit, aes(x = Gender, y = Balance, fill = Gender)) +
  geom_boxplot() +
  labs(title = "Średni bilans w zależności od płci",
       x = "Płeć",
       y = "Bilans") +
  theme_minimal()

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.9639640 0.9349349 0.0000000 0.8688689

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?

Wnioski z testu ANOVA i ANOVA z bootstrapem:

Wyniki testu ANOVA bez bootstrapu (\(F = 0.0463\), \(p = 0.9548\)) oraz z bootstrapem (\(p = 0.954955\)) wskazują na brak istotnych różnic w średnim bilansie na karcie kredytowej między grupami etnicznymi. Włączenie bootstrapu, który uwzględnia potencjalne odstępstwa od normalności danych i zmienności wariancji między grupami, nie wpłynęło znacząco na wyniki, co wskazuje na stabilność metody. Bootstrap pozwala na bardziej elastyczne oszacowanie p-wartości, szczególnie w przypadku małych prób lub danych o niestandardowym rozkładzie, jednak w tej analizie, gdzie próbki są duże i wariancje równe, wyniki klasycznej ANOVA i ANOVA z bootstrapem są niemal identyczne. Wnioski końcowe potwierdzają, że średni bilans nie różni się istotnie między grupami etnicznymi, a choć bootstrap jest bardziej uniwersalny, w tym przypadku jego zastosowanie nie zmienia ogólnego wniosku.