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

# Overview

## Learning material

* Textbook: Fahrmeir et al (2013): Chapter 5.2, 5.3.
* [Classnotes 27.09.2018](https://www.math.ntnu.no/emner/TMA4315/2018h/TMA4315M4H20180927.pdf) 
* [Classnotes 04.10.2018](https://www.math.ntnu.no/emner/TMA4315/2018h/TMA4315M4H20181004.pdf) 

---

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

* examples of count data
* the Poisson distribution
* regression with count data
* Poisson regression with log-link
* parameter estimation (ML): log-likelhood, score vector, information matrix to give iterative calculations 
* asymptotic MLE properties
* confidence intervals and hypothesis tests (Wald, score and LRT)

Jump to [interactive (week 1)](#interactivew1)

---

### [Second week](#secondweek)

* Count data with Poisson regression (continued)
* deviance, model fit and model choice
* overdispersion
* rate models and offset

* Modelling continuous response data: lognormal and gamma
* the gamma distribution
* the gamma GLM model
* gamma likelihood and derivations thereof
* dispersion parameter: scaled and unscaled deviance

---

**<a id="secondweek"> SECOND WEEK</a>**

# Poisson regression for count data

```{r, echo=FALSE}
library(ggplot2)
library(GGally)
crab=read.table("https://www.math.ntnu.no/emner/TMA4315/2018h/crab.txt")
colnames(crab)=c("Obs","C","S","W","Wt","Sa")
crab=crab[,-1] #remove column with Obs
crab$C=as.factor(crab$C)
crab$S=as.factor(crab$S)

```

What did we do last week?

- Examples 
— GLM model 
— loglikelihood, score function and Fisher information matrix 
— asympototic results for $\hat{\boldsymbol{\beta}}$ and Wald, score and LRT.

# Deviance

The deviance is the difference in model fit between 2 models. Three models are helpful:

- **Null model**: model with only an intercept
- **Candidate model:** The model that we are investigating. We maximize the likelihood and get $\hat{\beta}$
- **Saturated model:** A model with a parameter for every $\lambda_i$ by the observed count for observation $i$. 

Notation:

$\hat{\lambda}_i = \hat{y}_i$: prediction for the candidate model for observation $i$
$\tilde{\lambda}_i=y_i$: prediction for the saturated model for observation $i$

---

# Model assessment and model choice: using the deviance

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

The (log-) likelihood ratio between two models is

$$
\begin{aligned}
l(\hat{\lambda})-l(\tilde{\lambda})) &=\sum_{i=1}^n [y_i \ln(\hat{\lambda}_i)-\hat{\lambda}_i-\ln(y_i!)] - \sum_{i=1}^n [y_i \ln(\tilde{\lambda}_i)-\tilde{\lambda}_i-\ln(y_i!)] \\
&= \sum_{i=1}^n [y_i \left(\ln(\hat{\lambda})_i - \ln(\tilde{\lambda}_i) \right)-\hat{\lambda}_i-\tilde{\lambda}_i] \\
&= \sum_{i=1}^n [y_i \ln(\hat{\lambda}_i/\tilde{\lambda}_i)- (\hat{\lambda}_i-\tilde{\lambda}_i)] \\
\end{aligned}
$$

($\lambda_i=E(y_i)$ for the model of interest)

The Saturated Deviance is 

$$
D = -2 l(\lambda)=-2\sum_{i=1}^n [y_i \ln(y_i)-y_i-\ln(y!)] 
$$

because $\hat{\lambda}_i = y_i$ for this model

---

# Deviance test

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

---

## Likelihood Ratio Tests with Deviance

The likelihood ratio test can be performed using the difference between two deviances:

$$
\begin{aligned}
LRT &= D_1 - D_2 =-2(l_1(\hat{\lambda}_1)-l(\tilde{\lambda})) -
(-2(l_2(\hat{\lambda}_2)-l(\tilde{\lambda}))) \\
&=-2(l_1(\hat{\lambda}_1)-l_2(\tilde{\lambda}_2))
\end{aligned}
$$

This follows a $\chi^2_p$ distribution with $k$ equal to the difference in number of parameters

---

## Deviance test

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


<!-- http://thestatsgeek.com/2014/04/26/deviance-goodness-of-fit-test-for-poisson-regression/
We are thus not guaranteed, even when the sample size is large, that the test will be valid (have the correct type 1 error rate). One of the few places to mention this issue is Venables and Ripley's book, Modern Applied Statistics with S. Venables and Ripley state that one situation where the chi-squared approximation may be ok is when the individual observations are close to being normally distributed and the link is close to being linear. Pawitan states in his book In All Likelihood that the deviance goodness of fit test is ok for Poisson data provided that the means are not too small. When the mean is large, a Poisson distribution is close to being normal, and the log link is approximately linear, which I presume is why Pawitan's statement is true (if anyone can shed light on this, please do so in a comment!).
-->
$$ D=2 \sum_{i=1}^n [y_i\ln(\frac{y_i}{\hat{\lambda}_i})-(y_i-\hat{\lambda}_i)]$$

Remark: if $\sum_{i=1}^n y_i=\sum_{i=1}^n \hat{\lambda}_i$ then deviance will be equal to 
$$ D=2 \sum_{i=1}^n y_i\ln(\frac{y_i}{\hat{\lambda}_i})$$

The deviance statistic approximately follows a $\chi^2_{n-p}$, at least when the counts are not low. 


---

## Pearson test

The Pearson $\chi^2$-goodness of fit statistic is given as the sum of the squared Pearson residuals

$$ X_P^2=\sum_{i=1}^n r_i^2=\sum_{i=1}^n  \frac{(y_j-\hat{\lambda}_i)^2}{\hat{\lambda}_i}$$
where $\hat{\lambda}_i=\exp({\bf x}_i^T\hat{\beta})$.
The Pearson $\chi^2$ statistic is asymptotically equivalent to the deviance statistic and thus is asymptotically $\chi^2_{n-p}$ (proof: do a Taylor series expansion of the deviance).

---

## Remarks

The asymptotic distribution of both statistics (deviance and Pearson) are questionable when there are many low counts. Agresti (1996, page 990) suggest analysing grouped data, for example by grouping by width in the horseshoe crab example.

The Pearson statistic is also used for testing independence in contingency tables - we will do that in Compulsory Exercise 2.

---

## Example: goodness of fit with female horseshoe crabs

Comment on the analysis. Is this a good fit? What might a bad fit be due to? 

```{r}
model3=glm(Sa~W+C,family=poisson(link=log),data=crab,contrasts=list(C="contr.sum"))
# summary(model3)
1-pchisq(model3$deviance,model3$df.residual)
Xp=sum(residuals(model3,type="pearson")^2)
Xp
1-pchisq(Xp,model3$df.residual)
```


---

## AIC

Identical to Module 3 - we may use the Akaike informations criterion. Let $p$ be the number of regression parameters in our model.
$$\text{AIC} =-2 \cdot l(\hat{\beta})+2p$$
A scaled version of AIC, standardizing for sample size, is sometimes preferred. And, we may also use the BIC, where $2p$ is replaced by $\log(n)\cdot p$.

---

## Analysis of deviance
Identical to Module 3 we may also sequentially compare models, and use analysis of deviance for this.

---

# Residuals

Two types of residuals are popular: _deviance_ and _Pearson_. These are based on the deviance and $\chi^2$ statistics: basically they are the contribution of each data point to that statistic. 

But the sign has to be included in the deviance residual


---

## Deviance residuals


The deviance residuals are given by a signed version of each element in the sum for the deviance, that is

$$
d_i=\text{sign}(y_i-\hat{\lambda}_i)\cdot \left\{ 
2[y_i\ln(\frac{y_i}{\hat{\lambda}_i})-(y_i-\hat{\lambda}_i)]\right\}^{1/2}
$$

where the term $\text{sign}(y_i-\hat{\lambda}_i)$ makes negative residuals possible - and we get the same sign as the _Pearson residuals_

---

## Pearson residuals

The Pearson residuals are given as 
$$ r_i=\frac{o_i-e_i}{\sqrt{e_i}}$$
where $o_i$ is the observed count for observation $i$ and $e_i$ is the estimated expected count for observation $i$. We have 
that $o_i=y_i$ and $e_i==\hat{\lambda}_i=\exp({\bf x}_i^T\hat{\beta})$.

Remark: A standardized version scales the Pearson residuals with $\sqrt{1-h_{ii}}$ similar to the standardized residuals for the normal model. Here $h_{ii}$ is the diagonal element number $i$ in the hat matrix ${\bf H}={\bf X}({\bf X}^T{\bf X})^{-1}{\bf X}^T$.

---


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

Normality of residuals is not assumed, but for large counts can be reasonable, and can be checked using qq-plots as for the MLR in Module 2.

Below - notice the trend in the residuals, this is due to the discrete nature of the response. The plot with different shades of blue shows that the structures are for equal values of $y$.

---

# R Code

```{r}
model3=glm(Sa~W+C,family=poisson(link=log),
           data=crab,
           contrasts=list(C="contr.sum",S="contr.sum"))
df=data.frame("Sa"=crab$Sa,"fitted"=model3$fitted.values,
              "dres"=residuals(model3,type="deviance"),
              "pres"=residuals(model3,type="pearson"))

library(ggplot2)
# create the plot
gg1=ggplot(df)+geom_point(aes(x=fitted,y=dres,color="deviance"))+
     geom_point(aes(x=fitted,y=pres,color="pearson"))
```

---


```{r,echo=TRUE}
gg1   
```

---

## Pearson Residuals

```{r,echo=showsol}
gg2=ggplot(df)+geom_point(aes(x=fitted,y=pres,color=Sa))
gg2

```

---

## Normal Prabability Plots


```{r,echo=TRUE, fig.height=5}
dff=data.frame("devres"=residuals(model3,type="deviance"),"pearsonres"=residuals(model3,type="pearson"))
ggplot(dff, aes(sample = devres)) +
  stat_qq(pch = 19) +
  geom_abline(intercept = 0, slope = 1, linetype = "dotted") +
  labs(x = "Theoretical quantiles", y = "Deviance residuals", title = "Normal Q-Q", subtitle = deparse(model3$call))
```

---

```{r,echo=TRUE, fig.height=5}
ggplot(dff, aes(sample = pearsonres)) +
  stat_qq(pch = 19) +
  geom_abline(intercept = 0, slope = 1, linetype = "dotted") +
  labs(x = "Theoretical quantiles", y = "Pearson residuals", title = "Normal Q-Q", subtitle = deparse(model3$call))
```

<!-- 
https://theoreticalecology.wordpress.com/2016/08/28/dharma-an-r-package-for-residual-diagnostics-of-glmms/

Diagnostics in Poisson regression: <https://www.youtube.com/watch?v=T-UPWa0RVM0>
snakker om hvor mange utenfor 1.96?-->
<!--
library(faraway)
modl <- glm(doctorco ~ sex + age + agesq + income + levyplus + freepoor + 
            freerepa + illness + actdays + hscore + chcond1 + chcond2,
            family=poisson, data=dvisits)
plot(modl)
This is the appearance you expect of such a plot when the dependent variable is discrete.

Each curvilinear trace of points on the plot corresponds to a fixed value kk of the dependent variable yy. Every case where y=ky=k has a prediction ŷ y^; its residual--by definition--equals k−ŷ k−y^. The plot of k−ŷ k−y^ versus ŷ y^ is obviously a line with slope −1−1. In Poisson regression, the x-axis is shown on a log scale: it is log(ŷ )loga(y^). The curves now bend down exponentially. As kk varies, these curves rise by integral amounts. Exponentiating them gives a set of quasi-parallel curves. (To prove this, the plot will be explicitly constructed below, separately coloring the points by the values of yy.)

We can reproduce the plot in question quite closely by means of a similar but arbitrary model (using small random coefficients):

# Create random data for a random model.
set.seed(17)
n <- 2^12                       # Number of cases
k <- 12                         # Number of variables
beta = rnorm(k, sd=0.2)         # Model coefficients
x <- matrix(rnorm(n*k), ncol=k) # Independent values
y <- rpois(n, lambda=exp(-0.5 + x %*% beta + 0.1*rnorm(n)))

# Wrap the data into a data frame, create a formula, and run the model.
df <- data.frame(cbind(y,x))    
s.formula <- apply(matrix(1:k, nrow=1), 1, function(i) paste("V", i+1, sep=""))
s.formula <- paste("y ~", paste(s.formula, collapse="+"))
modl <- glm(as.formula(s.formula), family=poisson, data=df)

# Construct a residual vs. prediction plot.
b <- coefficients(modl)
y.hat <- x %*% b[-1] + b[1]     # *Logs* of the predicted values
y.res <- y - exp(y.hat)         # Residuals
colors <- 1:(max(y)+1)          # One color for each possible value of y
plot(y.hat, y.res, col=colors[y+1], main="Residuals v. Fitted")
-->

---

# Overdispersion

Count data might show greater variability in the response counts than we would expect if the response followed a Poisson distribution. This is called _overdispersion_. 

Example: newspaper sales with tourist bus.

Our model states that the variance $\text{Var}(Y_i)=\lambda_i$. If we change the model to $\text{Var}(Y_i)=\phi\lambda_i$ we may allow for an increased variance due to heterogeneity among subjects.

Or, we can miss several covariate, then any data point is a mixture of several Poisson populations, each with its own mean for the response. 

This heterogeniety may give an overall response distribution where the variance is greater than the standard Poisson variance.
 
---

The overdispersion parameter can be estimated from the Pearson statistic or deviance

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

where $D$ is the deviance. Note that similarity to $\hat{\sigma^2} = 1/(n-p)\cdot\text{SSE}$ in the MLR.

Cov$(\hat{\beta})$ can then be changed to $\hat{\phi}F^{-1}(\hat{\beta})$, so we multiply the standard error by the square root of $\hat{\phi}_D$.

---

# Estimating Overdispersion

(more on quasipoisson later in the course)

```{r}
m.od=glm(Sa~W, family=poisson(link=log), 
             data=crab)
m.qp=glm(Sa~W, family=quasipoisson(link=log), data=crab)
# summary.glm(m.od)
summary.glm(m.qp)$dispersion

(OverDisp.Dev <- m.od$deviance/m.od$df.residual)
Chi2 <- sum((crab$Sa-fitted(m.od))^2/fitted(m.od))
(OverDisp.Chi2 <- Chi2/m.od$df.residual)

```

---

# Rate models with Counts

In the Poisson process we might analyse an event that occurs within a time interval or region in space, and therefore it is often of interest to model the _rate_ at which events occur.

Examples:

* crime rates in cities
* death rate for smokers vs. non-smokers
* rate of auto thefts in cities 

We model the rates by using an *offset* to convert to counts.

---


We don't want a model for $Y_i$ but for $Y_i/t_i$:

* Let $t_i$ denote the index (population size in the example) associated with observation $i$. 
* We still assume that $Y_i$ follows a Poisson distribution, but we now include the index in the modelling and focus on $Y_i/t_i$.
* The expected value of $Y_i/t_i$ would then be $\text{E}(Y_i)/t_i=\lambda_i/t_i$.

A log-linear model would be
$$ \log(\lambda_i/t_i)={\bf x}_i^T \beta$$
We may equivalently write the model as

$$
\log(\lambda_i)-\log(t_i)={\bf x}_i^T \beta
$$
This adjustment term is called an _offset_ and is a known quantity. Equivalently we have $\log(\lambda_i)={\bf x}_i^T \beta + \log(t_i)$


The expected number of outcomes will then satisfy
$$ \text{E}(Y_i)=\lambda_i=t_i \exp({\bf x}_i^T \beta).$$

---

## Example: British doctors and rate models

British doctors  sent a questionnaire (in 1951) about whether they smoked tobacco, and later information about their deaths were collected. 

Research questions: 

1) Is the death rate higher for smokers than for non-smokers? 
2) If so, by how much? 
3) How is this related to age?

