Yes, they are included, but only ever for the observed levels of the random factor. You can turn this using the by variable smooth trick however.
Consider the following example taken from ?gam.models:
dat <- gamSim(1,n=400,scale=2) ## simulate 4 term additive truth
## Now add some random effects to the simulation. Response is
## grouped into one of 20 groups by `fac' and each groups has a
## random effect added....
fac <- as.factor(sample(1:20,400,replace=TRUE))
dat$X <- model.matrix(~fac-1) #$ rendering bug
b <- rnorm(20)*.5
dat$y <- dat$y + dat$X%*%b
rm1 <- gam(y ~ s(fac, bs="re") + s(x0) + s(x1) + s(x2) + s(x3),
data = dat, method = "ML")
Now lets get the additive term contributions from the model and compare them with the full blown model predictions:
p <- predict(rm1, type = "terms")
head(rowSums(p) + attr(p, "constant"))
head(predict(rm1, type = "response"))
which gives
> head(rowSums(p) + attr(p, "constant"))
1 2 3 4 5 6
14.265260 6.433342 2.766193 12.864771 5.296381 7.341790
> head(predict(rm1, type = "response"))
1 2 3 4 5 6
14.265260 6.433342 2.766193 12.864771 5.296381 7.341790
So we are convinced now that the two ways of generating the predicted values are equivalent. now look at p the additive term contributions to the fitted values:
> head(p)
s(fac) s(x0) s(x1) s(x2) s(x3)
1 -0.03786017 -0.1683648 3.868927 2.6485134 0.157054343
2 0.21328630 0.5304765 -1.902366 -0.1972856 -0.007759325
3 -0.36501307 0.1058000 -1.661677 -3.0955348 -0.014372627
4 -0.12519987 0.5474540 2.554656 2.2189534 -0.128083342
5 -0.12519987 -0.3720668 -1.817144 -0.1451364 -0.041061989
6 -0.05481148 0.3490905 -1.216908 0.3411783 0.126251294
The first column is the s(fac) which was a random effect spline in the fitted GAM.
I will add that the gamm() function also in mgcv can give the within-group predictions (fitted values):
m2 <- gamm(y ~ s(x0) + s(x1) + s(x2) + s(x3),
data = dat, method = "ML",
random = list(fac = ~ 1))
head(predict(m2$lme))
> head(predict(m2$lme))
1/1/1/1/14 1/1/1/1/6 1/1/1/1/19 1/1/1/1/9 1/1/1/1/9 1/1/1/1/15
14.265259 6.433345 2.766196 12.864770 5.296383 7.341790
> head(predict(rm1, type = "response"))
1 2 3 4 5 6
14.265260 6.433342 2.766193 12.864771 5.296381 7.341790
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
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.
Best Answer
From version 1.8.8 of mgcv
predict.gamhas gained anexcludeargument which allows for the zeroing out of terms in the model, including random effects, when predicting, without the dummy trick that was suggested previously.Older approach
Simon Wood has used the following simple example to check this is working:
Which works for me. Likewise:
also works.
So I would check the data you are supplying in
newdatais what you think it is as the problem may not be withVesselID— the error is coming from the function that would have been called by thepredict()calls in the examples above, andRailis a factor in the first example.