The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
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.
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:
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)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
)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
)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.
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.