11
$\begingroup$

I have a binomial GAM with a smooth-by-random-factor interaction. From this I am able to predict and visualize the smooth term for any level of my random effects:

#Simulate data
set.seed(0)
means = rnorm(5, mean=0, sd=2)
df = data_frame(group = as.factor(rep(1:5, each=100)),
                x = rep(seq(-3,3, length.out =100), 5),
                y=as.numeric(dnorm(x, mean=means[group]) > 0.4*runif(10)))

#Fit model
library(mgcv)
gam_model = gam(y ~ te(x, group, bs=c("ts", "re")), data=df, family = binomial)

#Visualize
df2 = predict(gam_model, type="response", se.fit=TRUE)
df2 = cbind(df, response = df2$fit, lwr = df2$fit-2*df2$se.fit, upr = df2$fit+2*df2$se.fit)

library(ggplot2)
ggplot() +
  geom_ribbon(data = df2, mapping=aes(x=x, ymin=lwr, ymax=upr, fill=group), alpha=0.25) +
  geom_line(data = df2, mapping=aes(x=x, y=response, col=group)) +
  geom_point(data = df, mapping=aes(x=x, y=y, col=group)) +
  facet_wrap(~group)

enter image description here

How do I predict the mean smooth and confidence intervals around it?

$\endgroup$
1
  • 6
    $\begingroup$ Update: The top answerer and the questioner (with friends) now have an entire paper on this topic: peerj.com/articles/6876 $\endgroup$
    – Noam Ross
    Jun 17, 2020 at 13:10

2 Answers 2

11
$\begingroup$

The solution suggested by Simon Wood to the simpler problem of predicting the population level effect from a model with random intercepts represented as a smooth is to use a by variable in the random effect smooth. See this Answer for some detail.

You can't do this dummy trick directly with your model as you have the smooth and random effects all bound up in the 2d spline term. As I understand it, you should be able to decompose your tensor product spline into "main effects" and the "spline interaction". I quote these as the decomposition will be to split out the fixed effects and random effects parts of the model.

Nb: I think I have this right but it would be helpful to have people knowledgeable with mgcv give this a once over.

## load packages
library("mgcv")
library("ggplot2")
set.seed(0)
means <- rnorm(5, mean=0, sd=2)
group <- as.factor(rep(1:5, each=100))

## generate data
df <- data.frame(group = group,
                 x = rep(seq(-3,3, length.out =100), 5),
                 y = as.numeric(dnorm(x, mean=means[group]) > 
                       0.4*runif(10)),
                 dummy = 1) # dummy variable trick

This is what I came up with:

gam_model3 <- gam(y ~ s(x, bs = "ts") + s(group, bs = "re", by = dummy) + 
                  ti(x, group, bs = c("ts","re"), by = dummy),
                  data = df, family = binomial, method = "REML")

Here I've broken out the fixed effects smooth of x, the random intercepts and the random - smooth interaction. Each of the random effect terms includes by = dummy. This allows us to zero out these terms by switching dummy to be a vector of 0s. This works because by terms here multiply the smooth by a numeric value; where dummy == 1 we get the effect of the random effect smooth but when dummy == 0 we are multiplying the effect of each random effect smoother by 0.

To get the population level we need just the effect of s(x, bs = "ts") and zero out the other terms.

newdf <- data.frame(group = as.factor(rep(1, 100)), 
                    x = seq(-3, 3, length = 100),
                    dummy = rep(0, 100)) # zero out ranef terms
ilink <- family(gam_model3)$linkinv      # inverse link function
df2 <- predict(gam_model3, newdf, se.fit = TRUE)
ilink <- family(gam_model3)$linkinv
df2 <- with(df2, data.frame(newdf,
                            response = ilink(fit),
                            lwr = ilink(fit - 2*se.fit),
                            upr = ilink(fit + 2*se.fit)))

(Note that all this was done on the scale of the linear predictor and only backtransformed at the end using ilink())

Here's what the population-level effect looks like

theme_set(theme_bw())
p <- ggplot(df2, aes(x = x, y = response)) +
geom_point(data = df, aes(x = x, y = y, colour = group)) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1) +
geom_line()
p

enter image description here

And here are the group level smooths with the population level one superimposed

df3 <- predict(gam_model3, se.fit = TRUE)
df3 <- with(df3, data.frame(df,
                            response = ilink(fit),
                            lwr = ilink(fit - 2*se.fit),
                            upr = ilink(fit + 2*se.fit)))

and a plot

p2 <- ggplot(df3, aes(x = x, y = response)) +
geom_point(data = df, aes(x = x, y = y, colour = group)) +
geom_ribbon(aes(ymin = lwr, ymax = upr, fill = group), alpha = 0.1) +
geom_line(aes(colour = group)) +
geom_ribbon(data = df2, aes(ymin = lwr, ymax = upr), alpha = 0.1) +
geom_line(data = df2, aes(y = response))
p2

From a cursory inspection this looks qualitatively similar to the result from Ben's answer but it is smoother; you don't get the blips where the next group's data is not all zero.

enter image description here

