logLik.lm(): 为什么 R 使用 (p + 1) 而不是 p 作为自由度?

logLik.lm(): Why does R use (p + 1) instead of p for degree of freedom?

我试图理解 R 中 AIC/BIC 的结果。出于某种原因,R 将要估计的参数数量加 1。因此,R 使用与 2 * p - 2 * logLik 不同的公式(在高斯情况下 logLik 是残差平方和)。实际上它使用:2 * (p + 1) - 2 * logLik.

经过研究,我发现问题与stats:::logLik.lm()有关。

> stats:::logLik.lm  ## truncated R function body
## ...
##     attr(val, "df") <- p + 1
## ...

作为真实示例(使用 R 的内置数据集 trees),请考虑:

x <- lm(Height ~ Girth, trees)  ## a model with 2 parameters
logLik(x)
## 'log Lik.' -96.01663 (df=3)

这真是百思不得其解。有人知道为什么吗?


Edit1:glm @crayfish44 的例子

model.g <- glm(dist ~ speed, cars, family=gaussian)
logLik(model.g) # df=3
model.p <- glm(dist ~ speed, cars, family=poisson)
logLik(model.p) #df=2
model.G <- glm(dist ~ speed, cars, family=Gamma)
logLik(model.G) #df=3

Edit2:logLik

的方法
> methods(logLik)
[1] logLik.Arima*  logLik.glm*  logLik.lm*  logLik.logLik* logLik.nls*

当我们决定检查 stats:::logLik.lm 时,我们实际上非常接近答案。如果我们进一步检查 stats:::logLik.glm (感谢@crayfish44 的 glm 示例:伙计,你真棒。再一次你给了我灵感,因为上次 post 关于 persp()trans3d()。谢谢!),我们本来可以解决问题的。

使用:::的缺点是我们无法查看代码的注释。所以我决定检查一下 R-3.3.0 的源文件。您可以打开文件 R-3.3.0/src/library/stats/R/logLik.R 查看通用函数 logLik.**.

的注释代码
## log-likelihood for glm objects
logLik.glm <- function(object, ...)
{
    if(!missing(...)) warning("extra arguments discarded")
    fam <- family(object)$family
    p <- object$rank
    ## allow for estimated dispersion
    if(fam %in% c("gaussian", "Gamma", "inverse.gaussian")) p <- p + 1
    val <- p - object$aic / 2
    ## Note: zero prior weights have NA working residuals.
    attr(val, "nobs") <- sum(!is.na(object$residuals))
    attr(val, "df") <- p
    class(val) <- "logLik"
    val
 }

注意以下几行:

p <- object$rank
## allow for estimated dispersion
if(fam %in% c("gaussian", "Gamma", "inverse.gaussian")) p <- p + 1

p为rank-detection后模型系数的影响数

  • 当我们有"gaussian()""Gamma()""inverse.gaussian()"响应时,自由度加1,因为我们需要估计指数分布的离散参数。
  • 对于“binomial()”和“poisson()”响应,已知色散参数为 1,因此无需估计。

也许?logLik应该考虑解释一下,以防有像我们一样愚蠢的人!