---
title: "TMA4315 Generalized linear models H2018"
subtitle: "Module 4: Count and continuous positive response data (Poisson and gamma regression)"
author: "Mette Langaas, Department of Mathematical Sciences, NTNU -- with contributions from Ingeborg Hem"
date: " 27.09.2018 and 04.10.2018 [PL], 28.09.2018 and 05.10.2018 [IL]"
output: #3rd letter intentation hierarchy
  html_document:
    toc: true
    toc_float: true
    toc_depth: 2
    css: "../TMA4315RMarkdown.css"
 # pdf_document:
 #   toc: true
 #   toc_depth: 2
 #   keep_tex: yes
#  beamer_presentation:
#    keep_tex: yes
#    fig_caption: false
#    latex_engine: xelatex
---

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

<details><summary>Answer/Hint (not yet added)</summary>
42
</details>

# <a id="interactivew2">Interactive session - second week</a>

## Problem 1: Exam 2007 (Problem 1, a bit modified) - Smoking and lung cancer

(Permitted aids for the exam was "Tabeller og formler i statistikk", Matematisk formelsamling (Rottmann), one A5 sheet with your own handwritten notes, and a simple calculator.)

The dataset given in `smoking.txt` consists of four variables:

* `deaths`: number of lung cancer deaths over a period of six years [remark: incorrectly 1 year in exam question]
* `population`: the number of people [remark: incorrectly in 100 000 people in exam question]
* `age`: in five-year age groups (`40-44, 45-49, 50-54, 55-59, 60-64, 65-69, 70-74, 75-79, 80+`)
* `ageLevel`: age group as numbers from 1 to 9 (1 corresponds to 40-44, 2 to 45-59, and so on)
* `smoking status`: doesn't smoke (`no`), smokes cigars or pipe only (`cigarPipeOnly`), smokes cigarettes and cigar or pipe (`cigarettePlus`), and smokes cigarettes only (`cigaretteOnly`)

You can look at the dataset [here](https://www.math.ntnu.no/emner/TMA4315/2018h/smoking.txt). The data set is probably taken from [Best & Walker (1964)](https://www.jstor.org/stable/41983444?seq=1#page_scan_tab_contents), who say that the people studied were males who have contributed in wars before 1956 and who answered a questionaire. The authors point out that the dataset is not representative for the whole population.

We are interested in studying if the mortality rate due to lung cancer (the number of deaths due to lung cancer per individual during six years) controlled for age group and smoking status. Assume that the number of deaths for each set of covariate values, $Y_i$, can be considered Poisson distributed, $Y_i \sim \text{Poisson}(\lambda_i)$. We fit the following model:

```{r}
# load data:
smoking <- read.table(file = "https://www.math.ntnu.no/emner/TMA4315/2018h/smoking.txt")
# Correct typo. There are a few ways of doing this
smoking$age <- gsub("45-5", "45-4", smoking$age)
head(smoking)
nrow(smoking)
model1 <- glm(dead ~ age + smoke, family = "poisson", data = smoking, offset = log(pop)) 
# note that the size of the population for each combination of the covariates is the offset here
summary(model1)
```


<span class="question">a) What is an offset?</span>

<details><summary>Answer</summary>
In a model, an offset can be seen as a covariate with a fixed effect (usually 1). It is used because the fundamental structure of the problem means the effect size is known.

The problem here is a typical example: we expect the number of deaths to be proportional to the population size, i.e. if the probability of an individual dying is $p$, the expected number of deaths, $\text{E}(n)$ should be proportional to $N$, the population size. So $\text{E}(n) = \mu = Np$, and hence in a GLM with a log link we get $\eta = \log N + \log p$, so we can use $\log N$ as an offset.
</details>

<span class="question">For 53-year old non-smokers, what is the estimated number of deaths due to lung cancer (per person over 6 years)?  </span>

<details><summary>Answer</summary>
First by hand (as you would have to do in an exam)

$\mu = \exp(-3.63222 + 0.98039 -0.04781  + \log(249)) = \exp(2.817813) = 16.7402$

Then, in R. The first line predicts the mean for all of the data, and the second extracts the number we want.

```{r}
smoking$predict <- predict(model1, type="response")
smoking$predict[smoking$age == "50-54" & smoking$smoke=="no"]
```

</details>

<span class="question">Why is the number of degrees of freedom for the deviance of this model 24?</span>

