2  Parametric statistical inference

2.1 Introduction

Probability distributions like the binomial, Poisson and normal, enable us to calculate probabilities, and other quantities of interest (e.g. expectations) for a probability model of a random process. Therefore, given the model, we can make statements about possible outcomes of the process.

Statistical inference is concerned with the inverse problem. Given outcomes of a random process (observed data), what conclusions (inferences) can we draw about the process itself?

We assume that the $n$ observations of the response $\mathit{y}=(y_1,\dots ,y_n{)}^T$ are observations of random variables $\mathit{Y}=(Y_1,\dots ,Y_n{)}^T$, which have joint p.d.f. $f_{\mathit{Y}}$ (joint p.f. for discrete variables). We use the observed data $\mathit{y}$ to make inferences about $f_{\mathit{Y}}$.

We usually make certain assumptions about $f_{\mathit{Y}}$. In particular, we often assume that $y_1,\dots ,y_n$ are observations of independent random variables. Hence $$f_{\mathit{Y}}(\mathit{y})=f_{Y_1}(y_1)f_{Y_2}(y_2)\cdots f_{Y_n}(y_n)=\prod _{i=1}^n{f}_{Y_i}(y_i).$$

In parametric statistical inference, we specify a joint distribution $f_{\mathit{Y}}$, for $\mathit{Y}$, which is known, except for the values of parameters ${\theta }_1,{\theta }_2,\dots ,{\theta }_p$ (sometimes denoted by $\mathit{\theta }$). Then we use the observed data $\mathit{y}$ to make inferences about ${\theta }_1,{\theta }_2,\dots ,{\theta }_p$. In this case, we usually write $f_{\mathit{Y}}$ as $f_{\mathit{Y}}(\mathit{y};\mathit{\theta })$, to make explicit the dependence on the unknown $\mathit{\theta }$.

2.2 The likelihood function

We often think of the joint density $f_{\mathit{Y}}(\mathit{y};\mathit{\theta })$ as a function of $\mathit{y}$ for fixed $\mathit{\theta }$, which describes the relative probabilities of different possible values of $\mathit{y}$, given a particular set of parameters $\mathit{\theta }$. However, in statistical inference, we have observed $y_1,\dots ,y_n$ (values of $Y_1,\dots ,Y_n$). Knowledge of the probability of alternative possible realisations of $\mathit{Y}$ is largely irrelevant. What we want to know about is $\mathit{\theta }$.

Our only link between the observed data $y_1,\dots ,y_n$ and $\mathit{\theta }$ is through the function $f_{\mathit{Y}}(\mathit{y};\mathit{\theta })$. Therefore, it seems sensible that parametric statistical inference should be based on this function. We can think of $f_{\mathit{Y}}(\mathit{y};\mathit{\theta })$ as a function of $\mathit{\theta }$ for fixed $\mathit{y}$, which describes the relative likelihoods of different possible (sets of) $\mathit{\theta },$ given observed data $y_1,\dots ,y_n$. We write $$L(\mathit{\theta };\mathit{y})=f_{\mathit{Y}}(\mathit{y};\mathit{\theta })$$ for this likelihood, which is a function of the unknown parameter $\mathit{\theta }$. For convenience, we often drop $\mathit{y}$ from the notation, and write $L(\mathit{\theta })$.

The likelihood function is of central importance in parametric statistical inference. It provides a means for comparing different possible values of $\mathit{\theta }$, based on the probabilities (or probability densities) that they assign to the observed data $y_1,\dots ,y_n$.

Notes

1. Frequently it is more convenient to consider the log-likelihood function $\ell (\mathit{\theta })=\log L(\mathit{\theta })$.
2. Nothing in the definition of the likelihood requires $y_1,\dots ,y_n$ to be observations of independent random variables, although we shall frequently make this assumption.
3. Any factors which depend on $y_1,\dots ,y_n$ alone (and not on $\mathit{\theta }$) can be ignored when writing down the likelihood. Such factors give no information about the relative likelihoods of different possible values of $\mathit{\theta }$.

Example 2.1 (Bernoulli) $y_1,\dots ,y_n$ are observations of $Y_1,\dots ,Y_n$, independent identically distributed (i.i.d.) Bernoulli$(p)$ random variables. Here $\theta =(p)$ and the likelihood is $$L(p)=\prod _{i=1}^n{p}^{y_i}(1-p{)}^{1-y_i}=p^{\sum _{i=1}^n{y}_i}(1-p{)}^{n-\sum _{i=1}^n{y}_i}.$$ The log-likelihood is $$\ell (p)=\log L(p)=n\overline{y}\log p+n(1-\overline{y})\log (1-p).$$

Example 2.2 (Normal) $y_1,\dots ,y_n$ are observations of $Y_1,\dots ,Y_n$, i.i.d. $N(\mu ,{\sigma }^2)$ random variables. Here $\mathit{\theta }=(\mu ,{\sigma }^2)$ and the likelihood is $$\begin{array}{rl}L(\mu ,{\sigma }^2)& =\prod _{i=1}^n\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp (-\frac{1}{2{\sigma }^2}(y_i-\mu {)}^2)\\ & =(2\pi {\sigma }^2{)}^{-\frac{n}{2}}\exp (-\frac{1}{2{\sigma }^2}\sum (y_i-\mu {)}^2)\\ & \propto ({\sigma }^2{)}^{-\frac{n}{2}}\exp (-\frac{1}{2{\sigma }^2}\sum (y_i-\mu {)}^2).\end{array}$$ The log-likelihood is $$\ell (\mu ,{\sigma }^2)=\log L(\mu ,{\sigma }^2)=-\frac{n}{2}\log (2\pi )-\frac{n}{2}\log ({\sigma }^2)-\frac{1}{2{\sigma }^2}\sum (y_i-\mu {)}^2.$$