```{r}
library(boot)
data(breslow)
#n=person-year, ns=smoker-years, 
#y=number of deaths due to cad, 
breslow$age <- factor(breslow$age) #age=midpoint 10 year age group, 
breslow$smoke <- factor(breslow$smoke) # smoke=smoking status
```

---

# Writing an offset

Here our count depends on $n$ 

(we are actually using a Poisson approximation to the binomial)

There are 2 ways of coding an offset in `R`:

```{r}
# first age + smoke (but no interaction)
fit1 <- glm(y~age+smoke,offset=log(n),family=poisson, data=breslow)
fit1a<- glm(y~age+smoke + offset(log(n)),family=poisson, data=breslow)
```

--- 

# Do we need an interaction?

```{r}
# do we need interaction? 
fit2<- update(fit1,. ~. +smoke*age)
anova(fit1,fit2,test="Chisq")

```

---

# The final model

Number of deaths per 1000 doctors. 

```{r}
# year 40 nonsmokers should only be the intercept
exp(fit2$coefficients[1])*1000
# 80 year olds who smoke
1000*exp(sum(fit2$coefficients[c(1,5,6,10)]))
```

(you can also use `predict()` to do this)

# Modelling continuous positive response data

## Examples of continuous positive responses

* Insurance: Claim sizes
* Medicine: Time to blood coagulation (main example)
* Biology: Time in various development stages for fruit fly
* Meteorology: Amount of precipitation (interactive session - exam question 2012)