<details><summary>Answer</summary>
Because we have `r nrow(smoking)` data points, and 12 parameters (an intercept, 8 age contrasts to age `40-44`, and 3 contrasts to smoke level `cigarPipeOnly`).

As an aside, I would make the `no` smoking level the reference level.
</details>

<span class="question">Does the model give a good fit?</span>

<details><summary>Answer</summary>
Yes. We can see that the residual deviance is 21.5. If the model gives a good fit, this should follow a chi-squared distribution with 24 degrees of freedom. Without calculating the p-value, we can see it will not be significant, because the expected value is the degrees of freedom, and the residual deviance would have to be larger than this to reject the null hypothesis of a good fit.

Formally, we can calculate the p-value:

```{r}
pchisq(model1$deviance, model1$df.residual, lower.tail = FALSE)
```

And, yes, it is not significant.

</details>

b) Let $\lambda(a,s)$ denote the expected number of lung cancer deaths per person in age group $a$ with smoking status $s$. For two different smoking statuses $s_1$ and $s_2$, define

$$r(a, s_1, s_2) = \frac{\lambda(a, s_1)}{\lambda(a, s_2)}.$$

<span class="question">Explain why $r(a, s_1, s_2)$ does **not** vary as a function of $a$ in `model1`.  </span>

<details><summary>Hint</summary>
The most difficult part of this is understanding the question!
</details>
<details><summary>Answer</summary>
Because there is no interaction in the model.
</details>

<span class="question">For $s_1 =$ `cigarPipeOnly` and $s_2 =$ `cigaretteOnly`, find an estimate value for $r(a, s_1, s_2)$ and an approximate 90 % confidence interval.</span>

<details><summary>Hint</summary>
This is a bit simpler than it appears. Work out how you would calculate predictions for each level.
</details>
<details><summary>Answer</summary>
The long version starts by noting (with a slight abuse of notation) that $\lambda(a, s_1) = \exp(\beta_0 + \beta_a + \beta_{s_1})$, so 

$$
\begin{aligned}
r(a, s_1, s_2) &= \frac{\lambda(a, s_1)}{\lambda(a, s_2)}\\
&= \frac{\exp(\beta_0 + \beta_a + \beta_{s_1})}{\exp(\beta_0 + \beta_a + \beta_{s_2})} \\
&= \exp((\beta_0 + \beta_a + \beta_{s_1}) -(\beta_0 + \beta_a + \beta_{s_2})) \\
&=\exp(\beta_{s_1} - \beta_{s_2})
\end{aligned}
$$

This is easier than it appears: `cigarPipeOnly` is the reference class, so $\beta_{s_1}=0$ and we only need the `cigaretteOnly` effect. From above this is 0.36915, with a standard error of 0.03791.

So $r(a, s_1, s_2) = e^{-0.36915} = 0.69$. And the the 90% confidence interval is 
$e^{-0.36915 \mp 1.645\times0.03791} = (0.65, 0.74)$

We can compare this with the values from R:

```{r}
exp(-coef(model1)["smokecigarretteOnly"])
exp(-confint(model1, level=0.9)["smokecigarretteOnly",])

```
</details>

<span class="question">Is there a significant difference in the expected number of lung cancer deaths for individuals that smoke cigarettes `cigaretteOnly` versus those that smoke cigar/pipe `cigarPipeOnly`?</span>

<details><summary>Answer</summary>
This depends on what level $\alpha$ is used! If we use 10%, then from the previous question we can say no, because the confidence intervals do not overlap 0. This is the same as using a Wald test.

If we read the z test from the summary table we see that the z value is 9.737, with an associated p-value which is < 2e-16 (in fact it is about $10^{-22}$). So for any sane $\alpha$ there is a significant difference.
</details>

**c)**
We will now consider two alternative models, `model2` and `model3`:

```{r}
model2 <- glm(dead ~ smoke, family = "poisson", data = smoking, offset = log(pop))
model3 <- glm(dead ~ ageLevel + smoke, family = "poisson", data = smoking, offset = log(pop))

summary(model2)
summary(model3)
```


<span class="question">In `model3` we have `ageLevel`, which is continuous, as described above it is the age group as numbers from 1 to 9 (1 corresponds to 40-44, 2 to 45-59, and so on). So how do you interpret the estimated effect?</span>

<details><summary>Answer</summary>
The age classes are all 5 years wide, so going up 1 class is the same as aging 5 years (on average). Thus the estimated effect is the of of the increase in death rate for doctors 5 years of age apart. There are a few equivalent ways of expressing this, so your answer might be a bit different, but still be correct.

