1
$\begingroup$

I want to make sure my R code below is accurate. Because generally -2*logLik(lm_object) should equal deviance(lm_object).

But in my R code below, these two numbers don't match. Do I have a bug I'm missing?

set.seed(2)
n = 300
x <- rnorm(n)
e <- rnorm(n, 0, 3)
B0 = .5 ; B1 = 2

y <- B0 + B1*x + e
m <- lm(y~x)

-2*logLik(m) #   1500.681 !!

deviance(m)  #  2612.792  !!
$\endgroup$
  • $\begingroup$ The deviance formula is sum(weighted.residuals(object)^2, na.rm = TRUE) and you can check the function getAnywhere('logLik.lm'). It should be a different one $\endgroup$ – akrun Oct 4 at 22:52
  • $\begingroup$ @akrun, thanks! But see the equality that I mentioned applies here: library(lme4);hsb <- read.csv('https://raw.githubusercontent.com/rnorouzian/e/master/hsb.csv'); f1 <- lmer(math ~ ses + (1 | sch.id), data = hsb); -2*logLik(f1); deviance(f1) $\endgroup$ – rnorouzian Oct 4 at 22:59
  • 1
    $\begingroup$ arguably the deviance calculation in lme4 is wrong. I had thought that deviance==-2L for a continuous probability density, but that's not true ... $\endgroup$ – Ben Bolker Oct 5 at 0:58
2
$\begingroup$

As this CrossValidated answer points out, the deviance (which is the difference between -2*log(L) for the model and -2*log(L) for the saturated model, multiplied by the dispersion) for the linear model is equal to the sum of squared residuals (weighted, if necessary).

Here is the computational analogue of the algebra done in the linked answer:

deviance(m) ## 2612.792
sum(residuals(m)^2) ## 2612.792
nll2_m <- -2*sum(dnorm(y, predict(m), sd = sigma(m), log=TRUE))  ## 1500.688
nll2_s <- -2*sum(dnorm(y, y, sd=sigma(m), log=TRUE))  ## 1202.688
(nll2_m - nll2_s)*sigma(m)^2  ## 2612.792
| cite | improve this answer | |
$\endgroup$
0
$\begingroup$

You could try testing with the original statistical definition of log likelihood:

#Data
set.seed(2)
n = 300
x <- rnorm(n)
e <- rnorm(n, 0, 3)
B0 = .5 ; B1 = 2
#Model
y <- B0 + B1*x + e
m <- lm(y~x)

What you did is correct:

#Code 1
-2*logLik(m) 

Output:

'log Lik.' 1500.681 (df=3)

And the statistical perspective:

#Code 2 below statistical definition
sigma <- summary(m)$sigma
-2* sum(log(dnorm(x = y, mean = predict(m), sd = sigma)))

Output:

[1] 1500.688
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ But shouldn't -2*logLik(m) equal deviance(m)? $\endgroup$ – rnorouzian Oct 4 at 22:46
  • $\begingroup$ It depends of you compare the full saturated model and the estimated model. $\endgroup$ – Duck Oct 4 at 23:00
  • $\begingroup$ But Duck, it applies in here: library(lme4);hsb <- read.csv('https://raw.githubusercontent.com/rnorouzian/e/master/hsb.csv'); f1 <- lmer(math ~ ses + (1 | sch.id), data = hsb); -2*logLik(f1); deviance(f1) $\endgroup$ – rnorouzian Oct 4 at 23:02
  • $\begingroup$ It is a generalized linear model with different degrees of freedom. That might be the reason why values are equal. In the case of a linear model the results can differ! $\endgroup$ – Duck Oct 4 at 23:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.