2.3 Maximum likelihood estimation

One of the primary tasks of parametric statistical inference is estimation of the unknown parameters ${\theta }_1,\dots ,{\theta }_p$. Consider the value of $\mathit{\theta }$ which maximises the likelihood function. This is the ‘most likely’ value of $\mathit{\theta }$, the one which makes the observed data ‘most probable’. When we are searching for an estimate of $\mathit{\theta }$, this would seem to be a good candidate.

We call the value of $\mathit{\theta }$ which maximises the likelihood $L(\theta )$ the maximum likelihood estimate (MLE) of $\mathit{\theta }$, denoted by $\hat{\mathit{\theta }}$. $\hat{\mathit{\theta }}$ depends on $\mathit{y}$, as different observed data samples lead to different likelihood functions. The corresponding function of $\mathit{Y}$ is called the maximum likelihood estimator and is also denoted by $\hat{\mathit{\theta }}$.

Note that as $\mathit{\theta }=({\theta }_1,\dots ,{\theta }_p)$, the MLE for any component of $\mathit{\theta }$ is given by the corresponding component of $\hat{\mathit{\theta }}=({\hat{\theta }}_1,\dots ,{\hat{\theta }}_p{)}^T$. Similarly, the MLE for any function of parameters $g(\mathit{\theta })$ is given by $g(\hat{\mathit{\theta }})$.

As $\log$ is a strictly increasing function, the value of $\mathit{\theta }$ which maximises $L(\mathit{\theta })$ also maximises $\ell (\mathit{\theta })=\log L(\mathit{\theta })$. It is almost always easier to maximise $\ell (\mathit{\theta })$. This is achieved in the usual way; finding a stationary point by differentiating $\ell (\mathit{\theta })$ with respect to ${\theta }_1,\dots ,{\theta }_p$, and solving the resulting $p$ simultaneous equations. It should also be checked that the stationary point is a maximum.

Example 2.3 (Bernoulli) $y_1,\dots ,y_n$ are observations of $Y_1,\dots ,Y_n$, i.i.d. Bernoulli$(p)$ random variables. Here $\mathit{\theta }=(p)$ and the log-likelihood is $$\ell (p)=n\overline{y}\log p+n(1-\overline{y})\log (1-p).$$ Differentiating with respect to $p$, $$\frac{\partial }{\partial p}\ell (p)=\frac{n\overline{y}}{p}-\frac{n(1-\overline{y})}{1-p}$$ so the MLE $\hat{p}$ solves $$\frac{n\overline{y}}{\hat{p}}-\frac{n(1-\overline{y})}{1-\hat{p}}=0.$$ Solving this for $\hat{p}$ gives $\hat{p}=\overline{y}$. Note that $$\frac{{\partial }^2}{\partial p^2}\ell (p)=-n\overline{y}/p^2-n(1-\overline{y})/(1-p{)}^2<0$$ everywhere, so the stationary point is clearly a maximum.

Example 2.4 (Normal) $y_1,\dots ,y_n$ are observations of $Y_1,\dots ,Y_n$, i.i.d. $N(\mu ,{\sigma }^2)$ random variables. Here $\mathit{\theta }=(\mu ,{\sigma }^2)$ and and the log-likelihood is $$\ell (\mu ,{\sigma }^2)=-\frac{n}{2}\log (2\pi )-\frac{n}{2}\log ({\sigma }^2)-\frac{1}{2{\sigma }^2}\sum (y_i-\mu {)}^2.$$ Differentiating with respect to $\mu$ $$\frac{\partial }{\partial \mu }\ell (\mu ,{\sigma }^2)=\frac{1}{{\sigma }^2}\sum (y_i-\mu )=\frac{n(\overline{y}-\mu )}{{\sigma }^2}$$ so $(\hat{\mu },{\hat{\sigma }}^2)$ solve

$$\begin{array}{}\text{(2.1)}& \frac{n(\overline{y}-\hat{\mu })}{{\hat{\sigma }}^2}=0.\end{array}$$

Differentiating with respect to ${\sigma }^2$ $$\frac{\partial }{\partial {\sigma }^2}\ell (\mu ,{\sigma }^2)=-\frac{n}{2{\sigma }^2}+\frac{1}{2({\sigma }^2{)}^2}\sum (y_i-\mu {)}^2,$$ so

$$\begin{array}{}\text{(2.2)}& -\frac{n}{2{\hat{\sigma }}^2}+\frac{1}{2({\hat{\sigma }}^2{)}^2}\sum (y_i-\hat{\mu }{)}^2=0\end{array}$$

Solving (Equation 2.1) and (Equation 2.2), we obtain $\hat{\mu }=\overline{y}$ and $${\hat{\sigma }}^2=\frac{1}{n}\sum (y_i-\hat{\mu }{)}^2=\frac{1}{n}\sum (y_i-\overline{y}{)}^2.$$

Strictly, to show that this stationary point is a maximum, we need to show that the Hessian matrix (the matrix of second derivatives with elements $[\mathit{H}(\mathit{\theta }){]}_{ij}=\frac{{\partial }^2}{\partial {\theta }_i\partial {\theta }_j}\ell (\theta )$) is negative definite at $\mathit{\theta }=\hat{\mathit{\theta }}$, that is ${\mathit{a}}^T\mathit{H}(\hat{\mathit{\theta }})\mathit{a}<0$ for every $\mathit{a}\ne \mathbf{0}$. Here $$\mathit{H}(\hat{\mu },{\hat{\sigma }}^2)=\left(\begin{array}{cc}-\frac{n}{{\hat{\sigma }}^2}& 0\\ 0& -\frac{n}{2({\hat{\sigma }}^2{)}^2}\end{array}\right)$$ which is clearly negative definite.

