Serial correlation AR(1) model for residuals: how to generalize to irregular times

Question

I am working on a Bayesian serial correlation model for binary and ordinal logistic models (proportional odds model). I am modeling the serial correlation structure on the random effects of the model (on the logit scale) so that it is easy to handle ordinal outcome variables. If all subjects are measured at regular, say integer-valued, times, the following thinking works:

• Let $T$ be the integer maximum observation time over all subjects
• Suppose that observation times are $t=1, \dots, T$
• Let $\gamma_i$ be a $n(0, \sigma_\gamma$) random effect for the $i$th subject.
• Let $\epsilon_{i,1}, ... \epsilon_{i,T}$ be the within-subject white noise for the $i$th subject that is $n(0, \sigma_w)$, where $T$ is the maximum follow-up time (we may only use the first few of these for a given subject)
• $\epsilon_{i,1} = \gamma_i$ so the usual random effect as generated for a hierarchical (compound symmetric correlation pattern) repeated measures model is the starting white noise for a given subject (this random effect may have a different standard deviation $\sigma_\gamma$)
• Then the random effect for subject $i$ at time $t$ is $r_{i,t} = \rho^{t}\gamma_i + \rho^{t-1}\epsilon_{i,1} + \rho^{t-2}\epsilon_{i,2} + ... \epsilon_{i,T}$
• The random effects must all be defined regardless of which observations are actually observed, so we can re-write the model as using a matrix of random effects $r_{i,1} = \gamma_i, r_{i,2} = \rho r_{i,1} + \epsilon_{i,2}, r_{i,3} = \rho r_{i, 2} + \epsilon_{i,3}, ...$.

I want to generalize this to an irregular-spaced continuous time AR(1) structure. For that purpose I would rather not envision an $N \times T$ white noise (or random effects) matrix but would like to develop an irregular continuous-time AR(1) structure or something that is similar to AR(1) in inducing higher correlation between two nearby measurements when compared to the correlation between two distant measurements. Thoughts for how to envision and develop this would be welcomed.

Answer 1

The Ornstein-Uhlenbeck (OU) process can be considered as it is a natural extension of the (V)AR(1).

I here skip the mathematical derivation. I first present the multivariate version, which I worked with. The OU process can be specified as follows:

$Y(0)\sim N(\mu, \Omega)$,

$Y(t+\Delta t)|Y(t) \sim N(\mu + e^{-\Gamma \Delta t}\left(Y(t)-\mu\right), \Omega-e^{-\Gamma \Delta t}\Omega e^{-\Gamma^T \Delta t}),$

where $\mu, \Omega, \Gamma$ are the parameters satisfying the following constraints:

• The real part of each eigenvalue of $\Gamma$ is positive.
• $\Gamma\Omega+\Omega\Gamma^T$ is a covariance matrix,
• $\Omega$ is a covariance matrix.

You can see now that if we reduce to one response:

$Y(0)\sim N(\mu, \omega^2)$,

