Mesures répétées anova: lm vs lmer

10

J'essaie de reproduire plusieurs tests d'interaction entre les deux lmet lmersur des mesures répétées (2x2x2). La raison pour laquelle je veux comparer les deux méthodes est que le GLM de SPSS pour les mesures répétées donne exactement les mêmes résultats que l' lmapproche présentée ici, donc à la fin je veux comparer SPSS vs R-lmer. Jusqu'à présent, je n'ai réussi à reproduire (de près) que certaines de ces interactions.

Vous trouverez ci-dessous un script pour mieux illustrer mon propos:

library(data.table)
library(tidyr)
library(lmerTest)
library(MASS)

set.seed(1)

N     <- 100 # number of subjects
sigma <- 1   # popuplation sd
rho   <- .6  # correlation between variables

# X1:   a  a  a  a  b  b  b  b
# X2:   a  a  b  b  a  a  b  b
# X3:   a  b  a  b  a  b  a  b
mu <- c(5, 3, 3, 5, 3, 5, 5, 3) # means

# Simulate the data
sigma.mat <- rep(sigma, length(mu))
S <- matrix(sigma.mat, ncol = length(sigma.mat), nrow = length(sigma.mat))
Sigma <- t(S) * S * rho  
diag(Sigma) <- sigma**2
X <- data.table( mvrnorm(N, mu, Sigma) )
setnames(X, names(X), c("aaa", "aab", "aba", "abb", "baa", "bab", "bba", "bbb"))
X[, id := 1:.N]

# Long format
XL <- data.table( gather(X, key, Y, aaa:bbb) )
XL[, X1 := substr(key, 1, 1)]
XL[, X2 := substr(key, 2, 2)]
XL[, X3 := substr(key, 3, 3)]

# Recode long format (a = +1; b = -1)
XL[, X1c := ifelse(X1 == "a", +1, -1)]
XL[, X2c := ifelse(X2 == "a", +1, -1)]
XL[, X3c := ifelse(X3 == "a", +1, -1)]


### Composite scores to be used with lm
# X2:X3 2-way interaction (for half the data; i.e. when X1 == "a")
X[, X1a_X2.X3 := (aaa - aab) - (aba - abb)]

# X2:X3 2-way interaction (for all the data)
X[, aa := (aaa + baa) / 2]
X[, ab := (aab + bab) / 2]
X[, ba := (aba + bba) / 2]
X[, bb := (abb + bbb) / 2]
X[, X2.X3 := (aa - ab) - (ba - bb)]

# X1:X2:X3 3-way interaction (for all the data)
X[, X1.X2.X3 := ( (aaa - aab) - (aba - abb) ) - ( (baa - bab) - (bba - bbb) )]


### Fit models
# X2:X3 2-way interaction (for half the data; i.e. when X1 == "a")
summary( lm(X1a_X2.X3 ~ 1, X) ) # t = 34.13303
summary( lmer(Y ~ X2c*X3c + (X2c+X3c|id), XL[X1 == "a"]) ) # t = 34.132846  close match
summary( lmer(Y ~ X2c*X3c + (X2c+X3c||id), XL[X1 == "a"]) ) # t = 34.134624  close match

# X2:X3 2-way interaction (for all the data) 
summary( lm(X2.X3 ~ 1, X) ) # t = 0.3075025
summary( lmer(Y ~ X2c*X3c + (X2c+X3c|id), XL) ) # t = 0.1641932
summary( lmer(Y ~ X2c*X3c + (X2c+X3c||id), XL) ) # t = 0.1640710
summary( lmer(Y ~ X2c*X3c + (X2c*X3c|id), XL) ) # t = 0.1641765
anova(   lmer(Y ~ X2c*X3c + (X2c*X3c|id), XL), ddf = "Kenward-Roger" ) # t = 0.1643168
summary( lmer(Y ~ X2c*X3c + (X2c*X3c|id), XL, REML = FALSE) ) # t = 0.1645303
summary( lmer(Y ~ X2c*X3c + (X2c*X3c||id), XL) ) # t = 0.1640704

