8
$\begingroup$

Recently I'm trying to model daily counts using negative binomial distribution within mgcv's GAM model in order to obtain RR's. The data comes from 20 years of daily time series with weekly and monthly seasonality. Given prior series decomposition I assumed that seasonality is evolving with each year. Hence I used tensor products to model it:

library(mgcv)

gam_model <- gam(counts ~ te(year, day_of_week, k=c(20,7), bs=c("tp","cc") +
                          te(year, day_of_year, k=c(20,12), bs=c("tp", "cc"),
                          knots=list(day_of_year=c(0.5, 366.5), day_of_week=c(0.5, 7.5)),
                          data=daily_ts,
                          method="REML",
                          family=nb())

By accident I changed the order in equation to:

library(mgcv)

gam_model <- gam(counts ~ te(year, day_of_year, k=c(20,12), bs=c("tp","cc") +
                          te(year, day_of_week, k=c(20,7), bs=c("tp", "cc"),
                          knots=list(day_of_year=c(0.5, 366.5), day_of_week=c(0.5, 7.5)),
                          data=daily_ts,
                          method="REML",
                          family=nb())

And results are different, for example the tensor for year and day_of_year changes from (1st eq): enter image description here

to (2nd eq): enter image description here

Why this situation happen? I've checked for concurvity, there's none. I thought that the order of equation's components doesn't matter. Which solution I should choose or how to correct it?

$\endgroup$

1 Answer 1

9
$\begingroup$

The order shouldn't matter, but your model is not well specified.

I suspect it's because the main effects of year are in both smooths and hence you have the same effects in the model twice. mgcv tries to fix this by identifying these redundant columns from the model matrix and I can well imagine that this process is dependent upon the ordering of the columns in the model matrix.

You should refit your model using ti() terms as follows

m <- gam(counts ~ ti(year, k = 20, bs = "tp") +
    ti(day_of_week, k = 7, bs = "cc") +
    ti(day_of_year, k = 12, bs = "cc") +
    ti(year, day_of_week, k = c(10, 7), bs = c("tp", "cc")) +
    ti(year, day_of_year, k = c(10, 12), bs = c("tp", "cc")),
  knots = list(day_of_year = c(0.5, 366.5),
    day_of_week = c(0.5, 7.5)),
  data = daily_ts,
  method = "REML",
  family = nb())

(you can use s() for the univariate terms, but you might want to be explicit about the basis type there as s() uses "tp" by default and ti() use "cr" and I'm not sure how much this is an issue or not.)

The point of the ti() bivariate terms is that these smooths are pure interactions, having had the main effects removed from the basis. With this formulation, the model has three smooths for the main effects, and then two bivariate tensor product interactions that only include the smooth interaction between pairs of covariates. The only downside is ghat this model form requires more smoothing parameters to be estimated.

$\endgroup$
2
  • $\begingroup$ That explains a lot, thank you. Let me ask what should I do in case when I want to introduce two variables with interaction effect into "by" arguments for each smooth? I want to look for different effects for this variables. I specify fixed effects and then add interaction into each "by" in ti()'s components right? m <- gam(counts ~ age + sex + ti(year, by=interaction(age, sex)(...) $\endgroup$
    – Tom
    Commented Jan 19 at 10:40
  • 2
    $\begingroup$ I think you should probably ask that as a separate question, as I'm not 100% clear on what you are asking just from your comment and the snippet of code $\endgroup$ Commented Jan 19 at 14:41

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.