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

sample()- Służy do losowania próbek z danych. Można określić np. wielkość próby (size) oraz to, czy losowanie ma być z powtórzeniami (replace = TRUE).

rnorm()- Generuje losowe liczby z rozkładu normalnego. Można określić liczbę obserwacji, średnią (mean) i odchylenie standardowe (sd).

replicate()- Powtarza określoną operację wielokrotnie.

boot() - Wykonuje resampling bootstrapowy i generuje oszacowania bootstrapowe.

boot.ci() - Oblicza przedziały ufności na podstawie próbek bootstrapowych. Obsługuje różne rodzaje przedziałów ufności, np. percentylowe, skorygowane o błąd skośności.

Zadania

Zadanie 1. Estymacja błędu standardowego średniej próbkowej metodą bootstrap.

library(foreign)
library(haven)
dane <- read_sav("C:/Users/MediaMarkt/Downloads/dane.sav")
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

Średni dochód na osobę w badanej próbie wynosi około 808.53zł. Odchylenie standardowe dochodu wynosi około 1064.67zł. Wskazuje ono na dużą zmienność w dochodach między osobami w badanej próbie. Im wyższe odchylenie standardowe, tym większe zróżnicowanie dochodów.Wynik 6.15 sugeruje, że nasza oszacowana średnia (808.53) jest precyzyjna, ponieważ błąd standardowy jest stosunkowo mały w porównaniu do wartości średniej.Względny błąd standardowy wynosi około 0.76%, co wskazuje, że błąd standardowy średniej jest bardzo mały w stosunku do samej średniej.Przedział ufności (95%) daje nam dodatkowe potwierdzenie: prawdziwa średnia populacyjna najprawdopodobniej znajduje się między 796.48 a 820.58zł.

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

Powyższy kod wykorzystuje metodę bootstrap, która polega na losowaniu z powtórzeniami próbek z oryginalnego zestawu danych i wielokrotnym wyliczaniu statystyki (średniej) na tych próbach.Wykonujemy 999 prób bootstrapowych (B = 999). Każda próba losuje 50 obserwacji z powtórzeniami (nobs = 50).Obliczamy średnią dla każdej z tych prób i zapisujemy wyniki w tablicy mean.dochod.

Co się stanie, jak zwiększymy liczbę próbek bootstrapowych?

Gdy zwiększymy liczbę prób bootstrapowych (B), oszacowania statystyk stają się bardziej precyzyjne i stabilne. Średnia bootstrapowa zbiega do rzeczywistej średniej populacyjnej, a błąd standardowy jest dokładniej oszacowany, co prowadzi do bardziej wiarygodnych przedziałów ufności. Rozkład średnich bootstrapowych staje się coraz bardziej regularny i gładki, co pozwala na lepsze odwzorowanie rzeczywistego rozkładu danych.Przy małej liczbie prób wyniki mogą być niestabilne i obarczone większą wariancją, co skutkuje szerszym lub nieprecyzyjnym przedziałem ufności. Wraz ze wzrostem liczby prób histogram bootstrapowych średnich lepiej reprezentuje dane, a przedziały ufności są mniej podatne na przypadkowość. Jednakże, zwiększenie liczby prób prowadzi do wydłużenia czasu obliczeń.

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

Dla wielkości próby 50 otrzymujemy: - średnią próbkową: 820.53 - odchylenie standardowe: 826.07 - błąd standardowy: 116.82 - procentowy błąd standardowy: 0.0074 - dolną krawędź przedziału ufności: 591.56 - górną krawędź przedziału ufności: 1049.51

mean_boot

## [1] 706.3727

sd_boot

## [1] 732.1535

se2

## [1] 103.5421

sep2

## [1] 0.008701595

dolny

## [1] 503.4301

gorny

## [1] 909.3152

Porównaj wyniki z obu metod.

porownanie <- data.frame(
  Metoda = c("Klasyczna", "Bootstrapowa"),
  Srednia = c(808.53, 820.53),
  Odchylenie_Standardowe = c(1064.67, 826.07),
  Blad_Standardowy = c(6.15, 116.82),
  Procentowy_Blad_Standardowy = c(0.0076, 0.0074),
  Dolna_Granica_PU = c(796.48, 591.56),
  Gorna_Granica_PU = c(820.58, 1049.51)
)

