A toy example in cmdstanr. I had to put a pretty strong prior on sigma1 to keep it from getting stuck at zero, since in this toy example the residuals and random effects aren’t well identified otherwise.

```
data {
int<lower=1> N;
int<lower=1> J;
real y1[N];
int y2[N];
}
parameters {
// variance parameters
real<lower=0> sigma1;
real<lower=0> phi2;
// latent variables
vector[J] L[N];
cholesky_factor_corr[J] cholcor;
// scale of latent variables
vector<lower=0>[J] sigmaL;
}
model {
phi2 ~ exponential(2); // negative binomial scale
sigma1 ~ gamma(6,10); // lognormal scale
sigmaL ~ exponential(1); // scales of latent variables
cholcor ~ lkj_corr_cholesky(2); // correlations of latent variables
L ~ multi_normal_cholesky(rep_vector(0,J),cholcor); // latent variables
y1 ~ lognormal(to_vector(L[1:N,1])*sigmaL[1],sigma1); // y1 is lognormal
y2 ~ neg_binomial_2(exp(to_vector(L[1:N,2])*sigmaL[2]),phi2); // y2 is negative binomial
}
generated quantities {
real y1out[N] = lognormal_rng(to_vector(L[1:N,1])*sigmaL[1],sigma1);
int y2out[N] = neg_binomial_2_rng(exp(to_vector(L[1:N,2])*sigmaL[2]),phi2);
}
```

```
library(cmdstanr)
library(bayesplot)
library(posterior)
library(tidyverse)
# Generate some continuous and integer-valued correlated data
cormat <- matrix(c(1,0.8,0.8,1),ncol=2)
y <- MASS:: mvrnorm(50, c(0,0), cormat)
y <- exp(y)
y[,2] <- round(y[,2])
datlist <- list(
y1 = y[,1],
y2 = as.integer(y[,2]),
N = NROW(y),
J = NCOL(y)
)
model <- cmdstan_model("scripts/latent_re.stan")
fit <- model$sample(data = datlist, parallel_chains = 4, adapt_delta = 0.99, max_treedepth = 15)
```