Eu tenho modelos como este:
require(nlme)
set.seed(123)
n <- 100
k <- 5
cat <- as.factor(rep(1:k, n))
cat_i <- 1:k # intercept per kategorie
x <- rep(1:n, each = k)
sigma <- 0.2
alpha <- 0.001
y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma)
plot(x, y)
m1 <- lm(y ~ x)
summary(m1)
m2 <- lm(y ~ cat + x)
summary(m2)
m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit)
summary(m3)
Agora, estou tentando avaliar se o efeito aleatório deve estar presente no modelo. Então, comparo os modelos usando AIC ou anova e recebo o seguinte erro:
> AIC(m1, m2, m3)
df AIC
m1 3 1771.4696
m2 7 -209.1825
m3 4 -154.0245
Warning message:
In AIC.default(m1, m2, m3) :
models are not all fitted to the same number of observations
> anova(m2, m3)
Error in anova.lmlist(object, ...) :
models were not all fitted to the same size of dataset
Como você pode ver, nos dois casos eu uso o mesmo conjunto de dados. Encontrei dois remédios, mas não os considero satisfatórios:
- Adicionando
method = "ML"
à chamada lme () - não tenho certeza se é uma boa ideia alterar o método. - Usando em
lmer()
vez disso. Surpreendentemente, isso funciona, apesar de lmer () usar o método REML. No entanto, eu não gosto desta solução, porque almer()
não mostra valores de p para coeficientes - eu gosto de usar mais velhoslme()
.
Você tem alguma idéia se isso é um bug ou não e como podemos contornar isso?
fonte
Isso é peculiar definitivamente. Como um primeiro pensamento: ao fazer a comparação de modelos em que os modelos estão tendo diferentes estruturas de efeitos fixos (ML y k kX=0
m2
e,m3
por exemplo), é melhor para nós o pois "mudará" . (Ele será multiplicado por , onde ) Isso é interessante porque funciona usando o que me faz acreditar que pode não ser um bug. Parece quase que aplica "boas práticas".y k k X = 0REML
method="ML"
Dito isto, vamos olhar sob o capô:
onde no seu caso você pode ver que:
e obviamente
logLik
está fazendo algo (talvez?) inesperado ...? não, na verdade, se você olhar para a documentação delogLik
,?logLik
verá que está explicitamente declarado:o que nos leva de volta ao nosso ponto original, você deveria estar usando
ML
.Para usar um ditado comum no CS: "Não é um bug; é um recurso (real)!"
EDIT : (Apenas para endereçar seu comentário :) Suponha que você se enquadre em outro modelo usando
lmer
este tempo:e você faz o seguinte:
O que parece ser uma clara discrepância entre os dois, mas realmente não é como Gavin explicou. Apenas para afirmar o óbvio:
Há uma boa razão pela qual isso acontece em termos de metodologia, eu acho.
lme
tenta entender a regressão LME para você enquanto,lmer
ao fazer comparações de modelos, ela volta imediatamente aos resultados de ML. Eu acho que não há bugs sobre esse assuntolme
elmer
apenas diferentes razões por trás de cada pacote.Veja também o comentário de Gavin Simposon sobre uma explicação mais perspicaz do que aconteceu
anova()
(a mesma coisa acontece praticamente comAIC
)fonte
lmer
está usando REML (consulte o resumo do modelo) e funciona bem na AIC? Portanto, existem duas possibilidades: 1) a mensagem de erro é * um recurso , não um bug, e o fato de que ele funcionalmer
é um bug. Ou 2) a mensagem de erro é um bug , não um recurso.lmer()
não usa REML quando você solicita comparações. No IIRC, eles incluíram um pouco de açúcar sofisticado,lmer()
então você não precisou reajustar o modeloML
apenas para comparar ajustes quando desejarREML
em ajustes individuais para obter melhores estimativas dos parâmetros de variação. Veja?lmer
, execute o primeiro exemplo LME até e incluindo aanova(fm1, fm2)
chamada. Veja as probabilidades de log relatadas poranova()
e as relatadas anteriormente na saída impressa. Aanova()
está ficando ML estima para você.lmer
consegui os dois ao mesmo tempo (ele usa o PLS, portanto calcula apenas um de cada vez). Eu esqueci de verificar o que você mencionou.anova()
é a baseada em ML. O AIC relatado apenasAIC()
é aquele baseado em REML.