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.

  • $\begingroup$ Would weighting the rho values in ri,t by the length of an irregular interval of time adequately represent a decay function? Isn't the main concern the distribution of the impact of errors over a time span? $\endgroup$ Dec 1, 2020 at 22:20
  • $\begingroup$ Could you use a Latent variable Gaussian process? $\endgroup$
    – JTH
    Dec 2, 2020 at 4:35
  • 2
    $\begingroup$ Matérn covariance kernels with shape $\nu = p - 1/2$ with $p = 1$, $2$, $\dots$ define some special AR(p) which work fine in continuous time, via Kalman filering. The case $p=1$ (hence $\nu = 1/2$) corresponds to the Ornstein-Ulhenbeck process. The Kalman filter is used in discrete-time, looping over the ordered observation times. For all AR odrers $p$ only one parameter is used. $\endgroup$
    – Yves
    Dec 2, 2020 at 9:44

3 Answers 3


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}).$

  • 1
    $\begingroup$ This is very helpful and will address many modeling problems including prediction. But for my particular problem, which is comparing two treatments in a randomized clinical trial, I require an approach that provides causal inference on the effect of treatment. For that purpose, the multivariate process has to be a marginal one with respect to $Y(t)$, i.e., $Y(t)$ cannot condition on $Y(s)$ where $s < t$ since that would be conditioning on an earlier treatment effect. $\endgroup$ Dec 2, 2020 at 10:56
  • $\begingroup$ I make be mistaken about that. The $Y(t)$ refer to residuals in my model, not to raw response values. I think it is OK for residuals to depend on earlier residuals. $\endgroup$ Dec 2, 2020 at 13:34
  • $\begingroup$ Is there a form for which the $Y(t + \Delta t)$ distribution can ignore $Y(t)$? I'd like to be able to deal with "tall and thin" datasets where information from the current record (one time for one subject) can stand alone. $\endgroup$ Dec 2, 2020 at 13:53
  • $\begingroup$ On the other hand I can just add a lagged time to my dataset. $\endgroup$ Dec 2, 2020 at 14:46
  • 1
    $\begingroup$ Please see my newest post in the Stan forum (discourse.mc-stan.org/t/…) $\endgroup$
    – TrungDung
    Dec 12, 2020 at 23:20

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





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


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


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


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


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

  • 2
    $\begingroup$ This is a fascinating approach that I had not seen before. The one downside seems to be a possible need for more control over the correlation as a function of the time lag, and I do favor isotropic structures. $\endgroup$ Dec 2, 2020 at 10:52
  • $\begingroup$ I think we could use alternate decaying structures like $e^{-a \delta t}$ for the correlation to get more control over decay if this one is too steep. Could you please let me know what you meant by isotropic structures in this context? $\endgroup$ Dec 2, 2020 at 15:45
  • 1
    $\begingroup$ Isotropic = correlation is a function of $|t_{i} - t_{j}|$ and not the individual times. $\endgroup$ Dec 2, 2020 at 16:15
  • 1
    $\begingroup$ Hmm.. I guess then the only way is to get true isotropic structures is to sample the vector of random effects for ith individual $r_{ij}$'s from $N(\boldsymbol{0},\boldsymbol(\Sigma(\delta T_i))$. Which would then be similar to your $N \times T$ structure but not a matrix as we have different observations. Alternatively, we can try to find a multiplication factor to reduce the non-isotropy of the correlation as low as possible as a tradeoff to avoid sampling a vector. $\endgroup$ Dec 2, 2020 at 17:56

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))

#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
  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',
  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
  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

#fit with integration (faster, linearised approximation)
ro <- ctStanFit( datalong = d,
  ctstanmodel = m,cores=cores,
  intoverstates = T,nopriors=F,
# 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

#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,
  chains = cores,plot=F,
  intoverstates = FALSE,
# s
  • $\begingroup$ Very nice work and R package. My understanding is that the continuous time AR(1) generalization you are using requires lots of parameters because of lots of possible $\Delta t$ values. Is there an alternative formulate that is an even more continuous AR(1) with few parameters? $\endgroup$ Dec 3, 2020 at 13:54
  • 1
    $\begingroup$ The number of $\Delta t$ values does not influence the number of parameters, no -- for a univariate AR1 process there is just the single 'auto effect' parameter, in this case a 1x1 matrix 'A'. The discrete time covariance with an earlier time point is then $e^{A \Delta t}$ $\endgroup$
    – Charlie
    Dec 3, 2020 at 14:00
  • $\begingroup$ Great. Is there a fundamental reference for this continuous time AR(1) formulation? $\endgroup$ Dec 3, 2020 at 14:16
  • 1
    $\begingroup$ For stochastic differential equation fundamentals, Gardiner, C. W. (1985). Handbook of stochastic methods (Vol. 3, pp. 2-20). Berlin: springer. For the most up to date / best explained version of what I'm specifically using, Driver, C. C., & Voelkle, M. C. (2018). Hierarchical Bayesian continuous time dynamic modeling. Psychological Methods, 23(4), 774. $\endgroup$
    – Charlie
    Dec 3, 2020 at 15:41
  • 1
    $\begingroup$ Oh - then I misunderstood. That would be expected , so far as I follow -- if you nonlinearly transform a variable that has a certain linear relationship originally, the resulting linear relationship will be different, and only capture part of the 'true' relationship that was in place before the transform. $\endgroup$
    – Charlie
    Dec 23, 2020 at 17:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.