2.4 Score

Let $$u_i(\mathit{\theta })\equiv \frac{\partial }{\partial {\theta }_i}\ell (\theta ),i=1,\dots ,p$$ and $\mathit{u}(\mathit{\theta })\equiv [u_1(\mathit{\theta }),\dots ,u_p(\mathit{\theta }){]}^T$. Then we call $\mathit{u}(\mathit{\theta })$ the vector of scores or score vector. Where $p=1$ and $\mathit{\theta }=(\theta )$, the score is the scalar defined as $$u(\theta )\equiv \frac{\partial }{\partial \theta }\ell (\theta ).$$ The maximum likelihood estimate $\hat{\mathit{\theta }}$ satisfies $$u(\hat{\mathit{\theta }})=\mathbf{0},$$ that is, $$u_i(\hat{\mathit{\theta }})=0,i=1,\dots ,p.$$ Note that $u(\mathit{\theta })$ is a function of $\mathit{\theta }$ for fixed (observed) $\mathit{y}$. However, if we replace $y_1,\dots ,y_n$ in $u(\mathit{\theta })$, by the corresponding random variables $Y_1,\dots ,Y_n$ then we obtain a vector of random variables $U(\mathit{\theta })\equiv [U_1(\mathit{\theta }),\dots ,U_p(\mathit{\theta }){]}^T$.

An important result in likelihood theory is that the expected score at the true (but unknown) value of $\mathit{\theta }$ is zero:

Theorem 2.1 We have $E[U(\mathit{\theta })]=\mathbf{0}$, i.e. $E[U_i(\mathit{\theta })]=0,$ $i=1,\dots ,p,$ provided that

1. The expectation exists.
2. The sample space for $\mathit{Y}$ does not depend on $\mathit{\theta }$.

Proof. Our proof is for continuous $\mathit{y}$ – in the discrete case, replace $\int$ by $\sum$. For each $i=1,\dots ,n$ $$\begin{array}{rl}E[U_i(\mathit{\theta })]& =\int U_i(\mathit{\theta })f_{\mathit{Y}}(\mathit{y},\mathit{\theta })d\mathit{y}\\ & =\int \frac{\partial }{\partial {\theta }_i}\ell (\theta )f_{\mathit{Y}}(\mathit{y};\mathit{\theta })d\mathit{y}\\ & =\int \frac{\partial }{\partial {\theta }_i}\log f_{\mathit{Y}}(\mathit{y};\mathit{\theta })f_{\mathit{Y}}(\mathit{y};\mathit{\theta })d\mathit{y}\\ & =\int \frac{\frac{\partial }{\partial {\theta }_i}{f}_{\mathit{Y}}(\mathit{y};\mathit{\theta })}{f_{\mathit{Y}}(\mathit{y};\mathit{\theta })}{f}_{\mathit{Y}}(\mathit{y};\mathit{\theta })d\mathit{y}\\ & =\int \frac{\partial }{\partial {\theta }_i}{f}_{\mathit{Y}}(\mathit{y};\mathit{\theta })d\mathit{y}\\ & =\frac{\partial }{\partial {\theta }_i}\int f_{\mathit{Y}}(\mathit{y};\mathit{\theta })d\mathit{y}\\ & =\frac{\partial }{\partial {\theta }_i}1=0,\end{array}$$ as required.

Example 2.5 (Bernoulli) $y_1,\dots ,y_n$ are observations of $Y_1,\dots ,Y_n$, i.i.d. Bernoulli$(p)$ random variables. Here $\mathit{\theta }=(p)$ and $$u(p)=n\overline{y}/p-n(1-\overline{y})/(1-p).$$ Since $E[U(p)]=0$, we must have $E[\overline{Y}]=p$ (which we already know is correct).

Example 2.6 (Normal) $y_1,\dots ,y_n$ are observations of $Y_1,\dots ,Y_n$, i.i.d. $N(\mu ,{\sigma }^2)$ random variables. Here $\mathit{\theta }=(\mu ,{\sigma }^2)$ and $$\begin{array}{rl}{u}_1(\mu ,{\sigma }^2)& =n(\overline{y}-\mu )/{\sigma }^2\\ u_2(\mu ,{\sigma }^2)& =-\frac{n}{2{\sigma }^2}+\frac{1}{2({\sigma }^2{)}^2}\sum _{i=1}^n(y_i-\mu {)}^2\end{array}$$ Since $E[\mathit{U}(\mu ,{\sigma }^2)]=\mathbf{0}$, we must have $E[\overline{Y}]=\mu$ and $E[\frac{1}{n}\sum _{i=1}^n(Y_i-\mu {)}^2]={\sigma }^2.$

2.5 Information

Suppose that $y_1,\dots ,y_n$ are observations of $Y_1,\dots ,Y_n$, whose joint p.d.f. $L(\theta )$ is completely specified except for the values of $p$ unknown parameters $\mathit{\theta }=({\theta }_1,\dots ,{\theta }_p{)}^T$. Previously, we defined the Hessian matrix $H(\mathit{\theta })$ to be the matrix with components $$[H(\mathit{\theta }){]}_{ij}\equiv \frac{{\partial }^2}{\partial {\theta }_i\partial {\theta }_j}\ell (\theta )i=1,\dots ,p;j=1,\dots ,p.$$ We call the matrix $-H(\mathit{\theta })$ the observed information matrix. Where $p=1$ and $\mathit{\theta }=(\theta )$, the observed information is a scalar defined as $$-H(\theta )\equiv -\frac{\partial }{\partial {\theta }^2}\ell (\theta ).$$

