Recently I am reading a paper where the authors use the GAM to make predictions. In brief, the data looks like following:

y    i    j     x    weekend
5.6    1    1   4.6    Mon.
6.5    1    2   5.6    Mon.
...
4.6    2    1   6.7    Sta.
2.4    2    2   1.2    Sta.
...

where y, x1, x2 are continuous numbers, weekend is the day of the week. In the paper, the authors use the following formula:

$$y_{ij} = \beta_0 + b_{0i} + \beta_1{\rm weekend}_i + f_1(x_{ij}, {\rm weekend}_i) + \varepsilon_{ij}$$

In the formula, $\beta_0$ is the overall mean, $b_{0i}$ is the random intercept, ${\rm weekend}_i$ determines whether it is weekday or weekend. Ans so I transform ${\rm weekend}$ from {Mon., Thu., .., Sun.} into {0, 1}. And $f_1$ is cubic regression function with 17 spline knots, and in fact will generate two smooth functions one for weekday, another for weekend.

I want to use following code:

gam(y~ s(i,bs="re") + weekend + s(x, by=weekend, bs="cr", k=17))

But I'm not sure whether it fits the formula or not. My questions are:

1. gam will automatically generate the mean of the model, so there is no need to specify a $\beta_0$ in the code?
2. Is it right that by using s(i,bs="re"), the gam will calculate different random effect with distribution $N(0, \delta_i)$ for every $i$ specifically?
3. Is it good to transform weekend into 0-1 value? and in the code s(x, by=weekend, bs="cr", k=17), does the by keyword mean that it will generate different smooth functions of x for different weekend value?
4. The last question is that without specifying knots=list(), as in the above code, the default behaviour of the model is to put knot points evenly of the range of value?
• I'm going to be the dissenting voice here: I think this is marginal as regards being OT here. The OP is trying to relate the model specification to their statistical understanding of the problem. This wouldn't be on topic for Stack Overflow and I don't know of a better place for it within the Stack Exchange family. Feb 18, 2015 at 16:46
• That's a reasonable position, @GavinSimpson. I've retracted my close vote. Feb 18, 2015 at 20:53

## Q1

Not really; in R's formula system, the intercept is implied and created when R parses the formula and builds the model matrix. If you want to suppress it, you need to add -1 or + 0 to the formula.

## Q2

No, assuming $i$ is a grouping variable? mgcv::gam() will fit a spline that is equivalent to a random intercept in the variable i, i.e. the intercepts are drawn from a mean zero, Gaussian distribution with a single unknown variance to be estimated from the data.

## Q3

Yes, you will get a spline for the data where weekend == 0 and a different spline for weekend == 1. You don't have to recode this as 0 or 1, just make sure that weekend is a factor variable. it may help for example to have weekend be a factor with levels c("weekday", "weekend"), corresponding to your 0 and 1 respectively as that will help you recall the coding.

## Q4

Yes, boundary knots are placed at the minimum and maximum of the observed data for x and then the remaining knots are spread evenly over the interval of the data. For some spline bases it makes no sense to fiddle with knots, such as the p spline bases (bs = "ps") and some bases don't even use knots, like the thin-plate splines (bs = "tprs") that mgcv::gam() defaults to using.

• does grouping variable mean c(1,2,3...)? sorry, I can't understand answer to Q2. Two ideas come to me. One is that the gam will fit only one variance of random intercept for the whole data. The other is if i is a grouping variable, then gam will fit variances for every different i. If the latter one is wrong, then how to achieve it ?
– 宇宙人
Feb 19, 2015 at 10:25
• A grouping variable might be patient, or lake, or plot. With bs = "re", gam() will fit a standard random intercept (random effect) where the mean of the response within each level of the group (i.e. each patient, or each lake) is assumed to be a mean zero Gaussian random variable with variance to be estimated. In other words, the variability in the mean response between groups is what is then estimated. This is a bog-standard random intercept. It doesn't fit a separate variance for each level of the group (each patient or lake, etc). I didn't think random effects were used for that. Feb 19, 2015 at 17:01
• The gamm() function will fit via the nlme package (in your case with a Gaussian GAM) and you can use the weights argument and a variance function to estimate a difference per group. See ?varClasses after doing library("nlme") Feb 19, 2015 at 17:03