logLik.lm(): Почему R использует (p + 1) вместо p для степени свободы?

Question

Я пытаюсь понять результаты AIC/BIC в R. По какой-то причине 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* 

Answer 1

Мы были на самом деле очень близко к отве когда мы решили осмотреть stats:::logLik.lm. Если бы мы дополнительно осмотрели stats:::logLik.glm(Спасибо за пример glm от @ crayfish44: Mate, вы потрясающий. Еще раз вы дадите мне вдохновение, поскольку последнее сообщение касается 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 является эффект число коэффициентов модели после ранга-обнаружения.

• Когда мы имеем "gaussian()", "Gamma()" и "inverse.gaussian()" ответ, степень свободы добавляется 1, как мы должны оценку дисперсии параметра экспоненциального распределения.
• Для ответа «binomial()» и «poisson()» параметр дисперсии известен как 1, поэтому его не нужно оценивать.

Возможно, ?logLik следует рассмотреть объяснение этого, в случае, если некоторые из них глупы, как мы!