ggplot-glm-interactions.Rmd

---
title: "Customise a ggplot for a dataset with interactions between variables"
author: JR Ferrer-Paris (https://github.com/jrfep)
date: 2021-06-21
output:
  pdf_document: 
    fig_width: 4
    fig_height: 3
  html_document: default
---

I was given a dataset and a glm model with interactions, and the user wants to show several interactions in a single plot

## Libraries

First load the libraries

```{r}
require(ggplot2)
require(ggnewscale)
```

## Dataset

```{r data}
n.species <- c(20, 12, 8, 15, 9, 15, 15, 8, 6, 3, 12, 9, 19, 3, 2, 36, 13, 27, 13, 16, 15, 7, 12, 36, 30, 7, 17, 7, 3, 11, 17, 4, 15, 7, 12, 24, 11, 9, 22, 11, 22, 22, 17, 11, 14, 24, 14, 13, 18, 24, 4, 19, 17, 3, 10, 8, 16, 5, 40, 20)
##n.species <- c(4, 4, 9, 7, 7, 6, 11, 7, 18, 16, 9, 14, 2, 14, 15, 5, 3, 11, 4, 9, 8, 7, 6, 8, 7, 4, 4, 9, 8, 10, 15, 11, 8, 7, 9, 5, 9, 14, 4, 6, 9, 10, 13, 8, 9, 7, 7, 6, 12, 11, 6, 6, 5, 5, 8, 8, 10, 11, 11, 0)
log.area <- c(4.600786168, 4.600786168, 2.587377991, 0.690311286, 0.123721006, 1.475014166, 0.93399822, 1.725134595, 2.240236394, 0.611585246, 1.239139249, -0.117247844, 2.301647774, 1.820698544, 1.820698544, 1.657655225, 1.657655225, 2.341799021, 1.820698544, 1.350418601, 2.156923017, 2.156923017, 1.602848605, 0.549125927, 1.818928908, 1.790386511, 1.818928908, 0.991226076, 1.209675835, 1.457836498, 1.457836498, 2.996815091, 2.996815091, 2.996815091, 2.471614231, 2.471614231, 2.996815091, 2.996815091, 2.996815091, 1.347113302, 4.600786168, 0.951837592, 0.901534602, 0.509283189, 1.65374203, 1.65374203, 2.546547722, 2.546547722, 2.546547722, 2.546547722, 0.660751674, 0.594547233, 2.244172185, 2.244172185, 2.804492179, 4.235323083, 2.281941707, 2.281941707, 4.322433684, 4.322433684)
rain.diff <- c(98.3, 111.6, -172.8, 102.1, -126.8, -126.8, -29.8, 64.6, 72.5, 133.4, 66.9, 155.5, -119.5, -222.7, -222.7, 45.7, 45.7, 430.6, -40.7, 50.9, 61.9, -39.8, -31.8, 123.8, 130.8, -248.2, 130.8, 201.4, 170.4, 158.5, 160.1, 184.5, 212.3, 212.3, 429.6, 430, 131.1, 131.1, 131.1, 363.6, 274.1, 52.9, 62.9, 76.9, 6.1, 6.1, -75.9, -75.9, -81.9, -293.1, -278.8, -164.8, -3.3, -3.3, -3.3, -214.3, -188, -191.6, -18.8, -18.8)
blocks <- factor(c(2, 2, 2, 2, 1, 1, 1, 1, 3, 3, 1, 2, 3, 3, 3, 1, 1, 3, 3, 1, 1, 1, 2, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 3, 3, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 1, 3, 3, 3, 3),labels=c("G1","G2","G3"))

obs.data <- data.frame(n.species,log.area,rain.diff,blocks)

```

## Model

This model has many interactions, not all are important, but we want to visualise the differences:

```{r model}
mdl <- glm(n.species ~ (log.area*blocks) * (rain.diff*log.area), data=obs.data, family=quasipoisson())
summary(mdl)
```

## Plot

We have four variables that we want to display, three are continuous, one categorical. If we use the two axes for two variables, how do we show the other two?

Size and shape:

```{r plot-size-shape}
ggplot(data=obs.data) +
  geom_point(aes(x=log.area,y=n.species, size=rain.diff, shape=blocks))
```

Colour and shape:
```{r plot-colour-shape}
ggplot(data=obs.data) +
  geom_point(aes(x=log.area,y=n.species, colour=rain.diff, shape=blocks))
```

Or better size and colour and shape:
```{r plot-size-colour-shape}
ggplot(data=obs.data) +
  geom_point(aes(x=log.area,y=n.species, colour=rain.diff, size=rain.diff,shape=blocks))
```

That is not so bad, but we want to improve on this, you can save the plot and add elements to it:
```{r add-points}

(plot_with4vars <- ggplot(data=obs.data) +
  geom_point(aes(x=log.area,y=n.species, colour=rain.diff, size=rain.diff,shape=blocks)))

```

The previous plot is still too cluttered, we can divide it in facets:

```{r add-facet}
plot_with4vars + facet_grid(blocks~.)
```

This is more readable, but we want to join the legends of size and colour, and maybe we don't need the legend for blocks, since we already have them in different panels:

```{r unify-legends}
range(obs.data$rain.diff)
var3.limits <- c(-300,450)
var3.breaks <- seq(-300,450,by=150)

plot_with4vars + facet_grid(blocks~.) +
scale_size_continuous(limits=var3.limits, breaks=var3.breaks) +
scale_colour_continuous(limits=var3.limits, breaks=var3.breaks) +
guides(size=guide_legend(), colour=guide_legend(), shape='none')

```

We further split the plot by the sign of the difference in rainfall, left size are the negative values, right side the positive values.

To do this we need to add a new variable to the data set and we better repeat the call to the original plot (otherwise we will get warnings and errors for missing faceting variables)

```{r split-by-sign}
obs.data$rain.sign <- sign(obs.data$rain.diff)

plot_with4vars <- ggplot(data=obs.data) +
  geom_point(aes(x=log.area,y=n.species, colour=rain.diff, size=rain.diff,shape=blocks))

plot_with4vars + facet_grid(blocks~rain.sign) +
scale_size_continuous(limits=var3.limits, breaks=var3.breaks) +
scale_colour_continuous(limits=var3.limits, breaks=var3.breaks) +
guides(size=guide_legend(), colour=guide_legend(), shape='none')

```

Now we want to add the results of the model. We create a new data frame with new observations. We will use expand grid to create sequence of values along the range log.area and all levels of the blocks variables, but using only two quantiles of the rain.diff variable:

```{r prediction-grid}
pred.data <- expand.grid(log.area=seq(min(obs.data$log.area), max(obs.data$log.area), length=20), rain.diff=quantile(obs.data$rain.diff,c(.25,.75)), blocks=levels(obs.data$blocks))

pred.data$rain.sign <- factor(sign(pred.data$rain.diff))
```

We will add the predictions to this object. First predict the expected values and standard errors in the scale of the link function. Then, assuming assymptotic normality, we calculate the lower and upper intervals. Finally, we use the inverse link function to transform these values to the scale of the response variable.

```{r prediction}
prd <- predict(mdl,newdata=pred.data,type='link',se.fit=T)
ilink <- family(mdl)$linkinv

pred.data$prediction <- ilink(prd$fit)
pred.data$pred.lower <- ilink(prd$fit-1.96*prd$se.fit)
pred.data$pred.upper <- ilink(prd$fit+1.96*prd$se.fit)

```


Now we can add the lines and a ribbon with the confidence intervals
```{r add-lines}
plot_observed_data <- plot_with4vars + facet_grid(blocks~rain.sign) +
scale_size_continuous(limits=var3.limits, breaks=var3.breaks) +
scale_colour_continuous(limits=var3.limits, breaks=var3.breaks) +
guides(size=guide_legend(), colour=guide_legend(), shape='none')

plot_observed_data +
  geom_ribbon(data=pred.data, aes(x=log.area, ymin = pred.lower, ymax = pred.upper), alpha = 0.2, color = NA) +
      geom_line(data=pred.data,aes(x=log.area,y=prediction) )

```

Well that looks... awful!

We have very large confidence intervals (I warned you, the model is not optimal). We need to limit our y-axis to the range of the observed values:

```{r set-ylim}
plot_obs_and_pred_data <- plot_observed_data +
  geom_ribbon(data=pred.data, aes(x=log.area, ymin = pred.lower, ymax = pred.upper), alpha = 0.2, color = NA) +
      geom_line(data=pred.data,aes(x=log.area,y=prediction) )

plot_obs_and_pred_data + coord_cartesian(ylim = range(obs.data$n.species))
```

Well that looks better. It shows the estimated trend and the uncertainty of the model (very important!).

Maybe we can make this look even nicer with different colours of the predictions for negative and positive values of rain.diff:

```{r colour-sign, error=TRUE, plot=FALSE}

plot_observed_data +
  geom_ribbon(data=pred.data, aes(x=log.area, ymin = pred.lower, ymax = pred.upper, fill=rain.sign), alpha = 0.2, color = NA) +
      geom_line(data=pred.data,aes(x=log.area,y=prediction,colour=rain.sign) ) + coord_cartesian(ylim = range(obs.data$n.species)) + scale_fill_discrete(labels=c("min","max"))
```

Upps! that throws an error. We are mixing two different colour scales, one for the points and one for the line.

We can try to fix it declaring a new scale `new_scale_color`


```{r mix-colours}

plot_observed_data + new_scale_color() +
  geom_ribbon(data=pred.data, aes(x=log.area, ymin = pred.lower, ymax = pred.upper, fill=rain.sign), alpha = 0.2, color = NA) +
      geom_line(data=pred.data,aes(x=log.area,y=prediction,colour=rain.sign) ) + coord_cartesian(ylim = range(obs.data$n.species)) + scale_fill_discrete(labels=c("min","max"))
```

That looks... confusing!?. The plot is not that bad, but the legends are a mess. We get again separate size and colour legends for the points, and separate legends for lines and ribbon. You can spend hours trying to fix this, or maybe we should simplify again the plot and reduce the aesthetics that are redundant.

Let's get back to the beginning and simplify the observed point data first. Shapes are redundant with the facets, and colour and size are also redundant:

```{r simplify-plot}
plot_with4vars <- ggplot(data=obs.data) +
  geom_point(aes(x=log.area,y=n.species, size=rain.diff)) + facet_grid(blocks~rain.sign)

(plot_observed_data <- plot_with4vars  +
scale_size_continuous(limits=var3.limits, breaks=var3.breaks) )

```

We can now add the lines and ribbons without conflict with previous colour scales


```{r simple-lines-points}
(plot_obs_and_pred_data <-
  plot_observed_data +
    geom_ribbon(data=pred.data, aes(x=log.area, ymin = pred.lower, ymax = pred.upper, fill=rain.sign), alpha = 0.2, color = NA) +
      geom_line(data=pred.data,aes(x=log.area,y=prediction,colour=rain.sign) ) + coord_cartesian(ylim = range(obs.data$n.species)) + scale_fill_discrete(labels=c("min","max")) + scale_colour_discrete(labels=c("min","max")) )
```

This looks good, but we suddently realize the scale in the x axes is a logarithm of area, but we want to show the real value of area (in hectars). Should we go back and fix the data and model and run all the code again?... or maybe we just need to add one line of code to fix the plot:

```{r log-scale}
plot_obs_and_pred_data + scale_x_continuous( breaks=seq(0,4), labels=c(1,10,100,1000,10000))
```

***Important***: this fixes the looks in the plot, but the model is still using the log.area as a variable, this is ok for this example but might not be so in other contexts. Make sure you add accurate descriptions in methods and figure legends when you work in your research project.


Now the plot looks good enough, all the information we want is in there, and we discarded some redundant features.

As a final step we can fix the legends and labels to make them more informative, notice the use of several functions to add or overwrite labels:

```{r final}
plot_obs_and_pred_data +
  facet_grid(blocks~rain.sign, labeller = labeller(rain.sign = c(`-1`='Negative rainfall diff.',`1`='Positive rainfall diff.'))) +
  scale_x_continuous( breaks=seq(0,4), labels=c(1,10,100,1000,10000))  +
  labs(title='This plot is Amazing!',caption='Prediction from a GLM with Quasipoisson error distribution and log-link, confidence interval are approximate.',x='Area (ha) in log-scale',y='# of species') +
  labs(colour='Model prediction \n+ 95% C.I.', fill='Model prediction \n+ 95% C.I.',size='Difference\nin rainfall (mm)') +
  scale_colour_discrete(labels=c("with neg. rain diff.","with pos. rain diff.")) +
  scale_fill_discrete(labels=c("with neg. rain diff.","with pos. rain diff."))

```

## Useful links
- https://fromthebottomoftheheap.net/2017/05/01/glm-prediction-intervals-i/
- https://www.datanovia.com/en/blog/how-to-change-ggplot-facet-labels/


## Debugging info

Check the R session info:

```{r}
sessionInfo()
```