---
title: "Generalized Linear Models"
author: "Mathieu RIBATET"
format: html
editor: visual
---

## 📝 Newton-Raphson

In this exercise, we will play a bit with the Newton-Raphson and Fisher's scoring algorithms. Althgouh this exercise is a bit theoretical, we will end with our own software implementation.

Let $f\colon \mathbb{R} \to \mathbb{R}$ be a smooth enough function, e.g., differentiable and co... For short, if you need a specific property for $f$, it is satisfied! (Yes, I'm a bit harsh here but I don't want to bother you with technical details here since it is not a lecture about optimization). We aim at finding a root $x_*$ of $f(x) = 0$.

1.  Let $x_* = x_0 + h$ (\$h\$ 'small' obviously). Give a Taylor expansion (to the 1st order) of $f$ around $x=x_*$. By the way, another name for this, is to give a tangent line approximation.

2.  Using the above approximation, give an approximation for $h$ and hence a 'better' approximation for $x_*$. (We will obviously assume that $f'(x_0) \neq 0$)

3.  Write a pseudo code that iterate this procedure... this is the Newton-Raphson algorithm! For those interested in the analysis of the convergence of this algorithm, Wikipedia is a good resource. Although very interesting, it is beyond the scope of this lecture.

4.  Write a Newton-Raphson algorithm to solve the following equations

    $$
    x^3 - 2 x - 5 = 0, \qquad \cos x = x^3,
    $$

    where you will use as 'initial guess' $x_0 = 2$ and $x_0 = 0.5$ respectively. (You can play a bit with these starting values if you feel adventurous 😉)

5.  Let $\{f(x; \theta)\colon x \in \mathbb{R}, \theta \in \Theta\}$ be a regular parametric statistical model. Given a sample $X_1, \ldots, X_n$, write down the score equation. Write a pseudo code for the Newton-Raphson algorithm in this context. What is $-f'$ in this situation?

6.  Application. Assume that our statistical model is the Poisson model and that we observed

    $$
    X_1 = 1, X_2 = 5, X_3 = 5, X_4 = 5, X_5 = 4, X_6 = 6, X_7 = 2, X_8 = 4, X_9 = 2, X_{10} = 5.
    $$

    Give a numerical approximation of the maximum likelihood estimator and compare it to its closed form solution.

7.  Try another implementation where you substitute $f'$ with its expected value. What is the name of this algorithm?

## 📝 The glm function

The objective of this exercise is to get you familiarized with the *glm* function. Of course, you should by now be familiar with the theory of GLMs. If not, reread the lecture notes and/or ask questions!

1.  Read the documentation of the *glm* function. Pay particular attention to the *family* argument. You probably have to look at the *family* function too.

2.  Run the following piece of code and describe what model is fitted here. What is the link function used here?

    ```{r}
    counts <- c(18,17,15,20,10,20,25,13,12)
    outcome <- gl(3,1,9)
    treatment <- gl(3,3)
    df <- data.frame(treatment, outcome, counts)
    fit <- glm(counts ~ outcome + treatment, family = poisson, data = df)
    ```

3.  Comment the output of the following code

    ```{r}
    summary(fit)
    ```

4.  Is the *`treatment`* variable statistically significant? What about *`outcome`*?

5.  Read the documentation of the *`predict`* (well actually *`predict.glm`*) function. Make sure you understand all the different types of prediction, i.e., when you set *`type`* to either `link`, `response` or `terms`.

## 📝 Clotting times

For this exercise we are interested in modelling clotting times of blood for normal plasma diluted to nine percentage concentrations with prothrombin-free plasma; clotting was induced by two lots of thromboplastin.

1.  Read [Section 8.4.2](https://mribatet.perso.math.cnrs.fr/CentraleNantes/SDR/ClottingTimes.pdf) of the McCullagh and Nelder's book that talks about modelling clotting times of blood and make sure you understand it!

2.  Do the initial plot mentioned in the text and comment.

3.  Using R, fit the models described at page 301. Are your results in agreement with what the book says?

4.  The text says "Clearly al the terms are necessary". How would you assert that the above statement is correct?

5.  Do the plots described in the paragraph starting with "The plot of the Pearson residuals..."

6.  Do you understand why the authors state that for a GLM with the Gamma family we have $\mbox{Var}(Y) \propto \mu^2$ while when fitting a usual linear model on $1/Y$ we have $\mbox{Var}(Y) \propto \mu^4$?

## 📝 Oesophageal cancer in Ille-et-Vilaine

For this exercise we will analyse the dataset *`esoph`* included directly within the R software.

1.  Familiarize yourself with the dataset.

2.  Run the following piece of code and explain what model is fitted here. Why in the R formula the left hand side has two entries *`ncases`* and *`ncontrols`*?

    ```{r}
    fit <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp + alcgp, data = esoph, family = binomial())
    ```

    **Give comments here.**

3.  Perform a residual analysis.

    ```{r}
    ## Insert code here
    ```

    **Give comments here.**

4.  Interpret results.

    ```{r}
    ## Insert code here
    ```

    **Give comments here.**

5.  Fit a `quasibinomial` GLM and spot any differences between this new model and the previous one.

    ```{r}
    ## Insert code here
    ```

    **Give comments here.**
