1: include `rm.ranef=TRUE` to zero all random effects.
2: include `all.terms=TRUE` to visualize parametric terms.
# 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`.
### 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))`.
### 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
### 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')
```
### 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
### 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)`.
### 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
### 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)
```
### 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