Hi,
I first want to say that these forums have been very helpful for me, who is a very novice user to both stan, and statistical modeling. I had a slow running model (40 days) that with some simple tweeks went down to 4 days. With that being said, given my inexperience, my question is likely extremely basic. I want to ask the devs rather than self-teach due to a very tight work deadline.
I have a 3 level hierarchical mixed effects model that represents disease progression of Alzheimer’s. There is inter-study variability (19 studies) and inter-subject variability (4485 individuals) with about 25,0000 total observations. I need to estimate the covariance between the inter-study variability on the baseline component of the model and the rate component of the model. The same goes for the inter-subject variability. I’m not sure how to do that. Can someone advise going forward? Thanks.
Ps. If there are other noticeable features of my code that is not very well optimized or blatantly wrong, feel free to critique.
Here is my code:
data{
int N; // number of total Observations
int P; // number of subjects
int M; // number of studies
int IDp[N]; // Patient ID
int IDs[N]; // Study ID
int SEX[N];
int AGE[N];
int COMED[N];
real APOE4[N];
vector[N] time; // time of observation (years)
real S[N]; // measured ADAScog score
}
parameters{
//Fixed effects and covariates
real<lower=0,upper=1> theta_S0; // Mean baseline score for average individual
real theta_r; // Mean progression rate for average individual
real<lower=0> theta_SEX;
real theta_AGE;
real<lower=0> theta_APOE4_b;
real theta_APOE4_r;
real theta_COMED;
real<lower=0> tau;
real<lower=0> beta;
//Inter-patient re
vector[P] eta_pb; // re for patient baseline
vector[P] eta_pr; // re for patient rate
vector[M] eta_sb;
vector[M] eta_sr;
real<lower=0> omega_pb; // std for patient baseline **Want to estimate covariance between
omega_pb and omega_pr as well as omega_sb and omega_sr**
real<lower=0> omega_pr; // std for patient rate
real<lower=0> omega_sb;
real<lower=0> omega_sr;
}
transformed parameters{
vector[N] baseline_cov;
vector[N] rate_cov;
vector[N] S0;
vector[N] r;
vector<lower=0,upper=1>[N] muS;
vector[N] ones;
for(i in 1:N) ones[i] = 1;
for(i in 1:N) {
baseline_cov[i] = theta_S0*(1+theta_SEX*(SEX[i]-1))*(1+theta_APOE4_b*(APOE4[i]-.72));
rate_cov[i] = theta_r*(1+theta_AGE*(AGE[i]-75))*(1+theta_APOE4_r*(APOE4[i]-.72))*(1+theta_COMED*(1-COMED[i]));
S0[i] = baseline_cov[i]*exp(eta_pb[IDp[i]] + eta_sb[IDs[i]]);
r[i] = rate_cov[i] + eta_pr[IDp[i]] + eta_sr[IDs[i]];
muS[i] = S0[i]/(S0[i]^beta +(1-S0[i]^beta)*exp(-beta*r[i]*time[i]))^(1/beta);
}
}
model{
//Priors
omega_pb~normal(0,1); //orig uniform(0,1)
omega_pr~normal(0,1); //orig uniform(0,1)
omega_sb~normal(0,1); //orig uniform(0,1)
omega_sr~normal(0,1); //orig uniform(0,1)
eta_pb~normal(0,omega_pb);
eta_pr~normal(0,omega_pr);
eta_sb~normal(0,omega_sb);
eta_sr~normal(0,omega_sr);
theta_S0~normal(.324,1);
theta_r~normal(0,1);
theta_SEX~normal(0,1);
theta_AGE~normal(0,1);
theta_APOE4_b~normal(0,1);
theta_APOE4_r~normal(0,1);
theta_COMED~normal(0,10);
tau~normal(80,8); //orig uniform(0,100)
beta~normal(10,2);
// Likelihood
S ~beta(muS*tau, (1-muS)*tau);
}
[edit: escaped code]