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
应该考虑解释一下,以防有像我们一样愚蠢的人!
我试图理解 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
应该考虑解释一下,以防有像我们一样愚蠢的人!