---
title: "TMA4315 Generalized linear models H2018"
subtitle: "Module 3: BINARY REGRESSION"
author: "Mette Langaas, Department of Mathematical Sciences, NTNU - with contributions from Øyvind Bakke, Thea Bjørnland and Ingeborg Hem"
date: "13.09 and 20.09 [PL],  14.09 and 21.09 [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)
showsol<-TRUE
```

**Hints and reminders are in bold**
 
<span class="question">Questions appear in blue.</span>

<details><summary>Some hints and answers are hidden away in a fold like this</summary>
Well done, you found me!
</details>

The exam style questions are optional, but will be helpful for, well, the exam. They will also help you now to see if you understand the content.

# Interactive lecture - second week

We will use a data set on contraceptive use in Fiji (data from 1975). The data is to be analysed with 
"current use of contraceptive" as response and some (or all of) "age", "level of education", "desire for more children" as covariates.

The data is available at <https://grodri.github.io/datasets/cuse.dat> with the following description:
  
* Contraceptive use: yes (`using`) or no (`notUsing`)
* `age`: categorical variable with 5 levels: "<25", "25-29", "30-39", "40-49"
* `education`: categorical variable with 2 levels giving highest level of education obtained: Lower and Upper 
* `wantsMore`: Desires more children: yes or no


```{r,eval=TRUE}
ds=read.table("https://grodri.github.io/datasets/cuse.dat",header=TRUE)
names(ds)
summary(ds)
dim(ds)
head(ds)
```

We will study binary regression using the logit model, and we will work with grouped data. So the proportion using contraception is `using`/(`using` + `notUsing`). For grouped data as a response, R wants 2 columns, with success in the fist column, and failure in the second. We can write this as `cbind(success, failure) ~ X1 + X2`.

<!-- Useful functions: `qlogis(p)` calculates the logit logit(p) = log(p/(1-p)) and the function `plogis(x)` is the inverse logit, that is exp(x)/(1-exp(x)). -->

**Plan: Start with Problem 2, then move to 3 and 4, and if you have time you look at Problem 1. **

## Problem 1: The null model - no covariates included.

(This is the most theoretical of the problems and rather technical - but with some cool results on the null model.)

<span class="question">**a)** Fit the null model and call it `fit0`. Explain what you see. </span>

<details><summary>Answer</summary>
```{r}
fit0=glm(cbind(using,notUsing)~1,data=ds,family=binomial(link=logit))
summary(fit0)
```

We see some output. The estimate of the intercept is `r round(coef(fit0),2)`: this is less than 0.5, so the mean probability that a woman uses contraception is <0.5 (the response is written as `cbind(success, failure)`: there are other ways of writing it, see the help page for `glm`).

The null and residual deviances are the same, because this is the null model.

</details>

Observe: When we only have an intercept in the model, the model matrix will be an $n \times 1$-matrix with only ones. 
Then ${\bf x}_j^T\boldsymbol{\beta}=\beta_0$. The log-likelihood can be written as (let $j=1,\ldots,G$)

\begin{align} l(\boldsymbol{\beta}) &=\sum_{j=1}^G (y_j{\bf x}_j^T \boldsymbol{\beta} - n_j \log(1+\exp({\bf x}_j^T\boldsymbol{\beta})))= \sum_{j=1}^G (y_j\beta_0 - n_j \log(1+\exp(\beta_0)))\\
&= \beta_0 \sum_{j=1}^G y_j - \log(1+\exp(\beta_0))\sum_{j=1}^G n_j =\beta_0 N_1 - \log(1+\exp(\beta_0))N
\end{align}

where $N_1=\sum_{j=1}^G y_j$ is the total number of successes and $N=\sum_{j=1}^G n_j$ is the total number of trials (over all covariates patterns, that is, here the number of individuals in the data set). Also $N_2=N-N_1$ is the total number of failures.

We will use this loglikelihood in the next question.

<span class="question">**b)** What is the relationship between your estimated coefficient and the proportion in the data set using contraception ($N_1/N$)?</span>

<details><summary>Hint 0</summary>
What would be your intuitive estimator for $\pi$ (common to all groups). </details>
<details><summary>Hint 1</summary>
Find the derivative of the log-likelihood with respect to $\beta_0$, and use this to find the MLE for $\beta_0$.
</details>
<details><summary>Hint 2</summary>
maybe easier to see what $\hat{\pi}$ is, where $\hat{\pi}=\frac{\exp{\hat{\beta}_0}}{1+\exp{\hat{\beta}_0}}$ (so `plogis`), and then $\hat{\boldsymbol{\beta}}_0=\text{logit}(\hat{\pi})$ (so `qlogis`). 
</details>
<details><summary>Hint 3</summary>
You can verify by using the estimate from `glm`.
</details>
<details><summary>Answer</summary>

The intuitive estimate of $\hat{\pi}$ is $N_1/N$. From above, $l(\boldsymbol{\beta}) = \beta_0 N_1 - \log(1+\exp(\beta_0))N$, so we can differentiate (to get the score), set that to 0 and solve...

$$
\begin{align} \frac{dl(\boldsymbol{\beta})}{d\beta_0} &=  N_1 - \frac{\log(1+\exp(\beta_0))N}{d\beta_0}
 \\
&= N_1 - N \frac{\log(u)}{du}\frac{(1+\exp(\beta_0))}{d\beta_0} && \text{using } u=1+\exp(\beta_0) \\
&= N_1 - N \frac{1}{1+\exp(\beta_0)}\exp(\beta_0) = 0 \\
 \hat{\pi} &= \frac{N_1}{N} = \frac{\exp(\beta_0)}{1+\exp(\beta_0)}
\end{align}
$$

Recognise this? it's the logistic function, so the inverse is the logit link: $\beta_0 = \log(\hat{\pi}/(1-\hat{\pi}))$.

We can verify this in R:

```{r}
N = sum(ds$using + ds$notUsing)
N1 = sum(ds$using)
N2 = N - N1

