J'essaie de reproduire plusieurs tests d'interaction entre les deux lm
et lmer
sur 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' lm
approche 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 lm
estimations ne correspond exactement à lmer
celles. 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 lm
résultats).
Prime:
J'ai remarqué que l' t value
associé à 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
lm
modèle) aveclmer
, et aussi savoir quelles sont les analyses correcteslmer
pour ce type de données.lm
; Je soupçonne que c'est pourquoi la statistique t est environ deux fois plus petitelmer
. 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.Réponses:
Étrange, quand j'utilise votre dernier modèle, je trouve une correspondance parfaite, pas une correspondance étroite:
la source
summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c|id), XL, control=lmerControl(check.nobs.vs.nRE="ignore")) )$coefficients
revientt = 46.5148684
pour moi. Pourrait être un problème de version? J'utiliseR version 3.5.3 (2019-03-11)
etlmerTest 3.1-0
.