Comment définir des contrastes personnalisés avec lmer dans R

9

J'utilise lmer dans R pour vérifier l'effet de condition ( cond) sur un résultat. Voici quelques données composées, où s est l'identifiant du sujet et a, bet csont des conditions.

library("tidyr")
library("dplyr")
set.seed(123)
temp <- data.frame(s = paste0("S", 1:30), 
                   a = rnorm(30, -2, 1), 
                   b = rnorm(30, -3, 1), 
                   c = rnorm(30, -4, 1)) 

Je voudrais comparer

  1. niveau aà la moyenne des niveaux bet cet
  2. niveau bà niveau c.

Ma question est la suivante: comment puis-je définir les contrastes de manière à ce que l'ordonnée à l'origine reflète la moyenne des trois conditions et que les deux estimations calculées reflètent directement les différences définies en 1. et 2.?

J'ai essayé avec

c1 <- cbind(c(-0.5, 0.25, 0.25), c(0, -0.5, 0.5))
gather(temp, cond, result, a, b, c) %>%
  lmer(result ~ cond + (1|s), data = ., contrasts = list(cond = c1))

cond2semble être OK, mais cond1ne l'est pas.

Suite Comment interpréter ces contrastes personnalisés? , J'ai essayé d'utiliser l'inverse généralisé à la place, mais ces estimations n'ont pas de sens non plus.

c2 <- t(ginv(c1))
gather(temp, cond, result, a, b, c) %>%
  lmer(result ~ cond + (1|s), data = ., contrasts = list(cond = c2))

J'ai aussi essayé les contrastes Helmert, mais les moyens ne correspondent toujours pas.

gather(temp, cond, result, a, b, c) %>%
  mutate(cond = factor(cond, levels = c("c", "b", "a"))) %>%
  lmer(result ~ cond + (1|s), data = ., contrasts = list(cond = contr.helmert))

Quelle est la bonne façon de procéder?

M4RT1NK4
la source
Cela ressemble à un contraste Helmert (c est le premier niveau, puis b, puis a).
Michael M
J'ai essayé Helmert aussi, mais les chiffres ne sont pas le moyen que je recherche. J'ai modifié la question pour inclure les contrastes Helmert, merci.
M4RT1NK4

Réponses:

13

Pour les étapes suivantes, nous avons besoin de la trame de données au format long. La trame de données datcontient la variable dépendante result, le prédicteur catégorique cond(niveaux: a, b, et c), et le facteur aléatoire s.

library(tidyr)
dat <- gather(temp, cond, result, a, b, c)

Dans ce qui suit, je vais illustrer deux approches pour créer une matrice de contraste correspondant aux conditions que vous souhaitez comparer:

  1. une-b+c2
  2. b-c

Contrastes personnalisés

La matrice matcorrespond aux différences de niveau.

mat <- rbind(c(1, -0.5, -0.5),     # a vs. (b + c) / 2
             c(0, 1, -1))          # b vs. c

Pour créer la matrice de contraste réelle, nous calculons l'inverse généralisé avec ginv(à partir de MASS).

library(MASS)
cMat <- ginv(mat)
#            [,1]          [,2]
# [1,]  0.6666667 -7.130169e-17
# [2,] -0.3333333  5.000000e-01
# [3,] -0.3333333 -5.000000e-01

Cette matrice de contraste cMatpeut être utilisée dans lmer.

library(lme4)
res <- lmer(result ~ cond + (1|s), data = dat, 
            contrasts = list(cond = cMat))
coef(summary(res))    
#              Estimate Std. Error    t value
# (Intercept) -2.948115  0.0946025 -31.163182
# cond1        1.351517  0.2006822   6.734612
# cond2        1.153918  0.2317279   4.979625

Comme vous pouvez le voir, les estimations à effet fixe correspondent aux différences spécifiées ci-dessus. De plus, l'ordonnée à l'origine représente la moyenne globale.

Helmert contraste avec contr.helmert

Vous pouvez également utiliser la contr.helmertfonction intégrée pour créer la matrice de contraste.

cHelmert <- contr.helmert(3)
#   [,1] [,2]
# 1   -1   -1
# 2    1   -1
# 3    0    2

Cependant, la commande ne correspond pas à celle que vous avez spécifiée dans la question. Par conséquent, nous devons inverser l'ordre des colonnes et des lignes. La première colonne correspond à bvs aet la seconde correspond à cvs la moyenne de bet a.

cHelmert2 <- cHelmert[c(3:1), 2:1]
#   [,1] [,2]
# 3    2    0
# 2   -1    1
# 1   -1   -1

Comparez la matrice de contraste cHelmert2à cMat. Vous remarquerez que les colonnes sont des versions à l'échelle de l'autre matrice.

Le résultat lmerest de:

library(lme4)
res2 <- lmer(result ~ cond + (1|s), data = dat, 
             contrasts = list(cond = cHelmert2))
coef(summary(res2))    
#               Estimate Std. Error    t value
# (Intercept) -2.9481150 0.09460250 -31.163182
# cond1        0.4505056 0.06689407   6.734612
# cond2        0.5769590 0.11586393   4.979625

t

Sven Hohenstein
la source
Grand merci! Juste pour m'assurer que je comprends cela maintenant - si je voulais comparer le premier niveau au reste des niveaux dans une variable à 4 niveaux, ce matserait c(1, -1/3, -1/3, -1/3)? J'ai donc toujours défini les nombres tels qu'ils seraient dans la formule (a + (b + c + d) / 3), puis je les ai ginvajustés de manière appropriée afin que les coefficients reflètent directement la différence. Et lorsque vous avez changé l'ordre dans l'exemple Helmert, c'était juste pour répondre à la question? Sinon, les résultats devraient être les mêmes, quel que soit l'ordre des contrastes, non?
M4RT1NK4
@ M4RT1NK4 Votre formule et le contraste correspondant sont corrects. L'ordre des colonnes vient d'être modifié pour correspondre à l'ordre des colonnes dans la question. L'ordre des lignes est cependant important, car le premier niveau est le niveau de référence. Dans votre exemple, le niveau de référence est le troisième niveau.
Sven Hohenstein
@SvenHohenstein J'ai eu une question connexe basée sur cette réponse, ça vous dérange? stats.stackexchange.com/questions/357781/…
mat