The fastglm package provides regression models for continuous, binary, count, and ordered outcomes, with an emphasis on efficient computation and numerical stability. Its scope includes model-specific mean bias reduction, coefficient and model-assumption tests, probability and mean-response summaries, and methods for repeated fitting or data too large to retain as a dense design. This overview identifies the available models and methods before developing a common workflow for specification, diagnostics, uncertainty, and repeated outcomes.
| Outcome or modeling need | Main interface | Estimation and worked guide |
|---|---|---|
| Generalized linear models (GLMs): continuous responses, binary outcomes, or counts | fastglm(); also
glm(method=fastglm_fit) |
Likelihood or quasi-likelihood fitting; mean bias reduction for supported likelihood family/link combinations. GLM specification and bias reduction |
| Counts with variation beyond the Poisson model | fastglm_nb() |
Maximum likelihood (ML) estimation of regression coefficients and negative-binomial size. Count models |
| Counts sampled conditional on being positive | fastglm_truncated() |
Zero-truncated Poisson or negative-binomial likelihood. Truncation and hurdles |
| Separate models for occurrence and positive counts | fastglm_hurdle() |
Binary occurrence and zero-truncated Poisson or negative-binomial components, fitted by ML. Hurdle models |
| Zeros arising from an additional component and a count distribution | fastglm_zi() |
Zero-inflated Poisson or negative-binomial likelihood. Zero-inflated models |
| A point mass at zero and a positive continuous amount | fastglm_twopart() |
Binary and Gamma components, with ML or mean bias reduction. Two-part models |
| Ordered outcome categories | fastglm_ordinal() |
Common or selected threshold-specific slopes; logit, probit, and complementary log-log links; ML or mean bias reduction. Ordinal models, inference, and computation |
Mean bias reduction is requested with firth=TRUE in the
interfaces that support it. For GLMs, support depends on the family and
link; a fixed-size negative-binomial family supplied through
negbin(theta) differs from jointly estimating size with
fastglm_nb(). The truncated, hurdle, and zero-inflated
count interfaces above currently use ML. The bias-reduction guide lists the supported
GLM combinations, while the ordinal and two-part guides describe their
respective adjustments.
| Analysis task | Interface | Scope |
|---|---|---|
| Coefficient intervals and joint tests | confint() and
fastglm_wald_test() |
Model-based Wald inference for supported fitted classes; formula terms and named contrasts can combine several coefficients |
| Profile intervals and likelihood-ratio tests | confint(method="profile"),
fastglm_profile_test(), and
fastglm_confint_theta() |
Dense binomial-logit Firth fits, dense ordinal ML fits with common or varying slopes, and negative-binomial size profiles, respectively |
| Score tests of coefficients or the common-slope restriction | fastglm_score_test() and
fastglm_po_test() |
Retained dense ordinal models for coefficient
constraints; fastglm_po_test() tests complete terms in
dense or streamed common-slope and partial models |
| Robust or clustered standard errors | fastglm_vcov() and supported
vcov(type=...) methods |
GLMs, count and two-part models, and dense or streamed ordinal fits, subject to each class’s requirements |
| Response or probability summaries | predict(se.fit=TRUE),
fastglm_average(), and
fastglm_marginal_effects() |
Predictions and averages for supported classes; automatic marginal effects target common-slope or partial ordinal probabilities |
| Refitting uncertainty and model assessment | fastglm_bootstrap(),
fastglm_diagnostics(), residuals, and
fastglm_ordinal_assess() |
Bootstrap support depends on retained data and model class; ordinal prediction scores, calibration, and optional boundary-row diagnostics support dense or streamed fits |
These methods use different assumptions and are not interchangeable across model classes. Dense partial ordinal fits support robust covariance, ML profiles, and constrained score inference. Streamed ordinal fits provide model-based and robust covariance, common-slope tests, prediction scores, and calibration on the original source. General constrained coefficient inference requires retained observations. The ordinal support guide states these distinctions explicitly. The sections below explain coefficient uncertainty and population averages, joint tests, and readable result displays.
| Computational setting | Available approach | Scope and guide |
|---|---|---|
| Ordinary in-memory analysis | Formula and dense-design interfaces | Model-specific fitting with coefficient, prediction, and diagnostic methods |
| Predictors with many zeros or stored outside ordinary R matrices | Sparse matrices and bigmemory::big.matrix
inputs |
GLMs; selected count interfaces also accept sparse designs. Large-data guide |
| Data read repeatedly from files or databases | fastglm_source(), Arrow and DBI adapters,
and streaming fitters |
GLMs, negative-binomial, hurdle, two-part, and common-slope or partial ordinal models. Streaming computation |
| Several responses with the same predictors | fastglm_prepare() /
fastglm_batch() and their ordinal counterparts |
GLMs and common-slope or partial ordinal models; supported calculations are reused across responses. Repeated outcomes |
Data representation changes the amount of retained information and the cost of rereading observations. It can therefore affect both fitting time and which later procedures are available. The large-data vignette explains these requirements. Benchmark examples examine speed together with numerical agreement.
The remaining sections demonstrate the shared workflow with a Poisson model: specify the mean, choose a numerical decomposition and assess convergence, quantify uncertainty, and reuse a design. Generalized linear model fitting uses iteratively reweighted least squares, which updates coefficients through a sequence of weighted regressions.
A formula specifies the response and predictors using the same
notation as glm(). The fitted object records factor levels,
the coding of categories into predictor columns, and offsets so that
subsequent predictions use the same model specification. An offset
represents a known contribution to the linear predictor. In the
following example, exposure multiplies the expected count, so its
logarithm enters the linear predictor with coefficient one.
library(fastglm)
set.seed(41)
d <- data.frame(x=rnorm(500), group=factor(rep(letters[1:4],125)),
exposure=runif(500,.5,2))
d$y <- rpois(500, d$exposure*exp(.2+.4*d$x))
fit <- fastglm(y~x+group+offset(log(exposure)), data=d, family=poisson(),
method="auto")
summary(fit)
#> Generalized linear model
#>
#> Call:
#> fastglm::fastglm(data = d, family = poisson(), method = "auto",
#> formula = y ~ x + group + offset(log(exposure)))
#> Family: poisson (log link)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 0.04085 0.07493 0.545 0.5856
#> x 0.41174 0.03465 11.884 <1e-04 ***
#> groupb 0.16419 0.09947 1.651 0.0988 .
#> groupc 0.15478 0.10065 1.538 0.1241
#> groupd 0.14228 0.10240 1.389 0.1647
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Dispersion Null deviance Null df Deviance Residual df AIC
#> 1 688.8 499 548.6 495 1504
#> Convergence: converged (5 iterations)
predict(fit, newdata=d[1:4,], type="response", se.fit=TRUE)
#> Predictions with standard errors
#> Scale: response mean
#> Residual scale: 1
#>
#> Row Estimate Std. Error
#> 1 0.8993 0.07737
#> 2 2.2479 0.15343
#> 3 2.1857 0.15936
#> 4 2.6889 0.21223For computations that already provide a design matrix,
fastglm(X, y, family) avoids constructing a model frame.
Prediction from a matrix fit requires the same column order as the
fitted design. New offsets can be supplied through the
offset argument to predict(). The fitting
method can also be used through
glm(..., method=fastglm_fit).
The optional broom::tidy() and
broom::glance() methods provide coefficient and model
summaries for matrix and formula fits. Standard update()
calls can revise a formula specification.
Fitted models and their summaries use a common display: a model
heading, labeled coefficient tables, and fitting status. The same
conventions apply to count, ordinal, streaming, and bootstrap results.
print(summary(fit),digits=3) controls displayed precision
while retaining the original numerical values. Predictions requested
with se.fit=TRUE place estimates beside their standard
errors. Their as.data.frame() method provides a numeric
table for reporting; the original vectors or matrices remain accessible
as $fit and $se.fit. For matrix predictions,
the table identifies both the observation row and the component, such as
an outcome category. The display shows at most 20 rows by default, with
max.rows=Inf displaying the full table.
Hypothesis-test results from fastglm_wald_test(),
fastglm_profile_test(), fastglm_score_test(),
and fastglm_po_test() use a compact statistic,
degrees-of-freedom, and p-value table. Their
as.data.frame() and optional broom::tidy()
methods provide full-precision numeric results for reports, with
details=TRUE adding method descriptions and explanations
for unavailable tests. Profile confidence intervals also print as a
compact table while retaining endpoint and fit diagnostics as
attributes.
For grouped binomial observations, the response can be a two-column matrix containing successes and failures. Trial totals enter the initialized prior weights and are retained for likelihood calculations. Observations with zero prior weight do not contribute to the residual degrees of freedom.
Each coefficient update requires solving a weighted regression. A matrix decomposition rewrites the design or its crossproduct in a form that makes this calculation possible. The available decompositions differ in speed and in their ability to cope with predictors that are redundant or nearly redundant.
| Name | Integer | Decomposition |
|---|---|---|
qr |
0 | Column-pivoted QR |
qr_unpivoted |
1 | Unpivoted QR, requiring full rank |
llt |
2 | Cholesky |
ldlt |
3 | LDLT |
qr_full |
4 | Full-pivoted QR |
svd |
5 | Singular value decomposition |
For ordinary matrices, method="auto" first rescales
predictor columns to reduce differences in their numerical magnitude,
then attempts Cholesky. If the calculation detects redundant or nearly
redundant predictors, it uses column-pivoted QR, which reorders columns
to identify a set that can be estimated separately. The returned
solver and solver.fallback entries record the
choice.
A coefficient is aliased when its effect cannot be distinguished from
a combination of other predictor effects. Pivoted QR reports such
coefficients as NA; fitted means use the remaining
coefficients. Predictions for new rows require the same dependence
between predictors as in the fitted design. Sparse matrices and matrices
accessed through bigmemory support the Cholesky methods and
require that all coefficient effects can be distinguished.
fastglm_diagnostics(fit)
#> Numerical diagnostics
#> Convergence: converged (5 iterations)
#> Solver: Cholesky (LLT)
#> Estimable coefficients (rank): 5
#> Boundary flag: no
#>
#> Check Value
#> Standardized ordinary score 2.131e-11
#> Covariance condition number 33.62A convergence flag means that the algorithm met its stopping rule. It
does not establish that every coefficient is well estimated.
fastglm_diagnostics() reports the score, which measures how
much the fitting objective would change with a small coefficient change,
and a standardized version that accounts for predictor scale. A large
standardized score can indicate incomplete numerical convergence. Fitted
probabilities near zero or one and nearly redundant predictors also
warrant examination. In logistic regression, predictors that perfectly
separate successes from failures can prevent finite maximum likelihood
estimates; the bias-reduction vignette describes an alternative for that
setting.
The rank is the number of coefficient effects that can be estimated
separately. A large covariance condition number indicates that some
combinations of coefficient estimates are much less precisely determined
than others; it depends on the predictor units as well as their
dependence. The boundary flag records whether the fit encountered a
parameter boundary. An unavailable check is labeled explicitly.
as.data.frame(fastglm_diagnostics(fit)) retains numeric and
logical fields for further analysis, with vector-valued details, such as
the individual coefficient scores, stored in list columns.
For an ordinal model,
fastglm_diagnostics(fit,boundary=TRUE) also identifies the
rows and outcome categories with the smallest fitted probabilities and
adjacent cumulative-predictor gaps. Their observed frequencies help
distinguish a rare fitted outcome from an unobserved category at a
particular predictor pattern. This optional calculation rereads a
streamed source once and retains only a few extreme rows. Small
probabilities can help explain numerical difficulty, but do not by
themselves establish that a solution is absent.
The coefficient covariance describes uncertainty in each estimate and how estimation errors in different coefficients vary together. Its diagonal entries are squared standard errors. The model-based calculation uses the fitted response distribution and is evaluated at the reported coefficient estimates.
confint() constructs Wald intervals by adding and
subtracting a multiple of each standard error. The multiplier comes from
a normal reference distribution when dispersion is fixed and a Student t
distribution when dispersion is estimated. Dispersion describes the
response variation beyond the family’s mean-variance relationship. These
intervals use a local approximation around the estimate;
profile-likelihood intervals instead refit the other parameters while
varying the parameter of interest.
vcov(fit)
#> (Intercept) x groupb groupc
#> (Intercept) 0.0056142098 -0.0006491624 -0.0053535815 -0.0053766734
#> x -0.0006491624 0.0012004262 0.0001672107 0.0002099121
#> groupb -0.0053535815 0.0001672107 0.0098947441 0.0052923971
#> groupc -0.0053766734 0.0002099121 0.0052923971 0.0101307820
#> groupd -0.0053841673 0.0002237697 0.0052943274 0.0053022873
#> groupd
#> (Intercept) -0.0053841673
#> x 0.0002237697
#> groupb 0.0052943274
#> groupc 0.0053022873
#> groupd 0.0104862176
confint(fit, parm="x")
#> 2.5 % 97.5 %
#> x 0.3438348 0.4796492
vcov(fit, type="HC1")
#> (Intercept) x groupb groupc
#> (Intercept) 0.0051767554 -5.571796e-04 -0.0049406554 -0.0049673130
#> x -0.0005571796 1.005518e-03 0.0001305478 0.0001789602
#> groupb -0.0049406554 1.305478e-04 0.0092242203 0.0048914533
#> groupc -0.0049673130 1.789602e-04 0.0048914533 0.0095847217
#> groupd -0.0049029537 5.966093e-05 0.0048747259 0.0048792052
#> groupd
#> (Intercept) -4.902954e-03
#> x 5.966093e-05
#> groupb 4.874726e-03
#> groupc 4.879205e-03
#> groupd 9.933362e-03The type="HC1" calculation estimates coefficient
uncertainty from the observed variability in contributions to the
estimating equations, with a finite-sample rescaling. This is useful
when the fitted mean model is appropriate but the response variance
differs from the family specification. This form of robust covariance is
also called a sandwich covariance.
When observations are correlated within groups, such as repeated
measurements on patients, fastglm_vcov() can sum those
contributions within each group before estimating uncertainty. The
groups, called clusters, are assumed independent of one another.
Reliable approximation generally requires sufficiently many independent
clusters. The calculation reads predictor rows in blocks and retains a
contribution vector for each cluster, so its memory requirement depends
on both block size and the number of clusters.
cluster <- rep(seq_len(100), each=5)
robust <- fastglm_vcov(fit, type="HC1", cluster=cluster)
fastglm_average(fit, vcov.=robust)
#> estimate std.error conf.low conf.high
#> 1 1.614 0.05835132 1.499634 1.728366fastglm_average() estimates an average response over a
supplied covariate distribution. Its standard error accounts for
coefficient uncertainty while conditioning on those covariates. A
contrast function can alter the design matrix, yielding the
average difference between transformed and original predictions. Such a
contrast describes the fitted mean model; a causal interpretation
requires additional identification assumptions.
A categorical predictor generally contributes several coefficients.
Testing those coefficients together addresses whether that predictor
contributes to the specified mean model, conditional on the remaining
terms. fastglm_wald_test() selects all columns associated
with a formula term and combines their estimates and covariance in a
joint Wald statistic.
The same interface supports ordinal, truncated-count, hurdle,
zero-inflated, and two-part models. For models with two components,
component selects the regression being tested; its default
tests the term in both components where present. The count and
semicontinuous-outcome vignettes illustrate these choices and their
interpretation.
fastglm_wald_test(fit,terms="group")
#> Joint Wald test
#> Terms: group
#> Estimation: fitted model coefficients
#> Covariance: model-based
#>
#> Statistic df p-value
#> 3.487 3 0.3225
#>
#> Tested contrasts:
#> Estimate Std. Error Null
#> groupb 0.16419 0.09947 0.00000
#> groupc 0.15478 0.10065 0.00000
#> groupd 0.14228 0.10240 0.00000
#>
#> Note: First-order chi-square approximation; validity depends on coefficient
#> estimability and the covariance assumptions.
fastglm_wald_test(fit,terms="group",vcov.=robust)
#> Joint Wald test
#> Terms: group
#> Estimation: fitted model coefficients
#> Covariance: supplied coefficient covariance
#>
#> Statistic df p-value
#> 3.599 3 0.3081
#>
#> Tested contrasts:
#> Estimate Std. Error Null
#> groupb 0.16419 0.09088 0.00000
#> groupc 0.15478 0.10315 0.00000
#> groupd 0.14228 0.10451 0.00000
#>
#> Note: First-order chi-square approximation; validity depends on coefficient
#> estimability and the covariance assumptions.The second test uses the cluster covariance calculated above. Both tests use a first-order chi-square reference; the covariance must correspond to the sampling structure. Small samples or few independent clusters can make this approximation inaccurate. A test of a main-effect term in a model containing interactions depends on the predictor coding and reference values. It does not test every effect involving that variable.
A named contrast specifies a linear combination of coefficients. For example, the difference between the coefficients for groups c and b describes their log mean ratio at equal values of the other predictors in this model.
group_comparison <- fastglm_wald_test(fit,contrast=c(groupc=1,groupb=-1))
group_comparison
#> Joint Wald test
#> Estimation: fitted model coefficients
#> Covariance: model-based
#>
#> Statistic df p-value
#> 0.009368 1 0.9229
#>
#> Tested contrasts:
#> Estimate Std. Error Null
#> Contrast 1 -0.009405 0.097163 0.000000
#>
#> Note: First-order chi-square approximation; validity depends on coefficient
#> estimability and the covariance assumptions.
as.data.frame(group_comparison)
#> statistic df p.value null.estimation available
#> 1 0.009368425 1 0.9228927 <NA> TRUEMatrix contrasts permit several comparisons in a joint test, and
null supplies their hypothesized values. Omitted named
coefficients receive weight zero. The returned result retains the
contrast estimates and standard errors at full precision. Nonestimable
comparisons, unsuccessful fits, and unusable covariance matrices produce
an unavailable test with an explanation.