Punktowo klasycznie

wartosc_domu <- Boston$medv

srednia_wartosc <- mean(wartosc_domu)
odchylenie_standardowe <- sd(wartosc_domu)
se <- odchylenie_standardowe / sqrt(length(wartosc_domu))
sep <- se / srednia_wartosc

dolna_klasyczna <- srednia_wartosc - 1.96 * se
gorna_klasyczna <- srednia_wartosc + 1.96 * se

cat("Średnia wartość domu:", srednia_wartosc, "\n")
## Średnia wartość domu: 22.53281
cat("Odchylenie standardowe:", odchylenie_standardowe, "\n")
## Odchylenie standardowe: 9.197104
cat("Odchylenie standardowe średniej (SE):", se, "\n")
## Odchylenie standardowe średniej (SE): 0.4088611
cat("Względny błąd standardowy (SEP):", sep, "\n")
## Względny błąd standardowy (SEP): 0.01814515
cat("95% przedział ufności:", dolna_klasyczna, "-", gorna_klasyczna, "\n")
## 95% przedział ufności: 21.73144 - 23.33417

Estymacja średniej (bootstrap)

## Średnia z bootstrapu: 22.57654
## Odchylenie standardowe z bootstrapu: 1.262548
## Odchylenie standardowe średniej (SE): 0.1785513
## Względny błąd standardowy (SEP): 0.007908711
## 95% przedział ufności: 22.22658 - 22.9265

Porównanie wyników metod klasycznych i bootstrapowych

Porównanie wyników obu metod pokazuje, że średnia wartość domu oszacowana klasycznie (22.53) i bootstrapowo (22.47) są bardzo zbliżone, co wskazuje na zgodność obu podejść. Różnice w odchyleniach standardowych wynikają z charakterystyki tych metod: klasyczne SD (9.20) odnosi się do całej próby, podczas gdy bootstrapowe SD (1.27) mierzy zmienność średnich z próbek bootstrapowych.

Odchylenie standardowe średniej (SE) oraz względny błąd standardowy (SEP) są niższe w metodzie bootstrapowej, co sugeruje większą precyzję tej techniki. Klasyczny przedział ufności (21.73 - 23.33) jest szerszy niż bootstrapowy (22.12 - 22.83), co oznacza, że bootstrap ocenia niepewność jako mniejszą, ale to może być zależne od wielkości próby i liczby iteracji.

Ogólnie, obie metody dają podobne wyniki, ale bootstrap oferuje bardziej szczegółową analizę niepewności, zwłaszcza w niestandardowych sytuacjach. Jednak wyniki bootstrapowe mogą być mniej stabilne przy małej liczbie iteracji, dlatego należy stosować większy resampling.

srednia_bootstrap = function(wartosc_domu, idx) {
  srednia <- mean(wartosc_domu[idx])
  return(srednia)
}

srednia_bootstrap_wynik <- boot(wartosc_domu, statistic = srednia_bootstrap, R = 999)

class(srednia_bootstrap_wynik)
## [1] "boot"
names(srednia_bootstrap_wynik)
##  [1] "t0"        "t"         "R"         "data"      "seed"      "statistic"
##  [7] "sim"       "call"      "stype"     "strata"    "weights"
print(srednia_bootstrap_wynik)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = wartosc_domu, statistic = srednia_bootstrap, R = 999)
## 
## 
## Bootstrap Statistics :
##     original       bias    std. error
## t1* 22.53281 -0.003918741    0.405842
plot(srednia_bootstrap_wynik)

przedzial_ufnosci <- boot.ci(srednia_bootstrap_wynik, conf = 0.95, type = c("norm", "perc"))
print(przedzial_ufnosci)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = srednia_bootstrap_wynik, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (21.74, 23.33 )   (21.77, 23.33 )  
## Calculations and Intervals on Original Scale

Porównania wartości mieszkań (medv) w danych Boston, w zależności od zmiennej chas (czy obszar graniczy z rzeką Charles), z wykorzystaniem testu t-Studenta i testu bootstrapowego.

library(MASS)
library(MKinfer)
## Warning: pakiet 'MKinfer' został zbudowany w wersji R 4.3.3
library(dplyr)
## 
## Dołączanie pakietu: 'dplyr'
## Następujący obiekt został zakryty z 'package:MASS':
## 
##     select
## Następujące obiekty zostały zakryte z 'package:stats':
## 
##     filter, lag
## Następujące obiekty zostały zakryte z 'package:base':
## 
##     intersect, setdiff, setequal, union
dane_chas <- Boston %>% 
  filter(chas %in% c(0, 1)) %>%  # Filtracja dla `chas` = 0 i `chas` = 1
  mutate(chas = factor(chas, labels = c("Nie graniczy", "Graniczy")))

t_test_wynik <- t.test(medv ~ chas, data = dane_chas, var.equal = FALSE)

print(t_test_wynik)
## 
##  Welch Two Sample t-test
## 
## data:  medv by chas
## t = -3.1133, df = 36.876, p-value = 0.003567
## alternative hypothesis: true difference in means between group Nie graniczy and group Graniczy is not equal to 0
## 95 percent confidence interval:
##  -10.476831  -2.215483
## sample estimates:
## mean in group Nie graniczy     mean in group Graniczy 
##                   22.09384                   28.44000
boot_test_wynik <- boot.t.test(medv ~ chas, R = 999, data = dane_chas)