# X1:X2:X3 3-way interaction (for all the data)
summary( lm(X1.X2.X3 ~ 1, X) ) # t = 46.50177
summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c - X1c:X2c:X3c|id), XL) ) # t = 49.0317599
anova(   lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c - X1c:X2c:X3c|id), XL), ddf = "Kenward-Roger" ) # t = 49.03176
summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c - X1c:X2c:X3c|id), XL, REML = FALSE) ) # t = 49.2677606
summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c - X1c:X2c:X3c||id), XL) ) # t = 46.5193774 close match
summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c|id), XL) ) # unidentifiable
summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c|id), XL,
              control = lmerControl(check.nobs.vs.nRE="ignore")) ) # t = 46.5148684 close match

Comme vous pouvez le voir ci-dessus, aucune des lmestimations ne correspond exactement à lmercelles. Bien que certains résultats soient très similaires et peuvent différer uniquement pour des raisons numériques / informatiques. L'écart entre les deux méthodes d'estimation est particulièrement important pour le X2:X3 2-way interaction (for all the data).

Ma question est de savoir s'il existe un moyen d'obtenir les mêmes résultats exacts avec les deux méthodes, et s'il existe un moyen correct d'effectuer les analyses avec lmer(bien qu'il puisse ne pas correspondre aux lmrésultats).


Prime:

J'ai remarqué que l' t valueassocié à l'interaction à 3 voies est affecté par la façon dont les facteurs sont codés, ce qui me semble très étrange:

summary( lmer(Y ~ X1*X2*X3 + (X1*X2*X3 - X1:X2:X3||id), XL) ) # t = 48.36
summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c - X1c:X2c:X3c||id), XL) ) # t = 56.52
tapis
la source
1
+1 parce que cela semble intéressant mais je n'ai aucune idée de ce que vous faites ici :) Pouvez-vous expliquer en mots ou en mathématiques pourquoi ces appels lm et lmer devraient produire les mêmes coefficients? Et quelle est la logique derrière tout cet exercice?
amoeba
@amoeba J'ai mis à jour mon message pour clarifier le but de ce message. Fondamentalement, je veux reproduire les résultats de SPSS (qui peuvent être traduits en lmmodèle) avec lmer, et aussi savoir quelles sont les analyses correctes lmer pour ce type de données.
mat
La raison de l'écart important en cas d'interaction bidirectionnelle pour les données complètes est que vous avez 2 points de données par combinaison de paramètres. L'intuition est que la taille effective de l'échantillon pour un modèle mixte est 2x plus petite que pour lm; Je soupçonne que c'est pourquoi la statistique t est environ deux fois plus petite lmer. Vous seriez probablement en mesure d'observer le même phénomène en utilisant une conception 2x2 plus simple et en regardant les effets principaux, sans vous soucier des interactions 2x2x2 et compliquées.
amoeba

Réponses:

3

Étrange, quand j'utilise votre dernier modèle, je trouve une correspondance parfaite, pas une correspondance étroite:

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  3.91221    0.07242 99.00001  54.025   <2e-16 ***
X1c          0.03277    0.05006 99.00000   0.655    0.514    
X2c         -0.04836    0.04644 99.00000  -1.042    0.300    
X3c          0.04248    0.05009 99.00001   0.848    0.398    
X1c:X2c      0.08370    0.08747 98.99998   0.957    0.341    
X1c:X3c     -0.07025    0.08895 98.99994  -0.790    0.432    
X2c:X3c     -0.02957    0.09616 99.00000  -0.308    0.759    
X1c:X2c:X3c -8.14099    0.17507 99.00003 -46.502   <2e-16 ***
user244839
la source
1
Pour être clair, de quel modèle parlez-vous?
mat
résumé (lmer (Y ~ X1c X2c X3c + (X1c X2c X3c | id), XL, control = lmerControl (check.nobs.vs.nRE = "ignore"))))
user244839
C'est en effet très étrange! summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c|id), XL, control=lmerControl(check.nobs.vs.nRE="ignore")) )$coefficientsrevient t = 46.5148684pour moi. Pourrait être un problème de version? J'utilise R version 3.5.3 (2019-03-11)et lmerTest 3.1-0.
mat
J'ai les mêmes versions de R & lmerTest que @mat et j'obtiens les mêmes résultats qu'eux (bien qu'avec de nombreux avertissements - échec de la convergence, etc.).
mkt
1
@mat Peut-être que je n'étais pas clair - j'obtiens les mêmes résultats que vous! Je pense que vous avez probablement raison de dire que user244839 utilise une version différente de la nôtre.
mkt