Nieklasyczne metody statystyki

Bootstrap w estymacji i testowaniu

Bootstrap w estymacji punktowej i przedziałowej na przykładzie badania średniej.

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

2. Wygeneruj ciąg losowo wybranych liczb z zakresu od 1 do 10. Zastosuj losowanie ze zwracaniem i bez zwracania.

Zadania

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

library(haven)
dane <- dataset <- read_sav('dane.sav')
Dochod_na_osobe<-dane$Dochod_na_osobe

#punktowo klasycznie
mean(Dochod_na_osobe)

## [1] 808.5322

sd(Dochod_na_osobe)

## [1] 1064.67

se<-sd(Dochod_na_osobe)/sqrt(length(Dochod_na_osobe))
se

## [1] 6.146569

sep<-se/mean(Dochod_na_osobe)
sep

## [1] 0.007602133

#przedziałowo klasycznie
cat(paste('Przedziały ufności w danych surowych: \n'))

## Przedziały ufności w danych surowych:

mean(Dochod_na_osobe)-1.96*se   #dolna krawędź

## [1] 796.4849

mean(Dochod_na_osobe)+1.96*se   #górna krawędź

## [1] 820.5794

B=1e3
mean.Dochod_na_osobe=rep(0,B)
nobs=length(Dochod_na_osobe) 
#liczba obserwacji - ma być: 50, 250, 500, 1000, 10000

cat(paste('Przyjęto wielkość obserwacji z bootstrapa:', B))

## Przyjęto wielkość obserwacji z bootstrapa: 1000

 for (i in 1:B) 
{
 boot.data=sample(Dochod_na_osobe,nobs,replace=TRUE)
 mean.Dochod_na_osobe[i]=mean(boot.data)
 }

cat(paste('\nStatystyki opisowe bootstrap: \n'))

## 
## Statystyki opisowe bootstrap:

mean(boot.data)

## [1] 804.2442

sd(boot.data)

## [1] 1058.132

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

## [1] 6.108825

sep2<-se/mean(boot.data)
sep2

## [1] 0.007642665

hist(boot.data)

plot(density(boot.data))

cat(paste('Przedziały ufności w danych surowych: \n'))

## Przedziały ufności w danych surowych:

#quantile(boot.data,0.975)
#quantile(boot.data,0.025)

dolny<-mean(boot.data)-1.96*se2
gorny<-mean(boot.data)+1.96*se2
cbind(dolny,gorny)

##         dolny    gorny
## [1,] 792.2709 816.2175

cat(paste('Przedziałowy Efron: \n'))

## Przedziałowy Efron:

gorny<-quantile(boot.data,0.975)
dolny<-quantile(boot.data,0.025)
cbind(dolny,gorny)

##      dolny    gorny
## 2.5%     0 3047.567

cat(paste('Przedziałowy Hall: \n'))

## Przedziałowy Hall:

dolny<-2*mean(boot.data)-quantile(boot.data,0.025)
gorny<-2*mean(boot.data)-quantile(boot.data,0.975)
cbind(dolny,gorny)

##         dolny     gorny
## 2.5% 1608.488 -1439.079

A teraz z bootstrapem:

#install.packages('boot')
library(boot)

mean.boot=function(Dochod_na_osobe,idx) {
ans=mean(Dochod_na_osobe[idx])
ans
}

DOCHOD.mean.boot = boot(Dochod_na_osobe,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_na_osobe, statistic = mean.boot, R = 999)
## 
## 
## Bootstrap Statistics :
##     original    bias    std. error
## t1* 808.5322 0.5337111     5.97603

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.3, 819.7 )   (797.4, 820.6 )  
## Calculations and Intervals on Original Scale

Testy t studenta

Czy dochody na osobę różnią się istotnie w woj. pomorskim i podkarpackim?

#install.packages('MKinfer')
library(MKinfer) 

## Warning: pakiet 'MKinfer' został zbudowany w wersji R 4.3.3

library(tidyverse)

## Warning: pakiet 'tidyverse' został zbudowany w wersji R 4.3.3

## Warning: pakiet 'ggplot2' został zbudowany w wersji R 4.3.3

## Warning: pakiet 'dplyr' został zbudowany w wersji R 4.3.3

## Warning: pakiet 'lubridate' 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

## uruchamianie serwera httpd dla pomocy ... wykonano

dane2<- dane %>%
  filter(Wojewodztwo %in% c(16,22)) 

  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.6326 
## bootstrap difference of means (SE) = 18.81337 (53.75527) 
## 95 percent bootstrap percentile confidence interval:
##  -82.78772 127.37836
## 
## Results without bootstrap:
## t = 0.39855, df = 1363.8, p-value = 0.6903
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -84.90621 128.20264
## sample estimates:
## mean in group 16 mean in group 22 
##         865.8145         844.1663

  t.test(Dochod_na_osobe~Wojewodztwo, dane2)

## 
##  Welch Two Sample t-test
## 
## data:  Dochod_na_osobe by Wojewodztwo
## t = 0.39855, df = 1363.8, p-value = 0.6903
## alternative hypothesis: true difference in means between group 16 and group 22 is not equal to 0
## 95 percent confidence interval:
##  -84.90621 128.20264
## sample estimates:
## mean in group 16 mean in group 22 
##         865.8145         844.1663

wyniki<-boot.t.test(Dochod_na_osobe[as.numeric(dane$Wojewodztwo)==22],
               Dochod_na_osobe[as.numeric(dane$Wojewodztwo)==16],R=999,data=dane2)
wyniki

## 
##  Bootstrap Welch Two Sample t-test
## 
## data:  Dochod_na_osobe[as.numeric(dane$Wojewodztwo) == 22] and Dochod_na_osobe[as.numeric(dane$Wojewodztwo) == 16]
## number of bootstrap samples:  999
## bootstrap p-value = 0.6987 
## bootstrap difference of means (SE) = -21.26549 (53.55983) 
## 95 percent bootstrap percentile confidence interval:
##  -119.55926   76.39046
## 
## Results without bootstrap:
## t = -0.39855, df = 1363.8, p-value = 0.6903
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -128.20264   84.90621
## sample estimates:
## mean of x mean of y 
##  844.1663  865.8145

Test proporcji

#install.packages('ISLR')
library(ISLR)
data("Credit")
attach(Credit)

cat(paste('Test niezależności grup \n'))

## Test niezależności grup

tabelka<-table(Student,Gender)

#chisq.test(tabelka,simulate.p.value = TRUE, B = 2000)
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.3123

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

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

Różnice nie wyszły istotne. Zależność między resamplingiem a wynikami polega na tym, że im wyższa ilość resamplingu tym wynik zbliżony do tego bez bootstrapa. Test ANOVA względem testu ANOVA z bootstrapem nie różni się znacząco. Wyniki i istotności zmiennych są na podobnych poziomach.