Linear mixed effects model with varying intercept and slope in matrix notation

Modeling
3

Thanks for the reply. I do have to work directly with Stan code for several reasons.

I have modified the code found using brms and from the following paper

Here is the code:

Xmat = model.matrix(~ 1 + length)

data.list = list(N=nrow(Xmat),Npreds = ncol(Xmat),Nlev1 = nlevels(pop),
                 n_u = ncol(Xmat),lev1 = as.numeric(pop),Xmat = Xmat,Z_u = Xmat,y = as.vector(mass))

LME_randomSlopeAndInterceptCor.stan = '
data {
  int<lower=0> N;               //no trials
  int<lower=1> Npreds;               //no fixefs
  int<lower=0> Nlev1;               //no subjects
  int<lower=1> n_u;             //no subj ranefs
   int<lower=1,upper=Nlev1> lev1[N]; //subject indicator
   row_vector[Npreds] Xmat[N];           //fixef design matrix
  row_vector[n_u] Z_u[N];       //subj ranef design matrix
  vector[N] y;                 // response
}

parameters {
  vector[Npreds] beta;               //fixef coefs
  cholesky_factor_corr[n_u] L_u;  //cholesky factor of subj ranef corr matrix
  
  vector<lower=0>[n_u] sigma_u; //subj ranef std
  real<lower=0> sigma_e;        //residual std
  vector[n_u] z_u[Nlev1];           //spherical subj ranef
}

transformed parameters {
  vector[n_u] u[Nlev1];             //subj ranefs
 
  {
    matrix[n_u,n_u] Sigma_u;    //subj ranef cov matrix
    Sigma_u = diag_pre_multiply(sigma_u,L_u);
    for(j in 1:Nlev1)
      u[j] = Sigma_u * z_u[j];
   }
}

model {
  //priors
  L_u ~ lkj_corr_cholesky(2.0);
 for (j in 1:Nlev1)
    z_u[j] ~ normal(0,1);
 
  //likelihood
  for (i in 1:N)
  	y[i] ~ normal(Xmat[i] * beta + Z_u[i] * u[lev1[i]], sigma_e);
}
'

fit_matrix = stan(file = "LME_randomSlopeAndInterceptCor.stan",data = data.list,chains = 4,
           iter = ni, warmup = nb, thin = nt)

print(fit_matrix,pars=c("beta", "sigma_e", "sigma_u"),probs = c(0.025, 0.5, 0.975))

It return the expected parameters.