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

Pierwszy krok- zainstalowanie niezbędnych pakietów oraz wgranie danych.

library(foreign)
dane <- read.spss("C:/Users/macie/OneDrive/Pulpit/dane.sav", to.data.frame = TRUE)
attach(dane)
dochod<-Dochod_na_osobe

Wyniki klasycznej metody estymacji punktowej:

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

Wyniki klasycznej estymacji przedziałowej, przy czym 1.96 to kwantyl rozkładu normalnego dla 95% ufności. Podane poniżej wyniki to krawędź dolna oraz górna.

mean(dochod)-1.96*se
## [1] 796.4849
mean(dochod)+1.96*se
## [1] 820.5794

Teraz zrobimy to samą metodą bootstrap przy uwzględnieniu 50 próbek (nobs). Objaśnienia skrótów: mean_boot- średnia próbki, sd_boot- odchylenie standardowe dla próbki, se2- błąd standardowy dla próbki, sep2- błąd standardowy wyrażony w procentach, dolny oraz gorny- dolna oraz górna granica przedziału ufności.

B=999
mean.dochod=rep(0,B)
nobs=50 

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

mean_boot <- mean(boot.data)
mean_boot
## [1] 741.3161
sd_boot <- sd(boot.data)
sd_boot
## [1] 660.2896
se2<-sd(boot.data)/sqrt(length(boot.data))
se2
## [1] 93.37906
sep2<-se2/mean(boot.data)
sep2
## [1] 0.1259639
hist(boot.data)

plot(density(boot.data))

#przedziałowo bootstrapowo  95% ufności:
dolny<-mean(boot.data)-1.96*se2
dolny
## [1] 558.2932
gorny<-mean(boot.data)+1.96*se2
gorny
## [1] 924.3391

Porównanie wyników dla metody klasycznej oraz bootstrappingu przy próbie o wielkości 50:

Dla wielkości próby ‘r nobs=50’ otrzymujemy: - średnią próbkową: 792.2717 - odchylenie standardowe: 805.5613 - błąd standardowy: 113.9236 - procentowy błąd standardowy: 0.007758157 - dolną krawędź przedziału ufności: 568.9815 - górną krawędź przedziału ufności: 1015.562

Dla klasycznej metody wykorzystującej całość zbioru danych otrzymano następujące wartości statystyk: Średnia wartość dotycząca całości zbioru danych: 808.5322 Odchylenie standardowe: 1064.67 Błąd standardowy: 6.146569 Procentowy błąd standardowy: 0.007602133 Dolna krawędź przedziału ufności: 796.4849 Górna krawędź przedziału ufności: 820.5794

Wykorzystany zbiór danych obejmuje ponad 30 tysięcy obserwacji, wyniki bootstrappingu przy wykorzystaniu zaledwie 50 obserwacji należących do tego zbioru danych mogą zatem znacząco się różnić. Przy wykonanym pomiarze odnotowano stosunkowo zbliżoną do wartości faktycznej wartość oszacowanej średniej- wartość średniej prókowej wynosiła 792.2717, co stanowi około 98% średniej dla całej populacji. Do istotnych konsekwencji wykorzystania bootstrappingu można zaliczyć bardzo silny wzrost wartości błędu standardowego, a także istotne poszerzenie zakresu przedziału ufności. Występowanie takiej sytuacji można uzasadnić znacząco większym znaczeniem obserwacji odstających przy wykorzystaniu do ekstrapolacji bardzo małej części analizowanego zbioru danych, a uwzględnienie przy bootstrappingu większej liczby obserwacji może pomóc w zmniejszeniu dysproporcji pomiędzy wartośiami statystyk takich jak błąd standardowy pomiędzy wynikami bootstrappingu a wynikami dotyczącymi całości populacji. W tym celu wykorzystany bootstrapping uwzględniający 500 oraz 1000 próbek.

B=999
mean.dochod=rep(0,B)
nobs=500 

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

