---
title: "Quick overview of plot functions"
author: "Jacolien van Rij"
date: "`r Sys.Date()`"
bibliography: bibliography.bib
output:   
  rmarkdown::html_document:
    theme: readable
    highlight: default
vignette: >
  %\VignetteIndexEntry{Quick overview of plot functions}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---
<style>
samp    {font-style:900;}
.i      {color:cornflowerblue;}
th      {
  border-bottom-style: solid; 
  border-bottom-width:2pt; 
  border-top-style: solid; 
  border-top-width:2pt;
  padding: 5pt;
}
td      {
  border-bottom-style: solid; 
  border-bottom-width:1pt; 
  padding:5pt;
}
.sideborder {
  border-right-style:solid; 
  border-right-width: 1pt;
}
.noborder {
  border-right-style:solid; 
  border-right-width: 1pt;
  border-bottom-width:0pt; 
}

</style>
```{r, include=FALSE}
library(itsadug)
infoMessages("off")
```

<p></br></p>
The table present the different plot functions in the packages **`mgcv`** [@Wood_2006; @Wood_2011] and <span style="color:cornflowerblue">**`itsadug`**</span> for visualizing GAMM models.

<p></br></p>
<table>
<tr> 
<th align="left" class="sideborder"> </th> 
<th align="left" class="sideborder"> Partial effect </th> 
<th align="left" class="sideborder"> Sum of all effects </th> 
<th align="left"> Sum of "fixed" effects<sup>1</sup> </th> 
</tr>
<tr> 
<td align="left" rowspan="2" class="sideborder"> surface </td> 
<td align="center" class="noborder" style="border-right:solid 1pt;"> <samp>plot.gam()</samp></td> 
<td align="center" class="noborder"> <samp>vis.gam()</samp></td> 
<td align="center" style="border-width:0pt;" bgcolor="lightgray"> </td> 
</tr>
<tr> 
<td align="center" class="sideborder"> <samp class="i"><a href="#pvisgam">pvisgam()</a></samp> </td> 
<td align="center" colspan="2"> <samp class="i"><a href="#fvisgam">fvisgam()</a></samp> </td> 
</tr>
<tr> 
<td align="left" class="sideborder"> smooth </td> 
<td align="center" class="sideborder"> <a href="#plotgam1"><samp style="color:black;">plot.gam()</samp></a></td> 
<td align="center" colspan="2"> <samp class="i"><a href="#plotsmooth">plot_smooth()</a></samp></td> 
</tr>
<tr> 
<td align="left" class="sideborder"> group estimates </td> 
<td align="center" class="sideborder"><a href="#plotgam2"><samp style="color:black;">plot.gam()<sup>2</sup></samp></a></td> 
<td align="center" colspan="2"> <samp class="i"><a href="#plotparametric">plot_parametric()</a></samp></td> 
</tr>
<tr> 
<td align="left" class="sideborder"> random smooths </td> 
<td align="center" class="sideborder"> <samp class="i"><a href="#getrandom">get_random()</a></samp>, <samp class="i"><a href="#inspectrandom">inspect_random()</a></samp></td> 
<td align="center" colspan="2"> </td> 
</tr>
</table>

<sup>1</sup>: include `rm.ranef=TRUE` to zero all random effects.

<sup>2</sup>: include `all.terms=TRUE` to visualize parametric terms.

<p></br></p>
# Examples

```{r}
library(itsadug)
library(mgcv)
data(simdat)

## Not run: 
# Model with random effect and interactions:
m1 <- bam(Y ~ Group + te(Time, Trial, by=Group)
    + s(Time, Subject, bs='fs', m=1, k=5),
    data=simdat)
# Simple model with smooth:
m2 <-  bam(Y ~ Group + s(Time, by=Group)
    + s(Subject, bs='re')
    + s(Subject, Time, bs='re'),
    data=simdat)
```

Summary model `m1`:
```{r, results='asis'}
gamtabs(m1, type='html')
```

Summary model `m2`:
```{r, results='asis'}
gamtabs(m2, type='html')
```


## a. Surfaces

Plotting the partial effects of `te(Time,Trial):GroupAdults` and `te(Time,Trial):GroupChildren`. 


### <a id="pvisgam"></a> pvisgam()

```{r, fig.width=8, fig.height=4}
par(mfrow=c(1,2), cex=1.1)

pvisgam(m1, view=c("Time", "Trial"), select=1,
        main="Children", zlim=c(-12,12))
pvisgam(m1, view=c("Time", "Trial"), select=2,
        main="Adults", zlim=c(-12,12))
```

Notes:

- Plots same data as `plot(m1, select=1)`: *partial effects plot*, i.e., the plot does not include intercept or any other effects.

- Make sure to set the zlim values the same when comparing surfaces

- Use the argument `cond` to specify the value of other predictors in a more complex interaction. For example, for plotting a modelterm `te(A,B,C)` use something like `pvisgam(model, view=c("A", "B"), select=1, cond=list(C=5))`.


### <a id="fvisgam"></a> fvisgam()


Plotting the fitted effects of `te(Time,Trial):GroupAdults` and `te(Time,Trial):GroupChildren`. 


```{r, fig.width=8, fig.height=4}
par(mfrow=c(1,2), cex=1.1)

fvisgam(m1, view=c("Time", "Trial"), cond=list(Group="Children"),
        main="Children", zlim=c(-12,12), rm.ranef=TRUE)
fvisgam(m1, view=c("Time", "Trial"), cond=list(Group="Adults"),
        main="Adults", zlim=c(-12,12), rm.ranef=TRUE)

```

Notes:

- Plots the *fitted effects*, i.e., the plot includes intercept and estimates for the other modelterms. 

- Make sure to set the zlim values the same when comparing surfaces

