Scale of time covariate in `corCAR1()` matters?

Question

I'm running a bunch of GAM analyses on time series data (measurements of my own weight). Unlike many examples I see online, my data is not spaced evenly in time. In fact, time points can come minutes (if not seconds) apart, and the entire range of the data is about a year (currently about 500 data points).

Taking the helpful advice of Gavin Simpson's blog posts on the matter, I am looking to make my model aware of the autocorrelation in the data, using something like:

m <- gamm(Y ~ s(time) + s(...) + ..., correlation = corCAR1(form = ~time))

However, after struggling for a long time trying to make sense of things, I realized that the results of my model are wildly different depending on the scale of the time covariate.

I initially just converted the datetimes into EPOCH time, which resulted in very large numbers---for example 1569201817. Using this covariate, the AR1 component of my model did little to nothing.

But when I scaled the same variable (ie. as.vector(scale(time)), after converting it into a numeric), the model was completely different, with the AR1 component playing a major role.

Playing around with it, I found that whenever I scale the range of the time covariate to greater than 70, the behavior changes to disfavoring the AR1 model.

What gives? The example the corCAR1() docs use has a time covariate that ranges from -0.16 to 1.16, so that suggests I should scale my data, but the fact that everything changes with an arbitrary scale doesn't make much sense to me.

How do I know the proper scale to use?

Answer 1

I'm assuming that s(time) or something like it in the model? If so, you can get into situations where both the smooth of time and the CAR(1) process, which are mathematically similar, cannot both be uniquely identified; one of the two ends up winning out, but which one can depend on a number of things. As you've found out it can depend on the scale of the values used for time and how these get handled by the floating point arithmetic going on in the CPU. I've seen wiggly smooths with effectively zero CAR(1) estimated on one CPU and the opposite model estimated on another CPU.

To diagnose if this is the problem, I would look at how the s() terms differ between the two versions of the model, and look at intervals(m$lme) to see if there may be any problems with either fit. With the latter you will often you will get a warning:

Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance

which is a clear indication that there is something wrong with the model, usually that the model you are trying to fit is too complex for the data.