$Y(t+\Delta t)|Y(t) \sim N(\mu + e^{-\gamma \Delta t}\left(Y(t)-\mu\right), \omega^2\times (1-e^{-2\gamma \Delta t}).$

Now we consider equidistant time points, the process is reduced to AR(1):

$Y(0)\sim N(\mu, \omega^2)$,

$Y(t+1)|Y(t) \sim N(\mu + e^{-\gamma}\left(Y(t)-\mu\right), \omega^2\times (1-e^{-2\gamma}).$

Finally, if we set $\mu$ equal to 0, then we have:

$Y(0)\sim N(0, \omega^2)$,

$Y(t+1)|Y(t) \sim N(e^{-\gamma}Y(t), \omega^2\times (1-e^{-2\gamma}).$

Answer 2

Assuming the data to be at irregular intervals with different number of measurements for each subject:

$y_{t_{11}},y_{t_{12}},\ldots,y_{t_{1K_1}}$

$y_{t_{21}},y_{t_{22}},\ldots,y_{t_{2K_2}}$

$\vdots$

$y_{t_{N1}},y_{t_{N2}},\ldots,y_{t_{NK_N}}$

For sake of simplicity, I am writing the equations for linear models instead of a proportional odds model which can be modeled using a logit link function.

We could model the first measurement for the ith subject as:

$$y_{t_{i1}} = \beta_0 +\beta_1t_{i1}+\boldsymbol{\beta X}+\gamma_i+\epsilon_{i1}$$

The second measurement should be autocorrelated with the first measurement. The correlation should be 1 when $t_{i2} = t_{i1}$ and 0 when they are too far away. One choice is to choose the correlation to be inversely proportional to the time interval. $$corr(t1,t2) = \frac{\rho}{\rho+f(t_{i2}-t_{i1})} $$ where $f$ can be chosen depending on the decay needed. If $f$ is chosen as the identity function, it means that we have a correlation of 0.5 at a time interval of $\rho$ (which is flexible and $\rho$ can be estimated accordingly)

The second measurement then becomes:

$$y_{t_{i2}} = \beta_0 +\beta_1t_{i2}+\boldsymbol{\beta X}+\frac{\rho}{\rho+(t_{i2}-t_{i1})}\gamma_i+\epsilon_{i2}$$

Similarly,

$$y_{t_{i3}} = \beta_0 +\beta_1t_{i3}+\boldsymbol{\beta X}+\frac{\rho}{\rho+(t_{i3}-t_{i1})}\gamma_i+\epsilon_{i3}$$

$$\vdots$$

$$y_{t_{iK_i}} = \beta_0 +\beta_1t_{iK_i}+\boldsymbol{\beta X}+\frac{\rho}{\rho+(t_{iK_i}-t_{i1})}\gamma_i+\epsilon_{iK_i}$$

This correlation structure seems simple and natural, however, (I think) as we go to the measurements at the end, like the correlation between the last observation and last but one observation is small compared to the first and second observation even if they are equally separated in time.

To remedy this an alternative specification of the random effects to preserve autocorrelation between measurements is:

$$y_{t_{i2}} = \beta_0 +\beta_1t_{i2}+\boldsymbol{\beta X}+\frac{\rho}{\rho+(t_{i2}-t_{i1})}\gamma_i+\epsilon_{i2}$$

Similarly,

$$y_{t_{i3}} = \beta_0 +\beta_1t_{i3}+\boldsymbol{\beta X}+\frac{\rho}{\rho+(t_{i2}-t_{i1})}\frac{\rho}{\rho+(t_{i3}-t_{i2})}\gamma_i+\epsilon_{i3}$$

$$\vdots$$

$$y_{t_{iK_i}} = \beta_0 +\beta_1t_{iK_i}+\boldsymbol{\beta X}+\frac{\rho}{\rho+(t_{i2}-t_{i1})}\frac{\rho}{\rho+(t_{i3}-t_{i2})} \ldots \frac{\rho}{\rho+(t_{iK_i}-t_{iK_{i-1}})}\gamma_i+\epsilon_{iK_i}$$

There may be some attenuation in correlation here too, but probably much less than the previous specification.

This model could be estimated with a Bayesian hierarchical model with appropriate priors on the random effects and beta coefficients without any non-identifiability issues.

I just wrote down my thoughts, I am not an expert in the field. Please let me know if I missed any major concept or if am wrong somewhere.

Answer 3

When you generalize to a continuous time form, you only need the random effects relevant to your specific observations, because the covariance with earlier random effects will be conditioned on the time interval. This paper describes a continuous-time Rasch model https://psycnet.apa.org/record/2019-22131-001 , using the ctsem software https://cran.r-project.org/web/packages/ctsem/index.html for R (I am the author) . The software is designed to handle this sort of structure, though only limited development (and more limited documentation) have gone into the non-continuous data side of things and once multivariate situations are encountered, the sampling is very slow.

Here is R code that generates and fits what I believe to be the sort of structure you're talking about, using the much faster extended kalman filter approximation, as well as the more rigorous full sampling approach via Stan's HMC:

#install software
# source(file = 'https://github.com/cdriveraus/ctsem/raw/master/installctsem.R')

invlogit=function (x) exp(x)/(1 + exp(x))
n.manifest=7

#generate data
gm <- ctModel(DRIFT=-.3, DIFFUSION=.3, CINT=.1, #dynamic system pars
  TRAITVAR=diag(.3,1), #old approach to allow individual variation 
  LAMBDA= matrix(rep(1,each=n.manifest)), #factor loading
  TDpredNames = 'intervention',
  TDPREDMEANS = matrix(c(rep(0,10),1,rep(0,9))), TDPREDEFFECT = 1, #intervention timing / effect
  Tpoints=20,
  MANIFESTMEANS=c(0,rep(c(.5,-.5),each=(n.manifest-1)/2)), #measurement offset
  T0MEANS=-.3, #initial latent state
  T0VAR=.5 #initial latent variance
  )

d=ctGenerate(gm,n.subjects = 50,logdtsd=.2) #generate continuous data
#convert to binary
d[,gm$manifestNames] <- rbinom(nrow(d)*gm$n.manifest,size=1,prob=invlogit(d[,gm$manifestNames]))

#model to fit
m <- ctModel(
  manifestNames = gm$manifestNames, #observation variables
  TDpredNames = 'intervention',
  latentNames='eta1',
  MANIFESTMEANS = c(0,paste0('m',2:n.manifest,'|param|FALSE')), #set prior to N(0,1), disable individual variation
  LAMBDA = rep(1,n.manifest), #factor loading
  T0MEANS='t0m|param|TRUE|.1',
  CINT = 'b|param|TRUE|1', #use standard normal for mean prior, individual variation = TRUE (default), default scale for sd
  type = "stanct" )

# ctModelLatex(m) #shows general SDE / measurement structure but assuming continuous observation type

m$manifesttype[]=1 #set observation type to binary
cores=6



#fit with integration (faster, linearised approximation)
ro <- ctStanFit( datalong = d,
  ctstanmodel = m,cores=cores,
  plot=10,verbose=0,
  intoverstates = T,nopriors=F,
  optimize=T,intoverpop=T)#,optimcontrol=list(stochastic=F))
so=summary(ro)
# so

ctKalman(ro,plot=T,kalmanvec='etasmooth',subjects=1:3) #latent performance, conditioned on parameters and all data
ctKalman(ro,plot=T,kalmanvec='etaprior')#conditioned only on pars and previous data
ctKalman(ro,plot=T,kalmanvec='yprior') #observation predictions
ctKalman(ro,plot=T,kalmanvec='ysmooth')




#fit without kalman filter integration (much slower, using Stan's HMC sampler)
r <- ctStanFit( datalong = d,
  #fit=FALSE, #set this to skip fitting and just get the standata and stanmodel objects
  ctstanmodel = m,
  iter = 200,verbose=0,
  control=list(max_treedepth=4),
  nopriors=FALSE,
  chains = cores,plot=F,
  intoverstates = FALSE,
  optimize=FALSE,intoverpop=F)
s=summary(r)
# s