Nested Models and Likelihood Ratio Test

Suppose two models, $f_0(\cdot ,\overrightarrow{{\theta }^0})$ and $f_1(\cdot ,\overrightarrow{{\theta }^1})$ are under consideration. Then we say $f_0$ is nested in $f_1$ if the parameter space of $f_0$ is a subspace of the parameter space of $f_1$. Often, this can be demonstrated by showing that setting parameters in $f_1$ to particular values results in a model of the form $f_0$.

Nested model examples

• •

  The exponential model is nested within a Weibull model. The Weibull pdf is

  $$f(x;\lambda ,k)=\frac{k}{\lambda }{\left(\frac{x}{\lambda }\right)}^{k-1}{e}^{-{(x/\lambda )}^k}$$

  for $x\ge 0$, $\lambda >0$, $k>0$, where $\lambda$ is the scale parameter and $k$ is the shape parameter. However, if we set $k=1$ we recover the exponential pdf.

• •

  The exponential model is nested within a gamma model. The gamma pdf is

  $$f(x;\alpha ,\beta )=\frac{{\beta }^{\alpha }}{\mathrm{\Gamma }(\alpha )}{x}^{\alpha -1}{e}^{-\beta x},$$

  for $x\ge 0$,$\alpha >0$,$\beta >0$, where $\beta$ is the rate paramter and $\alpha$ is the shape parameter. Now, if we set $\alpha =1$ then we recover the exponential pdf (albeit with a slightly different parameterisation: let $\lambda ={\beta }^{-1}$).

• •

  The Uniform[0,1] distribution is nested within the beta distribution. The beta pdf is

  $$f(x;\alpha ,\beta )=\frac{x^{\alpha -1}{(1-x)}^{\beta -1}}{B(\alpha ,\beta )}$$

  for $0\le x\le 1$, $\alpha >0$, $\beta >0$, where $\alpha$ and $\beta$ are both shape parameters. If we set $\alpha =\beta =1$ this is the Uniform[0,1] distribution (NB $B(1,1)=1$).

• •

  The negative binomial distribution gives the number of failures $x$ observed until $r$ successes are observed. The geometric distribution is nested within the negative binomial. The negative binomial pmf is

  $$p(x;r,p)=\left(\genfrac{}{}{0pt}{}{x+r-1}{x}\right){(1-p)}^r{p}^x.$$

  Setting $r=1$ yields the geometric distribution.

• •

  The normal distribution with known variance is nested within a normal distribution with unknown variance.

If we are comparing two nested models for data $\overrightarrow{x}$ we can use the likelihood ratio test. The set-up is to view the choice between the two models as a hypothesis test. We take $f_0$ as a nested model within $f_1$ and take ${\mathrm{\Omega }}_i$ to be the parameter space of $f_i$. Then ${\mathrm{\Omega }}_0\subseteq {\mathrm{\Omega }}_1$.

Suppose we have the two models $f_0$ and $f_1$ as described above. Our hypotheses are

$H_0$: The simpler model adequately describes $\overrightarrow{x}$, i.e. ${\theta }_0\in {\mathrm{\Omega }}_0$,
$H_1$: The more complex model is required, i.e. ${\theta }_0\in \{{\mathrm{\Omega }}_1\setminus {\mathrm{\Omega }}_0\}$.

Then the deviance (likelihood ratio test) for model comparison is

$$D(f_1,f_0)=2\{{\mathrm{\ell }}_1({\widehat{\overrightarrow{\theta }}}^1)-{\mathrm{\ell }}_0({\widehat{\overrightarrow{\theta }}}^0)\}.$$

Under the usual regularity conditions, if $f_1$ has $p_1$ unknown parameters and $f_0$ has $p_0$ unknown parameters (of course $p_1>p_0$) then under $H_0$, asymptotically as $n\to \mathrm{\infty }$,

$$D(f_1,f_0)\sim {\chi }_{p_1-p_0}^2.$$

Remarks:

• •

  It should not surprise you to learn that this theorem is closely related to Theorems 2, 7 and 12. In fact, all three of these are corollaries to the above result. Convincing yourself of this would be a good way of enhancing your understanding.

• •

  It is always the case that ${\mathrm{\ell }}_1({\widehat{\overrightarrow{\theta }}}^1)>{\mathrm{\ell }}_0({\widehat{\overrightarrow{\theta }}}^0)$. To see this, remember that ${\widehat{\overrightarrow{\theta }}}^0$ is a valid parameter choice for model $f_1$, and ${\widehat{\overrightarrow{\theta }}}^1$ must have at least as high a likelihood under $f_1$ by virtue of it being the MLE.

Theorem 13 suggests that we should evaluate $D(f_1,f_0)$ and compare to the corresponding critical value, $z_c^2$ from the ${\chi }^2$ distribution. If we do not reject $H_0$ and take the simpler model $f_0$ as adequate.

Example 6.1:  Failure times, ctd.

We revisit the failure times example first encountered in Section 2.6. There we fitted an exponential model to the failure time data; should a gamma model be preferred?

