Zadania

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.

  1. Wygeneruj przykładową populację liczności 25 o takim rozkładzie (funkcja rnorm()).
  2. Wygeneruj 1000 prób o liczności 25 z otrzymanych danych (funkcja replicate()).
  3. 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)
  4. Porównaj średnią wygenerowanej próby wyjściowej ze średnią wektora średnich z prób bootstrapowych. Oceń ich różnicę.
  5. Przy użyciu funkcji boot (pakiet boot) wyznacz błąd standardowy dla średniej.

Zadanie 2. Estymacja przedziałowa średniej metodą bootstrap.

  1. 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.
  2. 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.
#punktowo klasycznie
library(haven)
dane <- read_sav("dane.sav")
dane$Dochod_na_osobe->dochod
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
#przedziałowo klasycznie
mean(dochod)-1.96*se   #dolna krawędź
## [1] 796.4849
mean(dochod)+1.96*se   #górna krawędź
## [1] 820.5794
B=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] 809.0562
sd(boot.data)
## [1] 1046.451
se2<-sd(boot.data)/sqrt(length(boot.data))
se2
## [1] 6.041386
sep2<-se/mean(boot.data)
sep2
## [1] 0.007597209
hist(boot.data)

plot(density(boot.data))

quantile(boot.data,0.975)
##    97.5% 
## 3077.201
quantile(boot.data,0.025)
## 2.5% 
##    0
#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,] 797.2151 820.8973

Średnia: Obie metody dostarczają zbliżonych wyników dla średniej dochodów. Różnica między wartościami (808.53 vs. 803.02) jest niewielka i można ją uznać za nieistotną statystycznie.

Odchylenie standardowe: Różnica między odchyleniami standardowymi (1064.67 vs. 1050.18) także jest marginalna. Wartość bootstrapowa jest nieco niższa, co wynika z losowości resamplingu.

Błąd standardowy: Błędy standardowe są niemal identyczne w obu metodach (6.1466 vs. 6.0629), co potwierdza spójność estymacji.

Relatywny błąd standardowy: Różnica w relatywnych błędach standardowych jest minimalna (0.0076 vs. 0.0077), co oznacza, że względna precyzja obu metod jest porównywalna.

#przedziałowo Efron
dolny<-quantile(boot.data,0.975)
gorny<-quantile(boot.data,0.025)
cbind(dolny,gorny)
##          dolny gorny
## 97.5% 3077.201     0
#przedziałowo Hall
library(boot)

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% 1618.112 -1459.089

Przedziały ufności:

Klasyczny: [796.48; 820.58] Bootstrapowy: [791.14; 814.91] Efrona: [0; 3099.95] Halla: [-1493.90; 1606.05]

Metoda klasyczna i bootstrap klasyczny dostarczają bardzo podobnych wyników, z niewielkimi różnicami w przedziałach. Przedziały uzyskane metodami Efrona i Halla są znacznie bardziej rozbieżne.

Wielkość resamplingu wpływa na precyzję wyników. Większa liczba obserwacji prowadzi do stabilniejszych oszacowań, co skutkuje mniejszym błędem standardowym i bardziej precyzyjnymi przedziałami ufności.

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.002925157    6.147553
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.5, 820.6 )   (796.4, 820.8 )  
## Calculations and Intervals on Original Scale

Widać, że rozkład wyników jest symetryczny i przypomina rozkład normalny, co sugeruje, że średnia dochodów jest stabilna i dobrze oszacowana.

Bootstrapowe wyniki średniej dochodów dobrze odpowiadają normalnemu rozkładowi, co uzasadnia użycie przedziałów ufności opartych na normalności

RAPORT 7

library(MKinfer) 
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
library(tidyverse)
## ── 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.4     ✔ 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
# Przykład 1. 
library(haven)
dane <- read_sav("C:/Users/Beata Wachowicz/Downloads/dane.sav")
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.6707 
## bootstrap difference of means (SE) = 24.38195 (54.17378) 
## 95 percent bootstrap percentile confidence interval:
##  -77.59424 134.61599
## 
## 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
# Czy samochody z automatyczną i manualną skrzynią biegów
# mają istotnie różne spalanie?
#wyniki<-boot.t.test(dochod[,dane$Wojewodztwo==22],
                    # dochod[,dane$Wojewodztwo==16],R=999,data=dane)

