Qual é o equivalente lme4 :: lmer de uma ANOVA de três medidas repetidas?

11

Minha pergunta é baseada nesta resposta que mostrou qual lme4::lmermodelo corresponde a uma ANOVA de medidas repetidas nos dois sentidos:

require(lme4)
set.seed(1234)
d <- data.frame(
    y = rnorm(96),
    subject = factor(rep(1:12, 4)),
    a = factor(rep(1:2, each=24)),
    b = factor(rep(rep(1:2, each=12))),
    c = factor(rep(rep(1:2, each=48))))

# standard two-way repeated measures ANOVA:
summary(aov(y~a*b+Error(subject/(a*b)), d[d$c == "1",]))

# corresponding lmer call:
anova(lmer(y ~ a*b+(1|subject) + (1|a:subject) + (1|b:subject), d[d$c == "1",]))

Minha pergunta agora é sobre como estender isso ao caso de uma ANOVA de três vias:

summary(aov(y~a*b*c+Error(subject/(a*b*c)), d))
## [...]
## Error: subject:a:b:c
##           Df Sum Sq Mean Sq F value Pr(>F)
## a:b:c      1  0.101  0.1014   0.115  0.741
## Residuals 11  9.705  0.8822

A extensão natural e suas versões não correspondem aos resultados da ANOVA:

anova(lmer(y ~ a*b*c +(1|subject) + (1|a:subject) + (1|b:subject) + (1|c:subject), d))
## [...]
## a:b:c  1 0.1014  0.1014  0.1500

anova(lmer(y ~ a*b*c +(1|subject) + (1|a:subject) + (1|b:subject) + (1|c:subject) + 
               (1|a:b:subject) + (1|a:c:subject) + (1|b:c:subject), d))
## [...]
## a:b:c  1 0.1014  0.1014  0.1539

Observe que uma pergunta muito semelhante foi feita antes . No entanto, estavam faltando dados de exemplo (que são fornecidos aqui).

r anova mixed-model repeated-measures lme4-nlme Henrik
fonte
Tem certeza de que não deseja que seu modelo de dois níveis seja y ~ a*b + (1 + a*b|subject), d[d$c == "1",]? Ou talvez eu esteja perdendo alguma coisa?
Rasmus Bååth
@ RasmusBååth Vá em frente e tente encaixá-lo, lmerirá reclamar porque os efeitos aleatórios não são mais identificados. Inicialmente, eu também pensei que esse é o modelo que eu quero, mas não é. Se você comparar o modelo mais recente que proponho para o caso de duas vias com a ANOVA padrão, você verá que os valores F correspondem exatamente . Como dito na resposta eu liguei.
Henrik
3
Para o problema de três vias, não se espera que o primeiro lmermodelo que você escreveu (que exclui as interações aleatórias de duas vias) seja equivalente a uma RM-ANOVA de três vias, mas o segundo que você escreveu (que inclui a aleatória interações bidirecionais) deve ser. Quanto ao motivo pelo qual há uma discrepância, mesmo com esse modelo, tenho um palpite sobre qual é o problema. Ir para o jantar e examinar o conjunto de dados de brinquedos um pouco mais.
Jake Westfall

Respostas:

18

A resposta direta à sua pergunta é que o último modelo que você escreveu,

anova(lmer(y ~ a*b*c +(1|subject) + (1|a:subject) + (1|b:subject) + (1|c:subject) + 
           (1|a:b:subject) + (1|a:c:subject) + (1|b:c:subject), d))

Eu acredito que é "em princípio" correto, embora seja uma parametrização estranha que nem sempre parece funcionar bem na prática real.

Quanto ao motivo pelo qual a saída que você obtém deste modelo é discrepante com a aov()saída, acho que há duas razões.

  1. Seu conjunto de dados simulado simples é patológico, pois o modelo de melhor ajuste é aquele que implica componentes de variação negativa, que os modelos mistos se encaixam lmer()(e a maioria dos outros programas de modelos mistos) não permitem.
  2. Mesmo com um conjunto de dados não patológicos, a maneira como você configurou o modelo, como mencionado acima, nem sempre parece funcionar bem na prática, embora eu deva admitir que realmente não entendo o porquê. Também é geralmente estranho na minha opinião, mas isso é outra história.

Deixe-me primeiro demonstrar a parametrização que eu prefiro no seu exemplo inicial de ANOVA de duas vias. Suponha que seu conjunto de dados desteja carregado. Seu modelo (observe que eu mudei de códigos fictícios para códigos de contraste) era:

