Warum stimmt lrtest () nicht mit anova überein (test = “LRT”)

Ich suchte nach Möglichkeiten, einen Likelihood-Ratio-Test in R durchzuführen, um Modellanpassungen zu vergleichen. Ich habe es zuerst selbst codiert, dann sowohl die Standardfunktion anova()als auch lrtest()das lmtestPaket gefunden. Bei der Überprüfung wird jedoch anova()immer ein geringfügig anderer p-Wert als bei den anderen beiden Werten erzeugt, obwohl der Parameter 'test' auf "LRT" gesetzt ist. Führt man anova()tatsächlich einen subtil anderen Test durch oder verstehe ich etwas nicht?

Plattform: R 3.2.0 unter Linux Mint 17, lmtestVersion 0.9-33

Beispielcode:

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)

Wenn ich es ausführe, anova()ergibt sich ein p-Wert von 0,6071, während die anderen beiden 0,60599 ergeben. Ein kleiner Unterschied, aber konsistent und zu groß, um die Speicherung von Gleitkommazahlen genau zu bestimmen. Kann jemand erklären, warum anova()eine andere Antwort gibt?

Die Teststatistik wird unterschiedlich abgeleitet. anova.lmlistVerwendet die skalierte Differenz der verbleibenden Quadratsumme:

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

$n-k$$n$

Der in implementierte Likelihood-Ratio-Test lrtest()verwendet den ML-Schätzer für jedes Modell separat, während anova(..., test = "LRT")alternativ der OLS-Schätzer verwendet wird.

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

Dann ist die zu lrtest()berechnende Statistik

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") auf der anderen Seite verwendet

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

Unter der Nullhypothese sind beide natürlich asymptotisch äquivalent, aber bei endlichen Stichproben gibt es einen kleinen Unterschied.