Medidas repetidas anova: lm vs lmer

10

Estoy tratando de reproducir varias pruebas de interacción entre ambas lmy lmeren medidas repetidas (2x2x2). La razón por la que quiero comparar ambos métodos es porque el GLM de SPSS para medidas repetidas produce exactamente los mismos resultados que el lmenfoque presentado aquí, por lo que al final quiero comparar SPSS vs R-lmer. Hasta ahora, solo he logrado reproducir (de cerca) algunas de estas interacciones.

A continuación encontrará un script para ilustrar mejor mi punto:

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

Como puede ver desde arriba, ninguna de las lmestimaciones coincide exactamente con las lmermismas. Aunque algunos de los resultados son muy similares y pueden diferir solo debido a razones numéricas / computacionales. La brecha entre ambos métodos de estimación es especialmente grande para X2:X3 2-way interaction (for all the data).

Mi pregunta es si hay una manera de obtener exactamente los mismos resultados con ambos métodos, y si hay una forma correcta de realizar los análisis lmer(aunque puede que no coincida con los lmresultados).


Prima:

Me di cuenta de que lo t valueasociado con la interacción de 3 vías se ve afectado por la forma en que se codifican los factores, lo que me parece muy extraño:

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
estera
fuente
1
+1 porque parece interesante, pero no tengo idea de lo que estás haciendo aquí :) ¿Puedes explicar con palabras o matemáticas por qué estas llamadas lm y lmer deberían producir los mismos coeficientes? ¿Y cuál es la lógica detrás de todo este ejercicio?
ameba
@amoeba Actualicé mi publicación para aclarar el propósito de esta publicación. Básicamente, quiero reproducir los resultados de SPSS (que se puede traducir en un lmmodelo) con lmer, y también saber cuáles son los análisis correctos lmer para este tipo de datos.
mat
La razón de la gran discrepancia en el caso de la interacción bidireccional para los datos completos es que tiene 2 puntos de datos por combinación de parámetros. La intuición es que el tamaño de muestra efectivo para un modelo mixto es 2 veces menor que para lm; Sospecho que es por eso que la estadística t es aproximadamente dos veces más pequeña lmer. Probablemente podría observar el mismo fenómeno utilizando un diseño 2x2 más simple y observando los efectos principales, sin molestarse con las interacciones complicadas y 2x2x2.
ameba

Respuestas:

3

Extraño, cuando uso su último modelo, encuentro una combinación perfecta, no una coincidencia cercana:

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 ***
usuario244839
fuente
1
Para ser claros, ¿a qué modelo te refieres?
mat
summary (lmer (Y ~ X1c X2c X3c + (X1c X2c X3c | id), XL, control = lmerControl (check.nobs.vs.nRE = "ignore")))
user244839
¡Esto es realmente muy extraño! summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c|id), XL, control=lmerControl(check.nobs.vs.nRE="ignore")) )$coefficientsvuelve t = 46.5148684por mi ¿Podría ser un problema de versión? Estoy usando R version 3.5.3 (2019-03-11)y lmerTest 3.1-0.
mat
Tengo las mismas versiones de R & lmerTest que @mat y obtengo los mismos resultados que ellas (aunque con muchas advertencias: fallo de convergencia, etc.).
mkt - Restablecer Mónica
1
@mat Quizás no estaba claro: ¡estoy obteniendo los mismos resultados que tú! Creo que probablemente tengas razón en que user244839 está usando una versión diferente a la nuestra.
mkt - Restablecer a Mónica el