mean_boot <- mean(boot.data)
mean_boot
## [1] 782.2754
sd_boot <- sd(boot.data)
sd_boot
## [1] 1191.358
se2<-sd(boot.data)/sqrt(length(boot.data))
se2
## [1] 53.27914
sep2<-se2/mean(boot.data)
sep2
## [1] 0.06810791
hist(boot.data)

plot(density(boot.data))

#przedziałowo bootstrapowo  95% ufności:
dolny<-mean(boot.data)-1.96*se2
dolny
## [1] 677.8483
gorny<-mean(boot.data)+1.96*se2
gorny
## [1] 886.7026
B=999
mean.dochod=rep(0,B)
nobs=1000 

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

mean_boot <- mean(boot.data)
mean_boot
## [1] 796.7097
sd_boot <- sd(boot.data)
sd_boot
## [1] 1111.789
se2<-sd(boot.data)/sqrt(length(boot.data))
se2
## [1] 35.15786
sep2<-se2/mean(boot.data)
sep2
## [1] 0.04412882
hist(boot.data)

plot(density(boot.data))

#przedziałowo bootstrapowo  95% ufności:
dolny<-mean(boot.data)-1.96*se2
dolny
## [1] 727.8003
gorny<-mean(boot.data)+1.96*se2
gorny
## [1] 865.6191

Dla wszystkich wyników bootstrappingu odnotowano stosunkowo podobne wartości dotyczące średniej dla próbki, oscylujące pomiędzy wartościami 790 oraz 800, zbliżone do średniej wartości dla całej populacji. Przy każdym zwiększaniu populacji przy bootstrappingu odnotowano zmniejszenie się wartości błędu standardowego- podczas gdy przy wielkości populacji wynoszącej 50 wynosił on 113.9236, przy wykorzystaniu 500 obserwacji było to 45.32871, a dla 1000 obserwacjach odnotowano wartość 27.9694, najbardziej zbliżoną do średniej wartości błędu standardowego dla całej populacji wynoszącego 6.146569. Podobny trend, gdzie wraz ze wzrostem odsetka populacji podlegającego bootstrappingowi wartości danej statystyki zbliżały się do tych cechujących całość populacji, odnotowano dla dolnych oraz górnych krawędzi przedziału ufności. Przy wykorzystaniu 50, 500 oraz 1000 wartości dolny próg wynosił odpowiednio 568.9815, 701.6668 oraz 741.192, zbliżając się do wartości dla całej populacji wynoszącej 796.4849, a dla górnego progu jego wartości zbliżające się do wartości dla całości danych wynoszącej 820.5794 wynosiły kolejno 1015.562, 879.3553 oraz 850.8321. Można zatem dojść do konkluzji, że wraz ze zwiększeniem się rozmiaru próbki wykorzystywanej do bootstrappingu zmniejsza się różnica pomiędzy wartościami takich statystyk jak błąd standardowy oraz progi przedziału ufności pomiędzy całą populacją a bootstrappingiem dla populacji badanej, jednak zwiększenie rozmiaru populacji bootsrappingowanej z 50 do 500 przyniosło dużo większy efekt niż zwiększenie jej z 500 do 1000.

Zadanie 2. Estymacja błędu standardowego średniej próbkowej metodą bootstrap dla różnych wielkości próby.

Pracę z funkcjami dotyczącymi bootstrappingu zaczniemy od uruchomienia procedury bootstrappingu.

?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"

Statystyki dotyczące bootstrappingu:

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.01606554    6.006097

Wykresy:

plot(DOCHOD.mean.boot)

Przedziały ufności:

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

Przy zastosowaniu funkcji boot, odnotowano wartości odpowiednio błędu standardowego, dolnego oraz górnego progu przedziału ufności wynoszące odpowiednio 6.047981, 796.6 oraz 820.3, zbliżone do wartości otrzymanych przy zastosowaniu klasycznej metody dla całości zbioru danych, które wynosiły odpowiednio 6.146569, 796.4849 oraz 820.5794, co wskazuje na bardzo dobre dostowanie tej metody bootstrappingu do wykorzystanych danych.

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
## Warning: pakiet 'dplyr' został zbudowany w wersji R 4.3.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ 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) = -195.3439 (30.77184) 
## 95 percent bootstrap percentile confidence interval:
##  -257.4711 -136.7949
## 
## 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