This is also covered in **survival analysis**, but that often sidesteps modelling the actual distributions

---

## Models for continuous positive responses

* Lognormal distribution on response
* Gamma distribution on response
* Inverse Gaussian distribution on response (we will not consider this here)

---

## Time to blood coagulation

The data is clotting time of blood (in seconds) `y` for normal plasma diluted to nine different percentage concentrations `u` with prothrombin-free plasma (whatever that is!). 

To induce the clotting a chemical called thromboplasting was used, and in the experiment two different lots of the chemical were used - denoted `lot`. Our aim is to investigate the relationship between the clotting time and the dilution percentage, and look at differences between the lots.

---

```{r}
clot = read.table("https://www.math.ntnu.no/emner/TMA4315/2018h/clot.txt",header=T)
clot$lot = as.factor(clot$lot)
summary(clot)
```

# Lognormal distribution 

Let $Y_i$ be the response on the original scale, where $Y_i>0$. 

Transform the response to a logaritmic scale: $Y^*_i=\ln(Y_i)$. Then, assume that transformed responses follow a normal distribution (or follows approximately) and use ordinary MLR. This means we have a GLM with normal response and identity link (on logarithmic scale of reponse).

1. $Y^*_i \sim N(\mu^{*}_i,\sigma^{*2})$
2. $\eta_i={\bf x}_i^T\beta$
3. $\mu^{*}_i=\eta_i$ (identity link)

