I have univariate time series data that I am fitting a gamm to. I am using 10-fold cross validation to choose the best temporal correlation structure in the data. The lowest RMSE value tends to come from an overly complicated ARMA structure (e.g.,(3,2)). I believe the reason may be that with lower order ARMA structures the day covariate smoother becomes overly complex (I think perhaps it is "soaking" in some of the correlation, for lack of a better term). I know a-priori that there is little variability in the data across days, so I think it may help to either fix the k value low (I am thinking 3), or set the max k low for the day covariate. I am not sure how to pick the specific k value, or if this sounds reasonable. Any help is much appreciated, thank-you!

The code I used for all models is shown below, the only change I made was varying the correlation structure for values of p=0-3, and q=0-3.

gam=gamm(sv~s(day,bs="tp")+s(range,bs="tp")+s(time,bs="cc"),data=data.all[1:336,],correlation=corARMA(p=3,q=2),gamma=1.4,control=list(maxIter = 10000, msMaxIter=10000,niterEM=0,msMaxEval=10000))

2 example results from the GAMMs are shown below, note: "day" is a Julian day count for a single month, "range" is the tidal range values, and "time" is a variable with values 1-24 for time of day. The results for the gamm with an AR-1 structure is shown below: enter image description here

The results from the ARMA 3-2 (which produced a lower RMSE) is shown below: enter image description here


If you want to fix the smoothness of terms in the model, choose k = 3 and set fx = TRUE. (With a thin plate spline k = 3 in the smallest you can go down to for the basis dimension.)

Alternatively you can turn day into a linear term by just allowing it to enter as day and not s(day) in the model formula.

Thinking about this, I'm not sure I would approach the problem they way you are doing: the correlation structure here isn't improving the fit of the model. Why don't you just using AIC or BIC and fit ARMA models to the residuals of a GAM that you are happy with (having fixed whatever terms you want to avoid overfitting). You can do this using auto.arima() in the forecast package. Then take whatever ARMA model it selects and refit your gamm() with this order of ARMA terms.

In refitting the model to different subsets of the data, what you are cross-validating is both the model fit and the correlation structure. If your interest in in choosing the latter, limit the complexity of the spline terms until you are happy with the fit and instead fit ARMA models to the residuals of this model. That way the data you are fitting isn't changing under you all the time.

| cite | improve this answer | |
  • 1
    $\begingroup$ @Hannah Ping me on email again if you want to chat about the latter couple of paragraphs. (sorry I haven't yet replied to your last email - I've been away on and off the past few weeks and nursing a poorly toddler [fine now, having passed on her cold to me...]) $\endgroup$ – Gavin Simpson Oct 28 '15 at 18:13
  • $\begingroup$ Sounds good, I totally understand! I hope the cold clears up.. $\endgroup$ – Hannah Oct 28 '15 at 18:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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