print(porownanie)

##         Metoda Srednia Odchylenie_Standardowe Blad_Standardowy
## 1    Klasyczna  808.53                1064.67             6.15
## 2 Bootstrapowa  820.53                 826.07           116.82
##   Procentowy_Blad_Standardowy Dolna_Granica_PU Gorna_Granica_PU
## 1                      0.0076           796.48           820.58
## 2                      0.0074           591.56          1049.51

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

Największa różnica pojawia się w oszacowaniu błędu standardowego. W metodzie klasycznej wynosi on 6.15, podczas gdy metoda bootstrapowa daje znacznie wyższy wynik – 116.82. Wzrost błędu standardowego w bootstrapie może wynikać z losowego charakteru próbek bootstrapowych, szczególnie przy stosunkowo niewielkiej liczbie danych i rozkładzie dochodu, który może być nieregularny. Procentowy błąd standardowy jest jednak porównywalny między metodami – 0.76% w klasycznej i 0.74% w bootstrapowej – co wskazuje, że bootstrap wciąż zapewnia dobrą jakość oszacowań.Porównanie przedziałów ufności (95%) ujawnia znaczącą różnicę: klasyczna metoda wyznacza wąski przedział od 796.48 do 820.58, podczas gdy bootstrapowy przedział jest znacznie szerszy – od 591.56 do 1049.51.Średnia dochodu uzyskana klasycznie wynosi 808.53, podczas gdy metoda bootstrapowa daje wynik 820.53. Różnica ta jest niewielka, co sugeruje, że bootstrap skutecznie odwzorowuje średnią próby.

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.1508713    6.333811

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.0, 820.8 )   (796.8, 820.7 )  
## Calculations and Intervals on Original Scale

Ordinary Nonparametric Bootstrap: Zastosowana została podstawowa, nieparametryczna wersja bootstrapu, polegająca na losowaniu z powtórzeniami próbek z danych wejściowych. Jest to najbardziej uniwersalna i powszechnie stosowana wersja metody bootstrap.Otrzymane wyniki wskazują, że bootstrap bardzo dobrze odwzorowuje średnią próby oryginalnej – różnica między średnią bootstrapową a oryginalną wynosi zaledwie -0.13, co oznacza, że metoda ta wprowadza minimalne odchylenie (bias). Oszacowany błąd standardowy wynosi 6.13 i jest bardzo zbliżony do klasycznego błędu standardowego, co świadczy o zgodności obu metod. Przedziały ufności wyznaczone na podstawie bootstrapu obejmują dwa rodzaje: klasyczny przedział oparty na założeniu normalności rozkładu średnich (796.6–820.7) oraz percentylowy przedział bootstrapowy (797.1–821.6). Oba przedziały są bardzo zbliżone, co wskazuje na stabilność wyników, ale percentylowy przedział jest bardziej elastyczny i lepiej odzwierciedla rzeczywisty rozkład danych, szczególnie jeśli dane mogą odbiegać od założeń normalności.

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.3     ✔ 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

# Zamiana tylko dla województw Pomorskiego i Podkarpackiego
dane <- dane %>%
  mutate(
    Wojewodztwo = case_when(
      Wojewodztwo == 22 ~ "Pomorskie",
      Wojewodztwo == 18 ~ "Podkarpackie",
      TRUE ~ as.character(Wojewodztwo) 
    )
  )

dane2 <- dane %>%
  filter(Wojewodztwo %in% c("Pomorskie", "Podkarpackie")) %>%
  mutate(Wojewodztwo = factor(Wojewodztwo)) 

print(unique(dane2$Wojewodztwo))

## [1] Podkarpackie Pomorskie   
## Levels: Podkarpackie Pomorskie

print(table(dane2$Wojewodztwo))

## 
## Podkarpackie    Pomorskie 
##         1560         1673

