thanks, you are right. I change the declaration for g and mn as vector[] , the message gone. Now I also want to sample
the row vector of the matrix U[200,3] from multi_normal distribution, I have tried several version. this is the last one,
but not work..Irene
protein_code <-'
data {
int ns; //number of records
int nobs; //number of observation
int nmiss; //number of missing
int ncensor; //number of censored
int<lower=0> subject[ns]; //subject id
int<lower=0> proteinid[ns]; //protein id
int<lower=0,upper=1> classno[ns]; //normal comparing with control
int<lower=0,upper=1> run1[ns]; //run 1
int<lower=0,upper=1> run2[ns];
int<lower=0,upper=1> run3[ns];
int<lower=0,upper=1> run4[ns];
int<lower=0,upper=1> run5[ns];
int<lower=0,upper=1> run6[ns];
int<lower=0,upper=1> run7[ns];
int<lower=0,upper=1> run8[ns];
int<lower=0,upper=1> tlab113[ns]; //label 113
int<lower=0,upper=1> tlab114[ns];
int<lower=0,upper=1> tlab115[ns];
int<lower=0,upper=1> tlab116[ns];
int<lower=0,upper=1> tlab117[ns];
int<lower=0,upper=1> tlab118[ns];
int<lower=0,upper=1> tlab119[ns];
int<lower=0,upper=1> tlab121[ns];
real logofAUC_obs[nobs];
real logofAUC_censor[ncensor];
real m_z_centered[ns]; //logofAUC
corr_matrix[3] R;
vector[3] Sigma_sd;
cov_matrix[3] invprec;
vector[3] mn;
int<lower=0,upper=1> miss[ns];
}
transformed data
{
matrix[3,3] T;
T<-diag_matrix(Sigma_sd) * R * diag_matrix(Sigma_sd);
}
parameters {
matrix[200,3] U;
vector[3] g;
real<lower=0> ita;
real beta2[64];
real beta3_1;
real beta3_2;
real beta3_3;
real beta3_4;
real beta3_5;
real beta3_6;
real beta3_7;
real beta3_8;
real beta4_113;
real beta4_114;
real beta4_115;
real beta4_116;
real beta4_117;
real beta4_118;
real beta4_119;
real beta4_121;
real logofAUC_miss[nmiss];
real alpha;
real alpha1;
real alpha2;
}
model
{
real mu[ns];
for (pep in 1:ns)
mu[pep]<-beta2[subject[pep]]+U[proteinid[pep],1]+U[proteinid[pep],2]*m_z_centered[pep]
+U[proteinid[pep],3]*classno[pep]+beta3_1*run1[pep]+beta3_2*run2[pep]+beta3_3*run3[pep]+beta3_4*run4[pep]+beta3_5*run5[pep]+
beta3_6*run6[pep]+beta3_7*run7[pep]+beta4_113*tlab113[pep]+beta4_114*tlab114[pep]+beta4_115*tlab115[pep]+beta4_116*tlab116[pep]+beta4_117*tlab117[pep]+beta4_118*tlab118[pep]+beta4_119*tlab119[pep]+beta4_121*tlab121[pep];
beta3_1~normal(0,10); //prior for fixed coefficient
beta3_2~normal(0,10);
beta3_3~normal(0,10);
beta3_4~normal(0,10);
beta3_5~normal(0,10);
beta3_6~normal(0,10);
beta3_7~normal(0,10);
beta3_8~normal(0,10);
beta4_113~normal(0,10);
beta4_114~normal(0,10);
beta4_115~normal(0,10);
beta4_116~normal(0,10);
beta4_117~normal(0,10);
beta4_118~normal(0,10);
beta4_119~normal(0,10);
beta4_121~normal(0,10);
ita~gamma(200,10);
#hyper prior
#for (i in 1:3)
#g[i]~normal(mn[i],invprec[i,i]);
g~multi_normal(mn,invprec);
#for (prot in 1:200)
#{U[prot,1]~normal(g[1],T[1,1]);
#U[prot,2]~normal(g[2],T[2,2]);
#U[prot,3]~normal(g[3],T[3,3]);
#}
for (prot in 1:200)
U[prot,1:3]~multi_normal(g,T);
for (sub in 1:64)
beta2[sub]~normal(0,1);
for (pep in 1:nobs) //nobs includes number of censored
{logofAUC_obs[pep]~normal(mu[pep],ita); //lp__ <- lp__ + ncensor*log1m(normal_p(-9.2,mu[pep],ita))
}
for (pep in 1:ncensor)
logofAUC_censor[pep]~normal(mu[pep],ita);
for (pep in 1:nobs)
miss[pep]~bernoulli(inv_logit(alpha + alpha1*m_z_centered[pep]+alpha2*logofAUC_obs[pep]));
for (pep in 1:nmiss)
{logofAUC_miss[pep]~normal(mu[pep],ita);
miss[pep]~bernoulli_logit(alpha + alpha1*m_z_centered[pep]+alpha2*logofAUC_miss[pep]);
}
}
'