options(contrasts=c("contr.sum","contr.poly"))
mod1 <- lmer(y ~ a*b+(1|subject) + (1|a:subject) + (1|b:subject),
         data = d[d$c == "1",])
anova(mod1)
# Analysis of Variance Table
#     Df  Sum Sq Mean Sq F value
# a    1 2.20496 2.20496  3.9592
# b    1 0.13979 0.13979  0.2510
# a:b  1 1.23501 1.23501  2.2176

que funcionou bem aqui, pois correspondia à aov()saída. O modelo que eu prefiro envolve duas alterações: codificar manualmente os fatores para não trabalhar com objetos de fator R (o que eu recomendo fazer em 100% dos casos) e especificar os efeitos aleatórios de maneira diferente:

d <- within(d, {
  A <- 2*as.numeric(paste(a)) - 3
  B <- 2*as.numeric(paste(b)) - 3
  C <- 2*as.numeric(paste(c)) - 3
})
mod2 <- lmer(y ~ A*B + (1|subject)+(0+A|subject)+(0+B|subject),
             data = d[d$c == "1",])
anova(mod2)
# Analysis of Variance Table
# Df  Sum Sq Mean Sq F value
# A    1 2.20496 2.20496  3.9592
# B    1 0.13979 0.13979  0.2510
# A:B  1 1.23501 1.23501  2.2176

logLik(mod1)
# 'log Lik.' -63.53034 (df=8)
logLik(mod2)
# 'log Lik.' -63.53034 (df=8)

As duas abordagens são totalmente equivalentes no problema simples de duas vias. Agora vamos passar para um problema de três vias. Mencionei anteriormente que o exemplo de conjunto de dados que você deu era patológico. Então, o que eu quero fazer antes de abordar seu exemplo de conjunto de dados é primeiro gerar um conjunto de dados a partir de um modelo de componentes de variação real (ou seja, onde componentes de variação diferentes de zero são incorporados ao modelo verdadeiro). Primeiro, mostrarei como minha parametrização preferida parece funcionar melhor do que a que você propôs. Em seguida, demonstrarei outra maneira de estimar os componentes de variância que não impõem que eles sejam não negativos. Então, estaremos em condições de ver o problema com o conjunto de dados de exemplo original.

O novo conjunto de dados será idêntico em estrutura, exceto que teremos 50 assuntos:

set.seed(9852903)
d2 <- expand.grid(A=c(-1,1), B=c(-1,1), C=c(-1,1), sub=seq(50))
d2 <- merge(d2, data.frame(sub=seq(50), int=rnorm(50), Ab=rnorm(50),
  Bb=rnorm(50), Cb=rnorm(50), ABb=rnorm(50), ACb=rnorm(50), BCb=rnorm(50)))
d2 <- within(d2, {
  y <- int + (1+Ab)*A + (1+Bb)*B + (1+Cb)*C + (1+ABb)*A*B +
    (1+ACb)*A*C + (1+BCb)*B*C + A*B*C + rnorm(50*2^3)
  a <- factor(A)
  b <- factor(B)
  c <- factor(C)
})

As relações F que queremos corresponder são:

aovMod1 <- aov(y ~ a*b*c + Error(factor(sub)/(a*b*c)), data = d2)
tab <- lapply(summary(aovMod1), function(x) x[[1]][1,2:4])
do.call(rbind, tab)
#                          Sum Sq Mean Sq F value
# Error: factor(sub)       439.48    8.97        
# Error: factor(sub):a     429.64  429.64  32.975
# Error: factor(sub):b     329.48  329.48  27.653
# Error: factor(sub):c     165.44  165.44  17.924
# Error: factor(sub):a:b   491.33  491.33  49.694
# Error: factor(sub):a:c   305.46  305.46  41.703
# Error: factor(sub):b:c   466.09  466.09  40.655
# Error: factor(sub):a:b:c 392.76  392.76 448.101

Aqui estão nossos dois modelos:

mod3 <- lmer(y ~ a*b*c + (1|sub)+(1|a:sub)+(1|b:sub)+(1|c:sub)+
  (1|a:b:sub)+(1|a:c:sub)+(1|b:c:sub), data = d2)
anova(mod3)
# Analysis of Variance Table
#       Df Sum Sq Mean Sq F value
# a      1  32.73   32.73  34.278
# b      1  21.68   21.68  22.704
# c      1  12.53   12.53  13.128
# a:b    1  60.93   60.93  63.814
# a:c    1  50.38   50.38  52.762
# b:c    1  57.30   57.30  60.009
# a:b:c  1 392.76  392.76 411.365