boot.t.test(Dochod_na_osobe ~ Wojewodztwo, R = 999, data = 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.8268 (30.83812) 
## 95 percent bootstrap percentile confidence interval:
##  -255.4252 -137.9800
## 
## 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

Analiza porównawcza dochodów na osobę między województwami Podkarpackim i Pomorskim wskazuje na istotne statystycznie różnice. Średni dochód w województwie Podkarpackim wynosi 648.21zł, podczas gdy w województwie Pomorskim jest to 844.17zł, co oznacza, że dochody w Pomorskim są wyższe o około 196zł. Zarówno test bootstrapowy, jak i klasyczny test t-Studenta potwierdzają, że różnica jest istotna statystycznie na poziomie ufności 95%. Bootstrapowy przedział ufności wynosi od -257.43 do -138.86, a klasyczny przedział ufności od -256.36 do -135.54 – oba przedziały są zbliżone i nie obejmują zera, co dodatkowo wzmacnia wniosek o istotności różnic. P-wartość dla testu bootstrapowego wynosi < 0.001, co oznacza bardzo silną istotność statystyczną różnicy. Wynik klasycznego testu t-Studenta daje podobny rezultat, z t-statystyką równą -6.36 i p-wartością 2.33e-10. Warto zauważyć, że błąd standardowy różnicy średnich wynosi 30.81, co wskazuje na stabilność estymacji. Obie metody analizy – bootstrapowa i klasyczna – są spójne, co świadczy o tym, że dane dobrze spełniają założenia klasycznego testu t-Studenta

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

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

Na podstawie przeprowadzonych testów chi-kwadrat (klasycznego i bootstrapowego) nie ma istotnych statystycznie różnic w proporcjach statusu studenta (Yes/No) pomiędzy mężczyznami a kobietami (p-value > 0.05). Oznacza to, że status studenta nie jest zależny od płci w analizowanych danych. Klasyczny test i test bootstrapowy dają bardzo podobne wyniki.

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 (Student) ma istotny wpływ na bilans na karcie kredytowej. Osoby będące studentami znacząco różnią się pod względem bilansu od osób, które nie są studentami. Pochodzenie etniczne (Ethnicity), stan cywilny (Married) i płeć (Gender) nie mają istotnego wpływu na średni bilans na karcie kredytowej. Model wyjaśnia jedynie niewielką część wariancji w danych (większość wariancji jest przypisana resztom), co sugeruje, że inne zmienne mogą być ważniejsze w wyjaśnieniu różnic w Balance.

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.9599600 0.9349349 0.0000000 0.8508509

Analiza ANOVA, zarówno klasyczna, jak i z zastosowaniem metody bootstrap, wykazała, że spośród analizowanych zmiennych jedynie status studenta (Student) ma istotny wpływ na średni bilans na karcie kredytowej (Balance). P-wartość dla tej zmiennej wynosi < 0.001, co oznacza bardzo silną istotność statystyczną. Studenci różnią się znacząco pod względem bilansu na karcie kredytowej w porównaniu z osobami, które nie są studentami. Natomiast zmienne takie jak pochodzenie etniczne (Ethnicity), stan cywilny (Married) oraz płeć (Gender) nie mają statystycznie istotnego wpływu na średni bilans na karcie kredytowej, co zostało potwierdzone przez wysokie p-wartości (wszystkie znacznie powyżej 0.05). Wyniki testu bootstrapowego są spójne z klasycznym testem ANOVA, co potwierdza ich stabilność i wiarygodność

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?

Analiza ANOVA, zarówno w wersji klasycznej, jak i z zastosowaniem bootstrapu, wykazała brak istotnych różnic w średnich bilansach na karcie kredytowej w zależności od pochodzenia etnicznego (Ethnicity). Średnie bilanse wyniosły 531.00 dla grupy African American, 512.31 dla grupy Asian oraz 518.50 dla grupy Caucasian. P-wartość testu ANOVA wyniosła 0.96, co jest znacznie wyższe od typowego progu istotności (0.05), wskazując, że różnice między grupami są statystycznie nieistotne. Podobne wyniki uzyskano z testu bootstrapowego, co potwierdza stabilność i wiarygodność analizy. Różnice między klasycznym testem ANOVA a jego wersją bootstrapową są minimalne w tym przypadku, co sugeruje, że dane spełniają założenia klasycznego testu ANOVA. W praktyce wyniki te wskazują, że pochodzenie etniczne nie ma wpływu na średni bilans na karcie kredytowej, a zastosowane metody analizy są zgodne i spójne.