$\endgroup$
8
  • $\begingroup$ Love all of these answers, the question, and the paper! Can we use an approach like this to predict response from a known x and an unknown group? It seems like zeroing out the other smooths in the model (leaving just s(x, bs = "ts") would get me the expectation of the response for an unknown group, but what are those confidence bands representing? $\endgroup$
    – mikoontz
    Jan 31, 2023 at 7:09
  • 2
    $\begingroup$ There is a better way to do this above now (I should write a new answer) as there are terms and exclude arguments to predict.gam() which allow you to specify either the terms to include in or the terms to exclude from the predicted values. I'm pretty sure the bands are only including the uncertainty from the "population" smooth, and not anything that is due to the random effect distribution. I wonder if fitting with gamm4() might allow what you want if glmer can do it for GLMMs? $\endgroup$ Jan 31, 2023 at 12:17
  • 1
    $\begingroup$ Great, thanks! I’ll explore the ‘terms’ and ‘exclude’ arguments to get the expectation for an unknown group based on the population smooth, and will check out ‘gamm4()’ to see if I can’t also incorporate the uncertainty from the random effects distribution. I’ll pose a new question if I get stuck, and report back here if I think my attempt with ‘terms’ and ‘exclude’ will be useful to others. Thanks so much! $\endgroup$
    – mikoontz
    Feb 1, 2023 at 14:23
  • $\begingroup$ I'd be interested in doing this with my model that has two continuous variables and one categorical (separate smooths and 3-way interaction, both implemented with the "by" argument), but I currently can't imagine how I might slip an additional dummy variable in there when the "by" argument is already used - does this call for a new question or can you elaborate here? $\endgroup$
    – mluerig
    May 11, 2023 at 10:42
  • 1
    $\begingroup$ @mluerig I should add another answer here, or edit this one, as things have changed in regards to {mgcv} now such that you don't need this dummy trick. See ?predict.gam and note the terms and exclude arguments; you use one or the other to include or exclude terms from the linear predictor when predicting. In this case, it is easiest to use exclude to ignore the random effect term when predicting, so exclude = c("s(group)", "ti(x,group)") and note that you have to get these exactly right - look at summary() to see how the smooths are labelled. $\endgroup$ May 11, 2023 at 10:49
4
$\begingroup$

It depends. There are a bunch of ways to define the "average" response. The answer here is based on the unweighted average across groups; with this simple, artificial example it doesn't make any difference, but in other cases you might want to take a population-weighted average.

n.b. there are several reasons the following is not quite right, although it's not an unreasonable start

  • for consistency, we ought to be taking the mean / combining the standard errors on the linear predictor (link) scale, not the response scale
  • the answer below essentially treats the groups as fixed effects. We know more about the distribution (at least the assumed distribution) of the conditional modes ... but it means there are a lot of possible definitions

I will update when I get a chance, but it's still a slightly useful answer

Repeating data generation for my convenience ...

library("dplyr") ## for data_frame
set.seed(0)
means = rnorm(5, mean=0, sd=2)
df = data_frame(group = as.factor(rep(1:5, each=100)),
            x = rep(seq(-3,3, length.out =100), 5),
            y=as.numeric(dnorm(x, mean=means[group]) > 0.4*runif(10)))

#Fit model
library(mgcv)
gam_model = gam(y ~ te(x, group, bs=c("ts", "re")),
                 data=df, family = binomial)
gam_avg = gam(y ~ s(x), data=df, family = binomial)

Tweak prediction step a tiny bit to retain se in the results (it would be nice to write a broom::augment() method for this case ...)

#Predict
pfun <- function(x, type="response") {
      pp <- predict(x, type="response", se.fit=TRUE)

df2 = predict(gam_model, type="response", se.fit=TRUE)
df2 = with(df2,data.frame(df,
                      response = fit,
                      se= se.fit, lwr = fit-2*se.fit,
                      upr = fit+2*se.fit))

Generate mean predictions by averaging at each value of x; construct confidence intervals by adding "in quadrature" (i.e. sqrt(sum(x^2))) (I don't know why c() is necessary, but it seems to be).

sumquad <- function(x) { sqrt(sum(c(x)^2)) }
dfsum <- df2 %>% group_by(x) %>%
   summarise(response=mean(c(response)),
             se=sumquad(se)) %>%
   mutate(lwr=response-2*se,upr=response+2*se)

Now visualize:

library("ggplot2"); theme_set(theme_bw())
gg1 <- ggplot(mapping=aes(x)) +
        geom_ribbon(data = df2,
          mapping=aes(ymin=lwr, ymax=upr, fill=group),
              alpha=0.25) +
       geom_line(data = df2, mapping=aes(y=response, col=group)) +
       geom_point(data = df, mapping=aes(y=y, col=group))

## add mean response + ribbon                                              
gg1 + geom_line(data=dfsum,aes(y=response))+
    geom_ribbon(data=dfsum,aes(ymin=lwr,ymax=upr),alpha=0.2)
$\endgroup$
2
  • $\begingroup$ To switch to the link scale for this answer, would it require more than the following? 1) changing to type="link" and 2) adding ` mutate(response=plogis(link), re_upr=plogis(upr), re_lwr=plogis(lwr))` for both df2 and dfsum $\endgroup$
    – jaimedash
    Jan 5, 2016 at 19:08
  • 2
    $\begingroup$ @jaimedash no (but plogis(response)), I don't think so (i.e. yes I think this is all you need). You're supposed to do all these things on the scale of the linear predictor and then back transform. More generally you can always grab the inverse of the link function via ilink <- family(mod)$linkinv and use ilink(foo). $\endgroup$ Jan 5, 2016 at 19:32

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.