Hi all,

I’m wondering if anyone can point me towards some resources for translating models into accurate statistical notation.

I’m working on translating an ordered-probit model (discussed in this thread) into statistical notation, but am not at all confident what I have is correct.

What I have so far:

Y_{sm} = \gamma_{sm} \Leftrightarrow \theta_{m,\gamma_{sm}} < \widetilde Y_{sm} \leq \theta_{m,\gamma_{sm}}, \gamma_{sm} \equiv {1,...,K_{m}}

where \gamma_{sm} is a category out of K_{m} ordered categories and \theta_{m} is a vector of threshold parameters for outcome m.

My brms model is:

```
Mod_fmla <- bf(mvbind(Valence, Arousal) | resp_thres(8, gr = Category) ~ Category +
Condition +
Category:Condition +
Intensity.c +
(1|o|Token) +
(1|p|Condition) +
(1|q|Condition:Subject)) +
lf(disc ~ 0 + Category, resp = "Valence", cmc = FALSE) +
lf(disc ~ 0 + Category, resp = "Arousal", cmc = FALSE)
## priors
priors <- c(
set_prior("normal(0,6)",
class = "Intercept",
coef = "",
group = "Neutral",
resp = ""),
set_prior("normal(0,6)",
class = "Intercept",
coef = "",
group = "Pleasant",
resp = ""),
set_prior("student_t(6,1,1.5)",
class = "b",
coef = "",
resp = ""),
set_prior("normal(0,2)",
class = "b",
coef = "CategoryNeutral",
resp = "",
dpar = "disc"),
set_prior("normal(0,2)",
class = "b",
coef = "CategoryPleasant",
resp = "",
dpar = "disc"),
set_prior("exponential(2)",
class = "sd",
coef = "Intercept",
group = "Condition",
resp = ""),
set_prior("exponential(3)",
class = "sd",
coef = "Intercept",
group = "Condition:Subject",
resp = ""),
set_prior("exponential(2)",
class = "sd",
coef = "Intercept",
group = "Token",
resp = "")
)
```

The stan code is:

```
parameters {
vector[K_Valence] b_Valence; // population-level effects
ordered[nthres_Valence[1]] Intercept_Valence_1; // temporary thresholds for centered predictors
ordered[nthres_Valence[2]] Intercept_Valence_2; // temporary thresholds for centered predictors
ordered[nthres_Valence[3]] Intercept_Valence_3; // temporary thresholds for centered predictors
vector[K_disc_Valence] b_disc_Valence; // population-level effects
vector[K_Arousal] b_Arousal; // population-level effects
ordered[nthres_Arousal[1]] Intercept_Arousal_1; // temporary thresholds for centered predictors
ordered[nthres_Arousal[2]] Intercept_Arousal_2; // temporary thresholds for centered predictors
ordered[nthres_Arousal[3]] Intercept_Arousal_3; // temporary thresholds for centered predictors
vector[K_disc_Arousal] b_disc_Arousal; // 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
matrix[M_2, N_2] z_2; // standardized group-level effects
cholesky_factor_corr[M_2] L_2; // cholesky factor of correlation matrix
vector<lower=0>[M_3] sd_3; // group-level standard deviations
matrix[M_3, N_3] z_3; // standardized group-level effects
cholesky_factor_corr[M_3] L_3; // 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_Valence_1;
vector[N_1] r_1_Arousal_2;
matrix[N_2, M_2] r_2; // actual group-level effects
// using vectors speeds up indexing in loops
vector[N_2] r_2_Valence_1;
vector[N_2] r_2_Arousal_2;
matrix[N_3, M_3] r_3; // actual group-level effects
// using vectors speeds up indexing in loops
vector[N_3] r_3_Valence_1;
vector[N_3] r_3_Arousal_2;
// compute actual group-level effects
r_1 = (diag_pre_multiply(sd_1, L_1) * z_1)';
r_1_Valence_1 = r_1[, 1];
r_1_Arousal_2 = r_1[, 2];
// compute actual group-level effects
r_2 = (diag_pre_multiply(sd_2, L_2) * z_2)';
r_2_Valence_1 = r_2[, 1];
r_2_Arousal_2 = r_2[, 2];
// compute actual group-level effects
r_3 = (diag_pre_multiply(sd_3, L_3) * z_3)';
r_3_Valence_1 = r_3[, 1];
r_3_Arousal_2 = r_3[, 2];
}
model {
vector[nmthres_Valence] merged_Intercept_Valence; // merged thresholds
// initialize linear predictor term
vector[N_Valence] mu_Valence = X_Valence * b_Valence;
// initialize linear predictor term
vector[N_Valence] disc_Valence = X_disc_Valence * b_disc_Valence;
vector[nmthres_Arousal] merged_Intercept_Arousal; // merged thresholds
// initialize linear predictor term
vector[N_Arousal] mu_Arousal = X_Arousal * b_Arousal;
// initialize linear predictor term
vector[N_Arousal] disc_Arousal = X_disc_Arousal * b_disc_Arousal;
merged_Intercept_Valence[Kthres_start_Valence[1]:Kthres_end_Valence[1]] = Intercept_Valence_1;
merged_Intercept_Valence[Kthres_start_Valence[2]:Kthres_end_Valence[2]] = Intercept_Valence_2;
merged_Intercept_Valence[Kthres_start_Valence[3]:Kthres_end_Valence[3]] = Intercept_Valence_3;
merged_Intercept_Arousal[Kthres_start_Arousal[1]:Kthres_end_Arousal[1]] = Intercept_Arousal_1;
merged_Intercept_Arousal[Kthres_start_Arousal[2]:Kthres_end_Arousal[2]] = Intercept_Arousal_2;
merged_Intercept_Arousal[Kthres_start_Arousal[3]:Kthres_end_Arousal[3]] = Intercept_Arousal_3;
for (n in 1:N_Valence) {
// add more terms to the linear predictor
mu_Valence[n] += r_1_Valence_1[J_1_Valence[n]] * Z_1_Valence_1[n] + r_2_Valence_1[J_2_Valence[n]] * Z_2_Valence_1[n] + r_3_Valence_1[J_3_Valence[n]] * Z_3_Valence_1[n];
}
for (n in 1:N_Arousal) {
// add more terms to the linear predictor
mu_Arousal[n] += r_1_Arousal_2[J_1_Arousal[n]] * Z_1_Arousal_2[n] + r_2_Arousal_2[J_2_Arousal[n]] * Z_2_Arousal_2[n] + r_3_Arousal_2[J_3_Arousal[n]] * Z_3_Arousal_2[n];
}
for (n in 1:N_Valence) {
// apply the inverse link function
disc_Valence[n] = exp(disc_Valence[n]);
}
for (n in 1:N_Arousal) {
// apply the inverse link function
disc_Arousal[n] = exp(disc_Arousal[n]);
}
// priors including all constants
target += student_t_lpdf(b_Valence | 6,1,1.5);
target += normal_lpdf(Intercept_Valence_2 | 0,6);
target += normal_lpdf(Intercept_Valence_3 | 0,6);
target += normal_lpdf(b_disc_Valence[1] | 0,2);
target += normal_lpdf(b_disc_Valence[2] | 0,2);
target += student_t_lpdf(b_Arousal | 6,1,1.5);
target += normal_lpdf(Intercept_Arousal_2 | 0,6);
target += normal_lpdf(Intercept_Arousal_3 | 0,6);
target += normal_lpdf(b_disc_Arousal[1] | 0,2);
target += normal_lpdf(b_disc_Arousal[2] | 0,2);
target += exponential_lpdf(sd_1[1] | 2);
target += exponential_lpdf(sd_1[2] | 2);
target += std_normal_lpdf(to_vector(z_1));
target += lkj_corr_cholesky_lpdf(L_1 | 1);
target += exponential_lpdf(sd_2[1] | 3);
target += exponential_lpdf(sd_2[2] | 3);
target += std_normal_lpdf(to_vector(z_2));
target += lkj_corr_cholesky_lpdf(L_2 | 1);
target += exponential_lpdf(sd_3[1] | 2);
target += exponential_lpdf(sd_3[2] | 2);
target += std_normal_lpdf(to_vector(z_3));
target += lkj_corr_cholesky_lpdf(L_3 | 1);
// likelihood including all constants
if (!prior_only) {
for (n in 1:N_Valence) {
target += cumulative_probit_merged_lpmf(Y_Valence[n] | mu_Valence[n], disc_Valence[n], merged_Intercept_Valence, Jthres_Valence[n]);
}
for (n in 1:N_Arousal) {
target += cumulative_probit_merged_lpmf(Y_Arousal[n] | mu_Arousal[n], disc_Arousal[n], merged_Intercept_Arousal, Jthres_Arousal[n]);
}
}
}
[details="Summary"]
This text will be hidden
[/details]
```