mod4 <- lmer(y ~ A*B*C + (1|sub)+(0+A|sub)+(0+B|sub)+(0+C|sub)+
  (0+A:B|sub)+(0+A:C|sub)+(0+B:C|sub), data = d2)
anova(mod4)
# Analysis of Variance Table
#       Df Sum Sq Mean Sq F value
# A      1  28.90   28.90  32.975
# B      1  24.24   24.24  27.653
# C      1  15.71   15.71  17.924
# A:B    1  43.56   43.56  49.694
# A:C    1  36.55   36.55  41.703
# B:C    1  35.63   35.63  40.655
# A:B:C  1 392.76  392.76 448.101

logLik(mod3)
# 'log Lik.' -984.4531 (df=16)
logLik(mod4)
# 'log Lik.' -973.4428 (df=16)

Como podemos ver, apenas o segundo método corresponde à saída de aov(), embora o primeiro método esteja pelo menos no estádio. O segundo método também atinge uma maior probabilidade de log. Não sei por que esses dois métodos dão resultados diferentes, pois, novamente, acho que são "em princípio" equivalentes, mas talvez seja por algumas razões numéricas / computacionais. Ou talvez eu esteja enganado e eles não sejam equivalentes nem em princípio.

Agora vou mostrar outra maneira de estimar os componentes de variação com base nas idéias tradicionais da ANOVA. Basicamente, tomaremos as equações quadradas médias esperadas para o seu projeto, substituiremos os valores observados dos quadrados médios e resolveremos os componentes de variância. Para obter os quadrados médios esperados, usaremos uma função R que escrevi há alguns anos, chamada EMS(), que está documentada AQUI . Abaixo, assumo que a função já esteja carregada.

# prepare coefficient matrix
r <- 1 # number of replicates
s <- 50 # number of subjects
a <- 2 # number of levels of A
b <- 2 # number of levels of B
c <- 2 # number of levels of C
CT <- EMS(r ~ a*b*c*s, random="s")
expr <- strsplit(CT[CT != ""], split="")
expr <- unlist(lapply(expr, paste, collapse="*"))
expr <- sapply(expr, function(x) eval(parse(text=x)))
CT[CT != ""] <- expr
CT[CT == ""] <- 0
mode(CT) <- "numeric"
# residual variance and A*B*C*S variance are confounded in
# this design, so remove the A*B*C*S variance component
CT <- CT[-15,-2]
CT
#        VarianceComponent
# Effect  e b:c:s a:c:s a:b:s a:b:c c:s b:s a:s b:c a:c a:b s   c   b   a
#   a     1     0     0     0     0   0   0   4   0   0   0 0   0   0 200
#   b     1     0     0     0     0   0   4   0   0   0   0 0   0 200   0
#   c     1     0     0     0     0   4   0   0   0   0   0 0 200   0   0
#   s     1     0     0     0     0   0   0   0   0   0   0 8   0   0   0
#   a:b   1     0     0     2     0   0   0   0   0   0 100 0   0   0   0
#   a:c   1     0     2     0     0   0   0   0   0 100   0 0   0   0   0
#   b:c   1     2     0     0     0   0   0   0 100   0   0 0   0   0   0
#   a:s   1     0     0     0     0   0   0   4   0   0   0 0   0   0   0
#   b:s   1     0     0     0     0   0   4   0   0   0   0 0   0   0   0
#   c:s   1     0     0     0     0   4   0   0   0   0   0 0   0   0   0
#   a:b:c 1     0     0     0    50   0   0   0   0   0   0 0   0   0   0
#   a:b:s 1     0     0     2     0   0   0   0   0   0   0 0   0   0   0
#   a:c:s 1     0     2     0     0   0   0   0   0   0   0 0   0   0   0
#   b:c:s 1     2     0     0     0   0   0   0   0   0   0 0   0   0   0
#   e     1     0     0     0     0   0   0   0   0   0   0 0   0   0   0

