I am fitting GAM models to check whether forest treatment (2 types of logging regimes) influence bird abundance across years. Abundance was counted on constant plots. Each plot is located in a constant treatment area. Bird abundance was surveyed in each plot 6 times between 2005-2020.

I created a simplified, reproducible example, which mirrors my dataset and fitted a GAM model:

plot <- rep(sprintf("p%s",seq(1:18)), each=6)
treatment <- rep(c("control", "treatment1", "treatment2"),each=36)
year <- rep(c(2005,2007,2008,2010,2012,2020), 18)
abundance <- c(sort(runif(36, min = 1, max = 40), decreasing = TRUE), sort(runif(36, min = 1, max = 35), decreasing = TRUE), sort(runif(36, min = 16, max = 32)))
piska_df <- as.data.frame(cbind(plot,treatment,year,abundance))
piska_df$plot <- as.factor(piska_df$plot)
piska_df$treatment <- as.factor(piska_df$treatment)
piska_df$abundance <- as.integer(piska_df$abundance)
piska_df$year <- as.integer(piska_df$year)

g1<-gam(abundance ~  treatment*year + s(plot,bs="re"), data=piska_df, family=poisson, method="REML")

The model works great! But i am now stuck on visualizing main results in a clear manner.

I decided to calculate predicted values for each treatment separately across years, while keeping random factor (“plot”) constraint. Afterwards I transformed my data using inv.logit to get true abundance values for birds. I calculated 95% CI based on SE. This is the code that I used:

year.pr <-seq(min(piska_df$year),max(piska_df$year), length.out = 100)
new_data_t2 <- as.data.frame(new_data_t2)
new_data_t1 <- as.data.frame(new_data_t1)
new_data_ctrl <- as.data.frame(new_data_ctrl)

ilink <- family(g1)$linkinv
g.pred.ctrl <- predict(g1,newdata=new_data_ctrl,
                           type="link",se.fit = TRUE)
g.pred.t1 <-predict(g1,newdata=new_data_t1,
                          type="link",se.fit = TRUE)
g.pred.t2 <-predict(g1,newdata=new_data_t2,
                          type="link",se.fit = TRUE)

g.pred.ctrl <- cbind(g.pred.ctrl, new_data_ctrl)
g.pred.ctrl <- transform(g.pred.ctrl, lwr_ci = ilink(fit - (2 * se.fit)),
                              upr_ci = ilink(fit + (2 * se.fit)),
                              fitted = ilink(fit))

g.pred.t1 <- cbind(g.pred.t1, new_data_t1)
g.pred.t1 <- transform(g.pred.t1, lwr_ci = ilink(fit - (2 * se.fit)),
                             upr_ci = ilink(fit + (2 * se.fit)),
                             fitted = ilink(fit))

g.pred.t2 <- cbind(g.pred.t2, new_data_t2)
g.pred.t2 <- transform(g.pred.t2, lwr_ci = ilink(fit - (2 * se.fit)),
                             upr_ci = ilink(fit + (2 * se.fit)),
                             fitted = ilink(fit))

g.pred.all <- rbind(g.pred.t2,g.pred.t1,g.pred.ctrl)

Then I plotted a ggplot graph, using the predicted values:

ggplot(g.pred.all, aes(x = year, y = fitted, colour = factor(treatment))) +
  theme_classic() +
  geom_ribbon(aes(ymin = lwr_ci, ymax = upr_ci, fill = factor(treatment)), alpha = 0.1) +
  geom_line(linewidth=1.5) +
  ggtitle("abundance") + xlab("year")

This is the graph:

enter image description here

And here comes my problem. I am interested in how treatment differs from control – I want it to be main focus of those graphs. I am not interested in general decrease/increase, I am interested in decrease/increase in relation to control.

Therefore I thought it would be a nice idea if I had control as a horizontal 0 line (with respective confidence intervals). Then my Y axis would become “abundance difference from control” instead of “abundance”.

My question would be: how to transform my predictions so that I can show control as straight line going through 0 while maintaining “true” mathematical relations between points and confidence intervals? Can I just calculate difference between all other values & mean control and plot this on the graph? Does it make sense mathematically speaking? Should the CI values be somehow recalculated?

