# Implementation of continuous AR(1) in nlme package

### 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
$$$$

• This is defined in Pinheiro and Bates (2000) Eq. 5.27: $h(s, \phi) = \phi^s, ~s\geq0,~\phi\geq0$. It is analogous to discrete AR(1), only the lag $s$ is now continuous. They cite Longitudinal Data with Serial Correlation: A Statespace Approach (Jones, 1993). If you need the actual implementation (i.e., how the correlation matrix is constructed), you probably need to study the source code, i.e. corStruct.c (search for CAR1_mat). Aug 19, 2020 at 6:42
• hmm... I got a hold of the book and checked this section. Do you @Roland have any intuition on how one would simulate these CAR(1) residuals? Whats the continuous-time version of the AR(1)'s $y_i = \phi y_{i-1} + \epsilon$? I'm assuming the variance of $\epsilon$ needs to vary with respect to the lag value in the continuous time version. Aug 19, 2020 at 12:31
• I have never done this but chapter 3.1 of this paper should be useful: ncbi.nlm.nih.gov/pmc/articles/PMC5655034 Aug 19, 2020 at 12:57
• Great! In the section 3.1 you referenced: Note that the CAR(1) model is also known as the Ornstein-Uhlenbeck process or model. Exactly what I was looking for! So it seems I'm already familiar with CAR(1) models, only knowing them by a different name. Aug 19, 2020 at 13:06
• If you implement the simulation as an R function, please make it available (e.g., as an answer here). Aug 19, 2020 at 13:32

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 = 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.