glm+.Rmd

---
title: "GLM with R" 
author: "D.-L. Couturier / R. Nicholls / C. Chilamakuri"
date: '`r format(Sys.time(), "Last modified: %d %b %Y")`'
output:
  html_document:
    theme: united 
    highlight: tango
    code_folding: show    
    toc: true           
    toc_depth: 2       
    toc_float: true     
    fig_width: 8
    fig_height: 6
---



<!--- rmarkdown::render("~/courses/cruk/LinearModelAndExtensions/git_linear-models-r/glm+.Rmd") --->
<!--- rmarkdown::render("/Volumes/Files/courses/cruk/LinearModelAndExtensions/git_linear-models-r/glm+.Rmd") --->

```{r message = FALSE, warning = FALSE, echo = FALSE} 
# change working directory: should be the directory containg the Markdown files:
# setwd("/Volumes/Files/courses/cruk/LinearModelAndExtensions/20200310/Practicals/")

# install gamlss package if needed
    # install.packages("gamlss")
```

# Section 1: Logistic regression

We will analyse the data collected by Jones (Unpublished BSc dissertation, University of Southampton, 1975). The aim of the study was to define if the probability of having Bronchitis is influenced by smoking and/or pollution.

The data are stored under data/Bronchitis.csv and contains information on 212 participants.


### Section 1.1: importation and descriptive analysis

Lets starts by

* importing the data set *Bronchitis* with the function `read.csv()`  
* displaying _bron_ (a dichotomous variable which equals 1 for participants having bronchitis and 0 otherwise) as a function of _cigs_, the number of cigarettes smoked daily.


```{r message = FALSE, warning = FALSE, echo = TRUE} 
Bronchitis = read.csv("data/Bronchitis.csv",header=TRUE)
plot(Bronchitis$cigs,Bronchitis$bron,col="blue4",
        ylab = "Absence/Presense of Bronchitis", xlab = "Daily number of cigarettes")
abline(h=c(0,1),col="light blue")
```

# Section 1.2: Model fit

Lets 

* fit a logistic model by means the function `glm()` and by means of the function `gamlss()` of the library `gamlss`.  
* display and analyse the results of the `glm` function : Use the function `summary()` to display the results of an R object of class `glm`.
```{r message = FALSE, warning = FALSE, echo = TRUE} 
fit.glm = glm(bron~cigs,data=Bronchitis,family=binomial)

library(gamlss)
fit.gamlss = gamlss(bron~cigs,data=Bronchitis,family=BI)

summary(fit.glm)
```

Let's now define the estimated probability of having bronchitis for any number of daily smoked cigarette and display the corresponding logistic curve on a plot:

```{r message = FALSE, warning = FALSE, echo = TRUE} 
plot(Bronchitis$cigs,Bronchitis$bron,col="blue4",
        ylab = "Absence/Presense of Bronchitis", xlab = "Daily number of cigarettes")
abline(h=c(0,1),col="light blue")

axe.x = seq(0,40,length=1000)
f.x = exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2])/(1+exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2]))
lines(axe.x,f.x,col="pink2",lwd=2)
```

## Section 1.3: Model selection

As for linear models, model selection may be done by means of the function `anova()` used on the glm object of interest. 

```{r message = FALSE, warning = FALSE, echo = TRUE} 
anova(fit.glm,test="LRT")
```

## Section 1.3: Model check

Lets assess is the model fit seems satisfactory by means 