As with the score, if we replace $y_1,\dots ,y_n$ in $H(\mathit{\theta })$, by the corresponding random variables $Y_1,\dots ,Y_n$, we obtain a matrix of random variables. Then, we define the expected information matrix or Fisher information matrix $$[\mathcal{I}(\mathit{\theta }){]}_{ij}=E(-[H(\mathit{\theta }){]}_{ij})i=1,\dots ,p;j=1,\dots ,p.$$

An important result in likelihood theory is that the variance-covariance matrix of the score vector is equal to the expected information matrix:

Theorem 2.2 We have $\text{Var}[U(\mathit{\theta })]=\mathcal{I}(\mathit{\theta })$, i.e. $$\text{Var}[U(\mathit{\theta }){]}_{ij}=[\mathcal{I}(\mathit{\theta }){]}_{ij},i=1,\dots ,p,j=1,\dots ,p$$ provided that

1. The variance exists.
2. The sample space for $\mathit{Y}$ does not depend on $\mathit{\theta }$.

Proof. Our proof is for continuous $\mathit{y}$ – in the discrete case, replace $\int$ by $\sum$.

For each $i=1,\dots ,p$ and $j=1,\dots ,p$, $$\begin{array}{rl}\text{Var}[U(\mathit{\theta }){]}_{ij}& =E[U_i(\mathit{\theta })U_j(\mathit{\theta })]\\ & =\int \frac{\partial }{\partial {\theta }_i}\ell (\theta )\frac{\partial }{\partial {\theta }_j}\ell (\theta )f_{\mathit{Y}}(\mathit{y};\mathit{\theta })d\mathit{y}\\ & =\int \frac{\partial }{\partial {\theta }_i}\log f_{\mathit{Y}}(\mathit{y};\mathit{\theta })\frac{\partial }{\partial {\theta }_j}\log f_{\mathit{Y}}(\mathit{y};\mathit{\theta })f_{\mathit{Y}}(\mathit{y};\mathit{\theta })d\mathit{y}\\ & =\int \frac{\frac{\partial }{\partial {\theta }_i}{f}_{\mathit{Y}}(\mathit{y};\mathit{\theta })}{f_{\mathit{Y}}(\mathit{y};\mathit{\theta })}\frac{\frac{\partial }{\partial {\theta }_j}{f}_{\mathit{Y}}(\mathit{y};\mathit{\theta })}{f_{\mathit{Y}}(\mathit{y};\mathit{\theta })}{f}_{\mathit{Y}}(\mathit{y};\mathit{\theta })d\mathit{y}\\ & =\int \frac{1}{f_{\mathit{Y}}(\mathit{y};\mathit{\theta })}\frac{\partial }{\partial {\theta }_i}{f}_{\mathit{Y}}(\mathit{y};\mathit{\theta })\frac{\partial }{\partial {\theta }_j}{f}_{\mathit{Y}}(\mathit{y};\mathit{\theta })d\mathit{y}.\end{array}$$

Now $$\begin{array}{rl}[\mathcal{I}(\mathit{\theta }){]}_{ij}& =E[-\frac{{\partial }^2}{\partial {\theta }_i\partial {\theta }_j}\ell (\theta )]\\ & =\int -\frac{{\partial }^2}{\partial {\theta }_i\partial {\theta }_j}\log f_{\mathit{Y}}(\mathit{y};\mathit{\theta })f_{\mathit{Y}}(\mathit{y};\mathit{\theta })d\mathit{y}\\ & =\int -\frac{\partial }{\partial {\theta }_i}\left[\frac{\frac{\partial }{\partial {\theta }_j}{f}_{\mathit{Y}}(\mathit{y};\mathit{\theta })}{f_{\mathit{Y}}(\mathit{y};\mathit{\theta })}\right]f_{\mathit{Y}}(\mathit{y};\mathit{\theta })d\mathit{y}\\ & =\int [-\frac{\frac{{\partial }^2}{\partial {\theta }_i\partial {\theta }_j}{f}_{\mathit{Y}}(\mathit{y};\mathit{\theta })}{f_{\mathit{Y}}(\mathit{y};\mathit{\theta })}+\frac{\frac{\partial }{\partial {\theta }_i}{f}_{\mathit{Y}}(\mathit{y};\mathit{\theta })\frac{\partial }{\partial {\theta }_j}{f}_{\mathit{Y}}(\mathit{y};\mathit{\theta })}{f_{\mathit{Y}}(\mathit{y};\mathit{\theta }{)}^2}]f_{\mathit{Y}}(\mathit{y};\mathit{\theta })d\mathit{y}\\ & =-\frac{{\partial }^2}{\partial {\theta }_i\partial {\theta }_j}\int f_{\mathit{Y}}(\mathit{y};\mathit{\theta })d\mathit{y}+\int \frac{1}{f_{\mathit{Y}}(\mathit{y};\mathit{\theta })}\frac{\partial }{\partial {\theta }_i}{f}_{\mathit{Y}}(\mathit{y};\mathit{\theta })\frac{\partial }{\partial {\theta }_j}{f}_{\mathit{Y}}(\mathit{y};\mathit{\theta })d\mathit{y}\\ & =\text{Var}[U(\mathit{\theta }){]}_{ij},\end{array}$$ as required.