All help would be very valuable. I am also open to any other simple and convenient ways to visualize those results (simple visualization of 3 treatments and their trends over years).

Thank you very much in advance.

  • 1
    $\begingroup$ You could use emmeans to calculate the difference between control and the treatment groups at some values of year. If you want the differences instead of ratios, it's a bit tricky. But this here works: contrast(emmeans(ref_grid(g1, at = list(year = seq(2005, 2020, length = 20)), regrid = "response"), "treatment", by = "year"), "trt.vs.ctrl", infer = c(TRUE, FALSE)). Convert this to a data frame to plot it using ggplot, for example. Caution: emmeans will average over the random effects, unlike predict. $\endgroup$ Feb 16 at 21:17
  • $\begingroup$ Thanks a lot for this awesome suggestion! I am okay with the average over the random effects. However, i run the analysis & made a ggplot - the CI of contrasts overlap 0 in both cases. The model, however, shows statistical significance of control-treatment2 difference at p<0.05 . Doesn't it mean that the CIs of contrasts should not overlap 0? Should i not worry about that? $\endgroup$ Feb 17 at 8:54
  • 1
    $\begingroup$ Because you included an interaction between treatment and year, the main effects cannot be intepreted as usual. The main effect of treatment compares the treatments to the control for a year of $0$, which makes little sense. Further, the contrasts at each year are adjusted according to Dunnett's test, which will enlarge the confidence intervals. In summary, the plot and CI are consistent with the model output. $\endgroup$ Feb 17 at 10:21
  • $\begingroup$ Thank you. One more small question: You say that "the main effects cannot be intepreted as usual". Does it mean that i should not interpret GAM model the way i intend to? Meaning: if there is statistical significance of control-treatment2 difference at p<0.05, can i say based on my GAM model that treatment2 differed from control? $\endgroup$ Feb 17 at 10:40
  • 1
    $\begingroup$ You could use anova(g1) to check whether the interaction is significant. If it is, there is no point in interpreting the significance of the main effects because you already have evidence now that the abundance is associated with treatment and year. In other words: There is evidence that the differences between treatments is not constant over the years (lines are not parallel on the scale of the linear predictor). It's possible that there is little difference in year 2005 but a large difference in 2020. The graph will show that. $\endgroup$ Feb 17 at 13:13

1 Answer 1


Comparing multiple treatments with a control can be done using Dunnett's test. The emmeans package makes this convenient. Some minor tweaks are required because by default, emmeans will compare the groups using ratios instead of differences (i.e. differences on the log-scale). Furthermore, emmeans will adjust $p$-values and confidence intervals using Dunnett's test within each year. Note: You currently estimate the model using restricted maximum likelihood (REML). This is fine if you want unbiased tests for the random effects but suboptimal for comparison of fixed effects. I suggest refitting the model using maximum likelihood (ML) before doing these comparisons. The figure below is for your model fitted with maximum likelihood.

Here is the code (I assume your code has been run before):


# Setting up the reference grid using 20 values for year
refgrid <- ref_grid(g1, at = list(year = seq(2005, 2020, length = 20)), regrid = "response")

# Marginal means
em <- emmeans(refgrid, "treatment", by = "year")

# Dunnett's test for each year ("trt.vs.ctrl")
contrs <- contrast(em, "trt.vs.ctrl", infer = c(TRUE, FALSE))

# Convert to data frame for plotting
plot_dat <- as.data.frame(contrs)


  • 2
    $\begingroup$ Is it okay to fit the model with REML for these comparisons? And for the ANOVA, I think one needs to use ML instead? $\endgroup$
    – dipetkov
    Feb 17 at 18:29
  • 1
    $\begingroup$ Thank you for the advice but I already know not to compare models with different fixed effects fitted with REML. I thought the OP might want to know too. See here: REML or ML to compare two mixed effects models with differing fixed effects, but with the same random effect? $\endgroup$
    – dipetkov
    Feb 17 at 19:05
  • 2
    $\begingroup$ @dipetkov Good point, I missed this in your first comment (I thought you are OP). I updated the answer and the graph accordingly. $\endgroup$ Feb 17 at 21:38
  • 1
    $\begingroup$ Thank you both, i am going to rerun my models using ML then. $\endgroup$ Feb 20 at 8:42

Your Answer

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

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