# Help with notation for hierarchical gam in BRMS

I have been asked to include the notation for a model in the supplement of a paper we are submitting. I have come up with my best attempt at notating what we fit in brms, but I would greatly appreciate if anyone is able to either verify that what I have is correct or suggest how I can improve what I have.

The brms syntax to fit the model is

dn15_result_PREV_KCAL<-brm(dn15_result~s(fai) + class_focal + s(prev_day_total_kcal_using_ap_low_fermentation) + s(isotope_time_collected) + (1 +  fai + prev_day_total_kcal_using_ap_low_fermentation + isotope_time_collected  | name_focal),
data = data[complete.cases(data[,c("dn15_result","isotope_time_collected","fai","prev_day_total_kcal_using_ap_low_fermentation")]),],
family=gaussian(),iter = 2000,
prior=c(prior(normal(0,10),class="Intercept"),
prior(normal(0,10),class="b"),
prior(normal(0,10),class="sigma"),
prior(normal(0,10),class="sd")),
control = list(max_treedepth = 11, adapt_delta = .999))


The stan code from this model is:

// generated with brms 2.16.3
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);
}
}
data {
int<lower=1> N;  // total number of observations
vector[N] Y;  // response variable
int<lower=1> K;  // number of population-level effects
matrix[N, K] X;  // population-level design matrix
// data for splines
int Ks;  // number of linear effects
matrix[N, Ks] Xs;  // design matrix for the linear effects
// data for spline s(fai)
int nb_1;  // number of bases
int knots_1[nb_1];  // number of knots
// basis function matrices
matrix[N, knots_1[1]] Zs_1_1;
// data for spline s(prev_day_total_kcal_using_ap_low_fermentation)
int nb_2;  // number of bases
int knots_2[nb_2];  // number of knots
// basis function matrices
matrix[N, knots_2[1]] Zs_2_1;
// data for spline s(isotope_time_collected)
int nb_3;  // number of bases
int knots_3[nb_3];  // number of knots
// basis function matrices
matrix[N, knots_3[1]] Zs_3_1;
// 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_1;
vector[N] Z_1_2;
vector[N] Z_1_3;
vector[N] Z_1_4;
int<lower=1> NC_1;  // number of group-level correlations
int prior_only;  // should the likelihood be ignored?
}
transformed data {
int Kc = K - 1;
matrix[N, Kc] Xc;  // centered version of X without an intercept
vector[Kc] means_X;  // column means of X before centering
for (i in 2:K) {
means_X[i - 1] = mean(X[, i]);
Xc[, i - 1] = X[, i] - means_X[i - 1];
}
}
parameters {
vector[Kc] b;  // population-level effects
real Intercept;  // temporary intercept for centered predictors
vector[Ks] bs;  // spline coefficients
// parameters for spline s(fai)
// standarized spline coefficients
vector[knots_1[1]] zs_1_1;
real<lower=0> sds_1_1;  // standard deviations of spline coefficients
// parameters for spline s(prev_day_total_kcal_using_ap_low_fermentation)
// standarized spline coefficients
vector[knots_2[1]] zs_2_1;
real<lower=0> sds_2_1;  // standard deviations of spline coefficients
// parameters for spline s(isotope_time_collected)
// standarized spline coefficients
vector[knots_3[1]] zs_3_1;
real<lower=0> sds_3_1;  // standard deviations of spline coefficients
real<lower=0> sigma;  // dispersion parameter
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 {
// actual spline coefficients
vector[knots_1[1]] s_1_1;
// actual spline coefficients
vector[knots_2[1]] s_2_1;
// actual spline coefficients
vector[knots_3[1]] s_3_1;
matrix[N_1, M_1] r_1;  // actual group-level effects
// using vectors speeds up indexing in loops
vector[N_1] r_1_1;
vector[N_1] r_1_2;
vector[N_1] r_1_3;
vector[N_1] r_1_4;
// compute actual spline coefficients
s_1_1 = sds_1_1 * zs_1_1;
// compute actual spline coefficients
s_2_1 = sds_2_1 * zs_2_1;
// compute actual spline coefficients
s_3_1 = sds_3_1 * zs_3_1;
// compute actual group-level effects
r_1 = scale_r_cor(z_1, sd_1, L_1);
r_1_1 = r_1[, 1];
r_1_2 = r_1[, 2];
r_1_3 = r_1[, 3];
r_1_4 = r_1[, 4];
}
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[N] mu = Intercept + rep_vector(0.0, N) + Xs * bs + Zs_1_1 * s_1_1 + Zs_2_1 * s_2_1 + Zs_3_1 * s_3_1;
for (n in 1:N) {
// add more terms to the linear predictor
mu[n] += r_1_1[J_1[n]] * Z_1_1[n] + r_1_2[J_1[n]] * Z_1_2[n] + r_1_3[J_1[n]] * Z_1_3[n] + r_1_4[J_1[n]] * Z_1_4[n];
}
target += normal_id_glm_lpdf(Y | Xc, mu, b, sigma);
}
// priors including constants
target += normal_lpdf(b | 0, 10);
target += normal_lpdf(Intercept | 0, 10);
target += normal_lpdf(bs | 0, 10)
- 3 * normal_lccdf(0 | 0, 10);
target += student_t_lpdf(sds_1_1 | 3, 0, 2.5)
- 1 * student_t_lccdf(0 | 3, 0, 2.5);
target += std_normal_lpdf(zs_1_1);
target += student_t_lpdf(sds_2_1 | 3, 0, 2.5)
- 1 * student_t_lccdf(0 | 3, 0, 2.5);
target += std_normal_lpdf(zs_2_1);
target += student_t_lpdf(sds_3_1 | 3, 0, 2.5)
- 1 * student_t_lccdf(0 | 3, 0, 2.5);
target += std_normal_lpdf(zs_3_1);
target += normal_lpdf(sigma | 0, 10)
- 1 * normal_lccdf(0 | 0, 10);
target += normal_lpdf(sd_1 | 0, 10)
- 4 * normal_lccdf(0 | 0, 10);
target += std_normal_lpdf(to_vector(z_1));
target += lkj_corr_cholesky_lpdf(L_1 | 1);
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept - dot_product(means_X, b);
// 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];
}
}
}