- Use the argument `cond` to specify the value of other predictors in the model. 

- The argument `rm.ranef` cancels the contribution of random effects.

- The argument `transform` accepts a function for transforming the fitted values into the original scale.


## b. Smooths

### <a id="plotgam1"></a> plot.gam()

Plotting the partial effects of `s(Time):GroupAdults` and `s(Time):GroupChildren`. 

```{r, fig.width=8, fig.height=4}
par(mfrow=c(1,2), cex=1.1)

plot(m2, select=1, shade=TRUE, rug=FALSE, ylim=c(-15,10))
abline(h=0)
plot(m2, select=2, shade=TRUE, rug=FALSE, ylim=c(-15,10))
abline(h=0)
```

Notes:

- Not possible to combine different smooth terms in a single plot.

Alternatively:

```{r, fig.width=4, fig.height=4}
par(mfrow=c(1,1), cex=1.1)

# Get model term data:
st1 <- get_modelterm(m2, select=1)
st2 <- get_modelterm(m2, select=2)

# plot model terms:
emptyPlot(2000, c(-15,10), h=0,
    main='s(Time)', 
    xmark = TRUE, ymark = TRUE, las=1)
plot_error(st1$Time, st1$fit, st1$se.fit, shade=TRUE)
plot_error(st2$Time, st2$fit, st2$se.fit, shade=TRUE, col='red', lty=4, lwd=2)

# add legend:
legend('bottomleft',
  legend=c('Children', 'Adults'),
  fill=c(alpha('black'), alpha('red')),
  bty='n')
```



### <a id="plotsmooth1"></a> plot_smooth()


Plotting the fitted effects of `te(Time,Trial):GroupAdults` and `te(Time,Trial):GroupChildren` i.e., the plot includes intercept and estimates for the other modelterms.  


```{r, fig.width=8, fig.height=4}
par(mfrow=c(1,2), cex=1.1)

plot_smooth(m1, view="Time", cond=list(Group="Children"),
    rm.ranef=TRUE, ylim=c(-6,10))
plot_smooth(m1, view="Time", cond=list(Group="Adults"),
    col="red", rug=FALSE, add=TRUE,
    rm.ranef=TRUE)

# or alternatively:
plot_smooth(m1, view="Time", plot_all="Group",
    rm.ranef=TRUE)
```

Notes:

- Use the argument `cond` to specify the value of other predictors in the model. 

- The argument `rm.ranef` cancels the contribution of random effects.

- The argument `transform` accepts a function for transforming the fitted values into the original scale.

- The argument `plot_all` plots all levels for the given predictor(s).


## c. Group estimates

### <a id="plotgam2"></a> plot.gam()

Plotting the partial effect of grouping predictors such as `Group`:

```{r, fig.width=4, fig.height=4}
par(mfrow=c(1,1), cex=1.1)

plot.gam(m1, select=4, all.terms=TRUE, rug=FALSE)
```


Alternatively, use `get_coefs()` to extract the coefficients and plot these:

```{r, fig.width=8, fig.height=4, results='hold'}
coefs <- get_coefs(m1)
coefs

par(mfrow=c(1,2), cex=1.1)

b <- barplot(coefs[,1], beside=TRUE, 
             main="Parametric terms",
             ylim=c(0,5))
errorBars(b, coefs[,1], coefs[,2], xpd=TRUE)

# Note that the effect of Group is a *difference* estimate
# between intercept (=GroupChildren) and Group Adults

b2 <- barplot(coefs[1,1], beside=TRUE, 
             main="Estimate for Group",
             ylim=c(0,5), xlim=c(0.1,2.5))
mtext(row.names(coefs), at=b, side=1, line=1)
abline(h=coefs[1,1], lty=2)
rect(b[2]-.4, coefs[1,1], b[2]+.4, coefs[1,1]+coefs[2,1],
     col='gray')
errorBars(b, coefs[,1]+c(0,coefs[1,1]), coefs[,2], xpd=TRUE)
```

Notes:

- For large models `get_coefs()` is faster than `summary(model)`.


### <a id="plotparametric"></a> plot_parametric()

Plotting the fitted effects of grouping predictors such as `Group`:

```{r, fig.width=4, fig.height=4}
pp <- plot_parametric(m1, pred=list(Group=c("Children", "Adults")) )
pp
```


## Random effects

### <a id="getrandom"></a>get_random()

For extracting the random effects coefficients (random adjustments of intercept and slopes):

```{r, fig.width=8, fig.height=4}
par(mfrow=c(1,2), cex=1.1)
plot(m2, select=3)
plot(m2, select=4)
```


Or alternatively:

```{r, fig.width=4, fig.height=4}
pp <- get_random(m2)

emptyPlot(range(pp[[1]]), range(pp[[2]]), h=0,v=0,
     xlab='Random intercepts', ylab='Random slopes',
     main='Correlation')

text(pp[[1]], pp[[2]], labels=names(pp[[1]]), 
     col='steelblue', xpd=TRUE)
```


### <a id="inspectrandom"></a>inspect_random()

For plotting and extracting the random smooths:

```{r, fig.width=8, fig.height=4}
par(mfrow=c(1,2), cex=1.1)

inspect_random(m1, select=3, main='s(Time, Subject)')

children <- unique(simdat[simdat$Group=="Children", "Subject"])
adults   <- unique(simdat[simdat$Group=="Adults", "Subject"])

inspect_random(m1, select=3, main='Averages', 
      fun=mean, 
      cond=list(Subject=children))
inspect_random(m1, select=3, 
      fun=mean, cond=list(Subject=adults),
      add=TRUE, col='red', lty=5)

# add legend:
legend('bottomleft',
  legend=c('Children', 'Adults'),
  col=c('black', 'red'), lty=c(1,5),
  bty='n')
```

# References

