¿Por qué lrtest() no coincide con anova(test="LRT")?

Question

Estaba buscando la manera de hacer una prueba de razón de verosimilitud en R para comparar los ajustes del modelo. Primero lo codifiqué yo mismo, luego encontré que tanto el método por defecto anova() y también lrtest() en el lmtest paquete. Sin embargo, cuando lo comprobé, anova() siempre produce un valor p ligeramente diferente de los otros dos, aunque el parámetro 'test' esté ajustado a "LRT". Es anova() ¿Realizando alguna prueba sutilmente diferente, o no estoy entendiendo algo?

Plataforma: R 3.2.0 corriendo en Linux Mint 17, lmtest versión 0.9-33

Código de ejemplo:

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)

Cuando lo ejecute, anova() da un valor p de 0,6071, mientras que los otros dos dan 0,60599. Una diferencia pequeña, pero consistente, y demasiado grande para ser imprecisión en cómo se almacenan los números en coma flotante. ¿Puede alguien explicar por qué anova() da una respuesta diferente?

Answer 1

Como se menciona en la respuesta anterior, la diferencia se reduce a una diferencia de escala, es decir, a diferentes estimadores de la desviación típica de los errores. Las fuentes de la diferencia son (1) el escalado por $n-k$ (el estimador MCO insesgado) frente al escalado por $n$ (el estimador ML sesgado), y (2) utilizar el estimador bajo la hipótesis nula o alternativa.

La prueba del cociente de probabilidades aplicada en lrtest() utiliza el estimador ML para cada modelo por separado, mientras que anova(..., test = "LRT") utiliza el estimador MCO bajo la alternativa.

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

Entonces la estadística que lrtest() calcula es

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 otra parte utiliza

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

Bajo la hipótesis nula ambas son asintóticamente equivalentes, por supuesto, pero en muestras finitas hay una pequeña diferencia.

Answer 2

La estadística de prueba se obtiene de forma diferente. anova.lmlist utiliza la diferencia escalada de la suma de cuadrados residuales:

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