Generalized linear models relate an outcome mean to predictors while
allowing response distributions such as binomial or Poisson. The
fastglm package performs the fitting calculations in C++ and
provides several ways to solve the weighted regressions used during
fitting. fastglm is compatible with R’s
family objects (see ?family). The package can
be installed via
and loaded via:
The main function, fastglm(), accepts either a model
formula and data frame or a predictor matrix and response vector. A
matrix is convenient when the predictor columns have already been
constructed. The following example illustrates that interface:
data(esoph)
x <- model.matrix(~ agegp + unclass(tobgp) + unclass(alcgp),
data = esoph)
y <- cbind(esoph$ncases, esoph$ncontrols)
gfit1 <- fastglm(x = x, y = y, family = binomial(link = "cloglog"))
summary(gfit1)
#> Generalized linear model
#>
#> Call:
#> fastglm.default(x = x, y = y, family = binomial(link = "cloglog"))
#> Family: binomial (cloglog link)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -4.29199 0.29777 -14.414 <1e-04 ***
#> agegp.L 3.30677 0.63456 5.211 <1e-04 ***
#> agegp.Q -1.35525 0.57311 -2.365 0.018 *
#> agegp.C 0.20296 0.43334 0.468 0.640
#> agegp^4 0.15056 0.29318 0.514 0.608
#> agegp^5 -0.23425 0.17855 -1.312 0.190
#> unclass(tobgp) 0.33276 0.07285 4.568 <1e-04 ***
#> unclass(alcgp) 0.80384 0.07538 10.664 <1e-04 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Dispersion Null deviance Null df Deviance Residual df AIC
#> 1 368 87 96.95 80 228
#> Convergence: converged (6 iterations)A formula can be passed directly to fastglm(). The
fitting routine can also be used inside glm() by supplying
method = fastglm_fit, which is useful when existing
analysis code already uses glm():
gfit2 <- glm(cbind(ncases, ncontrols) ~ agegp + unclass(tobgp) +
unclass(alcgp), data = esoph,
family = binomial(link = "cloglog"),
method = fastglm_fit)
summary(gfit2)
#> Generalized linear model
#>
#> Call:
#> glm(formula = cbind(ncases, ncontrols) ~ agegp + unclass(tobgp) +
#> unclass(alcgp), family = binomial(link = "cloglog"), data = esoph,
#> method = fastglm_fit)
#> Family: binomial (cloglog link)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -4.29199 0.29777 -14.414 <1e-04 ***
#> agegp.L 3.30677 0.63456 5.211 <1e-04 ***
#> agegp.Q -1.35525 0.57311 -2.365 0.018 *
#> agegp.C 0.20296 0.43334 0.468 0.640
#> agegp^4 0.15056 0.29318 0.514 0.608
#> agegp^5 -0.23425 0.17855 -1.312 0.190
#> unclass(tobgp) 0.33276 0.07285 4.568 <1e-04 ***
#> unclass(alcgp) 0.80384 0.07538 10.664 <1e-04 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Dispersion Null deviance Null df Deviance Residual df AIC
#> 1 368 87 96.95 80 228
#> Convergence: converged (6 iterations)The coefficient summary uses the covariance calculated during fitting, evaluated at the reported coefficient estimates. It therefore obtains standard errors without fitting the model again.
A GLM is usually fitted through a sequence of weighted regression updates, called iteratively reweighted least squares (IRLS). An update can increase the discrepancy between the fitted and observed outcomes or produce means outside the distribution’s allowed range. The fastglm safeguards shorten such steps until an acceptable update is found. This step-halving procedure can improve convergence, although difficult models can still fail.
The following example deliberately fits an inappropriate link, the function relating the response mean to the linear predictor. It illustrates numerical difficulties that can arise when proposed means leave the allowed range. Different starting values and step safeguards can lead to different convergence behavior. Numerical convergence alone does not make the chosen response distribution or link scientifically appropriate.
set.seed(1)
x <- matrix(rnorm(10000 * 100), ncol = 100)
y <- exp(0.25 * x[,1] - 0.25 * x[,3] + 0.5 * x[,4] - 0.5 * x[,5] + rnorm(10000)) + 0.1
# fastglm
system.time(gfit1 <- glm(y ~ x, family = Gamma(link = "sqrt"),
method = fastglm_fit))
#> user system elapsed
#> 0.299 0.007 0.307
# glm
system.time(gfit2 <- glm(y ~ x, family = Gamma(link = "sqrt")))
#> user system elapsed
#> 1.566 0.035 1.603
# glm2
system.time(gfit3 <- glm(y ~ x, family = Gamma(link = "sqrt"),
method = glm2::glm.fit2))
#> user system elapsed
#> 1.125 0.029 1.155
## Note that fastglm() returns estimates with the
## largest likelihood
logLik(gfit1)
#> 'log Lik.' -16046.66 (df=102)
logLik(gfit2)
#> 'log Lik.' -16704.05 (df=102)
logLik(gfit3)
#> 'log Lik.' -16046.66 (df=102)
coef(gfit1)[1:5]
#> (Intercept) x1 x2 x3 x4
#> 1.426850e+00 1.242357e-01 -7.475901e-05 -1.254684e-01 2.361046e-01
coef(gfit2)[1:5]
#> (Intercept) x1 x2 x3 x4
#> 1.431168e+00 1.251936e-01 -6.896739e-05 -1.281857e-01 2.366473e-01
coef(gfit3)[1:5]
#> (Intercept) x1 x2 x3 x4
#> 1.426864e+00 1.242616e-01 -9.860241e-05 -1.254873e-01 2.361301e-01
## check convergence of fastglm
gfit1$converged
#> [1] TRUE
## number of IRLS iterations
gfit1$iter
#> [1] 16
## now check convergence for glm()
gfit2$converged
#> [1] FALSE
gfit2$iter
#> [1] 25
## check convergence for glm2()
gfit3$converged
#> [1] TRUE
gfit3$iter
#> [1] 19
## increasing number of IRLS iterations for glm() does not help that much
system.time(gfit2 <- glm(y ~ x, family = Gamma(link = "sqrt"),
control = list(maxit = 100)))
#> user system elapsed
#> 5.923 0.131 6.063
gfit2$converged
#> [1] FALSE
gfit2$iter
#> [1] 100
logLik(gfit1)
#> 'log Lik.' -16046.66 (df=102)
logLik(gfit2)
#> 'log Lik.' -16068.27 (df=102)