Zaróno przy teście t-studenta, jak i bootstrapowym na podstawie zbliżonych do zera wartości p-value należy przyjąć hipotezę alternatywną, zgodnie z którą średnie dochody na osobę istotnie różnią się w obu wybranych województw, a różnicę tą potwierdzają ujemne wartości dotyczące przedziałów ufności, które są w przypadku obu testów zbliżone do siebie oraz ujemne.

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)
## Warning: pakiet 'ISLR' został zbudowany w wersji R 4.3.2
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.3043
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

W przypadku testu Chi-kwadrat zarówno przy wykorzystaniu bootstrapu, jak i przy niewykorzystaniu go nie można odrzucić hipotezy zerowej mówiącej w tym przypadku, że status studenta nie różni się w zależności od płci w istotnym statystycznie stopniu. Na podstawie tego zmienne Student oraz Gender można zatem uznać za niepowiązane od tego, niezależnie od tego czy wykorzystujemy bootstrapping.

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.9479479 0.9329329 0.0000000 0.8888889

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.

Wykres zależności pomiędzy bilansem na karcie a pochodzeniem:

library(ggstatsplot)
## Warning: pakiet 'ggstatsplot' został zbudowany w wersji R 4.3.2
## 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
)

Wykres zależności pomiędzy bilansem na karcie a stanem małżeńskim:

library(ggstatsplot)
ggbetweenstats(data=Credit,
  y=Balance,
  x=Married,
  nboot=999  #liczba prób bootstrapowych
)

Wykres zależności pomiędzy bilansem na karcie a statusem studenta:

library(ggstatsplot)
ggbetweenstats(data=Credit,
  y=Balance,
  x=Student,
  nboot=999  #liczba prób bootstrapowych
)

Wykres zależności pomiędzy bilansem na karcie a płcią:

library(ggstatsplot)
ggbetweenstats(data=Credit,
  y=Balance,
  x=Gender,
  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?

Wartości p-value dla testu Anova zarówno przy wykorzystaniu bootstrappingu, jak i przy niekorzystaniu z niego wskazują na to, że spośród analizowanych zmiennych mogących potencjalnie wpływać na poziom środków na karcie kredytowej jedynie status studenta ma istotny wpływ na poziom środków na karcie. Jak wskazują wartości p-value dotyczące testu Anova dla zmiennych dotyczących pochodzenia, stanu cywilnego oraz płci, a także przedstawione wykresy, średnie wartości dochodów można uznać za niepowiązane z tymi zmiennymi. Dla zmiennych Ethnicity, Gender oraz Married odnotowano nieznacznie wyższe wartości p-value przy włączonym bootstrapingu.

Zdecydowano się także na porównanie wyników przy zmianie liczby prób bootstrappingowych na 99 oraz 9999 w celu dokonania porównania pomiędzy pierwotnym pomiarem wykorzystującym 999 pomiarów, gdzie p-value dotyczące zmiennych Ethnicity, Married, Student oraz Gender wynosiły odpowiednio 0.9639640, 0.9459459, 0.0000000 oraz 0.8728729, jednak nie zaobserwowano jednoznacznego powiązania pomiędzy wzrostem liczby prób a zmianą p-value.

library(lmboot)
anova_boot<-ANOVA.boot(Balance ~ Ethnicity + Married + Student + Gender,data=Credit,B=99)
## Warning in ANOVA.boot(Balance ~ Ethnicity + Married + Student + Gender, : This function has only been fully tested for one-way and two-way ANOVA.
## Warning in ANOVA.boot(Balance ~ Ethnicity + Married + Student + Gender, : Number of bootstrap samples is recommended to be more than the number of observations.
anova_boot$`p-values`
## [1] 0.9696970 0.9090909 0.0000000 0.8585859
library(lmboot)
anova_boot<-ANOVA.boot(Balance ~ Ethnicity + Married + Student + Gender,data=Credit,B=9999)
## 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.9571957 0.9379938 0.0000000 0.8730873