log(N1/N2)
fit0$coefficients
```

</details>


**c)** We know that the (asymptotic) estimated covariance matrix of the ML estimate is the inverse of the expected Fisher information matrix, here the matrix is only a scalar and the covariance matrix is just the variance. 

<span class="question">Find the mathematical expression for the estimated variance of our estimated intercept.</span>

<details><summary>Hint 1</summary>
We have $F(\boldsymbol{\beta}) = \sum_{j=1}^G x_jx_j^T n_j \pi_j(1-\pi_j)$, and then insert $x_j = 1$ and $\pi_j(1-\pi_j)=\pi(1-\pi)$, and hopefully you found above that $\hat{\pi}=N_1/N$ in our model with only intercept. 
</details>
<details><summary>Hint 2</summary>
$\frac{N}{N_1 \cdot N_2}=\frac{1}{N_1}+\frac{1}{N_2}$ to make things prettier.
</details>
<details><summary>Answer</summary>

$Var(\boldsymbol{\beta}) = F^{-1}(\boldsymbol{\beta})$, so find $F(\boldsymbol{\beta})$:

$$
\begin{aligned}
F(\boldsymbol{\beta}) &= \sum_{j=1}^G x_jx_j^T n_j \pi_j(1-\pi_j) \\
&= N \pi(1-\pi) \\
&= N \frac{N_1}{N}\left(\frac{N_2}{N}\right) \\
&= \frac{N_1 N_2}{N}
\end{aligned}
$$

So, $F^{-1}(\boldsymbol{\beta}) = Var(\boldsymbol{\beta}) = N/(N_1 N_2) = \frac{1}{N_1}+\frac{1}{N_2}$
</details>


<span class="question">**d)** What is the estimated (numerical value) standard deviation of the parameter estimate? Did your calculation above gives the same result?</span>

<details><summary>Hint</summary>
`vcov(fit0)`
</details>
<details><summary>Answer</summary>

```{r}
vcov(fit0)
```

Hopefully, "yes" (up to rounding error), especially if `1/sum(ds$notUsing) + 1/sum(ds$using) = ` `r 1/sum(ds$notUsing) + 1/sum(ds$using)`
</details>

<span class="question">**e)** What is the asymptotic distribution of the estimated regression coefficient? Use this to write down the formula for the 95% confidence interval for the intercept, and calculate the interval in R. Compare numerically to `confint.default` and `confint`.</span>

<details><summary>Answer</summary>
The coefficient is asymptotically normally distributed, so the (approximate) 95% confidence interval is

$\hat{\beta}_0 \pm 1.96 \sqrt{\sigma^2_{\beta_0}} =$ `r round(coef(fit0),2)` $\pm 1.96 \times$ `r round(sqrt(vcov(fit0)),2)` = (`r round(coef(fit0)- 1.96*sqrt(vcov(fit0)),2)`, `r round(coef(fit0)+ 1.96*sqrt(vcov(fit0)),2)`).

```{r}
c(coef(fit0) - 1.96*sqrt(vcov(fit0)), coef(fit0) + 1.96*sqrt(vcov(fit0)))
(ci=confint.default(fit0))
confint(fit0)

