1

I am modeling fishery CPUE as a function of a number of a number of covariates using a GAM approach that includes fixed and random effects.

I understand that there are limitations with regards to predicting random effects (predict function only addresses fixed effects) with gamm4. How does the predict function in the basic gam (mgcv, using bs="re") deal with the random effects? Are they included in predictions? Any thoughts would be much appreciated...

Peterson
  • 11
  • 1
  • 3
  • I should add that, unless there is wording to the contrary somewhere, `gamm4` returns the full `lmer` fit and according to the documentation for the **lme4** package, predictions can be at the population level only if requested. `predict` on the `$mer` component of the returned object should return predictions conditional upon the random effects (i.e. including them). – Gavin Simpson Jun 30 '15 at 22:28
  • Thanks Gavin for all your input. Very helpful. I have tried to predict on the $mer component using gamm4; however, I have been unable to get around the following error message. "Error: couldn't evaluate grouping factor Xr within model frame: try adding grouping factor to data frame explicitly if possible". My random effects are factors, so I am unsure where my problem is. It appears that the gam option may be a workaround. – Peterson Jul 01 '15 at 17:29

1 Answers1

3

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
Gavin Simpson
  • 37,567
  • 5
  • 110
  • 153