Repeated measure hierarchical reinforcement learning with groups (2 x 2 x 2 design)

[martinmodrak]

So, there is a lot to unpack… I’ll admit I don’t follow all the details, but for the sake of getting you at least some hints sooner rather than later, I’ll start with some surface-level observations and hope they’ll help you move forward. I’ll admit that I’ve never used hbayesdm or built any models like the one you describe…

@carinaufer

However the model does not fit very well resulting in Rhat > 3 for some parameters. (I read in the forum that I can ignore Rhat for predictions from the generate quantities part. So is it correct that I can ignore them? Because some show me NAs there?)

Yes, that looks bad. Rhat for generated quantities can sometimes be safely ignored - if it is NA, it often means that all posterior samples from all chains have the same value. To what extent this is to be expected and to what extent it may signal a bug is a model-specific question.

Anyway, I would take this one step at a time and first check, whether there isn’t some fundamental problem unrelated to your tweaks to the model - if you take the hbayesdm model as-is (ignoring for a while that it doesn’t capture some aspects of the data), are you able to fit without problems (divergences, treedepth warnings, reasonable Rhat)? If so, are you then able to use the same Stan code to fit the same model yourself? (or maybe that’s what you did already?) Are you able to fit some example dataset from hbayesdm with your custom code?

@carinaufer

I have been reading the manual and the forum, but both at the same time causes a lot more confusion than it clarifies questions. Often I feel like the manual is to sparse explanation in order that I can really understand it. I try to structure my questions below:

Yes, that happens, been there some time ago myself… Unfortunately once you move beyond simple stuff, there really is a lot of background knowledge one needs to build up to be able to succesfully run complex models and lot of it is not yet captured in accessible tutorials/classes/… So yes, the first steps can be overwhelming, and Bayesian modelling can be very hard, but I believe you will be able to conquer it :-) Just don’t be afraid to ask.

(I would recommend McElreath’s Statistical Rethinking as a good primer that builds a lot of knowledge from the ground up - there is even an online version of course materials and talks: GitHub - rmcelreath/statrethinking_winter2019: Statistical Rethinking course at MPI-EVA from Dec 2018 through Feb 2019, but I think it still doesn’t cover some of the more advanced stuff. The case studies and courses by Mike Betancourt Writing - betanalpha.github.io do cover some of the more advanced stuff, but can be hard to digest)

@carinaufer

Since the model does not fit well, I wanted to include an additional hierarchical layer accounting for the repeated measure.

It is very rare that you are able to get rid of fitting issues by adding more parameters/structure/… So I would first debug the simpler model. One thing to look at would be the pairs plot for the parameters (possibly subsetting to show only some of the A_pr and tau_pr and keep the plot reasonable). See also Visual MCMC diagnostics using the bayesplot package • bayesplot for some more discussion on using visualisations to debug.

@carinaufer

I am not sure if and how I should limit the parameters. I read in the forum that a truncation through initialization as vector<lower=0, upper=1>[x] does not really limit the parameter since there are two different spaces one for the estimation itself (not limited still moves through the unbounded space) and one for the parameter that is returned (which is limited than). Did I understand it correct?

If you have a constrained parameter - e.g. vector<lower=0, upper = 1>[5] x; then all elements of x will be constrained between 0 and 1. Stan creates an unconstrained version of the parameter behind the scenes, but for the most part, you don’t need to worry about it too much.

@carinaufer

3. How can I test my hypothesis adequately? I thought I just estimate the parameters and than run classical statistical analysis where my model would be something like:

lme((alpha)~S_O*Admin*Observation+blockN, random= ~S_O|subject, contrasts=list(Admin="contr.sum",Observation="contr.sum",S_O="contr.poly"),
                  corAR1(0,form=~1|subject), data = parameterList_con, method = "ML",na.action=na.exclude, control=list(opt="optim"))

However, it would not be a bayesian way and would not account for the uncertainty and take the point estimates.

Hypothesis testing is a bit of an unnatural approach in Bayesian inference (though it can be done). The most natural approach is to inteprepte either the posterior distribution of the parameters, or predictions of the model.

If you want to think about some sort of Bayesian alternative to classical point hypothesis, I’ve written about this elsewhere: Bayesian Hypothesis testing - #15 by martinmodrak