I’m working with two data sets from the same group of subjects. They completed a survey of 4 questions on a 5-point Likert scale .
The group was divided in two with one group receiving an intervention. Everyone then completed a different Likert scale.
I would like to use the responses from the first scale as a predictor in the model for the second scale.
I started setting up a combined model with an ordinal probit model for each scale.
I’m calculating sum-scores from the first model per row, which is a sum score for each question for each person.
How can I calculate an average sum score per person (i.e., group by subject index and take an average of calculated sum scores)?
I’d like to then use that average sum-score as the predictor in the second model. But I’m unsure how I can or should index the average to line up with the data for the second model.
Here is my stan code created by combining two brms models:
functions {
/* compute correlated group-level effects
* Args:
* z: matrix of unscaled group-level effects
* SD: vector of standard deviation parameters
* L: cholesky factor correlation matrix
* Returns:
* matrix of scaled group-level effects
*/
matrix scale_r_cor(matrix z, vector SD, matrix L) {
// r is stored in another dimension order than z
return transpose(diag_pre_multiply(SD, L) * z);
}
/* cumulative-probit log-PDF for a single response
* Args:
* y: response category
* mu: latent mean parameter
* disc: discrimination parameter
* thres: ordinal thresholds
* Returns:
* a scalar to be added to the log posterior
*/
real cumulative_probit_lpmf(int y, real mu, real disc, vector thres) {
int nthres = num_elements(thres);
real p;
if (y == 1) {
p = Phi(disc * (thres[1] - mu));
} else if (y == nthres + 1) {
p = 1 - Phi(disc * (thres[nthres] - mu));
} else {
p = Phi(disc * (thres[y] - mu)) - Phi(disc * (thres[y - 1] - mu));
}
return log(p);
}
/* compute monotonic effects
* Args:
* scale: a simplex parameter
* i: index to sum over the simplex
* Returns:
* a scalar between 0 and rows(scale)
*/
real mo(vector scale, int i) {
if (i == 0) {
return 0;
} else {
return rows(scale) * sum(scale[1 : i]);
}
}
}
data {
// Comp
int<lower=1> N_Comp; // total number of observations
array[N_Comp] int Y_Comp; // response variable
int<lower=2> nthres_Comp; // number of thresholds
// data for group-level effects of ID 1
int<lower=1> N_1_Comp; // number of grouping levels
int<lower=1> M_1_Comp; // number of coefficients per level
array[N_Comp] int<lower=1> J_1_Comp; // grouping indicator per observation
// group-level predictor values
vector[N_Comp] Z_1_1_Comp;
// Aut and Con
int<lower=1> N; // total number of observations
int<lower=1> N_Aut; // number of observations
array[N_Aut] int Y_Aut; // response variable
int<lower=2> nthres_Aut; // number of thresholds
int<lower=1> K_Aut; // number of population-level effects
matrix[N_Aut, K_Aut] X_Aut; // population-level design matrix
int<lower=1> Ksp_Aut; // number of special effects terms
int<lower=1> Imo_Aut; // number of monotonic variables
array[Imo_Aut] int<lower=1> Jmo_Aut; // length of simplexes
array[N_Aut] int Xmo_Aut_1; // monotonic variable
vector[Jmo_Aut[1]] con_simo_Aut_1; // prior concentration of monotonic simplex
int<lower=1> K_disc_Aut; // number of population-level effects
matrix[N_Aut, K_disc_Aut] X_disc_Aut; // population-level design matrix
int<lower=1> N_Con; // number of observations
array[N_Con] int Y_Con; // response variable
int<lower=2> nthres_Con; // number of thresholds
int<lower=1> K_Con; // number of population-level effects
matrix[N_Con, K_Con] X_Con; // population-level design matrix
int<lower=1> Ksp_Con; // number of special effects terms
int<lower=1> Imo_Con; // number of monotonic variables
array[Imo_Con] int<lower=1> Jmo_Con; // length of simplexes
array[N_Con] int Xmo_Con_1; // monotonic variable
vector[Jmo_Con[1]] con_simo_Con_1; // prior concentration of monotonic simplex
int<lower=1> K_disc_Con; // number of population-level effects
matrix[N_Con, K_disc_Con] X_disc_Con; // population-level design matrix
// data for group-level effects of ID 1
int<lower=1> N_1; // number of grouping levels
int<lower=1> M_1; // number of coefficients per level
array[N_Aut] int<lower=1> J_1_Aut; // grouping indicator per observation
array[N_Con] int<lower=1> J_1_Con; // grouping indicator per observation
// group-level predictor values
vector[N_Aut] Z_1_Aut_1;
vector[N_Con] Z_1_Con_2;
int<lower=1> NC_1; // number of group-level correlations
// data for group-level effects of ID 2
int<lower=1> N_2; // number of grouping levels
int<lower=1> M_2; // number of coefficients per level
array[N_Aut] int<lower=1> J_2_Aut; // grouping indicator per observation
// group-level predictor values
vector[N_Aut] Z_2_Aut_1;
// data for group-level effects of ID 3
int<lower=1> N_3; // number of grouping levels
int<lower=1> M_3; // number of coefficients per level
array[N_Con] int<lower=1> J_3_Con; // grouping indicator per observation
// group-level predictor values
vector[N_Con] Z_3_Con_1;
int prior_only; // should the likelihood be ignored?
}
parameters {
// Comp
ordered[nthres_Comp] Intercept_Comp; // temporary thresholds for centered predictors
vector<lower=0>[M_1_Comp] sd_1_Comp; // group-level standard deviations
array[M_1_Comp] vector[N_1_Comp] z_1_Comp; // standardized group-level effects
//Aut and Con
vector[K_Aut] b_Aut; // population-level effects
ordered[nthres_Aut] Intercept_Aut; // temporary thresholds for centered predictors
simplex[Jmo_Aut[1]] simo_Aut_1; // monotonic simplex
vector[Ksp_Aut] bsp_Aut; // special effects coefficients
vector[K_disc_Aut] b_disc_Aut; // population-level effects
vector[K_Con] b_Con; // population-level effects
ordered[nthres_Con] Intercept_Con; // temporary thresholds for centered predictors
simplex[Jmo_Con[1]] simo_Con_1; // monotonic simplex
vector[Ksp_Con] bsp_Con; // special effects coefficients
vector[K_disc_Con] b_disc_Con; // population-level effects
vector<lower=0>[M_1] sd_1; // group-level standard deviations
matrix[M_1, N_1] z_1; // standardized group-level effects
cholesky_factor_corr[M_1] L_1; // cholesky factor of correlation matrix
vector<lower=0>[M_2] sd_2; // group-level standard deviations
array[M_2] vector[N_2] z_2; // standardized group-level effects
vector<lower=0>[M_3] sd_3; // group-level standard deviations
array[M_3] vector[N_3] z_3; // standardized group-level effects
real b_PCL; // beta for PCL sum-scores
}
transformed parameters {
//Comp
vector[N_Comp] mu_Comp_pred = rep_vector(0, N_Comp); // linear predictor
//vector[nmthres] merged_Intercept; // merged thresholds
real disc = 1; // discrimination parameters
vector[N_1_Comp] r_1_1_Comp; // actual group-level effects
vector[N_Comp] p_1; // probability of choice 1
vector[N_Comp] p_2; // probability of choice 2
vector[N_Comp] p_3; // probability of choice 3
vector[N_Comp] p_4; // probability of choice 4
vector[N_Comp] p_5; // probability of choice 5
vector[N_Comp] Sum_Score; // probability of choice 5
//Aut and Con
matrix[N_1, M_1] r_1; // actual group-level effects
// using vectors speeds up indexing in loops
vector[N_1] r_1_Aut_1;
vector[N_1] r_1_Con_2;
vector[N_2] r_2_Aut_1; // actual group-level effects
vector[N_3] r_3_Con_1; // actual group-level effects
// Comp merge thresholds
r_1_1_Comp = sd_1_Comp[1] * z_1_Comp[1];
// Comp response probabilities
mu_Comp_pred += r_1_1_Comp[J_1_Comp];
for (n in 1 : N_Comp) {
p_1[n] = inv_logit(Intercept_Comp[1] - (mu_Comp_pred[n]));
p_2[n] = inv_logit(Intercept_Comp[2] - mu_Comp_pred[n]) - p_1[n];
p_3[n] = inv_logit(Intercept_Comp[3] - mu_Comp_pred[n]) - (p_1[n] + p_2[n]);
p_4[n] = inv_logit(Intercept_Comp[4] - mu_Comp_pred[n]) - (p_1[n] + p_2[n] + p_3[n]);
p_5[n] = 1 - (p_1[n] + p_2[n] + p_3[n] + p_4[n]);
Sum_Score[n] = p_1[n]*1 + p_2[n]*2 + p_3[n]*3 + p_4[n]*4 + p_5[n]*5;
}
// for (n in 1 : N_1_Comp) {
// Sum_Score[n] = p_1[J_1_Comp[n]]*1 + p_2[J_1_Comp[n]]*2 + p_3[J_1_Comp[n]]*3 + p_4[J_1_Comp[n]]*4 + p_5[J_1_Comp[n]]*5;
// }
// Aut and Con
// compute actual group-level effects
r_1 = scale_r_cor(z_1, sd_1, L_1);
r_1_Aut_1 = r_1[ : , 1];
r_1_Con_2 = r_1[ : , 2];
r_2_Aut_1 = sd_2[1] * z_2[1];
r_3_Con_1 = sd_3[1] * z_3[1];
}
model {
// Comp
real lprior_Comp = 0; // prior contributions to the log posterior
// priors not including constants
lprior_Comp += normal_lupdf(Intercept_Comp | 0, 4);
// lprior_Comp += normal_lupdf(Intercept_Comp_1 | 0, 4);
// lprior_Comp += normal_lupdf(Intercept_Comp_2 | 0, 4);
// lprior_Comp += normal_lupdf(Intercept_Comp_3 | 0, 4);
// lprior_Comp += normal_lupdf(Intercept_Comp_4 | 0, 4);
lprior_Comp += normal_lupdf(sd_1_Comp | 0, 1);
target += lprior_Comp;
target += std_normal_lupdf(z_1_Comp[1]);
// Aut and Con
real lprior = 0; // prior contributions to the log posterior
// priors not including constants
lprior += normal_lupdf(b_Aut[1] | 0, 2);
lprior += normal_lupdf(b_Aut[3] | 0, 2);
lprior += normal_lupdf(b_Aut[4] | 0, 2);
lprior += normal_lupdf(Intercept_Aut | 0, 4);
lprior += dirichlet_lupdf(simo_Aut_1 | con_simo_Aut_1);
lprior += normal_lupdf(bsp_Aut[1] | 0, 2);
lprior += normal_lupdf(b_disc_Aut | 2, 2);
lprior += normal_lupdf(b_Con[1] | 0, 1);
lprior += normal_lupdf(b_Con[3] | 0, 1);
lprior += normal_lupdf(b_Con[4] | 0, 1);
lprior += normal_lupdf(Intercept_Con | 0, 4);
lprior += dirichlet_lupdf(simo_Con_1 | con_simo_Con_1);
lprior += normal_lupdf(bsp_Con[1] | 0, 1);
lprior += normal_lupdf(b_disc_Con | 2, 2);
lprior += normal_lupdf(sd_1[1] | 0, 1);
lprior += normal_lupdf(sd_1[2] | 0, 1);
lprior += lkj_corr_cholesky_lupdf(L_1 | 1);
lprior += normal_lupdf(sd_2[1] | 0, 1);
lprior += normal_lupdf(sd_3[1] | 0, 1);
target += lprior;
target += std_normal_lupdf(to_vector(z_1));
target += std_normal_lupdf(z_2[1]);
target += std_normal_lupdf(z_3[1]);
b_PCL ~ normal(0,1);
// Comp
// likelihood not including constants
if (!prior_only) {
// initialize linear predictor term
vector[N_Comp] mu_Comp = rep_vector(0.0, N_Comp);
for (n in 1 : N_Comp) {
// add more terms to the linear predictor
mu_Comp[n] += r_1_1_Comp[J_1_Comp[n]] * Z_1_1_Comp[n];
}
for (n in 1 : N_Comp) {
target += cumulative_probit_lpmf(Y_Comp[n] | mu_Comp[n], disc, Intercept_Comp);
}
}
// Aut and Con
// likelihood not including constants
if (!prior_only) {
// initialize linear predictor term
vector[N_Aut] mu_Aut = rep_vector(0.0, N_Aut);
// initialize linear predictor term
vector[N_Aut] disc_Aut = rep_vector(0.0, N_Aut);
// initialize linear predictor term
vector[N_Con] mu_Con = rep_vector(0.0, N_Con);
// initialize linear predictor term
vector[N_Con] disc_Con = rep_vector(0.0, N_Con);
mu_Aut += X_Aut * b_Aut;
disc_Aut += X_disc_Aut * b_disc_Aut;
mu_Con += X_Con * b_Con;
disc_Con += X_disc_Con * b_disc_Con;
for (n in 1 : N_Aut) {
// add more terms to the linear predictor
mu_Aut[n] += bsp_Aut[1] * mo(simo_Aut_1, Xmo_Aut_1[n])
//mu_Aut[n] += bsp_Aut[1] * mo(simo_Aut_1, Xmo_Aut_1[n]) + b_PCL*Sum_Score[J_1_Aut[n]]
+ r_1_Aut_1[J_1_Aut[n]] * Z_1_Aut_1[n]
+ r_2_Aut_1[J_2_Aut[n]] * Z_2_Aut_1[n];
}
for (n in 1 : N_Con) {
// add more terms to the linear predictor
mu_Con[n] += bsp_Con[1] * mo(simo_Con_1, Xmo_Con_1[n])
//mu_Con[n] += bsp_Con[1] * mo(simo_Con_1, Xmo_Con_1[n]) + b_PCL*Sum_Score[J_1_Con[n]]
+ r_1_Con_2[J_1_Con[n]] * Z_1_Con_2[n]
+ r_3_Con_1[J_3_Con[n]] * Z_3_Con_1[n];
}
disc_Aut = exp(disc_Aut);
disc_Con = exp(disc_Con);
for (n in 1 : N_Aut) {
target += cumulative_probit_lpmf(Y_Aut[n] | mu_Aut[n], disc_Aut[n], Intercept_Aut);
}
for (n in 1 : N_Con) {
target += cumulative_probit_lpmf(Y_Con[n] | mu_Con[n], disc_Con[n], Intercept_Con);
}
}
}
generated quantities {
// compute actual thresholds
vector[nthres_Aut] b_Aut_Intercept = Intercept_Aut;
// compute actual thresholds
vector[nthres_Con] b_Con_Intercept = Intercept_Con;
// compute group-level correlations
corr_matrix[M_1] Cor_1 = multiply_lower_tri_self_transpose(L_1);
vector<lower=-1, upper=1>[NC_1] cor_1;
// extract upper diagonal of correlation matrix
for (k in 1 : M_1) {
for (j in 1 : (k - 1)) {
cor_1[choose(k - 1, 2) + j] = Cor_1[j, k];
}
}
}
Do I need to index the sum-score calculations to J_1_Comp?
Right now, I get a Sum_Score for each observation at each iteration (8000 x 496), where groups of 4 columns are from 1 subject, so I’d like to average groups of 4 Sum_Scores per iteration (8000 x 124). Then, use the 1:124 as indexes for the second model.