Return to Revisions

If you had consulted the help page of the logLik function, you would see

For the gaussian, Gamma and inverse.gaussian families it assumed that the dispersion of the GLM is estimated and has been counted as a parameter in the AIC value, and for all other families it is assumed that the dispersion is known. Note that this procedure does not give the maximized likelihood for "glm" fits from the Gamma and inverse gaussian families, as the estimate of dispersion used is not the MLE.

The crux is essentially the same as this question: the null distribution of (double) log-likelihood ratio $\lambda$ itself does not follow $\chi_1^2$ distribution, it is its scaled version $\sigma^{-2}\lambda \sim \chi_1^2$, where $\sigma$ is the dispersion parameter of the error term $\varepsilon$. Also note that these two tests are strictly equivalent only when $\varepsilon$ is normally distributed, i.e., when $\varepsilon \sim N(0, \sigma^2)$. Otherwise neither the $F$-statistic would have $F$-distribution nor the scaled log-likelihood ratio statistic would have $\chi^2$ distribution.

Because you didn't reveal what $\sigma$ is in your data, it is actually infeasible to reconcile these two tests for the reason stated above. Also because the mechanism how logLik evaluates the log-likelihood (i.e., always estimate $\sigma$ by $\hat{\sigma}$. My judgment is from the verbal description quoted above as I don't want to spend too much time on digging the source code. Yourself may confirm), the likelihood ratio statistic 2 * (logLik(withCov) - logLik(withInt)) would not rigorously follow $\chi^2$ distribution, even when $\varepsilon$ are normally distributed (this is something you cannot guarantee if you started with observations Y and X directly). What you should expect is that 2 * (logLik(withCov) - logLik(withInt)) * sigma_hat^2 / sigma^2 follows $\chi^2$ distribution, where sigma_hat is the estimated residual standard error, which could be extracted by sigma() function after fitting the model.

Based on these discussions, a simulation to verify the equivalence between the $F$-test and the log-likelihood ration test should be designed as follows:

1. Simulate $\varepsilon_1, \ldots, \varepsilon_n \text{ i.i.d. } \sim N(0, \sigma^2)$ and $x_1, \ldots, x_n$.
2. Let $y_i = a + bx_i + \varepsilon_i, i = 1, 2, \ldots, n$.
3. Run lm to collect the $F$-test p-value.
4. Calculate the p-value of the likelihood ratio test by 1 - pchisq(2 * (logLik(withCov) - logLik(withInt)) * sigma_hat^2 /sigma^2, df = 1).

Here is a snippet I used to implement this simulation:

set.seed(24)
sigma <- 0.8
n <- 50
epsilon <- rnorm(n, 0, sigma)
x <- runif(n, -1, 1)
y <- 1 + epsilon

withCov <- lm(y ~ x)
withInt <- lm(y ~ 1)
sigma_hat <- sigma(withCov)

1 - pchisq(2 * (logLik(withCov) - logLik(withInt)) * sigma_hat /sigma, df = 1)
# returns 0.3300543 

summary(withCov)

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2707 -0.3257  0.1145  0.4200  1.5336 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.8465     0.1088   7.782 4.71e-10 ***
x            -0.1841     0.1843  -0.999    0.323    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7684 on 48 degrees of freedom
Multiple R-squared:  0.02036,   Adjusted R-squared:  -5.366e-05 
F-statistic: 0.9974 on 1 and 48 DF,  p-value: 0.323

The p-values are closer, but still not identical.