Example 2.7 (Bernoulli) $y_1,\dots ,y_n$ are observations of $Y_1,\dots ,Y_n$, i.i.d. Bernoulli$(p)$ random variables. Here $\mathit{\theta }=(p)$ and $$\begin{array}{rl}u(p)& =\frac{n\overline{y}}{p}-\frac{n(1-\overline{y})}{(1-p)}\\ -H(p)& =\frac{n\overline{y}}{p^2}+\frac{n(1-\overline{y})}{(1-p{)}^2}\\ \mathcal{I}(p)& =\frac{n}{p}+\frac{n}{(1-p)}=\frac{n}{p(1-p)}.\end{array}$$

Example 2.8 (Normal) $y_1,\dots ,y_n$ are observations of $Y_1,\dots ,Y_n$, i.i.d. $N(\mu ,{\sigma }^2)$ random variables. Here $\mathit{\theta }=(\mu ,{\sigma }^2)$ and $$\begin{array}{rl}{u}_1(\mu ,{\sigma }^2)& =\frac{n(\overline{y}-\mu )}{{\sigma }^2}\\ u_2(\mu ,{\sigma }^2)& =-\frac{n}{2{\sigma }^2}+\frac{1}{2({\sigma }^2{)}^2}\sum (y_i-\mu {)}^2.\end{array}$$ Therefore $$-\mathit{H}(\mu ,{\sigma }^2)=\left(\begin{array}{cc}\frac{n}{{\sigma }^2}& \frac{n(\overline{y}-\mu )}{({\sigma }^2{)}^2}\\ \frac{n(\overline{y}-\mu )}{({\sigma }^2{)}^2}& \frac{1}{({\sigma }^2{)}^3}\sum (y_i-\mu {)}^2-\frac{n}{2({\sigma }^2{)}^2}\end{array}\right)$$ $$\mathcal{I}(\mu ,{\sigma }^2)=\left(\begin{array}{cc}\frac{n}{{\sigma }^2}& 0\\ 0& \frac{n}{2({\sigma }^2{)}^2}\end{array}\right).$$

2.6 Asymptotic distribution of the MLE

Maximum likelihood estimation is an attractive method of estimation for a number of reasons. It is intuitively sensible and usually reasonably straightforward to carry out. Even when the simultaneous equations we obtain by differentiating the log-likelihood function are impossible to solve directly, solution by numerical methods is usually feasible.

Perhaps the most compelling reason for considering maximum likelihood estimation is the asymptotic behaviour of maximum likelihood estimators.

Suppose that $y_1,\dots ,y_n$ are observations of independent random variables $Y_1,\dots ,Y_n$, whose joint p.d.f. $f_{\mathit{Y}}(\mathit{y};\mathit{\theta })=\prod _{i=1}^n{f}_{Y_i}(y_i;\mathit{\theta })$ is completely specified except for the values of an unknown parameter vector $\mathit{\theta }$, and that $\hat{\mathit{\theta }}$ is the maximum likelihood estimator of $\mathit{\theta }$.

Then, as $n\to \mathrm{\infty }$, the distribution of $\hat{\mathit{\theta }}$ tends to a multivariate normal distribution with mean vector $\mathit{\theta }$ and variance covariance matrix $\mathcal{I}(\mathit{\theta }{)}^{-1}$.

Where $p=1$ and $\mathit{\theta }=(\theta )$, the distribution of the MLE $\hat{\theta }$ tends to $N[\theta ,1/\mathcal{I}(\theta )]$.

For ‘large enough $n$’, we can treat the asymptotic distribution of the MLE as an approximation. The fact that $E(\hat{\mathit{\theta }})\approx \mathit{\theta }$ means that the maximum likelihood estimator is approximately unbiased for large samples. The variance of $\hat{\mathit{\theta }}$ is approximately $\mathcal{I}(\mathit{\theta }{)}^{-1}$. It is possible to show that this is the smallest possible variance of any unbiased estimator of $\mathit{\theta }$ (this result is called the Cramér–Rao lower bound, which we do not prove here). Therefore the MLE is the ‘best possible’ estimator in large samples (and therefore we hope also reasonable in small samples, though we should investigate this case by case).

2.7 Quantifying uncertainty in parameter estimates

The usefulness of an estimate is always enhanced if some kind of measure of its precision can also be provided. Usually, this will be a standard error, an estimate of the standard deviation of the associated estimator. For the maximum likelihood estimator $\hat{\theta }$, a standard error is given by $$s.e.(\hat{\theta })=\frac{1}{\mathcal{I}(\hat{\theta }{)}^{\frac{1}{2}}},$$ and for a vector parameter $\mathit{\theta }$ $$s.e.({\hat{\theta }}_i)=[\mathcal{I}(\hat{\mathit{\theta }}{)}^{-1}{]}_{ii}^{\frac{1}{2}},i=1,\dots ,p.$$

An alternative summary of the information provided by the observed data about the location of a parameter $\theta$ and the associated precision is a confidence interval.

