Medidas repetidas anova: lm vs lmer

10

Estou tentando reproduzir vários testes de interação entre ambos lme lmerem medidas repetidas (2x2x2). O motivo pelo qual desejo comparar os dois métodos é porque o GLM do SPSS para medidas repetidas produz exatamente os mesmos resultados da lmabordagem apresentada aqui, então, no final, quero comparar o SPSS vs o R-lmer. Até agora, apenas consegui reproduzir (de perto) algumas dessas interações.

Abaixo, você encontrará um script para ilustrar melhor meu argumento:

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 você pode ver acima, nenhuma das lmestimativas corresponde exatamente lmeràquelas. Embora alguns dos resultados sejam muito semelhantes e possam diferir apenas devido a razões numéricas / computacionais. A diferença entre os dois métodos de estimativa é especialmente grande para o X2:X3 2-way interaction (for all the data).

Minha pergunta é se existe uma maneira de obter exatamente os mesmos resultados com ambos os métodos e se existe uma maneira correta de executar as análises lmer(embora isso possa não corresponder aos lmresultados).


Bônus:

Percebi que o t valueassociado à interação de três vias é afetado pela maneira como os fatores são codificados, o que me parece muito estranho:

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
esteira
fonte
11
+1 porque parece interessante, mas não tenho idéia do que você está fazendo aqui :) Você pode explicar em palavras ou matemática por que essas chamadas lm e lmer devem produzir os mesmos coeficientes? E qual é a lógica por trás de todo esse exercício?
Ameba
@amoeba Atualizei meu post para esclarecer o objetivo deste post. Basicamente, desejo reproduzir os resultados do SPSS (que podem ser traduzidos em um lmmodelo) lmere também saber quais são as análises corretas lmer para esse tipo de dados.
mat
O motivo da grande discrepância no caso da interação bidirecional para os dados completos é que você tem 2 pontos de dados por combinação de parâmetros. A intuição é que o tamanho efetivo da amostra para um modelo misto é 2x menor do que para lm; Eu suspeito que é por isso que a estatística t é aproximadamente duas vezes menor lmer. Você provavelmente seria capaz de observar o mesmo fenômeno usando um design 2x2 mais simples e observando os principais efeitos, sem se preocupar com 2x2x2 e interações complicadas.
Ameba

Respostas:

3

Estranho, quando uso seu último modelo, encontro uma combinação perfeita, não uma combinação próxima:

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
fonte
11
Só para deixar claro, a qual modelo você está se referindo?
mat
resumo (lmer (Y ~ X1c X2c X3c + (X1c X2c X3c | id), XL, control = lmerControl (check.nobs.vs.nRE = "ignore")))
user244839
Isso é realmente muito estranho! summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c|id), XL, control=lmerControl(check.nobs.vs.nRE="ignore")) )$coefficientsretorna t = 46.5148684para mim. Pode ser um problema de versão? Eu estou usando R version 3.5.3 (2019-03-11)e lmerTest 3.1-0.
esteira
Eu tenho as mesmas versões R & lmerTest que @mat e obtenho os mesmos resultados que elas (embora com muitos avisos - falha na convergência, etc).
mkt - Reinstala Monica
11
@mat Talvez eu não tenha esclarecido - estou obtendo os mesmos resultados que você! Acho que você provavelmente está certo que o user244839 está usando uma versão diferente da nossa.
mkt - Reinstala Monica