---
title: "Chapter 09: Models for the Gamma Family"
output: rmarkdown::html_vignette
bibliography: REFERENCES.bib
reference-section-title: References
vignette: >
  %\VignetteIndexEntry{Chapter 09: Models for the Gamma Family}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r setup}
library(glmbayes)
```

## 1. Introductory Discussion

Gamma regression models are used when the response variable is strictly positive and right‑skewed.  
Typical applications include:

- insurance claim severities  
- cost or expenditure data  
- reaction times and waiting times  
- positive continuous measurements with heteroskedasticity

Gamma regression is a standard GLM for positive, right‑skewed responses [@NelderWedderburn1972; @McCullagh1989; @Agresti2015].

In classical statistics, these models are fit using

glm(..., family = Gamma(link = ...)).

In glmbayes, the Bayesian analogue is

glmb(..., family = Gamma(link = ...), pfamily = dNormal(...))

or, when dispersion is unknown and modeled,

glmb(..., family = Gamma(link = ...), pfamily = dNormal_Gamma(...))

or related pfamilies.

This chapter:

1. reviews the Gamma GLM structure and link functions  
2. shows how to fit classical and Bayesian Gamma regression models  
3. uses a realistic insurance example (carinsca)  
4. demonstrates how to estimate dispersion using gamma.dispersion()  
5. compares classical and Bayesian summaries, including DIC


## 2. Gamma Likelihood and Model Structure (Log Link, Fixed Dispersion)

Gamma regression models strictly positive, right-skewed responses.  
In a Gamma GLM with fixed dispersion, the response satisfies

\[
Y_i > 0, 
\qquad 
Y_i \sim \text{Gamma}(\text{mean}=\mu_i,\ \text{dispersion}=\phi),
\]

where the variance is

\[
\mathrm{Var}(Y_i) = \phi\,\mu_i^2.
\]

The mean is linked to a linear predictor through

\[
\eta_i = x_i^\top \beta,
\qquad
\mu_i = g^{-1}(\eta_i).
\]

In this chapter we focus exclusively on the **log link**, which is the standard choice for Gamma GLMs and the only link supported by `glmb()` when dispersion is fixed.


### 2.1 Weighted Gamma Log-Likelihood (Log Link)

Let \(w_i\) denote observation weights (e.g., claim counts or exposures).  
Under the log link,

\[
\mu_i = e^{\eta_i},
\]

and the Gamma log-likelihood (up to constants) is

\[
\ell(\beta)
=
\sum_{i=1}^n
w_i\left[
-\frac{1}{\phi}\,\eta_i
-
\frac{1}{\phi}\,y_i e^{-\eta_i}
\right].
\]

Terms not involving \(\beta\) have been absorbed into the constant.  
This form is used by both `glm()` and the Bayesian functions `glmb()` and `rglmb()` when dispersion is fixed.


### 2.2 Exponential-Family Representation

The Gamma distribution belongs to the exponential family [@McCullagh1989; @Agresti2015] with canonical parameter

\[
\theta_i = -\frac{1}{\mu_i},
\]

and cumulant function

\[
b(\theta_i) = -\log(-\theta_i).
\]

The variance is

\[
\mathrm{Var}(Y_i) = \phi\,\mu_i^2.
\]

Under the **log link**, the canonical parameter is a monotone transformation of the linear predictor, and the resulting log-likelihood is **log-concave in \(\eta_i\)**.  
This property is essential for stable envelope construction and accept–reject sampling in *glmbayes*.


### 2.3 Log Link and Its Properties

The log link is defined by

\[
g(\mu_i) = \log(\mu_i),
\qquad
\mu_i = e^{\eta_i}.
\]

It is the preferred link for Gamma GLMs because:

- it automatically enforces \(\mu_i > 0\),
- it yields multiplicative interpretations for coefficients,
- it produces a log-concave likelihood in \(\eta_i\),
- it is numerically stable for skewed data,
- it is the only link supported by `glmb()` when dispersion is fixed.

Other links (inverse, identity) do not preserve log-concavity and are therefore not used in this chapter.  
Models with **estimated dispersion** or **non-log links** are covered in Chapter 11.


### 2.4 Likelihood, Prior, and Posterior (Normal Prior on \(\beta\))

With the log link, the weighted log-likelihood is

\[
\ell(\beta)
=
\sum_{i=1}^n
w_i\left[
-\frac{1}{\phi}\,\eta_i
-
\frac{1}{\phi}\,y_i e^{-\eta_i}
\right].
\]

A multivariate normal prior is specified as

\[
\log p(\beta)
=
-\tfrac{1}{2}(\beta - \mu_0)^\top \Sigma_0^{-1}(\beta - \mu_0)
+ \text{const}.
\]

Combining likelihood and prior yields the log-posterior:

\[
\log p(\beta \mid y)
=
\sum_{i=1}^n
w_i\left[
-\frac{1}{\phi}\,\eta_i
-
\frac{1}{\phi}\,y_i e^{-\eta_i}
\right]
-
\tfrac{1}{2}(\beta - \mu_0)^\top \Sigma_0^{-1}(\beta - \mu_0)
+ \text{const}.
\]

Because both terms are concave in \(\beta\), the posterior is **log-concave**, enabling efficient iid sampling via the envelope-based accept–reject sampler implemented in *glmbayes*.
## 3. Link Functions for Gamma GLMs

The link function relates the mean \(\mu\) to the linear predictor \(\eta = X\beta\):

\[
g(\mu) = \eta.
\]

The Gamma family in base R supports:

| Link      | Formula              | Notes                                   |
|-----------|----------------------|-----------------------------------------|
| inverse   | \(\eta = 1/\mu\)     | canonical link                          |
| identity  | \(\eta = \mu\)       | must ensure \(\mu > 0\)                 |
| log       | \(\eta = \log(\mu)\) | most common; ensures \(\mu > 0\)        |

In practice, the log link is usually preferred because:

- it enforces \(\mu > 0\) automatically  
- it gives multiplicative interpretations for coefficients  
- it tends to be numerically stable  
- it preserves log‑concavity of the likelihood

At present (due to the need for log-concavity), only the log link is implemented in the **glmb** function.

## 4. Gamma Regression Example: Insurance Claims

We now consider a Gamma regression example based on the carinsca dataset.  
The goal is to model average claim cost as a function of rating variables Merit and Class, using claim counts as weights.

### 4.1 Data Setup

We begin by preparing the data and setting appropriate factor levels and contrasts.


```{r carinsca}

    data(carinsca)
    carinsca$Merit <- ordered(carinsca$Merit)
    carinsca$Class <- factor(carinsca$Class)
    oldopt <- options(contrasts = c("contr.treatment", "contr.treatment"))

    Claims  <- carinsca$Claims
    Insured <- carinsca$Insured
    Merit   <- carinsca$Merit
    Class   <- carinsca$Class
    Cost    <- carinsca$Cost
