Bayesian Logistic Regression with MAP prior

Modeling
1

I am very new to Stan. My Stan code is below.
I am getting this error. Can you please help? Probably it is trivial and I am messing up with the syntax and looping.

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

Unknown variable in assignment; lhs variable=lp
error in ‘model33103996a7f_rstanforum’ at line 92, column 5

90: EFF[c] ~ binomial(N[c], q[c]);
91: 
92: lp = multi_normal_lpdf(log_ab[1] | mn1,cov1)+multi_normal_lpdf(log_ab[2] | mn2,cov2);
        ^
93:

PARSER EXPECTED: 

data {
  int<lower=1> C;              // number of cohorts
int<lower=1> N[C];           //number of Subjects by cohort
int<lower=0> EFF[C];         //number of EFFs by cohort
real<lower=0> dd[C,2];     //dd amount by cohort
real<lower=0> dd_ref[2];   //standard dd  
int<lower=1> D;              //number of all combinations
real<lower=0> dd_set[D,2]; //dd amount by combination

real Mean_Ma1;        //normal parameters for mu1 mean
real<lower=0> SD_Ma1; //normal parameters for mu1  sd
real Mean_Ma2;       //normal parameters for mu2 mean
real<lower=0> SD_Ma2; //normal parameters for mu2  sd
real Mean_Sa1;        // t distn for tau1 mean df is 4 default
real<lower=0> SD_Sa1; // t distn for tau1 sd
real Mean_Sa2;        // t distn for tau2 mean 
real<lower=0> SD_Sa2; // t distn for tau2 sd 

real Mean_Mb1;        //normal parameters for mu1 mean
real<lower=0> SD_Mb1; //normal parameters for mu1  sd
real Mean_Mb2;       //normal parameters for mu2 mean
real<lower=0> SD_Mb2; //normal parameters for mu2  sd
real Mean_Sb1;        // t distn for tau1 mean df is 4 default
real<lower=0> SD_Sb1; // t distn for tau1 sd
real Mean_Sb2;        // t distn for tau2 mean 
real<lower=0> SD_Sb2; // t distn for tau2 sd 
real<lower=-1, upper=1> Rho_min;
real<lower=-1, upper=1> Rho_max;
int<lower=0, upper=1> Eta_Dist;  //0:Normal, 1:LogNormal
real e; // uniform 
real f;
real Eta_Mean;
real<lower=0> Eta_SD;
}

parameters {
vector[2] log_ab[2]; // 4
real eta0;  //1

}

transformed parameters {
real eta;
real<lower=Rho_min, upper=Rho_max> rho1;
real<lower=Rho_min, upper=Rho_max> rho2;
real ma1; 
real ma2;
real<lower=0> sa1;
real<lower=0> sa2;
real mb1; 
real mb2;
real<lower=0> sb1;
real<lower=0> sb2;
real log_a[2];
real log_b[2];
real b[2];
real q0[C,2];
real q[C];
vector[2] mn1;
matrix[2,2] cov1;
vector[2] mn2;
matrix[2,2] cov2;

if (Eta_Dist == 0) {
eta = eta0; // Normal
} else {
eta = exp(eta0); // lognormal
}
for (k in 1:2) {
log_a[k] = log_ab[k,1];
log_b[k] = log_ab[k,2];
b[k] = exp(log_b[k]);
}
for (c in 1:C) {
for (k in 1:2)
q0[c,k] = dd[c,k] == 0 ? 0 : inv_logit(log_a[k] + b[k]*log(dd[c,k]/dd_ref[k]));
q[c] = inv_logit(logit(q0[c,1] + q0[c,2] - q0[c,1]*q0[c,2]) + eta*(dd[c,1]/dd_ref[1])*(dd[c,2]/dd_ref[2]));
mn1[1]= ma1; mn1[2]=ma2; 
cov1[1,1]= square(sa1); cov1[1,2]= sa1*sa2*rho1;
cov1[2,1]= sa1*sa2*rho1; cov1[2,2]= square(sa2);
mn2[1]= mb1; mn2[2]=mb2; 
cov2[1,1]= square(sb1); cov2[1,2]= sb1*sb2*rho2;
cov2[2,1]= sb1*sb2*rho2; cov2[2,2]= square(sb2);
}
}

model {

for (c in 1:C)
vector[C] lp;
real final;
lp[c]= binomial_lpmf(EFF[c]| N[c], q[c])+multi_normal_lpdf(log_ab[1] | mn1,cov1)+multi_normal_lpdf(log_ab[2] | mn2,cov2);

final= exp(lp);
ma1 ~ normal(Mean_Ma1,SD_Ma1);
ma2 ~ normal(Mean_Ma2,SD_Ma2);
sa1 ~ student_t(4, Mean_Sa1, SD_Sa1);
sa2 ~ student_t(4, Mean_Sa2, SD_Sa2);
mb1 ~ normal(Mean_Mb1,SD_Mb1);
mb2 ~ normal(Mean_Mb2,SD_Mb2);
sb1 ~ student_t(4, Mean_Sb1, SD_Sb1);
sb2 ~ student_t(4, Mean_Sb2, SD_Sb2);
rho1~ uniform(e,f)
rho2 ~ uniform(e,f)
log_ab[1] ~ multi_normal(mn1, cov1);
log_ab[2] ~ multi_normal(mn2,cov2);

eta0 ~ normal(Eta_Mean, Eta_SD);
}


generated quantities {
real pE[D,2];
real p[D];
for (d in 1:D)
for (k in 1:2)
pE[d,k] = dd_set[d,k] == 0 ? 0 : inv_logit(log_a[k] + b[k]*log(dd_set[d,k]/dd_ref[k]));
for (d in 1:D)
p[d] = inv_logit(logit(pE[d,1] + pE[d,2] - pE[d,1]*pE[d,2]) + eta*(dd_set[d,1]/dd_ref[1])*(dd_set[d,2]/dd_ref[2]));
}
2

lp is not declared I think. Puting vector[2] lp; at the start of the model block might solve the problem.

3

Yes. trying that.

4

There should be more error message here.

All that numbering of variables is an opportunity to use arrays instead. All of Stan’s distributions are vectorized so this will be much shorter that way.

Where are rho1 and rho2 defined? I only see declarations and uses of it, which won’t work. Variables in Stan need to be defined before they are used. I think rho1 and rho2 should be parameters, since you have sampling statements, but not definitions. Sampling statements only increment the log density. They don’t literally do random number generation.

Also, bounds shouldn’t be used for rejection if that was the intent.

You might want to start with a simpler model and work up to something like this after getting a feel for how Stan works.

5

Thanks for the reply.
Going through some videos of RStan and implementing it.