# In R package mgcv, is it valid to have a random effect smooth on two continuous variables?

In my model, I have two variables, distance and time, that influence the performance of different subjects (success or failure). The relationship between distance and time is nonlinear, not separable into single additive terms, and (probably) vary by subject. There are other variables involved in other terms in the model, but my focus at the moment is these two variables.

I don't have a parametric form that I can expect, but I do have thousands of observations across hundreds of subjects, and so can estimate the overall interaction between distance and time nonparametrically. I have been using mgcv::bam (the large-data version of mgcv::gam) to estimate a (simplified) model along the lines of:

mgcv::bam(success ~ s(distance, time), data=mydata)


However, what I'd like to explore is how this distance + time relationship varies across subjects, using a random effects framework. Something like:

mgcv::bam(success ~ s(distance, time, subject, bs="re"), data=mydata)


Running R code this way doesn't produce an error, but I am not sure it's doing what I think it is. I have not found any references to using multiple numeric variables in a mgcv random effect smooth this way, and wonder if there's a reason for that.

As a side note, one approach might be to combine distance and time into a single variable like "speed" (distance / time), and then use:

mgcv::bam(success ~ s(speed, subject, bs="re"), data=mydata)


... but this potentially loses some of the details in the complex way distance and time interact, leading to my question.

Consulting the documentation for the mgcv::gam function in R, you can specify a smooth on a factor variable to be a random slope of another numeric variable.

If g is a factor and x is numeric, then s(x,g,bs="re") produces an i.i.d. normal random slope relating the response to x for each level of g.

It's unclear to me whether you can specify more than one numeric variable in conjunction with a random effect smooth. i.e. s(x, y, g, bs="re"). Or even if this is a sensible think to ask.

Simplified example:

## adapted from example in mgcv docs
dat <- gamSim(1,n=400,scale=2) ## simulate 4 term additive truth

fac <- sample(1:20,400,replace=TRUE)
b <- rnorm(20)*.5
dat$y <- dat$y + b[fac]
dat$fac <- as.factor(fac) rm1 <- gam(y ~ s(x0, fac,bs="re") + s(x1, x2, fac, bs="re") + s(x3),data=dat,method="ML")  The smooth on x0 & fac is ok per the documentation, but I don't know if the smooth on x1, x2, & fac is valid and interpretable. • The other place to ask might be the R-SIG-mixed-models mailing list. – DWin May 2 '18 at 19:53 ## 2 Answers Huh... I made my post as a guest on SO because I am still on a suspension, but then the question got migrated here! So, if I understand you correctly, there's not really any similarity between the smooth s(x1, x2) and the random effect s(x1, x2, fac, bs = "re"), correct? Correct. The function name "s" does not mean "smooth function" when s() is used to construct a random effect. Broadly speaking, s() is just a model term constructor routine that constructs a design matrix and a penalty matrix. What I was envisioning was something smoothing in 2 dimensions like the former, but with some deviations from the average by factor level. You can get separate smooths per factor level using s(x1, x2, by=fac), but that completely separates the data for each factor level, rather than doing some partial pooling. s(x1, x2, by = fac) gives you something pretty close to what you want, except that as you said, data from different factor levels are treated independently. Technically, "close" means that s(x1, x2, by = fac) gives you the correct design matrix but not the correct penalty matrix. In this regard, you are probably aiming at te(x1, x2, fac, d = c(2, 1), bs = c("tp", "re")). I have never seen such model term before, but its construction is definitely possible in mgcv: library(mgcv) x1 <- runif(1000) x2 <- runif(1000) f <- gl(5, 200) ## "smooth.spec" object smooth_spec <- te(x1, x2, f, d = c(2, 1), bs = c("tp", "re")) ## "smooth" object sm <- smooth.construct(smooth_spec, data = list(x1 = x1, x2 = x2, f = f), knots = NULL)  You can check that this smooth term has 2 smoothing parameters as expected, one for the s(x1, x2, bs = 'tp') margin, the other for the s(f, bs = 're') margin. Specification of k turns out subtle. You need to explicitly pass nlevels(f) to the random effect margin. For example, if you want a rank-10 thin-plate regression spline, ## my example factor f has 5 levels smooth_spec <- te(x1, x2, f, d = c(2, 1), bs = c("tp", "re"), k = c(10, 5)) sapply(smooth_spec$margin, "[[", "bs.dim")
# [1] 10  5


At first I was thinking that perhaps we can simply pass NA to the random effect margin, but it turns out not!

smooth_spec <- te(x1, x2, f, d = c(2, 1), bs = c("tp", "re"), k = c(10, NA))
sapply(smooth_spec$margin, "[[", "bs.dim") # [1] 25 5 ## ?? why is it 25? something has gone wrong!  This might implies that there is a tiny bug... will have a check when available. • Thanks for the detailed explanation, and the suggestion of using te() as outlined. I've accepted your answer and will try that out! May 8 '18 at 14:05 • +1 I'll quibble a bit about "The function name "s" does not mean "smooth function"...". It really depends what you mean by "smooth function". Here one really means a spline with an associated wiggliness penalty in quadratic form. In that sense, the re spline is really a spline, just with certain basis functions and penalty that allows it to implement a random effect. It really is a spline, just not one you might easily visualise as a wiggly line. I have a couple of comments... May 9 '18 at 15:38 • (just one it turns out) ...I would not use a 2-d TPRS spline for distance and time in the OP's examples, unless you really want the same wiggliness in both distance and time, which seems unlikely. I think the only way to get separate smoothness parameters for distance and time would be to use a tensor product of three univariate marginals. I.e. te(distance, time, subject, bs = c("tp", "tp", "re")). May 9 '18 at 15:49 • @GavinSimpson Yes, I was aware of that, but did not comment it to OP. Also, I was wondering that he should not use Gaussian family if he works with binary data. I think you could produce another answer on such practical aspect? It might be better to use cubic regression spline for 1D spline margin. May 9 '18 at 15:50 • @李哲源 Who says a spline has to be be smooth in an english/language sense of the word? (I suppose you mean in the sense of second derivative?) You can construct splines that use different definitions of smoothness and basis functions. And I did say "quibble a bit" (emphasis on the "bit", i.e. only a little). May 9 '18 at 16:02 Consider the random effect part of your example toy <- gam(y ~ s(x0, fac, bs = "re") + s(x1, x2, fac, bs="re"), data = dat, method = "REML")  This is just the penalized version of the following linear regression model: toy.lm <- lm(y ~ x0:fac + x1:x2:fac, data = dat)  where a ridge penalty is applied to x0:fac and x1:x2:fac. The construction of simple random effect in mgcv::gam or mgcv::bam is fairly routine: • generate design matrix from X1 <- model.matrix(~x0:fac - 1, data = dat) X2 <- model.matrix(~x1:x2:fac - 1, data = dat)  • generate ridge penalty matrix S2 <- S1 <- diag(nlevels(dat$fac))

• Thanks for that explanation. So, if I understand you correctly, there's not really any similarity between the smooth s(x1, x2) and the random effect s(x1, x2, fac, bs ="re"), correct? What I was envisioning was something smoothing in 2 dimensions like the former, but with some deviations from the average by factor level. You can get separate smooths per factor level using s(x1, x2, by=fac), but that completely separates the data for each factor level, rather than doing some partial pooling. Sounds like that might be a little too ambitious a design. Is that a fair assessment? May 7 '18 at 21:37