Replicate t or F test from regression using regression likelihoods

Question

I've heard that the t-test and F-test we use to get the significance of our regression results are derived from the likelihood ratio test, but I'm having trouble replicating the p-value of the t/F tests with a likelihood ratio test on the regression likelihoods

Using the dataset at bottom, running these three regressions in R:

withCov<-lm(Y~X)
logLik(withCov) # 'log Lik.' -61.98043 (df=3)
withInt<-lm(Y~1)
logLik(withInt) # 'log Lik.' -63.18456 (df=2)
withNone<-lm(Y~0)
logLik(withNone) # 'log Lik.' -65.32909 (df=1)

When I do a likelihood ratio test for the significance of beta_1 for example I get

1-pchisq(2*(logLik(withCov)-logLik(withInt)),df = 1)
# 0.1206958
# compare to 
summary(withCov)
Call:
lm(formula = Y ~ X)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.98508 -0.44415 -0.02294  0.59907  1.66593 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -0.2568     0.1206  -2.129   0.0384 *
X            -0.2045     0.1329  -1.539   0.1304  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8531 on 48 degrees of freedom
Multiple R-squared:  0.04702,   Adjusted R-squared:  0.02717 
F-statistic: 2.369 on 1 and 48 DF,  p-value: 0.1304

While the values are close, 0.1304 is clearly not the same as 0.1206958. And for all 6 combinations of likelihood ratio tests I cannot recover any of the p-values from lm.

What am I doing wrong? thank you!

The below R code sets up the X and Y variables:

X<-c(0.147462983739098,-0.552822250273655,0.0791413008277721,1.64705442993914,0.68248797132517,-0.315633973636888,1.40047872456033,0.776033233883272,-0.343689753439294,-0.526351059575156,0.275611183793461,1.7751174954305,-0.932913434395753,0.163895795533468,-1.00807323166609,-0.434797856886559,0.112072003876197,-0.319445372479459,1.6989373340732,0.639738727242407,0.02980494131095,0.564230726111237,1.61191604859755,-0.474733281647485,-1.46213462226689,-1.33823497263023,0.0771907623203949,-0.698998017227839,-0.775816444324552,1.46468840793603,-0.0940659257837727,-1.85718512842644,0.109762500597346,0.293088069440979,-1.33774986808507,0.804321505460817,1.24638780351287,-1.66909878637454,-0.107871283787089,-0.286526054190293,-1.30268476505327,0.241186917275982,0.0941940655245403,0.426156461908492,-0.951908401332523,-0.782389908678191,0.436387212517629,0.491981905432585,0.863964246361868,-0.715853080383197)
Y<-c(0.123592053622774,-0.156626343170975,-1.00704515111936,-0.333485064105835,0.539121555715671,0.0827543989201612,-1.38138102448818,-0.510996824863877,1.40679987755037,-0.289587787513398,-1.00983646108717,0.20397559183913,-1.53795623798144,-0.374138120013244,0.753381985793875,-0.195106670583694,-0.395410227522805,0.314058948313302,-0.567075906102368,-0.823999734395967,0.334195737940288,-0.382748868492066,-1.1924589353617,-2.14478307786596,0.198528331051895,0.616636192164357,-1.66623239116232,0.906880778089087,0.663909698659448,-0.96026161021686,-0.0289513529973109,0.103477089026045,-0.517840007699842,-1.42770065687853,-1.18735568345121,0.441872906307648,-0.814014579874374,-0.96658546843308,0.474244931448731,-0.975319336891573,-0.0672523439350411,-0.743286641930158,0.788159757412595,1.23723123779866,-0.508581741865352,0.220065470600985,0.822051420773978,-0.383737198032438,-2.04890944754108,1.5555393999865)

Answer 1

The short answer is: you should NOT expect these two tests are strictly equivalent (or mutually derivable).

To corroborate my point above, we need to review some theory. To this end, let's denote the linear model in matrix form by $y = X\beta + \epsilon$, where $y, \epsilon \in \mathbb{R}^{n}, X \in \mathbb{R}^{n \times (p + 1)}$. We are interested in testing the hypothesis:

