---
title: "GLM CMCMC"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{GLM CMCMC}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE, purl=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = FALSE,
  purl = FALSE
)
```

## 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.
\]

```{r binomial-logit}
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:

```{r binomial-probit}
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.
\]

```{r poisson-log}
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)\).

```{r gaussian-identity}
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")`:

```{r 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\):

```{r gaussian-inverse}
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)\).

```{r gamma-log}
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:

```{r gamma-tail-priors}
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.

```{r output-methods}
coef(fit_binomial)
summary(fit_gaussian)
```