The asymptotic distribution of the maximum likelihood estimator can be used to provide approximate large sample confidence intervals. Asymptotically, ${\hat{\theta }}_i$ has a $N({\theta }_i,[\mathcal{I}(\mathit{\theta }{)}^{-1}{]}_{ii})$ distribution and we can find $z_{1-\frac{\alpha }{2}}$ such that $$P(-z_{1-\frac{\alpha }{2}}\le \frac{{\hat{\theta }}_i-{\theta }_i}{[\mathcal{I}(\mathit{\theta }{)}^{-1}{]}_{ii}^{\frac{1}{2}}}\le z_{1-\frac{\alpha }{2}})=1-\alpha .$$ Therefore $$P({\hat{\theta }}_i-z_{1-\frac{\alpha }{2}}[\mathcal{I}(\mathit{\theta }{)}^{-1}{]}_{ii}^{\frac{1}{2}}\le {\theta }_i\le {\hat{\theta }}_i+z_{1-\frac{\alpha }{2}}[\mathcal{I}(\mathit{\theta }{)}^{-1}{]}_{ii}^{\frac{1}{2}})=1-\alpha .$$ The endpoints of this interval cannot be evaluated because they also depend on the unknown parameter vector $\mathit{\theta }$. However, if we replace $\mathcal{I}(\mathit{\theta })$ by its MLE $\mathcal{I}(\hat{\mathit{\theta }})$ we obtain the approximate large sample $100(1-\alpha )\mathrm{\%}$ confidence interval $$[{\hat{\theta }}_i-z_{1-\frac{\alpha }{2}}[\mathcal{I}(\hat{\mathit{\theta }}{)}^{-1}{]}_{ii}^{\frac{1}{2}},{\hat{\theta }}_i+z_{1-\frac{\alpha }{2}}[\mathcal{I}(\hat{\mathit{\theta }}{)}^{-1}{]}_{ii}^{\frac{1}{2}}].$$ For $\alpha =0.1,0.05,0.01$, $z_{1-\frac{\alpha }{2}}=1.64,1.96,2.58$.

Example 2.9 (Bernoulli) If $y_1,\dots ,y_n$ are observations of $Y_1,\dots ,Y_n$, i.i.d. Bernoulli$(p)$ random variables then asymptotically $\hat{p}=\overline{y}$ has a $N(p,p(1-p)/n)$ distribution, and a large sample 95% confidence interval for $p$ is $$\begin{array}{rl}& [\hat{p}-1.96[\mathcal{I}(\hat{p}{)}^{-1}{]}^{\frac{1}{2}},\hat{p}+1.96[\mathcal{I}(\hat{p}{)}^{-1}{]}^{\frac{1}{2}}]\\ & =[\hat{p}-1.96[\hat{p}(1-\hat{p})/n{]}^{\frac{1}{2}},\hat{p}+1.96[\hat{p}(1-\hat{p})/n{]}^{\frac{1}{2}}]\\ & =[\overline{y}-1.96[\overline{y}(1-\overline{y})/n{]}^{\frac{1}{2}},\overline{y}+1.96[\overline{y}(1-\overline{y})/n{]}^{\frac{1}{2}}].\end{array}$$

2.8 Comparing statistical models

If we have a set of competing probability models which might have generated the observed data, we may want to determine which of the models is most appropriate. In practice, we proceed by comparing models pairwise. Suppose that we have two competing alternatives, $f_{\mathit{Y}}^{(0)}$ (model $M_0$) and $f_{\mathit{Y}}^{(1)}$ (model $M_1$) for $f_{\mathit{Y}},$ the joint distribution of $Y_1,\dots ,Y_n$. Often $H_0$ and $H_1$ both take the same parametric form, $f_{\mathit{Y}}(\mathit{y};\mathit{\theta })$ but with $\mathit{\theta }\in {\mathrm{\Theta }}^{(0)}$ for $H_0$ and $\mathit{\theta }\in {\mathrm{\Theta }}^{(1)}$ for $H_1$, where ${\mathrm{\Theta }}^{(0)}$ and ${\mathrm{\Theta }}^{(1)}$ are alternative sets of possible values for $\mathit{\theta }$. In the regression setting, we are often interested in determining which of a set of explanatory variables have an impact on the distribution of the response.

2.8.1 Hypothesis testing

A hypothesis test provides one mechanism for comparing two competing statistical models. A hypothesis test does not treat the two hypotheses (models) symmetrically. One hypothesis, $$H_0\text{: the data were generated from model }{M}_0,$$ is accorded special status, and referred to as the null hypothesis. The null hypothesis is the reference model, and will be assumed to be appropriate unless the observed data strongly indicate that $H_0$ is inappropriate, and that $$H_1\text{: the data were generated from model }{M}_1,$$ (the alternative hypothesis) should be preferred. The fact that a hypothesis test does not reject $H_0$ should not be taken as evidence that $H_0$ is true and $H_1$ is not, or that $H_0$ is better supported by the data than $H_1,$ merely that the data does not provide sufficient evidence to reject $H_0$ in favour of $H_1$.

A hypothesis test is defined by its critical region or rejection region, which we shall denote by $C$. $C$ is a subset of ${\mathbb{R}}^n$ and is the set of possible $\mathit{y}$ which would lead to rejection of $H_0$ in favour of $H_1$, i.e.

• If $\mathit{y}\in C$, $H_0$ is rejected in favour of $H_1$;
• If $\mathit{y}\notin C$, $H_0$ is not rejected.

As $\mathit{Y}$ is a random variable, there remains the possibility that a hypothesis test will produce an erroneous result. We define the size (or significance level) of the test $$\alpha =\underset{\mathit{\theta }\in {\mathrm{\Theta }}^{(0)}}{max}P(\mathit{Y}\in C;\mathit{\theta })$$ This is the maximum probability of erroneously rejecting $H_0$, over all possible distributions for $\mathit{Y}$ implied by $H_0$. We also define the power function $$\omega (\mathit{\theta })=P(\mathit{Y}\in C;\mathit{\theta })$$ It represents the probability of rejecting $H_0$ for a particular value of $\mathit{\theta }$. Clearly we would like to find a test with where $\omega (\mathit{\theta })$ is large for every $\mathit{\theta }\in {\mathrm{\Theta }}^{(1)}\setminus {\mathrm{\Theta }}^{(0)}$, while at the same time avoiding erroneous rejection of $H_0$. In other words, a good test will have small size, but large power.

