using hierarchical models with random slope resutling in sampling problem

[YuanboBQ]

Hi, I am currently working with the HSSM package (v 0.2.3) and have encountered an issue while handling my own data, including the "cavanagh_theta" data provided by the package. Specifically, when using hierarchical models with random slope (e.g., stim|participant_id), I am facing the following problem: the sampling process terminates quickly, but the obtained samples are mostly invalid, resulting in all NaN values for r_hat. However, if remove the random slope, the sampling process will be more normal and the results of the sampling will also be acceptable. Could you please investigate this issue? Can you check this issue? Currently, there are still some problems with HSSM when dealing with hierarchical models, so it is still difficult to use it in formal research.

hierarchical_ddm_model_reg_vazt = hssm.HSSM(
        model = "ddm",
        data = cav_data,
        prior_settings="safe",
        loglik_kind = "approx_differentiable", #"analytical"      
        include = [
            {
                "name": "v",
                "formula": "v ~ 0 + stim + (0 + stim|participant_id)",      
                "link": "identity",
            },
            {
                "name": "a",
                "formula": "a ~ 0 + stim + (0 + stim|participant_id)",      
                "link": "identity",
            },
            {
                "name": "z",
                "formula": "z ~ 0 + stim + (0 + stim|participant_id)",      
                "link": "identity",
            },      
            {
                "name": "t",
                "formula": "t ~ 0 + stim + (0 + stim|participant_id)",      
                "link": "identity",
            },                                                 
                     
        ],
    )
    
    print(hierarchical_ddm_model_reg_vazt)

or the following code

# Define a basic hierarchical model with trial-level covariates
hierarchical_ddm_model_reg_vazt = hssm.HSSM(
    model = "ddm",
    data = cav_data,
    loglik_kind = "approx_differentiable", #"analytical"      
    include = [
        {
            "name": "v",
            "formula": "v ~ 1 + stim + (1 + stim|participant_id)",
            "link": "identity",
        },               
        {
            "name": "a",
            "formula": "a ~ 1 + stim + (1 + stim|participant_id)",
            #"formula": "v ~ (1|subj_idx) + theta",
            "link": "identity",
        },  
        {
            "name": "z",
            "formula": "z ~ 1 + stim + (1 + stim|participant_id)",
            "link": "identity",
        }, 
        {
            "name": "t",
            "formula": "t ~ 1 + (1|participant_id)",
            #"formula": "v ~ (1|subj_idx) + theta",
            "link": "identity",
        },                                         
    ],
)

print(hierarchical_ddm_model_reg_vazt)

hierarchical_ddm_model_reg_vazt.sample(
    sampler="nuts_numpyro", 
    chains=2,  # how many chains to run
    cores=4,  # how many cores to use
    draws=1000,  # number of draws from the markov chain
    tune=500,  # number of burn-in samples
    #target_accept = 0.90,
    )
    

[frankmj]

