机器算法验证 - 通过卡方和对数似然比较嵌套 GLM - 吾爱随笔录

通过卡方和对数似然比较嵌套 GLM

机器算法验证 r 广义线性模型卡方检验

2022-03-12 11:06:54

我想将两个 GLM 与二项式因变量进行比较。结果是：

 m1 <- glm(symptoms ~ 1,         data=data2)
 m2 <- glm(symptoms ~ phq_index, data=data2)

模型测试给出以下结果：

 anova(m1, m2)
         no AIC    logLik   LR.stat df  Pr(>Chisq)   
 m1      1  4473.9 -2236.0                        
 m2      9  4187.3 -2084.7  302.62  8   < 2.2e-16 ***

我习惯于使用卡方值、卡方差和卡方差异检验来比较这些类型的模型。由于论文中的所有其他模型都以这种方式进行比较，并且我想将它们一起报告在一个表格中：为什么这个模型测试与我得到卡方值和差异测试的其他模型测试不同？我可以从这个测试中获得卡方值吗？

其他模型比较（例如 GLMER）的结果如下所示：

    #Df AIC     BIC     logLik  Chisq   Chi     DF diff Pr(>Chisq)
m3  13  11288   11393   -5630.9 392.16          
m4  21  11212   11382   -5584.9 92.02   300.14  8       0.001

1个回答

您正在寻找的“卡方值”是偏差（-2*（对数似然），至少达到一个对于推理而言无关紧要的附加常数。R 为您提供上面的对数似然( logLik) 和似然比统计量 ( LR.stat)：LR 统计量是对数似然差的两倍 (2*(2236.0-2084.7))。

我对您的anova结果有点困惑，因为它们似乎与我anova()在两次运行时得到的格式不匹配glm()：特别是，stats:::anova.glm（对象的anova方法glm）不打印 AIC m1...m2真的glm对象吗？（例如尝试class("m1")）

你能给我们一个可重复的例子吗？这是我得到的一个简单示例，修改自?glm：

## Dobson (1990) Page 93: Randomized Controlled Trial :
d.AD <- data.frame(counts=c(18,17,15,20,10,20,25,13,12),
                   outcome=gl(3,1,9),
                   treatment=gl(3,3))
glm1 <- glm(counts ~ outcome + treatment, family = poisson, data=d.AD)
glm0 <- update(glm1, . ~ 1)

模型比较给出了剩余偏差（您的“卡方值”）以及它们之间的差异......

anova(glm0,glm1,test="Chisq")
## Analysis of Deviance Table
## 
## Model 1: counts ~ 1
## Model 2: counts ~ outcome + treatment
##   Resid. Df Resid. Dev Df    Deviance Pr(>Chi)
## 1         8    10.5814                     
## 2         4     5.1291  4      5.4523    0.244

（如果我省略了test="Chisq"，我会得到以上所有但没有 p 值）

我从您在http://article.gmane.org/gmane.comp.lang.r.general/299377的交叉发布中看到您实际上正在使用ordinal::clm，在这种情况下，我们确实得到了看起来像您上面的输出（重要的是要准确）......结果并不相同，因为模型略有不同，并且给您的是对数似然（=-deviance/2）而不是偏差，但是偏差之间的差异（“ LR.stat") 类似。

library(ordinal)
clm1 <- clm(ordered(counts) ~
   outcome + treatment, family = poisson, data=d.AD)
clm0 <- update(clm1, . ~ 1)
anova(clm0,clm1)
##       no.par      AIC  logLik LR.stat df Pr(>Chisq)
## clm0      7    50.777 -18.389                      
## clm1     11    52.839 -15.419  5.9389  4     0.2038

其它你可能感兴趣的问题