10  Other GLMs

In the previous chapters, we covered important extensions of the ordinary least squares (OLS) model: logistic regression for binary responses and Poisson regression for count responses. It turns out that these two models are special cases of an important and powerful family of models called generalized linear models (GLMs) (Nelder and Wedderburn 1972). In this chapter, we describe GLMs and explain how to fit regularized versions of GLMs using glmnet.

10.1 What are GLMs?

A generalized linear model (GLM) consists of 3 parts:

• A linear predictor: \(\eta_i = \mathbf{x}_i^\top \mathbf{\beta}\),

• A link function: \(\eta_i = g(\mu_i)\), and

• A variance function as a function of the mean \(V = V(\mu_i)\).

The response \(y_i\) is assumed to come from a distribution¹ with mean \(\mu_i\) and variance \(\mbox{Var}(y_i)=V(\mu_i)\). The user gets to specify the link function \(g\) and the variance function \(V\). GLMs cover a wide range of models, including the ones that we’ve seen in previous chapters:

GLM family / Regression type	Link function	Variance function	Response distribution
Gaussian regression (Section 1.3)	Identity function \(g(\mu) = \mu\)	Constant function \(V(\mu) = \sigma^2\)	\(y \sim \mathcal{N}(\mu, \sigma^2)\)
Logistic regression (Section 7.1)	\(g(\mu) = \log \left( \dfrac{\mu}{1-\mu} \right)\)	\(V(\mu) = \mu (1 - \mu)\)	\(y \sim \text{Bern}(\mu)\)
Poisson regression (Section 9.1)	\(g(\mu) = \log \mu\)	\(V(\mu) = \mu\)	\(y \sim \text{Poi}(\mu)\)

Here is a longer list of GLM examples with their representations in R:

Examples of generalized linear models (GLMs) and their representations in R.

GLM family / Regression type	Response type	Representation in R
Gaussian	\(\mathbb{R}\)	gaussian()
Logistic	\(\{ 0, 1 \}\)	binomial()
Probit	\(\{ 0, 1 \}\)	binomial(link = "probit")
Quasi-Binomial	\(\{ 0, 1 \}\)	quasibinomial()
Poisson	\(\mathbb{N}_0 = \{ 0, 1, \dots \}\)	poisson()
Quasi-Poisson	\(\mathbb{N}_0\)	quasipoisson()
Negative binomial	\(\mathbb{N}_0\)	MASS::negative.binomial()
Gamma	\(\mathbb{R}_+ = [0, \infty)\)	Gamma()
Inverse Gaussian	\(\mathbb{R}_+\)	inverse.gaussian()
Tweedie	Depends on variance power parameter	statmod::tweedie()

10.2 Fitting a GLM

In OLS, we fit the model by minimizing the residual sum of squares. In GLMs, we fit the model by maximum likelihood estimation. Let \(\ell(y_i, \beta_0 + \mathbf{x}_i^\top \beta)\) denote the log-likelihood term associated with observation \(i\). In maximum likelihood estimation, we want to find \(\hat\beta_0\) and \(\hat\beta\) such that

\[(\hat\beta_0, \hat\beta) = \underset{(\beta_0, \beta) \in \mathbb{R}^{p+1}}{\text{argmin}} \; -\sum_{i=1}^n \ell (y_i, \beta_0 + \mathbf{x}_i^\top \beta). \tag{10.1}\]

This turns out to coincide with minimizing the residual sum of squares (OLS) when fitting Gaussian models. For the Gaussian family, \(\ell(y_i, \beta_0 + \mathbf{x}_i^\top \beta) = \frac{1}{2}\left[ y_i- (\beta_0 + \mathbf{x}_i^\top \beta) \right]^2\). Plugging this into Equation 10.1, we obtain Equation 1.1.

It is easy to fit GLMs in R using the glm function from the stats package; the user can specify the GLM to be fit using class family objects. These objects capture details of the GLM, such as the link function and the variance function.

As an example, probit regression is a model for binary responses that assumes that

\[y \sim \text{Bern}(\mu), \quad \mu = \Phi (\mathbf{x}^\top \beta),\]