#wyniki
  
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.6527 
## bootstrap difference of means (SE) = 21.76314 (53.67489) 
## 95 percent bootstrap percentile confidence interval:
##  -74.95652 134.20318
## 
## 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
# Przykład 2. Test Chi2
library(ISLR)
data("Credit")
attach(Credit)
# Czy status studenta różni się istotnie wg płci?
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.3168
ggbarstats(
  data=Credit,
  x=Gender,
  y=Student
)

ggpiestats(
  data=Credit,
  x=Gender,
  y=Student
)

Wniosek końcowy

Nie stwierdzono istotnej różnicy między płciami w kontekście statusu studenta.

Wynik testu statystycznego potwierdza, że rozkład statusu studenta jest niezależny od płci. Wykresy służą jako dodatkowe narzędzie ułatwiające interpretację wyników.

library(readr)
efektywność_systemy <- read_delim("efektywność_systemy.csv", 
    delim = ";", escape_double = FALSE, trim_ws = TRUE)
## New names:
## Rows: 20 Columns: 3
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: ";" chr
## (1): system dbl (2): ...1, efektywnosc
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`
#import danych efektywnosc_systemy
library(ggstatsplot)

library(readr)
efekt <- efektywność_systemy
dane_1_ <-efektywność_systemy

ggbetweenstats(efektywność_systemy, y=efektywnosc, x=system)

ggbetweenstats(efektywność_systemy, y=efektywnosc, x=system,nboot = 999)

#2. Sprawdzenie założenia o normalności
a=subset(efektywność_systemy$efektywnosc,efektywność_systemy$system=='A')
shapiro.test(a)
## 
##  Shapiro-Wilk normality test
## 
## data:  a
## W = 0.91538, p-value = 0.5006
b=subset(efektywność_systemy$efektywnosc,efektywność_systemy$system=='B')
shapiro.test(b)
## 
##  Shapiro-Wilk normality test
## 
## data:  b
## W = 0.92826, p-value = 0.5846
c=subset(efektywność_systemy$efektywnosc,efektywność_systemy$system=='C')
shapiro.test(c)
## 
##  Shapiro-Wilk normality test
## 
## data:  c
## W = 0.81741, p-value = 0.1115
d=subset(efektywność_systemy$efektywnosc,efektywność_systemy$system=='D')
shapiro.test(d)
## 
##  Shapiro-Wilk normality test
## 
## data:  d
## W = 0.94976, p-value = 0.7355

Wyniki testu Shapiro-Wilka sugerują, że dla wszystkich czterech systemów (A, B, C, D) można przyjąć założenie o normalności rozkładu zmiennej efektywność. Dzięki temu można bezpiecznie stosować testy parametryczne, takie jak ANOVA czy test t-Studenta.

#### 3. Sprawdzenie założenia o jednorodności wariancji, używamy testu Bartlett’a
bartlett.test(efektywność_systemy$efektywnosc ~ system, efekt)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  efektywność_systemy$efektywnosc by system
## Bartlett's K-squared = 0.13684, df = 3, p-value = 0.9871
#### 4. Porównanie średnich – Anova jednoczynnikowa
anova(aov(efektywność_systemy$efektywnosc ~ system, efekt))
## Analysis of Variance Table
## 
## Response: efektywność_systemy$efektywnosc
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## system     3 273.75  91.250  8.5882 0.001257 **
## Residuals 16 170.00  10.625                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hipoteza o jednorodności wariancji została przyjęta (wariancje są jednorodne). Istnieje istotna różnica między średnimi efektywności systemów, co sugeruje, że któryś z systemów jest wyraźnie lepszy lub gorszy w porównaniu do innych.

library(agricolae)
LSD.test(aov(efektywność_systemy$efektywnosc ~ system, efekt), "system",
p.adj="bonferroni", main="Efektywnosc czyszczenia terenu z
wyciekow oleju")

#### 6. Test post-hoc HSD Tukeya

post_model1 <- aov(efektywnosc ~ system, data = efektywność_systemy)
efekt_hsd <- TukeyHSD(post_model1)
print(efekt_hsd)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = efektywnosc ~ system, data = efektywność_systemy)
## 
## $system
##     diff        lwr        upr     p adj
## B-A   -7 -12.898143 -1.1018565 0.0174543
## C-A    3  -2.898143  8.8981435 0.4854618
## D-A   -3  -8.898143  2.8981435 0.4854618
## C-B   10   4.101857 15.8981435 0.0009174
## D-B    4  -1.898143  9.8981435 0.2510126
## D-C   -6 -11.898143 -0.1018565 0.0454680
plot(efekt_hsd)

Interpretacja: Istnieją istotne różnice w efektywności między: systemem B a A (p=0.0175), systemem C a B (p=0.0009), systemem D a C (p=0.0455). Nie stwierdzono istotnych różnic między: systemem C a A (p=0.4855), systemem D a A (p=0.4855), systemem D a B (p=0.2510).

ANOVA

#dodac anove (wedlug wojewodztwa, wydatki, dochod, czy roznie sie isttonie wzgledem wojewodztwa zrobic ten bokstrap z i bez ,test t-student

# ANOVA dla dochodu na osobę
anova_dochod <- aov(Dochod_na_osobe ~ factor(Wojewodztwo), data = dane)
summary(anova_dochod)
##                        Df    Sum Sq  Mean Sq F value Pr(>F)    
## factor(Wojewodztwo)    15 2.293e+08 15288460   13.57 <2e-16 ***
## Residuals           29987 3.378e+10  1126442                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Post-hoc test Tukeya
tukey_dochod <- TukeyHSD(anova_dochod)
print(tukey_dochod)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Dochod_na_osobe ~ factor(Wojewodztwo), data = dane)
## 
## $`factor(Wojewodztwo)`
##              diff         lwr           upr     p adj
## 4-2     -2.138444 -120.321786  116.04489708 1.0000000
## 6-2     -8.675223 -122.370158  105.01971140 1.0000000
## 8-2     50.850169 -104.106909  205.80724632 0.9992145
## 10-2    57.786341  -46.560938  162.13361950 0.8777119
## 12-2     9.903396  -93.551013  113.35780479 1.0000000
## 14-2   234.304768  142.571098  326.03843871 0.0000000
## 16-2   114.027701  -29.053748  257.10914937 0.3107886
## 18-2  -103.576642 -220.763672   13.61038812 0.1570308
## 20-2    29.433834 -112.505424  171.37309221 0.9999978
## 22-2    92.379486  -22.339000  207.09797127 0.2931209
## 24-2    80.435673  -13.167571  174.03891815 0.1934378
## 26-2   -53.087169 -192.692420   86.51808250 0.9958803
## 28-2   -43.347332 -173.930491   87.23582732 0.9991015
## 30-2    70.348011  -33.412997  174.10901805 0.6065380
## 32-2    37.739909  -84.504098  159.98391545 0.9996201
## 6-4     -6.536779 -134.518818  121.44525937 1.0000000
## 8-4     52.988613 -112.735736  218.71296211 0.9994219
## 10-4    59.924785  -59.830042  179.67961189 0.9435404
## 12-4    12.041840 -106.935799  131.01947966 1.0000000
## 14-4   236.443213  127.503182  345.38324315 0.0000000
## 16-4   116.166145  -38.511471  270.84376114 0.4176362
## 18-4  -101.438197 -232.532301   29.65590594 0.3628356
## 20-4    31.572278 -122.049390  185.19394670 0.9999980
## 22-4    94.517930  -34.374253  223.41011324 0.4622419
## 24-4    82.574118  -27.944801  193.09303644 0.4272605
## 26-4   -50.948724 -202.416510  100.51906143 0.9989525
## 28-4   -41.208887 -184.403948  101.98617326 0.9998394
## 30-4    72.486455  -46.757877  191.73078659 0.7766089
## 32-4    39.878353  -95.755202  175.51190769 0.9997906
## 8-6     59.525392 -103.028586  222.07937049 0.9972657
## 10-6    66.461564  -48.866033  181.78916084 0.8397207
## 12-6    18.578619  -95.941749  133.09898735 0.9999999
## 14-6   242.979992  138.926322  347.03366181 0.0000000
## 16-6   122.702924  -28.572986  273.97883407 0.2810663
## 18-6   -94.901418 -221.964001   32.16116446 0.4279144
## 20-6    38.109057 -112.086991  188.30510597 0.9999676
## 22-6   101.054709  -23.734838  225.84425557 0.2837589
## 24-6    89.110897  -16.594672  194.81646622 0.2207087
## 26-6   -44.411945 -192.404263  103.58037257 0.9997308
## 28-6   -34.672108 -174.185768  104.84155198 0.9999753
## 30-6    79.023234  -35.774182  193.82064968 0.5793959
## 32-6    46.415132  -85.325891  178.15615519 0.9982213
## 10-8     6.936172 -149.222761  163.09510444 1.0000000
## 12-8   -40.946773 -196.510497  114.61695148 0.9999480
## 14-8   183.454599   35.426481  331.48271813 0.0023095
## 16-8    63.177532 -121.130793  247.48585667 0.9986979
## 18-8  -154.426811 -319.442136   10.58851479 0.0970754
## 20-8   -21.416335 -204.839372  162.00670259 1.0000000
## 22-8    41.529317 -121.742201  204.80083475 0.9999666
## 24-8    29.585505 -119.608413  178.77942223 0.9999988
## 26-8  -103.937337 -285.560254   77.68557900 0.8469427
## 28-8   -94.197500 -268.980849   80.58584839 0.8994364
## 30-8    19.497842 -136.269948  175.26563145 1.0000000
## 32-8   -13.110260 -181.754427  155.53390611 1.0000000
## 12-10  -47.882945 -153.128995   57.36310543 0.9753182
## 14-10  176.518428   82.768853  270.26800176 0.0000000
## 16-10   56.241360  -88.140832  200.62355187 0.9946785
## 18-10 -161.362982 -280.134681  -42.59128405 0.0003651
## 20-10  -28.352507 -171.602880  114.89786642 0.9999988
## 22-10   34.593145  -81.743639  150.92992896 0.9997594
## 24-10   22.649333  -72.930392  118.22905752 0.9999868
## 26-10 -110.873509 -251.811590   30.06457111 0.3333355
## 28-10 -101.133672 -233.140785   30.87344027 0.3809092
## 30-10   12.561670  -92.985774  118.10911439 1.0000000
## 32-10  -20.046432 -143.810375  103.71751088 0.9999999
## 14-12  224.401372  131.646627  317.15611762 0.0000000
## 16-12  104.124305  -39.613921  247.86253016 0.4851560
## 18-12 -113.480038 -231.468073    4.50799752 0.0748946
## 20-12   19.530438 -123.070858  162.13173367 1.0000000
## 22-12   82.476090  -33.060517  198.01269614 0.5125852
## 24-12   70.532278  -24.071866  165.13642065 0.4312536
## 26-12  -62.990564 -203.268869   77.28773963 0.9781319
## 28-12  -53.250728 -184.553195   78.05173961 0.9919135
## 30-12   60.444615  -44.220199  165.10942799 0.8374862
## 32-12   27.836513  -95.175575  150.84860026 0.9999928
## 16-14 -120.277068 -255.823845   15.26971036 0.1522920
## 18-14 -337.881410 -445.739780 -230.02303983 0.0000000
## 20-14 -204.870934 -339.211476  -70.53039282 0.0000204
## 22-14 -141.925283 -247.096377  -36.75418870 0.0004233
## 24-14 -153.869095 -235.491368  -72.24682170 0.0000000
## 26-14 -287.391937 -419.264050 -155.51982387 0.0000000
## 28-14 -277.652100 -399.932881 -155.37131868 0.0000000
## 30-14 -163.956758 -257.053346  -70.86016955 0.0000002
## 32-14 -196.564860 -309.897226  -83.23249331 0.0000003
## 18-16 -217.604342 -371.522056  -63.68662848 0.0001459
## 20-16  -84.593867 -258.100684   88.91295081 0.9545385
## 22-16  -21.648215 -173.694898  130.39846738 1.0000000
## 24-16  -33.592027 -170.410996  103.22694117 0.9999789
## 26-16 -167.114869 -338.717576    4.48783708 0.0661166
## 28-16 -157.375032 -321.721691    6.97162599 0.0781809
## 30-16  -43.679690 -187.638745  100.27936439 0.9996921
## 32-16  -76.287792 -234.089756   81.51417211 0.9577035
## 20-18  133.010476  -19.846041  285.86699256 0.1766760
## 22-18  195.956127   67.976861  323.93539388 0.0000182
## 24-18  184.012315   74.559453  293.46517778 0.0000010
## 26-18   50.489473 -100.202225  201.18117114 0.9989989
## 28-18   60.229310  -82.144575  202.60319505 0.9876511
## 30-18  173.924652   55.667693  292.18161145 0.0000543
## 32-18  141.316550    6.550238  276.08286260 0.0288249
## 22-20   62.945652  -88.026683  213.91798595 0.9893076
## 24-20   51.001839  -84.622208  186.62588676 0.9963460
## 26-20  -82.521003 -253.172521   88.13051637 0.9576102
## 28-20  -72.781166 -236.134391   90.57205959 0.9796708
## 30-20   40.914177 -101.909706  183.73805925 0.9998483
## 32-20    8.306075 -148.460988  165.07313746 1.0000000
## 24-22  -11.943812 -118.749524   94.86189998 1.0000000
## 26-22 -145.466654 -294.246756    3.31344818 0.0634797
## 28-22 -135.726817 -276.075861    4.62222670 0.0709426
## 30-22  -22.031475 -137.842698   93.77974809 0.9999993
## 32-22  -54.639577 -187.264951   77.98579741 0.9905148
## 26-24 -133.522842 -266.702252   -0.34343185 0.0486096
## 28-24 -123.783005 -247.472498   -0.09351252 0.0495889
## 30-24  -10.087663 -105.026990   84.85166428 1.0000000
## 32-24  -42.695765 -157.546653   72.15512274 0.9967751
## 28-26    9.739837 -151.589482  171.06915630 1.0000000
## 30-26  123.435179  -17.069392  263.93975046 0.1644975
## 32-26   90.827077  -63.829915  245.48406891 0.8189919
## 30-28  113.695342  -17.848832  245.23951674 0.1855920
## 32-28   81.087240  -65.477150  227.65163014 0.8785417
## 32-30  -32.608102 -155.878154   90.66194965 0.9999446
plot(tukey_dochod)

# ANOVA dla całkowitych wydatków
anova_wydatki <- aov(Calkowite_wydatki_GD_GUS ~ factor(Wojewodztwo), data = dane)
summary(anova_wydatki)
##                        Df    Sum Sq  Mean Sq F value Pr(>F)    
## factor(Wojewodztwo)    15 8.940e+08 59598753   27.12 <2e-16 ***
## Residuals           29987 6.589e+10  2197245                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Post-hoc test Tukeya
tukey_wydatki <- TukeyHSD(anova_wydatki)
print(tukey_wydatki)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Calkowite_wydatki_GD_GUS ~ factor(Wojewodztwo), data = dane)
## 
## $`factor(Wojewodztwo)`
##               diff         lwr          upr     p adj
## 4-2    -51.2239112 -216.283628  113.8358051 0.9995955
## 6-2     59.1938515  -99.597172  217.9848745 0.9966766
## 8-2    107.2378679 -109.181561  323.6572967 0.9479291
## 10-2   100.4253717  -45.310330  246.1610738 0.5775058
## 12-2   196.8899984   52.401315  341.3786816 0.0003417
## 14-2   429.9508761  301.831858  558.0698939 0.0000000
## 16-2   493.5855187  293.752082  693.4189558 0.0000000
## 18-2   -35.1483866 -198.816613  128.5198399 0.9999965
## 20-2   -83.6053763 -281.843583  114.6328306 0.9880188
## 22-2   181.0246157   20.804059  341.2451719 0.0105262
## 24-2    51.7664774  -78.963665  182.4966195 0.9936799
## 26-2   -92.8691959 -287.847632  102.1092406 0.9628020
## 28-2  -136.5981440 -318.975956   45.7796676 0.4226888
## 30-2   240.8832810   95.966390  385.8001724 0.0000015
## 32-2     0.1377959 -170.593213  170.8688046 1.0000000
## 6-4    110.4177627  -68.327217  289.1627422 0.7548309
## 8-4    158.4617791  -72.995664  389.9192221 0.5892005
## 10-4   151.6492829  -15.605234  318.9038003 0.1279365
## 12-4   248.1139096   81.944844  414.2829750 0.0000364
## 14-4   481.1747873  329.024659  633.3249153 0.0000000
## 16-4   544.8094299  328.780309  760.8385513 0.0000000
## 18-4    16.0755246 -167.015893  199.1669419 1.0000000
## 20-4   -32.3814651 -246.935807  182.1728762 1.0000000
## 22-4   232.2485269   52.232402  412.2646522 0.0010773
## 24-4   102.9903886  -51.364879  257.3456566 0.6343565
## 26-4   -41.6452847 -253.191425  169.9008555 0.9999989
## 28-4   -85.3742328 -285.366345  114.6178792 0.9865027
## 30-4   292.1071922  125.565653  458.6487310 0.0000002
## 32-4    51.3617071 -138.069695  240.7931091 0.9999241
## 8-6     48.0440164 -178.985556  275.0735892 0.9999971
## 10-6    41.2315202 -119.839745  202.3027860 0.9999637
## 12-6   137.6961469  -22.247710  297.6400037 0.1909317
## 14-6   370.7570246  225.431389  516.0826606 0.0000000
## 16-6   434.3916673  223.113509  645.6698259 0.0000000
## 18-6   -94.3422381 -271.803068   83.1185918 0.9091536
## 20-6  -142.7992278 -352.569208   66.9707520 0.5993452
## 22-6   121.8307642  -52.455450  296.1169783 0.5514352
## 24-6    -7.4273741 -155.060120  140.2053723 1.0000000
## 26-6  -152.0630474 -358.755206   54.6291113 0.4561466
## 28-6  -195.7919955 -390.642512   -0.9414786 0.0474219
## 30-6   181.6894295   21.358637  342.0202225 0.0100728
## 32-6   -59.0560556 -243.050988  124.9388767 0.9993951
## 10-8    -6.8124962 -224.910485  211.2854928 1.0000000
## 12-8    89.6521305 -127.614566  306.9188272 0.9903625
## 14-8   322.7130082  115.970848  529.4551680 0.0000105
## 16-8   386.3476508  128.935063  643.7602389 0.0000317
## 18-8  -142.3862545 -372.853446   88.0809369 0.7546582
## 20-8  -190.8432443 -447.019404   65.3329151 0.4327138
## 22-8    73.7867478 -154.244971  301.8184662 0.9993341
## 24-8   -55.4713905 -263.841753  152.8989720 0.9999399
## 26-8  -200.1070638 -453.769100   53.5549724 0.3284519
## 28-8  -243.8360119 -487.945627    0.2736030 0.0506148
## 30-8   133.6454131  -83.906290  351.1971158 0.7624966
## 32-8  -107.1000720 -342.635452  128.4353081 0.9754438
## 12-10   96.4646267  -50.526336  243.4555897 0.6623010
## 14-10  329.5255044  198.590993  460.4600161 0.0000000
## 16-10  393.1601470  191.510039  594.8102554 0.0000000
## 18-10 -135.5737583 -301.455198   30.3076813 0.2686675
## 20-10 -184.0307481 -384.100112   16.0386157 0.1131907
## 22-10   80.5992440  -81.881494  243.0799822 0.9474525
## 24-10  -48.6588943 -182.149470   84.8316810 0.9974005
## 26-10 -193.2945676 -390.134488    3.5453530 0.0606881
## 28-10 -237.0235157 -421.390079  -52.6569520 0.0011550
## 30-10  140.4579093   -6.953993  287.8698120 0.0820304
## 32-10 -100.2875758 -273.141390   72.5662381 0.8325222
## 14-12  233.0608777  103.515785  362.6059705 0.0000001
## 16-12  296.6955204   95.944802  497.4462387 0.0000478
## 18-12 -232.0383850 -396.825328  -67.2514415 0.0001606
## 20-12 -280.4953747 -479.658210  -81.3325393 0.0001599
## 22-12  -15.8653827 -177.228560  145.4977947 1.0000000
## 24-12 -145.1235209 -277.251559  -12.9954831 0.0157496
## 26-12 -289.7591943 -485.677644  -93.8407443 0.0000470
## 28-12 -333.4881424 -516.870569 -150.1057154 0.0000001
## 30-12   43.9932826 -102.185901  190.1724664 0.9997212
## 32-12 -196.7522025 -368.555944  -24.9484607 0.0085991
## 16-14   63.6346426 -125.675563  252.9448485 0.9989606
## 18-14 -465.0992627 -615.738699 -314.4598261 0.0000000
## 20-14 -513.5562525 -701.181779 -325.9307260 0.0000000
## 22-14 -248.9262604 -395.812537 -102.0399842 0.0000008
## 24-14 -378.1843987 -492.181423 -264.1873745 0.0000000
## 26-14 -522.8200720 -706.998090 -338.6420544 0.0000000
## 28-14 -566.5490201 -737.331389 -395.7666509 0.0000000
## 30-14 -189.0675951 -319.090120  -59.0450705 0.0000728
## 32-14 -429.8130802 -588.097725 -271.5284353 0.0000000
## 18-16 -528.7339053 -743.701716 -313.7660945 0.0000000
## 20-16 -577.1908951 -819.517654 -334.8641363 0.0000000
## 22-16 -312.5609030 -524.915554 -100.2062516 0.0000533
## 24-16 -441.8190413 -632.906041 -250.7320418 0.0000000
## 26-16 -586.4547146 -826.122113 -346.7873159 0.0000000
## 28-16 -630.1836627 -859.716966 -400.6503591 0.0000000
## 30-16 -252.7022377 -453.761375  -51.6431003 0.0017826
## 32-16 -493.4477228 -713.840437 -273.0550082 0.0000000
## 20-18  -48.4569898 -261.942689  165.0287096 0.9999925
## 22-18  216.1730023   37.431894  394.9141106 0.0035562
## 24-18   86.9148640  -65.951506  239.7812344 0.8534342
## 26-18  -57.7208093 -268.183034  152.7414149 0.9999122
## 28-18 -101.4497574 -300.294982   97.3954671 0.9339533
## 30-18  276.0316676  110.869134  441.1942013 0.0000012
## 32-18   35.2861825 -152.933993  223.5063580 0.9999994
## 22-20  264.6299921   53.775819  475.4841646 0.0018314
## 24-20  135.3718538  -54.046270  324.7899772 0.5104454
## 26-20   -9.2638196 -247.602751  229.0751117 1.0000000
## 28-20  -52.9927677 -281.138602  175.1530664 0.9999899
## 30-20  324.4886574  125.014948  523.9623670 0.0000030
## 32-20   83.7431723 -135.204157  302.6905011 0.9956161
## 24-22 -129.2581383 -278.427389   19.9111128 0.1821930
## 26-22 -273.8938116 -481.686223  -66.1014005 0.0006977
## 28-22 -317.6227597 -513.640008 -121.6055116 0.0000034
## 30-22   59.8586653 -101.888053  221.6053834 0.9969299
## 32-22 -180.8868198 -366.116873    4.3432338 0.0642698
## 26-24 -144.6356733 -330.639516   41.3681694 0.3539879
## 28-24 -188.3646214 -361.114455  -15.6147879 0.0173965
## 30-24  189.1168036   56.520634  321.7129733 0.0001179
## 32-24  -51.6286815 -212.034156  108.7767927 0.9993747
## 28-26  -43.7289481 -269.048112  181.5902158 0.9999991
## 30-26  333.7524769  137.518013  529.9869410 0.0000007
## 32-26   93.0069918 -122.993325  309.0073085 0.9853205
## 30-28  377.4814250  193.761420  561.2014303 0.0000000
## 32-28  136.7359399  -67.961916  341.4337959 0.6324031
## 32-30 -240.7454851 -412.909510  -68.5814601 0.0001893
plot(tukey_wydatki)

T-student

# Filtracja danych dla dwóch województw
woj10 <- subset(dane$Dochod_na_osobe, dane$Wojewodztwo == 10)
woj12 <- subset(dane$Dochod_na_osobe, dane$Wojewodztwo == 12)

# Test T-Studenta
t_test_dochod <- t.test(woj10, woj12)
print(t_test_dochod)
## 
##  Welch Two Sample t-test
## 
## data:  woj10 and woj12
## t = 1.7642, df = 4616.9, p-value = 0.07777
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   -5.328738 101.094628
## sample estimates:
## mean of x mean of y 
##  809.5731  761.6902

Interpretacja: Test t-Studenta sugeruje, że różnica w średnich dochodach między województwem 10 a województwem 12 nie jest istotna statystycznie, ponieważ wartość p wynosi 0.07777, co jest większe niż 0.05. Choć średni dochód w województwie 10 (809.57) jest wyższy niż w województwie 12 (761.69), różnice te nie są wystarczająco silne, aby można było stwierdzić ich istotność statystyczną na poziomie 5%.

Podsumowanie: Na podstawie tego testu nie możemy odrzucić hipotezy zerowej, więc nie stwierdzamy, że średnie dochody w tych dwóch województwach różnią się istotnie.

Bootstrap: Estymacja średniej i błędu standardowego

library(boot)

# Funkcja dla średniej
boot_mean <- function(data, indices) {
  return(mean(data[indices]))
}

# Bootstrap dla dochodu
boot_dochod <- boot(data = dane$Dochod_na_osobe, statistic = boot_mean, R = 1000)

# Wyniki
print(boot_dochod)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = dane$Dochod_na_osobe, statistic = boot_mean, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original     bias    std. error
## t1* 808.5322 0.03051041    6.033991
plot(boot_dochod)

# Przedziały ufności
boot_ci_dochod <- boot.ci(boot_dochod, conf = 0.95, type = "perc")
print(boot_ci_dochod)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_dochod, conf = 0.95, type = "perc")
## 
## Intervals : 
## Level     Percentile     
## 95%   (797.2, 819.5 )  
## Calculations and Intervals on Original Scale
# Bootstrap dla wydatków
boot_wydatki <- boot(data = dane$Calkowite_wydatki_GD_GUS, statistic = boot_mean, R = 1000)

# Wyniki
print(boot_wydatki)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = dane$Calkowite_wydatki_GD_GUS, statistic = boot_mean, 
##     R = 1000)
## 
## 
## Bootstrap Statistics :
##     original      bias    std. error
## t1* 2068.324 -0.03149227    8.449473
plot(boot_wydatki)

# Przedziały ufności
boot_ci_wydatki <- boot.ci(boot_wydatki, conf = 0.95, type = "perc")
print(boot_ci_wydatki)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_wydatki, conf = 0.95, type = "perc")
## 
## Intervals : 
## Level     Percentile     
## 95%   (2052, 2086 )  
## Calculations and Intervals on Original Scale

Z wyników ANOVA widać, że różnice między systemami są istotne statystycznie (p < 0.05). Testy Tukeya wskazują, które grupy mają istotne różnice w średnich, np. różnice między niektórymi województwami.