GAM (mgcv): question on trend

Question

I am new to the exchange and self-teaching myself GAM (via mgcv package).

Context: I would like to estimate a panel regression via gam.

1. periodicity: monthly
2. balanced panel

Question: may I estimate a GAM to include a global trend and local trend and random? It would look like this:

model_1 = gam(Y ~ X + s(LOCATION, bs="re") + s(time, by=COUNTRY),
              data= dat, correlation = corAR1(form = ~ 1 | LOCATION), 
              method = "REML",...) 

The first smooth term (i.e., s(LOCATION, bs="re")) is a random intercept. The second is a non-parametric trend by location. Is this model appropriate, as I would like to model the panel in this way? I am intentionally simplifying in hopes it does strips away unneeded complexity. Very grateful for any advice!

Answer 1

You don't really provide much to go with here. For example you say the periodicity is monthly but there's nothing in the model to account for any seasonal effect, and you conflate LOCATION, COUNTRY, and "location". Also, what's X?

If LOCATION indexes the subjects that are measured over time, you could imagine having the smooth equivalent of random trends (slopes) for each LOCATION rather than having them grouped at the COUNTRY level as you have here.

gamm(Y ~ X + s(time) + s(time, LOCATION, bs = 'fs')
     data = dat, correlation = corAR1(form = ~ 1 | LOCATION), 
     method = "REML", ...)

Would give a similar model to your's, but with smooth trends for each LOCATION instead of COUNTRY, with the s(time) term representing the average trend over all observations. The random intercept is contained in the random smooth of time for LOCATIONs, so you don't need a separate random effect smooth just for LOCATION.

If you are going to use a factor by term then you should also include COUNTRY (the factor) in the model as a fixed (or random effect), so i would suggest you add + COUNTRY to your model if you want smooth trends at the COUNTRY level.

Also, to fit models with correlation structures via argument correlation and using corAR1() you have to use the gamm() function. You can fit a nested AR(1) with bam() but you have to specify the $\rho$ you want, as it won't be estimated.

In the corAR1() you should perhaps consider passing it an integer time variable rather than relying on the order of the data in dat.