where \(\Phi(\cdot)\) is the cumulative distribution function of the standard Gaussian random variable. The implied link function here is \(\eta = \Phi^{-1}(\mu)\). The code below shows the family object associated with the probit regression model:

class(binomial(link = "probit"))
#> [1] "family"
str(binomial(link = "probit"))
#> List of 13
#>  $ family    : chr "binomial"
#>  $ link      : chr "probit"
#>  $ linkfun   :function (mu)  
#>  $ linkinv   :function (eta)  
#>  $ variance  :function (mu)  
#>  $ dev.resids:function (y, mu, wt)  
#>  $ aic       :function (y, n, mu, wt, dev)  
#>  $ mu.eta    :function (eta)  
#>  $ initialize: language {     if (NCOL(y) == 1) { ...
#>  $ validmu   :function (mu)  
#>  $ valideta  :function (eta)  
#>  $ simulate  :function (object, nsim)  
#>  $ dispersion: num 1
#>  - attr(*, "class")= chr "family"

# link function for probit regression
binomial(link = "probit")$linkfun
#> function (mu) 
#> qnorm(mu)
#> <environment: namespace:stats>

# variance function for probit regression
binomial(link = "probit")$variance
#> function (mu) 
#> mu * (1 - mu)
#> <bytecode: 0x55bbe17285c0>
#> <environment: 0x55bbe25740f0>

By passing a class "family" object to the family argument of a glm call, glm has all the information it needs to fit the model. Here is an example of how to fit a probit regression model in R:

data("BinomialExample")
str(BinomialExample)
#> List of 2
#>  $ x: num [1:100, 1:30] -0.619 1.094 -0.357 -2.469 0.567 ...
#>  $ y: int [1:100] 0 1 1 0 1 0 0 0 1 1 ...
fit <- glm(y ~ x, data = BinomialExample, family = binomial(link = "probit"))

10.3 Regularized GLMs

To fit an elastic net model for a GLM, we simply add the elastic net penalty to the objective function to be minimized in Equation 10.1:

\[(\hat\beta_0, \hat\beta) = \underset{(\beta_0, \beta) \in \mathbb{R}^{p+1}}{\text{argmin}} \; \left[-\frac{1}{n} \sum_{i=1}^n \ell (y_i, \beta_0 + \mathbf{x}_i^\top \beta) + \lambda \left(\frac{1-\alpha}{2}\|\beta \|_2^2 + \alpha \|\beta \|_1 \right) \right].\]

Fitting such a regularized GLM in R is easy with glmnet: all one has to do is specify the family argument in the function call.

10.3.1 Specialized families

There are a few commonly used GLMs for which we have specialized, hyper-efficient routines for fitting. For these GLMs, we specify the family argument as a character string. The options here are:

• family="gaussian"

• family="binomial" (with default logit link)

• family="poisson"

If you are fitting one of these common GLM families, you should use the string family argument.

10.3.2 All other families

To fit any other regularized GLM, we can pass an object of class "family" to the family argument. For example, the code below fits a regularized quasi-Poisson model that allows for overdispersion:

glmnet(x, y, family = quasipoisson())

10.4 Additional notes

We discussed the "multinomial" family in the Multiclass Section. Although this family is closely aligned with GLMs, the response is not a scalar. Multinomial models are more similar in spirit to the models for multiple responses discussed in the following Multitask section.

For algorithmic details on how to fit unregularized GLMs as well as how glmnet fits regularized GLMs, see Tay et al. (2023).

Nelder, J. A., and R. W. M. Wedderburn. 1972. “Generalized Linear Models.” Journal of the Royal Statistical Society. Series A (General) 135 (3): 370. https://doi.org/10.2307/2344614.

Tay, J. Kenneth, Balasubramanian Narasimhan, and Trevor Hastie. 2023. “Elastic Net Regularization Paths for All Generalized Linear Models.” Journal of Statistical Software 106 (1). https://doi.org/10.18637/jss.v106.i01.

---

1. The distribution of the response \(y\) must come from a particular distribution in an exponential family. Exponential families include many of probability distributions we are familiar with, e.g. the normal, binomial, Poisson and gamma distributions.