$$H_0: \beta_1 = \cdots = \beta_p = 0 \text{ v.s. } H_a: \text{ at least one } \beta_i \text{ is non-zero.} \tag{$\dagger$}\label{0}$$

Call the reduced model when $H_0$ is true $M_0$, and the general model $M$. Furthermore, for a given dataset $(y, X)$, denote the residual sum of squares of $M_0$ and $M$ by $RSS(M_0)$ and $RSS(M)$ respectively.

Let's consider the following two cases: $\epsilon \sim N(0, \sigma^2 I_{(n)})$ and $\epsilon$ is arbitrary error.

Case 1: $\epsilon \sim N(0, \sigma^2 I_{(n)})$.

In this case, it is well-known that the likelihood-ratio test statistic for testing $\eqref{0}$ is given by \begin{align*} \Lambda = n\log\left(\frac{RSS(M_0)}{RSS(M)}\right), \tag{1}\label{1} \end{align*} and the F-test statistic for testing $\eqref{0}$ is given by \begin{align*} F = \frac{\frac{RSS(M_0) - RSS(M)}{p - 1}}{\frac{RSS(M)}{n - p - 1}}. \tag{2}\label{2} \end{align*} Under the normality assumption, we know that $F$ in $\eqref{2}$ has an exact $F_{p - 1, n - p - 1}$ distribution (check this answer for the proof). But does $\Lambda$ in $\eqref{1}$ have an exact $\chi_{p - 1}^2$ distribution? NO! It is just asymptotically $\chi_{p - 1}^2$ distributed. This asymptotic result is a consequence of the famous Wilk's theorem. At the end of this answer, I derived the exact distribution of $\Lambda$, from which we can also see the asymptotic distribution of $\Lambda$ is $\chi_{p - 1}^2$.

In conclusion, in Case 1, the p-value of the $F$-test is exact while the $p$-value of the likelihood-ratio test is approximated. Therefore, they are not directly comparable.

Case 2: make no assumption on the distributional form of $\epsilon$.

In this case, although we can still calculate statistics $\eqref{1}$ and $\eqref{2}$, we really cannot assert what their theoretical null distributions are. The reason that $\eqref{2}$ no longer follows $F_{p - 1, n - p - 1}$ can be understood by going through the proof of why $\eqref{2}$ has $F$-distribution under Case 1 (because without the normality assumption, we cannot guarantee the numerator and the denominator have $\chi^2$ distribution, as well as their independence). The reason that $\eqref{1}$ no longer follows $\chi_{p - 1}^2$ asymptotically because with errors that are not homoscedastic Gaussian, $\eqref{1}$ in general is no longer the log-likelihood ratio statistic for comparing $M_0$ and $M$.

The numerical example you presented clearly falls in Case 2. Based on the above analysis, there is no surprise to see these two tests give different $p$-values. A subtler case is Case 1, for which I did the simulation below to illustrate:

set.seed(1)
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)

1 - pchisq(2 * (logLik(withCov) - logLik(withInt)), df = 1)
# returns 0.04640223 

summary(withCov)

Call:
lm(formula = y ~ x)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.80682 -0.38006  0.09737  0.51327  1.05222 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.08038    0.09134  11.828 7.86e-16 ***
x            0.33067    0.16610   1.991   0.0522 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6459 on 48 degrees of freedom
Multiple R-squared:  0.07627,   Adjusted R-squared:  0.05703 
F-statistic: 3.963 on 1 and 48 DF,  p-value: 0.05221

It can be seen that the two tests even gave contradictory results at the significance level $\alpha = 0.05$. Clearly, the $F$-test did a better job (as expected), and this empirical result supports my theoretical "Case 1" analysis conclusion -- even under the much stronger (and usually unrealistic) normality assumption, the $F$-test and the likelihood-ratio test still do not concur! Therefore, the first sentence of your post "the t-test and F-test we use to get the significance of our regression results are derived from the likelihood ratio test" is a misconception -- although they are numerically related (e.g., if $\Lambda$ in $\eqref{1}$ is large, then $F$ in $\eqref{2}$ is large too because $\eqref{2}$ is a monotonic function of $\eqref{1}$, see $\eqref{3}$ below), it is by no means to say that (the null distribution of) the $F$-test is "derived" from the likelihood-ratio test!

