Hi All,
I am relatively new to Stan, and my problem concerns model parameter outcomes from Stan application.
Study:
For each patient we observe two clinical endpoints over three visits; no grouping factors. We specify two covariance matrices: 1) between the endpoints Sigma (2x2), and 2) within endpoint, Tigma (3x3); in both cases we want them to be unstructured. Our patient covariance matrix VCOV (6x6) is Kronecker product of these two.
Below you can find:
a) Part of simulated data: groups1:6 correspond to each biomarker readout row and the biormarker - visit pairs, with visit 1 as baseline.
b) The codes for model, data and sampling.
So far:
The model complied. For beta0 and Sigma the chains mixed relatively well with outcomes in line with (or acceptably close to) the simulation input. Unfortunately, for Tigma the chains did not mix well and outcome levels were rather far from the simulation input. It was not converging even if matrix T reflected exactly the ones used in simulations. As you can see Ive used Wishart priors for covariance matrices. Before Ive tried fitting the model with proc mcmc in SAS and JAGS but no success (I guess) due to the priors as Gibbs failed in recognizing conjugacy. Any ideas why Tigma fails? Is it the model specification itself, or nature of that problem due to the chosen priors? In SAS, proc mixed gave me estimates in line with the simulation inputs.
Any help or insight is much appreciated.
Robert
> head(DataSim, 10)
SUBJID CVORRSN VISIT CVTSTCD group1 group2 group3 group4 group5 group6
1 1 0.39027051 DAY 1 AK001 1 0 0 0 0 0
2 1 0.12891371 DAY 1 PASP 0 1 0 0 0 0
3 1 0.46329558 DAY 2 AK001 1 0 1 0 0 0
4 1 0.08407614 DAY 2 PASP 0 1 0 1 0 0
5 1 0.34342493 DAY 3 AK001 1 0 0 0 1 0
6 1 0.18086933 DAY 3 PASP 0 1 0 0 0 1
7 2 0.27029459 DAY 1 AK001 1 0 0 0 0 0
8 2 0.18741587 DAY 1 PASP 0 1 0 0 0 0
9 2 0.58242061 DAY 2 AK001 1 0 1 0 0 0
10 2 0.38295721 DAY 2 PASP 0 1 0 1 0 0
The model:
model{
data{
int<lower=1> N;
int<lower=1> NSUBJID;
int<lower=0, upper=1> group1[N];
int<lower=0, upper=1> group2[N];
int<lower=0, upper=1> group3[N];
int<lower=0, upper=1> group4[N];
int<lower=0, upper=1> group5[N];
int<lower=0, upper=1> group6[N];
int<lower=0> df_Sigma;
int<lower=0> df_Tigma;
int<lower=0> df_VCOV;
real y[N];
real mu1;
real s1;
real mu2;
real s2;
int<lower=1> ID[N];
cov_matrix[2] S;
cov_matrix[3] T;
}
parameters{
vector[6] beta[NSUBJID];
vector[6] beta0;
cov_matrix[2] Sigma;
cov_matrix[3] Tigma;
cov_matrix[6] VCOV;
real<lower=0> sigmaErr;
}
model{
real mu[N];
matrix[6,6] Omega;
for(i in 1:N){
mu[i] = beta[ID[i],1] * group1[i] +
beta[ID[i],2] * group2[i] +
beta[ID[i],3] * group3[i] +
beta[ID[i],4] * group4[i] +
beta[ID[i],5] * group5[i] +
beta[ID[i],6] * group6[i];
y[i] ~ normal(mu[i], sigmaErr);
}
sigmaErr ~ uniform(0, 2);
for(j in 1:NSUBJID){
beta[j, 1:6] ~ multi_normal(beta0[1:6], VCOV);
}
beta0[1] ~ normal(mu1, s1);
beta0[2] ~ normal(mu2, s2);
beta0[3] ~ normal(0, 0.1);
beta0[4] ~ normal(0, 0.1);
beta0[5] ~ normal(0, 0.1);
beta0[6] ~ normal(0, 0.1);
Sigma ~ inv_wishart(df_Sigma,S);
Tigma ~ inv_wishart(df_Tigma,T);
// Elementwise Kronecker product covariance matrix:
for (k in 1:2){
for (l in 1:2){
for (m in 1:3){
for (n in 1:3){
Omega[(k-1)*3+m, (l-1)*3+n] = Sigma[k,l] * Tigma[m,n];
}
}
}
}
VCOV ~ inv_wishart(df_VCOV, Omega);
}
Data and sampling (controls are subjective):
DataSim_Stan <- list(N=nrow(DataSim),
ID = DataSim $ SUBJID,
NSUBJID=length(unique(DataSim $ SUBJID)),
group1=DataSim $ group1,
group2=DataSim $ group2,
group3=DataSim $ group3,
group4=DataSim $ group4,
group5=DataSim $ group5,
group6=DataSim $ group6,
y = DataSim $ CVORRSN,
S = diag(0.1, 2),
T=diag(0.1,3),
mu1=0.5, s1=0.1, mu2=0.5, s2=0.1,
df_Sigma=2+1, df_Tigma=3+1, df_VCOV=6+1
)
theta_Stan = sampling(object=stan.compiled,
data=DataSim_Stan,
pars=c('beta0', 'Sigma', 'Tigma'),
warmup = 5e2,
iter = 1e3,
thin=1,
chains = 1,
init='random',
control=list(adapt_delta=0.9,
max_treedepth=15,
stepsize=1e-3,
epsilon=1e-8))
}