Por que o lrtest () não corresponde ao anova (test = "LRT")

Eu estava procurando maneiras de fazer um teste de razão de verossimilhança em R para comparar os ajustes do modelo. Eu primeiro o codifiquei, depois encontrei a anova()função padrão e também lrtest()no lmtestpacote. Quando verifiquei, no entanto, anova()sempre produz um valor p ligeiramente diferente dos outros dois, embora o parâmetro 'test' esteja definido como "LRT". anova()Realmente está realizando algum teste sutilmente diferente ou não estou entendendo alguma coisa?

Plataforma: R 3.2.0 em execução no Linux Mint 17, lmtestversão 0.9-33

Código de amostra:

set.seed(1) # Reproducibility
n=1000
y = runif(n, min=-1, max=1)
a = factor(sample(1:5, size=n, replace=T))
b = runif(n)

# Make y dependent on the other two variables
y = y + b * 0.1 + ifelse(a==1, 0.25, 0)
mydata = data.frame(y,a,b)

# Models
base = lm(y ~ a, data=mydata)
full = lm(y ~ a + b, data=mydata)

# Anova
anova(base, full, test="LRT")

# lrtest
library(lmtest)
lrtest(base, full)

# Homebrew log-likelihood test
like.diff = logLik(full) - logLik(base)
df.diff = base$df.residual - full$df.residual
pchisq(as.numeric(like.diff) * 2, df=df.diff, lower.tail=F)

Quando eu executo, anova()dá um valor de p de 0,6071, enquanto os outros dois dão 0,60599. Uma pequena diferença, mas consistente e grande demais para ser imprecisa na maneira como os números de ponto flutuante são armazenados. Alguém pode explicar por que anova()dá uma resposta diferente?

As estatísticas de teste são derivadas de maneira diferente. anova.lmlistusa a diferença escalada da soma residual dos quadrados:

anova(base, full, test="LRT")
#  Res.Df    RSS Df Sum of Sq Pr(>Chi)
#1    995 330.29                      
#2    994 330.20  1   0.08786   0.6071

vals <- (sum(residuals(base)^2) - sum(residuals(full)^2))/sum(residuals(full)^2) * full$df.residual 
pchisq(vals, df.diff, lower.tail = FALSE)
#[1] 0.6070549

Como mencionado na resposta anterior, a diferença se resume a uma diferença no dimensionamento, ou seja, estimadores diferentes para o desvio padrão dos erros. As fontes para a diferença são (1) dimensionamento por$n-k$ (o estimador OLS imparcial) vs. escala por $n$ (o estimador enviesado de ML) e (2) usando o estimador sob a hipótese nula ou alternativa.

O teste da razão de verossimilhança implementado lrtest()usa o estimador ML para cada modelo separadamente, enquanto anova(..., test = "LRT")usa o estimador OLS sob a alternativa.

sd_ols <- function(object) sqrt(sum(residuals(object)^2)/df.residual(object))
sd_mle <- function(object) sqrt(mean(residuals(object)^2))

Então a estatística que lrtest()calcula é

ll <- function(object, sd) sum(dnorm(model.response(model.frame(object)),
  mean = fitted(object), sd = sd, log = TRUE))
-2 * (ll(base, sd_mle(base)) - ll(full, sd_mle(full)))
## [1] 0.266047

anova(..., test = "LRT") por outro lado, usa

-2 * (ll(base, sd_ols(full)) - ll(full, sd_ols(full)))
## [1] 0.2644859

Sob a hipótese nula, ambos são assintoticamente equivalentes, é claro, mas em amostras finitas há uma pequena diferença.