I am new in GAMM and I am analyzing a data to find the effect of year, tree species and their interaction on richness of insects on the trees. I expect that year will have an overall nonlinear effect. Furthermore, the effect of year on richness will have different shapes for each tree species. I limit k as I do not need too many knots and otherwise the model does not work. My model is as fallows: m1 <- gam(richness ~ s(year_begin,k=6) + Tree + s(year_begin,Tree, by = Tree, k=6),random=~(1|Region/PlotID/TreeID) + (1|observation), family= quasipoisson, data)

The model works fine and results are OK. Then I try to plot the summed effects for each tree species (levels of the factor Tree) e.g. AH -maple:vis.gam(m1, view=c("year_begin", "Tree"), cond=list(Tree="AH"), plot.type='persp') vis.gam fails to subset as I used Tree factor coerced in the model (Tree) for view argument. At the end I got all 13 levels stacked in the plotsummed effect plot for all factors No labeling or transparency is possible to see which smooth is which tree. I also tried a model with only the interaction term m2 <- gam(richness ~ s(year_begin, Tree, bs="fs",k=6), random... and plotted it: plot(m2, select=1, rug=FALSE) But again there are no labels on the smooths so I cannot understand which one is which tree. enter image description here Id appreciate if anyone can tell me how to plot the summed effects of tree species separately or together but with labels. Thanks in advance

  • $\begingroup$ Once the model is fit, you can define new "predicting" data frames, in which you give values to the categorical variable Tree and then use the function "predict". For instance, make a loop over all levels of the variable Tree and make a prediction against all years, like this you get 13 plots, one for each factor level. $\endgroup$
    – nukimov
    Jun 18, 2019 at 12:53

1 Answer 1


Your model didn't work fine, not if that is the model you ran anyway, as the gam() function doesn't have an argument random. This looks like lme4 notation so did you use gamm4() instead of gam()? If you used gam(), your random effect was silently eaten by the ... argument and promptly ignored.

Using plot() on the model (the $gam component of the model in the case of gamm4() fits) will show the estimated smooths for factor by models. As you've seen, this doesn't work so well for the fs basis because it might make more sense with this formulation to not worry about the specific identify (these are random splines) but focus on the variability.

That said, the solution here is to predict from the model for each species.

pdat <- with(data,
             expand.grid(year_begin = seq(min(year_begin), max(year_begin), length = 100),
                         Tree = levels(Tree)))


pdat <- cbind(pdat, fitted = predict(model, newdata = pdat, type = "response"))

will get you predictions for 100 equally-spaced, ordered values over the range of year_begin for each species, which you can plot, and you'll have the indicator of which species each curve belongs in the pdat variable.

Or you can use functions from my gratia package to get a tidy data frame equivalent of what plot.gam() using evaluate_smooth(model, "smooth").

  • $\begingroup$ Thanks Gavin Simpson. Yes, I use gam4. It is OK with random design but not handy in plotting for my purposes. Your solution works fine. $\endgroup$
    – Centaur
    Jul 3, 2019 at 22:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.