R package ctsem - hierarchical continuous time dynamic models

Publicity
1

Hi all. As part of my phd I’ve been looking into dynamic models within developmental psychology, and working on some software to make it a bit easier for people to deal with. The software sets up a continuous (or discrete) time hierarchical state space model, using the Kalman filter (by default) to integrate over latent states. There is still some work needed, but it’s fairly functional as it stands, though definitely fails some of the posted guidelines - I’m not sure if they are all so applicable (avoiding compilation time when the models take much longer to fit anyway doesn’t seem so crucial) but will take a better look and change what I can sometime soon.

An intro to the package can be found at http://github.com/cdriveraus/ctsem/raw/master/vignettes/hierarchical.pdf
To get the very quick overview, looking at / running the sunspots example included in the intro vignette, starting bottom of page 11, might be simplest. If anyone is particularly interested I can also send through the soon to be submitted paper. I would of course be happy for any comments, on any aspect…

btw, the CRAN version is still only frequentist mixed effects.

5 Likes
2

@Charles_Driver
Hi, I am working on a time-series analysis in which the time intervals are not equal and the number of observations are not large for about 1200 subjects.
My outcome and predictors both are ordinal not continuous. I tried a syntax in Stan and shared it here. Then I realized I need to add time as a continuous dynamic and end up reading your paper. I am looking also to obtain individual-based estimates for the autoregressive and lagged models. I appreciate it if you can add your input or guide me to incorporate continuous time. I am also reading your package, but to familiarize myself with all required maths and address ordinal response and predictors.

https://discourse.mc-stan.org/t/time-series-in-stan-i-am-new-to-stan-and-need-hints-to-develop-the-model-thanks

Thanks

1 Like
3

I’m not sure whether to respond there or here, but there is too much going on in that thread for me to easily follow. Happy to provide some input but give me some bite sized questions not a sprawling labyrinth ;) from what I could follow, you have 24 ordinal measures at 4 time points with variable intervals and you want to model these. I don’t think you stated the goal anywhere – what are your trying to learn? If you just want to forecast then it may not be too bad to ignore the time differences, depending on how large they are. If you are really interested in inference on the process, then yes you’ll need to be more serious (and hopeful!). If it were me, I would be tempted to try approximate the ordinal likelihood as some transformation of a gaussian, and approximate this using the extended Kalman filter in ctsem. At least, I never had much luck getting Stan to sample the hidden states in coupled multivariate processes efficiently, integrating out the states via kalman filters was always more effective. ctsem will also allow you to generate stan code that is very close for the state sampling / proper ordinal approach, but that code might be a bit wild :)

Since the initial post is a few years old by now:
For a few examples see: https://cdriver.netlify.app/

4

Hello,
Thanks for your response. I plan to know the inference on the process than forecasting. I also need to obtain an individual-based estimate.

“I would be tempted to try to approximate the ordinal likelihood as some transformation of a gaussian,” we tried this in stan, indeed. I will try ctsem and get also stan code. What is the syntax to get stan code from ctstem?
I may get back to you here if I encounter questions.
Thanks

5

what is an individual based estimate – you want to estimate unique dynamic systems and measurement models for every subject? or you just want estimates of all the latent states returned, based on an overall system / measurement model? If you really want every parameter in the system to vary by subject, you will need very tight priors to get anywhere with 4 time points, I expect. I would start by getting estimates out of the simplest sensible model you can work with and building up with lots of sanity checks.

the stan code behind a fit object can always be found using
myfit$stanmodeltext
and you can create a skeleton fit object without actually fitting, by using the argument fit=FALSE .

there is a section on nonlinear filtering in this doc (but start linear!):

6

This a great advice. I am working on a data that not only the time-intervals are varied between visits but also the number of measurement times are varied. I chose 4 time points as everyone has at least 4 time-point measurements which still can be a good representative with more sample size.
As there are individual-based variations in the response and predictors I am hoping to get person-based estimates in addition to population-based estimates. I am also interested in the latent values. What I tried so far has been more focused on population-based estimate, and I got one estimate per time-series. Then I thought I should add the coefficients as a matrix and possibly do Kronecker product to get time specific estimates as well. Searching for all these I ended up reading your papers and have thought it can address most of my concerns.
I may get back here and ask some questions while working on this. Thanks again for your advice and sharing.

1 Like
7

I would still like to understand what you mean by a person based estimate. Predictors and responses usually vary across individuals, times, weather patterns, coffee consumption, whatever you like – whether or not you need to allow for differences in the model parameters due to these factors is not a straightforward decision! depends on exactly what you want to learn, and how close the resulting inferences are when we represent the complex infinite dimensional complexity of the ‘true’ parameters as single points…

8

I couldn’t enter to the link.

9

unfortunately I can’t edit the original post anymore. The intro paper can be found here: https://github.com/cdriveraus/ctsem/raw/master/vignettes/hierarchicalmanual.pdf
and the simpler / more visual quick start is on https://cdriver.netlify.app/
the discussion paper that describes the linear version of the hierarchical sde in some detail is here:
https://www.researchgate.net/profile/Charles_Driver/publication/310747801_Hierarchical_Bayesian_Continuous_Time_Dynamic_Modeling/links/5c4a0a4f458515a4c73da33a/Hierarchical-Bayesian-Continuous-Time-Dynamic-Modeling.pdf

