Running ctsem on cluster

[ballardtj]

Hi again Charles,

I'm running a ctsem model on a cluster using full Bayes (as opposed to optimisation). I'm wondering whether it's possible to make use of stan's within-chain parallelisation features to spread operations across cores (e.g., as implemented via the map_rect or reduce_sum functions)? Does ctsem have anything built in to make use of these features? Or are there any other ways of parallelising/speeding up models that run using full Bayes on a cluster? I see that the ctStanFit() function allows you set specify the number of cores you want to use, but presumably, this allows only one core per chain?

Cheers,
Tim

[cdriveraus]

Using the different optimizers ctsem offers allows full parallelization, but it's a pretty big rewrite to get maprect working so that the sampler can be properly parallel too. Now that rstan is finally getting up to date this is on the list for next year :)

[ballardtj]

Thanks, Charles. That's good to know! I'll let you know when I get around to attempting the GPU.

I'm having another issue when trying to generate data for a parameter recovery analysis. Here's the model:

ctmodel = ctModel(
type='stanct',
LAMBDA=diag(2),
n.latent = 2,
n.manifest = 2,
manifestNames = c("work_hours_scaled","mh"),
latentNames = c("HPW","MH"),
id = "subject",
time = "time",
DIFFUSION = matrix(c("sigma_change_hpw",0,0,"sigma_change_mh"),nrow=2),
DRIFT = matrix(c("tau_hpw","mh_on_hpw + mh_on_hpw2 * HPW","hpw_on_mh","tau_mh"),nrow=2),
CINT = c("mu_hpw","mu_mh"),
T0MEANS = c("init_hpw_mean","init_mh_mean"),
T0VAR = matrix(c("init_hpw_sd",0,0,"init_mh_sd"),nrow=2),
MANIFESTVAR = matrix(c("hpw_error_sd",0,0,"mh_error_sd"),ncol=2),
MANIFESTMEANS = c(0,0),
TRAITVAR = NULL,
T0TRAITEFFECT = NULL,
PARS = c("mh_on_hpw","mh_on_hpw2")
)

I'm looking to generate data by sampling from the prior. My understanding is that this can be done using the gendata=TRUE argument to the ctStanFit() function. I've tried this a few different ways, and am getting different errors. When I use full bayes, the analysis cuts out after the burnin phase.

ssfit <- ctStanFit(dat, ctmodel,optimize=F,nopriors=F,cores=1,gendata=T)

Output below:

`Bayesian estimation via Stan's NUTS sampler requested
Free T0VAR parameters as well as indvarying T0MEANS -- fixing T0VAR pars to diag matrix of 1e-3
Compiling model...
DIAGNOSTIC(S) FROM PARSER:
Info: integer division implicitly rounds to integer. Found int division: d * d - d / 2
Positive values rounded down, negative values rounded up or down in platform-dependent way.

Sampling...

SAMPLING FOR MODEL 'ctsem' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000127 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.27 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 1000 [ 2%] (Warmup)
Chain 1: Iteration: 40 / 1000 [ 4%] (Warmup)
Chain 1: Iteration: 60 / 1000 [ 6%] (Warmup)
Chain 1: Iteration: 80 / 1000 [ 8%] (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 1: Iteration: 120 / 1000 [ 12%] (Warmup)
Chain 1: Iteration: 140 / 1000 [ 14%] (Warmup)
Chain 1: Iteration: 160 / 1000 [ 16%] (Warmup)
Chain 1: Iteration: 180 / 1000 [ 18%] (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 1: Iteration: 220 / 1000 [ 22%] (Warmup)
Chain 1: Iteration: 240 / 1000 [ 24%] (Warmup)
Chain 1: Iteration: 260 / 1000 [ 26%] (Warmup)
Chain 1: Iteration: 280 / 1000 [ 28%] (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 1: Iteration: 320 / 1000 [ 32%] (Warmup)
Chain 1: Iteration: 340 / 1000 [ 34%] (Warmup)
Chain 1: Iteration: 360 / 1000 [ 36%] (Warmup)
Chain 1: Iteration: 380 / 1000 [ 38%] (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 1: Iteration: 420 / 1000 [ 42%] (Warmup)
Chain 1: Iteration: 440 / 1000 [ 44%] (Warmup)
Chain 1: Iteration: 460 / 1000 [ 46%] (Warmup)
Chain 1: Iteration: 480 / 1000 [ 48%] (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
[1] "Error in sampler$call_sampler(args_list[[i]]) : "
[2] " c++ exception (unknown reason)"
error occurred during calling the sampler; sampling not done
Stan model 'ctsem' does not contain samples.
Error in fit@inits[[1]] : subscript out of bounds`

