I copied the Stan code. Hopefully, that give you more eyeballs.
data {
int<lower=1> N; // Total number of data points
real<lower=0> time[N]; // Time points
vector[2] y[N]; // antiPT-antiPRN levels on log scale
int<lower=1> J; // Number of subjects
int<lower=1, upper=J> subj[N]; // Subject id
int<lower=0, upper=1> Land[N]; // Country indication
int<lower=0, upper=1> C2[N]; // Child indication
}
parameters {
real<lower=0> beta1; // Basic parms for antiPT
real<lower=0> gamma1;
real<lower=4,upper=8> h;
real<lower=0> alpha1;
real<lower=0> A01;
real<lower=0> beta2; // Basic parms for antiPT
real<lower=0> gamma2;
// real<lower=4,upper=8> h; // Assume the same h
real<lower=0> alpha2;
real<lower=0> A02;
vector[J] b1_raw; // To apply non-centered parameterization
vector[J] bbeta1_raw; // Random parts for antiPT
vector[J] bgamma1_raw;
vector[J] balpha1_raw;
vector[J] b2_raw; // To apply non-centered parameterization
vector[J] bbeta2_raw; // Random parts for antiPRN
vector[J] bgamma2_raw;
vector[J] balpha2_raw;
real<lower=0> sigma1_A0; // Variance part for antiPT
real<lower=0> sigma1_beta;
real<lower=0> sigma1_gamma;
real<lower=0> sigma1_alpha;
//real<lower=0> sigma1;
real<lower=0> sigma2_A0; // Variance part for antiPRN
real<lower=0> sigma2_beta;
real<lower=0> sigma2_gamma;
real<lower=0> sigma2_alpha;
//real<lower=0> sigma2;
real v1A0PT; real v2A0PT; // Covariate part for antiPT
real v1betaPT; real v2betaPT;
real<lower=-10,upper=10> v1gammaPT;real<lower=-10,upper=10> v2gammaPT;
real<lower=-10,upper=10> v1alphaPT;real<lower=-10,upper=10> v2alphaPT;
real<lower=-10,upper=10> v1h; real<lower=-10,upper=10> v2h; // h: only v1h, v2h, same for PT and PRN
real v1A0PRN; real v2A0PRN; // Covariate part for antiPRN
real v1betaPRN; real v2betaPRN;
real<lower=-10,upper=10> v1gammaPRN;real<lower=-10,upper=10> v2gammaPRN;
real<lower=-10,upper=10> v1alphaPRN;real<lower=-10,upper=10> v2alphaPRN;
//real<lower=0> omega; // Covariance term
//cholesky_factor_corr[2] L_corr;
//vector<lower=0,upper=1>[J] sigma;
cholesky_factor_corr[2] L_corr;
vector<lower=0>[2] sigma;
}
transformed parameters{
// For antiPT
vector[J] b1; vector[J] balpha1; vector[J] bgamma1; vector[J] bbeta1;
real betapart1[N];
real gammapart1[N];
real alphapart1[N];
real A0part1[N];
real mu1[N];
real diffPT[N]; real diff2PT[N];
real hpart[N]; //hpart the same for both PT and PRN
// For antiPRN
vector[J] b2; vector[J] balpha2; vector[J] bgamma2; vector[J] bbeta2;
real betapart2[N];
real gammapart2[N];
real alphapart2[N];
real A0part2[N];
real mu2[N];
real diffPRN[N]; real diff2PRN[N];
vector[2] mu[N]; // vector of the mean
// Calculate for antiPT & antiPRN
b1 = sigma1_A0*b1_raw; bgamma1=sigma1_gamma*bgamma1_raw; bbeta1=sigma1_beta*bbeta1_raw;
balpha1=sigma1_alpha*balpha1_raw; // For antiPT
b2 = sigma2_A0*b2_raw; bgamma2=sigma2_gamma*bgamma2_raw; bbeta2=sigma2_beta*bbeta2_raw;
balpha2=sigma2_alpha*balpha2_raw; // For antiPRN
for (i in 1:N) {
A0part1[i] = b1[subj[i]] +v1A0PT*C2[i] + v2A0PT*Land[i] ;
betapart1[i] = beta1 + bbeta1[subj[i]] +v1betaPT*C2[i] + v2betaPT*Land[i] ;
gammapart1[i] = gamma1 + bgamma1[subj[i]] +v1gammaPT*C2[i] + v2gammaPT*Land[i];
alphapart1[i] = alpha1 + balpha1[subj[i]] +v1alphaPT*C2[i] + v2alphaPT*Land[i]; // For antiPT
hpart[i] = h + v1h*C2[i] +v2h*Land[i] ; // hpart the same for both
A0part2[i] = b2[subj[i]] +v1A0PRN*C2[i] + v2A0PRN*Land[i] ;
betapart2[i] = beta2 + bbeta2[subj[i]] +v1betaPRN*C2[i] + v2betaPRN*Land[i] ;
gammapart2[i] = gamma2 + bgamma2[subj[i]] +v1gammaPRN*C2[i] + v2gammaPRN*Land[i];
alphapart2[i] = alpha2 + balpha2[subj[i]] +v1alphaPRN*C2[i] + v2alphaPRN*Land[i];
// The mean
if ( time[i] <= 2 ) {
mu1[i] = log(A01*exp(A0part1[i] -betapart1[i]*time[i])) ;
mu2[i] = log(A02*exp(A0part2[i] -betapart2[i]*time[i])) ; }
else if ( time[i] >2 && time[i] <=hpart[i] ) {
mu1[i] = log(A01*exp(A0part1[i] -2*betapart1[i] + gammapart1[i]*(time[i]-2) ) ) ;
mu2[i] = log(A02*exp(A0part2[i] -2*betapart2[i] + gammapart2[i]*(time[i]-2) ) ) ;}
else if ( time[i] > hpart[i]) {
mu1[i] = log(A01*exp(A0part1[i] - 2*betapart1[i] + gammapart1[i]*(hpart[i]-2) -alphapart1[i]*(time[i]-hpart[i])) );
mu2[i] = log(A02*exp(A0part2[i] - 2*betapart2[i] + gammapart2[i]*(hpart[i]-2) -alphapart2[i]*(time[i]-hpart[i])) );}
mu[i,1] = mu1[i];
mu[i,2] = mu2[i];
}
}
model {
// Priors for antiPT
b1_raw ~ normal(0, 1); //subj random effects - As we apply non-centered parameterization
bbeta1_raw ~ normal(0, 1);
bgamma1_raw~ normal(0, 1);
balpha1_raw~ normal(0, 1);
// Priors for antiPRN
b2_raw ~ normal(0, 1); //subj random effects - As we apply non-centered parameterization
bbeta2_raw ~ normal(0, 1);
bgamma2_raw~ normal(0, 1);
balpha2_raw~ normal(0, 1);
sigma ~ cauchy(0, 5);
L_corr ~ lkj_corr_cholesky(1);
// Likelihood
for (i in 1:N){
for (k in 1:2) {
target += multi_normal_cholesky_lpdf(y[N] | mu[k], diag_pre_multiply(sigma, L_corr) );
} ;
}
}
To figure out what is going wrong it could be a good idea to start from a simpler model. For instance, a non-hierarchical model and see whether you can make that work. I also typically debug with iter=200
. Stan is good at finding the typical set if your model is well specified.
In no particular order, I noticed the following.
-
mu1[i] = log(A01*exp(A0part1[i] -betapart1[i]*time[i]))
can be rewritten as mu1[i] = log(A01) + exp(A0part1[i] -betapart1[i]*time[i]))
which is going to be numerical stabler and requires less calculations. You could even work directly with log_A01 as a parameter if you correctly adjust your prior on A01.
-
You don’t need mu1 and mu2, you can directly assign to mu[i,j].
-
I think you could speed up things if you work with matrices and vectors for your parameters. Stan treats matrices in a special way under the hood to calculate the derivatives.
-
The following has some mistakes
for (i in 1:N){
for (k in 1:2) {
target += multi_normal_cholesky_lpdf(y[N] | mu[k], diag_pre_multiply(sigma, L_corr) );
} ;
};
It should be
for (i in 1:N){
target += multi_normal_cholesky_lpdf(y[i] | mu[i], diag_pre_multiply(sigma, L_corr) )
};
I think this is why you did not get any results.
You can speed this up by
{
L_cov matrix[2,2];
L_cov = diag_pre_multiply(sigma, L_corr); //only calculate this once
target += multi_normal_cholesky_lpdf(y | mu, L_cov);
}