0

我想使用 ldply() 从 GLM 模型列表中提取偏差

示例数据（来自 R 基础安装）：

    library(reshape2)
    library(plyr)
    mtcars.1 <- mtcars[, c("am",  "qsec" , "drat")  ]
    mtcars.m <- melt(mtcars.1, id= 1 ) 

    glm.cars <- dlply( mtcars.m ,  .(variable) ,  
    glm,  formula=  am ~ value , family=binomial )

走到这一步：

    ldply(  glm.cars  ,  summarise ,   "Null Deviance" = null.deviance , 
        "Residual Deviance" = deviance , "Deviance"= "??"    )

这给出了这个：

      variable  Null Deviance     Residual Deviance    Deviance
1     qsec      43.22973          41.46512             ??
2     drat      43.22973          21.65003             ??

偏差不见了！我该如何提取它？

那么如何提取上面示例中的偏差呢？

当然，我可以做 null.deviance + deviance ，但我只是不想那样做。我想我想更好地了解 G 统计数据的原因。我觉得我经历了提取、减去和做 chisqr 的步骤，我会学得更好。

PS 我很困惑地发现 glm.model$devinc

4

1 回答 1

3

正如你所说，你很困惑。对于每个模型，您都有两个偏差。有趣的统计测量是这两个偏差的差异（......不是它们的总和）。（我猜你是在类比残差平方和和模型平方和的加性性质，但如果是这样，那么你在语言类比洞中跟着错误的兔子走。）你需要将差异与95% 卡方值，其自由度与空模型和“残差模型”之间的自由度差异相同。如果您在模型上执行 str(.)，您可以向下滚动列表输出以发现名称是：

 str(glm(am~qsec, mtcars, family=binomial)  )
 .....
 $ deviance         : num 41.5
 $ aic              : num 45.5
 $ null.deviance    : num 43.2
 .....
 $ df.residual      : int 30
 $ df.null          : int 31
 .....

因此，您的 dlply 代码需要提取这些，然后您计算null.deviance-deviance并df.null -df.residual可能显示qchisq(0.95, df.null-df.residual). 如果您想了解它是如何被 R-Core 打包的，请查看：

 anova( glm(am~qsec, mtcars, family=binomial)  )
于 2013-01-25T16:19:06.800 回答