# Meta analysis with correlation structure between Y and X

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
)


**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?

1 Like

Hi,
I unfortunately find the model very hard to follow - presumably because you had trouble writing this in matrix format.

An array of N 2x2 matrices can be written (a bit unintuitively, following the C++ convention):

matrix[2,2] my_matrices[N];


and the order of indices for indexing would be as in:

for(n in 1:N) {
matrix[2,2] m = my_matrices[n];
}


i.e. the array index comes before the two matrix indices.

Additionally, the model as written does not make much sense to me - you specify two separate densities for \theta_i - Stan will happily do that for you, but it will usually result in weird likelihood/prior structure (which may then manifest as divergences). Could you elaborate a bit more on what is the motivation behind the model and what are you trying to achieve â€śin the real worldâ€ť? Maybe we can figure out a better way to encode your assumptions/needs into the model.

Best of luck!