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 &lt;- 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 &lt;- as.factor(sample(1:20,400,replace=TRUE))
dat$X &lt;- model.matrix(~fac-1)                                  #$ rendering bug
b &lt;- rnorm(20)*.5
dat$y &lt;- dat$y + dat$X%*%b

rm1 &lt;- 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 &lt;- predict(rm1, type = "terms")
head(rowSums(p) + attr(p, "constant"))
head(predict(rm1, type = "response"))


which gives

&gt; head(rowSums(p) + attr(p, "constant"))
        1         2         3         4         5         6 
14.265260  6.433342  2.766193 12.864771  5.296381  7.341790 
&gt; 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:

&gt; 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 &lt;- gamm(y ~ s(x0) + s(x1) + s(x2) + s(x3),
           data = dat, method = "ML",
           random = list(fac = ~ 1))

head(predict(m2$lme))

&gt; 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 
&gt; head(predict(rm1, type = "response"))
        1         2         3         4         5         6 
14.265260  6.433342  2.766193 12.864771  5.296381  7.341790


gam models with random effect R
I have rocks orbiting moons, moons orbiting planets, planets orbiting suns, and suns orbiting black holes, and the current system could have many many layers of orbitage.

the position of any object is a function of time and relative to the object it orbits.
(so far so good).

now I want to know for a given 2 objects(A,B), a start time and a speed, how can I work out the when and where to go. I can work out where A and B is given a time.
so i just need.

1: direction to travel in from A to B(remember B is moving(not in a straight line))
2: Time to get to b in a straight line. 

travel must be in a straight line with the shortest possible distance.

as an extension to this question, how will i know if its better to wait, EG is it faster to stay on object A and wait for a hour when the objects may be closer, than to set off from A to B at the start.

Cheers, it hurt my brain.