1 Like
10

Hello,
Thanks for asking this question. The model is really complex.
I have 24 variables that the fluctuation in each can impact on another.
What I like to know is
1- How a a time-variant factor (eye) at time 4 is related to eye at time 1, 2 and 3, and get population-and individual based- AR estimates and and latent values.
2- How the factor eye at time 4 is related to another time variant factor (hr) at time 1, 2, and 3 ( cross-lagged) and get population-based, individual-based estimates and latent values.
3- The interval between each time is variable for each person ranging from 1 year to 10 years. It would be great if I can get both population and individual-based estimates for time. Also, using ctsem allows me to use continuous time.
I went through the manual. In the data I am using, the factors are mainly time-dependent.
represent inputs to the system that vary over time and are independent
of fluctuations in the system.
4- the eventual goal is to know AR of each variable for the population and each individual and cross-lagged estimates for each time-point for the population and each individual and the role of time itself (estimate for time in association with each variable).
I have 800 individuals and 4 time-points measurements.

11

@Charles_Driver
I think the part that I am dealing with now is how to incorporate extended Kalman Filter in ctsem package for ordinal variables. I could not figure it out from the manual and the examples.
Thanks

12

@ssalimi yes, you need some differentiable function relating the latent processes to the expected mean and error std deviation of the indicators. the default is a linear factor model, in which case the measurement model doesn’t change depending on the latent state. I have never done much with ordinal variables, so far as I understand various approximations have been used. If I was going to start modelling such a thing I would probably start with some polynomial relation between the latent process and the prediction of y, and have the std dev of the measurement error a function of some baseline value plus the distance between the prediction of y and possible values of y. If you find something appropriate that others have done or have some specific ideas yourself, I can show how it should be included. If I find some time I’ll try it out and post an example.

13

Hi,
The suggestions have been to use Gaussian for the latent variable then logistic regression for the model. There is a case study here for ordered logistic regression in Stan.
https://betanalpha.github.io/assets/case_studies/ordinal_regression.html
I thought to change the normal_lpdf in your syntax to ordered_logistic and keep the latent model Gaussian. The distribution is Dirichlet.
But, I am not very familiar with all details of your Stan code and not sure how to modify it accordingly.

14

Here is some code from an example using binary data, showing the more appropriate but slower non integrated approach, as well as the linearised approximation.

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

invlog=function (x) exp(x)/(1 + exp(x))
n.manifest=21

#gen data
gm <- ctModel(DRIFT=-.3, DIFFUSION=.3, CINT=.1,
  TRAITVAR=diag(.3,1), #old approach to allow individual variation 
  LAMBDA= rep(1,each=n.manifest),
  n.latent=1,n.manifest=n.manifest,Tpoints=20,
  MANIFESTMEANS=c(0,rep(c(.5,-.5),each=(n.manifest-1)/2)),T0MEANS=-.3,T0VAR=.5)

d=ctGenerate(gm,n.subjects = 50,logdtsd=.2)
d[,gm$manifestNames] <- rbinom(nrow(d)*gm$n.manifest,size=1,prob=invlog(d[,gm$manifestNames]))

#model to fit
m <- ctModel( n.latent = 1,
  n.manifest = n.manifest,
  MANIFESTMEANS = c(0,paste0('m',2:n.manifest,'|param|FALSE')), #set prior to N(0,1), disable individual variation
  LAMBDA = rep(1,n.manifest),
  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" )

#plot(m)

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

#fit without integration
r <- ctStanFit( datalong = d,
  #fit=FALSE, #set this to skip fitting and just get the standata and stanmodel objects
  ctstanmodel = m,
  iter = 500,verbose=0,control=list(max_treedepth=8),nopriors=FALSE,
  chains = 2,plot=T,
  intoverstates = FALSE,
  optimize=FALSE,intoverpop=F)
s=summary(r)
s

#r$standata contains data structure
#r$stanmodeltext contains model text

#fit with integration (linear approximation)
ro <- ctStanFit( datalong = d,
  ctstanmodel = m,cores=2,
  intoverstates = T,nopriors=T,
  optimize=T,intoverpop=T)
so=summary(ro)
so

The important lines of code for you will be the measurement model (o indexes all observed data, o1 indexes indicators flagged as manifesttype 1, which is currently used for binary but you can just change it as neeed):

syprior[o] = sLAMBDA[o] * state[1:nlatent] + sMANIFESTMEANS[o,1];
syprior[o1] = to_vector(inv_logit(to_array_1d( syprior[o1] )));

and the likelihood:

if(nbinary_y[rowi] > 0) llrow[rowi] += sum(log(Y[rowi,o1d] .* (syprior[o1d]) + (1-Y[rowi,o1d]) .* (1-syprior[o1d])));

1 Like
15

@Charles_Driver
Thank you for the hints and code.