# Solution to GLM exercise


```{r setuppp, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,warning=FALSE,message = FALSE)
library(tidyverse)
```

```{r}
theme_set(theme_bw())
d<-read.csv("/home/aqm/course/data/HedisteCounts.csv")

```

## Visualise the data

The best way to look at these data quickly is by plotting histograms of the counts. I'll use ggplot and facet wrapping

```{r}
ggplot(d, aes(x=Count)) +geom_histogram(binwidth=3,fill="grey", col="black") +facet_wrap(~Substrate)
```

## Simple non parametric

So its very obvious that the data are **not** normally distributed and there are lots of zeros. So, an undergraduate would probably use a Wilcoxon rank sum test (Mann Whitney in SPSS) in order to get a p-value.

```{r}
wilcox.test(d$Count~d$Substrate)
```
Great! We've found a significant difference. However, where is the effect size?

Let's get the means anyway.


```{r}
d%>% group_by(Substrate) %>% summarise(mean(Count))

```
Looks like there are around four times as many ragworm in the mud than the sand. But we can't say that based on the non-parametric test.

## T test

Let's do a T test

```{r}
t.test(d$Count~d$Substrate)
```

Yippee! It's significant again. But my tutor is gong to tell me to put it into SPSS and SPSS will tell me that the test isn't valid.

## One way Anova

```{r}
mod<-lm(data=d,Count~Substrate)
anova(mod)
```

Test residuals for normality

```{r}
shapiro.test(residuals(mod))
```

Sugar. Not good. 

Test for homogeneity of variance?

```{r}
bartlett.test(residuals(mod)~d$Substrate)
```

Not valid. So, I'll have to use the wilcox test? 

## Poisson regression

Now I'm a master's student I'm more sophisticated. I'll use a GLM of the poisson family with a log link function. That will convince my tutor!


```{r}
mod<-glm(data=d,Count~Substrate,family="poisson")
summary(mod)
```

Oops. I've spotted over dispersion. Should have seen that one coming from the histograms.

## Negative binomial

But... I'm a very savvy MSC student because I took the AQM course. I have an answer to that. I'll use a glm of the negative binomial family.

```{r}
library(MASS)
mod<-glm.nb(data=d,Count~Substrate)
summary(mod)
```

That looks much better. 

Let's try to understand the coefficients.

```{r}
exp(coef(mod))

```

Aha. The first is the mean for mud. What is the second one?

```{r}
exp(coef(mod)[1]) * exp(coef(mod)[2])
```

Right. It's the value that I need to multiply the count in mud by in order to get the count in sand.

Effect size with confidence intervals?

```{r}
exp(confint(mod))
```


So the effect lies between 12% and 58%. We have an answer!




## Presence absence

Converting to presence absence.

```{r}
d$pres<-1*(d$Count>0)
```


Tabulate the data.

```{r}
table(d$pres,d$Substrate)
```




Plot as bar chart. (An undergrad would use Excel)

```{r}
d %>% group_by(Substrate,pres) %>% summarise(n=n()) %>%
  ggplot(aes(x=Substrate,y=n, fill=as.factor(pres))) + geom_col(position="dodge") 
```


Chi square test

It’s a 2x2 contingency table. I can do a chi squared!

```{r}
chisq.test( table(d$pres,d$Substrate))
```


Yippee! It’s significant. But … where is the effect size?
Fit a binomial glm

```{r}
mod<-glm(data=d,pres~Substrate,family="binomial")
summary(mod)
```



```{r}
exp(coef(mod))
```



Hang on… how do we interpret these exponentiated coeficients? They are odds.

```{r}
table(d$pres,d$Substrate)
```



So 27 to 23 for presence in mud gives ..

```{r}
27/23
```



Multiply this by the “relative risk” of finding ragworm in sand.

```{r}
1.173913 * 0.3097643 
```




Should be the same as

```{r}
16/44
```



Which it is.

So when we get confidence intervals the key value refers to the odds ratio.

```{r}
exp(confint(mod))
```



So we have got a confidence interval for a meaningful quantity, but we have unfortunately entered the interpretabilty minefield thrown up by relative odds ratios. Still better than using a chi squared!


We end up by being able to state that the relative odds of finding ragworm present in “sand” cores compared to “mud” cores (note that the variable itself is based on a binary split) is between 14% and 68%. Compare this to the difference between the result for the analysis of abundance. It is surprisingly similar. This illustrates that there is not really that much information in the actual counts. Presence vs absence is just about as good a measure (statistically). Of course, if you want to make an estimate of the overall number of ragworm then you do need to count them.

Note that if you have been following stories on vaccine trials in the news, this is equivalent to the stats being reported for a RCT (randomised control trial) between vaccinated and placebo groups. This is because analyses of drug trials also tend to use logistic regression (i.e glms of the binomial family)


## Chi squared test and poisson glms

As an aside, I mentioned that the chi squared test is essentially just a poisson glm.
This is how you can use a poisson glm to get the same p-value. The input data are the tabulated counts.

```{r}
d %>% group_by(pres,Substrate) %>% summarise(count=n()) -> tb
#tb<-as.data.frame(table(d$pres,d$Substrate)) ## base R tables need to be turned to df
#names(tb)<-c("pres","substrate","count")
tb$pres<-as.factor(tb$pres)
tb
```


If you remember your introductory stats the chisquared test using a contingency table is a test for association. So the element that is of interest in the model is an interaction between presence and substrate.

```{r}
mod<-glm(data=tb,count~pres*Substrate,family="poisson")
summary(mod)
```


The p value should be the same as the chi squared test, if the technical continuity correction used for contingency tables is not made.

```{r}
chisq.test( table(d$pres,d$Substrate),correct = FALSE)
```



Which it is. There is a very minor numerical difference due to rounding and minor differences in the actual mechanism of calculation, but the two numbers are practically equivalent. This equivalence applies to all situations, not just this data set. Notice it is also the same as the p-value for the binomial model.

The point being illustrated is that the assumption made by the chi squared test is that the counts follow a poisson distribution with a variance equal to the mean. This is used when testing for a difference between the observed number of counts and the expected number of counts under the null hypothesis of no association (interaction) between presence and substrate. So, the chi squared test is not even really a non parametric test!

Of the three possible ways of getting a p-value the most useful is clearly the binomial glm, as this provides confidence intervals. The poisson glm for the counts actually also could provide the same intervals, by looking at the interaction term, but this is rarely used in this context.


```{r}
exp(confint(mod))
```






