Bootstrap
Mikołaj Mogielnicki
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/Mikołaj/Documents/studia/NMS/dane.sav", to.data.frame = TRUE)
attach(dane)
dochod<-Dochod_na_osobe
#punktowo klasycznie
mean(dochod)## [1] 808.5322sd(dochod)## [1] 1064.67se<-sd(dochod)/sqrt(length(dochod))
se## [1] 6.146569sep<-se/mean(dochod)
sep## [1] 0.007602133#przedziałowo klasycznie
mean(dochod)-1.96*se #dolna krawędź## [1] 796.4849mean(dochod)+1.96*se #górna krawędź## [1] 820.5794# 1.96 to kwantyl rozkładu normalnego dla 95% ufnościTeraz 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*se2Wyniki 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’
#średnia próbkowa klasycznie
mean(dochod)## [1] 808.5322#średnia próbkowa bootstrap
mean_boot## [1] 806.7861Porównaj wyniki z obu metod. Próba n = 50 Średnia próbkowa Klasycznie: 808,5322 Bootstrap: 910,1497
Odchylenie standardowe Klasycznie: 1064.67 Bootstrap: 1051,86
błąd standardowy Klasycznie: 6,14 Bootstrap: 148,75
procentowy błąd standardowy Klasycznie: 0,0076 Bootstrap: 0,0067
dolną krawędź przedziału ufności Klasycznie: 796.4849 Bootstrap: 618,5871
górna krawędź przedziału ufności Klasycznie: 820.5794 Bootstrap: 1201,7123
Próba n = 10000 Średnia próbkowa Klasycznie: 808,5322 Bootstrap: 813,2188
Odchylenie standardowe Klasycznie: 1064.67 Bootstrap: 1089,7752
błąd standardowy Klasycznie: 6,14 Bootstrap: 6,14
procentowy błąd standardowy Klasycznie: 0,0076 Bootstrap: 10,8977
dolną krawędź przedziału ufności Klasycznie: 796.4849 Bootstrap: 791,8592
górna krawędź przedziału ufności Klasycznie: 820.5794 Bootstrap: 834,5784
Jakie są wnioski? Czy różnice są istotne? Jak wielkość resamplingu wpływa na wyniki?
Przy małej wartości próby różnice są większe, natomiast przy wartości n = 10000 różnice są dość małe. Można zatem wyciągnąć wniosek że więlkość resamplingu wpływa na wyniki, im większa próba tym wyniki będą lepsze.
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
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.1585516 6.069679plot(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.8, 820.6 ) (796.8, 820.0 )
## Calculations and Intervals on Original ScaleTesty 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)
?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.0922 (30.76422)
## 95 percent bootstrap percentile confidence interval:
## -253.5272 -129.5912
##
## 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.1663Zarówno według testu t-studenta oraz testu bootstrapowego dochody na osobę różnią się istotnie dla tych województw. ponieważ w obu testach p-value jest bardzo małe.W obu testach różnica między średnimi dochodami województwami wynosi około 195 złotych na korzyść pomorskiego. Przedział ufności dla testu t wyniósł [-256.3683, -135.5439] a dla testu bootstrapowego [-250.5698, -131.5983], więc są do siebie bardzo zbliżone. Oba testy zwracają zbliżone i spójne wyniki więc dzięki swojej elastyczności bootstrap może być dobrym rozwiązaniem w takim przypadku.
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.3353chisq.test(tabelka) # różnice???##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tabelka
## X-squared = 0.87218, df = 1, p-value = 0.3504Zarówno test bootstrapowy jak i ten zwykły mają p-value na zbliżonym poziomie (0,3273 i 0,3504). Oba testy wskazują na brak istotnych różnic statusu studenta w zależności 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 ' ' 1A teraz z włączonym bootstrapem:
library(lmboot)
anova_boot<-ANOVA.boot(Balance ~ Ethnicity + Married + Student + Gender,data=Credit,B=999)
anova_boot$`p-values`## [1] 0.9469469 0.9389389 0.0000000 0.8828829W 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)
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?
Zarówno ANOVA z bootstrapem oraz ANOVA bez niego wskazują na istotną różnice w bilansie na karcie kredytowej jedynie w zależności statusu studenta.Dla pozostałych czynników wartości p-value są bardzo wysokie więc nie ma istotnych różnic. Wielkość resamplingu nieznacznie zmienia wartości p-value jednak nawet dla bardzo małych ilości próbek istotnie bilans istotnie różni się dla statusu studenta.