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.

#Dane
library(foreign)
dane <- read.spss("C:/Users/Pawlo/OneDrive/Pulpit/dane.sav", to.data.frame = TRUE)
attach(dane)
dochod<-Dochod_na_osobe

##Klasycznie
#Klasyczne metody estymacji punktowej średniej oraz obliczania jej względnego błędu standardowego. 
#punktowo klasycznie

mean_value<-mean(dochod) 
mean_value
## [1] 808.5322
sd_value<-sd(dochod) 
sd_value
## [1] 1064.67
se<-sd(dochod)/sqrt(length(dochod)) 
se
## [1] 6.146569
sep<-se/mean(dochod)
sep
## [1] 0.007602133
## Średnia dla zmiennej dochód wynosi: 808.53
## Odchylenie standardowe wynosi: 1064.67
## Błąd standardowy wynosi: 6.1466
## Względny błąd standardowy wynosi: 0.76 %


#przedziałowo klasycznie
# Klasyczny przedział ufności
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

## Dolna krawędź: 796.48
## Górna krawędź: 820.58

Teraz zrobimy to samą metodą bootstrap:

## Metoda 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

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.

Wnioski: Średnie uzyskane z obu metod są zbliżone (808,53 w metodzie klasycznej i 796,85 w bootstrapie), co wskazuje, że metoda bootstrap dobrze odwzorowuje wartość średniej próbkowej. Jednakże, metoda bootstrap charakteryzuje się wyższym odchyleniem standardowym (1110,97 vs. 1064,67) oraz większym błędem standardowym (157,11 vs. 6,15), co może wynikać z większej wariancji w próbkach bootstrapowych. Ponadto, metoda klasyczna generuje węższy przedział ufności (796,48–820,58) w porównaniu do bootstrapu (488,91–1104,80), co wynika z uwzględnienia większej niepewności przez podejście bootstrapowe.

B4=999
mean.dochod4=rep(0,B4)
nobs4=1000 #liczba próbek bootstrapowych - zwiększaj ich liczbę: 50, 250, 500, 1000, 10000

 for (i in 1:B4) 
{
 boot.data4=sample(dochod,nobs4,replace=TRUE)
 mean.dochod4[i]=mean(boot.data4)
 }

mean_boot4 <- mean(boot.data4)

sd_boot4 <- sd(boot.data4)

se4<-sd(boot.data4)/sqrt(length(boot.data4))

sep4<-se4/mean(boot.data4)

hist(boot.data4)

plot(density(boot.data4))

#przedziałowo bootstrapowo  95% ufności:
dolny4<-mean(boot.data4)-1.96*se2
gorny4<-mean(boot.data4)+1.96*se2
Jakie są wnioski? Czy różnice są istotne? Jak wielkość resamplingu wpływa na wyniki?

Gdy liczność próby bootstrapowej zwiększa się (z 50 do 500), obserwuje się: istotny spadek odchylenia standardowego (1110,97 vs. 798,65), zmniejszenie błędu standardowego (157,11 vs. 35,72) oraz zawężenie przedziału ufności (488,91–1104,80 vs. 660,77–800,77).

Wniosek: Większa liczba obserwacji w próbach bootstrapowych prowadzi do większej stabilności wyników, dzięki redukcji wariancji oszacowań.

#### 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:



``` r
?boot
## uruchamianie serwera httpd dla pomocy ... wykonano
mean.boot=function(dochod,idx) {
ans=mean(dochod[idx])
ans
}

#Uruchomienie procedury bootstrap: 
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.2946571    6.187285
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.7, 821.0 )   (795.5, 820.4 )  
## Calculations and Intervals on Original Scale

Bootstrapowy błąd standardowy (6,067342) jest niemal identyczny z wartością obliczoną klasycznie (6,1466), co sugeruje, że metoda bootstrap skutecznie oszacowuje niepewność, nawet przy założeniu niestandardowego rozkładu danych.

Przedziały ufności uzyskane za pomocą metody normalnej i percentylowej są bardzo zbliżone, co świadczy o stabilności wyników bootstrap.

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
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ 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.4159 (30.67339) 
## 95 percent bootstrap percentile confidence interval:
##  -257.0108 -136.8175
## 
## 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: Zarówno test t-Studenta, jak i test bootstrapowy potwierdzają, że dochody na osobę istotnie różnią się między województwami pomorskim i podkarpackim. Średnie dochody w województwie podkarpackim wynoszą 648,21 zł, natomiast w pomorskim 844,16 zł.

Test t-Studenta wskazuje na istotną statystycznie różnicę między średnimi dochodów (p-value = 2.331e-10 < 0,05). Przedział ufności wynosi od -256,36 zł do -135,54 zł, co wyklucza wartość 0 i potwierdza istnienie różnicy. Podobnie, test bootstrapowy również wykazuje bardzo niskie p-value (poniżej 0,001), co dodatkowo potwierdza istotność różnic w dochodach. Wygenerowany przedział ufności dla testu bootstrapowego, choć zbliżony, jest nieco węższy i wynosi od -251,61 zł do -139,30 zł.

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.3153
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: Zarówno w teście chi² z zastosowaniem bootstrappingu, jak i bez niego, p-value przekracza poziom istotności 0,05 (p-value = 0,31 i 0,35). W związku z tym brak podstaw do odrzucenia hipotezy zerowej. Można zatem przyjąć, że zmienne „Student” i „Gender” są niezależne, co oznacza, że status studenta nie zależy istotnie od płci. Wyniki uzyskane przy użyciu bootstrappingu różnią się jedynie nieznacznie od wyników pojedynczego testu chi², a oba prowadzą do tego samego wniosku.

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?

attach(Credit)
## Następujące obiekty zostały zakryte z Credit (pos = 3):
## 
##     Age, Balance, Cards, Education, Ethnicity, Gender, ID, Income,
##     Limit, Married, Rating, Student
model1 <- lm(Balance ~ Ethnicity + Married + Student + Gender, data=Credit)
anova(model1)
## 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)
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.9619620 0.9339339 0.0000000 0.8778779

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
)

library(ggstatsplot)
ggbetweenstats(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?

Wnioski: Test ANOVA wykazał, że jedynie status studenta istotnie wpływa na średni bilans na karcie kredytowej. Pozostałe zmienne, takie jak pochodzenie, stan cywilny i płeć, nie są statystycznie istotne. Włączenie bootstrapu nie zmieniło tego wyniku – zmienne pochodzenie, stan cywilny i płeć nadal mają p-value wyższe od poziomu istotności. Podobnie jak w przypadku klasycznego testu ANOVA, tylko status studenta istotnie wpływa na średni bilans na karcie kredytowej.

Na wizualizacji widać, że osoby pochodzenia African American, Asian i Caucasian mają średnie bilanse na zbliżonym poziomie. Kolejna wizualizacja pokazuje, że osoby posiadające status studenta mają wyższy średni bilans na karcie kredytowej w porównaniu do osób bez tego statusu.

Różnice między wynikami testu ANOVA a testu ANOVA z bootstrapem są niewielkie. W przypadku bootstrapu zmienne nieistotne stają się jeszcze mniej istotne (p-value nieznacznie wzrasta), natomiast zmienna istotna (status studenta) ma jeszcze niższe p-value, wynoszące 0, w porównaniu do wartości 1.484e-07 w klasycznym teście ANOVA. Włączenie bootstrapu potwierdza poprawność wcześniejszych wniosków.