# 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:

1. 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}$$

2. 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.

3. 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.

• The first thing I am noticing is that you are using the observed Fisher information matrix instead of the Fisher information matrix. First you take expectations of the second derivatives of the log-likelihood, then you evaluate them at the parameter values for the individual tests. Secondly, for the score test you are using the formula for model's with one parameter when this model has 3 parameters. Jan 11, 2022 at 18:19

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}$$.