Dear Max (or anyone else). I didn’t had the opportunity earlier to take a next step in this. Last week I’ve been trying to figure it out how to do this in Stan code, but I’m stuck… Can’t work my head around it.
My model input looks like this in brms:
Model <- brm(
bf(Score_1 ~ 1 + Gender) +
bf(Score_2 ~ 1 + Gender) +
bf(Score_3 ~ 1 + Gender) +
set_rescor(rescor = TRUE),
data = Data)
make_stancode() results in
> make_stancode( bf(Score_1 ~ 1 + Gender) +
+ bf(Score_2 ~ 1 + Gender) +
+ bf(Score_3 ~ 1 + Gender) +
+ set_rescor(rescor = TRUE),
+ data = Data)
// generated with brms 2.15.0
functions {
}
data {
int<lower=1> N; // total number of observations
int<lower=1> N_Score1; // number of observations
vector[N_Score1] Y_Score1; // response variable
int<lower=1> K_Score1; // number of population-level effects
matrix[N_Score1, K_Score1] X_Score1; // population-level design matrix
int<lower=1> N_Score2; // number of observations
vector[N_Score2] Y_Score2; // response variable
int<lower=1> K_Score2; // number of population-level effects
matrix[N_Score2, K_Score2] X_Score2; // population-level design matrix
int<lower=1> N_Score3; // number of observations
vector[N_Score3] Y_Score3; // response variable
int<lower=1> K_Score3; // number of population-level effects
matrix[N_Score3, K_Score3] X_Score3; // population-level design matrix
int<lower=1> nresp; // number of responses
int nrescor; // number of residual correlations
int prior_only; // should the likelihood be ignored?
}
transformed data {
int Kc_Score1 = K_Score1 - 1;
matrix[N_Score1, Kc_Score1] Xc_Score1; // centered version of X_Score1 without an intercept
vector[Kc_Score1] means_X_Score1; // column means of X_Score1 before centering
int Kc_Score2 = K_Score2 - 1;
matrix[N_Score2, Kc_Score2] Xc_Score2; // centered version of X_Score2 without an intercept
vector[Kc_Score2] means_X_Score2; // column means of X_Score2 before centering
int Kc_Score3 = K_Score3 - 1;
matrix[N_Score3, Kc_Score3] Xc_Score3; // centered version of X_Score3 without an intercept
vector[Kc_Score3] means_X_Score3; // column means of X_Score3 before centering
vector[nresp] Y[N]; // response array
for (i in 2:K_Score1) {
means_X_Score1[i - 1] = mean(X_Score1[, i]);
Xc_Score1[, i - 1] = X_Score1[, i] - means_X_Score1[i - 1];
}
for (i in 2:K_Score2) {
means_X_Score2[i - 1] = mean(X_Score2[, i]);
Xc_Score2[, i - 1] = X_Score2[, i] - means_X_Score2[i - 1];
}
for (i in 2:K_Score3) {
means_X_Score3[i - 1] = mean(X_Score3[, i]);
Xc_Score3[, i - 1] = X_Score3[, i] - means_X_Score3[i - 1];
}
for (n in 1:N) {
Y[n] = transpose([Y_Score1[n], Y_Score2[n], Y_Score3[n]]);
}
}
parameters {
vector[Kc_Score1] b_Score1; // population-level effects
real Intercept_Score1; // temporary intercept for centered predictors
real<lower=0> sigma_Score1; // residual SD
vector[Kc_Score2] b_Score2; // population-level effects
real Intercept_Score2; // temporary intercept for centered predictors
real<lower=0> sigma_Score2; // residual SD
vector[Kc_Score3] b_Score3; // population-level effects
real Intercept_Score3; // temporary intercept for centered predictors
real<lower=0> sigma_Score3; // residual SD
cholesky_factor_corr[nresp] Lrescor; // parameters for multivariate linear models
}
transformed parameters {
}
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[N_Score1] mu_Score1 = Intercept_Score1 + Xc_Score1 * b_Score1;
// initialize linear predictor term
vector[N_Score2] mu_Score2 = Intercept_Score2 + Xc_Score2 * b_Score2;
// initialize linear predictor term
vector[N_Score3] mu_Score3 = Intercept_Score3 + Xc_Score3 * b_Score3;
// multivariate predictor array
vector[nresp] Mu[N];
vector[nresp] sigma = transpose([sigma_Score1, sigma_Score2, sigma_Score3]);
// cholesky factor of residual covariance matrix
matrix[nresp, nresp] LSigma = diag_pre_multiply(sigma, Lrescor);
// combine univariate parameters
for (n in 1:N) {
Mu[n] = transpose([mu_Score1[n], mu_Score2[n], mu_Score3[n]]);
}
target += multi_normal_cholesky_lpdf(Y | Mu, LSigma);
}
// priors including constants
target += student_t_lpdf(Intercept_Score1 | 3, 10, 3);
target += student_t_lpdf(sigma_Score1 | 3, 0, 3)
- 1 * student_t_lccdf(0 | 3, 0, 3);
target += student_t_lpdf(Intercept_Score2 | 3, 11, 4.4);
target += student_t_lpdf(sigma_Score2 | 3, 0, 4.4)
- 1 * student_t_lccdf(0 | 3, 0, 4.4);
target += student_t_lpdf(Intercept_Score3 | 3, 14, 4.4);
target += student_t_lpdf(sigma_Score3 | 3, 0, 4.4)
- 1 * student_t_lccdf(0 | 3, 0, 4.4);
target += lkj_corr_cholesky_lpdf(Lrescor | 1);
}
generated quantities {
// actual population-level intercept
real b_Score1_Intercept = Intercept_Score1 - dot_product(means_X_Score1, b_Score1);
// actual population-level intercept
real b_Score2_Intercept = Intercept_Score2 - dot_product(means_X_Score2, b_Score2);
// actual population-level intercept
real b_Score3_Intercept = Intercept_Score3 - dot_product(means_X_Score3, b_Score3);
// residual correlations
corr_matrix[nresp] Rescor = multiply_lower_tri_self_transpose(Lrescor);
vector<lower=-1,upper=1>[nrescor] rescor;
// extract upper diagonal of correlation matrix
for (k in 1:nresp) {
for (j in 1:(k - 1)) {
rescor[choose(k - 1, 2) + j] = Rescor[j, k];
}
}
}
Where to start to be able to fix the effect of Gender to 1 parameter estimate rather than 3? I would really appreciate any help here. Many thanks!