The general hypothesis testing procedure is to fix $\alpha$ to be some small value (often 0.05), so that the probability of erroneous rejection of $H_0$ is limited. In doing this, we are giving $H_0$ precedence over $H_1$. Given our specified $\alpha$, we try to choose a test, defined by its rejection region $C$, to make $\omega (\mathit{\theta })$ as large as possible for $\mathit{\theta }\in {\mathrm{\Theta }}^{(1)}\setminus {\mathrm{\Theta }}^{(0)}$.

2.8.2 Likelihood ratio tests for nested hypotheses

Suppose that $H_0$ and $H_1$ both take the same parametric form, $f_{\mathit{Y}}(\mathit{y};\mathit{\theta })$ with $\mathit{\theta }\in {\mathrm{\Theta }}^{(0)}$ for $H_0$ and $\mathit{\theta }\in {\mathrm{\Theta }}^{(1)}$ for $H_1$, where ${\mathrm{\Theta }}^{(0)}$ and ${\mathrm{\Theta }}^{(1)}$ are alternative sets of possible values for $\mathit{\theta }$. A likelihood ratio test of $H_0$ against $H_1$ has a critical region of the form

$$\begin{array}{}\text{(2.3)}& C=\{\mathit{y}:\frac{\underset{\mathit{\theta }\in {\mathrm{\Theta }}^{(1)}}{max}L(\mathit{\theta })}{\underset{\mathit{\theta }\in {\mathrm{\Theta }}^{(0)}}{max}L(\mathit{\theta })}>k\}\end{array}$$

where $k$ is determined by $\alpha$, the size of the test, so $$\underset{\mathit{\theta }\in {\mathrm{\Theta }}^{(0)}}{max}P(\mathit{y}\in C;\mathit{\theta })=\alpha .$$ Therefore, we will only reject $H_0$ if $H_1$ offers a distribution for $Y_1,\dots ,Y_n$ which makes the observed data much more probable than any distribution under $H_0$. This is intuitively appealing and tends to produce good tests (large power) across a wide range of examples.

In order to determine $k$ in Equation 2.3, we need to know the distribution of the likelihood ratio, or an equivalent statistic, under $H_0$. In general, this will not be available to us. However, we can make use of an important asymptotic result.

First we notice that, as $\log$ is a strictly increasing function, the rejection region is equivalent to $$C=\{\mathit{y}:2\log \left(\frac{\underset{\mathit{\theta }\in {\mathrm{\Theta }}^{(1)}}{max}L(\mathit{\theta })}{\underset{\mathit{\theta }\in {\mathrm{\Theta }}^{(0)}}{max}L(\mathit{\theta })}\right)>k^{\prime }\}$$ where $$\underset{\mathit{\theta }\in {\mathrm{\Theta }}^{(0)}}{max}P(\mathit{y}\in C;\mathit{\theta })=\alpha .$$ Write $$L_{01}\equiv 2\log \left(\frac{\underset{\mathit{\theta }\in {\mathrm{\Theta }}^{(1)}}{max}L(\mathit{\theta })}{\underset{\mathit{\theta }\in {\mathrm{\Theta }}^{(0)}}{max}L(\mathit{\theta })}\right)$$ for the log-likelihood ratio test statistic. Provided that $H_0$ is nested within $H_1$, the following result provides a useful large-$n$ approximation to the distribution of $L_{01}$.

Theorem 2.3 Suppose that $H_0$: $\mathit{\theta }\in {\mathrm{\Theta }}^{(0)}$ and $H_1$: $\mathit{\theta }\in {\mathrm{\Theta }}^{(1)}$, where ${\mathrm{\Theta }}^{(0)}\subset {\mathrm{\Theta }}^{(1)}$. Let $d_0=\dim ({\mathrm{\Theta }}^{(0)})$ and $d_1=\dim ({\mathrm{\Theta }}^{(1)})$. Under $H_0$, the distribution of $L_{01}$ tends towards ${\chi }_{d_1-d_0}^2$ as $n\to \mathrm{\infty }$.

Proof. First we note that in the case where $\mathit{\theta }$ is one-dimensional and $\mathit{\theta }=(\theta )$, a Taylor series expansion of $\ell (\theta )$ around the MLE $\hat{\theta }$ gives $$\ell (\theta )=\ell (\hat{\theta })+(\theta -\hat{\theta })U(\hat{\theta })+\frac{1}{2}(\theta -\hat{\theta }{)}^2{U}^{\prime }(\hat{\theta })+\dots$$ Now, $U(\hat{\theta })=0$, and if we approximate $U^{\prime }(\hat{\theta })\equiv H(\hat{\theta })$ by $E[H(\theta )]\equiv -\mathcal{I}(\theta )$, and also ignore higher order terms, we obtain $$2[\ell (\hat{\theta })-\ell (\theta )]=(\theta -\hat{\theta }{)}^2\mathcal{I}(\theta )$$ As $\hat{\theta }$ is asymptotically $N[\theta ,\mathcal{I}(\theta {)}^{-1}]$, $(\theta -\hat{\theta }{)}^2\mathcal{I}(\theta )$ is asymptotically ${\chi }_1^2$, and hence so is $2[\ell (\hat{\theta })-\ell (\theta )]$.

