4 The EF likelihood and ML estimation 4 The EF likelihood and ML estimation 4.2 The likelihood for a GLM

4.1 A Poisson example

Consider the AIDS data presented in Chapter 1. To write down the likelihood function we need to specify the pmf under the assumption of independent observations. We wish to fit the model:

Y_{i}|t_{i}\sim\mathrm{Pois}(\mu_{i}),\quad\mathrm{for}~{}~{}i=1,2,\dots,n% \quad\mathrm{with}\quad\log(\mu_{i})=a+bt_{i},

where $y_{i}$ is the number of observed deaths in the quarter $t_{i}$ . We denote the pmf by $f_{i}(y_{i}|\mu_{i})$ of a single observation $y_{i}$ . Here $\mu$ is a parametric function rather than a parameter itself. Furthermore, as the value of $\mu$ can change between observations due to the linear relation to the time covariate, so we do not assume that the observations are identically distributed.

 aids = read.table("aids.dat")
 attach(aids)
 plot(time, number)
 a = 0.340 ; b = 0.257 ; mu = exp(a + b * time)
 lines(time, mu)

This plot is presented in Section 1.3.2. The free parameters that determine the mean here are $a$ and $b$ .

Free parameters

There are several ways to specify which parameters are free and which are constrained. The moment parameter $\mu$ is always in one-to-one correspondence with the canonical parameter, e.g. $\theta=\log(\mu)$ . We may write the likelihood in terms of $\mu$ or in terms of $\theta$ as $f(y|\mu)=f(y|\theta)$ .

Baseline models

Taking the $n$ values of $\mu$ , $\mu_{i}$ , $i=1,2,\dots,n$ , to be free leads to the saturated model with $n$ free parameters. It is of theoretical interest as a baseline model. However, as it does not simplify the data, it is not usually a model to fit in its own right; $n$ unknowns for $n$ observations!

Setting the $n$ values $\mu_{i}$ , $i=1,2,\dots,n$ , to be equal gives a model with just $1$ free parameter, the null model.

Parsimony suggests a model with greater than 1 but less than $n$ free parameters. Here, constraining the moment parameters by $\log(\mu_{i})=a+bt_{i}$ , where $\beta=(a,b)$ leads to a model with two free parameters.

Exercise 4.37
The contribution to the likelihood of the $i$ th observation, $y_{i}$ , drawn from the Poisson distribution is one of the following:

\displaystyle\begin{array}[]{lrcl}(1)&\mathbb{P}(Y=y|\mu)&=&\mu^{y}\exp\{-\mu% \}\frac{1}{y!},\\ (2)&\mathbb{P}(Y=y_{i}|\mu)&=&\mu^{y_{i}}\exp\{-\mu\}\frac{1}{y_{i}!},\\ (3)&\mathbb{P}(Y_{i}=y_{i}|\mu_{i})&=&\mu_{i}^{y_{i}}\exp\{-\mu_{i}\}\frac{1}{% y_{i}!}.\end{array}

Which?

Model formula

With $\beta=(a,b)$ the log-likelihood of $\beta$ given observations $y_{1},y_{2},\dots,y_{n}$ is:

	$\displaystyle\ell(\beta)$	$\displaystyle=$	$\displaystyle\sum_{i=1}^{n}\{y_{i}\log\mu_{i}-\mu_{i}-\log(y_{i}!)\}$
		$\displaystyle=$	$\displaystyle\sum_{i=1}^{n}\{y_{i}(a+bt_{i})-\exp\{a+bt_{i}\}-\log(y_{i}!)\}.$

This is a function on $\beta\in\mathbb{R}^{2}$ .

Notation

We can consider this log-likelihood as a function of the vector of means $\boldsymbol{\mu}=(\mu_{1},\mu_{,}\dots,\mu_{n})$ , and write

\displaystyle\ell(\boldsymbol{\mu})

\displaystyle=

\displaystyle\sum_{i=1}^{n}\{y_{i}\log\mu_{i}-\mu_{i}-\log(y_{i}!)\}

where the means are constrained by

\log(\boldsymbol{\mu})\in\mathcal{M}=\{\log(\boldsymbol{\mu})=a\mathbf{1}+b% \mathbf{t},(a,b)\in\mathbb{R}^{2}\},

spanned by the $n$ -dimensional vectors $\mathbf{1}$ and $\mathbf{t}$ . Recall that $g=\log$ is the link function for the Poisson distribution.

This notation allows the log-likelihoods for each model to be simply expressed as

$\ell(\boldsymbol{\mu})$ constrained by $\mathcal{M}$ .

Here the null model is $\mathcal{M}=\{\log(\boldsymbol{\mu})=a\mathbf{1},a\in\mathbb{R}^{1}\}$ . The saturated model is $\mathcal{M}=\{\log(\boldsymbol{\mu})\in\mathbb{R}^{n}\}$ . And a 3-dimensional model is $\mathcal{M}=\{\log(\boldsymbol{\mu})=a\mathbf{1}+b\mathbf{t}+c\log(\mathbf{t})\}$ .

ML estimates

The partial derivative of the log-likelihood for each free parameter $a$ and $b$ are:

\frac{\partial{\ell(a,b)}}{\partial{a}}=\sum_{i=1}^{n}\left(y_{i}-\exp\{a+bt_{% i}\}\right)\quad\frac{\partial{\ell(a,b)}}{\partial{b}}=\sum_{i=1}^{n}\left(y_% {i}t_{i}-t_{i}\exp\{a+bt_{i}\}\right)

Except in special cases, it is not possible to analytically solve the above equations to determine the maximum likelihood estimates $\hat{a}$ and $\hat{b}$ . We therefore need to use a numerical technique to find the maximum likelihood estimate of the parameters. We shall examine this numerical technique in the next chapter. In the meantime, we can use this algorithm in R via the glm() command.

The necessary arguments for the glm() command are:

formula: A symbolic description of the model to be fitted. E.g. number ~ 1 + time.
family: A description of the error distribution and link function to be used in the model. E.g. poisson(link="log")
data: A data frame containing the variables of the model. E.g. aids.

We can estimate the free parameters by:

 AIDSModel <- glm(number~1+time, family=poisson(link="log"), data=aids)
 AIDSModel

This will return

Call:  glm(formula = number ~ 1 + time, family = poisson(link = "log"),
    data = aids)

Coefficients:
(Intercept)         time
     0.3396       0.2565

Degrees of Freedom: 13 Total (i.e. Null);  12 Residual
Null Deviance:      207.3
Residual Deviance: 29.65        AIC: 86.58

The ML estimates $\hat{a}=0.340$ and $\hat{b}=0.257$ maximise the Poisson regression model $\mathcal{M}=\{\log(\boldsymbol{\mu})=a\mathbf{1}+b\mathbf{t}\}$ .