```

So, the asymptotic answer is (almost) the same as the one `confint.default()` gives. `confint()` (which is actually `confint.glm()`) is different in the third decimal place. It uses a different approach, something called a profile likelihood.

</details>

<span class="question">**f)** Translate the 95% CI to probability scale</span>
<details><summary>Hint</summary>
`plogis` is the inverse logit. Or write your own function.

Conceptually, think about what probability mass is above the upper limit on the linear scale and after it has been transformed.
</details>
<details><summary>Answer</summary>

If a transformation is linear, then an x% quantile before transformation is also an x% (or (100-x)%) quantile on the transformed scale. Recall that this does not hold for other statistics (e.g. the mean, mode or standard deviation).

So, we can simply transform the `ci`:

```{r}
plogis(ci)
fit0$family$linkinv(ci)
```

</details>



## Problem 2: We then study the effect of the covariate "wantsMore"
(a little less maths)

Fit a regression with `wantsMore` as covariate, and call this `fit1`.
<span class="question">**a)**  Explain what the estimated coefficient ($\beta_1$) means. </span>

<details><summary>Hint</summary>
Interpretation using odds -- if `wantsMore´ goes from 0=no to 1=yes, then...

Also, check whether a "success" is using or not using contraceptives.
</details>

<details><summary>Answer</summary>
```{r,eval=TRUE}
fit1=glm(cbind(using,notUsing)~wantsMore,data=ds,family=binomial)
summary(fit1)
exp(fit1$coefficients)
```

$\beta_1$ is the difference between the two levels of `wantsMore`, i.e. it is the difference in the log odds between women who want more children and those who do not. The estimate is negative (`r round(fit1$coefficients[2],2)`), so a woman who wants more children is less likely to be using contraceptives. The odds for a woman who wants more children are `r round(exp(fit1$coefficients[2])*100)`\% lower than they are for those who want more children, . Put the other way around, they are `r round(exp(-fit1$coefficients[2]),1)` ($=1/0.35$) times higher for women who do not want more children.
</details>

<span class="question">**b)** Is this covariate significant?</span>

- Write down the null and alternative hypothesis to test. 
- Then test using the Wald test - write down the formula yourself. 
- Use `vcov(fit1)` to access the inverse of the Fisher information matrix. 
- What is the degrees of freedom for this test? 
- Compare to the print-out from `summary`. 

<details><summary>Answer</summary>

- Write down the null and alternative hypothesis to test. 

$H_0: \beta_1=0$
$H_1: \beta_1 \ne 0$

- Then test using the Wald test - write down the formula yourself. 

Formula:

$$
w = \frac{\beta_1 - 0}{\sqrt{a_{kk}(\hat{\boldsymbol{\beta}})}} 
$$


- Use `vcov(fit1)` to access the inverse of the Fisher information matrix. 

In code:

```{r}
(w <- coef(fit1)["wantsMoreyes"]^2/vcov(fit1)["wantsMoreyes", "wantsMoreyes"])
(z <- coef(fit1)["wantsMoreyes"]/sqrt(vcov(fit1)["wantsMoreyes", "wantsMoreyes"]))

```

- What is the degrees of freedom for this test? 

1

- Compare to the print-out from `summary`. 

```{r}
summary(fit1)
```

They are the same!

</details>



**c)** Alternatively, the likelihood ratio test can be used. 

<span class="question">Write down the formula for the likelihood ratio test statistic. Use this in combination with the `fit0` and  `fit1` objects to calculate the likelihood ratio statistic. What is the $p$-value?</span>

<details><summary>Hint</summary>
Write a function to calculate the log likelihoods
</details>
<details><summary>Answer</summary>