---

Exact distribution of $\Lambda$ in Case 1.

It follows by $\eqref{1}$ and $\eqref{2}$ that \begin{align*} F = \frac{n - p - 1}{p - 1}(e^{\Lambda/n} - 1), \tag{3}\label{3} \end{align*} whence the Jacobian of transformation $\eqref{3}$ is \begin{align*} J = \frac{n - p - 1}{n(p - 1)}e^{\lambda/n}. \tag{4}\label{4} \end{align*} Since the pdf of a $F_{p - 1, n - p - 1}$ r.v. is \begin{align*} f_F(x) = \frac{1}{B\left(\frac{p - 1}{2}, \frac{n - p - 1}{2}\right)} \left(\frac{p - 1}{n - p - 1}\right)^{\frac{p - 1}{2}}x^{\frac{p - 1}{2} - 1} \left(1 + \frac{p - 1}{n - p - 1}x\right)^{-\frac{n}{2} + 1}, \tag{5}\label{5} \end{align*} it follows by $\eqref{3}, \eqref{4}$ and $\eqref{5}$ that the pdf of $\Lambda$ is given by \begin{align*} f_\Lambda(\lambda) = \frac{1}{B\left(\frac{p - 1}{2}, \frac{n - p - 1}{2}\right)} e^{-\lambda/2} \left(e^{\lambda/n} - 1\right)^{\frac{p - 1}{2} - 1}\frac{1}{n}e^{\frac{2\lambda}{n}}. \tag{6}\label{6} \end{align*} Evidently, $\eqref{6}$ is different from the pdf of a $\chi_{p - 1}^2$ r.v., whose pdf is \begin{align*} g(\lambda) = \frac{1}{2^{\frac{p - 1}{2}}\Gamma\left(\frac{p - 1}{2}\right)}\lambda^{\frac{p - 1}{2} - 1}e^{-\lambda/2}. \tag{7}\label{7} \end{align*} On the other hand, as $n \to \infty$, $f_\Lambda(\lambda)$ does converge to $g(\lambda)$ point-wisely (the verification is not too hard but tedious, where we need to use $e^x \sim 1 + x$ as $x \to 0$ and $\Gamma(z) \sim \sqrt{2\pi}z^{z - 1/2}e^{-z}$ as $z \to \infty$), showing that $\Lambda \to_d \chi_{p - 1}^2$ as $n \to \infty$, as Wilk's theorem asserted.

Answer 2

I've heard that the t-test and F-test we use to get the significance of our regression results are derived from the likelihood ratio test, but I'm having trouble replicating the p-value of the t/F tests with a likelihood ratio test on the regression likelihoods.

You have been applying the likelihood ratio test based on a distribution for the likelihood ratio that is approximated based on Wilks' theorem. It is not the same as the F-test or t-test which are based on exact distributions.

The likelihood ratio test is based on the use of a (log) likelihood ratio. For OLS, this can be expressed in terms of residual sum of squares of the two models, and is directly related to the F-ratio $$\log \Lambda = -\frac{n}{2} \log \left( \frac{RSS_0}{RSS_1}\right) = \frac{n}{2} \log \left( 1+F \frac{k_1-k_0}{n-k_1}\right)$$

I am not sure whether the exact distribution for $\log\Lambda$ has been described, but using the statistic $F$ (which follows an F-distribution) instead of the $\log\Lambda$ works the same.

---

See below an R code demonstration of the equivalence between log likelihood $\log \Lambda$ and f ratio $F$

set.seed(1)
n = 10
x = rnorm(n)
y = rnorm(n)

mod1 = lm(y~x)
mod0 = lm(y~1)

logLik(mod1)-logLik(mod0) # returns 0.7652048
f = anova(mod1,mod0)[2,5]
n/2*log(1+f*(2-1)/(n-2))  # returns 0.7652048