I am fitting a 2D tensor product smooth to a data set that includes 3 years of sampling. I want to feed year as a factor to the "by" argument so that an independent 2D smooth will be created for each year. My model takes the form:

y ~ year + te(x, z, by = year, bs = c("tp", "tp"), k = c(10, 10), m = 1)

This results in this rather fetching plot (using mgcViz).

enter image description here

The problem is that my sampling in xz space was not identical across the three years, so the plots should have different null space arising from the "too.far" argument. Instead, they all have the same null space, presumably based on the sampling distribution pooled across the three years.

A deeper problem is that I don't quite understand how the "by" argument deals with uneven sampling across factor levels. Does it attempt to model, for each factor level, the xz space of the pooled factor levels? This would seem problematic when the factor levels differ in xz space. Yet this seems to be just what it is doing with my data, since the "by year" smooths I get extend beyond the xz space actually sampled in a given year.

In summary, is there a way to handle uneven sampling across factor levels using the "by" argument, or should I settle for fitting a separate model for each year of my data?

  • 1
    $\begingroup$ Thanks, @Gavin Simpson. Is there any difference, though, between (1) setting up a single model with separate smooths for three dummy variables and (2) setting up three separate models based on annual subsets of my data? $\endgroup$
    – dbspon
    Commented Jun 17, 2020 at 4:57
  • $\begingroup$ The flippant answer is yes; you have 1 model with (1) or 3 models with (2). While I said you lose easy forms of comparison or inference, there is always something you can do if you have 1 model, but if you choose 3 separate models you have no easy ways of comparison at all. Personally I thin (1) is always better than (3)... $\endgroup$ Commented Jun 17, 2020 at 5:09

1 Answer 1


The basis functions for each of the smooths are the same, and are set up using the entire range of x and (I think you meant) z in the data.

The coefficients for those basis functions, which is what gives you the estimated surfaces, are estimated using only the x and z values for the individual levels. If you look at the model matrix for one of these models, the columns for the basis functions not associated with the $i$th observation are all set to zero; this part of the model matrix has a block diagonal appearance (assuming the observations/rows are in factor level order).

The plots (in mgcv and gratia at least, I'm not as familiar with mgcViz, but I suspect it is the same) are produced by first predicting for a grid over the full range of the data, i.e.

expand.grid(x = seq(min(x), max(x), length = 100),
            z = seq(min(z), max(z), length = 100))

and then any covariate combination that is further than some distance (Euclidean distance greater than 0.1 when data are cast into the interval 0,1, is the default) from an observation in the entire data set are set to NA.

This implies that in your case there is some extrapolation beyond the support of the data for a particular level, but the penalty and basis functions controlling the individual fits are in effect out beyond the data for that level to stop the smooth going too wiggly, but the model has to do something for those basis functions and as such it is estimating values for them conditional upon the other covariates in the model.

This isn't quite so bad as the basis functions aren't entirely local so the values that get estimated in some sense contribute to the fit of the data points.

If you want to have smooths that are entirely set up using the data for individual levels of a factor, you'll need to create by hand a dummy for each level, with a 1 in the $i$th row if the $i$th observation belongs to that level or a 0 if not. Then you would explicitly write three by variable smooths into the formula. For example with pseudo-ish code as there no reproducible example here:

data <- transform(data,
                  year2010 = ifelse(year == 2010, 1, 0),
                  year2011 = ifelse(year == 2011, 1, 0),
                  year2012 = ifelse(year == 2012, 1, 0))

y ~ year + te(x, z, by = year2010, bs = c("tp", "tp"), k = c(10, 10)) + 
           te(x, z, by = year2011, bs = c("tp", "tp"), k = c(10, 10)) + 
           te(x, z, by = year2012, bs = c("tp", "tp"), k = c(10, 10))

It doesn't matter that we use year to account for the year-specific means but use different by variables in the smooths as the by variables are just used as indicators to identify which rows need to be set to zero as the don't contribute to particular years/levels.

Doing this means you don't have a formal statistical way to compare the surfaces etc, but you could use posterior simulation to do that it you needed to do it. You also don't have a simple way to compare a model with separate te() smooths by years and one with a single te() for all years with the alternative I mention because the coefficient represent different basis functions and thus the model I showed is not nested in the simpler model:

y ~ year + te(x, z, bs = c("tp", "tp"), k = c(10, 10))

If you don't want to go to all that trouble but just don't like that the plots don't have different areas set to be NA, you could just evaluate the estimated smooths yourself at a grid of points, and then apply the exclude.too.far() three times, once per year, to only those rows for the particular year.

Another solution would be to treat the surfaces as random, and fit an "global" smooth for all years and then smooth deviations from that overall surface. You can do that using:

y ~ year + te(x, z, bs = c("tp", "tp"), k = c(10, 10)) +
           te(x, z, by = year, bs = c("tp", "tp"), k = c(10, 10), m = 1)

if you expect each year to have a different wiggliness (possibly with select = TRUE), or if you want to have the same wiggliness you could use

y ~ year + te(x, z, bs = c("tp", "tp"), k = c(10, 10)) +
           te(x, z, year, bs = c("tp", "tp", "re"), k = c(10, 10))


y ~ year + te(x, z, bs = c("tp", "tp"), k = c(10, 10, NA)) +
           t2(x, z, year, bs = c("tp", "tp", "re"), 
              k = c(10, 10, NA), full = TRUE)

The general idea here is that we can learn something about the smooth relationships between the response y and the covariates x and z from all years, rather than throwing up our hands and saying we know nothing — which is what setting something to NA amounts to. These kind of model are very much in the spirit of mixed effects, or hierarchical, models.

For more on these latter types of models, colleagues I wrote a (open access) paper on them last year, with lots of R examples for mgcv: Pedersen et al (2019).


Pedersen, E.J., Miller, D.L., Simpson, G.L., Ross, N., 2019. Hierarchical generalized additive models in ecology: an introduction with mgcv. PeerJ 7, e6876. https://doi.org/10.7717/peerj.6876


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.