Hi stan users, I am new to stan and am trying to fit the following model.
\left(\begin{array}{c}
\hat{\theta}_{i} \\
\hat{\gamma}_{i}
\end{array}\right) \mid\left(\begin{array}{l}
\theta_{i} \\
\gamma_{i}
\end{array}\right) \sim N\left(\left(\begin{array}{l}
\theta_{i} \\
\gamma_{i}
\end{array}\right),\left(\begin{array}{cc}
\sigma_{i}^{2} & \rho_{i} \sigma_{i} \delta_{i} \\
\rho_{i} \sigma_{i} \delta_{i} & \delta_{i}^{2}
\end{array}\right)\right)
\theta_i \mid \gamma_i,a,b,\tau \sim N(a+b\gamma_i,\tau^2)
Joint prior on (a,b,\tau^{2}) : \pi_{a,b,\tau}(.)
\alpha,\beta,\tau^2 are the parameters of interest. My goal is to obtain posterior samples for these parameters.
Following is the data associated with my model:
Y <- matrix(c(-.083,0,-.416,-.342,-.211,.122,-.342,-.128,-.386,-.446,-.896,-.357,
-1.08,-1.66,-.994,-.223,-.342,-.386,-.386,-.315,0.049,0,-.342,-.844,
-.163,.077,-1.14,-.528,-.799,-.236,-1.17,-1.14,-1.51,-1.02,.03,.285,
-.868,-.363,-.799,-.308,-.139,.207,-.357,.039,-.994,-.446,.058,.255,-.03,.049),ncol=2,byrow = T)
Y1 <- cbind(Y[,2],Y[,1])
colnames(Y1) <- c("OS","PFS")
sigi <- c(.24,.263,.811,.261,.301,.376,.816,.807,.189,.186,.116,.342,.22,.141,.115,
.366,.354,.226,.165,.160,.425,.431,.397,.112,.118)
deltai <- c(.150,.147,.687,.213,.253,.269,.484,.479,.369,.178,.112,.27,.196,.156,.121,.304,.269,.188,
.164,.161,.408,.25,.438,.23,.181)
rho=0.05
TT <- array(dim=c(2,2,length(deltai)))
for(i in 1:length(deltai)){
cov <- sigi[i]*deltai[i] * rho
TT[,,i] <- matrix(c(sigi[i]^2, cov,cov, deltai[i]^2 ),ncol=2,byrow = T)
}
Here is my stan model
data {
int<lower=0> J; // number of studies
vector[2] y1; // estimated treatment effects (OS,PFS)
vector[2] y2;
vector[2] y3;
vector[2] y4;
vector[2] y5;
vector[2] y6;
vector[2] y7;
vector[2] y8; vector[2] y9; vector[2] y10;
vector[2] y11; vector[2] y12; vector[2] y13; vector[2] y14; vector[2] y15;
vector[2] y16; vector[2] y17; vector[2] y18; vector[2] y19; vector[2] y20;
vector[2] y21; vector[2] y22; vector[2] y23; vector[2] y24; vector[2] y25;
matrix[2,2] T1; // var-cov of effect estimates
matrix[2,2] T2;
matrix[2,2] T3;
matrix[2,2] T4;
matrix[2,2] T5;
matrix[2,2] T6; matrix[2,2] T7; matrix[2,2] T8; matrix[2,2] T9; matrix[2,2] T10;
matrix[2,2] T11; matrix[2,2] T12; matrix[2,2] T13; matrix[2,2] T14; matrix[2,2] T15;
matrix[2,2] T16; matrix[2,2] T17; matrix[2,2] T18; matrix[2,2] T19; matrix[2,2] T20;
matrix[2,2] T21; matrix[2,2] T22; matrix[2,2] T23; matrix[2,2] T24; matrix[2,2] T25;
}
parameters {
real alpha;
real beta;
real<lower=0> tau2;
vector[2] theta1;
vector[2] theta2;
vector[2] theta3;
vector[2] theta4;
vector[2] theta5;
vector[2] theta6; vector[2] theta7; vector[2] theta8; vector[2] theta9; vector[2] theta10;
vector[2] theta11; vector[2] theta12; vector[2] theta13; vector[2] theta14; vector[2] theta15;
vector[2] theta16; vector[2] theta17; vector[2] theta18; vector[2] theta19; vector[2] theta20;
vector[2] theta21; vector[2] theta22; vector[2] theta23; vector[2] theta24; vector[2] theta25;
}
model {
tau2 ~ normal( 0 , 1) T[0, ]; //inv_gamma(3,1);
alpha ~ normal( 0 , 1000);
beta ~ normal( 0 , 1000);
theta1[2] ~ normal( 0 , 1000); theta2[2] ~ normal( 0 , 1000); theta3[2] ~ normal( 0 , 1000); theta4[2] ~ normal( 0 , 1000); theta5[2] ~ normal( 0 , 1000);
theta6[2] ~ normal( 0 , 1000); theta7[2] ~ normal( 0 , 1000); theta8[2] ~ normal( 0 , 1000); theta9[2] ~ normal( 0 , 1000); theta10[2] ~ normal( 0 , 1000);
theta11[2] ~ normal( 0 , 1000); theta12[2] ~ normal( 0 , 1000); theta13[2] ~ normal( 0 , 1000); theta14[2] ~ normal( 0 , 1000); theta15[2] ~ normal( 0 , 1000);
theta16[2] ~ normal( 0 , 1000); theta17[2] ~ normal( 0 , 1000); theta18[2] ~ normal( 0 , 1000); theta19[2] ~ normal( 0 , 1000); theta20[2] ~ normal( 0 , 1000);
theta21[2] ~ normal( 0 , 1000); theta22[2] ~ normal( 0 , 1000); theta23[2] ~ normal( 0 , 1000); theta24[2] ~ normal( 0 , 1000); theta25[2] ~ normal( 0 , 1000);
theta1[1] ~ normal( alpha + beta * theta1[2] , tau2); // epsilon i ~ N(0,tau2)
y1 ~ multi_normal(theta1, T1);
theta2[1] ~ normal( alpha + beta * theta2[2] , tau2);
y2 ~ multi_normal(theta2, T2);
theta3[1] ~ normal( alpha + beta * theta3[2] , tau2);
y3 ~ multi_normal(theta3, T3);
theta4[1] ~ normal( alpha + beta * theta4[2] , tau2);
y4 ~ multi_normal(theta4, T4);
theta5[1] ~ normal( alpha + beta * theta5[2] , tau2);
y5 ~ multi_normal(theta5, T5);
theta6[1] ~ normal( alpha + beta * theta6[2] , tau2);
y6 ~ multi_normal(theta6,T6);
theta7[1] ~ normal( alpha + beta * theta7[2] , tau2);
y7 ~ multi_normal(theta7,T7);
theta8[1] ~ normal( alpha + beta * theta8[2] , tau2);
y8 ~ multi_normal(theta8,T8);
theta9[1] ~ normal( alpha + beta * theta9[2] , tau2);
y9 ~ multi_normal(theta9,T9);
theta10[1] ~ normal( alpha + beta * theta10[2] , tau2);
y10 ~ multi_normal(theta10,T10);
theta11[1] ~ normal( alpha + beta * theta11[2] , tau2);
y11 ~ multi_normal(theta11,T11);
theta12[1] ~ normal( alpha + beta * theta12[2] , tau2);
y12 ~ multi_normal(theta12,T12);
theta13[1] ~ normal( alpha + beta * theta13[2] , tau2);
y13 ~ multi_normal(theta13,T13);
theta14[1] ~ normal( alpha + beta * theta14[2] , tau2);
y14 ~ multi_normal(theta14,T14);
theta15[1] ~ normal( alpha + beta * theta15[2] , tau2);
y15 ~ multi_normal(theta15,T15);
theta16[1] ~ normal( alpha + beta * theta16[2] , tau2);
y16 ~ multi_normal(theta16,T16);
theta17[1] ~ normal( alpha + beta * theta17[2] , tau2);
y17 ~ multi_normal(theta17,T17);
theta18[1] ~ normal( alpha + beta * theta18[2] , tau2);
y18 ~ multi_normal(theta18,T18);
theta19[1] ~ normal( alpha + beta * theta19[2] , tau2);
y19 ~ multi_normal(theta19,T19);
theta20[1] ~ normal( alpha + beta * theta20[2] , tau2);
y20 ~ multi_normal(theta20,T20);
theta21[1] ~ normal( alpha + beta * theta21[2] , tau2);
y21 ~ multi_normal(theta21,T21);
theta22[1] ~ normal( alpha + beta * theta22[2] , tau2);
y22 ~ multi_normal(theta22,T22);
theta23[1] ~ normal( alpha + beta * theta23[2] , tau2);
y23 ~ multi_normal(theta23,T23);
theta24[1] ~ normal( alpha + beta * theta24[2] , tau2);
y24 ~ multi_normal(theta24,T24);
theta25[1] ~ normal( alpha + beta * theta25[2] , tau2);
y25 ~ multi_normal(theta25,T25);
}"
study_data_try <- list(
J = 25,
y1 = Y1[1,],y2 = Y1[2,],y3 = Y1[3,],y4 = Y1[4,],y5 = Y1[5,],
y6 = Y1[6,],y7 = Y1[7,],y8 = Y1[8,],y9 = Y1[9,],y10 = Y1[10,],
y11 = Y1[11,],y12 = Y1[12,],y13 = Y1[13,],y14 = Y1[14,],y15 = Y1[15,],
y16 = Y1[16,], y17 = Y1[17,], y18 = Y1[18,], y19 = Y1[19,], y20 = Y1[20,],
y21 = Y1[21,], y22 = Y1[22,], y23 = Y1[23,], y24 = Y1[24,], y25 = Y1[25,],
T1=TT[,,1],T2=TT[,,2],T3=TT[,,3],T4=TT[,,4],T5=TT[,,5],
T6=TT[,,6],T7=TT[,,7],T8=TT[,,8],T9=TT[,,9],T10=TT[,,10],
T11=TT[,,11],T12=TT[,,12],T13=TT[,,13],T14=TT[,,14],T15=TT[,,15],
T16 = TT[,,16], T17 = TT[,,17], T18 = TT[,,18], T19 = TT[,,19], T20 = TT[,,20],
T21 = TT[,,21], T22 = TT[,,22], T23 = TT[,,23], T24 = TT[,,24], T25 = TT[,,25]
)
fit.meta_try <- stan(
#file = "schools.stan", # Stan program
model_code= model.def_try,
data = study_data_try, # named list of data
chains = 2, # number of Markov chains
warmup = 10000, # number of warmup iterations per chain
iter = 20000, # total number of iterations per chain
cores = 1, # number of cores (could use one per chain)
refresh = 0 , # no progress shown
control = list(adapt_delta = 0.99)
)
**It seems my mcmc chains do not converge and the results are not reliable. **
I am getting following errors
'config' variable 'CPP' is deprecated
clang -E
Warning messages:
1: In system(paste(CPP, ARGS), ignore.stdout = TRUE, ignore.stderr = TRUE) :
error in running command
2: There were 887 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
3: There were 9988 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
4: There were 2 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low
5: Examine the pairs() plot to diagnose sampling problems
6: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-ess
7: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-ess
I would really appreciate any help. I tried to include adapt_delta=0.99 as mentioned in stan manual to get rid of #divergent-transitions-after-warmup. I also beleive that there must be a better way to write this code but can anyone suggest how can I define an array of 2x2 matrix?