Smíšený model a interakce
library(readxl)
library(tidyverse)
library(gridBase)
library(grid)
library(lme4)
data <- read_excel("data1.xlsx")
# data$Age <- data$Age - mean(data$Age) # centrování - na p-hodnoty nemá vlivGraf pro zprůměrované hodnoty v rámci každého mluvčího
dataAV <- data %>% group_by(ID) %>% summarise_at(vars(Trvani:TempoSl), mean) %>% left_join(select(data, ID, Sex, Age) %>%
distinct, by = "ID")
ggplot(dataAV, aes(x = Age, y = F0, color = Sex)) + geom_point() + geom_smooth(formula = y ~ x, method = "lm")
F0 vs. Age*Sex
fit <- lmer(F0 ~ Age*Sex + (1|ID), data) # pro interpretace modelu se REML nevypíná
plot(fit) # graf odchylek -> o.k.
summary(fit)## Linear mixed model fit by REML ['lmerMod']
## Formula: F0 ~ Age * Sex + (1 | ID)
## Data: data
##
## REML criterion at convergence: 1049.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5632 -0.5457 0.0062 0.5464 4.6027
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 5.7441 2.3967
## Residual 0.7335 0.8564
## Number of obs: 360, groups: ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 16.54913 2.24539 7.370
## Age -0.09975 0.04377 -2.279
## SexM -11.61332 3.08625 -3.763
## Age:SexM 0.04738 0.06068 0.781
##
## Correlation of Fixed Effects:
## (Intr) Age SexM
## Age -0.946
## SexM -0.728 0.688
## Age:SexM 0.682 -0.721 -0.955“Correlation of Fixed Effects” je věc, do jejíhož vysvětlování se mnoho lidí raději vůbec nepouští (nepletu-li se, tak v knize Bodo Wintera jsem interpretaci nenašel). Řeknu-li to uživatelsky, hlavní je, že to není přesně rovno 1, to by pak bylo velmi podezřelé! Korelace s Interceptem jinak není potřeba řešit snad vůbec. Korelace koeficientů “slope” mohou signalizovat potenciální kolinearitu (o které B. Winter v knize mluví, ale ne v kapitole lmer), ovšem také to není přímo její důkaz (jestli se nepletu).
Nejvíce mi pomohla diskuze na této adrese: https://stats.stackexchange.com/questions/57240/how-do-i-interpret-the-correlations-of-fixed-effects-in-my-glmer-output
The “correlation of fixed effects” output doesn’t have the intuitive meaning that most would ascribe to it. Specifically, is not about the correlation of the variables. It is in fact about the expected correlation of the regression coefficients. Although this may speak to multicollinearity it does not necessarily. In this case it is telling you that if you did the experiment again and it so happened that the coefficient 1 got smaller, it is likely that so too will would the coefficient 2.
Přeloženo do češtiny (jestli to dobře chápu a jestli má autor pravdu), tak se nejedná přímo o korelaci proměnných (tam by čísla, která nám vycházejí, znamenala skutečný průšvih - vysokou prokorelovanost proměnných). Jedná se o čísla, která se snaží odhadnout, co by se dělo s koeficienty modelu při replikaci experimentu (tedy stejný základní soubor, ale jiný náhodný výběr). Konkrétně pro náš případ ta vysoká záporná korelace Intercept a Age znamená, že lze očekávat, že když by vyšel vyšší koeficient Intercept, vyšel by nižší koeficient Age. Což dává smysl - pokud jsou celkově vyšší hlasy, mají více kam klesat, takže s rostoucím věkem mohou rychleji padat.
Tedy, jak já si to pro sebe vysvětluji, o provázanosti proměnných tyto korelace vypovídají, ale není to o vysvětlení celkové variability dat, jak jsme běžně zvyklí počítat r^2. Takže když tato čísla vycházejí třeba 0.7 apod., není to ještě žádný průšvih. A jak jsem psal, korelace s Interceptem vůbec neřeším. Pokud je fixních jevů více, je určitě dobré si ještě před voláním lmer() samostatně napočítat klasické korelace mezi příslušnými sloupci v našich datech. Tam by vyšší korelace určitě znamenaly problém kolinearity.
# Varying slopes z těchto dat nelze modelovat - 1 člověk je tam sice vícekrát, ale stále se stejným Age.
# Kdybychom dělali longitudiální studii, bylo by samozřejmě vše jinak...
coef(fit) # různé intercepty napříč mluvčími ID.## $ID
## (Intercept) Age SexM Age:SexM
## HC101 19.02953 -0.09974722 -11.61332 0.04737633
## HC103 14.74737 -0.09974722 -11.61332 0.04737633
## HC104 14.43228 -0.09974722 -11.61332 0.04737633
## HC105 17.92157 -0.09974722 -11.61332 0.04737633
## HC106 17.08391 -0.09974722 -11.61332 0.04737633
## HC107 17.43002 -0.09974722 -11.61332 0.04737633
## HC109 17.71729 -0.09974722 -11.61332 0.04737633
## HC110 13.67428 -0.09974722 -11.61332 0.04737633
## HC112 16.20836 -0.09974722 -11.61332 0.04737633
## HC113 18.98027 -0.09974722 -11.61332 0.04737633
## HC114 16.29761 -0.09974722 -11.61332 0.04737633
## HC115 11.99144 -0.09974722 -11.61332 0.04737633
## HC117 15.90071 -0.09974722 -11.61332 0.04737633
## HC118 14.37251 -0.09974722 -11.61332 0.04737633
## HC120 16.58171 -0.09974722 -11.61332 0.04737633
## HC121 15.55792 -0.09974722 -11.61332 0.04737633
## HC122 19.68960 -0.09974722 -11.61332 0.04737633
## HC123 21.57184 -0.09974722 -11.61332 0.04737633
## HC124 17.34178 -0.09974722 -11.61332 0.04737633
## HC125 18.74556 -0.09974722 -11.61332 0.04737633
## HC126 17.73116 -0.09974722 -11.61332 0.04737633
## HC127 13.22562 -0.09974722 -11.61332 0.04737633
## HC128 17.08301 -0.09974722 -11.61332 0.04737633
## HC301 16.07329 -0.09974722 -11.61332 0.04737633
## HC302 18.27563 -0.09974722 -11.61332 0.04737633
## HC303 12.55271 -0.09974722 -11.61332 0.04737633
## HC304 18.28014 -0.09974722 -11.61332 0.04737633
## HC305 13.25037 -0.09974722 -11.61332 0.04737633
## HC306 17.93835 -0.09974722 -11.61332 0.04737633
## HC307 16.78801 -0.09974722 -11.61332 0.04737633
##
## attr(,"class")
## [1] "coef.mer"P-hodnoty
fit_full <- lmer(F0 ~ Age*Sex + (1|ID), data, REML = FALSE) # teprve pro vyhodnocení likelihood ratio se REML vypíná
fit_0agesex <- lmer(F0 ~ Age+Sex + (1|ID), data, REML = FALSE)
anova(fit_full, fit_0agesex) # Má interakce význam? Zdá se, že ne.## Data: data
## Models:
## fit_0agesex: F0 ~ Age + Sex + (1 | ID)
## fit_full: F0 ~ Age * Sex + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## fit_0agesex 5 1053.1 1072.5 -521.53 1043.1
## fit_full 6 1054.4 1077.7 -521.18 1042.4 0.6953 1 0.4044 # Pozn. Pokud by měla, bralo by se to automaticky
# jako prokázání toho, že i Age a Sex dohromady jsou významné, což může znamenat 2 věci:
# a) byly by významné i bez interakce, ale interakce ještě výrazně sníží nevysvětlený rozptyl
# b) samy o sobě nevýznamné (např. zcela opačné trendy F0 ~ Age pro obě hodnoty Sex),
# ale dohromady teprve začnou dávat smysl a výrazně sníží rozptyl
# (viz pokus.R, kde modré a červené mají zcela opačný směr, takže x samo o sobě nevysvětlí nic,
# ale interakce s class ano, tedy je jasné, že x význam má, jen je potřeba ho řešit pro každé
# class samostatně)
# Když tedy interakce nevyšla, zkusíme porovnat model bez interakce s modelem bez Age
fit_0age <- lmer(F0 ~ Sex + (1|ID), data, REML = FALSE)
anova(fit_0agesex, fit_0age)## Data: data
## Models:
## fit_0age: F0 ~ Sex + (1 | ID)
## fit_0agesex: F0 ~ Age + Sex + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## fit_0age 4 1057.3 1072.8 -524.64 1049.3
## fit_0agesex 5 1053.1 1072.5 -521.53 1043.1 6.2262 1 0.01259 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # Závěr: Age vychází významně. Interakce v těchto datech už ale nic víc nezpřesní =>
# nelze vyloučit, že mezi ženami a muži ve sklonu není rozdíl.Konfidenční intervaly pro koeficienty modelu
library(emmeans)
prumery <- emmeans(fit, ~ Sex | Age) # porovnání středů F0 dle Sex pro průměrný Age
prumery # - to není až tak překvapivé (že jsou v průměru ženy úplně jinde než muži)## Age = 48.6:
## Sex emmean SE df lower.CL upper.CL
## F 11.70 0.726 26 10.21 13.19
## M 2.39 0.553 26 1.25 3.53
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95plot(prumery)
# Pozor, u interakcí jsou interpretace rozdílů mezi skupinami obecně složitější. Zde počítáme a zobrazujeme tzv.
# main-effects, tedy rozdíl F0 pro PRŮMĚRNÝ VĚK. Pro jiný věk budou rozdíly F0 jiné a v obecném případě se může stát,
# že někde kladné, někde i záporné (když by byly rozdíly mezi slope výraznější, tak se někde přímky klidně protnou)
sklony <- emtrends(fit, ~ Sex, var = "Age") # porovnání slopes Age dle Sex
sklony # - toto je zajímavé, jde o naši hypotézu, že s věkem se mění f0## Sex Age.trend SE df lower.CL upper.CL
## F -0.0997 0.0438 26 -0.190 -0.00979
## M -0.0524 0.0420 26 -0.139 0.03402
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95plot(sklony, comparisons = TRUE)
# Takže ano, je to ve shodě s p-hodnotou významu interakce. Nelze prokázat, že je rozdíl mezi oběma sklony, protože se
# konfidenční intervaly překrývají. Ale na druhou stranu, pokud to vezmeme pro jednotlivá pohlaví, tak u žen je sklon
# prokazatelně záporný, zatímco u mužů není. Takže zní-li naše hypotéza, že s Age se mění F0, tak pro Sex == F ano,
# což je důležité. Skupina žen tedy přispěla k celkovému vyhodnocení významnosti Age, ale (viz i výše) skrz nejistoty
# ještě nelze prokázat, že je rozdíl mezi sklonem žen a mužů v závislosti na Age.
# P-hodnota významu interakce ovšem řešila hypotézu rozdílnosti sklonů M vs. F a díky šířce šipek
# dochází k jejich překryvu a rozdíl nelze prokázat.# profesionální nástroj, který rovnou zjistí p-hodnoty všech fixních efektů
library(afex)# takhle to odpovídá tomu mému
# 3. řádek Age:Sex ... vyjmutí interakce z plného modelu
# 2. řádek Sex ... vyjmutí Sex z modelu bez interakce
# 1. řádek Age ... vyjmutí Age z modelu bez interakce
mixed(F0 ~ Age*Sex + (1|ID), data, method = "LRT", type = 2)## Contrasts set to contr.sum for the following variables: Sex, ID## Numerical variables NOT centered on 0: Age
## If in interactions, interpretation of lower order (e.g., main) effects difficult.## REML argument to lmer() set to FALSE for method = 'PB' or 'LRT'## Mixed Model Anova Table (Type 2 tests, LRT-method)
##
## Model: F0 ~ Age * Sex + (1 | ID)
## Data: data
## Df full model(s): 5, 5, 6
## Effect df Chisq p.value
## 1 Age 1 6.23 * .013
## 2 Sex 1 47.74 *** <.001
## 3 Age:Sex 1 0.70 .404
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1 # pro vysvětlení type II a III viz
# https://www.r-bloggers.com/2011/03/anova-%E2%80%93-type-iiiiii-ss-explained/
mixed(F0 ~ Age*Sex + (1|ID), data, method = "LRT", type = 3) # defaultní hodnota type, tedy není nutné zadávat## Contrasts set to contr.sum for the following variables: Sex, ID## Numerical variables NOT centered on 0: Age
## If in interactions, interpretation of lower order (e.g., main) effects difficult.## REML argument to lmer() set to FALSE for method = 'PB' or 'LRT'## Mixed Model Anova Table (Type 3 tests, LRT-method)
##
## Model: F0 ~ Age * Sex + (1 | ID)
## Data: data
## Df full model: 6
## Effect df Chisq p.value
## 1 Age 1 6.50 * .011
## 2 Sex 1 13.04 *** <.001
## 3 Age:Sex 1 0.70 .404
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1# Jak se dospělo k tomu 0.011?
# Sice si zapnul contr.sum (já v příkladech nahoře používám defaultní contr.treatment), ale to bych uměl, viz níže.
fit1 <- lmer(F0 ~ Age + Sex + Age:Sex + (1|ID), data, REML = FALSE, contrasts = list(Sex = "contr.sum"))
fit2 <- lmer(F0 ~ Sex + Age:Sex + (1|ID), data, REML = FALSE, contrasts = list(Sex = "contr.sum"))
anova(fit1, fit2)## Data: data
## Models:
## fit1: F0 ~ Age + Sex + Age:Sex + (1 | ID)
## fit2: F0 ~ Sex + Age:Sex + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## fit1 6 1054.4 1077.7 -521.18 1042.4
## fit2 6 1054.4 1077.7 -521.18 1042.4 0 0# jenže toto nevede na rozdílné rozptyly, vysvětlení asi viz následující informace:
# The model comparisons for main effects with SS type III cannot be done using anova() due to the violation of marginality
# principle - anova() automatically includes all lower-order terms for the interactions it finds.
# Nicméně: podle různých rad na fórech se má type 3 používat jen v případě, že interakce je významná.
# Když není, má se použít type 2, takže můj postup nahoře je asi správně a zdejší snaha o pochopení
# 0.011 je sice zajímavá, ale ne nezbytná.
# A můj názor je, že když interakce je významná, pak stačí jen p-hodnota pro tu interakci a na zbytek se
# podívám v emmeans, protože mě stejně zajímá efekt jednotlivých hodnot faktorů.“Ruční” konstrukce středových přímek z koeficientů obdrženého lmer modelu
Vzhledem k interakci jsou různé sklony. Přímky by měly odpovídat těm v ggplot ze zprůměrovaných dat.
summary(fit)## Linear mixed model fit by REML ['lmerMod']
## Formula: F0 ~ Age * Sex + (1 | ID)
## Data: data
##
## REML criterion at convergence: 1049.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5632 -0.5457 0.0062 0.5464 4.6027
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 5.7441 2.3967
## Residual 0.7335 0.8564
## Number of obs: 360, groups: ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 16.54913 2.24539 7.370
## Age -0.09975 0.04377 -2.279
## SexM -11.61332 3.08625 -3.763
## Age:SexM 0.04738 0.06068 0.781
##
## Correlation of Fixed Effects:
## (Intr) Age SexM
## Age -0.946
## SexM -0.728 0.688
## Age:SexM 0.682 -0.721 -0.955fixef(fit) # takto se dostaneme k jednotlivým koeficientům ve formě vektoru## (Intercept) Age SexM Age:SexM
## 16.54912894 -0.09974722 -11.61331863 0.04737633names(fixef(fit)) # s těmito jmény## [1] "(Intercept)" "Age" "SexM" "Age:SexM"summary(fit)$coefficients # nebo ve formě tabulky takto## Estimate Std. Error t value
## (Intercept) 16.54912894 2.24539493 7.3702531
## Age -0.09974722 0.04376579 -2.2791137
## SexM -11.61331863 3.08625417 -3.7629171
## Age:SexM 0.04737633 0.06067835 0.7807782Rovnice vycházejí
F: F0 = 16.5491289 + -0.0997472 \(*\)
Age + -11.6133186 \(*\) 0 + 0.0473763
\(*\) Age \(*\) 0 =
= 16.5491289 + -0.0997472 \(*\)
Age
M: F0 = 16.5491289 + -0.0997472 \(*\)
Age + -11.6133186 \(*\) 1 + 0.0473763
\(*\) Age \(*\) 1 =
= 4.9358103 + -0.0523709 \(*\) Age
par(mfrow = c(1, 3))
library(scales)
sex2col <- function(v) { return(ifelse(v == "F", hue_pal()(2)[1], hue_pal()(2)[2])) } # show_col(hue_pal()(2))
gg <- ggplot(dataAV, aes(x = Age, y = F0, color = Sex)) + geom_point() + geom_smooth(method = "lm")
plot.new(); vps <- baseViewports(); pushViewport(vps$figure); print(gg, vp = plotViewport(c(1, 0, 2, 0)))
plot(F0 ~ Age, data, pch = 16, cex = 1, col = sex2col(data$Sex), main = "Včetně interakce")
x <- c(min(data$Age), max(data$Age)) # dva body z Age pro konstrukci přímky
y <- fixef(fit)["(Intercept)"] + fixef(fit)["Age"] * x
lines(x, y, col = sex2col("F"), lwd = 2)
y <- fixef(fit)["(Intercept)"] + fixef(fit)["SexM"] + (fixef(fit)["Age"] + fixef(fit)["Age:SexM"]) * x
lines(x, y, col = sex2col("M"), lwd = 2)
par(mfrow = c(1, 1))
emmip(fit, Sex ~ Age, cov.reduce = range) # To samé pomocí emmeans package