```  
    
The response Cost/Claims represents the **average claim cost** (severity) per claim.  
The weights Claims reflect the number of claims on which each average is based, improving efficiency and aligning with standard GLM practice for rate or average models.


### 4.2 Classical Gamma Regression Using glm()

We fit a classical Gamma GLM with a log link. As `glm()` only crudely estimates the dispersion, we use `MASS::gamma.dispersion()` [@VenablesRipley2002] and pass the estimate to `summary()`.

```{r carinsca_classical}


    out <- glm(Cost/Claims ~ Merit + Class,
               family  = Gamma(link = "log"),
               weights = Claims,
               x       = TRUE)

    ## Estimate the dispersion using MLE
    disp <- gamma.dispersion(out)
    
    summary(out,dispersion=disp)

```  

The summary shows:

- regression coefficients on the log mean‑cost scale  
- standard errors and z‑statistics  
- deviance and AIC  
- an estimated dispersion parameter (phi)

Under the log link, exp(\(\beta_{j}\)) can be interpreted as a multiplicative factor on the expected cost ratio for a one‑unit change in the corresponding covariate (or a change in factor level).


## 5. Bayesian Gamma Regression (fixed dispersion) with glmb()

To fit the Bayesian analogue, we proceed in two steps:

1. estimate the dispersion parameter \( \phi \) from the classical Gamma GLM,  
2. construct a prior for the regression coefficients using \texttt{Prior\_Setup()},  
3. call \texttt{glmb()} with a \texttt{dNormal} prior that incorporates the fixed dispersion \( \phi \).

### 5.1 Estimating Dispersion via gamma.dispersion()

The gamma.dispersion() helper estimates the dispersion parameter from a fitted classical Gamma GLM:

```{r carinsca_disp}

    ## better than the crude estimate from the summary function
    disp <- gamma.dispersion(out)
    disp
