```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,tidy=TRUE,message=FALSE,warning=FALSE)
showsol<-FALSE
```


# Overview

## Learning material

* Textbook: Fahrmeir et al (2013): Chapter 2.3, 5.1, B4.1-3
* [Classnotes 13.09.2018](https://www.math.ntnu.no/emner/TMA4315/2018h/TMA4315M3H20180913.pdf) 
* [Classnotes 20.09.2018](https://www.math.ntnu.no/emner/TMA4315/2018h/TMA4315M3H20180920.pdf) 

---

## Topics
### [First week](#firstweek)

* aim of binary regression
* how to model a binary response 
* three ingredients of a GLM model
* the logit model: logistic regression
* interpreting the logit model - with odds
* grouped vs. individual data
* parameter estimation with maximum likelihood
    + likelihood, log-likelihood, 
    + score function


---

### [Second week](#secondweek)

* Parameter estimation
    + score function- and mean and covariance thereof, 
    + observed and expected information matrix
* comparison with the normal distribution - score function and Fisher information
* exponential family and canonical link
* iterative calculation of ML estimator (Newton-Raphson and Fisher scoring) - and in R with `optim`
* asymptotic properties of ML estimators - how to use in inference?
* statistical inference
    + confidence intervals
    + hypothesis testing: Wald, and likelihood ratio
* deviance: definition, analysis of deviance, deviance residuals 
* model fit and model choice
* overdispersion and estimating overdispersion parameter
* sampling strategy: cohort, but also case-control data good for logit model


---

## SECOND WEEK

Remember the [beetle](#beetle1) and [infant respitory disease](#infant1) examples?

First, we look back at the [model requirements](#binaryregassump) for the binary regression - and the [loglikelihood](#loglik) and [score function](#score). 

---

# Likelihood and derivations thereof - continued

Individual data (not grouped):

Loglikelihood:
$$l(\boldsymbol{\beta})=\sum_{i=1}^n[y_i \ln \pi_i-y_i\ln(1-\pi_i)+\ln(1-\pi_i)]$$

Score function:
$$s(\boldsymbol{\beta})=\sum_{i=1}^n {\bf x}_i (y_i-\pi_i)$$

---

## Properties of the score function

Since the score function depends on $Y_i=y_i$ we may regard the score function as a random vector. We will now calculate the mean and covariance matrix for the score function. 

---

## $E(s(\boldsymbol{\beta}))$

The expected value of the score function is

$$
\begin{aligned}
E(s(\boldsymbol{\beta})) &= E(\sum_{i=1}^n(Y_i-\pi_i){\bf x}_i) \\
&= \sum_{i=1}^nE((Y_i-\pi_i){\bf x}_i) \\
&= \sum_{i=1}^n(E(Y_i)-\pi_i){\bf x}_i = 0 \\
\end{aligned}
$$

as $E(Y_i) = \pi_i$. 

We also see that $E(s_i(\boldsymbol{\beta})) = E((Y_i-\pi_i){\bf x}_i) = 0, \forall i$.


---

## Fisher Information and Variances of Estimates

The "amount of information" that the data carry about the parameters, $\boldsymbol{\beta}$, can be summarised by the curvature in the likelihood surface

```{r Curvature, echo=FALSE, fig.height=4.5}
beta0 <- seq(-4,4,length=50)
data <- rbinom(100, 1, 0.45)
CalcBinLhood <- function(b0, n, r) {
  eta <- 1/(1+exp(-b0))
  sum(dbinom(r, n, eta, log=TRUE))
}
lhoodLots <- sapply(beta0, CalcBinLhood, n=100, r=40)
lhoodLittle <- sapply(beta0, CalcBinLhood, n=10, r=4)
plot(beta0, lhoodLots-max(lhoodLots), type="l", ylim=c(-100,0), 
     xlab=expression(beta[0]), ylab="log likelihood")
lines(beta0, lhoodLittle-max(lhoodLittle), col=2)
legend(-2, -30, c("n=100", "n=10"), lty=1, col=1:2)

```

The more curvature, the more certain we are = the lower the variance

---

## Fisher information matrix

We can formalise this with the observed and expected Fisher information matrices

Observed: $H(\boldsymbol{\beta}) = -\frac{\partial^2l(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}\partial\boldsymbol{\beta}^T}$ - curvature

Expected: $F(\boldsymbol{\beta}) = \text{Cov}(s(\boldsymbol{\beta})) = \sum_{i=1}^n \text{Cov}(s_i(\boldsymbol{\beta}))$ - variance

We will calculate these for the binomial with a logit link, and show they are the same, i.e.  under mild regularity conditions (so we can change the order of $\int$ and $\frac{\partial}{\partial \boldsymbol{\beta}}$):

$$\text{Cov}(s(\boldsymbol{\beta}))=F(\boldsymbol{\beta}) = E\left( -\frac{\partial^2l(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}\partial\boldsymbol{\beta}^T} \right) = E(-\text{Hessian matrix of }l)$$


---

## Observed Fisher information matrix $H(\boldsymbol{\beta})$

What is the observed Fisher information matrix?

$$H(\boldsymbol{\beta}) = -\frac{\partial^2l(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}\partial\boldsymbol{\beta}^T} = -\frac{\partial s(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}^T} = \frac{\partial}{\partial\boldsymbol{\beta}^T}\left[\sum_{i=1}^n (\pi_i-y_i){\bf x}_i \right] $$

because $s(\boldsymbol{\beta}) = \sum_{i=1}^n (y_i-\pi_i){\bf x}_i$, so $-s(\boldsymbol{\beta}) = \sum_{i=1}^n (\pi_i-y_i){\bf x}_i$. 

Note that $\pi_i = \pi_i(\boldsymbol{\beta})$.

---

## Calculating $H(\boldsymbol{\beta})$

$$H(\boldsymbol{\beta}) = \sum_{i=1}^n \frac{\partial}{\partial\boldsymbol{\beta}^T}[{\bf x}_i\pi_i-{\bf x}_iy_i] = \sum_{i=1}^n \frac{\partial}{\partial\boldsymbol{\beta}^T}{\bf x}_i\pi_i = \sum_{i=1}^n {\bf x}_i \frac{\partial \pi_i}{\partial \eta_i} \frac{\partial \eta_i}{\partial \boldsymbol{\beta}^T} $$

---

Use that

$$ \frac{\partial \eta_i}{\partial \boldsymbol{\beta}^T}=\frac{\partial {\bf x}_i^T\boldsymbol{\beta}}{\partial \boldsymbol{\beta}^T} = \left(\frac{\partial {\bf x}_i^T\boldsymbol{\beta}}{\partial \boldsymbol{\beta}}\right)^T = {\bf x}_i^T $$

and 

$$
\begin{aligned}
\frac{\partial \pi_i}{\partial \eta_i} &=  \frac{\partial\left(\frac{\exp(\eta_i)}{1+\exp(\eta_i)}\right)}{\partial \eta_i}\\ &= \frac{(1+\exp(\eta_i))\exp(\eta_i) - \exp(\eta_i)\exp(\eta_i)}{(1+\exp(\eta_i))^2} \\
&= \frac{\exp(\eta_i)}{1+\exp(\eta_i)} \frac{1+\exp(\eta_i)-\exp(\eta_i)}{1+\exp(\eta_i)} \\
&=\pi_i (1-\pi_i).
\end{aligned}
$$

---

And thus

$$
\begin{aligned}
H(\boldsymbol{\beta}) &=\sum_{i=1}^n {\bf x}_i \frac{\partial \pi_i}{\partial \eta_i} \frac{\partial \eta_i}{\partial \boldsymbol{\beta}^T} \\
&=\sum_{i=1}^n {\bf x}_i\pi_i(1-\pi_i) {\bf x}_i^T\\
&= \sum_{i=1}^n {\bf x}_i{\bf x}_i^T \pi_i(1-\pi_i).
\end{aligned}
$$


---

## The expected Fisher information matrix $F(\boldsymbol{\beta})$

The covariance of $s(\boldsymbol{\beta})$ is called the expected Fisher information matrix, $F(\boldsymbol{\beta})$ and is given by


$$
\begin{aligned} 
F(\boldsymbol{\beta}) &= \text{Cov}(s(\boldsymbol{\beta})) = \sum_{i=1}^n \text{Cov}(s_i(\boldsymbol{\beta})) \\
&= \sum_{i=1}^n E\left[\Big(s_i(\boldsymbol{\beta}) - E(s_i(\boldsymbol{\beta}))\Big)\Big(s_i(\boldsymbol{\beta})-E(s_i(\boldsymbol{\beta}))\Big)^T\right] \\
&= \sum_{i=1}^n E(s_i(\boldsymbol{\beta})s_i(\boldsymbol{\beta})^T) = \sum_{i=1}^n F_i(\boldsymbol{\beta}) \end{aligned}
$$


assuming that the responses $Y_i$ and $Y_j$ are independent

---

## $F_i(\boldsymbol{\beta})$

Remember that $s_i(\boldsymbol{\beta})=(Y_i-\pi_i){\bf x}_i$, then:

$$
\begin{aligned}
F_i(\boldsymbol{\beta}) &= E(s_i(\boldsymbol{\beta})s_i(\boldsymbol{\beta})^T) = E((Y_i-\pi_i){\bf x}_i(Y_i-\pi_i){\bf x}_i^T) \\
&= {\bf x}_i{\bf x}_i^T E((Y_i-\pi_i)^2)\\ 
&= {\bf x}_i{\bf x}_i^T \pi_i(1-\pi_i) 
\end{aligned}
$$
where $E((Y_i-\pi_i)^2)=\text{Var}(Y_i)=\pi_i(1-\pi_i)$ is the variance of $Y_i$. Thus

$$F(\boldsymbol{\beta}) = \sum_{i=1}^n {\bf x}_i{\bf x}_i^T \pi_i(1-\pi_i).$$

So the observed and the expected Fisher information matrix are equal., i.e. $Cov(s(\beta) = -\frac{\partial^2l(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}\partial\boldsymbol{\beta}^T}$

This is not the case for the probit or complementary log-log models.

---

## Overview of the results for individual and grouped data

* Individual data: $i=1,\ldots, n$, and pairs $({\bf x}_i,y_i)$.
* Grouped data: $j=1,\ldots, G$ with $n_j$ observations for group $j$, and $Y_j=\sum Y_i$ for all $i$ member of group $j$. In total $\sum_{j=1}^G n_j$ observations. For each pair $({\bf x}_j, y_j)$, where ${\bf x}_j$ the covariate pattern for group $j$.

NB: we keep that $\eta_i=\ln (\frac{\pi_i}{1-\pi_i})$ - not changed for grouped data (but now $\mu_j=n_j\pi_j$).

---

**Log-likelihood:**

Individual:
$$l(\boldsymbol{\beta})=\sum_{i=1}^n[y_i \ln \pi_i-y_i\ln(1-\pi_i)+\ln(1-\pi_i)]$$
Grouped:
$$l(\boldsymbol{\beta})=\sum_{j=1}^G[y_j \ln \pi_j-y_j\ln(1-\pi_j)+n_j\ln(1-\pi_j)+ \ln {n_j \choose y_j}]$$
The last part is usually not include in calculations.

---

**Score function:**

Individual:
$$s(\boldsymbol{\beta})=\sum_{i=1}^n {\bf x}_i (y_i-\pi_i)$$

Grouped:
$$s(\boldsymbol{\beta})=\sum_{j=1}^G {\bf x}_j (y_j-n_j\pi_j)$$

---

**Expected Fisher information matrix:**

Individual:

$$F(\boldsymbol{\beta})=\sum_{i=1}^n {\bf x}_i {\bf x}_i^T\pi_i (1-\pi_i)$$

Grouped:

$$F(\boldsymbol{\beta})=\sum_{j=1}^G {\bf x}_j {\bf x}_j^T n_j \pi_j (1-\pi_j)$$

The observed Fisher information matrix equals the expected Fisher information matrix - because the logit model is the _canonical link_ for the binomial distribution.

---

## Fitting the Models

There is no analytic solution, so we have to resort to numerics. Luckily, these models behave well enough

---


## Iterative gradient-based methods

Use a first order multivariate Taylor approximation for $s(\hat{\boldsymbol{\beta}})$, around some chosen reference value $\boldsymbol{\beta}^{(0)}$:

$$
s(\hat{\boldsymbol{\beta}})\approx s(\boldsymbol{\beta}^{(0)})+\frac{\partial s(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\big\rvert_{\boldsymbol{\beta}=\boldsymbol{\beta}^{(0)}} (\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{(0)})
$$

Let $H(\boldsymbol{\beta}^{(0)})=-\frac{\partial s(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\big\rvert_{\boldsymbol{\beta}=\boldsymbol{\beta}^{(0)}}$.
Setting $s(\hat{\boldsymbol{\beta}})=0$ solving for $\hat{\boldsymbol{\beta}}$ gives 

$$
\hat{\boldsymbol{\beta}}=\boldsymbol{\beta}^{(0)} + H(\boldsymbol{\beta}^{(0)})^{-1} s(\boldsymbol{\beta}^{(0)})
$$
where $H(\boldsymbol{\beta}^{(0)})^{-1}$ is the matrix inverse of $H(\boldsymbol{\beta}^{(0)})$.

---

## Enter Newton and Raphson

If we start with some value $\boldsymbol{\beta}^{(0)}$ and then find a new value $\boldsymbol{\beta}^{(1)}$ by applying this equation, and then continue applying the equation until convergence we have the _Newton-Raphson_ method:

$$\boldsymbol{\beta}^{(t+1)}=\boldsymbol{\beta}^{(t)} + H(\boldsymbol{\beta}^{(t)})^{-1} s(\boldsymbol{\beta}^{(t)})$$

---

Replacing the observed Fisher information matrix ${\bf H}$ with the expected Fisher information matrix ${\bf F}$ yields the _Fisher-scoring_ method.

For the logit model these two methods are the same since the observed and expected Fisher information matrix is the same for canonical link functions (like the logit is for binary regression).

This algorithm is run until the relative difference in Euclidean distance between two iterations "(new-old)/old" is smaller than some chosen constant. 

---

## Requirements for convergence

For the Newton-Raphson algorithm we see that the observed Fisher information matrix $H$ needs to be invertible for all $\boldsymbol{\beta}$, alternatively for the Fisher scoring algorithm the expected Fisher information matrix $F$ needs to be invertible.

---

## Proof of convergence

In our logit model $$F(\boldsymbol{\beta})=\sum_{j=1}^G {\bf x}_j {\bf x}_j^T n_j \pi_j (1-\pi_j)$$ 

Let ${\bf X}$ be the design matrix, where the rows are ${\bf x}_j^T$. Then ${\bf X}^T{\bf X}=\sum_{j=1}^G {\bf x}_j {\bf x}_j^T$. 

If we require that the design matrix has full rank ($G$) then also ${\bf X}^T{\bf X}$ will have full rank (it will also be positive definite) and and in addition $\pi_j(1-\pi_j)>0$ for all $\pi_j \in (0,1)$, so then $F(\boldsymbol{\beta})$ will be positive definite and all is good. 

---

## Why is $F(\boldsymbol{\beta})$ positive definite if we require that the design matrix has full rank? 

Formally, let ${\bf X}$ be a $n\times p$ matrix and $\Lambda$ a $n\times n$ diagonal matrix where all the diagonal elements are positive (like our $\pi_j(1-\pi_j)$, yes, put them on the diagonal). 
Let ${\bf X}$ be of full rank $\Leftrightarrow$ ${\bf X}^T\Lambda{\bf X}$ is positive definite.


*Proof:* 

$\Rightarrow$: Let ${\bf v}$ be a $p$ dimensional column vector. Assume $0={\bf v}^T {\bf X}^T \Lambda {\bf X} {\bf v}=(\Lambda^{1/2}{\bf X}{\bf v})^T(\Lambda^{1/2}{\bf X}{\bf v})=\sum_{i=1}^n w_i^2$ where ${\bf W}=\Lambda^{1/2}{\bf X}{\bf v}$. Then, $w$ must be 0, that is $\Lambda^{1/2}{\bf X}{\bf v}=0$ since multiplication with $\Lambda^{1/2}$ is to multiply each element in ${\bf X}{\bf v}$ with a number different from 0. That is, we must have ${\bf v}={\bf 0}$ since ${\bf X}$ has independent columns.

$\Leftarrow$: Assume that ${\bf X}{\bf v}={\bf 0}$. Then ${\bf v}^T {\bf X}^T \Lambda {\bf X} {\bf v}={\bf 0}$ so ${\bf v}={\bf 0}$ since ${\bf X}^T\Lambda{\bf X}$ is positive definite. This is, ${\bf X}$ has independent columns.

*End of proof*

---

## We need a full rank

Therefore, for GLMs we will also - as for the multiple linear regression model in Module 2 - assume that the design matrix has full rank!

We will see in Module 5 that this is the requirement needed for GLMs in general.

---

## Convergence

Convergence is still not guarnateed, especially for small samples. 

According to our text book, Fahrmeir et al (2013), page 284, the conditions for uniqueness and existence of ML estimators are very complex, and the authors suggest that the GLM user instead checks for convergence in practice by performing the iterations.

In practice, the logit model most often causes problems, when (for grouped data) $y_i=0$ or $y_i=n_i$, because $\hat{\pi}_i=0/1$, so $\hat{\eta}_i=\pm \infty$. Computers do not like infinity

---


# Asymptotic properties of ML estimates

## Results

We need some weak regularity conditions, including 

- $\boldsymbol{\beta}$ falls in the interior of the parameter space and 
- $p$ is fixed that $n$ increases

(Agresti (2015) page 125):

--- 

## The Results

Let $\hat{\boldsymbol{\beta}}$ be the maximum likelihood (ML) estimate in the GLM model. As the total sample size increases, $n\rightarrow \infty$:

1. $\hat{\boldsymbol{\beta}}$ exists
2. $\hat{\boldsymbol{\beta}}$ is consistent (convergence in probability, yielding asymptotically unbiased estimator, variances goes towards 0)
3. $\hat{\boldsymbol{\beta}} \approx N_p(\boldsymbol{\beta},F^{-1}(\hat{\boldsymbol{\beta}}))$ 

So asymptotically $\text{Cov}(\hat{\boldsymbol{\beta}})=F^{-1}(\hat{\boldsymbol{\beta}})$: the inverse of the expected Fisher information matrix evaluated at the ML estimate. 


---

The _proof_ (for the univariate case) is given in the course TMA4295 Statistical Inference course, Casella and Berger (2002):"Statistical inference", page 472. 

Here we will sketch the proof. The strategy:

- make a first order Taylor expansion of the score function around the true parameter, 
- use the fact that the maximum likelihood estimate is defined as the zero of the score function. 

---

## The sketch

Use $\boldsymbol{\theta}$ as the parameter of interest

(there is the connection to $\boldsymbol{\mu}$, then to $\boldsymbol{\eta}$ and finally to $\boldsymbol{\beta}$)

We start with the multivariate version of the first order Taylor expansion of the score around the true parameter value $\boldsymbol{\theta}$:

$$s({\boldsymbol{\theta}}) \approx s(\hat{\boldsymbol{\theta}}) + s'({\boldsymbol{\theta}})(\hat{\boldsymbol{\theta}}-{\boldsymbol{\theta}})$$
As $s'({\boldsymbol{\theta}}) = {\bf H}({\boldsymbol{\theta}})$, and $s(\hat{\boldsymbol{\theta}})=0$, 

$$s(\hat{\boldsymbol{\theta}})\approx s({\boldsymbol{\theta}})-{\bf H}({\boldsymbol{\theta}})(\hat{\boldsymbol{\theta}}-{\boldsymbol{\theta}}) = 0$$

$$s({\boldsymbol{\theta}})\approx {\bf H}({\boldsymbol{\theta}})(\hat{\boldsymbol{\theta}}-{\boldsymbol{\theta}})$$
And premultiplying with ${\bf H}^{-1}(\boldsymbol{\theta})$ gives
$$ (\hat{\boldsymbol{\theta}}-{\boldsymbol{\theta}})\approx {\bf H}^{-1}({\boldsymbol{\theta}})s({\boldsymbol{\theta}})$$

---

Then, to use the central limit theorem we need some smart manipulations with $n$, so we start by multiplying with $\sqrt{n}$ and split that into $n$ and $\frac{1}{\sqrt{n}}$.

$$ \sqrt{n}(\hat{\boldsymbol{\theta}}-{\boldsymbol{\theta}})\approx \sqrt{n}{\bf H}^{-1}({\boldsymbol{\theta}})s({\boldsymbol{\theta}})=(\frac{1}{n}{\bf H}({\boldsymbol{\theta}}))^{-1}\frac{1}{\sqrt{n}}s({\boldsymbol{\theta}})$$


From the central limit theorem:

1) $\frac{1}{n}{\bf H}(\boldsymbol{\theta})$ goes to the expected value which is ${\bf F}(\boldsymbol{\theta})$ (in probability),

2) the part $\frac{1}{\sqrt{n}}s({\boldsymbol{\theta}})$ asymptotically goes to a random variable $W$ that follows a multivariate normal with ${\bf W} \sim N({\bf 0},\frac{1}{n}{\bf F}(\boldsymbol{\theta}))$:
   + mean $\text{E}(\frac{1}{\sqrt{n}}s({\boldsymbol{\theta}}))={\bf 0}$ and the 
   + covariance matrix is $\text{Cov}(\frac{1}{\sqrt{n}}s({\boldsymbol{\theta}}))=\frac{1}{n}{\bf F}(\boldsymbol{\theta})$

---

$$
\sqrt{n}(\hat{\boldsymbol{\theta}}-{\boldsymbol{\theta}}) \approx \bf{F}^{-1}(\boldsymbol{\theta}){\bf W}
$$

On the right side here we have a multivariate normal distributed random variable $\bf{F}^{-1}(\boldsymbol{\theta}){\bf W}$ with mean ${\bf 0}$ and covariance matrix 
$$\text{Cov}({\bf F}^{-1}(\boldsymbol{\theta}){\bf W})={\bf F}^{-1}(\boldsymbol{\theta})\frac{1}{n}{\bf F}(\boldsymbol{\theta}){\bf F}^{-1}(\boldsymbol{\theta})=\frac{1}{n}{\bf F}^{-1}(\boldsymbol{\theta})$$

This leads to the wanted result:

$$ \hat{\boldsymbol{\theta}}\approx N(\boldsymbol{\theta},{\bf F}^{-1}(\boldsymbol{\theta}))$$

Due to the Slutsky theorem (from TMA4295 Statistical inference) this also holds when ${\bf F}^{-1}(\boldsymbol{\theta}))$ is replaced by ${\bf F}^{-1}(\hat{\boldsymbol{\theta}}))$.

---

# Parameter estimation

Parameter estimation can be based on grouped data - so now we use $Y_j \sim \text{bin}(n_j,\pi_j)$ from 1 above, but keep 2 and 3 unchanged. The number of groups is $G$ and the total number of observations is $\sum_{j=1}^G
n_j$.

* Likelihood=joint distribution, exponential family.
$$ f(y\mid \theta)=\exp \left( \frac{y \theta-b(\theta)}{\phi}\cdot w + c(y, \phi, w) \right) $$
where we have that $\theta=\ln(\frac{\pi}{1-\pi})$ for the binomial distribution, which means that our logit model is gives the canonical link (remember, good properties!).

* Log-likelihood
$$l(\boldsymbol{\beta})=\sum_{j=1}^G[y_j \ln \pi_j-y_j\ln(1-\pi_j)+n_j\ln(1-\pi_j)+ \ln {n_j \choose y_j}]$$

---

* Score function=vector of partial derivatives of log-likelihood. Find ML by solving $s(\hat{\boldsymbol{\beta})})=0$ - but no closed form solutions.
$$s(\boldsymbol{\beta})=\sum_{j=1}^G {\bf x}_j (y_j-n_j\pi_j)$$
* Expected Fisher information matrix
$$F(\boldsymbol{\beta})=\sum_{j=1}^G {\bf x}_j {\bf x}_j^T n_j \pi_j (1-\pi_j)$$
* $\hat{\boldsymbol{\beta}}$ found iteratively using Newton-Raphson or Fisher scoring
$$\boldsymbol{\beta}^{(t+1)}=\boldsymbol{\beta}^{(t)} + F(\boldsymbol{\beta}^{(t)})^{-1} s(\boldsymbol{\beta}^{(t)})$$

* $\hat{\boldsymbol{\beta}} \approx N_p(\boldsymbol{\beta},F^{-1}(\hat{\boldsymbol{\beta}}))$ 

# Further statistical inference

Our further statistical inference (confidence intervals and hypotheses tests) are based on the asymptotic distribution of the parameter estimates
$$\hat{\boldsymbol{\beta}} \approx N_p(\boldsymbol{\beta},F^{-1}(\hat{\boldsymbol{\beta}}))$$
where $F^{-1}(\hat{\boldsymbol{\beta}}))$ is the inverse of the expected Fisher information matrix inserted $\hat{\boldsymbol{\beta}}$.

For the logit model we found that 
$$F(\boldsymbol{\beta})=\sum_{j=1}^G {\bf x}_j {\bf x}_j^T n_j \pi_j (1-\pi_j)$$
So we would need to do $\pi_j=\frac{\exp(\eta_j)}{1+\exp(\eta_j)}$ and $\eta_j={\bf x}_j^T\boldsymbol{\beta}$ as "usual", and then replace $\boldsymbol{\beta}$ with $\hat{\boldsymbol{\beta}}$. 

The asymptotic distribution still holds when we replace $\boldsymbol{\beta}$ with $\hat{\boldsymbol{\beta}}$ in ${\bf F}$.

---

If we make a diagonal matrix ${\bf W}$ with $n_j \pi_j (1-\pi_j)$ on the diagonal, then we may write the matrix $F(\boldsymbol{\beta})$ in matrix notation. As before ${\bf X}$ is the $G\times p$ design matrix.

$$F(\boldsymbol{\beta})=\sum_{j=1}^G {\bf x}_j {\bf x}_j^T n_j \pi_j (1-\pi_j)={\bf X}^T {\bf W} {\bf X}.$$
which means that $\text{Cov}(\hat{\boldsymbol{\beta}})=({\bf X}^T {\bf W} {\bf X})^{-1}$ for the binomial model (remember that $\hat{\boldsymbol{\beta}}$ comes in with $\hat{\pi}_j$ in ${\bf W}$).

**Q:** How is this compared to the normal case? 

A: $F(\boldsymbol{\beta})=\frac{1}{\sigma^2}{\bf X}^T{\bf X}$, and the inverse $\text{Cov}(\hat{\boldsymbol{\beta}})=({\bf X}^T {\bf X})^{-1}\sigma^2$.

---

Let ${\bf A}(\boldsymbol{\beta})=F^{-1}(\hat{\boldsymbol{\beta}})$, and $a_{kk}(\hat{\boldsymbol{\beta}})$ is diagonal element number $k$.

For one element of the parameter vector:
$$ Z_k=\frac{\hat{\boldsymbol{\beta}}_k-\boldsymbol{\beta}_k}{\sqrt{\hat{a}_{kk}({\boldsymbol{\beta}})}}$$
is asymptotically standard normal. We will use this now!

---

<!---

**Q**: Where can you find $\hat{\boldsymbol{\beta}}$ and $F^{-1}(\hat{\boldsymbol{\beta}})$ in the print-out below?

---

```{r}
library(investr)
fitgrouped=glm(cbind(y, n-y) ~ ldose, family = "binomial", data = investr::beetle) 
summary(fitgrouped)
summary(fitgrouped)$cov.scaled
sqrt(diag(summary(fitgrouped)$cov.scaled))
```
Note: `cov.unscaled` is the estimated covariance matrix for the estimated coefficients, and `cov.scaled` is the the `cov.unscaled` scaled with the dispersion parameter. For the binomial the dispersion is equal to 1, so no difference between the two.
--->


## Confidence intervals

In addition to providing a parameter estimate for each element of our parameter vector $\boldsymbol{\beta}$ we should also report a $(1-\alpha)100$% confidence interval (CI) for each element.

We focus on element $k$ of $\boldsymbol{\beta}$, called $\boldsymbol{\beta}_k$. It is known that asympotically 

$$
Z_k=\frac{\hat{\boldsymbol{\beta}}_k-\boldsymbol{\beta}_k}{\sqrt{ {a}_{kk}(\hat{\beta})}} \sim N(0,1)
$$

We use that to form confidence intervals.

Let $z_{\alpha/2}$ be such that $P(Z_k>z_{\alpha/2})=\alpha/2$.

---

We then use 
$$ P(-z_{\alpha/2}\le Z_k \le z_{\alpha/2})=1-\alpha$$
insert $Z_k$ and solve for $\beta_k$ to get
$$ P(\hat{\beta}_k-z_{\alpha/2}\sqrt{a_{kk}(\hat{\beta})}
\le \beta_k \le \hat{\beta}_k-z_{\alpha/2}\sqrt{a_{kk}(\hat{\boldsymbol{\beta}})})=1-\alpha$$

A $(1-\alpha)$% CI for $\beta_k$ is when we insert numerical values for the upper and lower limits.

**Q:** We write $a_{kk}(\hat{\boldsymbol{\beta}})$. Why not $a_{kk}(\hat{\beta}_{kk})$?

---

### Example with the beetle data

Again, we study our beetle data - in the grouped version.

Here we calculate the upper and lower limits of the confidence interval using the formula.

```{r,echo=-1}
library(investr)
fitgrouped=glm(cbind(y, n-y) ~ ldose, family = "binomial", data = investr::beetle) 
coeff=fitgrouped$coefficients
sds=sqrt(diag(summary(fitgrouped)$cov.scaled))
alpha=0.05
lower=coeff-qnorm(1-alpha/2)*sds
upper=coeff+qnorm(1-alpha/2)*sds
cbind(lower,upper)
```

**Q:** Explain what is done in the `R`-print-out.

---

## Hypothesis testing

There are three methods that are mainly used for testing hypotheses in GLMs:

- Wald test, 
- likelihood ratio test and 
- score test. 

We will look at the first two.

---

First, look at linear hypotheses: We study a binary regression model with $p=k+1$ covariates, and refer to this as model A (the larger model). As for the multiple linear model we then want  to investigate the null and alternative hypotheses of the following type(s):

\begin{eqnarray*}
 H_0: \beta_{j}&=&0 \text{ vs. } H_1:\beta_j\neq 0\\
 H_0: \beta_{1}&=&\beta_{2}=\beta_{3}=0 \text{ vs. } H_1:\text{ at least one of these }\neq 0\\
 H_0: \beta_{1}&=&\beta_{2}=\cdots=\beta_{k}=0 \text{ vs. } H_1:\text{ at least one of these }\neq 0\\
 \end{eqnarray*}

---

We call the restricted model (when the null hypotesis is true) model B, or the smaller model.

These null hypotheses and alternative hypotheses can all be rewritten as a linear hypothesis 
\[H_0: {\bf C}{\bf \boldsymbol{\beta}}={\bf d} \text{ vs. } {\bf C}{\bf \boldsymbol{\beta}}\neq {\bf d} \]
by specifying ${\bf C}$ to be a $r \times p$ matrix and ${\bf d}$ to be a column vector of length $d$.

---

### The Wald test

The Wald test statistic is given as:

$$
w=({\bf C}\hat{ {\bf \boldsymbol{\beta}}}-{\bf d})^{\text T}[{\bf C}F^{-1}(\hat{\boldsymbol{\beta}}){\bf C}^{\text T}]^{-1}({\bf C}\hat{ {\bf \boldsymbol{\beta}}}-{\bf d})
$$
and measures the distance between the estimate ${\bf C}\hat{\boldsymbol{\beta}}$ and the value under then null hypothesis ${\bf d}$, weighted by the asymptotic covariance matrix of ${\bf C}\hat{\boldsymbol{\beta}}$.


Remember: $\text{Cov}({\bf C}\hat{\boldsymbol{\beta}})={\bf C}F^{-1}(\hat{\boldsymbol{\beta}}){\bf C}^{\text T}$.

Asymptotically, under the null hypothesis $w \sim \chi_r^2$ distribution with (where $r$ is the number of hypotheses tested).

$P$-values are calculated in the upper tail of the $\chi^2$-distribution.

Observe: to perform the test you only need to fit the larger model (and not the smaller).

---

For the special case that we only test one regression parameter, for example $\beta_k$:
$$ H_0: \beta_k=0 \text{ vs. } H_1: \beta_k\neq 0.$$
Now ${\bf C}\hat{\boldsymbol{\beta}}=\beta_k$ and ${\bf C}[F(\hat{\boldsymbol{\beta}})]^{-1}{\bf C}^{\text T}={\bf C}{\bf A}(\boldsymbol{\beta}){\bf C}^{\text T}=a_{kk}(\boldsymbol{\beta})$, 
and the Wald test becomes 
$$ (\hat{\beta}_k-\beta_k)[a_{kk}(\hat{\boldsymbol{\beta}})]^{-1}(\hat{\beta}_k-\beta_k)=\left(\frac{\hat{\beta}_k-\beta_k}{\sqrt{a_{kk}(\hat{\boldsymbol{\beta}})}}\right)^2=Z_k^2$$
so, asymptotically the square of the standard normal, which we know follows a $\chi^2$-distribution with 1 degree of freedom.

---

**Q**: Explain what you find in the columns named `z value` and `Pr(>|z|)` below, and which hypothesis tests these are related to. Are the hypothesis tests performed using the Wald test?


```{r, echo=TRUE}
library(investr)
fitgrouped=glm(cbind(y,n-y) ~ ldose, family="binomial", 
               data = investr::beetle) 
knitr::kable(summary(fitgrouped)$coefficients, digits=2)
```


<!-- ```{r} -->
<!-- library(statmod) -->
<!-- small.glm=glm(cbind(y, n-y) ~ 1, family = "binomial", data = investr::beetle)  -->
<!-- summary(small.glm) -->
<!-- z=glm.scoretest(small.glm, x2=investr::beetle$ldose) -->
<!-- z -->
<!-- 2*pnorm(abs(z),lower.tail=FALSE) -->
<!-- ``` -->

---


### The likelihood ratio test

The likelihood ratio test (LRT), which compares the likelihood of two models.

* First we maximize the likelihood for model A (the larger model) to get $L(\hat{\boldsymbol{\beta}}_A)$ and $\hat{\boldsymbol{\beta}}_A$.
* Then we maximize the likelihood for model B (the smaller model) to get $L(\hat{\boldsymbol{\beta}}_B)$ and $\hat{\boldsymbol{\beta}}_B$.

---

The likelihood of the larger model (A) will always be larger or equal to the likelihood of the smaller mode (B). (Why?)

The likelihood ratio statistic is defined as
$$- 2\ln \lambda=-2(\ln L(\hat{\boldsymbol{\beta}}_B)-\ln L(\hat{\boldsymbol{\beta}}_A)) $$ 
(so, $-2$ times small minus large).

---

Under weak regularity conditions the test statistic is approximately $\chi^2$-distributed with degrees of freedom equal the difference in the number of parameters in the large and the small model. This is general - and not related to the GLM! More in TMA4295 Statistical Inference!

$P$-values are calculated in the upper tail of the $\chi^2$-distribution.

Observe: to perform the test you need to fit both the small and the large model.

Notice: asymptotically the Wald and likelihood ratio test statistics have the same distribution, but the value of the test statistics might be different. How different? 

---

For the beetle data we compare model A=model with `ldose` as covariate with model B=model with only intercept. We use the loglikelihood-function that we made for the lecture session for week 2.

```{r}
library(investr)
fitgrouped=glm(cbind(y, n-y) ~ ldose, family = "binomial", data = investr::beetle) 
fitnull=glm(cbind(y, n-y) ~ 1, family = "binomial", data = investr::beetle) 

loglik <- function(par, args){
  y <- args$y; x <- args$x; n <- args$n
  res <- sum(y*x%*%par - n*log(1 + exp(x%*%par)))
  return(res)
}
```

---
```{r}
# call this with parameters estimated under model A=larger model
beetleargs = list(y = investr::beetle$y, 
            x = cbind(rep(1, nrow(investr::beetle)), investr::beetle$ldose), 
            n = investr:: beetle$n)

llA=loglik(matrix(fitgrouped$coefficients,ncol=1),args=beetleargs)

# then the smaller model, then we set the coeff for ldose to 0. B=smaller model
llB=loglik(matrix(c(fitnull$coefficients,0),ncol=1),args=beetleargs)
lrt=-2*(llB-llA)
lrt
pchisq(lrt,df=1, lower.tail = FALSE)
```

---

```{r}
anova(fitnull, fitgrouped, test="LRT")
```

**Q and A:** Here the small model is the model with only intercept and the large is the one with dose as covariate. This means that the null hypothesis is that "the small model is preferred" and our $p$-value is very small, so we reject the null hyptheses and stick with the model with dose as covariate. Observe that the LRT can be performed using `anova`.

---

# Deviance

The _deviance_ is used to assess model fit and also for model choice, and is based on the likelihood 
ratio test statistic. 

**Saturated model:**
One parmeter per group: estimate $\pi_j$ by the observed frequency for the group: $\tilde{\pi}_j=\frac{y_j}{n_j}$. Then $\tilde{\pi}$ is a $G$-dimensional column vector with the elements $\tilde{\pi}_j$.

This "imaginary model" is called the _saturated_ model.

**Candidate model:**
The model that we are investigated can be thought of as a _candidate_ model. Then we maximize the likelihood to get $\hat{\boldsymbol{\beta}}$ & thus $\hat{\pi}_j$. Then $\hat{\pi}$ is a $G$-dimensional column vector with the elements $\hat{\pi}_j$.

---

The _deviance_ is defined as the likelihood ratio statistic, where we put the saturated model in place of the larger model A and our candidate model in place of the smaller model B:

$$
\begin{aligned}
D&=-2(\ln L(\text{candidate model})-\ln L(\text{saturated model})) \\
&=-2(l(\hat{\pi})-l(\tilde{\pi}))=
-2\sum_{j=1}^G(l_j(\hat{\pi}_j)-l_j(\tilde{\pi}_j))
\end{aligned}
$$

For our logit model this can be written as (after some maths):
$$ 
D=2\sum_{j=1}^G [y_j\ln(\frac{y_j}{n_j\hat{\pi}_j})+(n_j-y_j)\ln(\frac{n_j-y_j}{n_j-n_j\hat{\pi}_j})]
$$

Verify this by yourself.

---

If our model is good, it should not be too far from the saturated model, and we measure this distance by the deviance.

If we want to investigate the null hypothesis that "our model fits the data well" to the negation, it is useful to know that asymptotically $D$ is distributed as $\chi^2$ with $G-p$ degrees of freedom (same reason as for the likelihood ratio test statistic).

This result depends on $n_j$ being large, hard to say how large (at least 5 is a rule of thumb).

---

The deviance is in `summary.glm` outputted as "Residual deviance", which we read off as `r summary(fitgrouped)$deviance`. Let's check for our beetle example by computing the formula for $D$ directly:

<!-- 
dev=function(yvec,nvec,pivec)
{ return(sum(yvec*log(yvec/(nvec*pivec))+(nvec-yvec)*log((nvec-yvec)/(nvec*(1-pivec)))))}
dev(investr::beetle$y,investr::beetle$n,fitgrouped$fitted.values)
-->

```{r}
D=deviance(fitgrouped)
D
G=dim(investr::beetle)[1]
G
p=2
1-pchisq(D,G-p)
```

So, do we have a good fit? 

---

The null hypothesis is that the candiate model is equally good as the saturated model. We do not reject this hypothesis at level 0.05. This means that we are satisfied with the candidate model. 

In the summary from glm also the socalled _NULL deviance_ is given. This is the deviance when the candicate model is the model with only intercept term present. This deviance asymptotically distributed as $\chi^2$ with $G-1$ degrees of freedom.  

---

\footnotesize

```{r}
summary(fitgrouped)
```

\normalsize

**Q**: where is the deviance(s) here and how do we use these?

---

## Analysis of deviance

In MLR we have seen that we may produce a sequential analysis of variance (Type I) by adding more and more terms to the model and comparing the scaled decrease in SSE by the scaled SSE of a full model.

For the binary regression we may adapt a similar strategy, but with using the scaled change in deviance instead of the SSE.

We use the infant respiratory disease data as an example

---

\footnotesize

```{r}
library(faraway)
fit=glm(cbind(disease, nondisease)~sex*food,family=binomial(link=logit),data=babyfood)
summary(fit)
anova(fit,test="Chisq")
```

\normalsize

**Q**: is it recommended (from the test) to add an interaction term to the model? What does it mean that the Residual deviance is 0 for the `sex*food` model?

---

## Deviance residuals
The deviance residuals are given by a signed version of each element in the sum for the deviance, that is
$$d_k=\text{sign}(y_k-n_k\hat{\pi}_k)\cdot \left\{ 
2[y_k\ln(\frac{y_k}{n_k\hat{\pi}_k})+(n_k-y_k)\ln(\frac{n_k-y_k}{n_k-n_k\hat{\pi}_k})]\right\}^{1/2}$$
where the term $\text{sign}(y_k-n_k\hat{\pi}_k)$ makes negative residuals possible. 


# Model assessment and choice

The fit of the model can be assessed based on goodness of fit statistics (and related tests) and by residual plots. Model choice can be made from analysis of deviance, or by comparing the AIC for different models.

## Deviance test for grouped data

We may use the deviance test presented before to test if the model under study is preferred compared to the saturated model.

---

## Pearson test and residuals
An alternative to the deviance test is the Pearson test. We will look in more detail at this test in a Module 4. The Pearson test statistic can be written as a function of the Pearson residuals, which for the binomial regression is given as:
$$ r_{j}=\frac{y_j-n_j\hat{\pi}_j}{\sqrt{n_j \hat{\pi}_j(1-\hat{\pi}_j)}}$$
Remark: A standardized version scales the Pearson residuals with $\sqrt{1-h_{kk}}$ similar to the standardized residuals for the normal model. Here $h_{kk}$ is the diagonal element number $k$ in the hat matrix ${\bf H}={\bf X}({\bf X}^T{\bf X})^{-1}{\bf X}^T$.

---

The Pearson $\chi^2$-goodness of fit statistic is given as
$$ X_P^2=\sum_{j=1}^G r_j^2=\sum_{j=1}^G  \frac{(y_j-n_j\hat{\pi}_j)^2}{n_j \hat{\pi}_j(1-\hat{\pi}_j)}$$
The Pearson $\chi^2$ statistic is asymptotically equivalent to the deviance statistic and thus is asymptotically $\chi^2_{G-p}$.

The Pearson $\chi^2$ statistic is not a good choice if any of the groups have a low expected number of observations, i.e. $n_j \hat{\pi}_j$ is small (below 1).

---

## Model assessment with continuous covariates

If data have continuous covariates it is possible to form groups based making intervals for continuous covariates.
Alternatively grouping on predicted probabilites can be done.

For continuous data the Hosmer Lemeshow test can be used - not on our reading list.

---

## Plotting residuals

Deviance and Pearson residuals can be used for checking the fit of the model, by plotting the residuals against fitted values and covariates. 

If $n_j$ is small for the covariate patterns the residual plots may be relatively uninformative.

Residual plots for the logistics regression - and for the GLM in general - is highly debated, and we will not put much emphasis on residual plots for this module.

---

```{r,echo=FALSE, fig.height=3}
library(ggplot2)
df=data.frame("fitted"=fitgrouped$fitted.values,"dres"=residuals(fitgrouped,type="deviance"),"ldose"=investr::beetle$ldose)
ggplot(df,aes(x=fitted,y=dres))+geom_point()
ggplot(df,aes(x=ldose,y=dres))+geom_point()
```

---

## Other plots

A useful plot is to show observed and fitted proportions (grouped data) plotted against the linear predictor or covariates.

```{r,echo=TRUE, fig.height=4}
library(ggplot2)
df=data.frame("fitted"=fitgrouped$fitted.values,"dres"=residuals(fitgrouped,type="deviance"),"ldose"=investr::beetle$ldose,"frac"=investr::beetle$y/investr::beetle$n)
ggplot(df,aes(x=ldose))+geom_point(aes(y=frac,colour="observed"))+geom_point(aes(y=fitted,colour="fitted"))
```

---

## AIC

It is known to us from multiple linear regression that if a model is chosen based on a goodness of fit statistic (like the SSE or $R^2$ in multiple linear regression) will in general result in us choosing a to big model (to many parameters fit). The Akaike informations criterion - that we studied for multiple linear regression - can also be used for binary regression:
Let $p$ be the number of regression parameters in our model.
$$\text{AIC} =-2 \cdot l(\hat{\boldsymbol{\beta}})+2p$$
A scaled version of AIC, standardizing for sample size, is sometimes preferred.

To use AIC for model selection you use the model with the _smallest_ AIC.

We may also use the BIC, where $2p$ is replaced by $\log(G)\cdot p$ or $\log(n)\cdot p$.

---

\footnotesize

```{r}
library(faraway)
fit1=glm(cbind(disease, nondisease)~1,family=binomial(link=logit),data=babyfood)
fit2=glm(cbind(disease, nondisease)~sex,family=binomial(link=logit),data=babyfood)
fit3=glm(cbind(disease, nondisease)~food,family=binomial(link=logit),data=babyfood)
fit4=glm(cbind(disease, nondisease)~food+sex,family=binomial(link=logit),data=babyfood)
fit5=glm(cbind(disease, nondisease)~food*sex,family=binomial(link=logit),data=babyfood)
AIC(fit1,fit2,fit3,fit4,fit5)
```

\normalsize

**Q**: Which of these 5 models would you prefer?


# Overdispersion and estimating overdispersion parameter
 
When we have grouped data: $Y_j \sim \text{Bin} (n_j, \pi_j)$ and Var$(Y_j) = n_j\pi_j(1-\pi_j)$.

It is possible to estimate the variance (within a group) by $n_j\bar{y_j}(1-\bar{y_j})$ where $\bar{y_j} = y_j/n_j$ (this is an estimate of $\pi_j$ for group $j$). We call this the *empirical variance*.

In a logistic regresson we estimate $\hat{\pi}_j = h({\bf x}_j^T\hat{\boldsymbol{\beta}})$ ($h(\cdot)$ is the inverse link function) which is 

$$\hat{\pi_j} = \frac{\exp(x_j^T \hat{\boldsymbol{\beta}})}{1+\exp(x_j^T \hat{\boldsymbol{\beta}})} $$

for a logistic regression. This would give the *estimated binomial variance* for $Y_j$ as $n_j\hat{\pi_j}(1-\hat{\pi_j})$.

---

Some times the empirical variance is much larger than the estimated binomial variance of the model. This is called *overdispersion* and may occur when the individual responses within the groups are correlated, or when the model could be improved upon (missing/unobserved covariates?).

Positively correlated binary variables will give a variance of the sum that is larger than for uncorrelated variables, e.g.

$$\text{Var}(\sum_{k=1}^K Y_k) = \sum_{k=1}^K\text{Var}(Y_k) + 2\sum_{k<l} \text{Cov}(Y_k, Y_l).$$

---

This can be handeled by including an *overdispersion parameter*, named $\phi$, in the variance formula:

$$ \text{Var}(Y_j) = \phi n_j \pi_j (1-\pi_j)$$

---

The overdispersion parameter can be estimated as the average Pearson statistic or average deviance

$$\hat{\phi}_D = \frac{1}{G-p} D$$

where $D$ is the deviance. Note that similarity to $\hat{\sigma^2} = 1/(n-p)\cdot\text{SSE}$ in the MLR. The Cov$(\hat{\boldsymbol{\beta}})$ can then be changed to $\hat{\phi}F^{-1}(\hat{\boldsymbol{\beta}})$.  

Remark: We are now moving from likelihood to quasi-likelihood theory, where only E$(Y_j)$ and Var$(Y_j)$ - and not the distribution of $Y_j$ - are used in the estimation. 

In Modules 7 and 8 we will look at using multilevel models to handle correlated observations.

---

```{r}
library(investr)
estpi=investr::beetle$y/investr::beetle$n
empvars=investr::beetle$n*estpi*(1-estpi)
fit=glm(cbind(y, n-y) ~ ldose, family = "binomial", data = investr::beetle)
modelestvar=investr::beetle$n*fit$fitted.values*(1-fit$fitted.values)
cbind(empvars,modelestvar)
est.dispersion=fit$deviance/fit$df.residual 
est.dispersion
summary(fit,dispersion=est.dispersion,correlation=TRUE)
fitquasi=glm(cbind(y, n-y) ~ ldose, family = "quasibinomial", data = investr::beetle)
# preferred method of estimation is to use quasilikelihood
summary(fitquasi)$dispersion
```

---


# References for further reading

* A. Agresti (2015): "Foundations of Linear and Generalized Linear Models." Wiley.
* A. J. Dobson and A. G. Barnett (2008): "An Introduction to Generalized Linear Models", Third edition. 
* J. Faraway (2015): "Extending the Linear Model with R", Second Edition. <http://www.maths.bath.ac.uk/~jjf23/ELM/>
* P. McCullagh and J. A. Nelder (1989): "Generalized Linear Models". Second edition.


---


If we have time

---

## Look back at MLR - what is $s(\boldsymbol{\beta})$ and $F(\boldsymbol{\beta})$ then?

1. $Y_i \sim \text{N}(\mu_i, \sigma^2)$

2. $\eta_j = x_i^T\boldsymbol{\beta}$

3. $\mu_i = \eta_i$ (identity response function and link function)

Likelihood:

$$L(\boldsymbol{\beta}) = \left(\frac{1}{2\pi}\right)^{n/2}\left(\frac{1}{\sigma^2}\right)^{n/2} \exp\left(-\frac{1}{2\sigma^2}(y-{\bf X}\boldsymbol{\beta})^T(y-{\bf X}\boldsymbol{\beta})\right)$$

 
Loglikelihood:
$$l(\boldsymbol{\beta}) = \ln L(\boldsymbol{\beta}) = -\frac{n}{2} \ln (2\pi) - \frac{n}{2} \ln (\sigma^2) - \frac{1}{2\sigma^2}(y-{\bf X}\boldsymbol{\beta})^T(y-{\bf X}\boldsymbol{\beta})$$

---


Since $(y-{\bf X}\boldsymbol{\beta})^T(y-{\bf X}\boldsymbol{\beta}) = Y^TY - 2Y^T{\bf X}\boldsymbol{\beta}  + \boldsymbol{\beta}^T{\bf X}^T{\bf X}\boldsymbol{\beta}$, then


$$s(\boldsymbol{\beta}) = \frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = -\frac{1}{2\sigma^2}(2{\bf X}^T{\bf X}\boldsymbol{\beta}-2{\bf X}^TY) = \frac{1}{\sigma^2}({\bf X}^TY - {\bf X}^T{\bf X}\boldsymbol{\beta})$$

and $s(\hat{\boldsymbol{\beta}}) = 0$ gives ${\bf X}^TY - {\bf X}^T{\bf X}\boldsymbol{\beta} = 0$ which can be solved on closed form giving $\hat{\boldsymbol{\beta}} = ({\bf X}^T{\bf X})^{-1}{\bf X}^TY$. So, no need for iterative methods.

---


Finally, observed Fisher information matrix.

$$H(\boldsymbol{\beta}) = \frac{\partial s(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^T} = -\frac{\partial}{\partial \boldsymbol{\beta}_T}(\frac{1}{\sigma^2}{\bf X}^TY - \frac{1}{\sigma^2}{\bf X}^T{\bf X}\boldsymbol{\beta}) = \frac{1}{\sigma^2}{\bf X}^T{\bf X}$$

which is independent on $\boldsymbol{\beta}$, and also we see that $F(\boldsymbol{\beta})=\text{E}(H(\boldsymbol{\beta}))=H(\boldsymbol{\beta})$ since no random variables are present. The identity link is also the canonical link. Finally, the (asymptotic) covariance of the ML estimate is $F^{-1}(\hat{\boldsymbol{\beta}})=({\bf X}^T{\bf X})^{-1}\sigma^2$ which we know as $\text{Cov}(\hat{\boldsymbol{\beta}})$.

---


# Exponential family - and canonical link

In Module 1 we introduced distributions of the $Y_i$, that could be written in the form of a _univariate exponential family_
\[ f(y_i\mid \theta_i)=\exp \left( \frac{y_i \theta_i-b(\theta_i)}{\phi}\cdot w_i + c(y_i, \phi, w_i) \right) \]
where

* $\theta_i$ is called the canonical parameter and is a parameter of interest

* $\phi$ is called a nuisance parameter (and is not of interest to us=therefore a nuisance (plage))

* $w_i$ is a weight function, in most cases $w_i=1$

* $b$ and $c$ are known functions.

It can be shown that $\text{E}(Y_i)=b'(\theta_i)$ and $\text{Var}(Y_i)=b''(\theta_i)\cdot \frac{\phi}{w_i}$.

---

In Module 1 we found that the binomial distribution $Y_i\sim \text{bin}(1,\pi_i)$ is an exponential family (derivation from Module 1: <https://www.math.ntnu.no/emner/TMA4315/2017h/Module1ExponentialFamily.pdf>)

and that

* $\theta_i=\ln( \frac{\pi_i}{1-\pi_i})$ is the canonical parameter
* $\phi=1$, no nuisance
* $w_i=1$
* $b(\theta_i)=\ln(1+\exp(\theta_i))$

---

Recall that in a GLM we choose a link function $g$, linking the linear predictor and the mean: $\eta_i=g(\mu_i)$. For the logit model we had that $\eta_i=\ln(\frac{\pi_i}{1-\pi_i})$.

Now (new to us) - every exponential family has a unique _canonical link function_ such that
$$\theta_i=\eta_i$$
Since $\eta_i=g(\mu_i)$ this means to us that we need 
$$ g(\mu_i)=\theta_i$$
to have a canonical link.

**Q:** Is the logit link the canonical link for the binary model?

```{r , echo=FALSE, eval=FALSE, warning=FALSE, message=FALSE}
cat("Answer \n")
cat("yes, since $\theta_i=\ln( \frac{\pi_i}{1-\pi_i})=g(\pi_i)$ then the logit link is the canonical link for the binary regression.\n")
```


**A**:

Yes, since $\theta_i=\ln( \frac{\pi_i}{1-\pi_i})=g(\pi_i)$ then the logit link is the canonical link for the binary regression.

---

##Properties of a GLM with canonical link

1. The log-likelihod is always concave so that the ML estimated is always unique (given that it exists).

2. The observed Fisher information matrix $H(\boldsymbol{\beta})$ _equals_ the expected Fisher information matrix $F(\boldsymbol{\beta})$.
That is, $$-\frac{\partial^2 l}{\partial \boldsymbol{\beta} \boldsymbol{\beta}^T}=\text{E}(-\frac{\partial^2 l}{\partial \boldsymbol{\beta} \boldsymbol{\beta}^T})$$

Proving this is beyond the scope of this course.

---

# Parameter estimation - in practise
To find the ML estimate $\hat{\boldsymbol{\beta}}$ we need to solve
$$s(\hat{\boldsymbol{\beta}})=0$$
We have that the score function for the logit model is:

$$
s(\boldsymbol{\beta})=\sum_{j=1}^G {\bf x}_j (y_j-n_j\pi_j)
$$

<!---
where $\pi_j=\frac{\exp({\bf x}_j^T\hat{\boldsymbol{\beta}})}{1+\exp({\bf x}_j^T\hat{\boldsymbol{\beta}})}$. Observe that this is a non-linear function in $\boldsymbol{\beta}$, and has no closed form solution. 
--->

---