# get mean squares
(MSmod <- summary(aov(y ~ a*b*c*factor(sub), data=d2)))
#                   Df Sum Sq Mean Sq
# a                  1  429.6   429.6
# b                  1  329.5   329.5
# c                  1  165.4   165.4
# factor(sub)       49  439.5     9.0
# a:b                1  491.3   491.3
# a:c                1  305.5   305.5
# b:c                1  466.1   466.1
# a:factor(sub)     49  638.4    13.0
# b:factor(sub)     49  583.8    11.9
# c:factor(sub)     49  452.2     9.2
# a:b:c              1  392.8   392.8
# a:b:factor(sub)   49  484.5     9.9
# a:c:factor(sub)   49  358.9     7.3
# b:c:factor(sub)   49  561.8    11.5
# a:b:c:factor(sub) 49   42.9     0.9
MS <- MSmod[[1]][,"Mean Sq"]

# solve
ans <- solve(CT, MS)
cbind(rev(ans[c(grep("e",names(ans)),grep("s",names(ans)))])/
        c(1,2,2,2,4,4,4,1))
# s     1.0115549
# a:s   1.5191114
# b:s   1.3797937
# c:s   1.0441351
# a:b:s 1.1263331
# a:c:s 0.8060402
# b:c:s 1.3235126
# e     0.8765093
summary(mod4)
# Random effects:
#  Groups   Name        Variance Std.Dev.
#  sub      (Intercept) 1.0116   1.0058  
#  sub.1    A           1.5191   1.2325  
#  sub.2    B           1.3798   1.1746  
#  sub.3    C           1.0441   1.0218  
#  sub.4    A:B         1.1263   1.0613  
#  sub.5    A:C         0.8060   0.8978  
#  sub.6    B:C         1.3235   1.1504  
#  Residual             0.8765   0.9362  
# Number of obs: 400, groups:  sub, 50

Ok, agora retornaremos ao exemplo original. Os índices F que estamos tentando combinar são:

aovMod2 <- aov(y~a*b*c+Error(subject/(a*b*c)), data = d)
tab <- lapply(summary(aovMod2), function(x) x[[1]][1,2:4])
do.call(rbind, tab)
#                       Sum Sq Mean Sq F value
# Error: subject       13.4747  1.2250        
# Error: subject:a      1.4085  1.4085  1.2218
# Error: subject:b      3.1180  3.1180  5.5487
# Error: subject:c      6.3809  6.3809  5.2430
# Error: subject:a:b    1.5706  1.5706  2.6638
# Error: subject:a:c    1.0907  1.0907  1.5687
# Error: subject:b:c    1.4128  1.4128  2.3504
# Error: subject:a:b:c  0.1014  0.1014  0.1149

Aqui estão nossos dois modelos:

mod5 <- lmer(y ~ a*b*c + (1|subject)+(1|a:subject)+(1|b:subject)+
  (1|c:subject)+(1|a:b:subject)+(1|a:c:subject)+(1|b:c:subject),
  data = d)
anova(mod5)
# Analysis of Variance Table
#       Df Sum Sq Mean Sq F value
# a      1 0.8830  0.8830  1.3405
# b      1 3.1180  3.1180  4.7334
# c      1 3.8062  3.8062  5.7781
# a:b    1 1.5706  1.5706  2.3844
# a:c    1 0.9620  0.9620  1.4604
# b:c    1 1.4128  1.4128  2.1447
# a:b:c  1 0.1014  0.1014  0.1539

mod6 <- lmer(y ~ A*B*C + (1|subject)+(0+A|subject)+(0+B|subject)+
  (0+C|subject)+(0+A:B|subject)+(0+A:C|subject)+
  (0+B:C|subject), data = d)
anova(mod6)
# Analysis of Variance Table
#       Df Sum Sq Mean Sq F value
# a      1 0.8830  0.8830  1.3405
# b      1 3.1180  3.1180  4.7334
# c      1 3.8062  3.8062  5.7781
# a:b    1 1.5706  1.5706  2.3844
# a:c    1 0.9620  0.9620  1.4604
# b:c    1 1.4128  1.4128  2.1447
# a:b:c  1 0.1014  0.1014  0.1539

logLik(mod5)
# 'log Lik.' -135.0351 (df=16)
logLik(mod6)
# 'log Lik.' -134.9191 (df=16)

Nesse caso, os dois modelos produzem basicamente os mesmos resultados, embora o segundo método tenha uma probabilidade logarítmica muito ligeiramente maior. Nenhum método corresponde aov(). Mas vamos ver o que obtemos quando resolvemos os componentes de variação como fizemos acima, usando o procedimento ANOVA que não restringe os componentes de variação a não serem negativos (mas que só podem ser usados ​​em projetos balanceados sem preditores contínuos e sem dados ausentes; as suposições clássicas da ANOVA).

