---
title: "Chapter 02-S04: Gamma–Poisson Conjugacy for One Count Rate"
author: "Kjell Nygren"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
bibliography: REFERENCES.bib
reference-section-title: References
vignette: >
  %\VignetteIndexEntry{Chapter 02-S04: Gamma–Poisson Conjugacy for One Count Rate}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r setup}
library(glmbayes)
if (requireNamespace("bayesrules", quietly = TRUE)) library(bayesrules)
```

---

## 1. The model

### 1.1 Likelihood

Suppose we observe \(n\) independent counts \(y_1, \dots, y_n\) drawn from a Poisson distribution with unknown rate \(\lambda > 0\):

\[
  Y_i \mid \lambda \;\sim\; \mathrm{Poisson}(\lambda), \qquad i = 1,\dots,n.
\]

The joint likelihood is

\[
  p(y_1,\dots,y_n \mid \lambda)
  \;\propto\; e^{-n\lambda}\,\lambda^{\sum_i y_i}.
\]

Everything relevant for updating beliefs about \(\lambda\) is captured by two sufficient statistics: \(n\) (number of observations) and \(\sum_i y_i\) (total observed count).

**With exposure.**  In many applications the counts are observed over different exposure periods \(e_i\) (time, area, number of opportunities).  The model becomes \(Y_i \mid \lambda \sim \mathrm{Poisson}(e_i \lambda)\), and the sufficient statistics become \(\sum_i e_i\) and \(\sum_i y_i\).  An exposure-weighted example appears in **Appendix A** ([@Albert2009]).

### 1.2 Prior

The rate \(\lambda > 0\) calls for a positive prior.  The **Gamma distribution** in shape–rate form is the natural conjugate choice:

\[
  p(\lambda) \;=\; \frac{\beta^\alpha}{\Gamma(\alpha)}\,\lambda^{\alpha-1}\,e^{-\beta\lambda},
  \qquad \lambda > 0,
\]

with shape \(\alpha > 0\) and rate \(\beta > 0\).  Its mean is \(\alpha/\beta\) and its variance is \(\alpha/\beta^2\).

**Interpreting the parameters.**  Think of \(\alpha - 1\) as a "prior total count" and \(\beta\) as a "prior number of observations".  The prior mean \(\alpha/\beta\) is your pre-data estimate of the rate, and larger \(\beta\) encodes a stronger commitment to that estimate.

Note: many textbooks use **shape–scale** form (\(\theta = 1/\beta\)); **glmbayes** uses the **shape–rate** convention throughout.

### 1.3 Posterior derivation

Multiply the likelihood kernel \(e^{-n\lambda}\lambda^{\sum y_i}\) by the prior kernel \(\lambda^{\alpha-1}e^{-\beta\lambda}\):

\[
  p(\lambda \mid y)
  \;\propto\;
  \lambda^{(\alpha + \sum y_i) - 1}\,e^{-(\beta + n)\lambda}.
\]

This is a Gamma kernel, so the conjugate posterior is:

\[
  \boxed{\lambda \mid y \;\sim\; \mathrm{Gamma}\!\left(\alpha + \textstyle\sum_i y_i,\;\; \beta + n\right).}
\]

**Update rule:** add the total count to the shape, add the sample size (or total exposure) to the rate.

The posterior mean is

\[
  \mathbb{E}[\lambda \mid y]
  \;=\; \frac{\alpha + \sum y_i}{\beta + n}
  \;=\;
  \frac{\beta}{\beta + n}\cdot\frac{\alpha}{\beta}
  \;+\;
  \frac{n}{\beta + n}\cdot\bar{y},
\]

again a weighted average of the prior mean and the data mean.

---

## 2. `bayesrules` illustration

The **`bayesrules`** package [@JohnsonLauEtAl2022] provides `plot_gamma_poisson()` and `summarize_gamma_poisson()` for visualizing the Gamma–Poisson update.  The same scenario appears in the Bayes Rules! conjugate vignette and in Chapter 12 of the book.

**Scenario.**  You observe daily bicycle rental counts \(y_1, \dots, y_9\) from the first nine days of a year.  The data are \(y = (1, 1, 1, 0, 0, 0, 0, 0, 0)\), giving \(\sum y = 3\) across \(n = 9\) days.  Your prior belief is a **\(\Gamma(3, 4)\)** distribution on the daily rate \(\lambda\) (prior mean = 0.75 rentals/day).

```{r gp-br-data}
br_y_df  <- data.frame(y = c(1, 1, 1, rep(0, 6)))  ## n = 9, sum(y) = 3
br_n     <- nrow(br_y_df)
br_shape <- 3
br_rate  <- 4