---

There are two ways of looking at this, 

1. either this is just a transformation to achieve approximate normality, or 
2. we assume that the original data follows a lognormal distribution. 

If we assume that the original data really comes from a lognormal distribution, then it can be shown that 

$$\text{E}(Y_i)=\exp(\mu^{*}_i) \cdot \exp(\sigma^{*2}/2)$$
$$\text{Var}(Y_i) =\exp(\sigma^{*2} -1)\cdot \mu_i^2$$

i.e. standard deviation proportional to expectation.

---

# Gamma regression

## The gamma distribution
We have seen that a gamma distributed variable may be the result of the time between events in a Poisson process. The well known $\chi^2_{\delta}$-distribution is a special case of the gamma distribution ($\frac{\nu}{\mu_i}=2$, $\nu=\frac{\delta}{2}$).

There are many parameterization for the gamma distribution, but we will stick with the one used in our textbook (page 643):

$Y_i \sim Ga(\mu_i,\nu)$ with density
$$ f(y_i)=\frac{1}{\Gamma(\nu)} (\frac{\nu}{\mu_i})^{\nu} y_i^{\nu-1}\exp(-\frac{\nu}{\mu_i}y_i) \text{ for }y_i>0$$

---

# Comparing the lognormal and gamma

```{r}
orgmu=1; orgsd=0.3 # normal mean and sd
mu= exp(orgmu +  orgsd^2/2) # = shape*scale
scale = (exp(orgsd^2)-1)*mu
shape = mu/scale
library(ggplot2)
xrange=range(0,10)
```

