GLM CMCMC

Overview

glm_cmcmc() is the formula interface for the built-in CMCMC GLM kernels. It uses the same model formula style as stats::glm(), builds the model frame and design matrix, checks the response, chooses the correct built-in GLM target kernel, and then calls cmcmc() internally.

The supported family-link combinations are:

Family	Supported links	Sampled parameters
Binomial	`logit`, `probit`	regression coefficients
Poisson	`log`	regression coefficients
Gaussian	`identity`, `log`, `inverse`	regression coefficients and `sigma2`
Gamma	`log`	regression coefficients and `alpha`

All GLM models use independent normal priors on the coefficients,

\[ \beta_j \sim N(0, \sigma_\beta^2), \]

where prior_var is \(\sigma_\beta^2\). Gaussian and Gamma models also sample one auxiliary scale/shape parameter, controlled by prior.override.

glm_cmcmc() uses backend = "auto" through cmcmc(): CUDA is used when it is available, otherwise the CPU backend is used. Set verbose = TRUE to print which backend was selected.

Binomial

Use the binomial family for binary responses. Numeric responses must contain only 0 and 1; two-level factors and logical responses are also accepted. The supported links are logit and probit.

For the logit link,

\[ y_i \mid \beta \sim \operatorname{Bernoulli}(p_i), \qquad \operatorname{logit}(p_i) = x_i^\top \beta. \]

library(CMCMC)

set.seed(1)
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)

eta <- -0.3 + 0.8 * x1 - 0.5 * x2
p <- 1 / (1 + exp(-eta))
y <- rbinom(n, size = 1, prob = p)

dat <- data.frame(y = y, x1 = x1, x2 = x2)

fit_binomial <- CMCMC::glm_cmcmc(
  y ~ x1 + x2,
  data = dat,
  family = binomial("logit"),
  k = 128,
  it = 100,
  covit = 10,
  prior_var = 100,
  seed = 1,
  verbose = TRUE
)

fit_binomial
summary(fit_binomial)
coef(fit_binomial)

For a probit model, change only the family:

fit_binomial_probit <- CMCMC::glm_cmcmc(
  y ~ x1 + x2,
  data = dat,
  family = binomial("probit"),
  k = 128,
  it = 100,
  covit = 10,
  seed = 1
)

Poisson

Use the Poisson family for non-negative integer counts. The supported link is log.

\[ y_i \mid \beta \sim \operatorname{Poisson}(\mu_i), \qquad \log(\mu_i) = x_i^\top \beta. \]

library(CMCMC)

set.seed(2)
n <- 120
x1 <- rnorm(n)
x2 <- rnorm(n)

eta <- 0.2 + 0.4 * x1 - 0.3 * x2
mu <- exp(eta)
y <- rpois(n, lambda = mu)

dat <- data.frame(y = y, x1 = x1, x2 = x2)

fit_poisson <- CMCMC::glm_cmcmc(
  y ~ x1 + x2,
  data = dat,
  family = poisson("log"),
  k = 128,
  it = 100,
  covit = 10,
  prior_var = 100,
  seed = 2,
  verbose = TRUE
)

fit_poisson
summary(fit_poisson)
coef(fit_poisson)

Gaussian

Use the Gaussian family for finite numeric responses. The supported links are identity, log, and inverse.

For the identity link,

\[ y_i \mid \beta,\sigma \sim N(\mu_i,\sigma^2), \qquad \mu_i = x_i^\top \beta. \]

The sampled auxiliary parameter is sigma2. The supported residual-scale priors are:

`sigma_prior`	Prior	Main fields
`"uniform_sigma2"`	bounded uniform prior on \(\sigma^2\)	`sigma2_lower`, `sigma2_upper`, `sigma2_init`
`"uniform_sigma"`	bounded uniform prior on \(\sigma\)	`sigma_lower`, `sigma_upper`, `sigma2_init`
`"exponential_sigma"`	exponential prior on \(\sigma\), calibrated by \(P(\sigma > U)=q\)	`sigma_upper`, `sigma_tail_prob`, `sigma2_init`
`"half_cauchy_sigma"`	half-Cauchy prior on \(\sigma\), calibrated by \(P(\sigma > U)=q\)	`sigma_upper`, `sigma_tail_prob`, `sigma2_init`

For the tail-calibrated priors, sigma_upper is \(U\) and sigma_tail_prob is \(q\). The exponential prior uses \(\lambda = -\log(q) / U\). The half-Cauchy prior uses scale \(A = U \tan(\pi q / 2)\).

library(CMCMC)

set.seed(3)
n <- 120
x1 <- rnorm(n)
x2 <- rnorm(n)

mu <- 0.5 + 1.2 * x1 - 0.7 * x2
sigma <- 0.4
y <- mu + rnorm(n, sd = sigma)

dat <- data.frame(y = y, x1 = x1, x2 = x2)

fit_gaussian <- CMCMC::glm_cmcmc(
  y ~ x1 + x2,
  data = dat,
  family = gaussian("identity"),
  k = 128,
  it = 100,
  covit = 10,
  prior_var = 100,
  prior.override = list(
    sigma_prior = "uniform_sigma",
    sigma_upper = 1
  ),
  seed = 3,
  verbose = TRUE
)

