# How to set hard constraints on a transformed parameter within a hierarchical model?

Hi there, I have developed a hierarchical model of behavioural data in a psychophysical experiment. The model defines a grand average prior on parameters ‘bias’, ‘sens’, and ‘n_exp’. The model also defines per-session and per-subject variations to these parameters, which is where the hierarchical aspect of the model comes in. See code below.

So for example, the N_EXP parameter for subject 1 on session 1 is set as:
N_EXP = n_exp + <deviation for subject 1> + <deviation for session 1>
In the code below these deviations are held in the matrices b_sess and b_subj

One problem I’m having is that the variation in the parameters across sessions and subjects occasionally sets this N_EXP value to be less than zero, even though the global parameter n_exp was constraints to be positive. This is a problem when computing the likelihood because the contrastLeft or contrastRight values (which are raised to the power N_EXP) are sometimes 0. This makes the term 0^(negative number) which is undefined.

Therefore is there a way to constrain N_EXP transformed parameter like this, given that it is computed differently subject to subject and session to session?

Thank you
Peter

``````data {
int<lower=1> numTrials;
int<lower=1> numSessions;
int<lower=1> numSubjects;
int<lower=1,upper=numSessions> sessionID[numTrials];
int<lower=1,upper=numSubjects> subjectID[numTrials];
vector<lower=0,upper=1>[numTrials] contrastLeft;
vector<lower=0,upper=1>[numTrials] contrastRight;
int<lower=1,upper=3> choice[numTrials]; // 1=Left, 2=Right, 3=NoGo
int<lower=1,upper=numSubjects> subjID_session[numSessions]; //the subject ID for each session in the dataset
}
parameters {
//global parameters
vector bias;
vector<lower=0> sens;
real<lower=0> n_exp;

//per session deviations
vector<lower=0> sd_sess;
matrix[5,numSessions] z_sess; //standard normal variable used to draw from the covariance matrix
cholesky_factor_corr rho_sess; //correlations of deviations

//per subject deviations
vector<lower=0> sd_subj;
matrix[5,numSubjects] z_subj;
cholesky_factor_corr rho_subj;
}
transformed parameters {
vector logOdds[numTrials]; // ln pL/pNG, ln pR/pNG, ln pNG/pNG
matrix[5,numSessions] b_sess;
matrix[5,numSubjects] b_subj;

//draw samples of sess and subj deviations, according to the covariance structure in rho_ & sd_
b_sess = diag_pre_multiply(sd_sess, rho_sess) * z_sess;
b_subj = diag_pre_multiply(sd_subj, rho_subj) * z_subj;

{
//temp variables
real BL;
real BR;
real SL;
real SR;
real N_EXP;

//compute (non)linear model
for (n in 1:numTrials)
{
BL = bias + b_sess[1,sessionID[n]] + b_subj[1,subjectID[n]];
BR = bias + b_sess[2,sessionID[n]] + b_subj[2,subjectID[n]];
SL = sens + b_sess[3,sessionID[n]] + b_subj[3,subjectID[n]];
SR = sens + b_sess[4,sessionID[n]] + b_subj[4,subjectID[n]];
N_EXP= n_exp + b_sess[5,sessionID[n]] + b_subj[5,subjectID[n]];

logOdds[n] = BL + SL*contrastLeft[n]^N_EXP;
logOdds[n] = BR + SR*contrastRight[n]^N_EXP;
logOdds[n] = 0;
}
}
}
model {
//priors on global parameters
bias ~ normal(0, 2);
sens ~ normal(5, 2);
n_exp ~ normal(0.5,0.25);

//make z_std_normal be standard normal (non centred parameterisation)
to_vector(z_sess) ~ normal(0, 1);
to_vector(z_subj) ~ normal(0, 1);

//prior on the variation of the per-session deviations
sd_sess ~ cauchy(0,1);
sd_subj ~ cauchy(0,1);

//prior on the cholesky factor of the covariance matrix
rho_sess ~ lkj_corr_cholesky(2.0); //penalises extreme correlations between the deviations
rho_subj ~ lkj_corr_cholesky(2.0); //penalises extreme correlations between the deviations

for (n in 1:numTrials) {
choice[n] ~ categorical_logit( logOdds[n] );
}
}
generated quantities {
corr_matrix corr_sess;
corr_matrix corr_subj;
vector[numTrials] log_lik;

//write correlation matrix
corr_sess = rho_sess * rho_sess';
corr_subj = rho_subj * rho_subj';

//write loglik
for (n in 1:numTrials){
log_lik[n] = categorical_logit_lpmf(choice[n] | logOdds[n] );
}
}
``````

I realize this is an old post, but I have the same question.

Does anyone have an answer to this question?

Thank you,