## Analytic posterior
post_shape_br <- br_shape + sum(br_y_df$y)   ## 3 + 3 = 6
post_rate_br  <- br_rate  + br_n              ## 4 + 9 = 13
```

```{r gp-br-plot, eval = requireNamespace("bayesrules", quietly = TRUE)}
library(bayesrules)
plot_gamma_poisson(
  shape      = br_shape,
  rate       = br_rate,
  sum_y      = sum(br_y_df$y),
  n          = br_n,
  prior      = TRUE,
  likelihood = TRUE,
  posterior  = TRUE
)
```

```{r gp-br-summary, eval = requireNamespace("bayesrules", quietly = TRUE)}
summarize_gamma_poisson(shape = br_shape, rate = br_rate,
                        sum_y = sum(br_y_df$y), n = br_n)
```

**Reading the output.**  The posterior is \(\Gamma(6, 13)\) with mean \(6/13 \approx 0.46\).  The data (mean \(\bar y = 0.33\)) have pulled the posterior mean down from the prior mean (0.75), and the posterior is noticeably tighter.

```{r gp-analytic}
gp_analytic <- data.frame(
  Example = "Bayes Rules! daily counts",
  n = br_n,
  `sum(y)` = sum(br_y_df$y),
  Posterior = sprintf("Gamma(%d, %d)", post_shape_br, post_rate_br),
  Mean = post_shape_br / post_rate_br,
  SD = sqrt(post_shape_br) / post_rate_br,
  check.names = FALSE
)
knitr::kable(gp_analytic, digits = 4,
  caption = "Conjugate Gamma--Poisson posterior (Bayes Rules! scenario)")
```

---

## 3. Fitting with **`glmb()`**

The scalar Gamma–Poisson model maps to **`glmb()`** with `family = poisson(link = "identity")` and `dGamma(shape, rate, beta, Inv_Dispersion = FALSE)`.  The single coefficient is the rate \(\lambda\).

```{r gp-glmb}
gp_beta <- matrix(br_shape / br_rate, 1L, 1L, dimnames = list(NULL, "(Intercept)"))
gp_pf   <- dGamma(shape = br_shape, rate = br_rate,
                  beta = gp_beta, Inv_Dispersion = FALSE)

set.seed(2026)
fit_gp <- glmb(
  n       = 20000,
  y ~ 1,
  data    = br_y_df,
  weights = rep(1L, br_n),
  family  = poisson(link = "identity"),
  pfamily = gp_pf
)
print(fit_gp)
```

```{r gp-compare}
gp_compare <- data.frame(
  Example = "Bayes Rules! daily counts",
  Posterior = gp_analytic$Posterior,
  `Analytic Mean` = gp_analytic$Mean,
  `Analytic SD`   = gp_analytic$SD,
  `glmb Post.Mean` = fit_gp$coef.means["(Intercept)"],
  `glmb Post.Sd`   = sd(fit_gp$coefficients[, "(Intercept)", drop = TRUE]),
  check.names = FALSE
)
knitr::kable(gp_compare, digits = 4,
  caption = "Analytic vs. glmb() posterior mean and SD")
```

The **`glmb Post.Mean`** and **`glmb Post.Sd`** columns should match the analytic **Mean** and **SD** to Monte Carlo error because `dGamma(Inv_Dispersion = FALSE)` implements the exact conjugate update.

---

## 4. Connection to glmbayes

The scalar Gamma–Poisson model is the conjugate prototype for **Poisson `glmb()` with `family = poisson(link = "identity")`**.  In that setting the single coefficient is the rate \(\lambda\), the prior is `dGamma(shape, rate, beta, Inv_Dispersion = FALSE)`, and the posterior update is the closed form above.

For general Poisson **regression** with covariates (log link), conjugacy is lost and `glmb()` uses envelope-based accept–reject sampling with a `dNormal()` prior on the regression coefficients (Chapter 10).  The Bayes Rules! **equality_index** Poisson regression example is worked in **Chapter 10, Appendix A**.

---

# Appendix A. Heart transplant mortality (@Albert2009)

The **`hearttransplants`** dataset [@Albert2020LearnBayes; @ChristiansenMorris1995] records, for each of 94 U.S. hospitals, the **exposure** `e` (expected deaths under the national baseline rate) and **observed deaths** `y` within 30 days of surgery.  The model is \(y_i \mid \lambda_i \sim \mathrm{Poisson}(e_i \lambda_i)\) where \(\lambda_i\) is the hospital's standardized mortality ratio.

```{r ht-load, eval = requireNamespace("LearnBayes", quietly = TRUE)}
library(LearnBayes)
data("hearttransplants")
cat(sprintf("%d hospitals  |  sum(y) = %d  |  sum(e) = %.0f\n",
            nrow(hearttransplants), sum(hearttransplants$y), sum(hearttransplants$e)))