And this is my attempt to notate the model:

{\delta^{15}N\ }_{IJ} \sim Normal\left({mu}_{ij},sigma\right)\\ {mu}_{ij}=\alpha_0+\alpha_{ji}+s\left({Nutrient}_i,\ \tau\ \right)+\beta_{1j}{Nutrient}_i+s\left({FAI}_i,\ \tau\right)+\beta_{2j}{FAI}_i+\\s\left({sampletime\ }_i,\ \tau\right)+β3_jsampletime_i+\boldsymbol{\beta}AgeClass_i\\ \boldsymbol{\beta}{AgeClass}_i=\sum_{k=1}^{K}{\boldsymbol{\beta}_{ki}{AgeClass}_{Ki}}\\ \begin{bmatrix} \alpha_j\\β1j\\β2j\\β3j\end{bmatrix} \sim MVNormal\begin{pmatrix}\begin{bmatrix} 0\\0\\0\\0\end{bmatrix} ,S\end{pmatrix},for\ individual\ j=1,\dots J\\ S=\left(\begin{matrix}\sigma_{\alpha_j}^2&\rho_{\beta_0\beta_{1j}}&\rho_{\beta_0\beta_{2j}}&\rho_{\beta_0\beta_{3j}}\\\rho_{\beta_{1j}\alpha_j}&\sigma_{\beta_1}^2&\rho_{\beta_{1j}\beta_{2j}}&\rho_{\beta_{1j}\beta_{3j}}\\\rho_{\beta_{2j}\alpha_j}&\rho_{\beta_{2j}\beta_{1j}}&\sigma_{\beta_2}^2&\rho_{\beta_{2j}\beta_{3j}}\\\rho_{\beta_{3j}\alpha_j}&\rho_{\beta_{3j}\beta_{1j}}&\rho_{\beta_{3j}\beta_{2j}}&\sigma_{\beta_3}^2\\\end{matrix}\right)\\ \alpha_0,\alpha_j~Normal(0,10)\\ \beta_{1j},\beta_{2j},\beta_{3j},\boldsymbol{\beta} \sim Normal(0,10)\\ \sigma_{\alpha_j},\sigma_{\beta_{1j}},\sigma_{\beta_{2j}},\sigma_{\beta_{3j}},\sigma_{\boldsymbol{\beta}_j} \sim Normal(0,10)\\ \tau \sim StudentT(3,0,2.5)\\ sigma \sim Normal(0,10)\\ \boldsymbol{R} \sim {LKJ}_{corr}\left(1\right)