This is the sort of weird things people do in analyses, for reasons which made sense at the time. So working out the interpretation can mean a bit of poking around.
</details>

<span class="question">Why does `model2` and `model3` have 32 and 31 degrees of freedom, respectively?</span>

<details><summary>Answer</summary>
For both of these, there are 36 data points. `model2` estimates 4 parameters (the intercept and 3 for the `smoke` factor), leaving 32 degrees of freedom. `model3` uses an additional parameter to estimate the ageLevel slope, so there are 31 left.

</details>

<span class="question">If we want to compare the three models `model1`, `model2` and `model3`, which model would you choose as the best? Justify your answer by formulating relevant hypotheses and perform hypothesis tests, based on the print-outs above.</span>

<details><summary>Hint</summary>
Remember that to compare two models then one model needs to be nested within the other model.
</details>

<details><summary>Answer</summary>
We could use AIC to compare the models:

```{r}
AIC(model1, model2, model3)
```

Lower is better, so `model1` is the best, by a long way. BIC would give the same result.

Although it might not look like it, we could also use Analysis of Deviance (i.e. ANOVA for GLMs). The trick is to see that `model2` is nested within `model1`, even though the age variables are different. The only difference is that `ageLevel` is ordered, as we assume the effect is linear, whereas `age` assumes no ordering. We could actually fit an 8th order polynomial to `ageLevel`and get the same model.

If we compare the models we see, again that `model1`is best:

```{r}
anova(model2, model3, model1)

```



</details>

New: For `model3`: <span class="question">Plot the regression line for the expected number of lung cancer deaths per person as a function of `age` for each of the four different smoking levels in the same plot. First use all coefficients, and then only the ones that are significant. Use `ggplot`! Discuss what you see.</span>

<details><summary>Answer/Hint (not yet added)</summary>
42
</details>

<!--
**d)** The (asymptotic) covariance matrix for the regression parameters (in a GLM) are given as $({\bf X}^T{\bf W}{\bf X})^{-1}$ where ${\bf X}$ is the design matrix for the problem.
Find the difference in the ${\bf W}$ matrix with log-link and identity (linear) link for a Poisson model. Also find the difference in the ${\bf W}$ matrix with logit link and identity (linear) link for binomial data.
Remark: We are talking about the inverse of the expected Fisher information matrix, and what is meant is that you find the diagonal elements of ${\bf W}$ for each of the four cases. And, we have not seen the general formula for ${\bf W}$ yet - but will do in Module 5. 
-->

<!--
**1. Explain what is done in the `R`-print-out.
2. What if we instead want a CI for $\beta_0+\beta_1 x_{i1}$?
3. What if we instead want a CI for $\lambda_i=\exp(\beta_0+\beta_i x_{i1})$?
**d)** Overdispersion
-->


## Problem 2: TMA4315 Exam 2012, Problem 3: Precipitation in Trondheim, amount 
Remark: the text is slightly modified from the original exam since we parameterized the gamma as in our textbook.

Remark: Problems 1 (binomial) and 2 (multinomial) at the 2012 exam also asked about precipitation.

Remark: There is no data provided in this problem: it is purely theoretical (unlike the rain in Trondheim).


We want to model the amount of daily precipitation given that it *is* precipitation, and denote this quantity $Y$. It is common to model $Y$ as a gamma distributed random variable, $Y \sim Gamma(\nu,\mu)$, with density

$$ f_Y(y) = \frac{1}{\Gamma(\nu)} \left(\frac{\nu}{\mu}\right)^{\nu} y^{\nu-1}\exp\left(-\frac{\nu}{\mu} y \right) $$

In this problem we consider $N$ observations, each gamma distributed with $Y_i \sim Gamma(\nu, \mu_i)$ (remark: common $\nu$). Here $\nu$ is considered to be a known nuisance parameter, and the $\mu_i$s are unknown.

**a)** <span class="question">Show that the gamma distribution function is member of the exponential family when $\mu_i$ is the parameter of interest.</span>


<details><summary>Hint</summary>

You need to be able to write it in this form

$$
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)
$$
</details>

<details><summary>Answer</summary>

We want the likelihood in this form:

$$
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)
$$
it's easier to just look at the parts in the brackets, i.e. look at the log likelihood:

$$
\begin{aligned}
\log f_Y(y) &= -\log\Gamma(\nu) +\nu \log \left(\frac{\nu}{\mu}\right) +  (\nu-1)\log y - \frac{\nu}{\mu} y \\
&= - \frac{\nu}{\mu} y + \nu \log \frac{1}{\mu} + \nu \log {\nu} -\log\Gamma(\nu)  +  (\nu-1)\log y  \\
&= \frac{ -y\frac{1}{\mu} + \log \frac{1}{\mu}}{\frac{1}{\nu}} + \nu \log {\nu} -\log\Gamma(\nu)  +  (\nu-1)\log y
\end{aligned}
$$
So

$$ 
\begin{aligned}
\theta &= -\frac{1}{\mu} \\
b(\theta_i) &= -\log \frac{1}{\mu} \\
\phi &= \frac{1}{\nu} \\
c(y_i, \phi, w_i) &= \nu \log {\nu} -\log\Gamma(\nu)  +  (\nu-1)\log y
\end{aligned}
$$

</details>

<span class="question">Use this to find expressions for the expected value and the variance of $Y_i$, in terms of $(\nu,\mu_i)$, and interpret $\nu$.</span>

<details><summary>Answer/Hint (not yet added)</summary>
42
</details>

New: <span class="question">What is the canonical link for the Gamma distribution?</span>

<details><summary>Answer/Hint (not yet added)</summary>
$\theta &= -\frac{1}{\mu}$

So the negative inverse link.

</details>

Hint: if you want to focus on discussing - you may look at the solutions from [Module 1 together](https://www.math.ntnu.no/emner/TMA4315/2017h/Module1ExponentialFamily.pdf).

**b)** <span class="question">Explain what a saturated model is.</span>

<details><summary>Answer</summary>
A model where there is one parameter per data point.
</details>

<span class="question">Set up the log-likelihood function expressed by $\mu_i$, and use it to find the maximum likelihood estimators for $\mu_i$s of the saturated model.  </span>

<details><summary>Hint</summary>
Hint: Do you see directly that $\hat{\mu}_i=y_i$? If not, you may look at the likelihood for one observation and solve that the derivative equals 0, i.e. $s(\mu_i)=0$: follow how this was done for other distributions.
</details>

<details><summary>Answer</summary>
We want to differentiate with respect to $\mu_i$ to get the score, and solve for when this equals 0. We can do this for any one data point (and not worry about summing likelihoods!):

$$
\begin{aligned}
s(\mu_i) &= \frac{\partial \log f_Y(y_i)}{\partial \mu_i} \\
&= \frac{\partial {\frac{ -y_i\frac{1}{\mu} + \log \frac{1}{\mu}}{\frac{1}{\nu}} + \nu \log {\nu} -\log\Gamma(\nu)  +  (\nu-1)\log y_i}}{\partial \mu_i} \\
&= \nu\left(-\frac{\partial y_i \frac{1}{\mu_i}}{\partial \mu_i} + \frac{\partial \log\frac{1}{\mu_i}}{\partial \mu_i} \right)
\end{aligned}
$$

We can use the chain rule, with $\theta_i = \mu_i^{-1}$, so 

$$
\begin{aligned}
\frac{\partial \theta_i}{\mu_i} &= -\mu_i^{-2} \\
\frac{\partial y \theta_i}{\partial \theta_i} &=y_i \\
\frac{\partial \log\theta_i}{\partial \theta_i} &=\frac{1}{\theta_i}
\end{aligned}
$$

$$
\begin{aligned}
s(\mu_i) &= \frac{\partial \log f_Y(y_i)}{\partial \theta_i}\frac{\partial \theta_i}{\mu_i} =0\\
&= \nu\left(-y_i + \mu_i \right) \left(-\frac{1}{\mu_i^2} \right) \\
&= -y_i + \mu_i =0 \\
y_i &= \mu_i
\end{aligned}
$$

</details>

<span class="question">Find the deviance (based on all $N$ observations).</span>

<details><summary>Answer</summary>
The For a model that gives us estimates of $\mu_i$, the deviance is 

$$
D({\bf y}|{\boldsymbol \mu}) = -2(\sum_{i=1}^nl(y_i|\mu_i) - l(y_i|y_i) )
$$

For one datum for the gamma distribution this is 