# prepare coefficient matrix
r <- 1 # number of replicates
s <- 12 # number of subjects
a <- 2 # number of levels of A
b <- 2 # number of levels of B
c <- 2 # number of levels of C
CT <- EMS(r ~ a*b*c*s, random="s")
expr <- strsplit(CT[CT != ""], split="")
expr <- unlist(lapply(expr, paste, collapse="*"))
expr <- sapply(expr, function(x) eval(parse(text=x)))
CT[CT != ""] <- expr
CT[CT == ""] <- 0
mode(CT) <- "numeric"
# residual variance and A*B*C*S variance are confounded in
# this design, so remove the A*B*C*S variance component
CT <- CT[-15,-2]

# get mean squares
MSmod <- summary(aov(y ~ a*b*c*subject, data=d))
MS <- MSmod[[1]][,"Mean Sq"]

# solve
ans <- solve(CT, MS)
cbind(rev(ans[c(grep("e",names(ans)),grep("s",names(ans)))])/
        c(1,2,2,2,4,4,4,1))
# s      0.04284033
# a:s    0.03381648
# b:s   -0.04004005
# c:s    0.04184887
# a:b:s -0.03657940
# a:c:s -0.02337501
# b:c:s -0.03514457
# e      0.88224787
summary(mod6)
# Random effects:
#  Groups    Name        Variance  Std.Dev. 
#  subject   (Intercept) 7.078e-02 2.660e-01
#  subject.1 A           6.176e-02 2.485e-01
#  subject.2 B           0.000e+00 0.000e+00
#  subject.3 C           6.979e-02 2.642e-01
#  subject.4 A:B         1.549e-16 1.245e-08
#  subject.5 A:C         4.566e-03 6.757e-02
#  subject.6 B:C         0.000e+00 0.000e+00
#  Residual              6.587e-01 8.116e-01
# Number of obs: 96, groups:  subject, 12

Agora podemos ver o que é patológico sobre o exemplo original. O modelo de melhor ajuste é aquele que implica que vários componentes de variação aleatória são negativos. Mas lmer()(e a maioria dos outros programas de modelos mistos) restringe as estimativas dos componentes de variação a não serem negativas. Isso geralmente é considerado uma restrição sensata, uma vez que as variações nunca podem, de fato, ser negativas. No entanto, uma conseqüência dessa restrição é que modelos mistos são incapazes de representar com precisão conjuntos de dados que apresentam correlações intraclasses negativas, ou seja, conjuntos de dados em que as observações do mesmo cluster são menos(em vez de mais) semelhante em média do que as observações extraídas aleatoriamente do conjunto de dados e, consequentemente, onde a variação dentro do cluster excede substancialmente a variação entre os conjuntos. Esses conjuntos de dados são perfeitamente razoáveis ​​que, ocasionalmente, são encontrados no mundo real (ou simulam acidentalmente!), Mas não podem ser descritos de forma sensata por um modelo de variação-componentes, porque implicam componentes de variação negativos. No entanto, eles podem ser "não sensatos" descritos por esses modelos, se o software permitir. aov()permite. lmer()não.

Jake Westfall
fonte
+1. Re I am not sure why these two methods give different results, as again I think they are "in principle" equivalent, but maybe it is for some numerical/computational reasons- você talvez entender isso melhor agora (dois anos depois)? Tentei descobrir qual é a diferença, mas também não a entendi ...
ameba diz Reinstate Monica
@amoeba Meu pensamento atual ainda é praticamente o mesmo de então: AFAIK, os dois modelos são estatisticamente equivalentes (no sentido de que eles fazem as mesmas previsões sobre os dados e implicam os mesmos erros padrão), mesmo que os efeitos aleatórios sejam parametrizados diferentemente. Penso que as diferenças observadas - que parecem ocorrer apenas algumas vezes - se devem apenas a problemas computacionais. Em particular, suspeito que você possa mexer nas configurações do otimizador (como variar pontos de partida ou usar critérios de convergência mais rigorosos) até que os dois modelos retornem exatamente a mesma resposta.
Jake Westfall
Obrigado pela sua resposta. Não estou convencido: tentei mexer nas configurações do otimizador e não consegui alterar nada nos resultados; minha impressão é que os dois modelos estão bem convergidos. Talvez eu faça uma pergunta separada em algum momento.
Ameba diz Reinstate Monica
Ak(1|A:sub)(0+A|sub)k−1k(k−1)/2k=2os dois métodos estimam um parâmetro, e ainda não tenho muita certeza do porquê de discordarem.
Ameba diz Reinstate Monica
Voltando a esse problema ... notei que, no caso de dois fatores em que duas lmerchamadas produzem anova()saída idêntica , as variações de efeito aleatório são, no entanto, bem diferentes: veja VarCorr(mod1)e VarCorr(mod2). Não entendo bem por que isso acontece; você? Para mod3e mod4, pode-se ver que quatro de sete variações para mod3são realmente iguais a zero (para mod4todas as sete são diferentes de zero); essa "singularidade" mod3é provavelmente o motivo pelo qual as tabelas anova diferem. Além disso, como você usaria o seu "caminho preferencial" se ae bteve mais de dois níveis?
Ameba diz Reinstate Monica
1