My questions are

1. Did I manage to get this right
and
2. Is this a sufficient minimum level of detail?

Any feedback will be much appreciated!

Seems like not only was I quite off, but there was an easier-to-read way to notate this. This is notation generated by ChatGPT based on the brms and stan code. I think it makes sense to me, but if anyone can confirm that seems fine that would be very helpful!

Let\ \textit{i} \ index \ the \ observations, \\ \textit{j} \ index \ the \ groups \ (name\_focal), \\ and \ \textit{k} \ index \ the \ spline \ basis \ functions. \\ The \ model \ can \ be \ written \ as: \\ \begin{align} y_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + X_{i}\boldsymbol{b} + Xs_{i}\boldsymbol{b_s} + Zs_{1_{i}}\boldsymbol{s_{1}} + Zs_{2_{i}}\boldsymbol{s_{2}} + Zs_{3_{i}}\boldsymbol{s_{3}} + \sum_{m=1}^{M} r_{m_{j}} Z_{m_{i}} \\ \boldsymbol{b} &\sim \text{Normal}(0, 10) \\ \alpha &\sim \text{Normal}(0, 10) \\ \boldsymbol{b_s} &\sim \text{Normal}(0, 10) \\ \boldsymbol{s_{1}} &= sds_{1_{1}} \cdot \boldsymbol{z_{s_{1}}} \\ \boldsymbol{s_{2}} &= sds_{1_{2}} \cdot \boldsymbol{z_{s_{2}}} \\ \boldsymbol{s_{3}} &= sds_{1_{3}} \cdot \boldsymbol{z_{s_{3}}} \\ sds_{1_{1}} &\sim \text{Student-t}(3, 0, 2.5) \\ sds_{1_{2}} &\sim \text{Student-t}(3, 0, 2.5) \\ sds_{1_{3}} &\sim \text{Student-t}(3, 0, 2.5) \\ \boldsymbol{z_{s_{1}}} &\sim \text{Normal}(0, 1) \\ \boldsymbol{z_{s_{2}}} &\sim \text{Normal}(0, 1) \\ \boldsymbol{z_{s_{3}}} &\sim \text{Normal}(0, 1) \\ \sigma &\sim \text{Normal}(0, 10) \\ \boldsymbol{r}_{j} &\sim \text{Multivariate Normal}(\boldsymbol{0}, \boldsymbol{\Sigma}) \\ \boldsymbol{\Sigma} &= \text{diag}(\boldsymbol{sd}) \cdot \boldsymbol{R} \cdot \text{diag}(\boldsymbol{sd}) \\ \boldsymbol{sd} &\sim \text{Normal}(0, 10) \\ \boldsymbol{R} &\sim \text{LKJ}(1) \\ \end{align} \\ where \ \alpha \ is \ the \ global \ intercept, \ \textbf{b} \ are \ the \ population-level \ effects, \ \\textbf{b_s} \ are \ the \ spline \ coefficients \ for \ the \ smooth \ functions, \\ \textbf{s_{1}}, \ \textbf{s_{2}}, \ and \ \textbf{s_{3}} \ are \ the \ actual \ spline coefficients \ for \ fai, \\ prev\_day\_total\_kcal\_using\_ap\_low\_fermentation, \\ and \ isotope\_time\_collected, \ \textbf{r}_{j} \ are \ the \ actual \ group-level \ effects, \\ \textbf{Sigma} \ is \ the \ group-level \ covariance \ matrix, \\ \textbf{sd} \ are \ the \ group-level \ standard \ deviations, \ and \\ \textbf{R} \ is \ the \ group-level \ z \ correlation \ matrix.