From earlier,

$$\mathrm{\ell }(\theta )=n\log \theta -\theta n\overline{x};\widehat{\theta }={\overline{x}}^{-1}.$$

Using R we evaluate this at the MLE:

fail<-c(90,255,40,143,30,239,484,28,39,15)
exploglhd<-function(theta,x){
 n<-length(x)
 return(n*log(theta) - theta*n*mean(x))
}

exploglhd(1/mean(fail),fail)

This gives $\mathrm{\ell }(\widehat{\theta })=-59.14858$, with $\widehat{\theta }=0.00734$.

Now we consider the likelihood under a gamma distribution. Recall from Example 3.7 that maximising the gamma distribution reduces to solving

$$n\log \widehat{\alpha }-n\log \overline{x}+\sum _{i=1}^n\log x_i-n\gamma (\widehat{\alpha }),$$

and then

$$\widehat{\beta }=\frac{\widehat{\alpha }}{\overline{x}}.$$

The first equation can only be solved numerically.

Using R,

f<-function(alpha,x){
   n<-length(x)
   return(n*log(alpha)-n*log(mean(x))+sum(log(x))-n*digamma(alpha))
}

uniroot(f,lower=0.001,upper=10,x=fail)

giving $\widehat{\alpha }=1.013$, and so $\widehat{\beta }=0.00743$. (Note that the digamma function in R is the derivative of the log of the gamma function.)

Again, we evaluate the gamma log likelihood at its MLE in R.

gamloglhd<-function(alpha,beta,x){
    n<-length(x)
    return(n*alpha*log(beta) + (alpha-1)*sum(log(x))
Ψ- beta*sum(x) - n*log(gamma(alpha)))
}

gamloglhd(alpha=1.012568,beta=0.007428966,x=fail)

giving $\mathrm{\ell }(\widehat{\alpha },\widehat{\beta })=-59.14808$, only very slightly larger than the exponential distribution. This gives a deviance of $2(-59.14808--59.14858)=0.001$.

This is much less than the critical value of ${\chi }_1^2$ distribution, 3.84, hence we retain the simpler exponential model. We should not be surprised at our conclusion because $\widehat{\alpha }=1.013$, which is very close to 1, the value that would imply an exponential distribution.

Example 6.2:  Surgical Mortality Rates at Hospitals The below table gives mortality levels for cardiac surgery on babies at 12 hospitals. $r/m$ means $r$ deaths in $m$ operations.

A  0/47	B  18/148	C  8/119	D  46/810	E  8/211	F  13/196
G  9/148	H  31/215	I  14/207	J  8/97	K  29/256	L  24/360

Let ${\theta }_s$ be the mortality rate for hospital $s$. We would like to know whether mortality rates vary between hospitals. This can be expressed as the hypothesis test

$H_0$: ${\theta }_s={\theta }_{*}$ for all $s=A,\mathrm{\dots },L$.
$H_1$: Each ${\theta }_s$ is allowed to be different.

Clearly the model suggested by $H_0$ is nested within the model for $H_1$.

Letting $r_s$ denote the number of deaths in hospital $s$ and $m_s$ the total number of operations carried out at hospital $s$, the likelihood is

$$L(\overrightarrow{\theta })=\prod _{s=A}^L\left(\genfrac{}{}{0pt}{}{m_s}{r_s}\right){\theta }_s^{r_s}{(1-{\theta }_s)}^{m_s-r_s}.$$

Under $H_1$ this breaks down into separate Binomial distributions for each hospital, so we have

$${\widehat{\theta }}_s=\frac{r_s}{m_s}$$

for $s=A,\mathrm{\dots },L$.

Under $H_0$ this collapses into a single Binomial distribution where operations over all hospitals are added together, yielding

$${\widehat{\theta }}_{*}=\frac{\sum r_s}{\sum m_s}.$$

We use R to evaluate the likelihood under each scenario as follows:

####hospital mortality example
r<-c(0,18,8,46,8,13,9,31,14,8,29,24)
m<-c(47,148,119,810,211,196,148,215,207,97,256,360)

hosploglhd<-function(theta,r,m){
    return(log(prod(dbinom(r,m,theta))))
}

#h0:
hosploglhd(theta=sum(r)/sum(m),r,m)
#h1:
hosploglhd(theta=r/m,r,m)

This gives a log-likelihood of -44.15 for $H_0$ and -24.88 for $H_1$. The likelihood ratio statistic is then

$$2(-24.88--44.15)=38.54.$$

The model under $H_1$ has 11 parameters extra to the model under $H_0$; this gives a critical ${\chi }^2$ value of qchisq(0.95,df=11)=19.68. Since $38.54>19.68$ we reject $H_0$ and conclude that mortality rates vary between hospitals.

What is the interpretation of this?

• •

  Different case mix encountered by hospitals (e.g. specialised hospital taking more complex cases).

• •

  Some hospitals not as good as others.

Monitoring surgical mortality between hospitals (and between surgeons) is important to ensure that those performing below par are identified quickly.