Manually calculate Rao's score statistics (Lagrange multiplier test) I learn statistics by writing blog posts. Here in this post I just wrote, I am trying to manually calculate the fisher information, wald test, and score test for linear regression:
$$y_i = \beta_0 + \beta_1 x_i + \epsilon$$
However, I encountered problems while I am calculating Rao's score test - I get negative variance for parameters. Here is what I did briefly:

*

*the null hypothesis is given by $H_0: \beta_1 = 0$, so I found the maximum likelihood estimator under $H_0$:
$$\tilde \beta_0 = \bar y \\\\\\ 
\tilde \beta_1 = 0 \\\\\\ 
\tilde \sigma^2 = \frac{\sum_i^n (y_i - \bar y)^2}{n}
$$


*Then I calculated the negative second derivative of log-likelihood function, evaluated at $\tilde \beta_0, \tilde \beta_1, \tilde \sigma^2$ (steps are in the post), which gives me the fisher information.


*Inverse the fisher information - this matrix could be interpreted as variance/covariance for our estimated parameters, which could be used to perform inference (asymptotically).
Now I am stuck at the third step: I plugged in some values, and found the variance to be negative. I am not sure if I miscalculated any quantities, or if I misunderstand the procedure of Rao's score test. Very much appreciate it if you can check my calculation here, and point out my mistakes.