```

This value will be treated as known in our Bayesian model, so that the Bayesian and classical fits are directly comparable in terms of how they treat phi.


### 5.2 Prior Setup Using Prior_Setup()

Next we use `Prior_Setup()` to construct a Zellner‑type g‑prior [@GriffinBrown2010] for the coefficients, aligned with the Gamma log‑link model and weighting structure:

```{r carinsca_prior}
    ps <- Prior_Setup(Cost/Claims ~ Merit + Class,
                      family  = Gamma(link = "log"),
                      weights = Claims)

    mu <- ps$mu
    V  <- ps$Sigma
```

The output from ps (if printed) will include:

- prior means (intercept from a null model, effects centered at zero by default)  
- prior standard deviations (scaled using the pwt argument and the Fisher information)  
- any additional hyperparameters relevant for dispersion modeling

With the default pwt = 0.01, the prior is weakly informative, and the posterior estimates typically remain close to the MLEs.


### 5.3 Fitting the Bayesian Gamma Model

We now call glmb(), specifying the same model and likelihood as in the classical fit, but adding a dNormal prior family that fixes dispersion at the value estimated above:

```{r carinsca_glmb}
    out_glmb <- glmb(Cost/Claims ~ Merit + Class,
                 family  = Gamma(link = "log"),
                 pfamily = dNormal(mu = mu, Sigma = V, dispersion = disp),
                 weights = Claims)
```

The arguments mirror the glm() call, with pfamily providing the prior:

- family = Gamma(link = "log") specifies the likelihood and link  
- pfamily = dNormal(...) specifies a multivariate normal prior for beta with known dispersion  
- weights = Claims carries over the claim counts

### 5.4 Summarizing the Bayesian Model

Bayesian summary:

```{r carinsca_sumamry_glmb}

    summary(out_glmb)
    options(oldopt)
```

The Bayesian summary will report:

- posterior modes and posterior means for each coefficient  
- posterior standard deviations and Monte Carlo errors  
- tail probabilities relative to the prior mean  
- posterior percentiles (credible intervals)  
- effective number of parameters pD  
- expected residual deviance D  
- DIC

Because dispersion is fixed at the same value in both models, differences between out and out_glmb reflect primarily the influence of the prior on \(\beta\), not differences in \(\phi\).


## 6. Interpretation and Model Diagnostics

### 6.1 Coefficients

Under the log link, coefficient interpretations are multiplicative:

- exp(\(\beta_{j}\)) is the multiplicative change in expected cost for a one‑unit change in the covariate (or for a change from baseline to a given factor level).

With a weak prior, posterior means in out_glmb will be close to the MLEs in out, but with slightly more regularization and smaller posterior standard deviations.

### 6.2 Dispersion

By fixing dispersion at disp from the classical model, the Bayesian and classical fits share the same assumed level of variability.  

This allows:

- direct comparison of log‑likelihood based quantities,  
- a clean assessment of how the prior affects coefficient estimates.

If desired, a more advanced model can estimate dispersion by switching to a pfamily such as dNormal_Gamma, but that is covered in more detail in the dispersion‑focused vignettes.


### 6.3 DIC and Model Fit

The Bayesian summary of out_glmb reports:

- the effective number of parameters pD,  
- the expected residual deviance D,  
- the DIC value [@Spiegelhalter2002].

These can be compared to AIC and deviance from the classical model, especially when considering alternative specifications (for example, reduced models, different covariate sets, or alternative link functions).


## 7. Concluding Discussion

Gamma regression provides a flexible framework for modeling positive, skewed response variables [@McCullagh1989; @Gelman2013].  
In this chapter, we have:

- introduced the Gamma GLM and its link functions,  
- shown how to fit a classical Gamma regression using glm(),  
- demonstrated how to estimate dispersion with gamma.dispersion(),  
- used Prior_Setup() to construct a default Zellner‑type prior,  
- fit a Bayesian Gamma regression with glmb() using a dNormal prior with fixed dispersion,  
- and compared classical and Bayesian summaries in a realistic insurance setting.

In later vignettes, we extend these ideas to:

- Gamma models with unknown dispersion and hierarchical priors,  
- Gamma regression combined with other components in larger GLM systems,  
- and specialized accept‑reject schemes for dispersion parameters.
