Implementation of continuous AR(1) in nlme package

Question

Background

I used the nlme package to implement a mixed effects model with an AR(1) correlation structure on the residuals.

In nlme, the grouping structure of the AR(1) process must be the same as the random effects structure, as far as I understand. While my AR(1) process was in fact discrete-time, because I had missing data within my random effects groups, I had to instead implement this as a continuous-time AR(1) model (?corCAR1). To my knowledge, I cannot simply specify that there are multiple AR(1) groups nested within my random effect groupings.

Question

I've calculated some confidence intervals, and I wanted to make a short parametric bootstrap to make sure that my math was correct (my calculated confidence bands from the models covariance matrix agree with the parametric bootstrapped confidence bands).

Trying to simulate new datasets, I realised it's not clear to me exactly what the continuous-time AR(1) process is. help(corCAR1) doesn't give any mathematical details, but only references to three textbooks.

• Does anyone know how continuous-time AR(1) process are implemented in the nlme package? I'm aware that an Ornstein-Uhlenbeck process is considered an continuous time analogue of the AR(1) process... but I'm not sure if there are other implementations I'm not aware of which nlme may be using.

Extra info

In case more information is needed, below I've put the model summary.

Linear mixed-effects model fit by REML
 Data: data 
        AIC       BIC   logLik
  -6472.216 -6415.398 3244.108

Random effects:
 Formula: ~1 | fish
        (Intercept)  Residual
StdDev:  0.09100886 0.2251129

Correlation Structure: Continuous AR(1)
 Formula: ~time.2 | fish 
 Parameter estimate(s):
      Phi 
0.6738598 
Fixed effects: distance_log ~ period + current_speed + tide_m 
                   Value  Std.Error   DF  t-value p-value
(Intercept)    1.2148390 0.02487684 8957 48.83414  0.0000
periodDuring   0.0228981 0.01208525 8957  1.89471  0.0582
periodAfter    0.0519890 0.01285010 8957  4.04581  0.0001
current_speed  0.1018782 0.01850602 8957  5.50514  0.0000
tide_m        -0.0088376 0.00490187 8957 -1.80290  0.0714
 Correlation: 
              (Intr) prdDrn prdAft crrnt_
periodDuring  -0.292                     
periodAfter   -0.249  0.451              
current_speed -0.432  0.178  0.092       
tide_m         0.083 -0.054  0.015 -0.197

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-6.30361758 -0.56513258 -0.01809098  0.46685852  4.81926579 

Number of Observations: 8980
Number of Groups: 19 
```

Answer 1

Thanks to a paper suggestion (section 3.1) by @Roland, I realised that CAR(1) processes are in fact Ornstein-Uhlenbeck (OU) processes.

For my question, I wanted to know how to simulate the CAR(1) process, so I could run a parametric bootstrap on a model which fitted CAR(1) residuals.

Theory

With a conventional discrete-time AR(1), we're used to seeing simulations being made with the following formula: $$x_t = \mu + \phi (x_{t-1} - \mu) + \epsilon ,$$ where the step lengths follow a Gaussian distribution $\epsilon = N(0, \sigma^2)$. $\phi$ is our auto correlation parameter and $\mu$ is the drift.

Johnson (2008) lays out the OU analogue to this with $$ x_{t} = \mu + e^{-\beta \Delta} \left[ x_{t-1} - \mu \right] +\zeta .$$ In this case we have a Gaussian step length distribution which has a variance dependant on our time step value $\Delta$ where $\zeta = N \left(0, \sigma^2 \left[ 1 - e^{-2 \beta \Delta} \right] / 2 \beta \right)$. The correlation coefficient of our OU process is $\phi = e^{-\beta}$, assuming unit time steps ($\Delta = 1$) so the coefficient is comparable to the AR(1) version.

The wikipedia page for OU processes doesn't explicitly show this formulation. However, you can derive it from the provided formulas for $E(x_t)$ and $\text{cov}(x_s, x_t)$ there.

Code

Now we'll code up that OU process shown above and make some example simulated data.

###### Helper function
OU <- function(x, mu, delta, beta, sigma){
  # A function for simulating the next step in an OU (continuous-time AR(1) 
  # process)
  
  # x: Current value
  # mu: drift parameter
  # delta: time step
  # beta: correlation coefficient
  # sigma: variance of the step length distribution
  
  sigma_scaled <- (sigma^2) * (1 - exp(-2*beta*delta))/(2*beta)
  
  x2 <- mu + exp(-beta*delta)*(x - mu) + 
    rnorm(mean = 0, sd = sqrt(sigma_scaled), n = 1)
}

############ Simulate CAR(1) process with low correlation
reps <- 100
mu <- 0
sigma <- 1
phi <- 0.8
beta <- -log(phi)

# simulate random time intervals between steps
delta <- runif(min = 0, max = 1, n = reps)
# Uncomment below to make discreete time equivalent
#delta <- rep.int(1, times = reps)  


# Init position vectors
x <- vector(mode = 'numeric', length = reps); x[1] = 0

# loop and fill positions
for(i in 2:reps){
  x[i] <- OU(x[i - 1], mu, delta[i], beta, sigma)
  y[i] <- OU(x[i - 1], mu, delta[i], beta, sigma)
}

# Show result
plot(x, type = 'line')

# For uniform time intervals, run this to calc correlation coefficient
acf(x)

The results of the simulations are shown in below.

When we set our time steps to unit length, the OU process is equivalent to an AR(1) process in discrete time, as seen in the resulting ACF plot below.

References

• Johnson, D. S., London, J. M., Lea, M.-A. and Durban, J. W. (2008). CONTINUOUS‐TIME CORRELATED RANDOM WALK MODEL FOR ANIMAL TELEMETRY DATA. Ecology 89, 1208–1215.