* of the analysis of deviance residuals (function `plot()` on an object of class `glm`,
* of the analysis of randomised normalised quantile residuals (function `plot()` on an object of class `gamlss`,


```{r message = FALSE, warning = FALSE, echo = TRUE}
# deviance
par(mfrow=c(2,2),mar=c(3,5,3,0))
plot(fit.glm)
# randomised normalised quantile residuals
plot(gamlss(bron~cigs,data=Bronchitis,family=BI))
```

## Section 1.4: Fun

```{r message = FALSE, warning = FALSE, echo = TRUE}
# long format:
long = data.frame(mi = rep(c("MI","No MI"),c(104+189,11037+11034)),
                  treatment = rep(c("Aspirin","Placebo","Aspirin","Placebo"),c(104,189,11037,11034)))
# short format: 2 by 2 table
table2by2 = table(long$treatment,long$mi)

#
chisq.test(table2by2)
prop.test(table2by2[,"MI"],apply(table2by2,1,sum))
summary(glm(I(mi=="MI")~treatment,data=long,family="binomial"))
```


# Section 2: Poisson regression

The dataset *students.csv* shows the number of high school students diagnosed with an infectious disease for each day from the initial disease outbreak. 

# Section 2.2: Importation

Lets 

* import the dataset by means of the function `read.csv()`
* display the daily number of students diagnosed with the disease (variable `cases`) as a function of the days since the outbreak (variable `day`). 

```{r message = FALSE, warning = FALSE, echo = TRUE} 
students = read.csv("data/students.csv",header=TRUE)
plot(students$day,students$cases,col="blue4",
        ylab = "Number of diagnosed students", xlab = "Days since initial outbreak")
abline(h=c(0),col="light blue")
```

# Section 2.2: Model fit

Lets 

* fit a poisson model by means the function `glm()` and by means of the function `gamlss()` of the library `gamlss`.  
* display and analyse the results of the `glm` function : Use the function `summary()` to display the results of an R object of class `glm`.

```{r message = FALSE, warning = FALSE, echo = TRUE} 
fit.glm = glm(cases~day,data=students,family=poisson)

library(gamlss)
fit.gamlss = gamlss(cases~day,data=students,,family=PO)

summary(fit.glm)
```

Let's now define the estimated probability of having bronchitis for any number of daily smoked cigarette and display the corresponding logistic curve on a plot:

```{r message = FALSE, warning = FALSE, echo = TRUE} 
plot(students$day,students$cases,col="blue4",
        ylab = "Number of diagnosed students", xlab = "Days since initial outbreak")
abline(h=c(0),col="light blue")

axe.x = seq(0,120,length=1000)
f.x = exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2])
lines(axe.x,f.x,col="pink2",lwd=2)
```

## Section 2.3: Model selection

As for linear models, model selection may be done by means of the function `anova()` used on the glm object of interest. 

```{r message = FALSE, warning = FALSE, echo = TRUE} 
anova(fit.glm,test="LRT")
```

## Section 2.3: Model check

Lets assess is the model fit seems satisfactory by means 

* of the analysis of deviance residuals (function `plot()` on an object of class `glm`,
* of the analysis of randomised normalised quantile residuals (function `plot()` on an object of class `gamlss`,


```{r message = FALSE, warning = FALSE, echo = TRUE}
# deviance
par(mfrow=c(2,2),mar=c(3,5,3,0))
plot(fit.glm)
# randomised normalised quantile residuals
plot(fit.gamlss)
```



# Section 6: Practicals


### (i) *Bronchitis.csv*
Analyse further the Bronchitis data of Jones (1975) by 

* first investigating if the probability of having bronchitis also depends on _pollution_ (variable `poll`),
* second investigating if there is an interaction between the variables `cigs` and `poll`.


Lets plot the data first. 
```{r message = FALSE, warning = FALSE, echo = TRUE} 
Bronchitis = read.csv("data/Bronchitis.csv",header=TRUE)
# plot
plot(Bronchitis$poll,Bronchitis$bron,col="blue4",
        ylab = "Absence/Presense of Bronchitis", xlab = "Pollution level")
abline(h=c(0,1),col="light blue")
```
No obvious relationship between pollution and bronchitis is visible by means of this plot.

Lets fit a model assuming that the probability of getting bronchitis is a function of the pollution level.

```{r message = FALSE, warning = FALSE, echo = TRUE} 
# fit1:
fit.glm = glm(bron~poll,data=Bronchitis,family=binomial)
summary(fit.glm)
```

The intercept of the previous fit allows to define the probability of getting bronchitis when the level of pollution equals 0.
As a zero level of pollution is (i) out of range (ii) not a realistic value, we will 

* create the variable `poll_centered` defined as the pollution level minus the mean so that the intercept corresponds to the probability of getting bronchitis for an average pollution level in Cardiff,
* refit the model

```{r message = FALSE, warning = FALSE, echo = TRUE} 
# fit2:
Bronchitis$poll_centered = Bronchitis$poll-mean(Bronchitis$poll)
fit.glm = glm(bron~poll_centered,data=Bronchitis,family=binomial)
library(gamlss)
fit.gamlss = gamlss(bron~poll_centered,data=Bronchitis,family=BI)
summary(fit.glm)
```

Lets 

* perform a model check but plotting the randomised quantile residuals of gamlss a few times
* plot the fitted probabilities

```{r message = FALSE, warning = FALSE, echo = TRUE}
# model check
plot(gamlss(bron~poll_centered,data=Bronchitis,family=BI))
plot(gamlss(bron~poll_centered,data=Bronchitis,family=BI))
plot(gamlss(bron~poll_centered,data=Bronchitis,family=BI))

# plot fit
plot(Bronchitis$poll_centered,Bronchitis$bron,col="blue4",
        ylab = "Absence/Presense of Bronchitis", xlab = "Pollution level")
abline(h=c(0,1),col="light blue")
axe.x = seq(-10,10,length=1000)
f.x = exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2])/(1+exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2]))
lines(axe.x,f.x,col="pink2",lwd=2)
```
Model check suggests a good fit. Lets finally check if the interaction is significant:
    

```{r message = FALSE, warning = FALSE, echo = TRUE} 
# interaction ?
fit.glm = glm(bron~poll_centered*cigs,data=Bronchitis,family=binomial)
summary(fit.glm)
anova(fit.glm,test="LRT")
```
Interaction is not significant.



### (ii) *myocardialinfarction.csv*

The file *myocardialinfarction.csv* indicates if a participant had a myocardial infarction attack (variable `infarction`) as well the participant's treatment (variable `treatment`).  

Does _Aspirin_ decrease the probability to have a myocardial infarction attack ?

Lets (i) import the dataset, (ii) change the levels of the factor `treatment` so that `Placebo` corresponds to the reference group, (iii) and finally plot the (sample) probabilities to get an attack by treatment group
```{r message = FALSE, warning = FALSE, echo = TRUE} 
# import
myocardialinfarction = read.csv("data/myocardialinfarction.csv")
# by default, Aspirin is the reference group as the alphabetic order is used
myocardialinfarction$treatment = factor(myocardialinfarction$treatment,
                                         levels=c("Placebo","Aspirin"))
# plot
par(mfrow=c(1,1),mar=c(3,4,3,1))
pi.group =  tapply(myocardialinfarction$infarction=="attack",myocardialinfarction$treatment,mean)
table.group = tapply(myocardialinfarction$infarction=="attack",myocardialinfarction$treatment,table)
temp = barplot(pi.group,plot=FALSE)
barplot(pi.group,ylab="Probability",xlab="",
        main = "Probability of myocardial infarction\n per treatment group",names=rep("",2),
        cex.axis=.6,axes=FALSE,cex.main=1.4)
axis(2,las=2,cex.axis=.8)
axis(1,temp[,1],names(pi.group),cex.axis=1.25,tick=FALSE)
for(gw in 1:length(pi.group)){
    text(temp[gw],pi.group[gw]/2,.p(table.group[[gw]]["TRUE"]," / ",sum(table.group[[gw]])),
         col="red",cex=1.5)
    }   
```
The barplot seems to suggest that the treatment (aspirin) reduces the risk of myocardial infarction. Lets fit a logistic model to assess if this difference is significant. 
Note that, in this case (a dichotomous outcome and a dichotomous predictor), a test of equality of proportions or an independence test could also do the job. 
With a logistic model, other predictors could easily be added to the model and the beta parameter corresponding to the treatment can be interpreted by means of odd ratios (or relative risk ratios when prevalences are *small*, as we will note at the end of this practical).


```{r message = FALSE, warning = FALSE, echo = TRUE} 
# model fit
fit.glm = glm(I(infarction=="attack")~treatment,data=myocardialinfarction,family=binomial)
summary(fit.glm)
# test of equality of (independent) proportions
prop.test(unlist(lapply(table.group,function(x)x[2])),unlist(lapply(table.group,sum)))
# test of independence
chisq.test(matrix(unlist(table.group),ncol=2))
```
The three methods lead to the same conclusion: there is a significant difference between the probabilities of having a myocardial infarction of the two treatment groups.
Note that the two last methods get exactly the same results (as they use the same test X-squared statsitic).

Lets define the fitted probabilities to get an attack and compare them to the sample probabilities: they should match:

```{r message = FALSE, warning = FALSE, echo = TRUE} 
# pi_placebo:
pi_placebo = exp(-4.04971)/(1+exp(-4.04971))
pi_placebo
pi.group[1]
# pi_aspirin
pi_aspirin = exp( -4.04971-0.60544)/(1+exp( -4.04971-0.60544))
pi_aspirin
pi.group[2]
```

Finally, note that when prevalence are small, the exponential of the logistic regression corresponding to the treatment **may** also be interpreted as relative risk ratios. Indeed:  
```{r message = FALSE, warning = FALSE, echo = TRUE} 
# interpreation of exp(beta1) when prevalences are small
pi_aspirin/pi_placebo
exp(-0.60544)
c(exp(-0.60544-qnorm(.975)*0.12284),exp(-0.60544+qnorm(.975)*0.12284))
```
Thus, aspirin strongly reduces the risk of myocardial infarction.




### (ii) *crabs.csv*

This data set is derived from Agresti (2007, Table 3.2, pp.76-77). It gives 6 variables for each of 173 female horseshoe crabs:

* Explanatory variables that are thought to affect this included the female crab’s color (C), spine condition (S), weightweight (Wt)
* C: the crab's colour,
* S: the crab's spine condition, 
* Wt: the crab's weight,
* W: the crab's carapace width, 
* Sa: the response outcome, i.e., the number of satellites. 

Check if the width of female's back can explain the number of satellites attached by fitting a Poisson regression model with width.


Lets import the datasset, fit a poisson loglinear model, plot the fit and perfom a model check :
```{r message = FALSE, warning = FALSE, echo = TRUE} 
crabs = read.csv("data/crab.csv",header=TRUE)

# plot:
plot(crabs$W,crabs$Sa,col="blue4",
        ylab = "Number of satellites", xlab = "width of female's back")
abline(h=c(0),col="light blue")


# fit
fit.glm = glm(Sa~W,data=crabs,family=poisson)
library(gamlss)


# plot fit:
plot(crabs$W,crabs$Sa,col="blue4",
        ylab = "Number of satellites", xlab = "width of female's back")
abline(h=c(0),col="light blue")

axe.x = seq(15,40,length=1000)
f.x = exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2])
lines(axe.x,f.x,col="pink2",lwd=2)
    # not a great fit...

# model check
plot(gamlss(Sa~W,data=crabs,family=PO))
plot(gamlss(Sa~W,data=crabs,family=PO))
plot(gamlss(Sa~W,data=crabs,family=PO))

    # confirm lack of fit -> bin the estimates
    
# 2 alternative models
plot(gamlss(Sa~W,data=crabs,family=ZIP))
plot(gamlss(Sa~W,data=crabs,family=NBI))
    # check ?ZIP and ?NBI for detail       
```
Reasonably, there is a lack of fit> the estimates are not to be trusted.