$$
- 2\ln \lambda=-2(\ln L(\hat{\boldsymbol{\beta}}_B)-\ln L(\hat{\boldsymbol{\beta}}_A))
$$
We can calculate the log likelihoods:

```{r}
calcloglik <- 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)
}

args0 <- list(y = ds$using, n = ds$using + ds$notUsing, x = cbind(rep(1, nrow(ds))))
args1 <- list(y = ds$using, n = ds$using + ds$notUsing, x = cbind(rep(1, nrow(ds)), 
    as.numeric(ds$wantsMore=="yes")))

betas0 <- fit0$coefficients
betas1 <- fit1$coefficients

(ll0 <- calcloglik(betas0, args0))
(ll1 <- calcloglik(betas1, args1))


(LRTstat <- -2*(ll0-ll1))

```

</details>


**d)** The likelihood ratio test statistic can alternatively be calculated using the residual deviance in `fit1` and `fit0`, and is given as `fit0$deviance-fit1$deviance`. 

<span class="question">Do you see why?</span>

<details><summary>Hint</summary>
"No" may be correct, but is not the answer we are looking for. 

Also, think about deviances being additive.
</details>
<details><summary>Answer</summary>
Yes, of course I see why! The total deviance is the sum of the "explained" and residual deviance, so we can look at either the increase in deviance or the decrease in residual deviance: they are the same.

```{r}
fit0$deviance-fit1$deviance
```

</details>




<span class="question">**e)** Are the two test statistics (Wald and LRT) equal? Do the two tests give the same conclusions?</span>


<details><summary>Answer</summary>
No, they are not quite equal, but lead to the same conclusion: there is strong evidence for an effect of wanting more children on the probability of using contraception.
</details>


## Problem 3: Now two covariates - deviance and model comparison
(no maths - only definitions and print-out)

Now we study the response together with `age` and `wantsMore`. We will consider the following 5 models. See R-code and print-out below. 

```{r}
ds$Y <- cbind(ds$using, ds$notUsing)
models <- list(
  null     = glm(Y ~ 1, family=binomial, data=ds),
  age      = glm(Y ~ age, family=binomial, data=ds),
  desire   = glm(Y ~ wantsMore, family=binomial, data=ds),
  additive = glm(Y ~ age + wantsMore, family=binomial, data=ds),
  interact = glm(Y ~ age*wantsMore, family=binomial, data=ds)
)
models
lapply(models,summary)
data.frame(deviance=round(unlist(lapply(models,deviance)),2),df=unlist(lapply(models,df.residual)),aic=round(unlist(lapply(models,AIC))))
```
Note: The function `lapply` (list apply) will apply a function to each element of a list.


<span class="question">**a)** Explain what each of these models include. </span>

<details><summary>Hint</summary>
We're interested in which covariate is in which model
</details>
<details><summary>Answer</summary>

- null: Constant only, no effects. This is the same as the model in Q1
- age: only has Age as an effect
- desire: Only has whether the woman wants more children. This is the model in Q2
- additive: has two main effects: from the age and desire models,
- interact: includes an interaction, i.e. the effect of age depends on whether the woman wants more children.

</details>

<span class="question">**b)** What is the definition of the deviance and df. What can we use the deviance for? </span>

<details><summary>Hint</summary>
df = degrees of freedom

Deviance is a function of likelihoods.
</details>
<details><summary>Answer</summary>

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}


The df is the difference in the number of parameters of the models being compared.

We can use it to (a) compare models, (b) test how goodness of fit of a model, and (c) test for overdispersion.

</details>

<span class="question">Optional: derive the formula for the deviance for the logit model (the derivation is not given on the module pages.)</span>

<details><summary>Hint</summary>
The deviance is 

$$
D= -2\sum_{j=1}^G(l_j(\hat{\pi}_j)-l_j(\tilde{\pi}_j)) = 
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})]
$$

</details>
<details><summary>Answer</summary>

Because $D = -2\sum_{j=1}^G(l_j(\hat{\pi}_j)-l_j(\tilde{\pi}_j))$, we can just work with $l_j(\hat{\pi}_j)-l_j(\tilde{\pi}_j)$, where $\hat{\pi}_j$ is the estimate at our model and $l_j(\tilde{\pi}_j) = y_j/n_j$ is the estimate at the saturated model. 

