GLMMs with Random Intercepts and Slopes for Various Response Distributions

Model

For observation (i=1,,n), let (g_i{1,,m}) denote its group and let the first column of (_i) be one. The Gaussian mixed model is

\[ y_i = \mathbf{x}_i^{\mathsf T}\boldsymbol\beta + u_{g_i,1}x_{i1} + u_{g_i,2}x_{i2} + \varepsilon_i, \qquad \varepsilon_i \sim \operatorname{Normal}(0,\sigma^2). \]

Equivalently,

\[ y_i \mid \boldsymbol\beta,\mathbf{u}_{g_i},\sigma \sim \operatorname{Normal}\!\left( \mathbf{x}_i^{\mathsf T}\boldsymbol\beta + \mathbf{z}_i^{\mathsf T}\mathbf{u}_{g_i}, \sigma^2 \right), \qquad \mathbf{z}_i=(1,x_{i2})^{\mathsf T}. \]

The group-specific random intercept and slope satisfy

\[ \mathbf{u}_g= \begin{pmatrix} u_{g,1}\\ u_{g,2} \end{pmatrix} \sim \operatorname{MVN}_2(\mathbf{0},\boldsymbol\Sigma_u), \qquad g=1,\ldots,m, \]

with

\[ \boldsymbol\Sigma_u = \begin{pmatrix} \tau_0^2 & \rho\tau_0\tau_1\\ \rho\tau_0\tau_1 & \tau_1^2 \end{pmatrix}, \qquad \tau_0=e^{\ell_{\tau_0}},\quad \tau_1=e^{\ell_{\tau_1}},\quad \rho=\tanh(r_{\rho}). \]

The priors used below are

\[ \beta_j\sim\operatorname{Normal}(0,10), \qquad \log\sigma\sim\operatorname{Normal}(0,2), \]

\[ \ell_{\tau_0}\sim\operatorname{Normal}(0,2), \qquad \ell_{\tau_1}\sim\operatorname{Normal}(0,2), \qquad r_{\rho}\sim\operatorname{Normal}(0,2). \]

library(hobbs)
library(coda)

set.seed(123)

n <- 100000L
p0 <- 5L
m <- 10000L

x <- matrix(rnorm(n * p0), nrow = n, ncol = p0)
x <- cbind(1, x)
p <- ncol(x)

gid <- rep(seq_len(m), each = n / m)

beta_true <- c(0.5, 1.0, -0.75, 0.5, 0.0, -0.25)

sigma_true <- 0.75

## Random intercept/slope covariance
tau0_true <- 1.25
tau1_true <- 0.60
rho_true <- -0.35

Sigma_u_true <- matrix(c(
    tau0_true^2, rho_true * tau0_true * tau1_true,
    rho_true * tau0_true * tau1_true, tau1_true^2
), nrow = 2, byrow = TRUE)

## Draw group-level random effects:
## u[,1] = random intercept
## u[,2] = random slope for x[,2]
u_true <- matrix(rnorm(m * 2), nrow = m, ncol = 2) %*% chol(Sigma_u_true)

y <- as.numeric(
    x %*% beta_true +
        rowSums(u_true[gid, , drop = FALSE] * x[, 1:2, drop = FALSE]) +
        rnorm(n, 0, sigma_true)
)

gstart <- seq.int(1L, by = n / m, length.out = m)
gend <- gstart + n / m - 1L

data <- list(
    n = n,
    p = p,
    m = m,
    x = x,
    y = y,
    gid = gid,
    gstart = gstart,
    gend = gend
)

mod = '
param beta(p) sampler=slice;
param u(m,2) save=mean;
param logsigma(1) sampler=slice;
param logtau0(1) sampler=slice;
param logtau1(1) sampler=slice;
param rho_raw(1) sampler=slice;

func llk() {
  double sigma = exp(logsigma(1));
  for (i = 1:n) {
    y(i) ~ dnorm(mu(i),sigma);
  }
}

func pllk() {
  double sigma = exp(logsigma(1));
  for (i = gstart(j):gend(j)) {
    y(i) ~ dnorm(mu(i),sigma);
  }
}

