Hi all,

First off, I’m fairly new to stan, so the answer to my question may be fairly obvious… to others. I’ve been working on this for a few days now, so I figured it was worth a shot to ask around for some help.

I’ve created a model for estimating two different item scores that each contain repeated measurements. This design is unbalanced, such that a given individual has a different number of measurements for the first and second items. I work with psychophysiological data, so unbalanced measurements for items (trial types) is very commonplace after rejecting “bad” measurements with too much noise/artifact.

Right now, this model estimates the within-person standard deviations (sig_err) for each item as well as the between-person standard deviations (sig_id) for each item.

I’ve been trying to build off this model to compute the between-person covariance for the mean item scores (u_1, u_2). I have been pretty unsuccessful. Does anyone have any recommendations about how I can go about computing covariance between u_1 and u_2 in this model?

Any help would be greatly appreciated!

Thanks,

Peter

stan_model <- "

data {

int<lower=1> NOB1; //total observations of first measurement

int<lower=1> NOB2; //total observations of second measurement

int<lower=1> NSUB; //total number of subjects

int<lower=1, upper=NSUB> id1[NOB1]; //grouping id for first meas

int<lower=1, upper=NSUB> id2[NOB2]; //grouping id for second meas

vector[NOB1] meas1; //scores for first measurement

vector[NOB2] meas2; //scores for second measurement

}

parameters {

vector[2] mu; //intercepts for meas1 and meas2

vector<lower=0>[2] sig_id; //bet-person std dev for given meas

vector<lower=0>[2] sig_err; //within-person std dev for given meas

vector[NSUB] id_raw1;

vector[NSUB] id_raw2;

}

transformed parameters {

vector[NSUB] u_1; //predicted meas1 score

vector[NSUB] u_2; //predicted meas2 score

u_1 = mu[1] + id_raw1*sig_id[1];
u_2 = mu[2] + id_raw2*sig_id[2];

}

model {

id_raw1 ~ normal(0,1);

id_raw2 ~ normal(0,1);

mu ~ normal(0,10);

sig_id ~ cauchy(0,10);

sig_err ~ cauchy(0,10);

meas1 ~ normal(u_1[id1],sig_err[1]);

meas2 ~ normal(u_2[id2],sig_err[2]);

}

"