$$
\begin{aligned}
D_i(y_i|\mu_i) &= -2(l(y_i|\mu_i) - l(y_i|y_i))\\
&= -2 \left(\frac{ -y_i\frac{1}{\mu_i} + \log \frac{1}{\mu_i}}{\frac{1}{\nu}} - \frac{ -y_i\frac{1}{y_i} + \log \frac{1}{y_i}}{\frac{1}{\nu}} \right) \\
&= - 2 \nu\left(-\frac{y_i}{\mu_i} + \frac{y_i}{y_i} +\log \frac{1}{\mu_i} - \log \frac{1}{y_i }\right) \\
&= - 2 \nu \left(1-\frac{y_i}{\mu_i} + \log \frac{y_i}{\mu_i} \right)
\end{aligned}
$$

So for all $N$ observations the deviance is 

$$
\begin{aligned}
D_i({\bf y}| {\boldsymbol \mu}) &= - 2 \nu \sum_{i=1}^{N}\left(1-\frac{y_i}{\mu_i} + \log \frac{y_i}{\mu_i} \right) \\
&= - 2 \nu N + 2 \nu \sum_{i=1}^{N}\frac{y_i}{\mu_i} - \sum_{i=1}^{N}\log \frac{y_i}{\mu_i}
\end{aligned}
$$

(you can rewrite this in a couple of ways depending on what you think looks nicest)

</details>


**c)** We now want to construct a model for amount of precipitation (given that there are occurrence) with precipitation forecast as explanatory variable.  
<span class="question">Let $Y_i$ be amount of precipitation for day $i$, and let $x_i$ be the precipitation forecast valid for day $i$. Set up a GLM for this, and argue for your choice of link function and linear predictor.</span>

<details><summary>Hint</summary>
Hint: explain why this fits into our GLM-gamma framework and set up our 3 equation model (random, systematic, link).
</details>

<details><summary>Answer</summary>
- **Random**: We assume $Y_i$ follows a Gamma distribution with mean $\mu_i$
- **Systematic**: We assume $g(\mu_i) = \eta_i = \beta_0 + \beta_1/x_i$
- **Link**: We assume $g(\mu_i) -1/\eta_i$

The link function is the canonical link, so is what makes `glm()` happiest. The linear predictor is not quite obvious: it's linear in $1/x_i$ to conter-act the inverse link. We would expect $\beta_0=0$ and $\beta_1 = -1$, but only if we trust the forecasters. 
</details>


---

## Problem 3:  Taken from UiO, STK3100, 2015, problem 2
(For reference: [here is the original exam](http://www.uio.no/studier/emner/matnat/math/STK3100/oppgaver/stk3100-f15eng_final.pdf)).

This is a problem on the logistic regression.

Do not look at the dataset before the end of the exercise! You should solve the exercise without using `R`, just as if you were at an exam.

In this problem you shall consider data of survivals from a study of treatment for breast cancer. The response is the number that survived for three years. The covariates were the four factors

* `app`: appearance of tumor, two levels, 1 = malignant, 2 = benign
* `infl`: inflammatory reaction, two levels, 1 = minimal, 2 = moderate or severe
* `age`: age of patients, three levels, 1 = under 50, 2 = 50 to 69, 3 = 70 or older
* `country`: hospital of treatment, three levels, 1 = Japan, 2 = US, 3 = UK

The dataset we have used differs slightly from the one they used at UiO (the number of survivals and non-survivals differ).

The number of survivors is modelled as a binomially distributed variable using a canonical logit link. Dummy variable coding is used for all factors, with the level coded as "1" as the reference category.

**a)** The output from fitting the model where only appearance and country are used as covariates, i.e., a model with predictor on the form

$$ \eta = \beta_0 + \beta_1 \texttt{ fapp} + \beta_2 \texttt{ fcountry2} + \beta_3 \texttt{ fcountry3} $$

is displayed below. 

<span class="question">What is the interpretation of the estimate of the coefficient of appearance, `fapp` (`f` means factor)? Explain also how this coefficient can be expressed in terms of an odds ratio.</span>

<details><summary>Answer/Hint (not yet added)</summary>
42
</details>

<details><summary>Answer/Hint (not yet added)</summary>
42
</details>

Hint: first explain what is the predicted odds of survivial for country $j$ for a benign tumor, then the same for a malignant tumor- and then make an odds ratio.

New: The main difference between our dataset and UiO's dataset is the number of degrees of freedom for the null model; we have 34, they have 35. 

<span class="question">How many observations do we have in the dataset? And how many does UiO have?</span>

<details><summary>Answer/Hint (not yet added)</summary>
42
</details>

<details><summary>Answer/Hint (not yet added)</summary>
42
</details>

Hint: All the covariates in the dataset are categorical, two with 3 levels, and two with 2 levels. How many possible combinations of observations does that make?

```
Call:
glm(formula = cbind(surv, nsurv) ~ fapp + fcountry, family = binomial, 
    data = brc)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8142  -0.7279   0.2147   0.7576   1.8715  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.0803     0.1656   6.522 6.93e-11 ***
fapp2         0.5157     0.1662   3.103 0.001913 ** 
fcountry2    -0.6589     0.1998  -3.298 0.000972 ***
fcountry3    -0.4945     0.2071  -2.388 0.016946 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 57.588  on 34  degrees of freedom
Residual deviance: 36.625  on 31  degrees of freedom

```

**b)** The output below is an analysis of deviance for comparing various model spesifications. <span class="question">Fill out the positions indicated `i-iv`.  </span>
Note: `y ~ x1*x2` gives both the linear components and the interaction, i.e., it is the same as `y ~ x1 + x2 + x1:x2` (`x1:x2` gives interaction only).  