func build_sig() {
  vec zero2(2);
  mat Sigma_u(2,2);
  
  double tau0 = exp(logtau0(1));
  double tau1 = exp(logtau1(1));
  double rho = tanh(rho_raw(1));

  Sigma_u(1,1) = tau0 * tau0;
  Sigma_u(2,2) = tau1 * tau1;
  Sigma_u(1,2) = rho * tau0 * tau1;
  Sigma_u(2,1) = Sigma_u(1,2);
}

func update_u() {
  for (j = 1:m) {
    u(j,1:2) ~ dmvn(zero2,Sigma_u);
  }
}

block beta(j) {
  beta(j) ~ dnorm(0,10);
  llk();
} cache mu(n) {
  for (i = 1:n) {
    for (k = 1:p) {
      mu(i) += beta(k) * x(i,k);
    }
    for(l = 1:2) {
      mu(i) += u(gid(i),l) * x(i,l);
    }
  }
} update mu(n) {
  for (i = 1:n) {
    mu(i) += (proposal(beta(j)) - current(beta(j))) * x(i,j);
  }
}

block u(j,l) {
  build_sig();
  u(j,1:2) ~ dmvn(zero2,Sigma_u);
  pllk();
} update mu(n) {
  for (i = gstart(j):gend(j)) {
    mu(i) += (proposal(u(j,l)) - current(u(j,l))) * x(i,l);
  }
}

block logsigma(1) {
  logsigma(1) ~ dnorm(0,2);
  llk();
}

