Lab6
Paulina Fereniec
2025-01-12
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).
library(haven)
data <- read_sav("dane.sav")
data## # A tibble: 30,003 × 9
## Wojewodztwo Powierzchnia_mieszkalna Oplaty_miesieczne_GD
## <dbl+lbl> <dbl> <dbl>
## 1 12 [Malopolskie] 80 112.
## 2 10 [Lodzkie] 18 741.
## 3 2 [Dolnoslaskie] 29 188
## 4 14 [Mazowieckie] 28 184.
## 5 32 [Zachodniopomorskie] 53 247.
## 6 30 [Wielkopolskie] 24 140.
## 7 12 [Malopolskie] 50 624
## 8 14 [Mazowieckie] 20 205.
## 9 24 [Slaskie] 46 235.
## 10 10 [Lodzkie] 65 393.
## # ℹ 29,993 more rows
## # ℹ 6 more variables: Calkowity_dochod_GD_GUS <dbl>,
## # Calkowite_wydatki_GD_GUS <dbl>, Dochod_na_osobe <dbl>,
## # Calkowity_dochod_GD_dane <dbl>, Wiek_glowy_rodziny <dbl>, Wiek <dbl>Zadanie 1. Estymacja błędu standardowego średniej próbkowej metodą bootstrap.
Załóżmy, że chcemy badać populację o rozkładzie normalnym ze średnią 3 i odchyleniem standardowym 1. Wygeneruj przykładową populację liczności 25 o takim rozkładzie (funkcja rnorm()). Wygeneruj 1000 prób o liczności 25 z otrzymanych danych (funkcja replicate()). Oblicz średnią każdej z 1000 wylosowanych prób bootstrapowych oraz przedstaw histogram rozkładu tych średnich. Porównaj ten rozkład z rozkładem cechy (rozkładem normalnym N (3;1)). Oceń, czy przybliżenie rozkładu metodą bootstrap jest adekwatne (czy kształty obu rozkładów są podobne) Porównaj średnią wygenerowanej próby wyjściowej ze średnią wektora średnich z prób bootstrapowych. Oceń ich różnicę. Przy użyciu funkcji boot (pakiet boot) wyznacz błąd standardowy dla średniej.
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.007602133przedziałowo klasycznie
mean(dochod)-1.96*se #dolna krawędź## [1] 796.4849mean(dochod)+1.96*se #górna krawędź## [1] 820.5794B=999
mean.dochod=rep(0,B)
nobs=length(dochod)liczba obserwacji - ma być: 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.data)## [1] 806.7525sd(boot.data)## [1] 1036.584se2<-sd(boot.data)/sqrt(length(boot.data))
se2## [1] 5.984421sep2<-se/mean(boot.data)
sep2## [1] 0.007618903hist(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
cbind(dolny,gorny)## dolny gorny
## [1,] 795.023 818.482Przedziałowo Efron
dolny<-quantile(boot.data,0.975)
gorny<-quantile(boot.data,0.025)
cbind(dolny,gorny)## dolny gorny
## 97.5% 3065.7 0Przedziałowo 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% 1613.505 -1452.195Bootstrap
library(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.3890116 5.890107plot(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.4, 820.5 ) (795.9, 820.0 )
## Calculations and Intervals on Original ScalePorównaj wyniki z obu metod.
Jakie są wnioski? Czy różnice są istotne? Jak wielkość resamplingu wpływa na wyniki?
Większa liczba prób (większy resampling) prowadzi do bardziej stabilnych i dokładnych estymacji. Wartość błędu standardowego jest mniejsza, a przedziały ufności są węższe.
Zadanie 2. Estymacja przedziałowa średniej metodą bootstrap.
Wygeneruj próbkę losową 50 obserwacji wybranej zmiennej z pliku dane.sav (funkcja sample()). Jako populację potraktuj cały zbiór N-obserwacji w pliku. Oszacuj przedziałowo (klasycznie – stosując aproksymację rozkładem normalnym) oraz nieklasycznie (stosując funkcję boot.ci()) średnią wybranej zmiennej w próbce. Oblicz oraz porównaj oba błędy standardowe (względny %) estymacji przedziałowej klasycznej i nieklasycznej. Podaj wnioski.
mean.boot=function(dochod,idx) {
ans=mean(dochod[idx])
ans
}DOCHOD.mean.boot = boot(dochod,statistic=mean.boot, R=49)
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 = 49)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 808.5322 -0.6824419 6.433158plot(DOCHOD.mean.boot)
boot.ci(DOCHOD.mean.boot,conf=0.95,type=c("norm","perc"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 49 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = DOCHOD.mean.boot, conf = 0.95, type = c("norm",
## "perc"))
##
## Intervals :
## Level Normal Percentile
## 95% (796.6, 821.8 ) (792.6, 822.2 )
## Calculations and Intervals on Original Scale
## Some percentile intervals may be unstable### TEST T studenta
library(MKinfer) ## Warning: pakiet 'MKinfer' został zbudowany w wersji R 4.4.2library(tidyverse)## Warning: pakiet 'tidyverse' został zbudowany w wersji R 4.4.2## Warning: pakiet 'ggplot2' został zbudowany w wersji R 4.4.2## Warning: pakiet 'tibble' został zbudowany w wersji R 4.4.2## Warning: pakiet 'tidyr' został zbudowany w wersji R 4.4.2## Warning: pakiet 'dplyr' został zbudowany w wersji R 4.4.2## ── 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?boot.t.test## uruchamianie serwera httpd dla pomocy ... wykonanodane<-dataTest Proporcji
Przykład 1.
Czy mają istotnie różne spalanie? 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?
dane2<-dane%>%
filter(Wojewodztwo %in% c(16,22))
wyniki<-boot.t.test(Dochod_na_osobe~Wojewodztwo,R=999,data=dane2)
wyniki##
## Bootstrap Welch Two Sample t-test
##
## data: Dochod_na_osobe by Wojewodztwo
## number of bootstrap samples: 999
## bootstrap p-value = 0.7227
## bootstrap difference of means (SE) = 24.5738 (53.93514)
## 95 percent bootstrap percentile confidence interval:
## -79.37684 139.17219
##
## 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.1663z wizualizacją:
library(ggstatsplot)## Warning: pakiet 'ggstatsplot' został zbudowany w wersji R 4.4.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.03167ggbetweenstats(
data=dane2,
x=Wojewodztwo,
y=Dochod_na_osobe,
nboot=999
)
Przykład 2. Test Chi2 dla dwóch zmiennych jakościowych
library(ISLR)## Warning: pakiet 'ISLR' został zbudowany w wersji R 4.4.2attach(dane2)## Następujące obiekty zostały zakryte z data:
##
## Calkowite_wydatki_GD_GUS, Calkowity_dochod_GD_dane,
## Calkowity_dochod_GD_GUS, Dochod_na_osobe, Oplaty_miesieczne_GD,
## Powierzchnia_mieszkalna, Wiek, Wiek_glowy_rodziny, WojewodztwoCzy status studenta różni się istotnie wg płci?
data("Credit")
?Credit
attach(Credit)Czy status studenta różni się istotnie wg płci?
tabelka<-table(Student,Gender)
?chisq.test()
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.3253ggbarstats(
data=Credit,
x=Gender,
y=Student
)
ggpiestats(
data=Credit,
x=Gender,
y=Student
)
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? Czy wyniki z włączonym bootstrapem różnią się od wyników pojedynczego testu ANOVA?
ANOVA bez bootstrapa:
## 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)## Warning: pakiet 'lmboot' został zbudowany w wersji R 4.4.2anova_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.9573957 0.9362936 0.0000000 0.8793879Wizualizacja
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)
ggbetweenstats(data=Credit,
y=Balance,
x=Ethnicity,
nboot=9999
)
library(ggstatsplot)
ggbetweenstats(data=Credit,
y=Balance,
x=Ethnicity,
nboot=999
)
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? Anova z bootstrapem daje bardziej stabilne wyniki, ale
wartości p są bliskie wartością ze zwykłej ANOVA. Wartość p dla testu
ANOVA z bootstrapem jest mniejsza niż dla testu ANOVA bez
bootstrapa.