Základní statistické modely a testy
Karel Rečka
31. 10. 2022
1 Balíčky
Nové balíčky:
- psych je obecný balíček pro analýzu dat v psychologii.
- broom je balíček složící k extrakci informací z odhadnutých modelů do přehledných tabulek.
- car poskytuje mnoho užitečných funkcí pro diagnostiku regresních modelů.
- effectsize poskytuje odhady standardizovaných velikostí účinků.
- corrplot slouží k vytváření korelogramů.
- rstatix přepracovává mnohé funkce z base R pro odhad modelů tak, aby poskytovali přehlednější output a byly “pipable”.
# Instalace balíčků, pokud nejsou instalovány
pckg <- c("psych", "broom", "car", "effectsize",
"corrplot", "rstatix", "lm.beta")
pckg <- pckg[!(pckg %in% installed.packages()[,"Package"])]
if(length(pckg)) install.packages(pckg)library(tidyverse)
library(ggmosaic)
source("https://is.muni.cz/el/fss/podzim2022/PSYn5320/um/r101_funs.R")2 Import dat
Pro cvičení použijeme opět data z Cest do dospělosti.
# Import dat
df <- readRDS(
url("https://is.muni.cz/el/fss/podzim2022/PSYn5320/um/datasets/eukids.rds",
"rb"))
# Přehled dat
glimpse(df)
tibble(
var_name = names(attributes(df)$variable.labels),
var_label = attributes(df)$variable.labels
)3 Korelace
Nejprve se podíváme na obyčejné bivariační korelace. Vytvoříme si
vektor s proměnnými, pro které chceme vypočíst korelace, např.
optimismus, životní spokojenost, psychosomatické symptoma a depresivita
a z datasetu tyto proměnné vybereme (vytvoříme si nový dataset
df_cor).
variables <- c("optim", "ziv_sp", "selfe", "effi", "zdravi", "deprese")
df_cor <- df %>% select(all_of(variables))
df_cor># # A tibble: 768 × 6
># optim ziv_sp selfe effi zdravi deprese
># <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
># 1 18 3.67 2.75 3 30 2.1
># 2 12 2.83 3.12 2 19 1.65
># 3 18 2.83 3 3 19 2.25
># 4 NA 2.67 2.67 3 34 2.06
># 5 20 3.67 4 3.25 19 1.25
># 6 17 3.33 3 1.67 NA 2.21
># 7 17 2.83 2.75 2 32 2.05
># 8 20 3.17 3.12 2.75 27 1.65
># 9 NA 2.67 2.75 2 26 1.45
># 10 13 3 3.12 2.75 16 1.16
># # … with 758 more rowsMůžeme si také vytvořit tabulku se všemi možnými páry proměných
pomocí funkce combn().
cor_pairs <- combn(variables, 2) %>% # Tento řádek vytvoří matici se dvěma řádky
t() %>% # Tímto ji transponujeme
as_tibble(.names_repair = TRUE) %>% # Konvertujeme na tibble()
setNames(c("x", "y")) # Nastavíme jméma sloupců
cor_pairs># # A tibble: 15 × 2
># x y
># <chr> <chr>
># 1 optim ziv_sp
># 2 optim selfe
># 3 optim effi
># 4 optim zdravi
># 5 optim deprese
># 6 ziv_sp selfe
># 7 ziv_sp effi
># 8 ziv_sp zdravi
># 9 ziv_sp deprese
># 10 selfe effi
># 11 selfe zdravi
># 12 selfe deprese
># 13 effi zdravi
># 14 effi deprese
># 15 zdravi depreseVždy je vhodné se podívat na univariační rozdělení proměnných, které chceme korelovat. Pokud není rozdělení univariačně normální, nemůže být splněna ani bivariační normalita. Můžeme si také všímat efektů podlahy/stropu, míry zešikmení, variability apod.
Skript r101_funs.R obsahuje funkci
plot_histogram_mult(), do které stačí zadat použitý dataset
a character vektor s názvy sloupců. Vytvoří se nám histogramy všech
proměnných.
# Histogramy
plot_histogram_mult(df_cor, var = variables)
Vždy je také vhodné se podívat na scatterploty analyzovaných vztahů,
abychom mohli ověřit předpoklady včetně linearity, bivariační normality
či absence extrémních případů. To můžeme udělat pomocí funkce
plot_scatter_mult().
# Scatterploty
plot_scatter_mult(df_cor, # dataset
xvar = cor_pairs$x, # character vektor s proměnnými na ose X
yvar = cor_pairs$y, # character vektor s proměnnými na ose Y
geom = "jitter", # geom jitter místo standardního point
alpha = 0.25) # průhlednost bodů
Samotné korelace můžeme vypočíst pomocí funkce cor() z base
R. Defaultně tato funkce vrací NA, pokud se v daném páru
proměnných nachází byť jen jedna chybějící hodnota. Argumentem
use to můžeme změnit. Když nastavíme
use = "pairwise", využijí se pro výpočet korelace všechny
případy, které mají pro daný pár proměnných validní hodnoty (pairwise
deletion). Pokud nastavíme use = "pairwise", použijí se
pouze případy, které mají ve všech korelovaných proměnných platné
hodnoty (listwise deletion)
cor(df_cor) # defaultně vyřadí proměnné s chybějícími hodnotami># optim ziv_sp selfe effi zdravi deprese
># optim 1 NA NA NA NA NA
># ziv_sp NA 1 NA NA NA NA
># selfe NA NA 1 NA NA NA
># effi NA NA NA 1 NA NA
># zdravi NA NA NA NA 1 NA
># deprese NA NA NA NA NA 1cor(df_cor, use = "pairwise")># optim ziv_sp selfe effi zdravi deprese
># optim 1.0000000 0.5443483 0.4618558 0.3762663 -0.1341719 -0.3383751
># ziv_sp 0.5443483 1.0000000 0.6445489 0.4000271 -0.1622375 -0.5039272
># selfe 0.4618558 0.6445489 1.0000000 0.4160053 -0.2156670 -0.5166076
># effi 0.3762663 0.4000271 0.4160053 1.0000000 -0.1463420 -0.2819470
># zdravi -0.1341719 -0.1622375 -0.2156670 -0.1463420 1.0000000 0.3993650
># deprese -0.3383751 -0.5039272 -0.5166076 -0.2819470 0.3993650 1.0000000cor(df_cor, use = "complete")># optim ziv_sp selfe effi zdravi deprese
># optim 1.0000000 0.5245830 0.4435352 0.3672163 -0.1308969 -0.3257702
># ziv_sp 0.5245830 1.0000000 0.6270288 0.3925454 -0.1620275 -0.4964127
># selfe 0.4435352 0.6270288 1.0000000 0.4207780 -0.2156612 -0.5175092
># effi 0.3672163 0.3925454 0.4207780 1.0000000 -0.1614927 -0.2676479
># zdravi -0.1308969 -0.1620275 -0.2156612 -0.1614927 1.0000000 0.4101335
># deprese -0.3257702 -0.4964127 -0.5175092 -0.2676479 0.4101335 1.0000000Defaultně funkce cor() počítá Pearsonovy korelace, ale
argumentem method si můžeme vyžádat Spearmanovo rhó nebo
Kendallovo tau-b.
cor(df_cor, use = "pairwise") %>%
round(digits = 2) ># optim ziv_sp selfe effi zdravi deprese
># optim 1.00 0.54 0.46 0.38 -0.13 -0.34
># ziv_sp 0.54 1.00 0.64 0.40 -0.16 -0.50
># selfe 0.46 0.64 1.00 0.42 -0.22 -0.52
># effi 0.38 0.40 0.42 1.00 -0.15 -0.28
># zdravi -0.13 -0.16 -0.22 -0.15 1.00 0.40
># deprese -0.34 -0.50 -0.52 -0.28 0.40 1.00cor(df_cor, use = "pairwise", method = "kendall") %>%
round(digits = 2)># optim ziv_sp selfe effi zdravi deprese
># optim 1.00 0.39 0.33 0.29 -0.10 -0.22
># ziv_sp 0.39 1.00 0.48 0.28 -0.13 -0.36
># selfe 0.33 0.48 1.00 0.32 -0.16 -0.36
># effi 0.29 0.28 0.32 1.00 -0.10 -0.19
># zdravi -0.10 -0.13 -0.16 -0.10 1.00 0.30
># deprese -0.22 -0.36 -0.36 -0.19 0.30 1.00cor(df_cor, use = "pairwise", method = "spearman") %>%
round(digits = 2)># optim ziv_sp selfe effi zdravi deprese
># optim 1.00 0.51 0.43 0.37 -0.14 -0.30
># ziv_sp 0.51 1.00 0.62 0.37 -0.18 -0.50
># selfe 0.43 0.62 1.00 0.41 -0.22 -0.50
># effi 0.37 0.37 0.41 1.00 -0.14 -0.26
># zdravi -0.14 -0.18 -0.22 -0.14 1.00 0.42
># deprese -0.30 -0.50 -0.50 -0.26 0.42 1.00Bohužel funkce cor() neposkytuje intervaly spolehlivosti
ani p-hodnoty. Funkce cor.test() je poskytuje, ale je
poněkud těžkopádná, protože umožňuje vkládat pouze dva numerické vektory
zároveň. Funkce corr.test() z balíčku
psych je v tomto ohledu uživatelsky
přívětivější. Ukáže nám i velikosti vzorku pro jednotlivé korelace
(pokud zvolíme pairwise deletion), p-hodnoty i intervaly spolehlivosti.
Navíc umožňuje specifikovat typ korekce p-hodnot (např. známou
Bonferroniho korekci či Holmuvu korekci)
r <- psych::corr.test(df_cor,
use = "pairwise", # Způsob zacházení s chybějícími daty
adjust = "holm", # jakou korekci p-hodnot použít
method = "spearman", # typ korelace
alpha = 0.005) # hladina spolehlivosti bude 1 - alpha
print(r, short = FALSE)># Call:psych::corr.test(x = df_cor, use = "pairwise", method = "spearman",
># adjust = "holm", alpha = 0.005)
># Correlation matrix
># optim ziv_sp selfe effi zdravi deprese
># optim 1.00 0.51 0.43 0.37 -0.14 -0.30
># ziv_sp 0.51 1.00 0.62 0.37 -0.18 -0.50
># selfe 0.43 0.62 1.00 0.41 -0.22 -0.50
># effi 0.37 0.37 0.41 1.00 -0.14 -0.26
># zdravi -0.14 -0.18 -0.22 -0.14 1.00 0.42
># deprese -0.30 -0.50 -0.50 -0.26 0.42 1.00
># Sample Size
># optim ziv_sp selfe effi zdravi deprese
># optim 715 714 714 714 663 708
># ziv_sp 714 762 761 758 705 756
># selfe 714 761 761 758 704 755
># effi 714 758 758 758 703 752
># zdravi 663 705 704 703 705 701
># deprese 708 756 755 752 701 756
># Probability values (Entries above the diagonal are adjusted for multiple tests.)
># optim ziv_sp selfe effi zdravi deprese
># optim 0 0 0 0 0 0
># ziv_sp 0 0 0 0 0 0
># selfe 0 0 0 0 0 0
># effi 0 0 0 0 0 0
># zdravi 0 0 0 0 0 0
># deprese 0 0 0 0 0 0
>#
># Confidence intervals based upon normal theory. To get bootstrapped values, try cor.ci
># raw.lower raw.r raw.upper raw.p lower.adj upper.adj
># optim-zv_sp 0.43 0.51 0.58 0 0.40 0.60
># optim-selfe 0.34 0.43 0.51 0 0.32 0.53
># optim-effi 0.28 0.37 0.46 0 0.26 0.48
># optim-zdrav -0.24 -0.14 -0.03 0 -0.24 -0.03
># optim-deprs -0.39 -0.30 -0.20 0 -0.41 -0.18
># zv_sp-selfe 0.56 0.62 0.68 0 0.54 0.70
># zv_sp-effi 0.28 0.37 0.46 0 0.26 0.47
># zv_sp-zdrav -0.28 -0.18 -0.08 0 -0.30 -0.07
># zv_sp-deprs -0.57 -0.50 -0.42 0 -0.59 -0.39
># selfe-effi 0.33 0.41 0.49 0 0.30 0.51
># selfe-zdrav -0.31 -0.22 -0.11 0 -0.33 -0.10
># selfe-deprs -0.57 -0.50 -0.42 0 -0.59 -0.39
># effi-zdrav -0.24 -0.14 -0.04 0 -0.25 -0.03
># effi-deprs -0.35 -0.26 -0.16 0 -0.37 -0.14
># zdrav-deprs 0.33 0.42 0.51 0 0.31 0.52Custom funkce cor_table() z r101_funs.R pak
produkuje tabulky, které by měly odpovídat formátu APA i s průměry,
směrodatnými odhchylkami a velikostí vzorku (napravo od diagonály).
cor_table(df_cor, # dataset
adjust = "none", # jakou korekci p-hodnot použít
conf = .995, # hladina spolehlivosti
sig_levels = c(0.005, 0.001, 0.0001)) # hvězdičkami značené p-hodnoty># var n M SD optim ziv_sp selfe
># 1 optim 715 16.31 3.12 714 714
># 2
># 3 ziv_sp 762 2.91 0.48 0.54*** 761
># 4 [ 0.47; 0.61]
># 5 selfe 761 3.06 0.48 0.46*** 0.64***
># 6 [ 0.38; 0.54] [ 0.58; 0.70]
># 7 effi 758 2.72 0.54 0.38*** 0.40*** 0.42***
># 8 [ 0.28; 0.46] [ 0.31; 0.48] [ 0.33; 0.50]
># 9 zdravi 705 25.43 4.82 -0.13** -0.16*** -0.22***
># 10 [-0.24;-0.03] [-0.26;-0.06] [-0.31;-0.11]
># 11 deprese 756 2.01 0.47 -0.34*** -0.50*** -0.52***
># 12 [-0.43;-0.24] [-0.58;-0.42] [-0.59;-0.44]
># effi zdravi deprese
># 1 714 663 708
># 2
># 3 758 705 756
># 4
># 5 758 704 755
># 6
># 7 703 752
># 8
># 9 -0.15*** 701
># 10 [-0.25;-0.04]
># 11 -0.28*** 0.40***
># 12 [-0.37;-0.19] [ 0.31; 0.48]Když pouze explorujeme data anebo a chceme si udělat rychlý přehled o
bivariačních vztazích mezi kvantitativními proměnnými, můžeme použít
balíček corrplot a jeho funkci
corrplot.mixed(), která nám vytvoří korelogram proměnných.
Pomocí argumentu order = "hclust" si můžeme nechat proměnné
seřadit na základě hierarchické klustrové analýzy. To jednoduše znamená,
že proměnné, které vykazují podobný vzorec korelací s jinými proměnnými,
budou umístěny blíže u sebe (čím jsou si proměnné podobnější z hlediska
vzorce korelací s jinými proměnnými, tím blíže u sebe budou
umístěny)
r <- cor(df_cor, use = "pairwise")
r # korelační matice jako input do funkce corrplot.mixed()># optim ziv_sp selfe effi zdravi deprese
># optim 1.0000000 0.5443483 0.4618558 0.3762663 -0.1341719 -0.3383751
># ziv_sp 0.5443483 1.0000000 0.6445489 0.4000271 -0.1622375 -0.5039272
># selfe 0.4618558 0.6445489 1.0000000 0.4160053 -0.2156670 -0.5166076
># effi 0.3762663 0.4000271 0.4160053 1.0000000 -0.1463420 -0.2819470
># zdravi -0.1341719 -0.1622375 -0.2156670 -0.1463420 1.0000000 0.3993650
># deprese -0.3383751 -0.5039272 -0.5166076 -0.2819470 0.3993650 1.0000000# Defaultní korelogram
corrplot::corrplot.mixed(r)
corrplot::corrplot.mixed(r,
lower = "number", # Hodnoty korelací dole
upper = "ellipse", # Nahoře elipsy
lower.col = "black", # Čísla dole černě
order = "hclust") # Seřadit na základě klustrové analýzy
4 Lineární regrese
V rámci cvičné analýzy budeme predikovat depresivitu na základě pohlaví, věku, průměru známek, důvěry v rodiče, důvěry ve vrstevníky, sebehodnocení a optimismu.
V rámci seriózní analýzy bychom se chybějícím datům věnovali podrobněji, ale protože tohle je jen cvičení, uplatníme vždy listwise deletion a vybereme jen případy s platnými hodnotami ve všech analyzovaných proměnných.
df_clean <- df %>%
select(id, skola, trida, deprese, pohlavi, vekr,
gpa, duv_v, duv_r, selfe, optim) %>%
drop_na()Opět se můžeme podívat na korelace, abychom zjistili, které kvantitativní prediktory souvisejí se závislou proměnnou nejsilněji a zda prediktory nekorelují mezi sebou příliš silně.
df_clean %>%
select(where(is.numeric), -c(id:trida)) %>%
cor() %>%
corrplot::corrplot.mixed(lower = "number", # Hodnoty korelací dole
upper = "ellipse", # Nahoře elipsy
lower.col = "black") # Seřadit na základě klustrové analýzy
Podívat se můžeme také na vztah kategorického prediktoru a závislé
proměnné, například pomocí boxplotu.
df_clean %>%
ggplot(aes(pohlavi, deprese)) +
geom_boxplot()
Pro odhad lineárně regresního modelu se používá funkce lm()
(lm jako linear model). Prvním argumentem je regresní rovnice ve formátu
závislá_proměnná ~ prediktor1 + prediktor2 atd.. Použitý
dataset specifikujeme argumentem data. Do prvního modelu
vložíme jako prediktory pouze pohlaví a věk
lmfit1 <- lm(deprese ~ pohlavi + vekr, data = df_clean) Když vzniklý objekt lmfit1 vyvoláme, ukáže se nám jen
přehled vstupních argumentů (regresní rovnice a dataset) a
b-koeficienty.
lmfit1>#
># Call:
># lm(formula = deprese ~ pohlavi + vekr, data = df_clean)
>#
># Coefficients:
># (Intercept) pohlavizenske vekr
># 1.72601 0.13663 0.01461Celkový přehled modelu získáme tím, že na nově vytvořený objekt
použijeme funkci summary()
summary(lmfit1)>#
># Call:
># lm(formula = deprese ~ pohlavi + vekr, data = df_clean)
>#
># Residuals:
># Min 1Q Median 3Q Max
># -0.99768 -0.34822 -0.04824 0.32199 1.36085
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 1.726013 0.126252 13.671 < 2e-16 ***
># pohlavizenske 0.136627 0.035471 3.852 0.000128 ***
># vekr 0.014608 0.008722 1.675 0.094397 .
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
># Residual standard error: 0.4613 on 699 degrees of freedom
># Multiple R-squared: 0.02457, Adjusted R-squared: 0.02178
># F-statistic: 8.804 on 2 and 699 DF, p-value: 0.0001675To nám ukáže (1) kvantily reziduí; (2) regresní koeficienty společně
s jejich standardními chybami, statistickým testem těchto koeficientů
(t-value) a p-hodnotami (Pr(>|t|)); (3)
informace o fitu modelu včetně standardní odchylky (chyby) reziduí,
podílu vysvětleného rozptylu R2 a
adjustovaného R2 a F-testu celkového modelu
(nulovo hypotézou je, že v populaci R2 = 0).
Také si můžeme všimnout, že R automaticky převede kategorické proměnné
(typu faktor) na dummy proměnné.
Rovnou si můžeme odhadnout druhý model, do kterého doplníme ostatní
prediktory: průměr známek (gpa), důvěru ve vrstevníky a
rodiče (duv_v, duv_r), sebehodnocení
(selfe) a optimismus (optim).
lmfit2 <- lm(deprese ~ pohlavi + vekr +
gpa + duv_v + duv_r + selfe + optim,
data = df_clean)
lmfit2>#
># Call:
># lm(formula = deprese ~ pohlavi + vekr + gpa + duv_v + duv_r +
># selfe + optim, data = df_clean)
>#
># Coefficients:
># (Intercept) pohlavizenske vekr gpa duv_v
># 3.295377 0.140132 0.003581 0.099101 0.037321
># duv_r selfe optim
># -0.049454 -0.441584 -0.015016Mohli bychom to udělat i pomocí funkce update(),
kdybychom chtěli zadat jenom nové prediktory
lmfit2 <- update(lmfit1, # Původní model
. ~ . + gpa + duv_v + duv_r + selfe + optim)
lmfit2>#
># Call:
># lm(formula = deprese ~ pohlavi + vekr + gpa + duv_v + duv_r +
># selfe + optim, data = df_clean)
>#
># Coefficients:
># (Intercept) pohlavizenske vekr gpa duv_v
># 3.295377 0.140132 0.003581 0.099101 0.037321
># duv_r selfe optim
># -0.049454 -0.441584 -0.015016summary(lmfit2)>#
># Call:
># lm(formula = deprese ~ pohlavi + vekr + gpa + duv_v + duv_r +
># selfe + optim, data = df_clean)
>#
># Residuals:
># Min 1Q Median 3Q Max
># -1.0601 -0.2831 -0.0211 0.2455 1.2256
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 3.295377 0.172021 19.157 < 2e-16 ***
># pohlavizenske 0.140132 0.031484 4.451 9.96e-06 ***
># vekr 0.003581 0.008131 0.440 0.65980
># gpa 0.099101 0.020944 4.732 2.70e-06 ***
># duv_v 0.037321 0.030299 1.232 0.21846
># duv_r -0.049454 0.028297 -1.748 0.08097 .
># selfe -0.441584 0.034779 -12.697 < 2e-16 ***
># optim -0.015016 0.005248 -2.861 0.00435 **
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
># Residual standard error: 0.3794 on 694 degrees of freedom
># Multiple R-squared: 0.3449, Adjusted R-squared: 0.3383
># F-statistic: 52.21 on 7 and 694 DF, p-value: < 2.2e-16Spočítat přírůstek vysvětleného rozptylu můžeme pomocí funkce
r2_diff() z r101_funs.R
r2_diff(lmfit1, lmfit2)># [1] 0.320373Statistické srovnání nested modelů umožňuje funkce
anova(). Nulovou hypotézou je v tomto případě nulový
přírůstek vysvětleného rozptylu v populaci.
anova(lmfit1, lmfit2)># Analysis of Variance Table
>#
># Model 1: deprese ~ pohlavi + vekr
># Model 2: deprese ~ pohlavi + vekr + gpa + duv_v + duv_r + selfe + optim
># Res.Df RSS Df Sum of Sq F Pr(>F)
># 1 699 148.722
># 2 694 99.875 5 48.847 67.884 < 2.2e-16 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Balíček broom pak poskytuje několik funkcí
k práci se statistickými modely včetně modelů vytvořenýn pomocí
lm():
- funkce
tidy()slouží k extrakci koeficientů modelu do přehlednější tabulky - funkce
glance()slouži k extrakci informaci o fitu modelu - funkce
augment()slouží k uložení informací o jednotlivých případech včetně:- pridikovaných hodnot (
.fitted), - standardní chyby predikce (predikovaných hodnot)
(
.se.fit), - hodnot reziduí (
.resid), - leverage čili hat values čili páka (
.hat), - standardní chyby reziduí po vyřazení daného případu
(
.sigma), - Cookových vzdáleností (
.cooksd), - standardizovaných reziduí (
.std.resid), - intervalu spolehlivosti anebo predikčního intervalu
(
.lowera.upper), podle toho, jestliinterval = "confidence"nebointerval = "prediction"
- pridikovaných hodnot (
# Koeficienty modelu v přehledné tabulce (tibble)
broom::tidy(lmfit2, conf.int = TRUE, conf.leve = 0.99)># # A tibble: 8 × 7
># term estimate std.error statistic p.value conf.low conf.high
># <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
># 1 (Intercept) 3.30 0.172 19.2 5.51e-66 2.85 3.74
># 2 pohlavizenske 0.140 0.0315 4.45 9.96e- 6 0.0588 0.221
># 3 vekr 0.00358 0.00813 0.440 6.60e- 1 -0.0174 0.0246
># 4 gpa 0.0991 0.0209 4.73 2.70e- 6 0.0450 0.153
># 5 duv_v 0.0373 0.0303 1.23 2.18e- 1 -0.0409 0.116
># 6 duv_r -0.0495 0.0283 -1.75 8.10e- 2 -0.123 0.0236
># 7 selfe -0.442 0.0348 -12.7 2.31e-33 -0.531 -0.352
># 8 optim -0.0150 0.00525 -2.86 4.35e- 3 -0.0286 -0.00146# Informace o fitu modelu
broom::glance(lmfit2)># # A tibble: 1 × 12
># r.squared adj.r.squa…¹ sigma stati…² p.value df logLik AIC BIC devia…³
># <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
># 1 0.345 0.338 0.379 52.2 8.58e-60 7 -312. 641. 682. 99.9
># # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
># # variable names ¹adj.r.squared, ²statistic, ³deviance# Informací o jednotlivých případech
df_aug <- broom::augment(lmfit2, se_fit = TRUE,
interval = "confidence")
df_aug <- df_aug %>%
mutate(id = row_number()) %>%
select(id, everything())
df_aug># # A tibble: 702 × 18
># id deprese pohlavi vekr gpa duv_v duv_r selfe optim .fitted .lower
># <int> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
># 1 1 2.1 zenske 12.1 2 2.33 3.5 2.75 18 2.11 2.04
># 2 2 1.65 zenske 12.6 1 3.5 3.67 3.12 12 1.97 1.87
># 3 3 2.25 muzske 11.9 2 2.83 2.33 3 18 1.93 1.86
># 4 4 1.25 zenske 12.1 1 3 3.83 4 20 1.43 1.35
># 5 5 2.21 zenske 11.9 2 2.6 2.83 3 17 2.05 2.00
># 6 6 2.05 zenske 12.7 1.67 2.33 3.33 2.75 17 2.10 2.04
># 7 7 1.65 zenske 12.1 1.67 2.83 2.33 3.12 20 1.95 1.88
># 8 8 1.16 muzske 12.7 1.67 1.8 3.8 3.12 13 1.81 1.72
># 9 9 2.45 zenske 12.1 2 2.83 2.83 3.12 16 2.02 1.97
># 10 10 1.7 muzske 12.7 1.33 2.33 3.33 3.38 18 1.63 1.57
># # … with 692 more rows, and 7 more variables: .upper <dbl>, .se.fit <dbl>,
># # .resid <dbl>, .hat <dbl>, .sigma <dbl>, .cooksd <dbl>, .std.resid <dbl>Pro základní grafickou diagnostiku modelu stačí na fit-objekt použít
funkce plot(). Ukážou se nám celkem čtyři grafy. V konzoli
vždy musíme stisktnou ENTER, aby se nám ukázal další graf.
- Residuals vs Fitted. Tento graf ukazuje predikované hodnoty a rezidua. Neměl by mezi nimi být žádný vztah, jinak to může svědčit o nonlinearitě. Rozptyl reziduí by také měl být nezávislný na predikovaných hodnotách (neměli bychom pozorovat žádný trychtýř).
- Normal Q-Q graf. Pomocí toho můžeme ověřit normalitu reziduí. Pokud je dodržena, všechny boty by měly ležet blízko vyznačené přímky.
- Scale-Location (Spread-Location). Jedná se o graf predikovaných hodnot a odmocniny ze čtverce standardizovaných reziduí. Pomocí něj můžeme ověřit homoskedasticitu (homogenitu rozptylu). Velikost reziduí by neměla záviset na predikovaných hodnotách.
- Residuals vs Leverage. Tento graf zobrazuje páku (leverage) standardizované reziduum. Slouží k identifikaci extrémních a vlívných případů. Hlavně by nás měly znepokojovat případy s vysokou pákou a vysokým reziduem, protože ty mohou ovlivnit regresní koeficienty nejvíce. Vysokou páku mají ty případy, které v jednom nebo více prediktorech mají extrémní hodnoty, ale nemusejí být problematické, pokud je model dobře predikuje (nemají vysoké reziduum).
plot(lmfit2)
Graf Cookových vzdáleností je z nějakého důvodu samostatně. Cookovy vzdálenosti jsou vlastně měřítkem toho, nakolik by se změnily predikované hodnoty všech ostatních případů, kdybychom daný případ vyřadili z modelu.
plot(lmfit2, 4,
id.n = 10) # vyznačit 10 případů s největší Cookovou vzdáleností
Dále balíček car nabízí množství funkcí pro
diagnostiku modelu. Například funkci residualPlots()pro
diagnostiku linearity. Na ose X jsou vždy jednotlivé prediktory a na ose
Y rezidua z modelu se všemi prediktory (jen poslední graf ukazuje opět
predikované hodnoty a rezidua). Pokud je model dobře specifikovaný,
žádný vztah mezi prediktory a rezidui by neměl existovat (jinak to
znamená, že mezi prediktory a závislou proměnnou ještě existuje nějaký
nemodelovaný vztah, např. kvůli nezohledněné nonlinearitě)
car::residualPlots(lmfit2)
># Test stat Pr(>|Test stat|)
># pohlavi
># vekr 0.7957 0.4265
># gpa -1.4436 0.1493
># duv_v -1.3424 0.1799
># duv_r 1.5464 0.1225
># selfe 0.8799 0.3792
># optim 1.0938 0.2744
># Tukey test 1.4589 0.1446Dále balíček car nabízí funkci
vif() (variance inflation factor) pro kontrolu
multikolinearity. Je to vlastně odhad toho, kolikanásobně jsou větší
čtverce standardních chyb regresních koeficientů ve srovnání se stavem,
kdy by prediktory vzájemně vůbec nekorelovaly.
car::vif(lmfit2) ># pohlavi vekr gpa duv_v duv_r selfe optim
># 1.164743 1.285046 1.217895 1.146962 1.215053 1.362514 1.304913Jedním z testů závislosti reziduí je Durbin-Watsonův test, který vlastně analyzuje autorelaci mezi soudedními rezidui (s lagem = 1). Přestavte si dva totožné sloupce s rezidui, přičemž jeden z nich je pouze o jeden řádek posunutý.
car::durbinWatsonTest(lmfit2)># lag Autocorrelation D-W Statistic p-value
># 1 -0.002274117 2.002334 0.99
># Alternative hypothesis: rho != 0Protože jsme si do df_aug uložili statistiky pro
jednotlivé případy. můžeme s nimi dále pracovat jakkoli uznáme za
vhodné. Například se můžeme podívat na 10 případů s nejvyšší Cookovou
vzdáleností.
df_aug %>%
slice_max(order_by = .cooksd, n = 10)># # A tibble: 10 × 18
># id deprese pohlavi vekr gpa duv_v duv_r selfe optim .fitted .lower
># <int> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
># 1 464 3.35 zenske 16.8 2.67 3 2.17 3.38 10 2.12 2.03
># 2 204 2.85 zenske 12.2 2.67 3 2 3.62 12 1.98 1.87
># 3 599 3.2 zenske 16.6 2 2.5 1.67 3.12 16 2.08 2.00
># 4 68 1.25 zenske 12.1 2.67 2 2 3 21 2.08 1.97
># 5 51 3.05 muzske 12.8 2.5 1.5 2 2.88 14 2.07 1.97
># 6 415 1.7 muzske 13 2 2.67 1.5 1.88 12 2.56 2.46
># 7 91 2.95 muzske 12.8 3 2.8 4 2.88 18 2.01 1.91
># 8 250 3 muzske 13.1 5 3 2.83 2.88 13 2.34 2.21
># 9 701 1.35 muzske 11.9 1 1.67 2.67 2.25 11 2.21 2.11
># 10 576 1.7 muzske 16 2.67 3.67 3.33 2.6 8 2.32 2.19
># # … with 7 more variables: .upper <dbl>, .se.fit <dbl>, .resid <dbl>,
># # .hat <dbl>, .sigma <dbl>, .cooksd <dbl>, .std.resid <dbl>Nebo s nejvyšším reziduem (po převodu na absolutní hodnoty), tedy případy, které model nejhůře predikuje.
df_aug %>%
slice_max(order_by = abs(.resid), n = 10)># # A tibble: 10 × 18
># id deprese pohlavi vekr gpa duv_v duv_r selfe optim .fitted .lower
># <int> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
># 1 464 3.35 zenske 16.8 2.67 3 2.17 3.38 10 2.12 2.03
># 2 599 3.2 zenske 16.6 2 2.5 1.67 3.12 16 2.08 2.00
># 3 95 1.25 zenske 12.1 2.67 2 2.67 2.57 16 2.31 2.23
># 4 85 2.7 muzske 11.8 2 2 3.33 3.25 20 1.71 1.64
># 5 362 3.05 zenske 12.6 1.33 2.5 3.17 2.75 18 2.06 2.00
># 6 51 3.05 muzske 12.8 2.5 1.5 2 2.88 14 2.07 1.97
># 7 91 2.95 muzske 12.8 3 2.8 4 2.88 18 2.01 1.91
># 8 196 2.55 zenske 12.7 1.67 2.83 4 3.75 17 1.64 1.57
># 9 292 3.2 zenske 12.5 2 2.33 2.33 2.62 13 2.30 2.23
># 10 143 1.55 zenske 12.2 1.67 3 2 2.25 14 2.45 2.37
># # … with 7 more variables: .upper <dbl>, .se.fit <dbl>, .resid <dbl>,
># # .hat <dbl>, .sigma <dbl>, .cooksd <dbl>, .std.resid <dbl>Také můžeme vytvářet jakékoli diagnostické či jiné grafy.
df_aug %>%
ggplot(aes(.hat, .std.resid, size = .cooksd)) +
geom_text(aes(label = id))
Také bychom mohli klidně provést “analýzu senzitivity” a podívat se,
jak by výsledné odhady modelu ovlivnilo vyřazení extrémních případů
(pomocí funkce filter() bychom si stanovili nějaký cut-off,
např. na základě Cookových vzdáleností).
5 Srovnání dvou nezávislých skupin
5.1 Nezávislý t-test
Zkusíme pomocí nezávislého t-testu ověřit rozdíl v průměru známek mezi chlapci a dívkami (technicky přesněji to je vlastně rozdíl mezi průměry průměrů známek chlapců a dívek)
Nejdříve se podíváme na rozdělení závislé proměnné u obou skupin graficky. Můžeme si vyzkoušet různé typy grafů.
# Histogram s relativní četností
# Pozn. hustato pravděpodobnosti * binwidth = relativní četnost
df_clean %>%
ggplot(aes(gpa, y = stat(density)*0.25)) + # Hustota pravděpodobnosti
geom_histogram(binwidth = 0.25) +
facet_wrap(~pohlavi)
# Graf hustoty pravděpodobnosti
df_clean %>%
ggplot(aes(gpa, color = pohlavi)) +
geom_density(size = 2)
# Boxplot
df_clean %>%
ggplot(aes(pohlavi, gpa)) +
geom_boxplot()
# Q-Q grafy
df_clean %>%
ggplot(aes(sample = gpa)) +
geom_qq() +
geom_qq_line() +
facet_wrap(~pohlavi)
Vypočteme si relevantní deskriptivní statistiky.
df_clean %>%
group_by(pohlavi) %>%
summarise(
n = n(), # počet případů
ci_mean(gpa, conf.level = .95), # průměr s intervalem spolehlivosti
SD = sd(gpa), # směrodatná odchylka
skew = psych::skew(gpa) # šikmost
)># # A tibble: 2 × 7
># pohlavi n M Mci_low Mci_upp SD skew
># <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
># 1 muzske 284 2.36 2.27 2.45 0.782 0.229
># 2 zenske 418 2.12 2.05 2.19 0.721 0.338Zkusit můžeme také test normality (Shapiro-Wilk).
df_clean %>%
group_by(pohlavi) %>%
dlookr::normality(gpa)># # A tibble: 2 × 5
># variable pohlavi statistic p_value sample
># <chr> <fct> <dbl> <dbl> <dbl>
># 1 gpa muzske 0.973 0.0000310 284
># 2 gpa zenske 0.958 0.00000000163 418Samotný nezávislý t-test provádíme pomocí funkce t.test
z base R.
t_res <- t.test(gpa ~ pohlavi, # Závislá proměnná ~ nezávislá proměnná
data = df_clean, # Použitý dataset
conf.level = .995, # úroveň spolehlivosti
var.equal = TRUE) # předpokládat shodu rozptylů?
t_res>#
># Two Sample t-test
>#
># data: gpa by pohlavi
># t = 4.231, df = 700, p-value = 2.636e-05
># alternative hypothesis: true difference in means between group muzske and group zenske is not equal to 0
># 99.5 percent confidence interval:
># 0.08117923 0.40427576
># sample estimates:
># mean in group muzske mean in group zenske
># 2.359155 2.116427Output této funkce obsahuje údaje o testové statistice, stupních
volnosti, p-hodnotě, intervalu spolehlivosti pro rozdíl mezi průměry a
nakonec o průměrech porovnávaných skupin. Defaultně funkce nepředpokládá
shodu rozptylů (argument var.equal je defaultně
FALSE), což znamená použití Welchova t-testu. Obvykle je to
robustnější alternativa, ale může mít menší statistickou sílu, pokud je
předpoklad homogenity rozptylu naplněn.
Balíček effectsize pak nabízí několik funkcí k výpočtu standardizovaných velikostí účinku včetně:
- cohens_d(). Rozdíl mezi průměry standardizovaný (dělený) směrodatnou odchylkou.
- hedges_g(). Podobné jako Cohenovo d, ale zahrnuje korekci pro nadhodnocení (bias) velikosti účinku, které hrozí hlavně u malých vzorků.
- glass_delta(). Používá ke standardizaci pouze směrodatnou odchylku druhé skupiny.
Argumentem pooled_sd můžeme specifikovat, zda se má ke
standardizaci použít společná směrodatná ochylka (pooled standard
deviation), což má smysl tehdy, je-li dodržena homogenita rozptylů.
Defaultně je pooled_sd = FALSE a ke standardizaci se
použije průměr směrodatných odchylek obou skupin.
effectsize::cohens_d(gpa ~ pohlavi, data = df_clean)># Cohen's d | 95% CI
># ------------------------
># 0.33 | [0.17, 0.48]
>#
># - Estimated using pooled SD.effectsize::cohens_d(gpa ~ pohlavi,
data = df_clean,
pooled_sd = FALSE, # Použít společnou směrodatnou odchylku?
ci = .995) # Úroveň spolehlivosti># Cohen's d | 99.5% CI
># ------------------------
># 0.32 | [0.10, 0.54]
>#
># - Estimated using un-pooled SD.5.2 Wilcoxon-Mann–Whitney (WMW) U-test
Pomocí WMW test testu můžeme ověřit rozdíl mezi chlapci a dívkami ve známce z češtiny, což je spíš ordinální proměnná. Tento test je vlastně neparametrická obdoba nezávislého t-testu, která ověřuje to, zda se skupiny liší v pravděpodobnosti superiority (nulová hypotéza je ta, že pravděpodobnost, že náhodně vybraný jedinec z populace 1 bude mít v dané proměnné vyšší skór, je totožná s pravděpodobností, že náhodně vybraný jedinec z populace 2 bude mít v dané proměnné vyšší skór, tj. p1sup = p2sup = 0,5).
df_clean <- df %>%
drop_na(pohlavi, zn_cj)Hlavní kvantily se příliš neliší, a to právě proto, že nemáme moc jedinečných hodnot.
df_clean %>%
group_by(pohlavi) %>%
summarise(
quantiles(zn_cj)
)># # A tibble: 2 × 6
># pohlavi Min Q2 Mdn Q3 Max
># <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
># 1 muzske 1 2 2 3 5
># 2 zenske 1 2 2 3 5Ale na základě četností to vypadá, že dívky dostávají častěji jedničky a dvojky než chlapci, zatímco chlapci častěji trojky a čtyřky než dívky
freq <- df_clean %>%
count(pohlavi, zn_cj) %>%
group_by(pohlavi) %>%
mutate(p = n/sum(n))
freq># # A tibble: 10 × 4
># # Groups: pohlavi [2]
># pohlavi zn_cj n p
># <fct> <dbl> <int> <dbl>
># 1 muzske 1 35 0.111
># 2 muzske 2 140 0.444
># 3 muzske 3 103 0.327
># 4 muzske 4 35 0.111
># 5 muzske 5 2 0.00635
># 6 zenske 1 95 0.212
># 7 zenske 2 219 0.489
># 8 zenske 3 110 0.246
># 9 zenske 4 21 0.0469
># 10 zenske 5 3 0.00670freq %>%
ggplot(aes(zn_cj, p, fill = pohlavi)) +
geom_col(position = "dodge")
Samotný WMW test provedeme takto.
wilcox.test(zn_cj ~ pohlavi, data = df_clean)>#
># Wilcoxon rank sum test with continuity correction
>#
># data: zn_cj by pohlavi
># W = 84400, p-value = 7.284e-07
># alternative hypothesis: true location shift is not equal to 0Kdychom chtěli znát velikost účinku s intervalem spolehlivosti, což
je mnohem informativnější než jen testová statistika s p-hodnotou,
můžeme použít funkci p_superiority() z balíčku
effectsize. Ten nabízí mnoho funkcí pro
výpočet různých jiných effect sizes, jak se můžete podívat zde.
effectsize::p_superiority(zn_cj ~ pohlavi,
data = df_clean,
iterations = 2000, # Bootstrap s 2000 iteracemi
parametric = FALSE) # Nepředpokládat normální rozdělení># Pr(superiority) | 95% CI
># ------------------------------
># 0.60 | [0.56, 0.64]
>#
># - Non-parametric CLES5.3 Párový t-test
Pomocí párového t-testu můžeme ověřit rozdíl mezi vnímanou vřelostí matky a otce. Párový t-test tedy ověřuje rozdíl mezi průměry dvou závislých měření. Také jej můžeme chápat jako test toho, že průměr rozdílových skórů je roven nule (alespoň taková je nejčastější nulová hypotéza). Nejdříve si vyčistíme data.
# Příprava dat
df_clean <- df %>%
drop_na(warm_m, warm_o)Poté si vypočteme i rozdílové skóry (warm_diff),
průměrnou vřelost pro oba rodiče (warm_both), celkový
průměr (grand_mean), mezisubjektový efekt
btw_eff (odchylku warm_both od
grand_mean) a adjustovaný skóry pro vřelost matky a otce (s
odečtením mezisubjektového efektu). Pro zopakování
warm_diff: rozdíl ve vřelosti matky – otce.warm_both: průměrná vřelost žáka pro oba rodiče.grand_mean: celkový průměr z hodnot warm_m i warm_obtw_eff: mezisubjektový efekt, o kolik se průměrná vřelost vůči oběma rodičůmwarm_bothodchyluje od celkového průměru grand_meanwarm_m_adjawarm_o_adj: adjustované hodnoty vřelosti pro odečtení mezisubjektového efektubtw_effodwarm_mawarm_o.
# Převod na delší formát
df_clean <- df_clean %>%
select(id, warm_m, warm_o) %>%
mutate(warm_diff = warm_m - warm_o, # rozdíl mezi vřelostí matky a otce
warm_both = (warm_m + warm_o)/2, # průměrná vřelost žáka pro oba rodiče
grand_mean = mean(warm_both), # celkový průměr obou typů vřelosti
btw_eff = warm_both - grand_mean, # odhad efektu žáka
warm_m_adj = warm_m - btw_eff,
warm_o_adj = warm_o - btw_eff
)
head(df_clean)># # A tibble: 6 × 9
># id warm_m warm_o warm_diff warm_both grand_mean btw_eff warm_m_adj warm_…¹
># <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
># 1 01-001 3.55 3.55 0 3.55 3.16 0.386 3.16 3.16
># 2 01-002 4 3.27 0.727 3.64 3.16 0.477 3.52 2.80
># 3 01-003 3.45 2.73 0.727 3.09 3.16 -0.0682 3.52 2.80
># 4 01-004 3.2 3.7 -0.5 3.45 3.16 0.291 2.91 3.41
># 5 01-005 4 3.73 0.273 3.86 3.16 0.705 3.30 3.02
># 6 01-006 3.6 3.7 -0.100 3.65 3.16 0.491 3.11 3.21
># # … with abbreviated variable name ¹warm_o_adjDále si můžeme vypočíst základní desriptivní statistiky.
df_clean %>%
skimr::skim(warm_m, warm_o, warm_diff)># ── Data Summary ────────────────────────
># Values
># Name Piped data
># Number of rows 743
># Number of columns 9
># _______________________
># Column type frequency:
># numeric 3
># ________________________
># Group variables None
>#
># ── Variable type: numeric ──────────────────────────────────────────────────────
># skim_variable n_missing complete_rate mean sd p0 p25 p50 p75
># 1 warm_m 0 1 3.23 0.470 1.18 3 3.27 3.55
># 2 warm_o 0 1 3.08 0.584 1 2.82 3.18 3.55
># 3 warm_diff 0 1 0.151 0.540 -2.03 -0.0909 0.0909 0.273
># p100 hist
># 1 4 ▁▁▃▇▇
># 2 4 ▁▁▃▇▆
># 3 2.64 ▁▂▇▁▁Také nás obvykle zajímá i korelace mezi párem proměnných.
cor.test(df_clean$warm_m, df_clean$warm_o, conf.level = .995)>#
># Pearson's product-moment correlation
>#
># data: df_clean$warm_m and df_clean$warm_o
># t = 15.382, df = 741, p-value < 2.2e-16
># alternative hypothesis: true correlation is not equal to 0
># 99.5 percent confidence interval:
># 0.4098709 0.5661465
># sample estimates:
># cor
># 0.4919613Můžeme si také vytvořit pár grafů rozdílových skórů.
# Histogram s křivkou normálního rozdělení
df_clean %>%
ggplot(aes(warm_diff)) +
geom_histogram(aes(y = stat(density))) +
stat_function(fun = dnorm, args = list(mean = mean(df_clean$warm_diff),
sd = sd(df_clean$warm_diff)),
color = "red", size = 2)
# Q-Q graf
df_clean %>%
ggplot(aes(sample = warm_diff)) +
stat_qq() +
stat_qq_line()
Kdybychom chtěli chybové úsečky s korektními intervaly spolehlivosti, použili bychom adjustované skóry a museli bychom si data převést do longer formátu.
df_long <- df_clean %>%
select(id, warm_m_adj, warm_o_adj) %>%
pivot_longer(cols = c(warm_m_adj, warm_o_adj), # Tyto sloupce sločit do jednoho
names_to = "warm", # Nový sloupec s identifikátorem původních sloupců
values_to = "value") # Nový sloupec z hodnotami
head(df_long)># # A tibble: 6 × 3
># id warm value
># <chr> <chr> <dbl>
># 1 01-001 warm_m_adj 3.16
># 2 01-001 warm_o_adj 3.16
># 3 01-002 warm_m_adj 3.52
># 4 01-002 warm_o_adj 2.80
># 5 01-003 warm_m_adj 3.52
># 6 01-003 warm_o_adj 2.80df_long %>%
ggplot(aes(warm, value)) +
geom_pointrange(stat = "summary",
fun.data = mean_cl_normal,
fun.args = list(conf.int= .995))
Samotná párový t-test můžeme provést takto.
t.test(df_clean$warm_m, df_clean$warm_o, # Porovnávaný pár proměnných
paired = TRUE, conf.level = .995)>#
># Paired t-test
>#
># data: df_clean$warm_m and df_clean$warm_o
># t = 7.6086, df = 742, p-value = 8.417e-14
># alternative hypothesis: true mean difference is not equal to 0
># 99.5 percent confidence interval:
># 0.09494551 0.20648662
># sample estimates:
># mean difference
># 0.1507161Jako standardizovanou velikost účinku můžeme vypočíst Cohenovo
d. Stačí opět použít funkci cohens_d() z balíčku
effectsize. Pro standardizaci se při
paired = TRUE použije směrodatná odchylka rozdílových
skórů.
effectsize::cohens_d(df_clean$warm_m, df_clean$warm_o,
paired = TRUE,
ci = .995)># Cohen's d | 99.5% CI
># ------------------------
># 0.28 | [0.17, 0.45]5.4 Wilcoxonův párový (sign-rank) test
Neparametrickou obdobou párového t-testu je Wilcoxonův párový test.
Podobně jako WMW test jej provedeme pomocí funkce
wilcox.test(), ale musíme nastavit
paired = TRUE.
wilcox.test(df_clean$warm_m, df_clean$warm_o,
paired = TRUE)>#
># Wilcoxon signed rank test with continuity correction
>#
># data: df_clean$warm_m and df_clean$warm_o
># V = 114689, p-value = 4.428e-14
># alternative hypothesis: true location shift is not equal to 0Jako velikost účinku opět můžeme vypočíst pravděpodobnost
superiority, ale opět nastavíme paired = TRUE.
effectsize::p_superiority(df_clean$warm_m, df_clean$warm_o,
data = df_clean,
paired = TRUE,
iterations = 2000, # Bootstrap s 2000 iteracemi
parametric = FALSE) # Nepředpokládat normální rozdělení># Pr(superiority) | 95% CI
># ------------------------------
># 0.68 | [0.64, 0.72]
>#
># - Non-parametric CLES6 Srovnání více nezávislých skupin
6.1 One-way ANOVA
V případě, že je závislá proměnná spojitá a má smysl srovnávat průměry, můžeme použít one-way ANOVA. Zkusíme si ověřit rozdíly v průměru známek mezi skupinami rozdělenými podle očekávaného vzdělání.
df_clean <- df %>%
select(id, skola, trida, ocek_vzd, gpa) %>%
drop_na()Podíváme se na četnosti skupina a základní deskriptivní statistiky.
stats <- df_clean %>%
group_by(ocek_vzd) %>%
summarise(
n = n(),
ci_mean(gpa),
SD = sd(gpa),
var = var(gpa),
skew = psych::skew(gpa),
kurt = psych::kurtosi(gpa)
)
stats># # A tibble: 6 × 9
># ocek_vzd n M Mci_low Mci_upp SD var skew kurt
># <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
># 1 základní 8 3.38 2.80 3.95 0.683 0.466 0.618 -0.941
># 2 vyučení 33 2.96 2.66 3.26 0.847 0.717 -0.659 -0.193
># 3 vyučení s maturitou 86 2.45 2.30 2.60 0.711 0.505 0.523 0.393
># 4 odborná střední škola 142 2.53 2.40 2.65 0.762 0.580 0.170 -0.360
># 5 gymnázium 42 2.21 2.01 2.41 0.647 0.419 0.473 0.432
># 6 vysoká škola 432 2.01 1.95 2.08 0.682 0.465 0.267 -0.663Pro představu o rozdělení v rámci skupin se můžeme podívat na boxploty.
df_clean %>%
ggplot(aes(ocek_vzd, gpa)) +
geom_boxplot()
Protože máme deskriptivní statistiky již vypočteny, můžeme na jejich základě vytvořit graf průměrů s intervaly spolehlivosti.
stats %>%
ggplot() +
geom_pointrange(aes(x = ocek_vzd,
y = M,
ymin = Mci_low,
ymax = Mci_upp))
Histogram a graf hustoty pravděpodobnosti by vypadaly takto.
df_clean %>%
ggplot(aes(gpa)) +
geom_histogram(binwidth = .25) +
facet_wrap(~ocek_vzd, scales = "free_y")
df_clean %>%
ggplot(aes(gpa)) +
geom_density() +
facet_wrap(~ocek_vzd)
Také se můžeme podívat na Q-Q grafy.
df_clean %>%
ggplot(aes(sample = gpa)) +
geom_qq() +
geom_qq_line() +
facet_wrap(~ocek_vzd)
Nebo provést Levenův test shody rozptylů.
car::leveneTest(gpa ~ ocek_vzd, data = df_clean)># Levene's Test for Homogeneity of Variance (center = median)
># Df F value Pr(>F)
># group 5 1.4868 0.1918
># 737A nakonec můžeme provést samotnou one-way ANOVA.
av1 <- aov(gpa ~ ocek_vzd, data = df_clean)
summary(av1)># Df Sum Sq Mean Sq F value Pr(>F)
># ocek_vzd 5 65.7 13.13 26.28 <2e-16 ***
># Residuals 737 368.4 0.50
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1broom::tidy(av1)># # A tibble: 2 × 6
># term df sumsq meansq statistic p.value
># <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
># 1 ocek_vzd 5 65.7 13.1 26.3 1.80e-24
># 2 Residuals 737 368. 0.500 NA NAbroom::glance(av1)># # A tibble: 1 × 6
># logLik AIC BIC deviance nobs r.squared
># <dbl> <dbl> <dbl> <dbl> <int> <dbl>
># 1 -794. 1601. 1634. 368. 743 0.151R nabízí i Welchovu ANOVA, kerá nepředpokládá shodu rozptylů.
av2 <- oneway.test(gpa ~ ocek_vzd, data = df_clean)
av2>#
># One-way analysis of means (not assuming equal variances)
>#
># data: gpa and ocek_vzd
># F = 23.029, num df = 5.000, denom df = 56.958, p-value = 1.441e-12broom::tidy(av2)># # A tibble: 1 × 5
># num.df den.df statistic p.value method
># <dbl> <dbl> <dbl> <dbl> <chr>
># 1 5 57.0 23.0 1.44e-12 One-way analysis of means (not assuming equa…Diagnostické grafy získáme stejně jako v případě lineární regrese
použitím funkce plot() na fit-objekt.
plot(av1)
Balíček effectsize() obsahuje i funkce pro výpočet
celkového efektu pro ANOVA, např. omega-squared.
effectsize::omega_squared(av1,
alternative = "two.sided",
ci = .995)># # Effect Size for ANOVA
>#
># Parameter | Omega2 | 99.5% CI
># ---------------------------------
># ocek_vzd | 0.15 | [0.08, 0.21]Dejte tomu, že ANOVA chceme doplnit o dvě plánovaná srovnání (ortogonální kontrasty), a to žáků, kteří
- očekávají vs. neočekávají, že budou mít maturitu,
- očekávají maturitu, ale očekávají vs. neočekávají, že budou mít vysokou školu.
Nejdříve si musíme stanovit váhy pro jednotlivé kontrasty. Pro jistotu si zobrazíme pořadí jednotlivých úrovní.
levels(df_clean$ocek_vzd)># [1] "základní" "vyučení" "vyučení s maturitou"
># [4] "odborná střední škola" "gymnázium" "vysoká škola"Pak si definujeme matici (matrix) s konstrastky takto
m <- matrix(
c(-1/2,-1/2, 1/4, 1/4, 1/4, 1/4, # úrovně bez maturity vs s maturitou
0, 0,-1/3,-1/3,-1/3, 1), # úrovně s maturitou, ale bez VŠ, vs. VŠ
nrow = 2, byrow = TRUE,
dimnames = list(
c("Maturita", "VŠ"), # rpwnames
levels(df_clean$ocek_vzd) # colnames
))
m># základní vyučení vyučení s maturitou odborná střední škola gymnázium
># Maturita -0.5 -0.5 0.2500000 0.2500000 0.2500000
># VŠ 0.0 0.0 -0.3333333 -0.3333333 -0.3333333
># vysoká škola
># Maturita 0.25
># VŠ 1.00K odhadu samotných konstrastů pak použijeme funkci
glht() z balíčku multcomp.
planned <- multcomp::glht(av1, linfct = m)
summary(planned)>#
># Simultaneous Tests for General Linear Hypotheses
>#
># Fit: aov(formula = gpa ~ ocek_vzd, data = df_clean)
>#
># Linear Hypotheses:
># Estimate Std. Error t value Pr(>|t|)
># Maturita == 0 -2.5578 0.2601 -9.833 < 1e-10 ***
># VŠ == 0 -0.3845 0.0593 -6.483 3.28e-10 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
># (Adjusted p values reported -- single-step method)confint(planned) # pro intervaly spolehlivosti>#
># Simultaneous Confidence Intervals
>#
># Fit: aov(formula = gpa ~ ocek_vzd, data = df_clean)
>#
># Quantile = 2.2401
># 95% family-wise confidence level
>#
>#
># Linear Hypotheses:
># Estimate lwr upr
># Maturita == 0 -2.5578 -3.1405 -1.9751
># VŠ == 0 -0.3845 -0.5173 -0.2516Kdybychom neměli apriorní hypotézy, ale chtěli provést post-hoc testy, balíček rstatix nabízí ty nejznámější. Např. Tukeyho HSD post-hoc test má dobrou statistickou sílu ale asi nejvíce předpokladů včetně normálního rozdělení v rámci skupin a shody rozptylů
df_clean %>%
rstatix::tukey_hsd(gpa ~ ocek_vzd) %>%
select(group1,group2, estimate, conf.low, conf.high, p.adj)># # A tibble: 15 × 6
># group1 group2 estimate conf.low conf.…¹ p.adj
># <chr> <chr> <dbl> <dbl> <dbl> <dbl>
># 1 základní vyučení -0.410 -1.21 0.386 6.82e-1
># 2 základní vyučení s maturitou -0.925 -1.67 -0.179 5.64e-3
># 3 základní odborná střední škola -0.848 -1.58 -0.114 1.29e-2
># 4 základní gymnázium -1.16 -1.94 -0.385 3.17e-4
># 5 základní vysoká škola -1.36 -2.08 -0.643 1.29e-6
># 6 vyučení vyučení s maturitou -0.515 -0.929 -0.101 5.32e-3
># 7 vyučení odborná střední škola -0.438 -0.828 -0.0473 1.77e-2
># 8 vyučení gymnázium -0.754 -1.22 -0.284 7.73e-5
># 9 vyučení vysoká škola -0.953 -1.32 -0.589 0
># 10 vyučení s maturitou odborná střední škola 0.0774 -0.199 0.353 9.67e-1
># 11 vyučení s maturitou gymnázium -0.239 -0.620 0.141 4.68e-1
># 12 vyučení s maturitou vysoká škola -0.438 -0.677 -0.200 2.93e-6
># 13 odborná střední škola gymnázium -0.317 -0.671 0.0381 1.11e-1
># 14 odborná střední škola vysoká škola -0.516 -0.711 -0.320 0
># 15 gymnázium vysoká škola -0.199 -0.526 0.127 5.04e-1
># # … with abbreviated variable name ¹conf.highGames-Howellův post-hoc test shodu rozptylů nepřepokládá, ale má menší stastickou sílu.
df_clean %>%
rstatix::games_howell_test(gpa ~ ocek_vzd) %>%
select(group1,group2, estimate, conf.low, conf.high, p.adj)># # A tibble: 15 × 6
># group1 group2 estimate conf.…¹ conf.…² p.adj
># <chr> <chr> <dbl> <dbl> <dbl> <dbl>
># 1 základní vyučení -0.410 -1.35 0.530 6.98e- 1
># 2 základní vyučení s maturitou -0.925 -1.84 -0.0130 4.7 e- 2
># 3 základní odborná střední škola -0.848 -1.76 0.0641 7 e- 2
># 4 základní gymnázium -1.16 -2.08 -0.249 1.2 e- 2
># 5 základní vysoká škola -1.36 -2.28 -0.450 6 e- 3
># 6 vyučení vyučení s maturitou -0.515 -1.01 -0.0230 3.5 e- 2
># 7 vyučení odborná střední škola -0.438 -0.916 0.0405 9 e- 2
># 8 vyučení gymnázium -0.754 -1.28 -0.230 1 e- 3
># 9 vyučení vysoká škola -0.953 -1.41 -0.499 4.11e- 6
># 10 vyučení s maturitou odborná střední škola 0.0774 -0.210 0.365 9.71e- 1
># 11 vyučení s maturitou gymnázium -0.239 -0.606 0.127 4.08e- 1
># 12 vyučení s maturitou vysoká škola -0.438 -0.680 -0.197 9.61e- 6
># 13 odborná střední škola gymnázium -0.317 -0.663 0.0298 9.3 e- 2
># 14 odborná střední škola vysoká škola -0.516 -0.722 -0.309 1.6 e-10
># 15 gymnázium vysoká škola -0.199 -0.510 0.112 4.17e- 1
># # … with abbreviated variable names ¹conf.low, ²conf.highBalíček rstatix má i vlastní funkci pro výpočet Cohenova d, která
zvládá i více skupin (Defaultně se vypočte Cohenovo d pro všechna
srovnání). O interval spolehlivosti bychom si mohli požádat pomocí
argumentu ci = TRUE, ale výpočet probíhá pomocí bootstrapu,
což může trvat delší dobu.
df_clean %>%
rstatix::cohens_d(gpa ~ ocek_vzd,
var.equal = TRUE) # Předpkládat shodu rozptylů?># # A tibble: 15 × 7
># .y. group1 group2 effsize n1 n2 magni…¹
># * <chr> <chr> <chr> <dbl> <int> <int> <ord>
># 1 gpa základní vyučení 0.501 8 33 modera…
># 2 gpa základní vyučení s maturitou 1.31 8 86 large
># 3 gpa základní odborná střední škola 1.12 8 142 large
># 4 gpa základní gymnázium 1.79 8 42 large
># 5 gpa základní vysoká škola 2.00 8 432 large
># 6 gpa vyučení vyučení s maturitou 0.686 33 86 modera…
># 7 gpa vyučení odborná střední škola 0.562 33 142 modera…
># 8 gpa vyučení gymnázium 1.02 33 42 large
># 9 gpa vyučení vysoká škola 1.37 33 432 large
># 10 gpa vyučení s maturitou odborná střední škola -0.104 86 142 neglig…
># 11 gpa vyučení s maturitou gymnázium 0.346 86 42 small
># 12 gpa vyučení s maturitou vysoká škola 0.639 86 432 modera…
># 13 gpa odborná střední škola gymnázium 0.429 142 42 small
># 14 gpa odborná střední škola vysoká škola 0.735 142 432 modera…
># 15 gpa gymnázium vysoká škola 0.293 42 432 small
># # … with abbreviated variable name ¹magnitude6.2 Kruskal-Wallisův Test
Neparametrickou alternativou one-way ANOVA je Kruskal-Wallisův Test. Abychom si věci nekomplikovali, zkusíme opět testovat vztah mezi očekávaných vzděláním a průměrem známek.
kw1 <- kruskal.test(gpa ~ ocek_vzd, data = df_clean)
kw1>#
># Kruskal-Wallis rank sum test
>#
># data: gpa by ocek_vzd
># Kruskal-Wallis chi-squared = 99.367, df = 5, p-value < 2.2e-16Jako souhrnnou velikost účinku se doporučuje uvádět epsilon-squred nebo eta-squared.
df_clean %>%
rstatix::kruskal_effsize(gpa ~ ocek_vzd,
ci = FALSE, # vypnuto, ať se nezdržujeme
conf.level = .95,
nboot = 1000)># # A tibble: 1 × 5
># .y. n effsize method magnitude
># * <chr> <int> <dbl> <chr> <ord>
># 1 gpa 743 0.128 eta2[H] moderateJde vlastně o odhad vysvětleného rozptylu po pořadové transformaci závislé proměnné, jak se můžeme přesvědčit níže (adjusted R-squred)
kw2 <- lm(rank(gpa) ~ ocek_vzd, data = df_clean)
summary(kw2)$r.squared># [1] 0.1339182summary(kw2)$adj.r.squared># [1] 0.1280425A ještě si ukážeme párová post-hoc srovnání pomocí Dunnova testu, což je neparametrická obdoba Tukeyho HSD testu.
df_clean %>%
rstatix::dunn_test(gpa ~ ocek_vzd,
p.adjust.method = "holm") %>%
select(group1, group2, statistic, p, p.adj) %>%
mutate(across(where(is.numeric), round, digits = 4))># # A tibble: 15 × 5
># group1 group2 statistic p p.adj
># <chr> <chr> <dbl> <dbl> <dbl>
># 1 základní vyučení -1.08 0.280 0.561
># 2 základní vyučení s maturitou -2.73 0.0063 0.0507
># 3 základní odborná střední škola -2.50 0.0123 0.0862
># 4 základní gymnázium -3.38 0.0007 0.0072
># 5 základní vysoká škola -4.40 0 0.0001
># 6 vyučení vyučení s maturitou -2.85 0.0044 0.0394
># 7 vyučení odborná střední škola -2.50 0.0123 0.0862
># 8 vyučení gymnázium -3.78 0.0002 0.0018
># 9 vyučení vysoká škola -6.34 0 0
># 10 vyučení s maturitou odborná střední škola 0.731 0.465 0.561
># 11 vyučení s maturitou gymnázium -1.57 0.117 0.396
># 12 vyučení s maturitou vysoká škola -4.76 0 0
># 13 odborná střední škola gymnázium -2.25 0.0246 0.123
># 14 odborná střední škola vysoká škola -6.84 0 0
># 15 gymnázium vysoká škola -1.65 0.0989 0.3967 Vztah mezi dvěma kategorickými proměnnými
Nakonec se podíváme na vztah mezi dvěma kategorickými proměnnými, např. známkou z češtiny a matematiky
df_clean <- df %>%
drop_na(zn_cj, zn_mat)Nejprve se podíváme na kontingenční tabulku pomocí funkce crosstabs z
r101_funs.R.
crosstabs(df, row_var = "zn_cj", col_var = "zn_mat")>#
># Pearson's Chi-squared test
>#
># data: freq
># X-squared = 567.21, df = 16, p-value < 2.2e-16># # A tibble: 20 × 7
># stat zn_cj `1` `2` `3` `4` `5`
># <chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
># 1 n 1 82 40 8 1 0
># 2 n 2 71 200 73 17 0
># 3 n 3 4 75 109 25 0
># 4 n 4 1 7 23 24 1
># 5 n 5 0 0 2 1 2
># 6 row % 1 62.6 30.5 6.11 0.76 0
># 7 row % 2 19.7 55.4 20.2 4.71 0
># 8 row % 3 1.88 35.2 51.2 11.7 0
># 9 row % 4 1.79 12.5 41.1 42.9 1.79
># 10 row % 5 0 0 40 20 40
># 11 col % 1 51.9 12.4 3.72 1.47 0
># 12 col % 2 44.9 62.1 34.0 25 0
># 13 col % 3 2.53 23.3 50.7 36.8 0
># 14 col % 4 0.63 2.17 10.7 35.3 33.3
># 15 col % 5 0 0 0.93 1.47 66.7
># 16 tot % 1 10.7 5.22 1.04 0.13 0
># 17 tot % 2 9.27 26.1 9.53 2.22 0
># 18 tot % 3 0.52 9.79 14.2 3.26 0
># 19 tot % 4 0.13 0.91 3 3.13 0.13
># 20 tot % 5 0 0 0.26 0.13 0.26Graficky bychom vztah mohli znázornit pomocí mozaikového grafu nebo pomocí bodového grafu s různou velikostí bodů podle zastoupení dané kombinace.
df_clean %>%
plot_mosaic(xvar = "zn_mat", yvar = "zn_cj") +
scale_fill_manual(values = palette.colors(n = 5, palette = "R4"))
df_clean %>%
ggplot(aes(zn_cj, zn_mat)) +
geom_count() +
scale_size_continuous(
range = c(1, 10)
)
A statisticky bychom mohli vztah mezi oběma proměnnými ověřit pomocí
testu nezávislosti chí-kvadrát
chisq.test(df$zn_cj, df$zn_mat)>#
># Pearson's Chi-squared test
>#
># data: df$zn_cj and df$zn_mat
># X-squared = 567.21, df = 16, p-value < 2.2e-16chisq.test(df$zn_cj, df$zn_mat,
simulate.p.value = TRUE, # p-hodnoty na základě simulace
B = 1e4) # počet simulací>#
># Pearson's Chi-squared test with simulated p-value (based on 10000
># replicates)
>#
># data: df$zn_cj and df$zn_mat
># X-squared = 567.21, df = NA, p-value = 9.999e-05Jako velikost účinku bývá nejčastěji v kontextu reportování chí-kvadrát testu uváděno Cramerovo V.
effectsize::cramers_v(df$zn_cj, df$zn_mat,
ci = .995, # hladina spolehlivosti
alternative = "two.sided", # typ intervalu
adjust = TRUE) # korekce biasu># Cramer's V (adj.) | 99.5% CI
># --------------------------------
># 0.43 | [0.37, 0.47]