---

# These should have the same mean and variance

```{r, fig.height=5}
ggplot(data.frame(x=xrange),aes(xrange))+
   xlab(expression(x))+ 
    stat_function(fun=dlnorm, args=list(meanlog=orgmu,sdlog=orgsd), geom="line", colour="red", n=1001)+
    stat_function(fun=dgamma, args=list(shape=shape,scale=scale), geom="line", colour="blue")
```

---

We found in Module 1 that the gamma distribution is an exponential family, with

* $\theta_i=-\frac{1}{\mu_i}$ is the canonical parameter
* $\phi=\frac{1}{\nu}$,
* $w_i=1$
* $b(\theta_i)=-\ln(-\theta_i)$
* $\text{E}(Y_i)=b'(\theta_i)=-\frac{1}{\theta_i}=\mu_i$
* $\text{Var}(Y_i)=b''(\theta_i)\frac{\psi}{w_i}=\frac{\mu_i^2}{\nu}$

(if you don't remember, work it out!)

---

For a GLM model we have canonical link if 
$$\theta_i=\eta_i$$
Since $\eta_i=g(\mu_i)$ this means to us that we need 
$$\theta_i=g(\mu_i)=-\frac{1}{\mu_i}$$
saying that with the canonical link is $-\frac{1}{\mu_i}$. 

However, the most commonly used link is $g(\mu_i)=\ln(\mu_i)$, and the identity link is also used.

**Q:** Discuss the implications on $\eta_i$ when using the canonical link. Why might the log-link be preferred?

---

Remark: often the inverse and not the negative inverse is used, and since
$$g(\mu_i)=-\frac{1}{\mu_i}={\bf x}_i^T\beta$$
then $$\frac{1}{\mu_i}=-{\bf x}_i^T\beta={\bf x}_i^T\beta^*$$
where $\beta^*=-\beta$.

---

## Gamma GLM model
1. $Y_i \sim Ga(\mu_i,\nu)$
2. $\eta_i={\bf x}_i^T\beta$

3. Popular link functions:
   + $\eta_i=\mu_i$ (identity link)
   + $\eta_i=\frac{1}{\mu_i}$ (inverse link)
   + $\eta_i=\ln(\mu_i)$ (log-link)

**Remark:** In our model the parameter $\mu_i$ varies with $i$ but $\nu$ is the same for all observations.

---

### Example: Time to blood coagulation
A simple model to start with is as follows (dosages often analysed on log scale):

\tiny
```{r}
fit1 = glm(time~lot+log(u),data=clot,family=Gamma(link=log))
summary(fit1)
```
\normalsize

**Q**: describe what you see in the print-out.


# Gamma regression: likelihood and derivations thereof

**Likelihood:**
$$L(\beta)=\prod_{i=1}^n \exp(-\frac{\nu y_i}{\mu_i}-\nu \ln \mu_i +\nu \ln \nu +(\nu-1)\ln y_i -\ln(\Gamma(\nu)))$$

**Log-likelihood:**
$$l(\beta)=\sum_{i=1}^n [-\frac{\nu y_i}{\mu_i}-\nu \ln \mu_i +\nu \ln \nu +(\nu-1)\ln y_i -\ln(\Gamma(\nu))]$$
Observe that we now- for the first time - have a nuisance parameter $\nu$ here. 

---

# Fitting the Model

To produce numerical estimates for the parameter of interest $\beta$ we may proceed to the score function, and solve using Newton Raphson or Fisher scoring. If we do not have the canonical link the observed and expected Fisher information matrix may not be equal.

What about $\phi=1/\nu$? Also estimated using maximum likelihood.

Further analyses: as before we use asymptotic distribution of parameter estimates, and of Wald, LRT and score test.

---

# Scaled and unscaled deviance

We have defined the deviance as
$$D=-2(\ln L(\text{candidate model})-\ln L(\text{saturated model}))$$

This is often called the _scaled deviance_.

The _unscaled deviance_ is then defined as $\phi D$, but is sadly sometimes also called the deviance - for example by `R`.

1. For the normal model the 

   + scaled deviance is $D=\frac{1}{\sigma^2}\sum_{i=1}^n (y_i-\hat{\mu}_i)^2$, while
   + unscaled deviance is $\phi D=\sum_{i=1}^n (y_i-\hat{\mu}_i)^2$

2. For the binomial and Poisson model $\phi=1$ so the scaled and unscaled deviance are equal.

3. What about the Gamma model?

---

Some calculations - see IL week 2, problem 2: 1b.

$$D=\frac{-2 \sum_{i=1}^n [\ln (\frac{y_i}{\hat{\mu}_i})-\frac{y_i-\hat{\mu}_i}{\hat{\mu}_i}]}{\phi}$$
and unscaled as $\phi D=-2 \sum_{i=1}^n [\ln (\frac{y_i}{\hat{\mu}_i})-\frac{y_i-\hat{\mu}_i}{\hat{\mu}_i}]$.

Compare to print-out from R: the deviance in R is the _unscaled deviance_.
```{r}
deviance(fit1)
(nu1=1/summary(fit1)$dispersion)
(D=-2*nu1*sum(log(fit1$y/fit1$fitted.values)-((fit1$y-fit1$fitted.values)/fit1$fitted.values)))
deviance(fit1)*nu1
```

# Comparing models 

## Comparing models based on deviance

```{r}
fit2=glm(time~lot+log(u)+lot:log(u),data=clot,family=Gamma(link=log))
anova(fit1,fit2)
```
The deviance table does not include $\phi$, so the unscaled deviance is reported.
If significance testing is done, the estimated $\phi$ from the largest model is used, and $p$-values are based on the scaled deviance.

---

```{r}
anova(fit1,fit2,test="Chisq")
1-pchisq((deviance(fit1)-deviance(fit2))/summary(fit2)$dispersion,fit1$df.residual-fit2$df.residual)
anova(fit1,fit2,test="F")
1-pf((deviance(fit1)-deviance(fit2))/summary(fit2)$dispersion,fit1$df.residual-fit2$df.residual,fit2$df.residual)
```

---

## Comparing models based on AIC
```{r}
AIC(fit1,fit2)
```
**Q**: would you prefer `fit1` or `fit2`?

---

AIC can also be used when we compare models with different link functions (models that are not nested).

The literature suggests to plot $y_i$ vs. each covariate to get a hint about which link function or transformation to use.

* Identity: Plot of $y_i$ vs $x_i$ should be close to linear
* $\ln:$ Plot of $\ln(y_i)$ vs $x_i$ should be close to linear
* Inverse (reciprocal): Plot of $1/y_i$ vs $x_i$ should be close to linear

---

# Compare link functions

```{r}
library(ggplot2)
library(ggpubr)
y=clot$time
x=clot$u

df=data.frame(y=y,x=x)
gg1=ggplot(df)+geom_point(aes(x=log(x),y=y)) + ggtitle("Identity")
gg2=ggplot(df)+geom_point(aes(x=log(x),y=log(y))) + ggtitle("Log")
gg3=ggplot(df)+geom_point(aes(x=log(x),y=1/y)) + ggtitle("Inverse")
gg4=ggplot(df)+geom_point(aes(x=sqrt(1/x),y=log(y))) + ggtitle("Log, sqrt x")
```

---

```{r}
ggarrange(gg1,gg2,gg3,gg4)
```

--- 

```{r}
fit4=glm(time~lot+sqrt(1/u),data=clot,
         family=Gamma(link=log))
AIC(fit1,fit4)
```

---

# Code for Residual Plots

```{r}
df4 = data.frame(fitted = fit4$fitted.values, dres = residuals(fit4, 
    type = "deviance"))
gg4 = ggplot(df4) + 
geom_point(aes(x = fitted, y = dres)) +
scale_color_discrete("") + ggtitle("time~lot+sqrt(1/u)")
df1 = data.frame(fitted = fit1$fitted.values, dres = residuals(fit1, 
    type = "deviance"))
gg1 = ggplot(df1) + 
geom_point(aes(x = fitted, y = dres))+
scale_color_discrete("") + ggtitle("time~lot+log(u)")
```

# The Plots

Are these good?

```{r, fig.height=5.5}
ggarrange(gg1,gg4)

```

---

# R packages

```{r, eval=FALSE}
install.packages(c("tidyverse", 
  "ggplot2", 
  "statmod",
  "corrplot", 
  "ggplot2", 
  "GGally",
  "boot"))
```

# Further reading
* A. Agresti (1996): "An Introduction to Categorical Data Analysis".
* 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.