Hi Yuanbo, try this (for the angle model but should also work for ddm if you just remove theta) data= hssm.load_data('cavanagh_theta'); model_angle_cav = hssm.HSSM( data= hssm.load_data('cavanagh_theta'), model = "angle", prior_settings = "safe", loglik_kind="approx_differentiable", include=[ { "name": "v", "formula": "v ~ 1 + stim + (1|participant_id)", "prior": { "Intercept": {"name": "Normal", "mu": 1, "sigma": 2, "initval": 1}, "stim": {"name": "Normal", "mu": 0, "sigma": 1},#, "initval": 0.1}, "1|participant_id": {"name": "Normal", "mu": 0, # using non-centered approach so mu's of indiv subject offsets should be 0 "sigma": {"name": "HalfNormal", "sigma": 1, }, "initval": 0 }, }, }, { "name": "t", "formula": "t ~ 1 + (1| participant_id)", "prior": { "Intercept": {"name": "Gamma", "mu": 0.4, "sigma": 0.2, "initval": 0.3}, "1|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": .3, "initval": .1 }, }, }, "link": "identity", }, { "name": "z", "formula": "z ~ 1 + (1 | participant_id)", "prior": { "Intercept": {"name": "Beta", "alpha": 10, "beta": 10, "initval": 0.5}, "1|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": 0.05, "initval": .01 }, }, }, }, { "name": "theta", "formula": "theta ~ 1 + (1 | participant_id)", "prior": { "Intercept": {"name": "HalfNormal", "sigma": 0.4, "initval": 0.2}, "1|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": .2, "initval": .05 }, }, }, "link": "identity", }, { "name": "a", #"formula": "a ~ 1 + stim + theta+ dbs + theta:dbs + (1 + theta +theta:dbs| participant_id)", "formula": "a ~ 1 + stim + (1 | participant_id)", "prior": { "Intercept": {"name": "Gamma", "mu": 0.5, "sigma": 1.75, "initval": 1}, "stim": {"name": "Normal", "mu": 0, "sigma": 1},#, "initval": 0.1}, "theta": {"name": "Normal", "mu": 0, "sigma": 1},#, "initval": 0.1}, "dbs": {"name": "Normal", "mu": 0, "sigma": 1},#, "initval": 0.1}, "theta:dbs": {"name": "Normal", "mu": 0, "sigma": 1},#, "initval": 0.1}, "1|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": 1, "initval": .3 }, }, "theta|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": .2, "initval": .05 }, }, }, }, ], noncentered = True, ) model_angle_cav Michael J Frank, PhD | Edgar L. Marston Professor Director, Carney Center for Computational Brain Science <https://www.brown.edu/carney/ccbs> Laboratory of Neural Computation and Cognition <https://www.lnccbrown.com/> Brown University website <http://ski.clps.brown.edu>

…

On Tue, Jul 16, 2024 at 11:07 PM yuanbo ***@***.***> wrote: Hi, I am currently working with the HSSM package (v 0.2.3) and have encountered an issue while handling my own data, including the "cavanagh_theta" data provided by the package. Specifically, when using hierarchical models with random slope (e.g., stim|participant_id), I am facing the following problem: the sampling process terminates quickly, but the obtained samples are mostly invalid, resulting in all NaN values for r_hat. However, if remove the random slope, the sampling process will be more normal and the results of the sampling will also be acceptable. Could you please investigate this issue? Can you check this issue? Currently, there are still some problems with HSSM when dealing with hierarchical models, so it is still difficult to use it in formal research. hierarchical_ddm_model_reg_vazt = hssm.HSSM( model = "ddm", data = cav_data, prior_settings="safe", loglik_kind = "approx_differentiable", #"analytical" include = [ { "name": "v", "formula": "v ~ 0 + stim + (0 + stim|participant_id)", "link": "identity", }, { "name": "a", "formula": "a ~ 0 + stim + (0 + stim|participant_id)", "link": "identity", }, { "name": "z", "formula": "z ~ 0 + stim + (0 + stim|participant_id)", "link": "identity", }, { "name": "t", "formula": "t ~ 0 + stim + (0 + stim|participant_id)", "link": "identity", }, ], ) print(hierarchical_ddm_model_reg_vazt) or the following code # Define a basic hierarchical model with trial-level covariates hierarchical_ddm_model_reg_vazt = hssm.HSSM( model = "ddm", data = cav_data, loglik_kind = "approx_differentiable", #"analytical" include = [ { "name": "v", "formula": "v ~ 1 + stim + (1 + stim|participant_id)", "link": "identity", }, { "name": "a", "formula": "a ~ 1 + stim + (1 + stim|participant_id)", #"formula": "v ~ (1|subj_idx) + theta", "link": "identity", }, { "name": "z", "formula": "z ~ 1 + stim + (1 + stim|participant_id)", "link": "identity", }, { "name": "t", "formula": "t ~ 1 + (1|participant_id)", #"formula": "v ~ (1|subj_idx) + theta", "link": "identity", }, ], ) print(hierarchical_ddm_model_reg_vazt) hierarchical_ddm_model_reg_vazt.sample( sampler="nuts_numpyro", chains=2, # how many chains to run cores=4, # how many cores to use draws=1000, # number of draws from the markov chain tune=500, # number of burn-in samples #target_accept = 0.90, ) — Reply to this email directly, view it on GitHub <#494>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAG7TFCG4OAANTCMDX5RDWLZMXUX5AVCNFSM6AAAAABK7ZHPHWVHI2DSMVQWIX3LMV43ASLTON2WKOZSGQYTENJZGEYTKNY> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[YuanboBQ]

Hi, Michael
I have try your provided code, and still get the same issue that the sampling process terminates quickly and resulting in all NaN values for r_hat. Here is my revised code

hierarchical_ddm_model_reg_vazt1 = hssm.HSSM(
    model = "ddm",
    data = cav_data,
    prior_settings="safe",
    loglik_kind = "approx_differentiable", #默认为"analytical"      
    include = [
        {
            "name": "v",
            "formula": "v ~ 1 + stim + (1 + stim|participant_id)",
            "prior": {
                "Intercept": {"name": "Normal", "mu": 1, "sigma": 2, "initval": 1},
                "stim": {"name": "Normal", "mu": 0, "sigma": 1}, # "initval": 0
                #"stim|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": 0.5}, "initval": 0.5},
                "1|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": 1}, "initval": 0}
                },            
            #"link": "identity",
        },    
        {
            "name": "a",
            "formula": "a ~ 1 + stim + (1 + stim|participant_id)",
            "prior": {
                "Intercept": {"name": "Gamma", "mu": 0.5, "sigma": 1.75, "initval": 1},
                "stim": {"name": "Normal", "mu": 0, "sigma": 1},# "initval": 0.1
                #"stim|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": 0.5}, "initval": 0.5},
                "1|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": 1}, "initval": 0.3}
                },            
            #"link": "identity",
        },                       
                 
    ],
   noncentered = True,   
)

print(hierarchical_ddm_model_reg_vazt1)

hierarchical_ddm_model_reg_vazt1.sample(
                              sampler="nuts_numpyro", 
                              chains=2,  # how many chains to run
                              cores=4,  # how many cores to use
                              draws=1000,  # number of draws from the markov chain
                              tune=500,  # number of burn-in samples
                              target_accept = 0.90,
                              )

[frankmj]

Ah yes there is an issue with categorical variables used as random effects, which leads to this error. If you recode your categorical variable as numeric it should work (e.g is stim is recoded as -1 and 1 or 1,2,3) but only worthwhile if you expect a linear effect of the different levels of your categorical variable. @paul Xu ***@***.***> this is actually an issue we have to circle back on

…

On Wed, Jul 17, 2024 at 8:04 AM yuanbo ***@***.***> wrote: Hi, Michael I have try your provided code, and still get the same issue that the sampling process terminates quickly and resulting in all NaN values for r_hat. Here is my revised code hierarchical_ddm_model_reg_vazt1 = hssm.HSSM( model = "ddm", data = cav_data, prior_settings="safe", loglik_kind = "approx_differentiable", #默认为"analytical" include = [ { "name": "v", "formula": "v ~ 1 + stim + (1 + stim|participant_id)", "prior": { "Intercept": {"name": "Normal", "mu": 1, "sigma": 2, "initval": 1}, "stim": {"name": "Normal", "mu": 0, "sigma": 1}, # "initval": 0 #"stim|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": 0.5}, "initval": 0.5}, "1|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": 1}, "initval": 0} }, #"link": "identity", }, { "name": "a", "formula": "a ~ 1 + stim + (1 + stim|participant_id)", "prior": { "Intercept": {"name": "Gamma", "mu": 0.5, "sigma": 1.75, "initval": 1}, "stim": {"name": "Normal", "mu": 0, "sigma": 1},# "initval": 0.1 #"stim|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": 0.5}, "initval": 0.5}, "1|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": 1}, "initval": 0.3} }, #"link": "identity", }, ], noncentered = True, ) print(hierarchical_ddm_model_reg_vazt1) hierarchical_ddm_model_reg_vazt1.sample( sampler="nuts_numpyro", chains=2, # how many chains to run cores=4, # how many cores to use draws=1000, # number of draws from the markov chain tune=500, # number of burn-in samples target_accept = 0.90, ) — Reply to this email directly, view it on GitHub <#494 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAG7TFBFYVDID36UHSSG4X3ZMZTWPAVCNFSM6AAAAABK7ZHPHWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEMZTGI3TKOBQGM> . You are receiving this because you commented.Message ID: ***@***.***>

[digicosmos86]

Hi @YuanboBQ,

Can you try adding categorical="stim" to your model definition? Like this:

hierarchical_ddm_model_reg_vazt = hssm.HSSM(
        model = "ddm",
        data = cav_data,
        prior_settings="safe",
        loglik_kind = "approx_differentiable", #"analytical"      
        categorical="stim" # a str or list of all categorical variables
        include = [
            {
                "name": "v",
                "formula": "v ~ 0 + stim + (0 + stim|participant_id)",      
                "link": "identity",
            },
            {
                "name": "a",
                "formula": "a ~ 0 + stim + (0 + stim|participant_id)",      
                "link": "identity",
            },
            {
                "name": "z",
                "formula": "z ~ 0 + stim + (0 + stim|participant_id)",      
                "link": "identity",
            },      
            {
                "name": "t",
                "formula": "t ~ 0 + stim + (0 + stim|participant_id)",      
                "link": "identity",
            },                                                 
        ],
    )

[AlexanderFengler]

Hi @YuanboBQ ,
If you use a regression formulation as v ~ 1 + (1 + C(stim)|participant_id) do you run into the same issue?

[YuanboBQ]

Hi, @frankmj and @digicosmos86, I have tried to add categorical="stim" or use the regression formulation as v ~ 1 + (1 + C(stim)|participant_id) or combine both of the two, and still got the same issue.

[AlexanderFengler]

@YuanboBQ ,
Could you take the regression out of the t parameter?
t is providing a bit of a pain point thus far, might be the actual culprit here too :/.

[YuanboBQ]

Hi, @AlexanderFengler, I did not include the regression of the t parameter. Here is my code

hierarchical_ddm_model_reg_vazt = hssm.HSSM(
    model = "ddm",
    data = cav_data,
    prior_settings="safe",
    loglik_kind = "approx_differentiable", #默认为"analytical"  
    #categorical="stim",    
    include = [
        {
            "name": "v",
            "formula": "v ~ 1 + stim + (1 + C(stim)|participant_id)",
            "prior": {
                "Intercept": {"name": "Normal", "mu": 1, "sigma": 2, "initval": 1},
                "stim": {"name": "Normal", "mu": 0, "sigma": 1}, # "initval": 0
                #"stim|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": 0.5}, "initval": 0.5},
                "1|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": 1}, "initval": 0}
                },            
            "link": "identity",
        },    
        {
            "name": "a",
            "formula": "a ~ 1 + stim + (1 + C(stim)|participant_id)",
            "prior": {
                "Intercept": {"name": "Gamma", "mu": 0.5, "sigma": 1.75, "initval": 1},
                "stim": {"name": "Normal", "mu": 0, "sigma": 1},# "initval": 0.1
                #"stim|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": 0.5}, "initval": 0.5},
                "1|participant_id": {"name": "Normal", "mu": 0, "sigma": {"name": "HalfNormal", "sigma": 1}, "initval": 0.3}
                },            
            "link": "identity",
        },                       
                 
    ],
   #noncentered = True,   
)

print(hierarchical_ddm_model_reg_vazt)