---
title: "Models, inference, and computation with fastglm"
author: "Jared Huling"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Models, inference, and computation with fastglm}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse=TRUE, comment="#>", message=FALSE)
```

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.

# Package scope and model selection {#package-scope}

## Models and estimation

| 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](#specifying-the-mean-model) and [bias reduction](firth-fastglm.html) |
| Counts with variation beyond the Poisson model | `fastglm_nb()` | Maximum likelihood (ML) estimation of regression coefficients and negative-binomial size. [Count models](count-firth-fastglm.html#negative-binomial-variation) |
| Counts sampled conditional on being positive | `fastglm_truncated()` | Zero-truncated Poisson or negative-binomial likelihood. [Truncation and hurdles](count-firth-fastglm.html#positive-counts-and-hurdle-models) |
| Separate models for occurrence and positive counts | `fastglm_hurdle()` | Binary occurrence and zero-truncated Poisson or negative-binomial components, fitted by ML. [Hurdle models](count-firth-fastglm.html#positive-counts-and-hurdle-models) |
| Zeros arising from an additional component and a count distribution | `fastglm_zi()` | Zero-inflated Poisson or negative-binomial likelihood. [Zero-inflated models](count-firth-fastglm.html#zero-inflated-counts) |
| 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](semicontinuous-fastglm.html) |
| 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](ordinal-fastglm.html#ordinal-scope) |

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](firth-fastglm.html) lists the supported GLM combinations, while the ordinal and two-part guides describe their respective adjustments.

## Inference, prediction, and assessment

| 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](ordinal-fastglm.html#ordinal-scope) states these distinctions explicitly. The sections below explain [coefficient uncertainty and population averages](#inference-and-population-averages), [joint tests](#joint-coefficient-tests), and [readable result displays](#specifying-the-mean-model).

## Data size and repeated fitting

| 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](large-data-fastglm.html#sparse-and-mapped-matrices) |
| 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](large-data-fastglm.html) |
| 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](#repeated-outcomes-on-a-shared-design) |

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](large-data-fastglm.html) explains these requirements. [Benchmark examples](benchmarks-fastglm.html) examine speed together with numerical agreement.

The remaining sections demonstrate the shared workflow with a Poisson model: [specify the mean](#specifying-the-mean-model), [choose a numerical decomposition and assess convergence](#choosing-a-decomposition), [quantify uncertainty](#inference-and-population-averages), and [reuse a design](#repeated-outcomes-on-a-shared-design). Generalized linear model fitting uses iteratively reweighted least squares, which updates coefficients through a sequence of weighted regressions.

# Specifying the mean model

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.

```{r}
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)
predict(fit, newdata=d[1:4,], type="response", se.fit=TRUE)
```

For 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.

# Choosing a decomposition

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.

```{r}
fastglm_diagnostics(fit)
```

A 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.

# Inference and population averages

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.

```{r}
vcov(fit)
confint(fit, parm="x")
vcov(fit, type="HC1")
```

The `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.

```{r}
cluster <- rep(seq_len(100), each=5)
robust <- fastglm_vcov(fit, type="HC1", cluster=cluster)
fastglm_average(fit, vcov.=robust)
```

`fastglm_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.

## Joint coefficient tests

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.

```{r}
fastglm_wald_test(fit,terms="group")
fastglm_wald_test(fit,terms="group",vcov.=robust)
```

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.

```{r}
group_comparison <- fastglm_wald_test(fit,contrast=c(groupc=1,groupb=-1))
group_comparison
as.data.frame(group_comparison)
```

Matrix 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.

# Repeated outcomes on a shared design

Repeated outcome analyses often use the same predictors and weights. `fastglm_prepare()` saves this common design and the information needed to construct predictors consistently. For Gaussian regression with an identity link, the expensive matrix decomposition depends only on the predictors and weights, so it can be reused across outcomes. Certain bias-reduced models also permit this reuse, as illustrated below. For models whose fitting weights change with the response and fitted mean, the preceding fit can still provide useful starting coefficients. Reusing those starting coefficients is called a warm start.

```{r}
X <- model.matrix(~x+group, d)
Y <- cbind(outcome1=rnorm(nrow(d))+.4*d$x,
           outcome2=rnorm(nrow(d))-.2*d$x)
design <- fastglm_prepare(X)
models <- fastglm_batch(design, Y)
vapply(models, function(z) coef(z)["x"], numeric(1))
```

Bias-reduced models can be fitted to several outcomes by specifying `firth=TRUE`. This setting can be saved in the prepared design or supplied to `fastglm_batch()`. For Gaussian-identity, Gamma-log, and Poisson-square-root models, the weights in the regression updates are constant across iterations when the predictors and prior weights are fixed. The default automatic solver and explicit Cholesky method `"llt"` reuse the common matrix decomposition and row-influence calculations across responses. Dispersion and response-dependent coefficient adjustments are still estimated separately for each outcome.

```{r}
set.seed(148)
mu1 <- exp(.2+.3*d$x)
mu2 <- exp(.5-.2*d$x)
Y_positive <- cbind(outcome1=rgamma(nrow(d),shape=3,scale=mu1/3),
                    outcome2=rgamma(nrow(d),shape=3,scale=mu2/3))
design_firth <- fastglm_prepare(X,firth=TRUE)
models_firth <- fastglm_batch(design_firth,Y_positive,family=Gamma("log"))
vapply(models_firth,function(z) coef(z)["x"],numeric(1))
vapply(models_firth,function(z) z$dispersion,numeric(1))
```

Changing the predictors or weights rebuilds the saved calculations. Offsets and outcomes may differ while the common matrix decomposition is retained because they do not determine that decomposition for these family-link combinations. Other supported Firth models also work through the batch interface, using their ordinary fitting calculations. Explicit QR, SVD, and LDLT choices retain their requested solver. The amount of time saved therefore depends on the model, the number of responses, and the cost of the shared setup.

`workers` sets the number of R processes used to fit outcomes concurrently and defaults to one. Using more processes can shorten elapsed time when individual fits are expensive enough to offset process startup and data-transfer costs. Each process needs its own working memory, and concurrent fits cannot use the immediately preceding result as a starting value. These costs should be included when comparing runtimes.

When only coefficients and fitted means are required, `inference="none"` skips the final covariance calculation. `inference="diagonal"` retains standard errors without retaining the full covariance; it does not avoid computing that covariance. `retain.data=FALSE` omits the design and response from the returned object. Subsequent procedures that require omitted information must be supplied with new data or require a refit.

```{r}
compact <- fastglm(X, d$y, poisson(), method="auto",
                   inference="none", retain.data=FALSE)
coef(compact)
```