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("C:/Users/czare/Desktop/NMS/Raport 6/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)
nobs1=50 #liczba próbek bootstrapowych - zwiększaj ich liczbę: 50, 250, 500, 1000, 10000

 for (i in 1:B) 
{
 boot.data1=sample(dochod,nobs1,replace=TRUE)
 mean.dochod[i]=mean(boot.data1)
 }

mean_boot1 <- mean(boot.data1)

sd_boot1 <- sd(boot.data1)

se1<-sd(boot.data1)/sqrt(length(boot.data1))

sep1<-se1/mean(boot.data1)

hist(boot.data1)

plot(density(boot.data1))

#przedziałowo bootstrapowo  95% ufności:
mean_boot1
## [1] 1103.475
sd_boot1
## [1] 1080.752
se1
## [1] 152.8415
sep1
## [1] 0.1385092
dolny<-mean(boot.data1)-1.96*se1
dolny
## [1] 803.9059
gorny<-mean(boot.data1)+1.96*se1
gorny
## [1] 1403.044

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’

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

 for (i in 1:B) 
{
 boot.data2=sample(dochod,nobs2,replace=TRUE)
 mean.dochod[i]=mean(boot.data2)
 }

mean_boot2 <- mean(boot.data2)

sd_boot2 <- sd(boot.data2)

se2<-sd(boot.data2)/sqrt(length(boot.data2))

sep2<-se2/mean(boot.data2)

hist(boot.data2)

plot(density(boot.data2))

#przedziałowo bootstrapowo  95% ufności:
mean_boot2
## [1] 824.6725
sd_boot2
## [1] 1152.874
se2
## [1] 72.91412
sep2
## [1] 0.08841585
dolny<-mean(boot.data2)-1.96*se2
dolny
## [1] 681.7608
gorny<-mean(boot.data2)+1.96*se2
gorny
## [1] 967.5842
B=999
mean.dochod=rep(0,B)
nobs3=500 #liczba próbek bootstrapowych - zwiększaj ich liczbę: 50, 250, 500, 1000, 10000

 for (i in 1:B) 
{
 boot.data3=sample(dochod,nobs3,replace=TRUE)
 mean.dochod[i]=mean(boot.data3)
 }

mean_boot3 <- mean(boot.data3)

sd_boot3 <- sd(boot.data3)

se3<-sd(boot.data3)/sqrt(length(boot.data3))

sep3<-se3/mean(boot.data3)

hist(boot.data3)

plot(density(boot.data3))

#przedziałowo bootstrapowo  95% ufności:
mean_boot3
## [1] 842.2532
sd_boot3
## [1] 1082.453
se3
## [1] 48.40875
sep3
## [1] 0.05747529
dolny<-mean(boot.data3)-1.96*se3
dolny
## [1] 747.3721
gorny<-mean(boot.data3)+1.96*se3
gorny
## [1] 937.1344
B=999
mean.dochod=rep(0,B)
nobs4=1000 #liczba próbek bootstrapowych - zwiększaj ich liczbę: 50, 250, 500, 1000, 10000

 for (i in 1:B) 
{
 boot.data4=sample(dochod,nobs4,replace=TRUE)
 mean.dochod[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:
mean_boot4
## [1] 831.9727
sd_boot4
## [1] 1188.711
se4
## [1] 37.59034
sep4
## [1] 0.04518218
dolny<-mean(boot.data4)-1.96*se4
dolny
## [1] 758.2957
gorny<-mean(boot.data2)+1.96*se4
gorny
## [1] 898.3496
B=999
mean.dochod=rep(0,B)
nobs5=10000 #liczba próbek bootstrapowych - zwiększaj ich liczbę: 50, 250, 500, 1000, 10000

 for (i in 1:B) 
{
 boot.data5=sample(dochod,nobs5,replace=TRUE)
 mean.dochod[i]=mean(boot.data5)
 }

mean_boot5 <- mean(boot.data5)

sd_boot5 <- sd(boot.data5)

se5<-sd(boot.data5)/sqrt(length(boot.data5))

sep5<-se5/mean(boot.data5)

hist(boot.data5)

plot(density(boot.data5))

#przedziałowo bootstrapowo  95% ufności:
mean_boot5
## [1] 819.3339
sd_boot5
## [1] 1123.24
se5
## [1] 11.2324
sep5
## [1] 0.01370919
dolny<-mean(boot.data5)-1.96*se5
dolny
## [1] 797.3184
gorny<-mean(boot.data2)+1.96*se5
gorny
## [1] 846.688
  1. Jakie są wnioski? Czy różnice są istotne? Jak wielkość resamplingu wpływa na wyniki?

Małe próby mają więszke błedy standartowe (np 50 - 109.1671 vs 10000 - 10.38394). Oraz mają większe przedziały ufnosci. W większych próbach po prostu ze względu na ilość danych mamy bardziej precyzyjne wyniki.

Porównojąc natomiast metode klasyczna do boostrapingi możemy dojść do wniosku, że w małcyh próbach typu 50 różnice są widoczne np bład standartowy 6,15 vs 109,17. Natomiast wraz ze wzrostem próby różnice są coraz mniejsze.

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.0557372    5.898963
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%   (797.0, 820.1 )   (797.1, 820.4 )  
## Calculations and Intervals on Original Scale

Bootstrapowy błąd standardowy - 6,067342 Klasyczny błąd standartowy - 6,15 Wyniki są bardzo zbliżone, więc można wywnioskowac, że metoda bootstrapu dobrze szacuje niepewność.

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) 
## Warning: pakiet 'MKinfer' został zbudowany w wersji R 4.3.3
library(tidyverse)
## Warning: pakiet 'ggplot2' został zbudowany w wersji R 4.3.3
## ── 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.5.1     ✔ 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) = -193.8267 (30.74618) 
## 95 percent bootstrap percentile confidence interval:
##  -254.8118 -135.7849
## 
## 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

Oba testy prowadzą do identycznych wniosków. Występują różnice w zarobkach pomiędzy województwem pomorskim a podkarpackim.

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

Wartość p w testach wynosi 0,3363 i 0,3504, zatem oba wyniki są większe od 0,05. Zatem mamy brak podstaw od odrzucenia hipotezy zerowej, więc status studenta nie zależy od płci.

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

A teraz z włączonym bootstrapem:

library(lmboot)
## Warning: pakiet 'lmboot' został zbudowany w wersji R 4.3.3
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.9579580 0.9179179 0.0000000 0.8958959

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?

Test ANOVA wykazał, że jedynie zmienna student jest istotna. Zatem można wnioskować, że status studenta wpływa na średni bilans na karcie kredytowej. Test z bootstrapem jedynie potwierdził powyższą tezę, wszystkie zmienne nieistotne za wyjątkiem statusu studenta.

Wizualizacja potwierdza, że nie ma różnic pomiędzy zmienną grupa etniczna. Mają bardzo zbliżone średnie bilansowe.