And here is the stancode:
// generated with brms 2.13.0
functions {
}
data {
int<lower=1> N; // number of observations
int<lower=2> ncat; // number of categories
int Y[N]; // response variable
vector<lower=0>[N] weights; // model weights
int<lower=1> K_muAS; // number of population-level effects
matrix[N, K_muAS] X_muAS; // population-level design matrix
vector[N] offsets_muAS;
int<lower=1> K_muICE; // number of population-level effects
matrix[N, K_muICE] X_muICE; // population-level design matrix
vector[N] offsets_muICE;
int<lower=1> K_muSACCZ; // number of population-level effects
matrix[N, K_muSACCZ] X_muSACCZ; // population-level design matrix
vector[N] offsets_muSACCZ;
int<lower=1> K_muNACCZ; // number of population-level effects
matrix[N, K_muNACCZ] X_muNACCZ; // population-level design matrix
vector[N] offsets_muNACCZ;
int<lower=1> K_muNSAF; // number of population-level effects
matrix[N, K_muNSAF] X_muNSAF; // population-level design matrix
vector[N] offsets_muNSAF;
// 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
int<lower=1> J_1[N]; // grouping indicator per observation
// group-level predictor values
vector[N] Z_1_muAS_1;
vector[N] Z_1_muICE_2;
vector[N] Z_1_muSACCZ_3;
vector[N] Z_1_muNACCZ_4;
vector[N] Z_1_muNSAF_5;
int<lower=1> NC_1; // number of group-level correlations
int prior_only; // should the likelihood be ignored?
}
transformed data {
int Kc_muAS = K_muAS - 1;
matrix[N, Kc_muAS] Xc_muAS; // centered version of X_muAS without an intercept
vector[Kc_muAS] means_X_muAS; // column means of X_muAS before centering
int Kc_muICE = K_muICE - 1;
matrix[N, Kc_muICE] Xc_muICE; // centered version of X_muICE without an intercept
vector[Kc_muICE] means_X_muICE; // column means of X_muICE before centering
int Kc_muSACCZ = K_muSACCZ - 1;
matrix[N, Kc_muSACCZ] Xc_muSACCZ; // centered version of X_muSACCZ without an intercept
vector[Kc_muSACCZ] means_X_muSACCZ; // column means of X_muSACCZ before centering
int Kc_muNACCZ = K_muNACCZ - 1;
matrix[N, Kc_muNACCZ] Xc_muNACCZ; // centered version of X_muNACCZ without an intercept
vector[Kc_muNACCZ] means_X_muNACCZ; // column means of X_muNACCZ before centering
int Kc_muNSAF = K_muNSAF - 1;
matrix[N, Kc_muNSAF] Xc_muNSAF; // centered version of X_muNSAF without an intercept
vector[Kc_muNSAF] means_X_muNSAF; // column means of X_muNSAF before centering
for (i in 2:K_muAS) {
means_X_muAS[i - 1] = mean(X_muAS[, i]);
Xc_muAS[, i - 1] = X_muAS[, i] - means_X_muAS[i - 1];
}
for (i in 2:K_muICE) {
means_X_muICE[i - 1] = mean(X_muICE[, i]);
Xc_muICE[, i - 1] = X_muICE[, i] - means_X_muICE[i - 1];
}
for (i in 2:K_muSACCZ) {
means_X_muSACCZ[i - 1] = mean(X_muSACCZ[, i]);
Xc_muSACCZ[, i - 1] = X_muSACCZ[, i] - means_X_muSACCZ[i - 1];
}
for (i in 2:K_muNACCZ) {
means_X_muNACCZ[i - 1] = mean(X_muNACCZ[, i]);
Xc_muNACCZ[, i - 1] = X_muNACCZ[, i] - means_X_muNACCZ[i - 1];
}
for (i in 2:K_muNSAF) {
means_X_muNSAF[i - 1] = mean(X_muNSAF[, i]);
Xc_muNSAF[, i - 1] = X_muNSAF[, i] - means_X_muNSAF[i - 1];
}
}
parameters {
vector[Kc_muAS] b_muAS; // population-level effects
real Intercept_muAS; // temporary intercept for centered predictors
vector[Kc_muICE] b_muICE; // population-level effects
real Intercept_muICE; // temporary intercept for centered predictors
vector[Kc_muSACCZ] b_muSACCZ; // population-level effects
real Intercept_muSACCZ; // temporary intercept for centered predictors
vector[Kc_muNACCZ] b_muNACCZ; // population-level effects
real Intercept_muNACCZ; // temporary intercept for centered predictors
vector[Kc_muNSAF] b_muNSAF; // population-level effects
real Intercept_muNSAF; // temporary intercept for centered predictors
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
}
transformed parameters {
matrix[N_1, M_1] r_1; // actual group-level effects
// using vectors speeds up indexing in loops
vector[N_1] r_1_muAS_1;
vector[N_1] r_1_muICE_2;
vector[N_1] r_1_muSACCZ_3;
vector[N_1] r_1_muNACCZ_4;
vector[N_1] r_1_muNSAF_5;
// compute actual group-level effects
r_1 = (diag_pre_multiply(sd_1, L_1) * z_1)';
r_1_muAS_1 = r_1[, 1];
r_1_muICE_2 = r_1[, 2];
r_1_muSACCZ_3 = r_1[, 3];
r_1_muNACCZ_4 = r_1[, 4];
r_1_muNSAF_5 = r_1[, 5];
}
model {
// initialize linear predictor term
vector[N] muAS = Intercept_muAS + Xc_muAS * b_muAS + offsets_muAS;
// initialize linear predictor term
vector[N] muICE = Intercept_muICE + Xc_muICE * b_muICE + offsets_muICE;
// initialize linear predictor term
vector[N] muSACCZ = Intercept_muSACCZ + Xc_muSACCZ * b_muSACCZ + offsets_muSACCZ;
// initialize linear predictor term
vector[N] muNACCZ = Intercept_muNACCZ + Xc_muNACCZ * b_muNACCZ + offsets_muNACCZ;
// initialize linear predictor term
vector[N] muNSAF = Intercept_muNSAF + Xc_muNSAF * b_muNSAF + offsets_muNSAF;
// linear predictor matrix
vector[ncat] mu[N];
for (n in 1:N) {
// add more terms to the linear predictor
muAS[n] += r_1_muAS_1[J_1[n]] * Z_1_muAS_1[n];
}
for (n in 1:N) {
// add more terms to the linear predictor
muICE[n] += r_1_muICE_2[J_1[n]] * Z_1_muICE_2[n];
}
for (n in 1:N) {
// add more terms to the linear predictor
muSACCZ[n] += r_1_muSACCZ_3[J_1[n]] * Z_1_muSACCZ_3[n];
}
for (n in 1:N) {
// add more terms to the linear predictor
muNACCZ[n] += r_1_muNACCZ_4[J_1[n]] * Z_1_muNACCZ_4[n];
}
for (n in 1:N) {
// add more terms to the linear predictor
muNSAF[n] += r_1_muNSAF_5[J_1[n]] * Z_1_muNSAF_5[n];
}
for (n in 1:N) {
mu[n] = [0, muAS[n], muICE[n], muSACCZ[n], muNACCZ[n], muNSAF[n]]';
}
// priors including all constants
target += student_t_lpdf(Intercept_muAS | 3, 0, 2.5);
target += student_t_lpdf(Intercept_muICE | 3, 0, 2.5);
target += student_t_lpdf(Intercept_muSACCZ | 3, 0, 2.5);
target += student_t_lpdf(Intercept_muNACCZ | 3, 0, 2.5);
target += student_t_lpdf(Intercept_muNSAF | 3, 0, 2.5);
target += student_t_lpdf(sd_1 | 3, 0, 2.5)
- 5 * student_t_lccdf(0 | 3, 0, 2.5);
target += std_normal_lpdf(to_vector(z_1));
target += lkj_corr_cholesky_lpdf(L_1 | 1);
// likelihood including all constants
if (!prior_only) {
for (n in 1:N) {
target += weights[n] * (categorical_logit_lpmf(Y[n] | mu[n]));
}
}
}
generated quantities {
// actual population-level intercept
real b_muAS_Intercept = Intercept_muAS - dot_product(means_X_muAS, b_muAS);
// actual population-level intercept
real b_muICE_Intercept = Intercept_muICE - dot_product(means_X_muICE, b_muICE);
// actual population-level intercept
real b_muSACCZ_Intercept = Intercept_muSACCZ - dot_product(means_X_muSACCZ, b_muSACCZ);
// actual population-level intercept
real b_muNACCZ_Intercept = Intercept_muNACCZ - dot_product(means_X_muNACCZ, b_muNACCZ);
// actual population-level intercept
real b_muNSAF_Intercept = Intercept_muNSAF - dot_product(means_X_muNSAF, b_muNSAF);
// 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;
// additionally draw samples from priors
real prior_Intercept_muAS = student_t_rng(3,0,2.5);
real prior_Intercept_muICE = student_t_rng(3,0,2.5);
real prior_Intercept_muSACCZ = student_t_rng(3,0,2.5);
real prior_Intercept_muNACCZ = student_t_rng(3,0,2.5);
real prior_Intercept_muNSAF = student_t_rng(3,0,2.5);
real prior_sd_1 = student_t_rng(3,0,2.5);
real prior_cor_1 = lkj_corr_rng(M_1,1)[1, 2];
// 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];
}
}
// use rejection sampling for truncated priors
while (prior_sd_1 < 0) {
prior_sd_1 = student_t_rng(3,0,2.5);
}
}
The offset is indeed specified in the stancode. But I don’t see any difference between a model with offset and one without. Is it because the offset variable is a proportion (between 0 and 1)?
Thanks!