\begin{aligned}
l_j({\hat{\pi}}_j)-l_j({\tilde{\pi}}_j) &= \left( y_j \ln {\hat{\pi}}_j + (1-y_j) \ln(n_j-{\hat{\pi}}_j) \right)
- \left( y_j \ln {\tilde{\pi}}_j + (1-y_j) \ln(n_j-{\tilde{\pi}}_j) \right) \\
&= y_j \ln \frac{{\hat{\pi}}_j}{{\tilde{\pi}}_j} + (1-y_j) \ln \frac{{n_j-\hat{\pi}}_j}{{1-\tilde{\pi}}_j} \\
&= y_j \ln \frac{{\hat{\pi}}_j}{{\tilde{\pi}}_j} + (1-y_j) \ln \frac{{n_j-\hat{\pi}}_j}{{1-\tilde{\pi}}_j} \\
&= y_j \ln \frac{n_j {\hat{\pi}}_j}{{y_j}} + (n_j-y_j) \ln \frac{n_j (1-{\hat{\pi}}_j)}{n_j-y_j}= y_j \ln \frac{n_j {\hat{\pi}}_j}{{y_j}} + (n_j-y_j) \ln \frac{n_j-n_j {\hat{\pi}}_j}{n_j-y_j}
\end{aligned}

So now we can plug this into $-2\sum D_j$ to get

$$
D = -2 \sum_{j=1}^G  y_j \ln \frac{n_j {\hat{\pi}}_j}{{y_j}} + (1-y_j) \ln \frac{n_j-n_j {\hat{\pi}}_j}{n_j-y_j}
$$


</details>


**c)** Perform a likelihood ratio test - based on the deviance results given in the data frame in the R chunk- to compare the additive and interact models. First write down the null and alternative hypothesis you are testing. 
<span class="question">Which of the two models would you prefer?</span>
<details><summary>Hint</summary>
The models are 

additive: `Y ~ age + wantsMore`
interact: `Y ~ age * wantsMore`
</details>

<details><summary>Answer</summary>

The null hypothesis is $H_0$: `Y ~ age + wantsMore` is the correct model, i.e. all of the interaction terms are 0.
The alternative hypothesis is $H_1$: `Y ~ age * wantsMore` is the correct model.

We could do this using `anova(additive, interact)`, but let's do it by hand!

The deviances are (using `summary(models$additive)$deviance`):

additive: `R summary(models$additive)$deviance`
interact: `R summary(models$interact)$deviance`

with `R summary(models$additive)$df.residual` and `R summary(models$interact)$df.residual` degrees of freedom respectively (using `summary(models$additive)$df.residual`). The test statistic is then:

```{r}
(test.dev <- summary(models$additive)$deviance - summary(models$interact)$deviance)
(test.df <-  summary(models$additive)$df.residual - summary(models$interact)$df.residual)

pchisq(test.dev, test.df, lower.tail = FALSE)

```

Hence, we reject the null hypothesis, and prefer the model with the interactions.


</details>


<span class="question">**d)** What if you use the AIC for model selection, which model (out of the 5) will you then use?</span>

<details><summary>Hint</summary>
Lower AIC is better. 
</details>
<details><summary>Answer</summary>

We can extract the AICs from each object (or just read the output above). Then we can use R to find the model with the lowest AIC:

```{r}
(AICs <- unlist(lapply(models, AIC))) # lapply() returns a list, so unlist() tidies this up

which(AICs ==min(AICs))
```

And we find the interaction model, once again, is best.

</details>


## Problem 4: Plotting
(yes, mainly R and trying to understand) 

Finally, we want to use the additive model with all three covariates: `Y~age+education+wantsMore` to study different plots.

```{r,eval=TRUE}
fit3add=glm(cbind(using,notUsing)~age+education+wantsMore,data=ds,family=binomial)
summary(fit3add)
exp(fit3add$coefficients)
```


<span class="question">**a)** Comment on the output from `summary`. Use the deviance to test if this is a good model. Which model do you then compare to (what is the saturated model here)?.</span>

<details><summary>Hint</summary>
Think about both the directions of effects and their size.

You should know how a saturates model is defined. For extra marks, can you write the R formula for it?

</details>
<details><summary>Answer</summary>
The summary suggests a lot is going on. In particular:

- low education and wanting more children are associated with a lower probability of using contraceptives: the effect of wanting more children is much stronger (perhaps not surprisingly).
- the age levels are contrasts to "<25", so the effect of age 25-29 is to increade the odds $e^{0.39} = 1.48$ times.
- we can see that as age increases, the age effect also increases, i.e. the older women are, the more likely they are to use contraception.

We can test the goodness of the model by using the deviance. This is, as you now know, based on the difference between the likelihood for the model, and the likelihood for a saturated model. The saturated model has one parameter per data point (i.e. row in the data). Because the data are all factors, and we have every comnination, this is the same as the model with all possible interactions, i.e. `cbind(using,notUsing) ~ age*education*wantsMore`.

We can test if this is a good model from this line:

`Residual deviance:  29.917  on 10  degrees of freedom` 

We can compare this to a $\chi^2$ distribution:

```{r}
pchisq(29.917, 10, lower.tail = FALSE)
pchisq(fit3add$deviance, fit3add$df.residual, lower.tail = FALSE) < 0.05
```

The second line is just to show you how to (a) extract the statistics from the model, and (b) reduce hte p-value to a yes-no, i.e. if it is less than 0.05.

We can see that the test is significant, i.e. the model does not appear to be good.

</details>

**b)** Look at the plot of fitted values vs. deviance residuals. 
<span class="question">Do you see a trend?</span>

<details><summary>Answer</summary>
I don't. With 16 data points, random can look weird in all sorts of ways.
</details>


```{r}
library(ggplot2)
plotdf <- data.frame(dres = fit3add$residuals, fitted = fit3add$fitted.values, 
    age = ds$age)
ggplot(plotdf, aes(x = fitted, y = dres)) + geom_point() + labs(x = "Fitted values", 
    y = "Deviance residuals")
```

**c)** Here is a code to produce a plot of fitted values for the 16 covariate patterns together with a saturated model (called frac) - on the logit scale. 

```{r,eval=TRUE}
library(ggplot2)
frac=ds$using/(ds$using+ds$notUsing)
logitdf=data.frame(lfitted=qlogis(fit3add$fitted),lfrac=qlogis(frac),age=ds$age,wantsMore=ds$wantsMore,education=ds$education)
ggplot(logitdf,aes(x=age))+geom_point(aes(y=lfrac,colour="saturated"))+geom_point(aes(y=lfitted,colour="fitted"))+facet_wrap(facets=~wantsMore*education) + labs(x = "age", y = "")
```

<span class="question">How can you see from the dots for fitted values that we are assuming an additive model on the logit scale (that is, the pattern for each panel for logit fitted vs age is the same for each panel if we have an additive model on the logit scale - you see that from the fitted points)? Do you also see that for the observed values? Implication? Any other observations?</span>

<details><summary>Hint</summary>
If it helps, write Age as continuous, and ask how the line changes when factor levels change.

For the observed values, are there clear patterns? Note that this is a subjective judgement!
</details>
<details><summary>Answer</summary>
If the model is additive on the scale here, the plots should look similar: the relative effects of age should be the same, i.e. the effect of going from high to low education should be to change the age effect by the same amount. So, with some obvious notation we can go from 

$$
\eta_i = \beta_0 + \beta_{age25-29} + \beta_{educationHigh}
$$

to 

$$
\eta_i = \beta_0 + \beta_{age25-29} + \beta_{educationLow} = \beta_1 + \beta_{age25-29}
$$

and changing age changes $\beta_{age}$, but not $\beta_1$.

For the observed values this is not so clear. For `no, low` (i.e. not wanting more children and low education) the residuals are close to the fitted line. But for `no, high` the line seems steeper, and `yes, low` has a lower value for ages 40-59. 