Similarly it can be shown that when $\mathit{\theta }\in \mathrm{\Theta }$, a multidimensional space, $2[\ell (\hat{\mathit{\theta }})-\ell (\mathit{\theta })]$ is asymptotically ${\chi }_p^2$, where $p$ is the dimension of $\mathrm{\Theta }$.

Now, suppose that $H_0$ is true and $\mathit{\theta }\in {\mathrm{\Theta }}^{(0)}$ and therefore $\mathit{\theta }\in {\mathrm{\Theta }}^{(1)}$. Furthermore, suppose that $\ell (\mathit{\theta })$ is maximised in ${\mathrm{\Theta }}^{(0)}$ by ${\hat{\mathit{\theta }}}^{(0)}$ and is maximised in ${\mathrm{\Theta }}^{(1)}$ by ${\hat{\mathit{\theta }}}^{(1)}$. Then

$$\begin{array}{rl}{L}_{01}& \equiv 2\log \left(\frac{\underset{\mathit{\theta }\in {\mathrm{\Theta }}^{(1)}}{max}L(\mathit{\theta })}{\underset{\mathit{\theta }\in {\mathrm{\Theta }}^{(0)}}{max}L(\mathit{\theta })}\right)\\ & =2\log L({\hat{\mathit{\theta }}}^{(1)})-2\log L({\hat{\mathit{\theta }}}^{(0)})\\ & =2[\log L({\hat{\mathit{\theta }}}^{(1)})-\log L(\mathit{\theta })]-2[\log L({\hat{\mathit{\theta }}}^{(0)})-\log L(\mathit{\theta })]\\ & =L_1-L_0.\end{array}$$

Therefore $L_1=L_{01}+L_0$ and we know that, under $H_0$, $L_1$ has a ${\chi }_{d_1}^2$ distribution and $L_0$ has a ${\chi }_{d_0}^2$ distribution. Furthermore, it is possible to show (although we will not do so here) that under $H_0$, $L_{01}$ and $L_0$ are independent. It can also be shown that under $H_0$ the difference $L_1-L_0$ can be expressed as a quadratic form of normal random variables. Therefore, it follows that under $H_0$, the log likelihood ratio statistic $L_{01}$ has a ${\chi }_{d_1-d_0}^2$ distribution. ◻

Example 2.10 $y_1,\dots ,y_n$ are observations of $Y_1,\dots ,Y_n$, i.i.d. Bernoulli$(p)$ random variables. Suppose that we require a size $\alpha$ test of the hypothesis $H_0$: $p=p_0$ against the general alternative $H_1$: ‘$p$ is unrestricted’ where $\alpha$ and $p_0$ are specified.

Here $\mathit{\theta }=(p)$, ${\mathrm{\Theta }}^{(0)}=\{p_0\}$ and ${\mathrm{\Theta }}^{(1)}=(0,1)$ and the log likelihood ratio statistic is $$L_{01}=2n\overline{y}\log \left(\frac{\overline{y}}{p_0}\right)+2n(1-\overline{y})\log \left(\frac{1-\overline{y}}{1-p_0}\right).$$ As $d_1=1$ and $d_0=0$, under $H_0$, the log likelihood ratio statistic has an asymptotic ${\chi }_1^2$ distribution. For a log likelihood ratio test, we only reject $H_0$ in favour of $H_1$ when the test statistic is too large (observed data are much more probable under model $H_1$ than under model $H_0$), so in this case we reject $H_0$ when the observed value of the test statistic above is ‘too large’ to have come from a ${\chi }_1^2$ distribution. What we mean by ‘too large’ depends on the significance level $\alpha$ of the test. For example, if $\alpha =0.05$, a common choice, then we should reject $H_0$ if the test statistic is greater than the 3.84, the 95% quantile of the ${\chi }_1^2$ distribution.

2.8.3 Information criteria for model comparison

It is more difficult to use the likelihood ratio test of Section 2.8.2 to compare two models if those models are not nested. An alternative approach is to record some criterion measuring the quality of the model for each of a candidate set of models, then choose the model which is the best according to this criterion.

When we were estimating the unknown parameters $\theta$ of a model, we chose the value which maximised the likelihood: that is, the value of $\theta$ that maximises the probability of observing the data we actually saw. It is tempting to use a similar system for choosing between two models, and to choose the model which has the greater likelihood, under which the probability of seeing the data we actually observed is maximised. However, if we do this we will always end up choosing complicated models, which fit the observed data very closely, but do not meet our requirement of parsimony.

For a given model depending on parameters $\theta \in {\mathbb{R}}^p$, let $\hat{\ell }=\ell (\hat{\theta })$ be the log-likelihood function for that model evaluated at the MLE $\hat{\theta }$. It is not sensible to choose between models by maximising $\hat{\ell }$ directly, and instead it is common to choose a model to maximise a criteria of the form $$\hat{\ell }-\text{penalty},$$ where the penalty term will be large for complex models, and small for simple models.

Equivalently, we may choose between models by minimising a criteria of the form $$-2\hat{\ell }+\text{penalty}.$$ By convention, many commonly-used criteria for model comparison take this form. For instance, the Akaike information criterion (AIC) is $$\text{AIC}=-2\hat{\ell }+2p,$$ where $p$ is the dimension of the unknown parameter in the candidate model, and the Bayesian information criterion (BIC) is $$\text{BIC}=-2\hat{\ell }+\log (n)p,$$ where $n$ is the number of observations.