fit_gaussian
summary(fit_gaussian)
coef(fit_gaussian)

For a Gaussian log-link model, generate positive means and use family = gaussian("log"):

set.seed(4)
n <- 120
x1 <- rnorm(n)
x2 <- rnorm(n)

eta <- 0.2 + 0.15 * x1 - 0.1 * x2
mu <- exp(eta)
y <- mu + rnorm(n, sd = 0.2)

dat <- data.frame(y = y, x1 = x1, x2 = x2)

fit_gaussian_log <- CMCMC::glm_cmcmc(
  y ~ x1 + x2,
  data = dat,
  family = gaussian("log"),
  k = 128,
  it = 100,
  covit = 10,
  prior.override = list(
    sigma_prior = "exponential_sigma",
    sigma_upper = 1,
    sigma_tail_prob = 0.05
  ),
  seed = 4
)

For the inverse link, the linear predictor must stay away from zero because \(\mu_i = 1 / \eta_i\):

set.seed(5)
n <- 120
x1 <- runif(n, -0.5, 0.5)
x2 <- runif(n, -0.5, 0.5)

eta <- 1.4 + 0.2 * x1 - 0.1 * x2
mu <- 1 / eta
y <- mu + rnorm(n, sd = 0.05)

dat <- data.frame(y = y, x1 = x1, x2 = x2)

fit_gaussian_inverse <- CMCMC::glm_cmcmc(
  y ~ x1 + x2,
  data = dat,
  family = gaussian("inverse"),
  k = 128,
  it = 100,
  covit = 10,
  prior.override = list(
    sigma_prior = "uniform_sigma",
    sigma_upper = 0.3
  ),
  seed = 5
)

Gamma

Use the Gamma family for positive continuous responses. The supported link is log.

CMCMC uses the Gamma GLM parameterisation

\[ y_i \mid \beta,\alpha \sim \operatorname{Gamma}\left( \operatorname{shape}=\alpha, \operatorname{rate}=\frac{\alpha}{\mu_i} \right), \qquad \log(\mu_i) = x_i^\top \beta. \]

This gives

\[ E(y_i) = \mu_i, \qquad \operatorname{Var}(y_i) = \frac{\mu_i^2}{\alpha}, \qquad CV = \frac{1}{\sqrt{\alpha}}. \]

The sampled auxiliary parameter is alpha, but the prior is specified on \(CV\). The supported CV priors are:

`cv_prior`	Prior	Main fields
`"uniform_cv"`	bounded uniform prior on \(CV\)	`cv_lower`, `cv_upper`, `alpha_init`
`"exponential_cv"`	exponential prior on \(CV\), calibrated by \(P(CV > U)=q\)	`cv_upper`, `cv_tail_prob`, `alpha_init`
`"half_cauchy_cv"`	half-Cauchy prior on \(CV\), calibrated by \(P(CV > U)=q\)	`cv_upper`, `cv_tail_prob`, `alpha_init`

For the tail-calibrated CV priors, cv_upper is \(U\) and cv_tail_prob is \(q\). The exponential prior uses \(\lambda = -\log(q) / U\). The half-Cauchy prior uses scale \(A = U \tan(\pi q / 2)\).

library(CMCMC)

set.seed(6)
n <- 120
x1 <- rnorm(n)
x2 <- rnorm(n)

eta <- 0.3 + 0.4 * x1 - 0.2 * x2
mu <- exp(eta)
alpha <- 8
y <- rgamma(n, shape = alpha, rate = alpha / mu)

dat <- data.frame(y = y, x1 = x1, x2 = x2)

fit_gamma <- CMCMC::glm_cmcmc(
  y ~ x1 + x2,
  data = dat,
  family = Gamma("log"),
  k = 128,
  it = 100,
  covit = 10,
  prior_var = 100,
  prior.override = list(
    cv_prior = "uniform_cv",
    cv_upper = 1
  ),
  seed = 6,
  verbose = TRUE
)

fit_gamma
summary(fit_gamma)
coef(fit_gamma)

For example:

fit_gamma_pc <- CMCMC::glm_cmcmc(
  y ~ x1 + x2,
  data = dat,
  family = Gamma("log"),
  k = 128,
  it = 100,
  covit = 10,
  prior.override = list(
    cv_prior = "exponential_cv",
    cv_upper = 0.75,
    cv_tail_prob = 0.05
  ),
  seed = 6
)

fit_gamma_half_cauchy <- CMCMC::glm_cmcmc(
  y ~ x1 + x2,
  data = dat,
  family = Gamma("log"),
  k = 128,
  it = 100,
  covit = 10,
  prior.override = list(
    cv_prior = "half_cauchy_cv",
    cv_upper = 0.75,
    cv_tail_prob = 0.05
  ),
  seed = 6
)

Output

glm_cmcmc() returns a data frame with class "glm_cmcmc". The first two columns are iter and particle; the remaining columns are posterior draws. coef() returns posterior means for regression coefficients only, while summary() reports marginal summaries for all sampled parameters.

coef(fit_binomial)
summary(fit_gaussian)