Enriching family objects: exponential family of distirbutions

Introduction

The enrichwith R package provides the enrich method to enrich list-like R objects with new, relevant components. The resulting objects preserve their class, so all methods associated with them still apply.

This vignette is a demo of the available enrichment options for family objects, focusing on objects that correspond to members of the exponential family of distributions.

Exponential family and family objects

family objects specify characteristics of the models used by functions such as glm. The families implemented in the stats package include binomial, gaussian, Gamma, inverse.gaussian, and poisson, which obvious corresponding distributions. These distributions are all special cases of the exponential family of distributions with probability mass or density function of the form $$ f_{Y}(y) = \exp\left\{\frac{y \theta - b(\theta) - c_1(y)}{\phi/m} - \frac{1}{2}a\left(-\frac{m}{\phi}\right) + c_2(y) \right\} $$ for some sufficiently smooth functions b(.), c1(.), a(.) and c2(.), and a fixed weight m. The expected value and the variance of Y are then where V(μ) and ϕ are the variance function and the dispersion parameter, respectively. Below we list the characteristics of the distributions supported by family objects.

Normal with mean μ and variance ϕ/m

θ = μ, $\displaystyle b(\theta) = \frac{\theta^2}{2}$, $\displaystyle c_1(y) = \frac{y^2}{2}$, a(ζ) = −log (−ζ), $\displaystyle c_2(y) = -\frac{1}{2}\log(2\pi)$

Binomial with index m and probability μ

$\displaystyle \theta = \log\frac{\mu}{1- \mu}$, b(θ) = log (1 + eθ), ϕ = 1, c1(y) = 0, a(ζ) = 0, $\displaystyle c_2(y) = \log{m\choose{my}}$

Poisson with mean μ

θ = log μ, b(θ) = eθ, ϕ = 1, c1(y) = 0, a(ζ) = 0, c2(y) = −log Γ(y + 1)

Gamma with mean μ and shape 1/ϕ

$\displaystyle \theta = -\frac{1}{\mu}$, b(θ) = −log (−θ), c1(y) = −log y, a(ζ) = 2log Γ(−ζ) + 2ζlog (−ζ), c2(y) = −log y

Inverse Gaussian with mean μ and variance ϕμ3

$\displaystyle \theta = -\frac{1}{2\mu^2}$, $\displaystyle b(\theta) = -\sqrt{-2\theta}$, $\displaystyle c_1(y) = \frac{1}{2y}$, a(ζ) = −log (−ζ), $\displaystyle c_2(y) = -\frac{1}{2}\log\left(\pi y^3\right)$

Components in family objects

family objects provide functions for the variance function (variance), a speficiation of deviance residuals (dev.resids) and the Akaike information criterion (aic). For example

inverse.gaussian()$dev.resids
## function (y, mu, wt) 
## wt * ((y - mu)^2)/(y * mu^2)
## <bytecode: 0x5586a2ea8750>
## <environment: 0x5586a2eac280>
inverse.gaussian()$variance
## function (mu) 
## mu^3
## <bytecode: 0x5586a2ea8910>
## <environment: 0x5586a2f0f180>
inverse.gaussian()$aic
## function (y, n, mu, wt, dev) 
## sum(wt) * (1 + log(dev/sum(wt) * 2 * pi)) + 3 * sum(log(y) * 
##     wt) + 2
## <bytecode: 0x5586a2ea8088>
## <environment: 0x5586a2f6b2b8>

Enrichment options for family objects

The enrichwith R package provides methods for the enrichment of family objects with a function that links the natural parameter θ with μ, the function b(θ), the first two derivatives of V(μ), a(ζ) and its first four derivatives, and c1(y) and c2(y). To illustrate, let’s write a function that reconstructs the densitities and probability mass functions from the components that result from enrichment

library("enrichwith")
dens <- function(y, m = 1, mu, phi, family) {
    object <- enrich(family)
    with(object, {
        c2 <- if (family == "binomial") c2fun(y, m) else c2fun(y)
        exp(m * (y * theta(mu) - bfun(theta(mu)) - c1fun(y))/phi -
            0.5 * afun(-m/phi) + c2)
    })
}

The following chunks test dens for a few distributions against the standard d* functions

## Normal
all.equal(dens(y = 0.2, m = 3, mu = 1, phi = 3.22, gaussian()),
          dnorm(x = 0.2, mean = 1, sd = sqrt(3.22/3)))
## [1] TRUE
## Gamma
all.equal(dens(y = 3, m = 1.44, mu = 2.3, phi = 1.3, Gamma()),
          dgamma(x = 3, shape = 1.44/1.3, 1.44/(1.3 * 2.3)))
## [1] TRUE
## Inverse gaussian
all.equal(dens(y = 0.2, m = 7.23, mu = 1, phi = 3.22, inverse.gaussian()),
          SuppDists::dinvGauss(0.2, nu = 1, lambda = 7.23/3.22))
## [1] TRUE
## Binomial
all.equal(dens(y = 0.34, m = 100, mu = 0.32, phi = 1, binomial()),
          dbinom(x = 34, size = 100, prob = 0.32))
## [1] TRUE