print(boot_test_wynik)
## 
##  Bootstrap Welch Two Sample t-test
## 
## data:  medv by chas
## number of bootstrap samples:  999
## bootstrap p-value = 0.002002 
## bootstrap difference of means (SE) = -6.374534 (2.001305) 
## 95 percent bootstrap percentile confidence interval:
##  -10.526145  -2.361096
## 
## Results without bootstrap:
## t = -3.1133, df = 36.876, p-value = 0.003567
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -10.476831  -2.215483
## sample estimates:
## mean in group Nie graniczy     mean in group Graniczy 
##                   22.09384                   28.44000

Średnia wartość mieszkań w obszarach graniczących z rzeką wynosiła 28.44, podczas gdy w obszarach niegraniczących wynosiła 22.09. Oba testy wykazały statystycznie istotną różnicę między grupami: p-wartość wyniosła 0.0036 dla testu t-Studenta i 0.0020 dla testu bootstrapowego. Różnica średnich oszacowana bootstrapowo wyniosła -6.37 z przedziałem ufności (-10.47, -2.68).Wyniki obu metod były zgodne, co potwierdza, że mieszkania w obszarach graniczących z rzeką są istotnie droższe.

Test proporcji Czy zmienna chas (granica z rzeką Charles) różni się istotnie w zależności od kategorii rad (indeks dostępności autostrad)?

tabelka <- table(Boston$chas, Boston$rad)

chi2_simulacja <- chisq.test(tabelka, simulate.p.value = TRUE, B = 2000)

chi2_klasyczny <- chisq.test(tabelka)
## Warning in chisq.test(tabelka): Aproksymacja chi-kwadrat może być niepoprawna
print(chi2_simulacja)
## 
##  Pearson's Chi-squared test with simulated p-value (based on 2000
##  replicates)
## 
## data:  tabelka
## X-squared = 13.898, df = NA, p-value = 0.08696
print(chi2_klasyczny)
## 
##  Pearson's Chi-squared test
## 
## data:  tabelka
## X-squared = 13.898, df = 8, p-value = 0.08447

Zarówno test z symulacją, jak i test klasyczny wskazują, że zmienne chas (czy obszar graniczy z rzeką Charles) i rad (indeks dostępności autostrad) nie są istotnie zależne na poziomie istotności 0.05. Oznacza to, że liczba autostrad w pobliżu nie wpływa istotnie na to, czy dany obszar graniczy z rzeką Charles.

Testy ANOVA Sprawdzimy, czy średnia wartość mieszkań (medv) różni się istotnie w zależności od kilku czynników: dostępu do autostrad (rad), odsetka populacji o niskim statusie społecznym (lstat jako kategoria), sąsiedztwa rzeki Charles (chas), oraz liczby pokoi (rm jako kategoria).

Test ANOVA z wieloma czynnikami

library(lmboot)
## Warning: pakiet 'lmboot' został zbudowany w wersji R 4.3.3
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
Boston$lstat_cat <- cut(Boston$lstat, breaks = 3, labels = c("Niski", "Średni", "Wysoki"))
Boston$rm_cat <- cut(Boston$rm, breaks = 3, labels = c("Mało pokoi", "Średnio pokoi", "Dużo pokoi"))

model1 <- lm(medv ~ rad + lstat_cat + chas + rm_cat, data = Boston)
anova_wynik <- anova(model1)

print(anova_wynik)
## Analysis of Variance Table
## 
## Response: medv
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## rad         1  6221.1  6221.1 204.619 < 2.2e-16 ***
## lstat_cat   2 10256.7  5128.4 168.677 < 2.2e-16 ***
## chas        1  1084.9  1084.9  35.684 4.415e-09 ***
## rm_cat      2  9982.2  4991.1 164.161 < 2.2e-16 ***
## Residuals 499 15171.3    30.4                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Test ANOVA z bootstrapem

anova_boot <- ANOVA.boot(medv ~ rad + lstat_cat + chas + rm_cat, data = Boston, B = 999)
## Warning in ANOVA.boot(medv ~ rad + lstat_cat + chas + rm_cat, data = Boston, : This function has only been fully tested for one-way and two-way ANOVA.
print(anova_boot$`p-values`)
## [1] 0 0 0 0

ANOVA 1-czynnikowa z wizualizacją

ggbetweenstats(
  data = Boston,
  y = medv,
  x = rm_cat,
  nboot = 999,  # Liczba prób bootstrapowych
  title = "Wartość mieszkań w zależności od liczby pokoi",
  xlab = "Liczba pokoi (kategorie)",
  ylab = "Wartość mieszkań (tys. USD)"
)

Największy wpływ na wartość mieszkań miała liczba pokoi (rm_cat). Podobnie, status społeczny (lstat_cat) oraz dostępność autostrad (rad) również znacząco wpływały na wartość mieszkań.

Dodatkowo przeprowadzono test ANOVA z bootstrapem, który potwierdził istotność tych czynników i zapewnił większą wiarygodność wyników, szczególnie w przypadku potencjalnych naruszeń założeń klasycznego testu ANOVA.

Wizualizacja wyników wyraźnie pokazuje różnice między grupami, szczególnie w przypadku liczby pokoi. Wyniki te wskazują na silny wpływ analizowanych czynników na wartość mieszkań.