When I use optimisation, the analysis runs, but it looks like it's only generating data for a single time series. In other words, there's the amount of data I would expect for one subject. I then tried adding the subjectMatrices=TRUE to the ctExtract() function, but got this error:

Computing quantities for 1000 samples...
starting worker pid=69538 on localhost:11705 at 15:04:14.063
starting worker pid=69552 on localhost:11705 at 15:04:14.261
Loading required package: data.table
Loading required package: data.table
Error in [email protected]$stan_fit_instance$constrain_pars(upars) :
c++ exception (unknown reason)
Error in [email protected]$stan_fit_instance$constrain_pars(upars) :
c++ exception (unknown reason)
Error in [email protected]$stan_fit_instance$constrain_pars(upars) :
c++ exception (unknown reason)
Error in [email protected]$stan_fit_instance$constrain_pars(upars) :
c++ exception (unknown reason)

Any help you can provide would be much appreciated!

[cdriveraus]

I may have to get rid of that argument, a nicer approach to generating data has been needed for some time now, does the ctStanGenerate function work for you perhaps?

[ballardtj]

Unfortunately not. I get a similar error:

ssfit2 <- ctStanFit(dat, ctmodel,optimize=F,nopriors=F,cores=2,iter=20)

ssfit_gen = ctStanGenerate(
cts=ssfit2,
cores=1
)

Fitting model to empty dataset...
Maximum a posteriori estimation requested
Compiling model...
DIAGNOSTIC(S) FROM PARSER:
Info: integer division implicitly rounds to integer. Found int division: d * d - d / 2
Positive values rounded down, negative values rounded up or down in platform-dependent way.

Using 1/8 logical CPU cores
Optimizing...
Finishing optimization...
Estimating Hessian
Computing posterior with 200 samples
Compilation needed -- compiling (usually ~ 1 min)
DIAGNOSTIC(S) FROM PARSER:
Info: integer division implicitly rounds to integer. Found int division: d * d - d / 2
Positive values rounded down, negative values rounded up or down in platform-dependent way.

Generating data from posterior
Using 1/8 logical CPU cores
Error in [email protected]$stan_fit_instance$constrain_pars(upars) :
c++ exception (unknown reason)
Error in [email protected]$stan_fit_instance$constrain_pars(upars) :
c++ exception (unknown reason)
Error in [email protected]$stan_fit_instance$constrain_pars(upars) :
c++ exception (unknown reason)
Error in [email protected]$stan_fit_instance$constrain_pars(upars) :
c++ exception (unknown reason)
Error in [email protected]$stan_fit_instance$constrain_pars(upars) :
c++ exception (unknown reason)
Error in [email protected]$stan_fit_instance$constrain_pars(upars) :
c++ exception (unknown reason)

...and a different error when I use a model and dataset instead of the fit object:

ssfit_gen = ctStanGenerate(
cts=ctmodel,
datastruct=dat,
cores=1
)`

Fitting model to empty dataset...
Maximum a posteriori estimation requested
Compiling model...
DIAGNOSTIC(S) FROM PARSER:
Info: integer division implicitly rounds to integer. Found int division: d * d - d / 2
Positive values rounded down, negative values rounded up or down in platform-dependent way.

Using 1/8 logical CPU cores
Optimizing...
Finishing optimization...
Estimating Hessian
Computing posterior with 200 samples
Error in makeNumericIDs(datalong, ctm$subjectIDname, ctm$timeName) :
(list) object cannot be coerced to type 'double'

A few other ideas I had:

1. I've tried compiling the model with cmdstanr, which has an option to sample from the prior (fixed_param = T). Full code below:

#Fit dummy model to generate model code
ssfit <- ctStanFit(dat, ctmodel,optimize=F,nopriors=F,cores=2,iter=2,chains=1,gendata=T,savesubjectmatrices = T)

#Compile with cmdstan
mod <- cmdstan_model(write_stan_file(ssfit$stanmodeltext))

#Sample from prior
fit <- mod$sample(
data = ssfit$standata,
iter_warmup=10,
iter_sampling=10,
chains = 2,
parallel_chains = 2,
fixed_param = T,
refresh = 1
)

#Extract draws
draws = fit$draws(format="matrix")

This works fine. And I've been able to find 'Ygen' among the variables in the draws object.

#Print unique variables in draws object
unique(gsub("\\[.*", "", colnames(draws)))

[1] "rawpopmeans" "rawpopsdbase" "sqrtpcov" "baseindparams" "rawpopsd"
[6] "rawpopcovbase" "rawpopcov" "rawpopcovchol" "rawpopcorr" "subset"
[11] "firstsub" "lastsub" "popmeans" "popsd" "popcov"
[16] "ll" "llrow" "etacova" "ycova" "etaa"
[21] "ya" "Ygen" "etaupdbasestates" "subj_T0MEANS" "subj_DRIFT"
[26] "subj_CINT" "subj_asymCINT" "pop_PARS" "pop_T0MEANS" "pop_LAMBDA"
[31] "pop_DRIFT" "pop_DIFFUSION" "pop_MANIFESTVAR" "pop_MANIFESTMEANS" "pop_CINT"
[36] "pop_T0VAR" "pop_DIFFUSIONcov" "pop_MANIFESTcov" "pop_T0cov" "pop_asymCINT"
[41] "pop_asymDIFFUSIONcov"

Based on the dimensions of the Ygen variable, there is 1 row for each datapoint in the dataframe and 1 column for each latent variable (though of course each dimension is stored as its own column with entries for each sample in the draws object).

It looks like some of the Ygen samples are NA:

draws[,str_detect(colnames(draws),"Ygen")]

A draws_matrix: 10 iterations, 2 chains, and 1816 variables
variable
draw Ygen[1,1] Ygen[2,1] Ygen[3,1] Ygen[4,1] Ygen[5,1] Ygen[6,1] Ygen[7,1] Ygen[8,1]
1 NaN NaN NaN NaN NaN NaN NaN NaN
2 NaN NaN NaN NaN NaN NaN NaN NaN
3 4.2 15 10.6 4.0 29.6 1.9 4.9 -24.1
4 NaN NaN NaN NaN NaN NaN NaN NaN
5 NaN NaN NaN NaN NaN NaN NaN NaN
6 8.2 -18 2.9 -12.2 -15.2 25.7 46.2 -4.6
7 NaN NaN NaN NaN NaN NaN NaN NaN
8 10.1 18 -7.6 1.5 -8.3 6.4 -7.2 -8.3
9 NaN NaN NaN NaN NaN NaN NaN NaN
10 NaN NaN NaN NaN NaN NaN NaN NaN
... with 10 more draws, and 1808 more variables

I think this happens because of the non-linear term in the model, which makes the system unstable in certain cases. When a run a linear model there are no NaNs.

2. Other idea I had -- can the ctKalman function be hacked to manually change the parameters that are used to generate the posterior predictives somehow? This might be too tricky....

[ballardtj]

Actually, as I'm playing around with the linear model a bit more, the ctStanGenerate() function works just fine. Is this something to do with the EKF in general or is it just the fact that the quadratic function in the non-linear model makes things too unstable?

Also, can I just confirm -- does the ctStanGenerate() generate data by sampling from the posterior or the prior?

Sorry for all the posts!

[cdriveraus]

I see. I think it's a combination of the EKF mathematics, various implementation choices I've made, and numeric issues in the Stan functions such as positive definite tolerances and the like. You occasionally see this kind of thing when fitting complex models, but often fitting to data constrains the parameters sufficiently that they never wander into the numerically difficult regions of parameter space. When you have a broad prior on a complex model (or just one prone to numerical problems), all the regions get 'explored' so to speak and the issues are more likely to arise. However, failing on one of the data generation samples shouldn't lead to complete failure, so I also worry that something else is going on. It might take me a month or two to look at this though -- in case you want to dive into it, try just getting one sample at a time and see if it at least works sometimes. An alternative approach is also to generate the data manually using an Euler-Maruyama (i.e. lots of tiny linear steps in time) approach -- can be a nice learning experience :) :)

edit: ok read the earlier post properly this time, I see that it's only failing on some samples. So you could experiment with tighter priors (probably on variances and drift terms) to see where things are breaking, or add some code to the ctStanGenerate function so that it outputs the sampled parameters as well, so you can see which parameter regions result in NA's.