Choosing the Right Link Function
Generalised Linear Models (GLM) offer a flexible framework for regression analysis. When building a GLM for a practical problem, the first and most important decision is selecting an appropriate link function — and this choice is driven by the nature of the response variable. A binary response calls for the logit link, count data calls for the log link, a proportion response calls for the probit or complementary log-log link, and so on. Once the link function is fixed, the model is ready to be fitted with data.
Fitting the Model: Iterative Optimisation Algorithms
GLMs are not fitted in one step. Since Maximum Likelihood Estimation for GLMs generally has no closed-form solution, parameter estimates are arrived at through numerical iterative optimisation algorithms. Starting from an initial guess, the algorithm updates the parameter estimates step by step, each iteration bringing the estimates closer to the values that maximise the likelihood.
Commonly used iterative algorithms in GLM fitting include:
- Newton-Raphson — uses both the gradient and the Hessian (matrix of second derivatives of the log-likelihood). Converges fast and accurately. Naturally handles predictors on different scales due to the curvature information in the Hessian.
- Fisher Scoring — a variant of Newton-Raphson that replaces the observed Hessian with the expected (Fisher) information matrix. Equivalent to Newton-Raphson for canonical link functions such as logit and log.
- Iteratively Reweighted Least Squares (IRLS) — the classical GLM fitting algorithm. At each step it solves a weighted least squares problem. Mathematically equivalent to Fisher Scoring for GLMs.
- Gradient Descent — uses only the first derivative (gradient). Computationally cheap per iteration but slow to converge and sensitive to feature scale.
- Stochastic Gradient Descent (SGD) — a variant that updates parameters using one or a small batch of observations at a time. Useful for very large datasets but introduces noise in convergence.
- L-BFGS (Limited-memory Broyden–Fletcher–Goldfarb–Shanno) — a quasi-Newton method that approximates the Hessian using gradient history from recent iterations. Memory-efficient but sensitive to differences in predictor scale.
- Conjugate Gradient (CG) — an iterative method that uses successive conjugate directions to minimise the objective. More efficient than plain gradient descent but still first-order in nature.
- Coordinate Descent — optimises one parameter at a time while keeping the rest fixed, cycling through all parameters. Very efficient for penalised regression problems such as Lasso and Elastic Net.
The essential aspect across all these algorithms is the stopping rule — the criterion that decides when the iteration ends. Two common approaches are:
- A predetermined number of steps — the algorithm runs for a fixed number of iterations and stops
- A tolerance limit — the algorithm compares parameter estimates from two consecutive steps and stops when the difference falls below a pre-set threshold, indicating convergence
R Libraries for GLM Fitting
R offers a rich set of packages for fitting GLMs. Each package is designed with a specific use case in mind — standard inference, penalised regression, bias correction, or handling large data. The following are the commonly used ones.
1. stats::glm() — Base R
The built-in glm() function in base R is the starting point for most GLM analysis. It uses IRLS (Fisher Scoring) as its fitting algorithm and covers the full range of exponential family distributions through the family argument.
model <- glm(y_response ~ age + balance + day + duration + campaign + pdays + previous,
data = mydata,
family = binomial(link = "logit"))
summary(model)
Output includes coefficients, standard errors, z-statistics, p-values, deviance, and AIC — no additional setup required. Convergence control is available through glm.control() where tolerance and maximum iterations can be set.
- CRAN / Documentation: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/glm.html
2. glmnet — Penalised GLM (Lasso, Ridge, Elastic Net)
glmnet fits GLMs with regularisation — L1 (Lasso), L2 (Ridge), or a combination (Elastic Net). It is particularly useful when the number of predictors is large or when multicollinearity is present. The fitting algorithm is Cyclical Coordinate Descent, which is highly efficient for penalised problems.
library(glmnet)
x <- model.matrix(y_response ~ age + balance + day + duration + campaign + pdays + previous,
data = mydata)[, -1]
y <- mydata$y_response
model <- glmnet(x, y, family = "binomial", alpha = 1) # alpha=1 for Lasso
plot(model)
glmnet fits the entire regularisation path across a sequence of lambda values in one call. Cross-validation for selecting the optimal lambda is available through cv.glmnet().
- CRAN: https://cran.r-project.org/package=glmnet
3. MASS::glm.nb() — Negative Binomial GLM
The MASS package extends base R GLM to handle Negative Binomial regression, which is appropriate for overdispersed count data where the Poisson assumption of equal mean and variance does not hold. It uses an alternating iteration between estimating the coefficients (via IRLS) and estimating the dispersion parameter.
library(MASS)
model <- glm.nb(count_response ~ x1 + x2 + x3, data = mydata)
summary(model)
- CRAN: https://cran.r-project.org/package=MASS
4. brglm — Bias-Reduced GLM for Binomial Response
Standard MLE in logistic regression can produce infinite estimates when complete or quasi-complete separation occurs in the data — a situation that is not uncommon in small or sparse datasets. brglm addresses this using Firth’s bias reduction method, which modifies the score equations to produce finite, second-order unbiased estimates.
library(brglm)
model <- brglm(y_response ~ age + balance + day + duration,
family = binomial(link = "logit"),
data = mydata)
summary(model)
The interface is nearly identical to glm(), making it an easy drop-in replacement when separation is suspected.
- CRAN: https://cran.r-project.org/package=brglm
5. brglm2 — Extended Bias Reduction for GLMs
brglm2 is the more recent and general successor to brglm. It supports bias reduction across a wider range of GLM families beyond binomial — including Poisson, Gamma, and others. It fits via a quasi Fisher Scoring algorithm and supports both mean bias reduction and median bias reduction.
library(brglm2)
model <- glm(y_response ~ age + balance + day + duration,
family = binomial(link = "logit"),
data = mydata,
method = "brglmFit")
summary(model)
- CRAN: https://cran.r-project.org/package=brglm2
6. speedglm — Fast GLM for Large Datasets
speedglm is designed for fitting GLMs efficiently to large datasets. It uses IRLS but avoids the overhead of base R’s glm() by working directly with matrix operations and optionally leveraging optimised BLAS routines. For datasets that exceed memory, it supports an iterative updating procedure through shglm().
library(speedglm)
model <- speedglm(y_response ~ age + balance + day + duration + campaign + pdays + previous,
data = mydata,
family = binomial(link = "logit"))
summary(model)
- CRAN: https://cran.r-project.org/package=speedglm
7. biglm — GLM for Data Too Large to Fit in Memory
biglm provides bigglm(), which fits GLMs on datasets that are too large to load into memory at once. It uses a chunking algorithm — reading and processing data in pieces, updating the sufficient statistics incrementally until convergence. The stopping rule checks the change in coefficients relative to their standard errors against a tolerance threshold.
library(biglm)
model <- bigglm(y_response ~ age + balance + day + duration,
data = mydata,
family = binomial(),
chunksize = 500)
summary(model)
- CRAN: https://cran.r-project.org/package=biglm
Summary
| Package | Algorithm | Best For |
|---|---|---|
stats::glm() | IRLS / Fisher Scoring | Standard GLM inference |
glmnet | Coordinate Descent | Penalised regression, high-dimensional data |
MASS::glm.nb() | IRLS + dispersion iteration | Overdispersed count data |
brglm | Modified Score (Firth) | Separation in binomial GLMs |
brglm2 | Quasi Fisher Scoring | Bias reduction across all GLM families |
speedglm | Optimised IRLS | Large datasets, fast fitting |
biglm | Chunked IRLS | Data too large to fit in memory |
R provides a package for virtually every practical GLM scenario — from standard inference to separation-prone sparse data to out-of-memory datasets. The base glm() function handles the majority of everyday situations cleanly, and the ecosystem of packages extends naturally from there without imposing any unnatural restrictions on the analyst.
References
- Nelder, J.A. and Wedderburn, R.W.M. (1972). Generalized Linear Models. Journal of the Royal Statistical Society, Series A, 135(3), 370–384.
- McCullagh, P. and Nelder, J.A. (1989). Generalized Linear Models (2nd ed.). Chapman and Hall.
- Friedman, J., Hastie, T. and Tibshirani, R. (2010). Regularization Paths for Generalized Linear Models via Coordinate Descent. Journal of Statistical Software, 33(1). https://cran.r-project.org/package=glmnet
- Firth, D. (1993). Bias reduction of maximum likelihood estimates. Biometrika, 80(1), 27–38.
- Kosmidis, I. et al. (2020). Mean and median bias reduction in generalized linear models. Statistics and Computing, 30, 43–59. https://cran.r-project.org/package=brglm2
- Enea, M. (2009). Fitting Linear Models and Generalized Linear Models with large data sets in R. https://cran.r-project.org/package=speedglm
- Lumley, T. (2006). biglm: Bounded Memory Linear and Generalized Linear Models. https://cran.r-project.org/package=biglm
All web links cited in this write-up were accessed and the content verified during the preparation of this document in May 2026.
Also refer: Transformation in GLM