```

**Albert's prior.**  Albert (Ch. 3.2) uses a \(\Gamma(16, 15174)\) prior on \(\lambda\), derived from historical national data (prior mean \(\approx 0.00105\)).  He focuses on two specific hospitals:

| Hospital | Exposure \(e\) | Deaths \(y\) | Analytic posterior | Posterior mean | Posterior SD |
|----------|-----------------|-------------|-------------------|----------------|--------------|
| A | 66 | 1 | \(\Gamma(17,\; 15240)\) | \(17/15240 \approx 0.001115\) | \(\sqrt{17}/15240 \approx 0.000271\) |
| B | 1767 | 4 | \(\Gamma(20,\; 16941)\) | \(20/16941 \approx 0.001181\) | \(\sqrt{20}/16941 \approx 0.000264\) |

```{r ht-prior}
alpha0_ht <- 16L;  beta0_ht <- 15174L
ex_A <- 66L;   yobs_A <- 1L
ex_B <- 1767L; yobs_B <- 4L
```

```{r ht-albert-analytic}
ht_albert <- data.frame(
  Hospital = c("A", "B"),
  e = c(ex_A, ex_B),
  y = c(yobs_A, yobs_B),
  Posterior = c("Gamma(17, 15240)", "Gamma(20, 16941)"),
  Mean = c(
    (alpha0_ht + yobs_A) / (beta0_ht + ex_A),
    (alpha0_ht + yobs_B) / (beta0_ht + ex_B)
  ),
  SD = c(
    sqrt(alpha0_ht + yobs_A) / (beta0_ht + ex_A),
    sqrt(alpha0_ht + yobs_B) / (beta0_ht + ex_B)
  ),
  check.names = FALSE
)
knitr::kable(ht_albert, digits = 6,
  caption = "Albert Ch. 3.2: analytic posterior mean and SD (shape--rate Gamma)")
```

### A.1 Fitting with **`glmb()`**

The exposure enters as `weights`.  The conjugate sampler updates \(\Gamma(\alpha_0 + y,\; \beta_0 + e)\) exactly.

```{r ht-glmb-A, eval = requireNamespace("LearnBayes", quietly = TRUE)}
ht_beta <- matrix(alpha0_ht / beta0_ht, 1L, 1L, dimnames = list(NULL, "(Intercept)"))
ht_pf   <- dGamma(shape = alpha0_ht, rate = beta0_ht,
                  beta = ht_beta, Inv_Dispersion = FALSE)

set.seed(2026)
fit_A <- glmb(n = 20000, y ~ 1, data = data.frame(y = yobs_A),
              weights = ex_A, family = poisson(link = "identity"), pfamily = ht_pf)
print(fit_A)
```

```{r ht-glmb-B, eval = requireNamespace("LearnBayes", quietly = TRUE)}
set.seed(2026)
fit_B <- glmb(n = 20000, y ~ 1, data = data.frame(y = yobs_B),
              weights = ex_B, family = poisson(link = "identity"), pfamily = ht_pf)
print(fit_B)
```

```{r ht-glmb-compare, eval = requireNamespace("LearnBayes", quietly = TRUE)}
ht_compare <- data.frame(
  Hospital = c("A", "B"),
  e = c(ex_A, ex_B),
  y = c(yobs_A, yobs_B),
  Posterior = c("Gamma(17, 15240)", "Gamma(20, 16941)"),
  `Albert Mean` = ht_albert$Mean,
  `Albert SD`   = ht_albert$SD,
  `glmb Post.Mean` = c(fit_A$coef.means["(Intercept)"],
                       fit_B$coef.means["(Intercept)"]),
  `glmb Post.Sd`   = c(sd(fit_A$coefficients[, "(Intercept)", drop = TRUE]),
                       sd(fit_B$coefficients[, "(Intercept)", drop = TRUE])),
  check.names = FALSE
)
knitr::kable(ht_compare, digits = 6,
  caption = "Albert analytic vs. glmb() posterior mean and SD")
```

Each fit is for **one hospital**: `(Intercept)` is the posterior draw for that hospital's rate \(\lambda\).  The **`glmb Post.Mean`** and **`glmb Post.Sd`** columns should match Albert's **Mean** and **SD** to Monte Carlo error.

---

## See also

- **Chapter-02-S01** — introduction to Bayesian inference and conjugacy.
- **Chapter-02-S02** — Normal–Normal conjugacy for one mean.
- **Chapter-02-S03** — Beta–Binomial conjugacy for one proportion.
- **Chapter-02-S05** — Gamma–Gamma conjugacy for a Gamma response rate.
- **Chapter 10, Appendix A** — Bayes Rules! Poisson regression (`equality_index`).
- **Chapter 18** — hierarchical Poisson mixed model (BikeSharing; block Gibbs with `rglmb`).
