R中的似然比检验

机器算法验证 r 物流 诊断
2022-02-15 00:10:39

假设我要对几个自变量进行单变量逻辑回归,如下所示:

mod.a <- glm(x ~ a, data=z, family=binominal("logistic"))
mod.b <- glm(x ~ b, data=z, family=binominal("logistic"))

我做了一个模型比较(似然比测试),看看这个命令模型是否比空模型好

1-pchisq(mod.a$null.deviance-mod.a$deviance, mod.a$df.null-mod.a$df.residual)

然后我建立了另一个包含所有变量的模型

mod.c <- glm(x ~ a+b, data=z, family=binomial("logistic"))

为了查看变量在多元模型中是否具有统计意义,我使用了lrtest来自的命令epicalc

lrtest(mod.c,mod.a) ### see if variable b is statistically significant after adjustment of a
lrtest(mod.c,mod.b) ### see if variable a is statistically significant after adjustment of b

我想知道pchisq方法和lrtest方法是否等效于进行对数似然检验?因为我不知道如何lrtest用于统一逻辑模型。

2个回答

基本上,是的,只要您使用对数似然的正确差异:

> library(epicalc)
> model0 <- glm(case ~ induced + spontaneous, family=binomial, data=infert)
> model1 <- glm(case ~ induced, family=binomial, data=infert)
> lrtest (model0, model1)
Likelihood ratio test for MLE method 
Chi-squared 1 d.f. =  36.48675 , P value =  0 
> model1$deviance-model0$deviance
[1] 36.48675

而不是在两种情况下都相同的空模型的偏差。df的个数是两个嵌套模型之间不同的参数个数,这里df=1。lrtest()顺便说一句,您只需键入即可查看源代码

> lrtest

在 R 提示符下。

另一种选择是lmtest包,它具有lrtest()接受单个模型的功能。以下是包中的示例?lrtestlmtest它适用于 LM,但也有一些适用于 GLM 的方法:

> require(lmtest)
Loading required package: lmtest
Loading required package: zoo
> ## with data from Greene (1993):
> ## load data and compute lags
> data("USDistLag")
> usdl <- na.contiguous(cbind(USDistLag, lag(USDistLag, k = -1)))
> colnames(usdl) <- c("con", "gnp", "con1", "gnp1")
> fm1 <- lm(con ~ gnp + gnp1, data = usdl)
> fm2 <- lm(con ~ gnp + con1 + gnp1, data = usdl)
> ## various equivalent specifications of the LR test
>
> ## Compare two nested models
> lrtest(fm2, fm1)
Likelihood ratio test

Model 1: con ~ gnp + con1 + gnp1
Model 2: con ~ gnp + gnp1
  #Df  LogLik Df  Chisq Pr(>Chisq)    
1   5 -56.069                         
2   4 -65.871 -1 19.605  9.524e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
>
> ## with just one model provided, compare this model to a null one
> lrtest(fm2)
Likelihood ratio test

Model 1: con ~ gnp + con1 + gnp1
Model 2: con ~ 1
  #Df   LogLik Df  Chisq Pr(>Chisq)    
1   5  -56.069                         
2   2 -119.091 -3 126.04  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1