Calculating per subject average from one model to use as predictor in another

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.