<details><summary>Answer/Hint (not yet added)</summary>
42
</details>

<details><summary>Answer/Hint (not yet added)</summary>
42
</details>

Remark: remember that `anova` gives the sequential comparison of deviances - that is - comparing each model to the previous model.


```
Analysis of Deviance Table

Model 1: cbind(surv, nsurv) ~ fapp + fage + fcountry
Model 2: cbind(surv, nsurv) ~ fapp + fage + finfl + fcountry
Model 3: cbind(surv, nsurv) ~ fapp + finfl + fage * fcountry
Model 4: cbind(surv, nsurv) ~ fapp * finfl + fage * fcountry
Model 5: cbind(surv, nsurv) ~ fapp * finfl + fapp * fage + fage * fcountry
Model 6: cbind(surv, nsurv) ~ fapp * finfl * fage * fcountry
  Resid. Df Resid. Dev Df Deviance
1        29     33.102            
2         i     33.095  1   0.0065
3        24     25.674 ii   7.4210
4        23     25.504  1      iii
5        21     22.021  2   3.4823
6         0      0.000 iv  22.0214
```

In the remaining parts of this problem we return to the model in part a) and consider the hypothesis

$$ H_0 : \beta_2 + \beta_3 = -1 \text{ versus } H_1 : \beta_2 + \beta_3 \neq -1 $$

Note: what are you testing now?

**c)** The estimated covariance matrix between the estimators of the coefficients $\beta_2$ and $\beta_3$ above is $\left(\begin{matrix} 0.040 & 0.021 \\ 0.021 & 0.043\end{matrix} \right)$. <span class="question">Use a Wald test to test the null hypothesis above.</span>

<details><summary>Answer/Hint (not yet added)</summary>
42
</details>

Note: remember - this was a written exam so you do this by hand!

**d)** <span class="question">Explain how the null hypothesis can be tested with a likelihood ratio test by fitting two suitable models.</span> No numerical calculations are necessary, but it must be specified how the predictors should be defined.  
Remark: rather technical - and hint: offset term?

<details><summary>Answer/Hint (not yet added)</summary>
42
</details>

<details><summary>Answer/Hint (not yet added)</summary>
42
</details>



New: <span class="question">Do this test in `R`.</span>

[Here is the data set!](https://www.math.ntnu.no/emner/TMA4315/2018h/breastcancer.txt)

<details><summary>Answer/Hint (not yet added)</summary>
42
</details>

# Work on your own: Exam questions

## December 2013 (Essay exam)

We will consider the following Poisson regression
$$Y_i \sim \text{Poisson}(\exp(\eta_i)), \text{ }i=1,\ldots,n$$
where the linear predictor is $\eta_i={\bf x}_i^T\beta$. Here ${\bf x}_i$ is a vector of the $p$ covariates for the $i$th observation $Y_i$ and $\beta$ is unknown $p$ dimensional column vector of unknown regression coefficients.

Write an introduction to Poisson regression and its practical usage, for a student
with a good background in statistics, but no knowledge about Generalized Linear
Models (GLM). Topics you may want to consider, are

* When to use it? Underlying assumptions.
* Parameter estimation, limiting results for the MLE, Fisher information and
observed Fisher information, confidence intervals and hypothesis testing.
* Output analysis, residual plots and interpretation of results.
* Deviance and its usage.
* What do we do when a covariate is a factor, and should the results be
interpreted? 
* The use of Poisson regression in the analysis of contingency tables.

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