For multiparameter models, $\boldsymbol{\theta} = (\theta_1, \cdots, \theta_n)$ the score test statistic is defined as
\begin{eqnarray*}
R = \boldsymbol{s}^{\prime}(\widehat{\boldsymbol{\theta}}_0) \boldsymbol{\mbox{I}}^{-1}(\widehat{\boldsymbol{\theta}}_0)\boldsymbol{s}(\widehat{\boldsymbol{\theta}}_0),
\end{eqnarray*}
where $\boldsymbol{s}(\widehat{\boldsymbol{\theta}}_0)$ denotes the vector of partial derivatives of the log-likelihood evaluated at the MLEs under the null  hypothesis and $\boldsymbol{\mbox{I}}^{-1}(\widehat{\boldsymbol{\theta}}_0)$ denotes the Fisher information matrix (not the observed Fisher information matrix) evaluated at the MLEs under the null  hypothesis.  Note that even if we test for the exclusion of a parameter, for instance $H_0: \theta_n = 0$, we do not immediately exclude it from the log-likelihood.  $\theta_n$ is still treated as a parameter appearing in the vector of first partial derivatives of the log-likelihood, $\boldsymbol{s}(\cdot)$, as well as the Fisher information matrix $\boldsymbol{\mbox{I}}(\cdot)$.  Now to the problem at hand.
Let $\boldsymbol{j}_n \in \mathbb{R}^n$ denote the vector of ones and $\boldsymbol{I}_n$ denote the $n\times n$ identity matrix.  Additionally define $\boldsymbol{P}_n = \boldsymbol{j}_n \boldsymbol{j}_n^{\prime}/n$ and $\boldsymbol{Q}_n = \boldsymbol{I}_n - \boldsymbol{P}_n$.  Define the design matrix to be $\boldsymbol{X} = \left[\boldsymbol{j}_n \,\,\, \boldsymbol{x}\right]$, where $\boldsymbol{x} \in \mathbb{R}^n$ represents the single covariate included in the linear regression model.  We shall also define $\boldsymbol{H} = \boldsymbol{X} \left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-1}\boldsymbol{X}^{\prime}$.  It may readily be shown that $\boldsymbol{H}\boldsymbol{x} = \boldsymbol{x}$ and $\boldsymbol{H}\boldsymbol{j}_n = \boldsymbol{j}_n$.  Furthermore, let $\boldsymbol{\beta} = (\beta_0 \,\,\, \beta_1)^{\prime}$ denote the vector of regression parameters and $\sigma^2$ denote the common variance of the response variable $\boldsymbol{y} \in \mathbb{R}^n$.  For convenience, let $\boldsymbol{\theta} = (\beta_0, \beta_1, \sigma^2)$.   Note that it is quicker to define a multivariate normal distribution instead of defining the product of $n$ independent univariate normal distributions.  Specifically, the simple linear regression model may be written as
\begin{eqnarray*}
\boldsymbol{y} \sim N_n (\boldsymbol{X}\boldsymbol{\beta}, \sigma^2 \boldsymbol{I}_n).
\end{eqnarray*}
The log-likelihood of the model, up to an additive constant, is
\begin{eqnarray*}
l(\boldsymbol{\theta} | \boldsymbol{y}) =-\frac{n}{2} \log (\sigma^2) -\frac{1}{2\sigma^2} \left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^{\prime}\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right).
\end{eqnarray*}
The first partial derivatives of the log-likelihood are
\begin{eqnarray*}
\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \beta_0} &=& \frac{\boldsymbol{j}_n^{\prime} \left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)}{\sigma^2}\\
\frac{\partial l(\boldsymbol{\theta} | \boldsymbol{y})}{\partial \beta_1} &=& \frac{\boldsymbol{x}^{\prime} \left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)}{\sigma^2}\\
\frac{\partial l(\boldsymbol{\theta} | \boldsymbol{y})}{\partial \sigma^2} &=& \frac{1}{2\sigma^4} \left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^{\prime}\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right) -\frac{n}{2\sigma^2}.
\end{eqnarray*}
I have shown the component-wise partial derivatives for the regression parameters instead of expressing it more compactly via
\begin{eqnarray*}
\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \boldsymbol{\beta}} &=& \frac{\boldsymbol{X}^{\prime} \left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)}{\sigma^2}.
\end{eqnarray*}
Before moving to the derivation of the Fisher information matrix, let
\begin{eqnarray*}
\boldsymbol{s} (\boldsymbol{\theta}) = \begin{pmatrix}
\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \beta_0} & \frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \beta_1} & \frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \sigma^2}
\end{pmatrix}^{\prime}.
\end{eqnarray*}
To compute the Fisher information matrix, we can either calculate the second partial derivatives and evaluate their negative expectations, or we can obtain the variance-covariance matrix of $\boldsymbol{s} (\boldsymbol{\theta})$.  I shall do the latter.  The Fisher information matrix will be written as
\begin{eqnarray*}
\boldsymbol{\mbox{I}} (\boldsymbol{\theta}) = \begin{pmatrix}
\mbox{I}_{\beta_0 \beta_0}(\boldsymbol{\theta}) & \mbox{I}_{\beta_0\beta_1}(\boldsymbol{\theta}) & \mbox{I}_{\beta_0 \sigma^2}(\boldsymbol{\theta}) \\
\cdot & \mbox{I}_{\beta_1 \beta_1}(\boldsymbol{\theta}) & \mbox{I}_{\beta_1 \sigma^2}(\boldsymbol{\theta}) \\ \cdot & \cdot & \mbox{I}_{\sigma^2\sigma^2}(\boldsymbol{\theta})
\end{pmatrix},
\end{eqnarray*}
where the missing values ($\cdot$) are filled in by symmetry and
\begin{eqnarray*}
\mbox{I}_{\beta_0 \beta_0}(\boldsymbol{\theta}) &=& \mbox{Var} \left[\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \beta_0}\right] &=& \frac{n} {\sigma^2} \\
\mbox{I}_{\beta_0 \beta_1}(\boldsymbol{\theta}) &=& \mbox{Cov} \left[\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \beta_0},\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \beta_1}\right] &=& \frac{\boldsymbol{j}_n^{\prime} \boldsymbol{x}}{\sigma^2} \\
\mbox{I}_{\beta_0 \sigma^2}(\boldsymbol{\theta}) &=& \mbox{Cov} \left[\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \beta_0},\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \sigma^2}\right] &=& 0 \\
\mbox{I}_{\beta_1 \beta_1}(\boldsymbol{\theta}) &=& \mbox{Var} \left[\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \beta_1}\right] &=& \frac{\boldsymbol{x}^{\prime}\boldsymbol{x}} {\sigma^2} \\
\mbox{I}_{\beta_1 \sigma^2}(\boldsymbol{\theta}) &=& \mbox{Cov} \left[\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \beta_1},\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \sigma^2}\right] &=& 0 \\
\mbox{I}_{\sigma^2 \sigma^2}(\boldsymbol{\theta}) &=& \mbox{Var} \left[\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \sigma^2}\right] &=& \frac{n} {4\sigma^4}.
\end{eqnarray*}
Inverting the Fisher information matrix, we find that
\begin{eqnarray*}
\boldsymbol{\mbox{I}}^{-1} (\boldsymbol{\theta}) = \frac{\sigma^2}{\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \boldsymbol{x}} \begin{pmatrix}
\frac{\boldsymbol{x}^{\prime}\boldsymbol{x} }{n} & \frac{-\boldsymbol{j}_n^{\prime}\boldsymbol{x} }{n} & 0 \\
\cdot & 1 & 0 \\ \cdot & \cdot & \frac{2 \sigma^2}{n}
\end{pmatrix}.
\end{eqnarray*}
Now suppose we wish to test $H_0: \beta_1 = 0$ vs. $H_1: \beta_1 \ne 0$ using Rao's score test.  Since we  have obtained $\boldsymbol{s}(\boldsymbol{\theta})$ and $\boldsymbol{\mbox{I}}^{-1} (\boldsymbol{\theta})$, we need only evaluate these quantities at the MLEs under the null hypothesis and use the equation given in the first equation.
Now, given $\beta_1 = 0$, $\boldsymbol{X}\boldsymbol{\beta} = \boldsymbol{j}_n \beta_0$.  Plugging this into the first partial derivatives and setting them to zero, we find that the MLE of $\beta_0$ is $\boldsymbol{j}^{\prime}_n \boldsymbol{y}/n$ and the MLE of $\sigma^2$ is $\boldsymbol{y}^{\prime} \boldsymbol{Q}_n \boldsymbol{y}/n$ under the null hypothesis.  Hence, the vector of MLEs under the null hypothesis is given by
\begin{eqnarray*}
\boldsymbol{\widehat{\theta}}_0 = \begin{pmatrix}
\frac{\boldsymbol{j}^{\prime}_n \boldsymbol{y}}{n} & 0 & \frac{\boldsymbol{y}^{\prime} \boldsymbol{Q}_n \boldsymbol{y}}{n}
\end{pmatrix}^{\prime}.
\end{eqnarray*}
Now clearly since the MLEs of $\beta_0$ and $\sigma^2$ were the zeroes of $\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \beta_0}$ and $\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \sigma^2}$, respectively, evaluating these components of the vector of first partial derivatives at $\boldsymbol{\widehat{\theta}}_0$ will be zero.  The only nonzero component will be associated with the component $\frac{\partial l(\boldsymbol{\theta}| \boldsymbol{y})}{\partial \beta_1}$.  Specifically,
\begin{eqnarray*}
\boldsymbol{s} (\boldsymbol{\widehat{\theta}}_0) = \begin{pmatrix}
0 & \frac{n\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \boldsymbol{y}}{\boldsymbol{y}^{\prime}\boldsymbol{Q}_n \boldsymbol{y}} & 0
\end{pmatrix}^{\prime}.
\end{eqnarray*}
Since the Rao score test statistic is a quadratic form, the above form for $\boldsymbol{s} (\boldsymbol{\widehat{\theta}}_0)$ implies that we need only use the component in the second row and second column of the inverse Fisher information matrix evaluated at $\boldsymbol{\widehat{\theta}}_0$, namely, $\frac{\boldsymbol{y}^{\prime}\boldsymbol{Q}_n \boldsymbol{y}}{n\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \boldsymbol{x}}$.
Therefore, Rao's score test statistic is
\begin{eqnarray*}
R &=& \left(\frac{n\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \boldsymbol{y}}{\boldsymbol{y}^{\prime}\boldsymbol{Q}_n \boldsymbol{y}}\right)^2 \frac{\boldsymbol{y}^{\prime}\boldsymbol{Q}_n \boldsymbol{y}}{n\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \boldsymbol{x}} \\
&=& \frac{n}{\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \boldsymbol{x}} \frac{\boldsymbol{y}^{\prime} \boldsymbol{Q}_n \boldsymbol{x}\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \boldsymbol{y}}{\boldsymbol{y}^{\prime}\boldsymbol{Q}_n \boldsymbol{y}}.
\end{eqnarray*}
Now since $\boldsymbol{Q}_n \boldsymbol{x}\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \boldsymbol{Q}_n = \boldsymbol{Q}_n \boldsymbol{x}\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \ne \boldsymbol{O}$, the numerator and denominator (scaled) chi-squared terms are not independent.  This may discourage one and lead them to think that one must use the asymptotic result; namely, $R \overset{d}{\rightarrow} \chi^2_1$ as $n \rightarrow \infty$.  However, it may be shown that this test statistic can be expressed as an increasing function of the Wald test statistic, which follows an $F_{1,n-2}$ distribution.  This will be shown now.
Let the test statistic $W$ satisfy
\begin{eqnarray*}
n\left(\boldsymbol{x}^{\prime} \boldsymbol{Q}_n \boldsymbol{x}\right) W = \frac{n\boldsymbol{y}^{\prime} \boldsymbol{Q}_n \boldsymbol{x}\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \boldsymbol{y}}{\left(\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \boldsymbol{x}\right) \boldsymbol{y}^{\prime}\left[\boldsymbol{I}_n - \boldsymbol{H}\right]\boldsymbol{y}}.
\end{eqnarray*}
Now the (scaled) chi-squared terms in the  numerator and denominator are independent since
\begin{eqnarray*}
\boldsymbol{Q}_n \boldsymbol{x}\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \left[\boldsymbol{I}_n - \boldsymbol{H}\right] &=& \boldsymbol{Q}_n \boldsymbol{x}\boldsymbol{x}^{\prime}\boldsymbol{Q}_n - \boldsymbol{Q}_n \boldsymbol{x}\boldsymbol{x}^{\prime}\left[\boldsymbol{I}_n - \frac{\boldsymbol{j}_n\boldsymbol{j}_n^{\prime}}{n}\right] \boldsymbol{H} \\
&=& \boldsymbol{Q}_n \boldsymbol{x}\boldsymbol{x}^{\prime}\boldsymbol{Q}_n - \boldsymbol{Q}_n \boldsymbol{x}\boldsymbol{x}^{\prime} + \boldsymbol{Q}_n \boldsymbol{x}\boldsymbol{x}^{\prime} \boldsymbol{P}_n \\ &=& \boldsymbol{O}.
\end{eqnarray*}
Furthermore, it may be shown that
\begin{eqnarray*}
\left(\boldsymbol{x}^{\prime} \boldsymbol{Q}_n \boldsymbol{x}\right)  \boldsymbol{H} &=& \boldsymbol{x}\boldsymbol{x}^{\prime} - \frac{\boldsymbol{j}_n^{\prime}\boldsymbol{x}}{n} \left(\boldsymbol{j}_n\boldsymbol{x}^{\prime} + \boldsymbol{x}\boldsymbol{j}_n^{\prime}\right) + \boldsymbol{x}^{\prime}\boldsymbol{x} \boldsymbol{P}_n
\end{eqnarray*}
by post-multipying both sides by $\boldsymbol{x}$ and simplifying.  Next note that
\begin{eqnarray*}
\boldsymbol{Q}_n \boldsymbol{x}\boldsymbol{x}^{\prime}\boldsymbol{Q}_n &=& \boldsymbol{x}\boldsymbol{x}^{\prime} - \frac{\boldsymbol{j}_n^{\prime} \boldsymbol{x}}{n}\left(\boldsymbol{j}_n\boldsymbol{x}^{\prime} + \boldsymbol{x}\boldsymbol{j}_n^{\prime}\right) + \left(\boldsymbol{x}^{\prime} \boldsymbol{P}_n\boldsymbol{x}\right)\boldsymbol{P}_n.
\end{eqnarray*}
The above relations are useful in order to simplify the following expression
\begin{eqnarray*}
1 + \left(\boldsymbol{x}^{\prime} \boldsymbol{Q}_n \boldsymbol{x}\right) W &=& \frac{\boldsymbol{y}^{\prime}\left[\left(\boldsymbol{x}^{\prime} \boldsymbol{Q}_n \boldsymbol{x}\right) \boldsymbol{I}_n - \left(\boldsymbol{x}^{\prime} \boldsymbol{Q}_n \boldsymbol{x}\right) \boldsymbol{H} + \boldsymbol{Q}_n \boldsymbol{x}\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \right]\boldsymbol{y}}{\left(\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \boldsymbol{x}\right) \boldsymbol{y}^{\prime}\left[\boldsymbol{I}_n - \boldsymbol{H}\right]\boldsymbol{y}} \\ &=& \frac{\boldsymbol{y}^{\prime}\left[\left(\boldsymbol{x}^{\prime} \boldsymbol{Q}_n \boldsymbol{x}\right) \boldsymbol{I}_n - \left[\boldsymbol{x}^{\prime}\boldsymbol{x} - \boldsymbol{x}^{\prime}\boldsymbol{P}_n \boldsymbol{x}\right]\boldsymbol{P}_n \right]\boldsymbol{y}}{\left(\boldsymbol{x}^{\prime}\boldsymbol{Q}_n \boldsymbol{x}\right) \boldsymbol{y}^{\prime}\left[\boldsymbol{I}_n - \boldsymbol{H}\right]\boldsymbol{y}} \\
&=& \frac{\boldsymbol{y}^{\prime} \boldsymbol{Q}_n \boldsymbol{y}}{\boldsymbol{y}^{\prime}\left[\boldsymbol{I}_n - \boldsymbol{H}\right]\boldsymbol{y}}
\end{eqnarray*}
Thus we have shown that
\begin{eqnarray*}
R = \frac{n\left(\boldsymbol{x}^{\prime} \boldsymbol{Q}_n \boldsymbol{x}\right) W}{1 + \left(\boldsymbol{x}^{\prime} \boldsymbol{Q}_n \boldsymbol{x}\right) W}.
\end{eqnarray*}
Finally, let $W^{\ast} = (n-2)\left(\boldsymbol{x}^{\prime} \boldsymbol{Q}_n \boldsymbol{x}\right)W$.  Now solving for $W$ and replacing it in the equation for $R$ we have
\begin{eqnarray*}
R = \frac{n W^{\ast}}{n-2 +  W^{\ast}}.
\end{eqnarray*}
Clearly $R$ is a strictly increasing function of $W^{\ast}$ and $W^{\ast}$ may be shown to follow an $F_{1,n-2}$ distribution.  Hence, we can conduct an exact test based on the transformation to $W^{\ast}$.