block logtau0(1) {
  logtau0(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block logtau1(1) {
  logtau1(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block rho_raw(1) {
  rho_raw(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}
'

fit = hobbs(
    model = mod,
    data = data,
    samples = 50000,
    burnin = 2000,
    out = "glmm.bin"
)

draws = read_hobbs("glmm.bin")
draws_mean = read_hobbs("glmm.mean.bin")

# -------------------------
# Posterior summaries
# -------------------------

beta_cols = paste0("beta[", seq_len(p), "]")

missing_beta = setdiff(beta_cols, names(draws))
if (length(missing_beta) > 0L) {
  stop(
    "Missing beta columns: ",
    paste(missing_beta, collapse = ", ")
  )
}

posterior_summary = function(x, truth) {
  q = quantile(
    x,
    probs = c(0.025, 0.500, 0.975),
    names = FALSE
  )

  c(
    truth = truth,
    mean = mean(x),
    q025 = q[1L],
    q500 = q[2L],
    q975 = q[3L],
    bias = mean(x) - truth,
    covered95 = as.numeric(truth >= q[1L] && truth <= q[3L])
  )
}

# Fixed effects
beta_draws = as.matrix(draws[, beta_cols, drop = FALSE])

beta_summary = do.call(
  rbind,
  lapply(seq_len(p), function(j) {
    posterior_summary(beta_draws[, j], beta_true[j])
  })
)

beta_summary = data.frame(
  parameter = paste0("beta", seq_len(p)),
  beta_summary,
  row.names = NULL,
  check.names = FALSE
)

cat("\nPosterior summaries for fixed effects:\n")
print(beta_summary, row.names = FALSE)

# Residual and random-effect covariance parameters
sigma_draws = exp(draws$`logsigma[1]`)
tau0_draws = exp(draws$`logtau0[1]`)
tau1_draws = exp(draws$`logtau1[1]`)
rho_draws = tanh(draws$`rho_raw[1]`)

covariance_draws = list(
  sigma = sigma_draws,
  tau0 = tau0_draws,
  tau1 = tau1_draws,
  rho = rho_draws
)

covariance_truth = c(
  sigma = sigma_true,
  tau0 = tau0_true,
  tau1 = tau1_true,
  rho = rho_true
)

covariance_summary = do.call(
  rbind,
  lapply(names(covariance_draws), function(nm) {
    posterior_summary(
      covariance_draws[[nm]],
      covariance_truth[[nm]]
    )
  })
)

covariance_summary = data.frame(
  parameter = names(covariance_draws),
  covariance_summary,
  row.names = NULL,
  check.names = FALSE
)

cat("\nPosterior summaries for sigma and random-effect covariance parameters:\n")
print(covariance_summary, row.names = FALSE)

cat("\n95% interval coverage of generating global parameters:\n")
print(c(
  fixed_effects = sum(beta_summary$covered95),
  fixed_effects_total = nrow(beta_summary),
  covariance_parameters = sum(covariance_summary$covered95),
  covariance_parameters_total = nrow(covariance_summary)
))

# -------------------------
# Random-effect recovery
# -------------------------

# u is declared with save=mean, so the companion mean file contains
# posterior means rather than the full m x 2 random-effect chain.
u_cols = grep("^u\\[", names(draws_mean), value = TRUE)

u_match = regexec(
  "^u\\[([0-9]+),([0-9]+)\\]$",
  u_cols
)
u_parts = regmatches(u_cols, u_match)
u_valid = lengths(u_parts) == 3L

if (sum(u_valid) == m * 2L) {
  u_post_mean = matrix(NA_real_, nrow = m, ncol = 2L)

  for (k in which(u_valid)) {
    g = as.integer(u_parts[[k]][2L])
    l = as.integer(u_parts[[k]][3L])
    u_post_mean[g, l] = as.numeric(draws_mean[[u_cols[k]]][1L])
  }

  random_effect_summary = data.frame(
    effect = c("random intercept", "random slope"),
    truth_sd = c(sd(u_true[, 1L]), sd(u_true[, 2L])),
    posterior_mean_sd = c(sd(u_post_mean[, 1L]), sd(u_post_mean[, 2L])),
    correlation = c(
      cor(u_true[, 1L], u_post_mean[, 1L]),
      cor(u_true[, 2L], u_post_mean[, 2L])
    ),
    rmse = c(
      sqrt(mean((u_post_mean[, 1L] - u_true[, 1L])^2)),
      sqrt(mean((u_post_mean[, 2L] - u_true[, 2L])^2))
    ),
    mae = c(
      mean(abs(u_post_mean[, 1L] - u_true[, 1L])),
      mean(abs(u_post_mean[, 2L] - u_true[, 2L]))
    )
  )

  cat("\nRecovery of simulated group-level random effects:\n")
  print(random_effect_summary, row.names = FALSE)
} else {
  warning(
    "Could not identify all saved random-effect means in glmm.mean.bin; ",
    "skipping random-effect recovery summaries."
  )
  u_post_mean = NULL
}

# -------------------------
# ESS summaries
# -------------------------

ess_draws = data.frame(
  beta_draws,
  sigma = sigma_draws,
  tau0 = tau0_draws,
  tau1 = tau1_draws,
  rho = rho_draws,
  check.names = FALSE
)
names(ess_draws)[seq_len(p)] = beta_cols

cat("\nEffective sample sizes for global parameters:\n")
print(effectiveSize(as.mcmc(ess_draws)))

# -------------------------
# Diagnostic plots
# -------------------------

op = par(
  mfrow = c(3, 2),
  mar = c(3, 4, 2, 1)
)

for (j in seq_len(p)) {
  plot(
    beta_draws[, j],
    type = "l",
    main = paste0("beta[", j, "] trace"),
    xlab = "Saved draw",
    ylab = paste0("beta[", j, "]")
  )
  abline(h = beta_true[j], col = 2)
}

par(op)

op = par(
  mfrow = c(2, 2),
  mar = c(3, 4, 2, 1)
)

plot(
  sigma_draws,
  type = "l",
  main = "sigma trace",
  xlab = "Saved draw",
  ylab = "sigma"
)
abline(h = sigma_true, col = 2)

plot(
  tau0_draws,
  type = "l",
  main = "tau0 trace",
  xlab = "Saved draw",
  ylab = "tau0"
)
abline(h = tau0_true, col = 2)

plot(
  tau1_draws,
  type = "l",
  main = "tau1 trace",
  xlab = "Saved draw",
  ylab = "tau1"
)
abline(h = tau1_true, col = 2)

plot(
  rho_draws,
  type = "l",
  main = "rho trace",
  xlab = "Saved draw",
  ylab = "rho"
)
abline(h = rho_true, col = 2)

par(op)

if (!is.null(u_post_mean)) {
  op = par(
    mfrow = c(1, 2),
    mar = c(4, 4, 2, 1)
  )

  plot(
    u_true[, 1L],
    u_post_mean[, 1L],
    pch = 16,
    cex = 0.35,
    xlab = "True random intercept",
    ylab = "Posterior mean",
    main = "Random-intercept recovery"
  )
  abline(a = 0, b = 1, col = 2)

  plot(
    u_true[, 2L],
    u_post_mean[, 2L],
    pch = 16,
    cex = 0.35,
    xlab = "True random slope",
    ylab = "Posterior mean",
    main = "Random-slope recovery"
  )
  abline(a = 0, b = 1, col = 2)

  par(op)
}

rm(beta_draws, ess_draws)
gc(verbose = FALSE)

unlink("glmm.bin")
unlink("glmm.bin.metadata.rds")
unlink("glmm.bin.param_names.rds")
unlink("glmm.bin.adaptation.csv")
unlink("glmm.bin.adapt_covariance.csv")

Other GLMM response distributions

The following sections keep the same random-intercept/random-slope structure and cached linear predictor,

\[ \eta_i = \mathbf{x}_i^{\mathsf T}\boldsymbol\beta + u_{g_i,1}x_{i1} + u_{g_i,2}x_{i2}. \]

Each section first shows only the code that changes from the Gaussian model and then gives the complete hobbs model string. These examples are not evaluated when the document is knitted.

Lognormal response

\[ \log y_i\sim\operatorname{Normal}(\eta_i,\sigma_{\log}^2). \]

exp(eta(i)) is the conditional median, not the arithmetic mean.

What changes

model_changes = '
param logsdlog(1);

func llk() {
  double sdlog = exp(logsdlog(1));
  for (i = 1:n) {
    y(i) ~ dlnorm(eta(i),sdlog);
  }
}

func pllk() {
  double sdlog = exp(logsdlog(1));
  for (i = gstart(j):gend(j)) {
    y(i) ~ dlnorm(eta(i),sdlog);
  }
}

block logsdlog(1) {
  logsdlog(1) ~ dnorm(0,2);
  llk();
}
'

Full model

mod = '
param beta(p);
param u(m,2) save=mean;
param logtau0(1);
param logtau1(1);
param rho_raw(1);
param logsdlog(1);

func llk() {
  double sdlog = exp(logsdlog(1));
  for (i = 1:n) {
    y(i) ~ dlnorm(eta(i),sdlog);
  }
}

func pllk() {
  double sdlog = exp(logsdlog(1));
  for (i = gstart(j):gend(j)) {
    y(i) ~ dlnorm(eta(i),sdlog);
  }
}

func build_sig() {
  vec zero2(2);
  mat Sigma_u(2,2);
  
  double tau0 = exp(logtau0(1));
  double tau1 = exp(logtau1(1));
  double rho = tanh(rho_raw(1));

  Sigma_u(1,1) = tau0 * tau0;
  Sigma_u(2,2) = tau1 * tau1;
  Sigma_u(1,2) = rho * tau0 * tau1;
  Sigma_u(2,1) = Sigma_u(1,2);
}

func update_u() {
  for (j = 1:m) {
    u(j,1:2) ~ dmvn(zero2,Sigma_u);
  }
}

block beta(j) {
  beta(j) ~ dnorm(0,10);
  llk();
} cache eta(n) {
  for (i = 1:n) {
    for (k = 1:p) {
      eta(i) += beta(k) * x(i,k);
    }
    for (l = 1:2) {
      eta(i) += u(gid(i),l) * x(i,l);
    }
  }
} update eta(n) {
  for (i = 1:n) {
    eta(i) += (proposal(beta(j)) - current(beta(j))) * x(i,j);
  }
}

block u(j,l) {
  build_sig();
  u(j,1:2) ~ dmvn(zero2,Sigma_u);
  pllk();
} update eta(n) {
  for (i = gstart(j):gend(j)) {
    eta(i) += (proposal(u(j,l)) - current(u(j,l))) * x(i,l);
  }
}

block logtau0(1) {
  logtau0(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block logtau1(1) {
  logtau1(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block rho_raw(1) {
  rho_raw(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block logsdlog(1) {
  logsdlog(1) ~ dnorm(0,2);
  llk();
}
'

Beta regression

\[ y_i\sim\operatorname{Beta}(\mu_i\phi,(1-\mu_i)\phi),\qquad \operatorname{logit}(\mu_i)=\eta_i. \]

This requires 0 < y(i) < 1.

What changes

model_changes = '
param logphi(1);

func llk() {
  double phi = exp(logphi(1));
  for (i = 1:n) {
    double mean_i = 1.0 / (1.0 + exp(-eta(i)));
    y(i) ~ dbeta(mean_i * phi,(1.0 - mean_i) * phi);
  }
}

func pllk() {
  double phi = exp(logphi(1));
  for (i = gstart(j):gend(j)) {
    double mean_i = 1.0 / (1.0 + exp(-eta(i)));
    y(i) ~ dbeta(mean_i * phi,(1.0 - mean_i) * phi);
  }
}

block logphi(1) {
  logphi(1) ~ dnorm(0,2);
  llk();
}
'

Full model

mod = '
param beta(p);
param u(m,2) save=mean;
param logtau0(1);
param logtau1(1);
param rho_raw(1);
param logphi(1);

func llk() {
  double phi = exp(logphi(1));
  for (i = 1:n) {
    double mean_i = 1.0 / (1.0 + exp(-eta(i)));
    y(i) ~ dbeta(mean_i * phi,(1.0 - mean_i) * phi);
  }
}

func pllk() {
  double phi = exp(logphi(1));
  for (i = gstart(j):gend(j)) {
    double mean_i = 1.0 / (1.0 + exp(-eta(i)));
    y(i) ~ dbeta(mean_i * phi,(1.0 - mean_i) * phi);
  }
}

func build_sig() {
  vec zero2(2);
  mat Sigma_u(2,2);
  
  double tau0 = exp(logtau0(1));
  double tau1 = exp(logtau1(1));
  double rho = tanh(rho_raw(1));

  Sigma_u(1,1) = tau0 * tau0;
  Sigma_u(2,2) = tau1 * tau1;
  Sigma_u(1,2) = rho * tau0 * tau1;
  Sigma_u(2,1) = Sigma_u(1,2);
}

func update_u() {
  for (j = 1:m) {
    u(j,1:2) ~ dmvn(zero2,Sigma_u);
  }
}

block beta(j) {
  beta(j) ~ dnorm(0,10);
  llk();
} cache eta(n) {
  for (i = 1:n) {
    for (k = 1:p) {
      eta(i) += beta(k) * x(i,k);
    }
    for (l = 1:2) {
      eta(i) += u(gid(i),l) * x(i,l);
    }
  }
} update eta(n) {
  for (i = 1:n) {
    eta(i) += (proposal(beta(j)) - current(beta(j))) * x(i,j);
  }
}

block u(j,l) {
  build_sig();
  u(j,1:2) ~ dmvn(zero2,Sigma_u);
  pllk();
} update eta(n) {
  for (i = gstart(j):gend(j)) {
    eta(i) += (proposal(u(j,l)) - current(u(j,l))) * x(i,l);
  }
}

block logtau0(1) {
  logtau0(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block logtau1(1) {
  logtau1(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block rho_raw(1) {
  rho_raw(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block logphi(1) {
  logphi(1) ~ dnorm(0,2);
  llk();
}
'

Student-t residuals

\[ y_i\sim t_{\nu}(\eta_i,\sigma). \]

What changes

model_changes = '
param logsigma(1);
param logdf(1);

func llk() {
  double sigma = exp(logsigma(1));
  double df = exp(logdf(1));
  for (i = 1:n) {
    y(i) ~ dt(df,eta(i),sigma);
  }
}

func pllk() {
  double sigma = exp(logsigma(1));
  double df = exp(logdf(1));
  for (i = gstart(j):gend(j)) {
    y(i) ~ dt(df,eta(i),sigma);
  }
}

block logsigma(1) {
  logsigma(1) ~ dnorm(0,2);
  llk();
}

block logdf(1) {
  logdf(1) ~ dnorm(1,1);
  llk();
}
'

Full model

mod = '
param beta(p);
param u(m,2) save=mean;
param logtau0(1);
param logtau1(1);
param rho_raw(1);
param logsigma(1);
param logdf(1);

func llk() {
  double sigma = exp(logsigma(1));
  double df = exp(logdf(1));
  for (i = 1:n) {
    y(i) ~ dt(df,eta(i),sigma);
  }
}

func pllk() {
  double sigma = exp(logsigma(1));
  double df = exp(logdf(1));
  for (i = gstart(j):gend(j)) {
    y(i) ~ dt(df,eta(i),sigma);
  }
}

func build_sig() {
  vec zero2(2);
  mat Sigma_u(2,2);
  
  double tau0 = exp(logtau0(1));
  double tau1 = exp(logtau1(1));
  double rho = tanh(rho_raw(1));

  Sigma_u(1,1) = tau0 * tau0;
  Sigma_u(2,2) = tau1 * tau1;
  Sigma_u(1,2) = rho * tau0 * tau1;
  Sigma_u(2,1) = Sigma_u(1,2);
}

func update_u() {
  for (j = 1:m) {
    u(j,1:2) ~ dmvn(zero2,Sigma_u);
  }
}

block beta(j) {
  beta(j) ~ dnorm(0,10);
  llk();
} cache eta(n) {
  for (i = 1:n) {
    for (k = 1:p) {
      eta(i) += beta(k) * x(i,k);
    }
    for (l = 1:2) {
      eta(i) += u(gid(i),l) * x(i,l);
    }
  }
} update eta(n) {
  for (i = 1:n) {
    eta(i) += (proposal(beta(j)) - current(beta(j))) * x(i,j);
  }
}

block u(j,l) {
  build_sig();
  u(j,1:2) ~ dmvn(zero2,Sigma_u);
  pllk();
} update eta(n) {
  for (i = gstart(j):gend(j)) {
    eta(i) += (proposal(u(j,l)) - current(u(j,l))) * x(i,l);
  }
}

block logtau0(1) {
  logtau0(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block logtau1(1) {
  logtau1(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block rho_raw(1) {
  rho_raw(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block logsigma(1) {
  logsigma(1) ~ dnorm(0,2);
  llk();
}

block logdf(1) {
  logdf(1) ~ dnorm(1,1);
  llk();
}
'

Laplace residuals

\[ y_i\sim\operatorname{Laplace}(\eta_i,b). \]

What changes

model_changes = '
param logscale(1);

func llk() {
  double scale = exp(logscale(1));
  for (i = 1:n) {
    y(i) ~ dlaplace(eta(i),scale);
  }
}

func pllk() {
  double scale = exp(logscale(1));
  for (i = gstart(j):gend(j)) {
    y(i) ~ dlaplace(eta(i),scale);
  }
}

block logscale(1) {
  logscale(1) ~ dnorm(0,2);
  llk();
}
'

Full model

mod = '
param beta(p);
param u(m,2) save=mean;
param logtau0(1);
param logtau1(1);
param rho_raw(1);
param logscale(1);

func llk() {
  double scale = exp(logscale(1));
  for (i = 1:n) {
    y(i) ~ dlaplace(eta(i),scale);
  }
}

func pllk() {
  double scale = exp(logscale(1));
  for (i = gstart(j):gend(j)) {
    y(i) ~ dlaplace(eta(i),scale);
  }
}

func build_sig() {
  vec zero2(2);
  mat Sigma_u(2,2);
  
  double tau0 = exp(logtau0(1));
  double tau1 = exp(logtau1(1));
  double rho = tanh(rho_raw(1));

  Sigma_u(1,1) = tau0 * tau0;
  Sigma_u(2,2) = tau1 * tau1;
  Sigma_u(1,2) = rho * tau0 * tau1;
  Sigma_u(2,1) = Sigma_u(1,2);
}

func update_u() {
  for (j = 1:m) {
    u(j,1:2) ~ dmvn(zero2,Sigma_u);
  }
}

block beta(j) {
  beta(j) ~ dnorm(0,10);
  llk();
} cache eta(n) {
  for (i = 1:n) {
    for (k = 1:p) {
      eta(i) += beta(k) * x(i,k);
    }
    for (l = 1:2) {
      eta(i) += u(gid(i),l) * x(i,l);
    }
  }
} update eta(n) {
  for (i = 1:n) {
    eta(i) += (proposal(beta(j)) - current(beta(j))) * x(i,j);
  }
}

block u(j,l) {
  build_sig();
  u(j,1:2) ~ dmvn(zero2,Sigma_u);
  pllk();
} update eta(n) {
  for (i = gstart(j):gend(j)) {
    eta(i) += (proposal(u(j,l)) - current(u(j,l))) * x(i,l);
  }
}

block logtau0(1) {
  logtau0(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block logtau1(1) {
  logtau1(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block rho_raw(1) {
  rho_raw(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block logscale(1) {
  logscale(1) ~ dnorm(0,2);
  llk();
}
'

Logistic residuals

\[ y_i\sim\operatorname{Logistic}(\eta_i,s). \]

What changes

model_changes = '
param logscale(1);

func llk() {
  double scale = exp(logscale(1));
  for (i = 1:n) {
    y(i) ~ dlogis(eta(i),scale);
  }
}

func pllk() {
  double scale = exp(logscale(1));
  for (i = gstart(j):gend(j)) {
    y(i) ~ dlogis(eta(i),scale);
  }
}

block logscale(1) {
  logscale(1) ~ dnorm(0,2);
  llk();
}
'

Full model

mod = '
param beta(p);
param u(m,2) save=mean;
param logtau0(1);
param logtau1(1);
param rho_raw(1);

param logscale(1);

func llk() {
  double scale = exp(logscale(1));
  for (i = 1:n) {
    y(i) ~ dlogis(eta(i),scale);
  }
}

func pllk() {
  double scale = exp(logscale(1));
  for (i = gstart(j):gend(j)) {
    y(i) ~ dlogis(eta(i),scale);
  }
}

func build_sig() {
  vec zero2(2);
  mat Sigma_u(2,2);
  
  double tau0 = exp(logtau0(1));
  double tau1 = exp(logtau1(1));
  double rho = tanh(rho_raw(1));

  Sigma_u(1,1) = tau0 * tau0;
  Sigma_u(2,2) = tau1 * tau1;
  Sigma_u(1,2) = rho * tau0 * tau1;
  Sigma_u(2,1) = Sigma_u(1,2);
}

func update_u() {
  for (j = 1:m) {
    u(j,1:2) ~ dmvn(zero2,Sigma_u);
  }
}

block beta(j) {
  beta(j) ~ dnorm(0,10);
  llk();
} cache eta(n) {
  for (i = 1:n) {
    for (k = 1:p) {
      eta(i) += beta(k) * x(i,k);
    }
    for (l = 1:2) {
      eta(i) += u(gid(i),l) * x(i,l);
    }
  }
} update eta(n) {
  for (i = 1:n) {
    eta(i) += (proposal(beta(j)) - current(beta(j))) * x(i,j);
  }
}

block u(j,l) {
  build_sig();
  u(j,1:2) ~ dmvn(zero2,Sigma_u);
  pllk();
} update eta(n) {
  for (i = gstart(j):gend(j)) {
    eta(i) += (proposal(u(j,l)) - current(u(j,l))) * x(i,l);
  }
}

block logtau0(1) {
  logtau0(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block logtau1(1) {
  logtau1(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block rho_raw(1) {
  rho_raw(1) ~ dnorm(0,2);
  build_sig();
  update_u();
}

block logscale(1) {
  logscale(1) ~ dnorm(0,2);
  llk();
}
'

Other available distributions

The same construction can be adapted to probability-scale Bernoulli and binomial models (dbern, dbinom), direct-rate Poisson models (dpois), size-probability negative-binomial models (dnbinom), and specialized positive or bounded responses using dunif, dinvgamma, dcauchy, dchisq, dpareto, dhalfnorm, or dhalfcauchy.

The multivariate distributions dbvn, dmvn, dwish, dinvwish, and dlkjcorr2 are normally used for multivariate responses or covariance priors, rather than as direct replacements for the scalar observation likelihood.