AIC, erro anova: os modelos não são todos ajustados ao mesmo número de observações, os modelos não foram todos ajustados ao mesmo tamanho do conjunto de dados

9

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:

  1. Adicionando method = "ML"à chamada lme ​​() - não tenho certeza se é uma boa ideia alterar o método.
  2. 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 a lmer()não mostra valores de p para coeficientes - eu gosto de usar mais velhos lme().

Você tem alguma idéia se isso é um bug ou não e como podemos contornar isso?

Curioso
fonte

Respostas:

6

Uma pesquisa rápida mostra que é possível (embora eu deva admitir que pensei que não era) e que não é um bug ... apenas outro caso em que os métodos em R estão ocultos e resultam em coisas que parecem 'inesperadas ', mas a multidão do RTFM diz: "Está na documentação". Enfim ... sua solução é fazer anovacom lmeque o primeiro argumento e os lmmodelos como o segundo (e terceiro, se você quiser) argumento (s). Se isso parece estranho, é porque é um pouco estranho. O motivo é que, quando você chama anova, o anova.lmemétodo é chamado apenas se o primeiro argumento for um lmeobjeto. Caso contrário, ele chama anova.lm(que por sua vez chama anova.lmlist; se você se interessar anova.lm, verá o porquê). Para detalhes sobre como você deseja ligaranovaNesse caso, solicite ajuda anova.lme. Você verá que é possível comparar outros modelos com lmemodelos, mas eles precisam estar em uma posição diferente do primeiro argumento. Aparentemente, também é possível usar anovamodelos ajustados usando a glsfunção sem se preocupar muito com a ordem dos argumentos do modelo. Mas não conheço os detalhes suficientes para determinar se é uma boa ideia ou não, ou o que exatamente isso implica (provavelmente parece bom, mas a sua decisão). A partir desse link, comparar lmcom lmeparece estar bem documentado e citado como um método, então eu errei nessa direção, se você fosse você.

Boa sorte.

russellpierce
fonte
1
Oh, resposta e das user11852 em relação ao AIC com adendo de Gavin stands, não há AIC.lme especial ou qualquer coisa para abordar essa questão e a coisa toda começa a escorregar além do meu salário grau
russellpierce
3

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 ( m2e, m3por 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 = 0MLREMLykkX=0method="ML"

Dito isto, vamos olhar sob o capô:

 methods(AIC)  
 getAnywhere('AIC.default')

 A single object matching AIC.default was found
 It was found in the following places
   registered S3 method for AIC from namespace stats
   namespace:stats with value

 function (object, ..., k = 2) 
 {
     ll <- if ("stats4" %in% loadedNamespaces()) 
         stats4:::logLik
     else logLik
     if (!missing(...)) {
         lls <- lapply(list(object, ...), ll)
         vals <- sapply(lls, function(el) {
             no <- attr(el, "nobs") #THIS IS THE ISSUE!
             c(as.numeric(el), attr(el, "df"), if (is.null(no)) NA_integer_ else no)
         })
         val <- data.frame(df = vals[2L, ], ll = vals[1L, ])
         nos <- na.omit(vals[3L, ])
         if (length(nos) && any(nos != nos[1L])) 
             warning("models are not all fitted to the same number of observations")
         val <- data.frame(df = val$df, AIC = -2 * val$ll + k * val$df)
             Call <- match.call()
             Call$k <- NULL
         row.names(val) <- as.character(Call[-1L])
         val
     }
     else {
         lls <- ll(object)
         -2 * as.numeric(lls) + k * attr(lls, "df")
     }     
 }

onde no seu caso você pode ver que:

  lls <- lapply(list(m2,m3), stats4::logLik)
  attr(lls[[1]], "nobs")
  #[1] 500
  attr(lls[[2]], "nobs")
  #[1] 498

e obviamente logLikestá fazendo algo (talvez?) inesperado ...? não, na verdade, se você olhar para a documentação de logLik, ?logLikverá que está explicitamente declarado:

 There may be other attributes depending on the method used: see
 the appropriate documentation.  One that is used by several
 methods is "nobs"’, the number of observations used in estimation
 (after the restrictions if REML = TRUE’)

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 lmereste tempo:

m3lmer <- lmer(y ~ x + 1|cat)

e você faz o seguinte:

lls <- lapply(list(m2,m3, m3lmer), stats4::logLik)
attr(lls[[3]], "nobs")
#[1] 500
 attr(lls[[2]], "nobs")
#[1] 498

O que parece ser uma clara discrepância entre os dois, mas realmente não é como Gavin explicou. Apenas para afirmar o óbvio:

 attr( logLik(lme(y ~ x, random = ~ 1|cat, na.action = na.omit, method="ML")),
 "nobs")
#[1] 500

Há uma boa razão pela qual isso acontece em termos de metodologia, eu acho. lmetenta entender a regressão LME para você enquanto, lmerao fazer comparações de modelos, ela volta imediatamente aos resultados de ML. Eu acho que não há bugs sobre esse assunto lmee lmerapenas 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 com AIC)

usεr11852
fonte
"você deveria estar usando o ML" - mas como você pode explicar que lmerestá 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 funciona lmeré um bug. Ou 2) a mensagem de erro é um bug , não um recurso.
Curioso
Veja post atualizado (eu tive que incluir algum código). Eu mesmo notei seu ponto válido ao escrever sua resposta original, mas originalmente optei por mantê-lo fora, para que a lógica por trás da minha resposta seja estritamente baseada em computação.
usεr11852
3
O @Tomas 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 modelo MLapenas para comparar ajustes quando desejar REMLem ajustes individuais para obter melhores estimativas dos parâmetros de variação. Veja ?lmer, execute o primeiro exemplo LME até e incluindo a anova(fm1, fm2)chamada. Veja as probabilidades de log relatadas por anova()e as relatadas anteriormente na saída impressa. A anova()está ficando ML estima para você.
Gavin Simpson
Bom ponto Gavin, eu esqueço que lmerconsegui os dois ao mesmo tempo (ele usa o PLS, portanto calcula apenas um de cada vez). Eu esqueci de verificar o que você mencionou.
usεr11852
2
@ rpierce: A AIC relatada dentroanova() é a baseada em ML. O AIC relatado apenas AIC()é aquele baseado em REML.
usεr11852