Mengapa lrtest () tidak cocok dengan anova (test = “LRT”)

15

Saya mencari cara untuk melakukan uji rasio kemungkinan dalam R untuk membandingkan model yang cocok. Saya pertama kali mengkodekannya sendiri, kemudian menemukan anova()fungsi default dan juga lrtest()dalam lmtestpaket. Ketika saya memeriksa, anova()selalu menghasilkan nilai p yang sedikit berbeda dari dua lainnya meskipun parameter 'test' diatur ke "LRT". Apakah anova()sebenarnya melakukan beberapa tes yang agak berbeda, atau apakah saya tidak memahami sesuatu?

Platform: R 3.2.0 berjalan pada Linux Mint 17, lmtestversi 0.9-33

Kode sampel:

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)

Ketika saya menjalankannya, anova()memberikan nilai p 0,6071, sedangkan dua lainnya memberikan 0,60599. Perbedaan kecil, tetapi konsisten, dan terlalu besar untuk menjadi tidak akurat dalam bagaimana angka floating point disimpan. Adakah yang bisa menjelaskan mengapa anova()memberi jawaban yang berbeda?

r anova likelihood-ratio Jason
sumber

Jawaban:

7

Statistik pengujian diturunkan secara berbeda. anova.lmlistmenggunakan perbedaan skala jumlah residu kuadrat:

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
Roland
sumber
16

n-kn

Tes rasio kemungkinan diterapkan dalam lrtest()menggunakan estimator ML untuk setiap model secara terpisah sementara anova(..., test = "LRT")menggunakan estimator OLS di bawah alternatif.

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

Maka statistik yang lrtest()menghitung adalah

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") di sisi lain menggunakan

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

Di bawah hipotesis nol keduanya tentu saja setara secara asimptotik, tetapi dalam sampel terbatas terdapat perbedaan kecil.

Achim Zeileis
sumber
1
Terima kasih atas jawabannya. Jadi, dapatkah kita mengatakan bahwa satu varian lebih baik daripada yang lain? Bisakah saya menggunakan anova-test tanpa khawatir?
Julian
1
Saya tidak tahu hasil teoretis mengenai pertanyaan ini, tetapi saya tidak akan terkejut jika varian OLS berkinerja lebih baik dalam sampel kecil dengan kesalahan Gaussian. Tetapi dalam sampel yang cukup besar perbedaannya harus diabaikan.
Achim Zeileis