São a, b, cfixos ou efeitos aleatórios? Se eles forem corrigidos, sua sintaxe será simplesmente

summary(aov(y~a*b*c+Error(subject), d))
Error: subject
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 11  13.47   1.225               

Error: Within
          Df Sum Sq Mean Sq F value  Pr(>F)   
a          1   1.41   1.408   1.730 0.19235   
b          1   3.12   3.118   3.829 0.05399 . 
c          1   6.38   6.381   7.836 0.00647 **
a:b        1   1.57   1.571   1.929 0.16889   
a:c        1   1.09   1.091   1.339 0.25072   
b:c        1   1.41   1.413   1.735 0.19168   
a:b:c      1   0.10   0.101   0.124 0.72518   
Residuals 77  62.70   0.814  

library(lmerTest)
anova(lmer(y ~ a*b*c+(1|subject), data=d))
Analysis of Variance Table of type 3  with  Satterthwaite 
approximation for degrees of freedom
      Sum Sq Mean Sq NumDF  DenDF F.value   Pr(>F)   
a     1.4085  1.4085     1 76.991  1.7297 0.192349   
b     3.1180  3.1180     1 76.991  3.8291 0.053995 . 
c     6.3809  6.3809     1 76.991  7.8363 0.006469 **
a:b   1.5706  1.5706     1 76.991  1.9289 0.168888   
a:c   1.0907  1.0907     1 76.991  1.3394 0.250716   
b:c   1.4128  1.4128     1 76.991  1.7350 0.191680   
a:b:c 0.1014  0.1014     1 76.991  0.1245 0.725183
Masato Nakazawa
fonte
Eles são efeitos fixos. No entanto, o modelo ANOVA que você se encaixa não é o modelo que parece ser o modelo ANOVA de medidas repetidas clássicas, veja, por exemplo, aqui . Veja os estratos de erro no seu caso e no meu.
Henrik
1
Na verdade, como eles fazem isso é incorreto. Se você tiver um planejamento fatorial totalmente cruzado de medidas repetidas (ou planejamento fatorial de blocos randomizados), deverá obter apenas 1 termo de erro, além de subject, para todos os efeitos (por exemplo, Within). Veja Design Experimental: Procedimentos para Ciências Comportamentais (2013) por Kirk, capítulo 10 (p.458) ou meu post aqui
Masato Nakazawa
Vamos contornar esta questão por enquanto e supor que o modelo que eu instalei seria o modelo correto. Como você ajustaria isso usando lmer? No entanto, vou receber minha cópia de Kirk (apenas na 2ª edição) e ver o que ela diz.
Henrik
Estou curioso para saber o que você acha do capítulo de Kirk. Eu acho que o número do capítulo na 2ª ed. é diferente. Enquanto isso, tentarei ajustar lmermodelos diferentes . A melhor maneira de verificar o ajuste do modelo é verificar seus dfs usando, lmerTestporque a aproximação KR deve fornecer a você exactdfs e, portanto, valores de p.
Masato Nakazawa
1
Eu tenho a 2ª edição de Kirk, onde acredito que a discussão relevante está nas páginas 443-449, que discute um exemplo de mão dupla (e não de mão dupla). Os quadrados médios esperados, assumindo aditividade de A e B ou não, são dados na p. 447. Assumindo que A e B são fixos e os sujeitos / blocos são aleatórios, podemos ver nos quadrados médios esperados listados por Kirk sob "modelo não aditivo" que os testes de A, B e AB envolvem termos de erro diferentes, a saber: as interações relevantes com o bloco / assunto. O mesmo princípio se estende ao presente exemplo tríplice.
Jake Westfall