Powtarzane miary anova: lm vs lmer

10

Próbuję odtworzyć kilka testów interakcji między obiema lmi lmerpowtarzanymi pomiarami (2x2x2). Powodem, dla którego chcę porównać obie metody, jest to, że GLM SPSS dla powtarzanych pomiarów daje dokładnie takie same wyniki, jak lmprzedstawione tutaj podejście, więc na koniec chcę porównać SPSS z R-lmer. Do tej pory udało mi się tylko odtworzyć (ściśle) niektóre z tych interakcji.

Poniżej znajduje się skrypt, który lepiej zilustruje mój punkt widzenia:

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

Jak widać z góry, żadna z tych lmszacunków nie jest dokładnie taka sama lmer. Chociaż niektóre wyniki są bardzo podobne i mogą się różnić tylko ze względu na przyczyny numeryczne / obliczeniowe. Różnica między obiema metodami szacowania jest szczególnie duża w przypadku X2:X3 2-way interaction (for all the data).

Moje pytanie brzmi, czy istnieje sposób na uzyskanie dokładnie takich samych wyników za pomocą obu metod i czy istnieje prawidłowy sposób wykonywania analiz lmer(chociaż może nie być zgodny z lmwynikami).


Premia:

Zauważyłem, że na t valuesposób interakcji z trójdrożnym wpływa sposób kodowania czynników, co wydaje mi się bardzo dziwne:

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
mata
źródło
1
+1, ponieważ wygląda interesująco, ale nie mam pojęcia, co tutaj robisz :) Czy możesz wyjaśnić słowami lub matematyką, dlaczego te wywołania lm i lmer powinny dawać te same współczynniki? A jaka jest logika tego całego ćwiczenia?
ameba
@amoeba Zaktualizowałem swój post, aby wyjaśnić cel tego postu. Zasadniczo chcę odtworzyć wyniki z SPSS (które można przetłumaczyć na lmmodel) lmer, a także wiedzieć, jakie są prawidłowe lmer analizy dla tego rodzaju danych.
mat
Powodem dużej rozbieżności w przypadku dwukierunkowej interakcji dla pełnych danych jest to, że masz 2 punkty danych na kombinację parametrów. Intuicja jest taka, że ​​efektywny rozmiar próbki dla modelu mieszanego jest 2x mniejszy niż dla lm; Podejrzewam, że właśnie dlatego statystyka t jest około dwa razy mniejsza lmer. Prawdopodobnie byłbyś w stanie zaobserwować to samo zjawisko, stosując prostszą konstrukcję 2x2 i patrząc na główne efekty, bez zawracania głowy 2x2x2 i skomplikowanymi interakcjami.
ameba

Odpowiedzi:

3

Dziwne, kiedy używam twojego ostatniego modelu, znajduję idealne dopasowanie, a nie ścisłe dopasowanie:

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 ***
użytkownik244839
źródło
1
Żeby było jasne, do którego modelu się odwołujesz?
mat
podsumowanie (lmer (Y ~ X1c X2c X3c + (X1c X2c X3c | id), XL, control = lmerControl (check.nobs.vs.nRE = „ignore”)))
user244839 16.04.19
To jest naprawdę bardzo dziwne! summary( lmer(Y ~ X1c*X2c*X3c + (X1c*X2c*X3c|id), XL, control=lmerControl(check.nobs.vs.nRE="ignore")) )$coefficientswraca t = 46.5148684po mnie. Może to być problem z wersją? Używam R version 3.5.3 (2019-03-11)i lmerTest 3.1-0.
mat
Mam te same wersje R & lmerTest co @mat i otrzymuję takie same wyniki jak one (choć z wieloma ostrzeżeniami - brak zbieżności itp.).
mkt - Przywróć Monikę
1
@mat Być może nie byłem pewien - otrzymuję takie same wyniki jak Ty! Myślę, że prawdopodobnie masz rację, że user244839 używa innej wersji niż my.
mkt - Przywróć Monikę