9

I'm interested in modelling a time series of temperature data across several years. The data are on the level of hourly observations, so I have variables for year, month, day, and time.

I found a great example of doing this by Gavin Simpson (found here). The blog only considers correlation within year, where as I have to deal with correlation within year and within day.

How can I best account for this correlation with gamm? Gavin uses the following code

modar2 <- gamm(apparentTemperature ~ s(month, bs = "cc", k = 12) + s(time, k = 20),data = timetemp, correlation = corARMA(form = ~ 1|year, p = 2),control = ctrl)

Where should I pass variables to account for correlation within day?

For reference, here is a sample of my data:

tibble::tribble(
               ~created_at,            ~time, ~month, ~year, 
~apparentTemperature,
     "2014-01-03 09:30:28",              9.5,      1,  2014,               -17.87,
     "2014-01-03 10:13:43", 10.2166666666667,      1,  2014,               -17.87,
     "2014-01-03 12:19:32", 12.3166666666667,      1,  2014,               -16.14,
     "2014-01-03 12:44:04", 12.7333333333333,      1,  2014,               -20.24,
     "2014-01-03 13:09:38",            13.15,      1,  2014,               -20.24,
     "2014-01-03 13:39:00",            13.65,      1,  2014,               -20.44
     )
Demetri Pananos
  • 36,121

1 Answers1

6

Depends how you want nest the autocorrelation, within days?

modar2 <- gamm(apparentTemperature ~ s(year) + s(month, bs = "cc", k = 12) +
               s(time, k = 20), data = timetemp,
               correlation = corARMA(form = ~ 1|day, p = 2),
               control = ctrl)

would have smooth long term trend, smooth seasonal effect, smooth time of day effect, with autocorrelation nested within days (for which you'd need to create a new variable day which generates the day of year from the date time variable.

If you have a lot of data, you really don't want to use form = ~ obs_seq for the correlation structure, where obs_seq is a sequence 1, 2, ..., number of observations, as that will create a massive covariance matrix that lme() will need to invert at each iteration. Having fitted such a model to high frequency data, it took gamm() a week to converge on powerful multicore workstation.

The reason I nested the correlation within year in that example was partly for this reason; that's a long monthly record and fitting a full ARMA function across all timepoints is not quick.

Gavin Simpson
  • 47,626
  • 1
    Hmm, my data is nearly 30,000 observations, so I guess that is far too large. I could summarize my data by day (e.g. mean apparent temperature for April 18 2018) which reduces my data to nearly 1,500 observations. Then I could just use your approach from the blog, yes? I'd have Year, month, and day of month, smoothing year and day with a seasonal effect of month. – Demetri Pananos Apr 20 '18 at 02:10
  • 1
    Why do you want to have the autocorrelation operate at beyond the daily level? Do you have evidence of longer scale autocorrelation? The main issue is that you're likely to run out of RAM unless you have a lot of it available. – Gavin Simpson Apr 20 '18 at 02:26
  • 1
    Start with the simpler model (nest the AR within day-of-year within year), use an AR(1) initially (corAR1()) rather than an ARMA as that is much more efficient. Then if that fits, look at the normalized residuals to see if you still have remaining autocorrelation. You could also fit without the AR and check that model's residuals. If you go that route, see bam() in mgcv. – Gavin Simpson Apr 20 '18 at 03:09