The problem with interpreting these is that it is easy to over-interpret random noise (this is called [Pareidolia](https://en.wikipedia.org/wiki/Pareidolia)). But we know from the test of the residual deviance that this model does not fit the data well enough, so perhaps there is something going on. Or perhaps this extra variation is random noise (Robert Abelson [in his book](https://www.routledge.com/Statistics-As-Principled-Argument/Abelson/p/book/9780805805284?srsltid=AfmBOorgdpIjsLBH9mVBbOE2WUZ5JI3EKTQHbX9lTqNT7otypBln3TOK) summarised this as "chance is lumpy").
 
</details>

<!---
These links are dead. *sniff*

To see what Rodríguez have found to be the best model, look at the bottom of <http://data.princeton.edu/wws509/R/c3s6.html>.
  
**Note:** This exercise is based on the excellent notes of Germán Rodríguez at Princeton, see in particular the R logs at <http://data.princeton.edu/wws509/R/c3s1.html>.
--->

<!--
We will use a data set on embryogenic anters of the plant species Dutura innoxia (as described in Dobson and Barnett, 2008), given as
```{r}
storage=c(1,1,1,2,2,2)
centrifuge=c(40,150,350,40,150,350)
anters=data.frame(y=c(55,52,57,55,50,50),n=c(102,99,108,76,81,90),storage=factor(storage),lncentrifuge=log(centrifuge))
knitr::kable(anters)
```
For each combination of `storage` (1= control, 2=treatment=at 3 degrees for 48 hours) and `lncentrifuge`(centrifuging force, three continuous values chosen on the natural log scale) a total of `n` anthers (pollenknapper) were prepared and resulted in `y` embryogenic anthers (germinated). As covariates we will used `storage` as a factor (you may choose either dummy or effect coding), and `lncentrifuge` will be used as a continuous covariate.

The aim is to build a model to predict the probability of germination, and to understand the possible influence of covariates storage and centrifuge. 

1. Explore the data set. Plot each covariate against fraction of anters germinated for the different groups. Does the plot hint to a relationship between the possible covariates and the response?

2. Why may a binominal regression be suitable for these data? Write down assumptions we make if we want to fit a logit model.

We will continue to look at the logit model.

3. Let $x$ be the `lncentrifuge`. We start fitting Model 1: "linear effect of lncentrifuge, but different intercept and different slope for control and treatment". Hint: this can be given as model formula `cbind(y,n-y)~storage*lncentrifuge`. 
-->

<!-- ## Team mode kahoot! -->

<!-- [Beamer presentation with the 8 questions (answers on last slide))](https://www.math.ntnu.no/emner/TMA4315/2017h/M3quiz.pdf)  -->


# Exam questions

For this module the following are exam questions to work with

* 2012 – [Problem 1](https://www.math.ntnu.no/emner/TMA4315/Exam/GLMeksamen2012E.pdf)
* 2011 – [Problem 1](https://www.math.ntnu.no/emner/TMA4315/Exam/eksDes11e.pdf)
<!-- * 2009 – [Problem 3](https://www.math.ntnu.no/emner/TMA4315/Exam/e2009.pdf) -->
<!-- * 2008 – Problem 2 -->
<!-- * 2005 – Problem 1 and 2 -->

In addition these essay-type exam questions are closely related to this module.

## December 2014
There are two main asymptotic results that are behind essentially everything we
have done with respect to inference for generalized linear models. These are

1. asymptotic normality of maximum likelihood estimators (MLE), and

2. asymptotic result for likelihood ratio tests (LRT).

State, describe and discuss the main assumptions behind these two
asymptotic results, how these results are used in practice to do inference
for parameters, testing various hypotheses and comparing nested
models, for logistic regression.

## December 2016
We will consider (binomial or binary) logistic regression, where we have independent observations
$$Y_i \sim \text{Binomial}(n_i,\pi_i) \text{ } i=1,\ldots,n$$
so that 
$$P(Y_i)=\binom{n_i}{y_i}\pi_i^{y_i}(1-\pi_i)^{n_i-y_i}, \text{ } y_i=0,1,\ldots, n_i$$
The linear predictor is 
$$ \eta_i={\bf x}_i^T\boldsymbol{\beta}$$
and 
$$ \pi_i=\frac{\exp(\eta_i)}{1+\exp(\eta_i)}$$
or 
$$ \text{logit}(\pi_i)=\eta_i.$$
Here, ${\bf x}_i$ is a vector of the $p$ covariates for the $i$th observation $y_i$ with size (number of trials) $n_i$, and $\boldsymbol{\beta}$ is the vector of $p$ unknown regression coefficients.

Write an introduction to logistic 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 (when it is possible) and interpretation of the $\boldsymbol{\beta}$-coefficients
* Deviance and its usage.

# R packages

```{r, eval=FALSE}
install.packages(c("tidyverse",
                   "investr",
                   "knitr",
                   "kableExtra",
                   "faraway",
                   "viridis",
                   "statmod"))
```


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


