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/Ewa/Desktop/python/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] 812.9525
sd_boot
## [1] 979.0649
se2
## [1] 9.790649
sep2
## [1] 0.007560797
dolny
## [1] 793.7629
gorny
## [1] 832.1422

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

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, takie jak próba o rozmiarze 50, charakteryzują się większymi błędami standardowymi oraz szerszymi przedziałami ufności, co oznacza, że wyniki są mniej precyzyjne. W takich przypadkach wyniki mogą być bardziej podatne na wpływ szumu w danych, co może prowadzić do trudności w dokładnej interpretacji i wyciąganiu wiarygodnych wniosków. Mała liczba obserwacji w próbie sprawia, że każdy pojedynczy pomiar ma większy wpływ na ostateczny wynik, co zwiększa niepewność całej analizy.

Z kolei duże próby, takie jak próbki liczące 10 000 jednostek, prowadzą do uzyskania bardziej stabilnych i precyzyjnych wyników. W takich przypadkach błąd standardowy maleje, a przedziały ufności stają się węższe, co oznacza, że oszacowania są bardziej wiarygodne i dokładne. Większa liczba obserwacji pozwala na lepsze uchwycenie prawdziwej struktury danych, minimalizując wpływ ewentualnych zakłóceń, a wyniki stają się bardziej reprezentatywne.

Różnice między wynikami uzyskanymi na małych próbach a wynikami uzyskanymi za pomocą klasycznych metod punktowych są widoczne, ale w miarę jak rozmiar próby rośnie, te różnice stają się coraz mniej istotne. Dla większych prób, takich jak próby o rozmiarze 1000 czy 10 000, wyniki uzyskane przy użyciu metod resamplingu są bardzo zbliżone do tych uzyskanych klasycznymi technikami. To wskazuje na to, że przy większych próbach metody resamplingu i klasyczne techniki dają spójne wyniki, a różnice między nimi stają się minimalne. Taki trend podkreśla znaczenie rozmiaru próby w zapewnianiu precyzyjnych i stabilnych oszacowań w analizach statystycznych.

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.08959405    6.233259
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.4, 820.8 )   (797.1, 820.9 )  
## 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.

#install.packages("MKinfer")
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
?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.9394 (30.77902) 
## 95 percent bootstrap percentile confidence interval:
##  -256.8053 -128.9856
## 
## 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ę wykazują istotne różnice między województwami pomorskim i podkarpackim. Zarówno test t-Studenta, jak i metoda bootstrapowa wskazują, że średni dochód w województwie podkarpackim jest niższy niż w województwie pomorskim. Obie techniki prowadzą do zbliżonych wniosków, co zwiększa pewność co do rzetelności wyników analizy. Dodatkowo, zastosowanie bootstrapu wzmacnia te wyniki, uwzględniając ewentualne naruszenia założeń testu t-Studenta, co czyni całą analizę bardziej odporną na potencjalne błędy związane z tymi założeniami.

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
#install.packages("ISLR")
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.3333
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 zaobserwowano istotnych różnic w rozkładzie statusu studenta między płciami. Wartości p uzyskane w obu testach, wynoszące odpowiednio 0.3293 i 0.3504, są wyższe od przyjętego poziomu istotności 0.05, co oznacza, że nie ma wystarczających dowodów na to, aby odrzucić hipotezę zerową o braku różnic. Oznacza to, że status studenta rozkłada się w podobny sposób niezależnie od płci. Warto jednak zaznaczyć, że w przypadku małych prób lub sytuacji, gdy liczności w tabelach są niskie, wyniki testu mogą być mniej wiarygodne. W takich przypadkach, bardziej wiarygodne i solidne mogą być wyniki uzyskane z testu z symulacjami Monte Carlo, który nie opiera się na klasycznych założeniach rozkładu i może lepiej radzić sobie z ograniczeniami związanymi z małymi próbami i niskimi częstościami w tabelach.

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

Status studenta jest jedynym czynnikiem, który ma istotny wpływ na średni bilans na karcie kredytowej. Różnice te są wyraźne, ponieważ studenci mają znacznie niższy średni bilans w porównaniu do osób, które nie są studentami. Wykres potwierdza, że status studenta odgrywa kluczową rolę w różnicowaniu średniego salda na karcie kredytowej, co sugeruje, że bycie studentem wiąże się z innym wzorcem korzystania z kredytu. Z kolei inne zmienne, takie jak pochodzenie, stan cywilny czy płeć, nie wykazują istotnego wpływu na wysokość średniego bilansu na karcie kredytowej. W kontekście tych wyników można stwierdzić, że to właśnie status studenta jest głównym czynnikiem, który wpływa na saldo na karcie kredytowej.

A teraz z włączonym bootstrapem:

#install.packages("lmboot")
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.9549550 0.9319319 0.0000000 0.8638639
#install.packages("ggstatsplot")
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?

Nie zaobserwowano istotnych różnic w średnim saldzie na karcie kredytowej (Balance) między grupami etnicznymi, co sugeruje, że pochodzenie etniczne nie ma znaczącego wpływu na tę zmienną. Przeprowadzenie analizy przy użyciu metody bootstrap pozwoliło na dodatkowe wzmocnienie wyników, zwiększając ich stabilność, choć w tym przypadku wpływ bootstrapu na wyniki był stosunkowo niewielki. Warto jednak zaznaczyć, że zarówno klasyczny test ANOVA, jak i jego wersja z bootstrapem dały zgodne wyniki, co świadczy o tym, że założenia klasycznego testu ANOVA zostały spełnione. Oznacza to, że wyniki uzyskane z klasycznej analizy statystycznej można uznać za wiarygodne, a stosowanie bootstrapu w tym przypadku nie wprowadziło zasadniczych zmian w interpretacji danych.