Issue when trying to fit a linear mixed effects model using Rstan

Starting to get out of my expertise here, but sigma and tau still seem to have very large values in the posterior, and very large variance between chains on each. Is there a way to rescale the data so that these parameters are approximately on the ±10 scale?

1 Like

Actually started looking at the data and I think maybe I can provide some expertise in querying your model choice (rather than attempting to diagnose why your current choice isn’t working). The data seem to be structured such that you have multiple individuals, each with a unique set of predictors (age, menopause, ppfev, ppdlco), and multiple measures of the outcome FEV1 obtained at time intervals specified by timeSince_t0. You model the influence of the predictors as a simple linear effect on the mean outcome and you model the influence of time as a Brownian stochastic process. I’m not familiar with the latter; is it common in your field? Does it encode deep knowledge of the data generating process? If not, then I’d suggest maybe taking a look at the more common (to my admittedly limited exposure) approach of a Gaussian Process (indeed, I suspect that there is probably a mathematical connection between Brownian stochastic process and a Gaussian process). GPs have been decently explored by Stan folks and you may be able to leverage more expertise by using that approach than with something you coded by hand yourself.

You might take a look at this example on how to do a non-hierarchical GP, and then this example on how to allow variability in the GP from participant to participant.

Oh, and note that the high correlation between the age and menopause variables might induce inefficient sampling that you could avoid by using the QR trick.

2 Likes

I actually thought about re-scaling/centering data. But I can only center age, ppfev1 and ppdlco. I will scale that and see whether it changes anything.

Yes, you are correct. My model is a random intercept model with a Brownian motion stochastic process. Random intercept is used to describe between-subject variation and Brownian motion is used to capture the evolution of outcome variable through time within each subject.

In biostatistics we commonly use BM and related stochastic processes to capture the within-subject variations. According to a preliminary study I did using the r package “lmenssp”, integrated Ornstein Uhlenbeck (IOU); which is related to scaled BM captures, within-subject variation well. So my plan was to use a more generalized version of BM, which is fractional BM for this data.

Yes, you are correct. BM and fractional BM are both (rather can be) gaussian processes. I will try to use the links you attached for my model. Thank you! :)

2 Likes


Now my model is much much better! It takes almost 22 hrs to run the model. I know it’s not the best. But now I can work on improving what I have. :)

Thank you @Ara_Winter and @mike-lawrence for your help! I’m very new to stan (2-3 weeks) and I have little to no experience with it. So I really appreciate your help. :)

2 Likes

Hey that’s looking way better. Glad it’s up and running.

I am not the best at optimizing models for speed. Hopefully someone else can jump in. I can look around and see who tackles those kinds of questions.

1 Like

If GPUs are available (you’d want one per chain, so likely run on EC2 or GCE) there’s GPU-based parallelism that would speed up the GP part. And the new reduce_sum should pretty flexibly speed up the final likelihood part. If it’s critical to get things going faster for multiple repeated applications of this model that is. Might not be worth working it all out if this is a one-off sampling scenario. (Though if you did work it out I imagine lots of other folks would love a thorough how-to on setting up a one-chain-per-multi-core-plus-GPU-instance on one of the cloud providers for models with simultaneous use of GPU-acceleration and reduce-sum )

Double-checking: no divergences and good rhats now?

1 Like

All R hats were <= 1. :) But there were 3 divergent transitions and I’m using adapt_delta = 0.999. I got the following two warnings too.

2: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-ess 
3: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-ess 

I thought running the model for more iterations. Any advice or solutions you could think of besides that???

I think that’s a great idea. But, I don’t have GPUs available right now.

Thanks!

Does H need to be this?

or can it be H ~ normal(0.5, 0.5) ? I didn’t think that will solve your ESS error but it’s a little cleaner.

There might be a need to parameterize the model? Looks like in the pairs plot some of your alphas are correlated. I am not as familiar with these kinds of models so maybe it’s how they are.

Yes, H is in (0, 1) by definition. Which is an open interval. But in the constraints for H, I added <lower=0, upper=1>. I don’t know whether that’s causing a problem.

Re-parameterization is a good idea. I will try that.

Thanks!

1 Like