Package {fastglm}

Type:	Package
Title:	Fast and Stable Fitting of Generalized Linear Models using 'RcppEigen'
Version:	0.1.3
Description:	Fits generalized linear models efficiently using 'RcppEigen'. The iteratively reweighted least squares implementation utilizes the step-halving approach of Marschner (2011) <doi:10.32614/RJ-2011-012> to help safeguard against convergence issues. Provides formula and matrix interfaces for count and cumulative-link ordinal models, mean bias reduction, and fitting from data read in blocks of rows for supported model classes.
BugReports:	https://github.com/jaredhuling/fastglm/issues
URL:	https://github.com/jaredhuling/fastglm, https://jaredhuling.org/fastglm/
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
Encoding:	UTF-8
Depends:	R (≥ 4.1.0)
Imports:	Rcpp (≥ 0.12.13), bigmemory, methods, Matrix, Formula
LinkingTo:	Rcpp, RcppEigen, BH, bigmemory
Suggests:	glm2, knitr, rmarkdown, testthat (≥ 3.1.7), sandwich, microbenchmark, MASS, statmod, logistf, brglm2, pscl, speedglm, tweedie, ggplot2, broom, DBI, RSQLite, numDeriv, arrow
Config/testthat/edition:	3
RoxygenNote:	7.3.3
VignetteBuilder:	knitr
NeedsCompilation:	yes
Packaged:	2026-11-08 01:09:43 UTC; huling
Author:	Jared Huling [aut, cre], Douglas Bates [cph], Dirk Eddelbuettel [cph], Romain Francois [cph], Yixuan Qiu [cph], Noah Greifer [ctb]
Maintainer:	Jared Huling <jaredhuling@gmail.com>
Repository:	CRAN
Date/Publication:	2026-11-09 10:50:09 UTC

big.matrix prod

Description

big.matrix prod

big.matrix prod

Usage

## S4 method for signature 'big.matrix,vector'
x %*% y

## S4 method for signature 'vector,big.matrix'
x %*% y

Arguments

Confidence intervals for GLM coefficients

Description

Returns Wald intervals by default. Binomial-logit Firth fits can instead use penalized profile likelihood intervals; see [fastglm_profile_test()] for assumptions, numerical controls, and diagnostics.

Usage

## S3 method for class 'fastglm'
confint(
  object,
  parm,
  level = 0.95,
  ...,
  method = c("wald", "profile"),
  control = list()
)

Arguments

Value

A matrix with lower and upper confidence limits. Profile intervals additionally carry convergence and endpoint diagnostics as attributes.

Fast generalized linear model fitting

Description

Fast generalized linear model fitting

bigLm default

Usage

fastglm(x, ...)

## Default S3 method:
fastglm(
  x,
  y,
  family = gaussian(),
  weights = NULL,
  offset = NULL,
  start = NULL,
  etastart = NULL,
  mustart = NULL,
  method = 0L,
  tol = 1e-08,
  maxit = 100L,
  firth = FALSE,
  inference = c("full", "diagonal", "none"),
  retain.data = TRUE,
  ...
)

Arguments

Value

A list with the elements

See Also

[fastglm_fit()]

Examples

x <- matrix(rnorm(10000 * 100), ncol = 100)
y <- 1 * (0.25 * x[,1] - 0.25 * x[,3] > rnorm(10000))

system.time(gl1 <- glm.fit(x, y, family = binomial()))

system.time(gf1 <- fastglm(x, y, family = binomial()))

system.time(gf2 <- fastglm(x, y, family = binomial(), method = 1))

system.time(gf3 <- fastglm(x, y, family = binomial(), method = 2))

system.time(gf4 <- fastglm(x, y, family = binomial(), method = 3))

max(abs(coef(gl1) - gf1$coef))
max(abs(coef(gl1) - gf2$coef))
max(abs(coef(gl1) - gf3$coef))
max(abs(coef(gl1) - gf4$coef))


## Not run: 
nrows <- 50000
ncols <- 50
bkFile <- "bigmat2.bk"
descFile <- "bigmatk2.desc"
bigmat <- filebacked.big.matrix(nrow=nrows, ncol=ncols, type="double",
                                backingfile=bkFile, backingpath=".",
                                descriptorfile=descFile,
                                dimnames=c(NULL,NULL))
for (i in 1:ncols) bigmat[,i] = rnorm(nrows)*i
y <- 1*(rnorm(nrows) + bigmat[,1] > 0)

system.time(gfb1 <- fastglm(bigmat, y, family = binomial(), method = 3))

## End(Not run)

Display profile confidence intervals

Description

Profile intervals print a compact table while retaining endpoint and fit diagnostics as attributes. Matrix indexing and numeric extraction are unchanged. Failed limits remain missing rather than being replaced by Wald limits.

Usage

## S3 method for class 'fastglm_confint'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

Arguments

Value

Invisibly returns x with all numerical values and attributes intact.

Display and extract numerical diagnostics

Description

Diagnostics use the same headings and numeric formatting as fitted models. Scores and remaining updates assess numerical convergence; they are not hypothesis tests. For bias-reduced GLMs the ordinary likelihood score can remain nonzero even when the adjusted equations have been solved.

Usage

## S3 method for class 'fastglm_diagnostics'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

## S3 method for class 'fastglm_diagnostics'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

Arguments

Value

The print method invisibly returns x. as.data.frame() returns one row with numeric and logical scalar columns at their original precision. Coefficient score vectors, aliased coefficients, collapsed categories, and rejected-evaluation details are retained as list columns.

Display and extract predictions with standard errors

Description

Predictions requested with se.fit=TRUE retain their fit and se.fit components and print aligned estimates and standard errors. For matrix predictions, each row and component (such as an outcome category) occupies one row of the displayed or extracted table. Missing standard errors remain missing, including for collapsed ordinal categories.

Usage

## S3 method for class 'fastglm_prediction'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

## S3 method for class 'fastglm_prediction'
print(x, digits = max(3L, getOption("digits") - 3L), max.rows = 20L, ...)

Arguments

Value

The print method invisibly returns x. as.data.frame() returns row identifiers, component labels for matrix predictions, and numeric fit and se.fit columns at full precision. Original matrix dimensions, names, and any model-specific metadata remain in x.

Empirical estimating-function and bread methods for 'fastglm' objects

Description

Methods for 'sandwich::estfun()' and 'sandwich::bread()', which let 'sandwich::vcovCL()', 'sandwich::vcovBS()', and the rest of the *sandwich* machinery work directly on 'fastglm' and 'fastglmFit' objects. Load 'sandwich' ('library(sandwich)') before calling them.

Usage

estfun.fastglm(x, ...)

estfun.fastglmFit(x, ...)

bread.fastglm(x, ...)

bread.fastglmFit(x, ...)

Arguments

Details

'estfun(x)' returns the per-observation empirical estimating-function contributions, an 'n x p' matrix whose 'i'-th row is '(y_i - mu_i) * mu.eta_i / variance(mu_i) * x_i / dispersion'. 'bread(x)' returns ‘(X’ W X)^(-1) * n * dispersion' (the *sandwich* convention). For poisson, binomial, and negative-binomial families the dispersion is fixed at 1, exactly matching 'sandwich::estfun.glm()' / 'sandwich::bread.glm()'.

Value

'estfun()' returns an 'n x p' matrix; 'bread()' returns a 'p x p' matrix.

Display and extract hypothesis tests

Description

Test objects print a compact table while retaining full-precision numeric results. Data-frame and optional broom methods provide a single row for reporting, including unavailable tests and their explanations.

Usage

## S3 method for class 'fastglm_test'
as.data.frame(x, row.names = NULL, optional = FALSE, details = FALSE, ...)

## S3 method for class 'fastglm_test'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

tidy.fastglm_test(x, ...)

## S3 method for class 'fastglm_test_list'
as.data.frame(x, row.names = NULL, optional = FALSE, details = FALSE, ...)

## S3 method for class 'fastglm_test_list'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

tidy.fastglm_test_list(x, ...)

Arguments

Value

Print methods invisibly return the original object. Extraction methods return a data frame with numeric statistics and p-values.

Tidy summaries for fastglm models

Description

Optional methods for the broom generics. Formula fits also retain the standard glm model structure for tools that consume coefficients, terms, prediction, and covariance methods.

Usage

tidy.fastglm(x, conf.int = FALSE, conf.level = 0.95, exponentiate = FALSE, ...)

tidy.fastglmFit(
  x,
  conf.int = FALSE,
  conf.level = 0.95,
  exponentiate = FALSE,
  ...
)

glance.fastglm(x, ...)

glance.fastglmFit(x, ...)

Arguments

Value

A data frame with one row per coefficient ('tidy') or one row per model ('glance').

Fit a generalized linear model from a formula

Description

Constructs the model frame using the same conventions as [stats::glm()]. Factors, contrasts, offsets, missing observations, and prediction terms are retained for subsequent inference and prediction.

Usage

## S3 method for class 'formula'
fastglm(
  x,
  data,
  family = gaussian(),
  subset,
  na.action,
  weights,
  offset,
  contrasts = NULL,
  method = 0L,
  tol = 1e-08,
  maxit = 100L,
  firth = FALSE,
  ...
)

Arguments

Value

A fitted model inheriting from 'fastglmFit', 'glm', and 'lm'.

Fast generalized linear model fitting

Description

Fast generalized linear model fitting

Usage

fastglmPure(
  x,
  y,
  family = gaussian(),
  weights = rep(1, NROW(y)),
  offset = rep(0, NROW(y)),
  start = NULL,
  etastart = NULL,
  mustart = NULL,
  method = 0L,
  tol = 1e-07,
  maxit = 100L,
  firth = FALSE,
  inference = c("full", "diagonal", "none")
)

Arguments

Value

A list with the elements

See Also

[fastglm_fit()]

Examples

set.seed(1)
x <- matrix(rnorm(1000 * 25), ncol = 25)
eta <- 0.1 + 0.25 * x[,1] - 0.25 * x[,3] + 0.75 * x[,5] -0.35 * x[,6] #0.25 * x[,1] - 0.25 * x[,3]
y <- 1 * (eta > rnorm(1000))

yp <- rpois(1000, eta ^ 2)
yg <- rgamma(1000, exp(eta) * 1.75, 1.75)

# binomial
system.time(gl1 <- glm.fit(x, y, family = binomial()))

system.time(gf1 <- fastglmPure(x, y, family = binomial(), tol = 1e-8))

system.time(gf2 <- fastglmPure(x, y, family = binomial(), method = 1, tol = 1e-8))

system.time(gf3 <- fastglmPure(x, y, family = binomial(), method = 2, tol = 1e-8))

system.time(gf4 <- fastglmPure(x, y, family = binomial(), method = 3, tol = 1e-8))

max(abs(coef(gl1) - gf1$coef))
max(abs(coef(gl1) - gf2$coef))
max(abs(coef(gl1) - gf3$coef))
max(abs(coef(gl1) - gf4$coef))

# poisson
system.time(gl1 <- glm.fit(x, yp, family = poisson(link = "log")))

system.time(gf1 <- fastglmPure(x, yp, family = poisson(link = "log"), tol = 1e-8))

system.time(gf2 <- fastglmPure(x, yp, family = poisson(link = "log"), method = 1, tol = 1e-8))

system.time(gf3 <- fastglmPure(x, yp, family = poisson(link = "log"), method = 2, tol = 1e-8))

system.time(gf4 <- fastglmPure(x, yp, family = poisson(link = "log"), method = 3, tol = 1e-8))

max(abs(coef(gl1) - gf1$coef))
max(abs(coef(gl1) - gf2$coef))
max(abs(coef(gl1) - gf3$coef))
max(abs(coef(gl1) - gf4$coef))

# gamma
system.time(gl1 <- glm.fit(x, yg, family = Gamma(link = "log")))

system.time(gf1 <- fastglmPure(x, yg, family = Gamma(link = "log"), tol = 1e-8))

system.time(gf2 <- fastglmPure(x, yg, family = Gamma(link = "log"), method = 1, tol = 1e-8))

system.time(gf3 <- fastglmPure(x, yg, family = Gamma(link = "log"), method = 2, tol = 1e-8))

system.time(gf4 <- fastglmPure(x, yg, family = Gamma(link = "log"), method = 3, tol = 1e-8))

max(abs(coef(gl1) - gf1$coef))
max(abs(coef(gl1) - gf2$coef))
max(abs(coef(gl1) - gf3$coef))
max(abs(coef(gl1) - gf4$coef))

Average predictions and contrasts with delta-method uncertainty

Description

Computes averages without retaining individual predictions. For a streaming fit, a contrast function can transform each design chunk before prediction. Uncertainty conditions on the supplied covariate distribution. For ordinal fits, returns an average for each category or cumulative probability. With collapsed categories, inference refers to the reduced category model and removed-category standard errors are unavailable. Nonparallel ordinal fits include the joint uncertainty of thresholds, common slopes, and threshold-specific slopes and check prediction validity at all rows with positive averaging weight.

Usage

fastglm_average(
  object,
  newdata = NULL,
  contrast = NULL,
  weights = NULL,
  vcov. = NULL,
  level = 0.95,
  type = "response"
)

Arguments

Value

A data frame with estimate, standard error, and confidence limits.

Fit multiple responses using a prepared design

Description

The same design can be used for many outcomes. Gaussian identity fits share the weighted crossproduct factor. With 'firth = TRUE', Gaussian identity, Gamma log, and Poisson square-root fits share their factor and observation leverages when using automatic or LLT fitting. Changes to predictors or weights rebuild this setup; offsets and responses can vary independently. Other families and decompositions use the ordinary fitting routine. Sequential fits may use the preceding converged fit as a starting value. Independent fits can run in separate R processes.

Usage

fastglm_batch(
  design,
  y,
  family = gaussian(),
  weights = NULL,
  offset = NULL,
  method = "auto",
  tol = 1e-08,
  maxit = 100L,
  warm.start = TRUE,
  workers = 1L,
  inference = c("full", "diagonal", "none"),
  retain.data = TRUE,
  firth = NULL
)

Arguments

Value

A named list of fitted models.

Refit bootstrap uncertainty for GLMs, count models, and ordinal models

Description

Repeatedly resamples observations or independent clusters, or generates new responses from the fitted model, and refits the regression. A statistic such as a coefficient or average predicted response is then evaluated after each refit. This includes uncertainty from estimating both model components and negative-binomial size when applicable.

Usage

fastglm_bootstrap(
  object,
  R = 199L,
  statistic = NULL,
  type = c("pairs", "cluster", "parametric"),
  cluster = NULL,
  seed = NULL,
  workers = 1L,
  level = 0.95,
  control = list()
)

## S3 method for class 'fastglm_bootstrap'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

Arguments

Details

Pair resampling treats rows as independent sampling units. Cluster resampling treats the supplied clusters as independent and retains all rows from each sampled cluster. The original observation weights and offsets are retained on every sampled copy. Parametric resampling assumes the fitted response distributions; it does not model unspecified within-cluster dependence. For a numeric binomial response, positive prior weights are interpreted as integer trial counts, as in 'stats::simulate()'; for a two-column response, the column sum supplies the trial count and any additional case weights are held fixed. The formula design is held fixed, including factor levels and spline bases, so data-driven preprocessing or model selection is not repeated. These procedures require regular, identifiable fits; they do not establish valid inference at separation or mixture boundaries. Streaming parametric refits regenerate responses one block at a time with fixed block seeds. Pair and cluster resampling also retain row indices or multiplicities, requiring memory proportional to the number of rows; neither procedure collects the full streaming design matrix.

Ordinal bootstrap support currently requires common slopes. For ordinal fits, pair and cluster resampling copy stored rows, including their frequency weights and grouped counts; they do not resample the individuals represented by a frequency row separately. Parametric resampling draws category counts with each row's effective integer total, then refits those counts with unit weights. The declared outcome levels, design columns, and offsets are preserved. Formula metadata is retained for statistics that predict on new data. Ordinal optimizer arguments may be changed through 'control', but the estimator, design, response, weights, and category order may not be overridden there. Stored zero-frequency rows remain sampling units in pairs and cluster resampling. Ordinal cluster identifiers may match retained rows or rows before missing-value omission, after any subset; every retained sampling row must have a cluster identifier. Ordinal results record 'rng.kind' alongside replicate seeds for reproducing the random-number stream.

Ordinal point replicates and category diagnostics are retained when a refit has empty categories or fails to converge, provided the statistic can be evaluated. For mean bias reduction with 'type = "parametric"', empty endpoint categories are allowed in the original fit and its refits because the fitted thresholds remain finite. The original fit and every refit must retain all interior categories, and every refit and statistic must succeed. An empty interior category ties adjacent thresholds and changes the fitted model, so uncertainty summaries remain unavailable. The same empty-category restriction applies to maximum likelihood and to pairs or cluster resampling. These rules also apply to custom statistics; point replicates remain available for inspection.

Every failed or nonconverged refit is recorded. If any replicate fails, standard errors and intervals are returned as missing rather than silently conditioning inference on the successful fits. Inspect 'diagnostics' and resolve the fitting problem before interpreting bootstrap uncertainty.

Value

A 'fastglm_bootstrap' list with the original 'estimate', a matrix of 'replicates', 'std.error', percentile 'conf.int', per-refit 'diagnostics', replicate 'seeds', and resampling details.

Profile-likelihood interval for negative-binomial size

Description

Refits regression coefficients at fixed theta and inverts the likelihood-ratio test on the log-theta scale. The interval may have an infinite upper endpoint when the Poisson boundary cannot be excluded. Bias-reduced fits are excluded because their estimates do not maximize this likelihood.

Usage

fastglm_confint_theta(object, level = 0.95, tol = 1e-05, maxit = 200L)

Arguments

Value

A named interval for theta with the evaluated profile attached.

Numerical diagnostics for a fitted GLM

Description

Reports whether the iteration converged, which coefficients can be estimated separately, and whether fitted parameters are close to a boundary. The score is the derivative of the log likelihood with respect to the coefficients; a large value after accounting for predictor scale can indicate incomplete convergence. Covariance conditioning describes how strongly nearly redundant coefficient effects can amplify numerical error. For bias-reduced fits the ordinary score need not vanish. The stored adjusted coefficient and dispersion scores instead assess the equations used to fit the model at the reported estimates.

Usage

fastglm_diagnostics(object, boundary = FALSE, max.rows = 5L)

Arguments

Value

A fastglm_diagnostics list with a compact printed summary. Existing components remain accessible with $ or [[ at full precision. as.data.frame() returns one row, with list columns for coefficient scores, category labels, and other nonscalar details. score omits division by the dispersion, and standardized.score uses the correspondingly unscaled information diagonal, retaining the ordinary GLM score convention. For a bias-reduced fit, adjusted.score contains the adjusted likelihood scores, including dispersion scaling, and standardized.adjusted.score is the largest absolute adjusted score divided by the square root of its expected information diagonal. adjusted.dispersion.score and standardized.dispersion.score describe the dispersion equation; they are missing for families whose dispersion is fixed. convergence.reason records the solver's stopping reason when available. Ordinary scores require retained data or a readable source. Stored adjusted scores remain available after data are omitted. With inference = "none", adjusted scores are retained only when the existing constant-information calculation supplies the required leverages, including prepared Gamma-log fits. Other compact fits and older objects without stored adjusted scores report them as unavailable rather than computing an additional covariance matrix. QR source fits also report their remaining predictor change, joint information-scaled score per positive-weight row, and dispersion checks. For binomial probit and complementary log-log bias reduction, step.reductions counts scoring updates reduced to limit their largest fitted linear-predictor change to one. A reduced step alone cannot establish convergence. These safeguards do not guarantee convergence for every design. For these bias-reduced binomial fits, scores and information use stable probability-tail calculations; rounded fitted probabilities are not used to reconstruct the estimating equations. Ordinal diagnostics report scores in the effective category model after any empty interior categories have been collapsed; their standardized score is calculated in the centered/scaled fitting coordinates. raw.score.maximum is the largest absolute score in the effective model's reported predictor units. Its magnitude therefore depends on predictor scaling. With boundary=TRUE, boundary.details identifies the smallest fitted probabilities and gaps, their row and category labels, observed frequencies, offsets, and original predictor values. A displayed zero probability can be numerical underflow; its log probability remains informative. These observations help diagnose a fit but do not prove that an interior solution is absent. Source rows and values must remain unchanged since fitting. Unavailable row diagnostics retain a reason.

Fast generalized linear model fitting

Description

'fastglm_fit()' is a fitting method for [glm()]. It works like 'glm.fit()', i.e., by being supplied to the 'method' argument of 'glm()'.

Usage

fastglm_fit(
  x,
  y,
  weights = rep(1, NROW(y)),
  start = NULL,
  etastart = NULL,
  mustart = NULL,
  offset = rep(0, NROW(y)),
  family = gaussian(),
  control = list(),
  intercept = TRUE,
  singular.ok = TRUE,
  firth = FALSE
)

fastglm_control(fastmethod = 0L, tol = 1e-07, maxit = 100L, firth = FALSE, ...)

## S3 method for class 'fastglmFit'
vcov(object, ...)

## S3 method for class 'fastglmFit'
summary(object, ...)

Arguments

Details

The purpose of the functions documented on this page is to facilitate integration with existing [glm()] utilities in base R. 'fastglm_fit()' is just a wrapper for [fastglmPure()] with some additional quality-of-life features. The 'vcov()' and 'summary()' methods use the unscaled coefficient covariance matrix returned directly from the C++ solver, so no refit is required.

Examples

set.seed(1234)
n <- 1e4
x <- matrix(rnorm(n * 25), ncol = 25)
eta <- 0.1 + 0.25 * x[,1] - 0.25 * x[,3] + 0.75 * x[,5] -0.35 * x[,6]
dat <- as.data.frame(x)

# binomial
dat$y <- rbinom(n, 1, pnorm(eta))

system.time({
    gl <- glm(y ~ ., data = dat,
              family = binomial)
})

system.time({
    gf0 <- glm(y ~ ., data = dat,
               family = binomial,
               method = fastglm_fit)
})

system.time({
    gf1 <- glm(y ~ ., data = dat,
               family = binomial,
               method = fastglm_fit,
               fastmethod = 1)
})

# poisson
dat$y <- rpois(n, eta^2)

system.time({
    gl <- glm(y ~ ., data = dat,
              family = poisson)
})

system.time({
    gf0 <- glm(y ~ ., data = dat,
               family = poisson,
               method = fastglm_fit)
})

system.time({
    gf1 <- glm(y ~ ., data = dat,
               family = poisson,
               method = fastglm_fit,
               fastmethod = 1)
})

# gamma
dat$y <- rgamma(n, exp(eta) * 1.75, 1.75)

system.time({
    gl <- glm(y ~ ., data = dat,
              family = Gamma(link = "log"))
})

system.time({
    gf0 <- glm(y ~ ., data = dat,
               family = Gamma(link = "log"),
               method = fastglm_fit)
})

system.time({
    gf1 <- glm(y ~ ., data = dat,
               family = Gamma(link = "log"),
               method = fastglm_fit,
               fastmethod = 1)
})

# Different (equivalent) ways of supplying
# control arguments:
gf1 <- glm(y ~ ., data = dat,
           family = Gamma(link = "log"),
           method = fastglm_fit,
           fastmethod = 1)

gf1 <- glm(y ~ ., data = dat,
           family = Gamma(link = "log"),
           method = fastglm_fit,
           control = list(fastmethod = 1))

gf1 <- glm(y ~ ., data = dat,
           family = Gamma(link = "log"),
           method = fastglm_fit,
           control = list(method = 1))

Hurdle (two-part) Poisson / Negative-Binomial regression

Description

Fits a hurdle model for count data with excess zeros: a binary (zero / non-zero) regression on the full sample combined with a zero- truncated count regression on the positive subset. The two parts have independent likelihoods and are estimated jointly. With 'dist = "negbin"' the negative-binomial dispersion 'theta' is estimated by an inner Brent-1D MLE inside the C++ driver, alternating with the count-side IRLS.

Usage

fastglm_hurdle(
  formula,
  data,
  subset,
  na.action,
  weights,
  offset,
  dist = c("poisson", "negbin"),
  zero.dist = "binomial",
  link = c("logit", "probit", "cloglog", "log"),
  init.theta = NULL,
  tol = 1e-08,
  maxit = 100L,
  outer.tol = 1e-07,
  outer.maxit = 50L,
  theta.tol = 1e-08,
  theta.maxit = 100L,
  model = TRUE,
  y = TRUE,
  x = FALSE,
  response = NULL,
  z = NULL,
  method = 2L,
  block.size = 8192L
)

Arguments

Details

Formula fits use the native count-model driver. Sparse designs and data supplied in blocks use the same likelihood with blockwise score updates. The latter permits repeated passes over data that need not fit in memory.

Value

A list of class 'c("fastglm_hurdle", "fastglm")' with elements 'coefficients' (a list with 'count' and 'zero'), 'se' (likewise), 'vcov' (block-diagonal), 'loglik', 'loglik_count', 'loglik_zero', 'theta', 'SE.theta', 'iter_*', and 'converged'.

Examples

set.seed(1)
n <- 300
x1 <- rnorm(n); x2 <- rnorm(n)
eta_count <- 1.0 + 0.4*x1 - 0.2*x2
lam <- exp(eta_count)
y_count <- rpois(n, lam)
eta_zero  <- -0.5 + 0.6*x1
p_pos     <- plogis(eta_zero)
is_pos    <- rbinom(n, 1, p_pos)
# zero-truncated: resample any rare zero from positive Poisson
for (i in seq_len(n)) {
  while (y_count[i] == 0) y_count[i] <- rpois(1, lam[i])
}
y <- ifelse(is_pos == 1, y_count, 0)
df <- data.frame(y = y, x1 = x1, x2 = x2)
fit <- fastglm_hurdle(y ~ x1 + x2, data = df, dist = "poisson")
coef(fit, model = "count")
coef(fit, model = "zero")

Average marginal effects for ordinal probabilities

Description

Estimates the average rate of change in each category probability as a continuous predictor changes. Uncertainty uses the delta method, conditional on the supplied predictor values and averaging weights.

Usage

fastglm_marginal_effects(
  object,
  variable,
  newdata = NULL,
  weights = NULL,
  vcov. = NULL,
  level = 0.95,
  type = "response",
  eps = NULL,
  offset = NULL
)

Arguments

Details

A marginal effect is a derivative, measured in probability units per unit of the continuous predictor. It describes a local rate of change; a finite change can give a different probability difference. Use fastglm_average(contrast = ...) for discrete comparisons, including factors. A numeric binary indicator is treated as continuous if requested.

Formula design derivatives use centered differences, with the fitted prediction basis and factor coding retained. Terms must be smooth functions of each row's predictors, or use a prediction basis frozen during fitting (for example, poly() or splines::ns()). Expressions that recompute sample summaries such as I(x - mean(x)) are unsupported. Rows and perturbed values with positive averaging weight must be finite and within the domain of the formula transformations. A smaller eps may be needed near a transformation's domain boundary.

Derivatives of probabilities and their parameter gradients are analytic conditional on the computed design derivative. The gradient includes uncertainty in both thresholds and slopes. Computation uses blocks of rows and retains the base prediction design, without storing a full array of observation-by-parameter-by-category derivatives. Normal intervals are first-order approximations; the covariance and reduced-category conventions match fastglm_average(). For nonparallel models, each cumulative probability uses its own slope. Evaluation rows must lie within the fitted model's valid predictor domain. Supplied covariance, including cluster covariance, retains covariance among thresholds, common slopes, and threshold-specific slopes.

Value

A data frame with variable, category, estimate, std.error, conf.low, and conf.high. Category effects sum to zero. Removed categories have zero effects and unavailable standard errors; inference otherwise refers to the reduced category model.

Examples

set.seed(91)
d <- data.frame(x=rnorm(150))
d$y <- ordered(cut(.5*d$x + rlogis(150), c(-Inf,-.5,.7,Inf)))
fit <- fastglm_ordinal(y ~ x, data=d, firth=TRUE)
fastglm_marginal_effects(fit, "x")

Fit a negative-binomial GLM with simultaneous (beta, theta) MLE

Description

'fastglm_nb()' fits a negative-binomial regression model, jointly maximising the likelihood over the regression coefficients 'beta' and the NB2 dispersion 'theta'. It is the fastglm analogue of [MASS::glm.nb()], Dense fits use native coefficient and size updates. Sparse and streamed fits alternate the existing GLM solver with size updates accumulated over row blocks, without loading the full design into a dense matrix.

Usage

fastglm_nb(
  x,
  y = NULL,
  weights = NULL,
  offset = NULL,
  start = NULL,
  init.theta = NULL,
  link = c("log", "sqrt", "identity"),
  method = 2L,
  tol = 1e-08,
  maxit = 100L,
  outer.maxit = 50L,
  outer.tol = 1e-07,
  theta.tol = 1e-08,
  theta.maxit = 100L,
  block.size = 8192L
)

Arguments

Value

A list of class 'c("fastglm_nb", "fastglm")' with the usual fastglm components plus 'theta', 'SE.theta', 'iter.theta', and 'twologlik' (twice the maximized NB log-likelihood). 'SE.theta' uses the joint observed information for coefficients and theta. The default coefficient covariance retains the GLM expected-information convention; 'vcov(fit, full=TRUE)' returns the joint observed covariance.

Examples

set.seed(1)
n <- 500
x <- cbind(1, matrix(rnorm(n * 2), n, 2))
eta <- x %*% c(0.3, 0.5, -0.2)
mu  <- exp(eta)
if (requireNamespace("MASS", quietly = TRUE)) {
  y   <- MASS::rnegbin(n, mu = mu, theta = 2)
  fit <- fastglm_nb(x, y)
  c(theta = fit$theta, fit$coefficients)
}

Fit a cumulative-link model for ordered outcomes

Description

Fits a cumulative-link model with common slopes and freely estimated thresholds, by maximum likelihood or mean bias reduction. Selected terms can have threshold-specific slopes in a dense fit with either estimation method. Positive slopes shift probability toward higher outcome categories.

Usage

## S3 method for class 'fastglm_ordinal_design'
fastglm_ordinal(
  x,
  y,
  weights = NULL,
  offset = NULL,
  levels = NULL,
  firth = NULL,
  tol = 1e-08,
  maxit = 200L,
  block.size = 2048L,
  control = list(),
  link = NULL,
  nonparallel = NULL,
  ...
)

fastglm_ordinal(x, ...)

## Default S3 method:
fastglm_ordinal(
  x,
  y,
  weights = NULL,
  offset = NULL,
  levels = NULL,
  firth = FALSE,
  tol = 1e-08,
  maxit = 200L,
  block.size = 2048L,
  control = list(),
  link = "logit",
  nonparallel = NULL,
  ...
)

## S3 method for class 'formula'
fastglm_ordinal(
  x,
  data,
  weights,
  offset,
  subset,
  na.action,
  contrasts = NULL,
  levels = NULL,
  firth = FALSE,
  tol = 1e-08,
  maxit = 200L,
  block.size = 2048L,
  control = list(),
  model = TRUE,
  link = "logit",
  nonparallel = NULL,
  ...
)

Arguments

Details

By default, thresholds are estimated in increasing order. The same location slopes apply to all cumulative probabilities. Full-rank location designs are required. Sparse matrices and repeatable sources are not accepted by this interface.

With nonparallel, the model is ⁠link{Pr(Y <= j)} = alpha[j] - x_common beta - x_nonparallel gamma[,j] - offset⁠. A selected predictor has a separate slope at each threshold and no additional common slope. Fitting requires every declared category to have positive frequency. The cumulative predictors must remain ordered at the observed predictor values, and predictions check this condition at each supplied row. Ordering at training values does not establish validity under extrapolation. A fully varying logit mean-bias-reduced model that is saturated across predictor patterns, with constant offsets within patterns and an empty intermediate category in a pattern, has no interior adjusted-score solution and is rejected before iteration. Other boundary failures retain their convergence diagnostics. See ordinal-partial-methods.

In common-slope models, an empty interior category is removed from the numerical fit and then restored to the declared scale with zero fitted probability and coincident thresholds. Convergence and information refer to the reduced model. Standard errors for the tied thresholds are unavailable. Bias-reduced fits retain empty endpoint categories; maximum likelihood with empty endpoints requires infinite thresholds and is rejected.

The covariance uses inverse expected Fisher information at the fitted estimate. It gives a first-order model-based approximation; this interface does not provide ordinal penalized-profile inference. Robust covariance is available through fastglm_vcov() or vcov(type="HC0") for converged retained dense or streamed fits with every declared category observed. Mean bias reduction solves the adjusted equations for the original thresholds and slopes. It does not maximize a log-determinant penalized ordinal likelihood.

logLik() evaluates the ordinary categorical frequency likelihood, omitting multinomial grouping constants. Consequently, grouping identical covariate patterns preserves this likelihood convention. For bias-reduced fits it is the ordinary likelihood evaluated at the bias-reduced estimate. deviance() compares this kernel with a saturated fit to the supplied response rows. It therefore uses a grouped multinomial saturated model for a count-matrix response and depends on the chosen grouping.

Value

An object of class fastglm_ordinal, with thresholds alpha, location coefficients beta, the full coefficients vector, category mapping, covariance, convergence diagnostics, and retained design. With selected nonparallel effects, the additional class fastglm_ordinal_partial provides threshold-specific slope methods; gamma stores these slopes as a predictor-by-threshold matrix, and beta contains only the common slopes. nobs() gives the total frequency, while nobs.rows counts rows with positive total frequency.

References

Kosmidis, I. (2014). Improved estimation in cumulative link models. Journal of the Royal Statistical Society, Series B, 76, 169–196. doi:10.1111/rssb.12025.

Examples

set.seed(91)
x <- rnorm(150)
y <- ordered(cut(.5*x + rlogis(150), c(-Inf,-.5,.7,Inf),
                 labels=c("low","middle","high")))
fit <- fastglm_ordinal(y ~ x, firth=TRUE)
coef(fit, component="location")
head(predict(fit))

Assess ordinal prediction and the common-slope assumption

Description

Compares predicted cumulative probabilities with observed frequencies, scores the full ordered prediction distribution, and optionally compares predictor slopes across binary divisions of the outcome scale or tests the common-slope restriction using an analytic score statistic.

Usage

fastglm_ordinal_assess(
  object,
  newdata = NULL,
  y = NULL,
  weights = NULL,
  offset = NULL,
  bins = 10L,
  proportional.odds = FALSE,
  po.control = list(),
  po.test = c("none", "score"),
  po.terms = NULL,
  po.test.control = list(),
  po.estimation = c("ML", "meanBR"),
  keep.rows = NULL,
  sink = NULL
)

## S3 method for class 'fastglm_ordinal_assessment'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'fastglm_ordinal_assessment'
plot(x, type = c("calibration", "proportional.odds"), level = 0.95, ...)

Arguments

Details

The ranked probability score (rps) is \sum_{j=1}^{K-1}\{F_j-I(Y\leq j)\}^2; rps_normalized divides by K-1, and lies between zero and one. log_loss is -\log p_Y. All are losses: smaller values are better. Grouped rows average the individual-outcome losses, rather than squaring differences from grouped outcome proportions. Overall scores average over the included individual frequency. Log probabilities are evaluated directly without clipping. A positive count in a category with structural probability zero has infinite log loss; a zero count contributes zero. Extremely small positive probabilities can have a finite log loss even when the displayed probability underflows to zero.

Missing outcomes, weights, predictors, or offsets, and zero-frequency rows, are retained in rows with missing scores and an exclusion reason. Infinite predictors or offsets are errors. Empty calibration bins remain in the table with zero frequency and missing means. With no usable outcomes, overall scores are missing. Bins include their lower boundary; the last bin also includes one. Fitting rows omitted under na.exclude are restored when assessing the original response. Explicit responses with no newdata must match the retained design rows.

Scores on the fitting outcomes describe apparent performance. Supplying held-out outcomes evaluates the frozen model; the function does not establish that those outcomes were excluded during fitting. The optional binary fits use the supplied assessment outcomes and design when supplied. They affect only the proportional-odds display, never predictions or scores. The dashed reference slopes always come from the original ordinal model. Positive slopes correspond to higher outcomes, matching the ordinal model. The binary estimates are correlated, and bias reduction can produce different finite-sample changes at each threshold. Their separate standard errors are descriptive; neither interval overlap nor this plot is a formal test of proportional odds. Thresholds with no outcome variation or unsuccessful fits retain their status and have unavailable estimates. A threshold is also unavailable when binary count conversion would round a positive minority frequency away (unrepresentable_binary_counts). Calibration alone is also not a test of proportional odds.

With po.test="score", a common-slope null on the fitted link is refitted to the included assessment outcomes, predictors, frequencies, and offsets. Maximum likelihood is the default, including when object used mean bias reduction. Set po.estimation="meanBR" to request a mean-bias-reduced null explicitly. The resulting ordinary efficient score retains its nuisance projection and has only first-order justification; it can be conservative in small samples. The null refit affects neither the original model nor the prediction scores. With supplied outcomes, the test concerns the common-slope restriction in those outcomes, rather than calibration of the original model's coefficients. The score uses expected Fisher information, averaged analytically over possible outcomes under this null, and projects out the threshold, common-slope, and any untested varying-effect nuisance scores. No resampling is used for its p-value.

If s design columns are tested across K categories, the asymptotic chi-square reference has s*(K-2) degrees of freedom. For a common-slope fit, other predictors retain common slopes in both the null and the local alternative. For a partial fit, other varying effects remain unrestricted. The latter uses the term-equality construction in fastglm_po_test. The descriptive binary-slope comparison is available only for common-slope fits. Independent observations and a correctly specified, identifiable interior null are required. The approximation can be inaccurate in small or sparse samples; a large p-value does not establish proportional odds. The test is not a robust, cluster-adjusted, or higher-order bias-corrected test.

The score test uses dense data or a repeatable source, at least three declared categories, at least one tested predictor, and positive frequency in every category. Empty categories, deficient rank, failed null fits, and unavailable information produce missing statistics and p-values with an explicit reason. Categories and degrees of freedom are not silently reduced, and the calculation does not switch to a mean bias-reduced fit. The descriptive slope comparison and formal test can be requested independently. Formula term selection requires retained design assignments from a current ordinal fit, including when model=FALSE.

Scoring and calibration use blocks of rows and require work proportional to the number of rows times the number of categories. The optional binary comparison fits one model per threshold without expanding frequency rows. The score test also processes frequencies without expanding observations; its information matrix grows with the predictor and category dimensions.

For an ordinal source fit with no supplied evaluation data, scoring and calibration reread the recorded source once. The same frequency, offset, category, and prediction calculations as dense assessment are used. The source must return unchanged rows and values. Weights and offsets cannot be replaced in this mode. Summary storage depends on the number of categories and bins; row calculations use the current source block. keep.rows=TRUE collects the row table, while sink can save it block by block. Source assessment describes apparent performance on the fitting outcomes. The optional common-slope score test requires additional source passes for refitting and information accumulation. Descriptive binary-slope comparisons still require explicit evaluation data.

The plot method uses ggplot2 when it is installed and base graphics otherwise. With ggplot2, main, xlab, ylab, xlim, ylim, pch, col, and cex are honored, and main titles the whole figure.

Value

An object of class fastglm_ordinal_assessment, containing scores, the overall frequency-weighted losses; rows, individual row averages and exclusions (NULL when keep.rows=FALSE); calibration, predicted and observed cumulative probabilities by threshold and bin; and proportional.odds, the optional slope comparison table. proportional.odds.test is NULL unless requested; otherwise it records statistic, df, p.value, available, tested terms and columns, null convergence, frequency, warnings, and an unavailability reason and explanation. It also records the evaluation source, category order, and exclusion counts. Plot methods use these saved tables.

References

Epstein, E. S. (1969). A scoring system for probability forecasts of ranked categories. Journal of Applied Meteorology, 8, 985–987.

Gneiting, T. and Raftery, A. E. (2007). Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 102, 359–378. doi:10.1198/016214506000001437.

Brant, R. (1990). Assessing proportionality in the proportional odds model for ordinal logistic regression. Biometrics, 46, 1171–1178. doi:10.2307/2532457.

SAS Institute Inc. (2013). SAS/STAT 12.3 User's Guide. The LOGISTIC Procedure, Score Statistics and Tests, pp. 4246–4248. Expected and observed information define distinct finite-sample score calculations; this implementation uses expected information.

Fit several ordinal responses using shared predictors

Description

Fits responses sequentially using a prepared ordinal design. The design setup is reused whenever effective row frequencies match. Each response uses its own empirical starting thresholds and category mapping; estimates from an earlier response are not used as starting values.

Usage

fastglm_ordinal_batch(
  design,
  y,
  weights = NULL,
  offset = NULL,
  levels = NULL,
  firth = NULL,
  tol = 1e-08,
  maxit = 200L,
  block.size = 2048L,
  control = list(),
  link = NULL,
  nonparallel = NULL
)

Arguments

Value

A list of fastglm_ordinal fits, with the response list's names.

Prepare predictors for repeated ordinal fits

Description

Retains a dense location design and its prediction coding for use with fastglm_ordinal() or fastglm_ordinal_batch(). Repeated fits reuse weighted centering, scaling, rank checks, and the transformed design when the predictors and effective row frequencies agree. Model scores, information, category mappings, and starting thresholds are recalculated for each fit.

Usage

fastglm_ordinal_prepare(
  x,
  data = NULL,
  weights = NULL,
  offset = NULL,
  levels = NULL,
  firth = FALSE,
  link = "logit",
  nonparallel = NULL
)

Arguments

Details

The cache contains ordinary R objects and can be saved with saveRDS(). Changing the design, effective row totals, selected nonparallel terms, or the standardize optimizer control rebuilds the relevant setup. For grouped counts, effective totals are counts times their fitting weights. Changing category allocations while retaining those totals permits reuse.

Fitting with a prepared design preserves the formula's factor coding and transformations for prediction. A replacement fitting offset is a total location offset. After changing prepared formula offsets, predictions on new data require an explicit replacement offset.

Value

A fastglm_ordinal_design object. Its cache environment is shared by repeated fits; each fitted model retains its own response and the common original design, without copying the cached transformed design into the fit.

Fit an ordinal model from repeatable blocks of observations

Description

Fits the same cumulative-link likelihood and mean bias-reduction equations as fastglm_ordinal() while retaining one source block and matrices whose dimensions depend on the number of coefficients. The source is read again for successive fitting evaluations. Mean bias reduction first combines the information from every block, then rereads the blocks to evaluate its adjustment using that global information.

Usage

fastglm_ordinal_streaming(
  source,
  n_chunks = NULL,
  levels = NULL,
  firth = FALSE,
  link = "logit",
  tol = 1e-08,
  maxit = 200L,
  block.size = 2048L,
  control = list(),
  nonparallel = NULL
)

Arguments

Details

The rows and values returned for each block must remain unchanged during fitting and later prediction or residual calculations. Empty blocks are not accepted; blocks with zero total frequency are allowed. Category reduction and rank checks use the complete source, so individual blocks may omit categories or have deficient predictor rank.

Storage depends on the largest source block, the number of coefficients, and block-size metadata. The source itself must also obtain blocks without retaining the complete dataset if bounded memory is required. Formula sources use their recorded predictor schema and fixed transformation basis.

Coefficients, model-based covariance, summaries, diagnostics, and Wald intervals use the usual ordinal methods. Robust and clustered covariance through fastglm_vcov() rereads the source and accumulates the same ordinary-score contributions as retained dense fits. Predictions, residuals, and weights can be returned one block at a time through a callback. Collecting their complete output requires memory proportional to the number of rows. The fitted design and response are not retained. fastglm_ordinal_assess() accumulates prediction scores and calibration over the original source; fastglm_po_test() refits common-slope hypotheses through source passes. Optional boundary rows from fastglm_diagnostics() also use one source pass. Profile intervals, general constrained coefficient-score tests, and refitting bootstrap calculations require retained observations. Average prediction contrasts can receive explicit evaluation data.

Value

A fastglm_ordinal_stream object inheriting from fastglm_ordinal, with source and reader references, source pass counts, and the usual ordinal estimates and convergence diagnostics.

Examples

set.seed(104)
x <- matrix(rnorm(120),ncol=1,dimnames=list(NULL,"exposure"))
y <- ordered(cut(.4*x[,1]+rlogis(120),c(-Inf,-.5,.7,Inf)))
reader <- function(k) {
    i <- ((k-1)*30+1):(k*30)
    list(X=x[i,,drop=FALSE],y=y[i])
}
fit <- fastglm_ordinal_streaming(reader,4,firth=TRUE)
summary(fit)
head(predict(fit))

Test common slopes in a cumulative-link model

Description

Computes a common-slope test directly, without calculating prediction losses or calibration summaries. Print the result for a compact test table, or use as.data.frame() or broom::tidy() for reporting.

Usage

fastglm_po_test(
  object,
  terms = NULL,
  null.estimation = c("ML", "meanBR"),
  newdata = NULL,
  y = NULL,
  weights = NULL,
  offset = NULL,
  control = list(),
  by.term = FALSE,
  method = c("score", "lr")
)

Arguments

Details

Both null-estimation choices for score tests use expected information and an asymptotic chi-square reference. Mean bias reduction gives a first-order equivalent score test under a regular interior null with fixed predictor and category dimensions. The nuisance projection is essential because the ordinary nuisance score is nonzero at a bias-reduced root. This is not an adjusted-score test with established higher-order accuracy. It can be conservative in small samples. The method never switches null estimators automatically. The fitted link is retained. For logit models the hypothesis is proportional odds; for probit and complementary log-log models it is the common-slope restriction on that link's scale.

For a partial model, a selected predictor with slopes gamma[,j] is tested through gamma[,j] - gamma[,1] = 0 for every threshold after the first. Other varying effects remain unrestricted in both the null and alternative. If s selected design columns have varying slopes across K categories, this imposes s*(K-2) independent restrictions. Binary outcomes and selections containing only already-common terms return an unavailable test with zero degrees of freedom and an explanation.

For partial score tests, null.estimation="meanBR" fits bias reduction in the reduced null model, embeds that estimate in the full partial model, and projects its ordinary score using full-model expected information. This preserves the common-slope test convention. It differs from fastglm_score_test() with null.estimation="meanBR", which uses the full-model adjustment in constrained fitting and in the score statistic. Likelihood ratio tests use ordinary ML likelihoods and first-order chi-square references; an unsuccessful full or constrained fit is reported as unavailable. No finite-sample or higher-order calibration is claimed.

Empty categories and deficient predictor rank remain unsupported. A sufficient separation check identifies two-valued predictors with a level observed only in an endpoint category. It does not detect all forms of separation; other unsuccessful null fits retain their convergence details. See fastglm_ordinal_assess() for the hypothesis, assumptions, and frequency conventions, and fastglm-test-methods for result extraction. With by.term=TRUE, term-specific p-values are unadjusted for multiple testing. The joint test remains the single test of all selected terms.

An ordinal source fit is tested on its original source when newdata and y are omitted. The source is replayed for refitting and information accumulation without collecting its design or outcomes. Recorded frequency weights and offsets are used; overrides require explicit evaluation data. Score tests retain the same null-estimation and term-selection conventions as dense fits. Partial-model likelihood ratio tests fit the equivalent reduced common-slope model and the full model by streamed ML. Memory is bounded by one source block and coefficient-sized matrices, in addition to the fitted source's existing block metadata. The reader must reproduce the same observations on every pass, as required for source fitting itself.

Value

A fastglm_po_test object inheriting from fastglm_test and htest, with numeric statistic, degrees of freedom, p-value, availability, null estimation method, tested terms, and detailed fitting diagnostics. With by.term=TRUE, a fastglm_test_list containing the joint and term-specific tests in ⁠$tests⁠, with print and data-frame methods.

References

Kosmidis, I. (2014). Improved estimation in cumulative link models. Journal of the Royal Statistical Society, Series B, 76, 169–196. doi:10.1111/rssb.12025. The present meanBR option uses the ordinary projected score at a bias-reduced null, rather than the larger-model adjusted-score statistic discussed in Section 12 of that article.

Examples

fit <- fastglm_ordinal(Sat ~ Infl + Type + Cont, MASS::housing,
                       weights = Freq, firth = TRUE)
fastglm_po_test(fit)
as.data.frame(fastglm_po_test(fit, terms = "Type"))

Prepare a reusable design for repeated GLM fits

Description

Records the predictors and weights shared by several outcome models. Formula inputs also retain the factor coding and transformations needed for prediction. Gaussian identity fits can reuse their matrix decomposition, the expensive calculation that depends only on predictors and weights. Dense Firth fits with Gaussian identity, Gamma log, or Poisson square-root links also reuse their common decomposition and bias-adjustment setup with the automatic or LLT solver. Saved calculations are rebuilt when predictors or weights change. Mean-bias-reduced Gaussian identity, Gamma log, and Poisson square-root fits also reuse the observation leverages needed for the adjusted score. A leverage measures how strongly a row contributes to the fitted regression. The shared setup is used by 'method = "auto"' and 'method = "llt"' (or '2').

Usage

fastglm_prepare(x, data = NULL, weights = NULL, offset = NULL, firth = FALSE)

Arguments

Value

A 'fastglm_design' object for [fastglm_batch()].

Profile likelihood inference

Description

Tests one or more coefficient values by comparing an unrestricted maximum with the maximum under those constraints. The remaining coefficients are re-estimated. Binomial-logit Firth fits use their penalized likelihood; ordinal maximum-likelihood fits use their ordinary likelihood. For intervals, use confint(object, method = "profile"); the default confint() method remains Wald inference.

Usage

fastglm_profile_test(object, parm, null = 0, control = list(), contrast = NULL)

Arguments

Details

For binomial-logit Firth models the objective is the binomial log likelihood plus one half of the log determinant of the full Fisher information matrix. Tested columns remain in that determinant when their coefficients are constrained. Grouped binomial responses, frequency weights, and offsets retain their fitted meanings. Effective successes and failures must be nonnegative integer counts, allowing arithmetic round-off. Fractional power or survey weights are unsupported. Zero-weight observations make no contribution.

For ordinal maximum-likelihood fits, the objective is the joint cumulative-link log likelihood, including partial proportional-odds and other nonparallel cumulative-link models. Thresholds and untested slopes are re-estimated subject to the selected coefficient constraints. Logit, probit, and complementary log-log links, grouped category counts, frequency weights, and offsets retain their fitted meanings. Ordinal mean bias reduction does not generally maximize a penalized likelihood and is unsupported by this profile method; fastglm_score_test() provides an analytic score procedure with an explicit null-estimation convention.

The likelihood-ratio statistic uses an asymptotic chi-square reference with one degree of freedom per constraint. Profile intervals invert its one-parameter version. These are approximations under an independent binomial or ordinal sampling model, including frequency-weighted observations; they are not exact finite-sample or cluster-robust procedures.

Sparse, mapped, and streaming designs are currently unsupported rather than converted to dense matrices. Native evaluation uses additional memory of order the squared number of coefficients plus a fixed block of rows. Ordinal constrained inference also creates a standardized dense design, requiring additional storage proportional to observations times predictors; it does not construct an observation-by-category-by-parameter array.

Every constrained optimization checks the ordinary or adjusted score of the free coefficients. The unrestricted fit is checked and refined when needed. Failed tests return missing statistics and p-values with a warning; failed confidence limits are missing, with endpoint diagnostics attached. No Wald replacement is made. Optimization starts from the fitted solution; convergence diagnostics do not establish that a global maximum has been found in a multimodal problem.

Value

fastglm_profile_test() returns a fastglm_test object inheriting from htest, with a compact printed table and as.data.frame() and optional broom::tidy() extraction. It retains the statistic, degrees of freedom, p-value, coefficient estimates, hypothesized values, the full and constrained log likelihoods (penalized for Firth fits, excluding constants independent of coefficients), and diagnostics containing convergence, scaled scores, optimizer codes, evaluation counts, and messages. confint(..., method = "profile") returns a two-column fastglm_confint matrix with diagnostics, fit.diagnostics, and method attributes. Each diagnostic row identifies the coefficient and endpoint, convergence, likelihood-ratio statistic and target, root error, scaled nuisance score, optimizer code, and failure message.

References

Heinze, G. and Schemper, M. (2002). A solution to the problem of separation in logistic regression. Statistics in Medicine, 21, 2409–2419. doi:10.1002/sim.1047.

Examples

d <- data.frame(y = rep(0:1, each = 10), x = rep(0:1, each = 10))
fit <- fastglm(y ~ x, d, family = binomial(), firth = TRUE)
confint(fit, parm = "x", method = "profile")
fastglm_profile_test(fit, parm = "x")

Score tests for ordinal model terms and coefficient contrasts

Description

Evaluates ordinary or mean-bias-adjusted likelihood scores under linear constraints, after re-estimating the nuisance parameters. The information matrix and, when requested, the bias adjustment refer to the full model.

Usage

fastglm_score_test(
  object,
  terms = NULL,
  contrast = NULL,
  null = 0,
  null.estimation = c("ML", "meanBR"),
  control = list()
)

Arguments

Details

Let B span the nuisance directions allowed by the constraints, I be full-model expected information, and U the selected score vector. The test removes the nuisance component using U_e=U-IB(B^TIB)^{-1}B^TU, and reports U_e^T I^{-1}U_e. With null.estimation="meanBR", constrained fitting solves B^TU^*=0 using the adjustment of the full model; it does not recompute a bias adjustment for a smaller model with the tested columns removed.

For independent ordinal observations or integer frequency counts, fixed parameter dimension, an interior true parameter, and nonsingular limiting information, with the usual local smoothness and consistent constrained estimation, the statistic has a first-order chi-square reference with one degree of freedom per independent constraint. The mean-bias adjustment is order one, while the likelihood score is order square root of sample size. This establishes the same first-order reference, not a higher-order calibration result or exact small-sample inference. No resampling is used.

Infeasible starts, failed nuisance fits, singular information, or exhausted iterations return an unavailable test with an explanation. A finite bias-reduced estimate does not guarantee a successful constrained fit. For nonparallel models, constraints and all nuisance updates retain ordered cumulative predictors at each positive-frequency training row. Equality contrasts between a predictor's threshold-specific slopes test its common slope restriction while allowing other predictors to remain nonparallel. Tests of a term's zero effects are different from tests of equal slopes. Streaming fits and collapsed or empty categories are currently unsupported.

Value

A fastglm_test object with the package's common printed table, numeric as.data.frame() extraction, and constrained-fit diagnostics.

References

Kosmidis, I. (2014). Improved estimation in cumulative link models. Journal of the Royal Statistical Society: Series B, 76, 169–196. doi:10.1111/rssb.12025. The nuisance projection and full-model constrained adjustment used here specify the implementation beyond the paper's general discussion of adjusted-score tests.

Examples

fit <- fastglm_ordinal(Sat ~ Infl + Type + Cont, MASS::housing, weights=Freq)
fastglm_score_test(fit, terms="Type")
fastglm_score_test(fit, terms="Type", null.estimation="meanBR")

Read a large dataset in repeatable blocks of rows

Description

A source describes how to obtain a block of observations when a fitting function needs it. This permits repeated model-fitting passes without loading the complete dataset into memory. The supplied function, called a callback, receives the block number and returns its observations.

Usage

fastglm_source(
  callback,
  n_chunks,
  reset = NULL,
  formula = NULL,
  prototype = NULL,
  weights = NULL,
  offset = NULL
)

Arguments

Details

When a formula is supplied, 'prototype' is an example data frame used to establish the design columns, factor levels, and transformations. This recorded specification, called the schema, is reused for every block so that a column has the same meaning throughout the analysis. For splines and other data-dependent transformations, the prototype must establish the intended basis, including knots and boundary knots.

Value

A 'fastglm_source' object for [fastglm_streaming()], [fastglm_ordinal_streaming()], or the count and two-part fitting functions.

Create a source from Parquet partitions

Description

Each file is one bounded chunk. The formula schema is supplied explicitly so that factor contrasts and transformations remain fixed across files.

Usage

fastglm_source_arrow(files, formula, prototype, weights = NULL, offset = NULL)

Arguments

Value

A 'fastglm_source'.

Create a source from indexed database queries

Description

Each query must return a stable, disjoint chunk in a deterministic order. The database snapshot must remain unchanged throughout fitting and inference.

Usage

fastglm_source_dbi(
  connection,
  queries,
  formula,
  prototype,
  weights = NULL,
  offset = NULL
)

Arguments

Value

A resettable 'fastglm_source'.

Fit a GLM by reading blocks of observations as needed

Description

fastglm_streaming() fits a GLM while keeping only a block of design rows and coefficient-sized summary matrices in memory. Reading successive blocks this way is called streaming. It permits fitting data that cannot be loaded into memory together, at the cost of rereading the observations as the coefficient estimates are updated.

Usage

fastglm_streaming(
  chunk_callback,
  n_chunks = NULL,
  family = gaussian(),
  start = NULL,
  method = 2L,
  tol = 1e-08,
  maxit = 100L,
  firth = FALSE
)

Arguments

Details

Iteratively reweighted least squares (IRLS) updates the coefficients by solving a sequence of weighted regressions. The updates and safeguards that shorten unsuccessful steps run in C++. For supported standard families, the distribution calculations also run in C++; the R function is called to supply each block of observations.

Cholesky methods accumulate the weighted crossproduct ⁠X' W X⁠. The QR method instead updates a triangular factor as blocks arrive, avoiding crossproducts in initialization, coefficient updates, and covariance calculation. It can retain precision for nearly redundant predictors at the cost of more arithmetic. Rank is assessed after combining all blocks; individual blocks may contain fewer rows than columns or lack factor levels. Globally redundant predictors cause an error.

Standard errors and cov.unscaled use the final weighted design. For QR fits, prediction uncertainty, leverage, robust covariance, and default average-prediction inference use the triangular factor directly. Supplying a covariance matrix explicitly to fastglm_average() uses that matrix. Storage includes the current block and matrices with dimensions determined by the number of coefficients, including a quadratic storage requirement in that number. Neither method retains all design rows.

QR bias-reduced convergence checks evaluate the actual adjusted equations using current working weights and leverage values. Verification requires additional source passes. For poorly conditioned designs, the calculation forms and refactors a better-scaled basis before checking these equations. A small score alone cannot certify convergence: the remaining predictor change and dispersion checks must also pass. Failed step searches retain a failure flag.

Value

An object with classes "fastglm_stream" and "fastglm" containing coefficient estimates and model summaries as in fastglm(), including coefficients, cov.unscaled, deviance, iter, and converged. The object retains the function used to read the observations. Prediction and fastglm_vcov() call that function again, so the data must remain accessible and unchanged. QR fits also retain a qr list containing a triangular matrix R, column permutation pivot, and column scales scale for the final weighted design. QR bias-reduced fits additionally save remaining.predictor.step, adjusted.score.rms, remaining.dispersion.step, and adjusted.dispersion.score.rms. The dispersion diagnostics are missing for fixed-dispersion families. refined.convergence.verification records whether verification required additional precision for the design.

Files and databases

fastglm_source_arrow() reads one Parquet file per block, and fastglm_source_dbi() reads the rows returned by one database query per block. fastglm_source() records the factor levels and transformations needed to give the design columns the same meaning in each block.

References

Zietkiewicz, P. and Kosmidis, I. (2024). Bounded-memory adjusted scores estimation in generalized linear models with large data sets. https://arxiv.org/abs/2307.07342.

Examples

# Simulate a "data source" that yields the design matrix in 4 row-blocks.
set.seed(1)
n <- 1000; p <- 5
X <- cbind(1, matrix(rnorm(n * (p - 1)), n, p - 1))
y <- rbinom(n, 1, plogis(X %*% c(0.2, 0.5, -0.3, 0.4, -0.2)))
chunk_size <- 250
chunks <- function(k) {
  idx <- ((k - 1) * chunk_size + 1):(k * chunk_size)
  list(X = X[idx, , drop = FALSE], y = y[idx])
}

fit_stream <- fastglm_streaming(chunks, n_chunks = 4, family = binomial())
fit_full   <- fastglm(X, y, family = binomial(), method = 2)
max(abs(coef(fit_stream) - coef(fit_full)))

Zero-truncated count regression

Description

Fits Poisson or negative-binomial regression conditional on a strictly positive response. The linear predictor describes the mean of the underlying untruncated distribution; response predictions describe the conditional mean among positive observations.

Usage

fastglm_truncated(
  x,
  y = NULL,
  data = NULL,
  dist = c("poisson", "negbin"),
  weights = NULL,
  offset = NULL,
  theta = NULL,
  start = NULL,
  tol = 1e-08,
  maxit = 100L,
  outer.tol = 1e-08,
  outer.maxit = 100L
)

Arguments

Value

A 'fastglm_truncated' model with joint covariance when theta is estimated.

Two-part regression for nonnegative continuous outcomes

Description

Fits a binomial model for a positive outcome and a Gamma model with a log link among positive outcomes. The marginal mean is the product of the probability of a positive outcome and the conditional positive mean.

Usage

fastglm_twopart(
  formula,
  data,
  subset,
  na.action,
  weights,
  offset,
  link = "logit",
  method = "auto",
  tol = 1e-08,
  maxit = 100L,
  firth = FALSE,
  response = NULL,
  z = NULL,
  block.size = 8192L
)

Arguments

Value

A 'fastglm_twopart' model supporting prediction, delta-method standard errors, simulation, and component summaries.

Robust covariance with bounded design memory

Description

Computes the sandwich meat in row blocks, optionally accumulating scores by cluster. Cluster storage is proportional to the number of distinct clusters times the number of coefficients. Streaming fits reread the source.

Usage

fastglm_vcov(
  object,
  type = c("HC3", "HC2", "HC1", "HC0"),
  cluster = NULL,
  cadjust = TRUE,
  block.size = 8192L,
  full = FALSE
)

Arguments

Details

For count and two-part models, the sandwich covariance uses joint estimating scores and the observed derivative of their sum. Cluster scores include both model components, so within-cluster dependence can produce nonzero covariance between them. The Gamma mean equation does not depend on dispersion; removing its common dispersion multiplier preserves the fitted mean parameters and their sandwich covariance. These approximations require independent observations, or independent clusters when supplied, and identifiable interior parameters. For GLM, count, and two-part models, sandwich covariance for bias-adjusted scores is not implemented.

Ordinal robust covariance supports dense, prepared, and repeatable source fits, with common and threshold-specific slopes jointly. Source fits reread the observations once without retaining the full design or response. The source determines the size of the loaded block; block.size limits the temporary calculations within that block. Cluster calculations retain one score vector and identifier per positive-frequency cluster. For ordinal fits, the covariance uses the inverse observed negative Jacobian of the ordinary likelihood score, rather than inverse expected Fisher information. Independent-observation meat sums individual score outer products linearly in the effective category frequencies. Cluster meat instead sums individual scores within each cluster before taking outer products. Thus a frequency row represents repeated independent individuals unless a cluster identifier groups them together. With total frequency N, parameter count d, and positive-frequency cluster count G, HC1 multiplies independent covariance by N/(N-d) and cluster covariance by (N-1)/(N-d); cadjust = TRUE additionally applies G/(G-1) to cluster covariance. Every declared category must have positive frequency, and the fit must have converged.

For mean-bias-reduced ordinal fits, this is explicitly a first-order ordinary-score sandwich evaluated at the bias-reduced estimate, not the covariance of the exact adjusted equations. Its asymptotic interpretation requires an interior regular limit, fixed parameter dimension, nonsingular limiting score derivative, an appropriate independent-observation or independent-cluster central limit theorem, and a score adjustment negligible relative to the square root of the sample size. Finite estimates under separation do not establish those conditions. This convention does not promise finite-sample calibration or bias reduction under misspecification. See fastglm_bootstrap() for refit-based uncertainty from retained-data fits; ordinal source bootstrapping remains unavailable.

Value

A coefficient covariance matrix.

Joint Wald tests for model terms and coefficient contrasts

Description

Tests specified linear combinations of regression coefficients in a fitted model. A factor or interaction can be tested as a complete formula term, and a named contrast can compare coefficients without refitting.

Usage

fastglm_wald_test(
  object,
  terms = NULL,
  contrast = NULL,
  null = 0,
  vcov. = NULL,
  component = "all"
)

Arguments

Details

For contrast matrix L, coefficient estimate b, and covariance V, the statistic is ⁠(L b - null)' solve(L V L', L b - null)⁠. Its reference is chi-square with one degree of freedom per independent contrast. This is a first-order approximation requiring asymptotically normal, estimable coefficients and a consistent covariance for the sampling design. It is not a finite-sample F test or a likelihood-ratio test.

Mean-bias-reduced estimates use the same first-order convention. Sparse outcomes, extreme effects, and few independent clusters can make the approximation inaccurate. Tests of regression coefficients do not assess the common-slope restriction; use fastglm_po_test() for that hypothesis.

Nonconvergence, collapsed ordinal categories, nonestimable contrasts, and missing or numerically singular contrast covariance return an unavailable result with an explanation. Degrees of freedom are not reduced silently. For aliased GLM coefficients, estimability is checked against the retained dependence between design columns. No generalized inverse is used. Tests concern regression coefficients, excluding dispersion or size parameters. An estimated negative-binomial size at its numerical boundary returns an unavailable result. Truncated, hurdle, and zero-inflated fits use regression covariance accounting for estimated size in regular fits. For fastglm_nb(), the default retains its conventional fixed-size GLM covariance; to use joint observed information, supply the regression block of vcov(object, full=TRUE). Supplied covariance must contain all regression coefficients, including both components when applicable, even for a selected-component test. Cross-component covariance is retained; it is not replaced with independent component calculations. The covariance restrictions of fastglm_vcov() continue to apply, including its restriction on adjusted two-part fits.

Value

A fastglm_wald_test inheriting from fastglm_test and htest, with statistic, degrees of freedom, p-value, availability, contrast estimates, standard errors, null values, and the contrast matrix and its covariance. Print the object directly or use as.data.frame() for numeric extraction.

Examples

fit <- fastglm_ordinal(Sat ~ Infl + Type + Cont, MASS::housing,
                      weights=Freq, firth=TRUE)
fastglm_wald_test(fit, terms="Type")
fastglm_wald_test(fit, contrast=c(InflHigh=1, InflMedium=-1))

Zero-inflated Poisson / Negative-Binomial regression

Description

Fits a zero-inflated count model: a binary inflation component combined with a Poisson or NB count component, where the count component is the *original* (non-truncated) distribution. Excess zeros above what the count model implies are absorbed by the inflation component.

Usage

fastglm_zi(
  formula,
  data,
  subset,
  na.action,
  weights,
  offset,
  dist = c("poisson", "negbin"),
  link = c("logit", "probit", "cloglog", "log"),
  init.theta = NULL,
  em.tol = 1e-08,
  em.maxit = 100L,
  tol = 1e-09,
  maxit = 100L,
  theta.tol = 1e-08,
  theta.maxit = 100L,
  model = TRUE,
  y = TRUE,
  x = FALSE
)

Arguments

Details

Estimation uses an EM algorithm with all numerical loops in C++. The E-step computes the posterior 'tau_i = P(Z_i = 1 | y_i)' analytically; M-steps fit the inflation logit/probit/cloglog/log and the Poisson/NB count regression via the existing native fastglm IRLS (with weights 'w_i (1 - tau_i)' on the count side). For NB, an inner Brent MLE re-estimates 'theta' after the count-side beta step. Final vcov comes from the numerical Jacobian of the analytical observed score at the EM fixed point (block-structured for '(gamma, beta, theta)').

This is the fastglm analogue of [pscl::zeroinfl()] – coefficients agree to 1e-5 on standard datasets (slightly looser than hurdle to allow for the EM iteration), with all numerical loops in C++.

Value

A list of class 'c("fastglm_zi", "fastglm")' with elements 'coefficients' (a list with 'count' and 'zero'), 'se' (likewise), 'vcov' (full, including theta if NB), 'loglik', 'theta', 'SE.theta', 'tau' (posterior P(Z=1|y)), 'em_iter', 'converged'.

Examples

set.seed(1)
n <- 400
x  <- rnorm(n)
eta_count <- 0.8 + 0.4*x
lam <- exp(eta_count)
eta_zero <- -0.6 + 0.5*x
p_inflate <- plogis(eta_zero)
z <- rbinom(n, 1, p_inflate)
y <- ifelse(z == 1, 0, rpois(n, lam))
df <- data.frame(y = y, x = x)
fit <- fastglm_zi(y ~ x, data = df, dist = "poisson")
coef(fit, model = "count")
coef(fit, model = "zero")

Negative binomial family with known dispersion

Description

A built-in negative-binomial family for use with [fastglm()] and [fastglmPure()] when 'theta' (the NB2 dispersion) is known. Equivalent to 'MASS::negative.binomial(theta, link)' but without taking a hard dependency on the **MASS** package. The variance is 'mu + mu^2 / theta'; as 'theta -> Inf' it reduces to Poisson.

Usage

negbin(theta, link = "log")

Arguments

Details

For joint estimation of 'theta' together with 'beta', use the dedicated [fastglm_nb()] function. The 'negbin()' family here is intended for use cases where 'theta' has been pre-specified or estimated separately.

Value

A 'family' object with class '"family"' and 'family$family == "Negative Binomial(theta)"'. The object also carries the slot 'family$theta', which retains the supplied size parameter.

Examples

set.seed(1)
n <- 500
x  <- cbind(1, matrix(rnorm(n * 2), n, 2))
mu <- exp(x %*% c(0.3, 0.4, -0.2))
y  <- rpois(n, lambda = mu * rgamma(n, shape = 2, rate = 2))   # NB(theta=2)

fit <- fastglm(x, y, family = negbin(theta = 2, link = "log"))
coef(fit)

Methods for ordinal models with threshold-specific slopes

Description

Partial proportional-odds models allow selected predictor effects to differ across cumulative thresholds. The same parameter layout is available for probit and complementary log-log links, with interpretation on their link scales. Common slopes and threshold-specific slopes are displayed separately.

Usage

## S3 method for class 'fastglm_ordinal_partial'
coef(
  object,
  component = c("all", "location", "threshold", "common", "nonparallel"),
  ...
)

## S3 method for class 'fastglm_ordinal_partial'
vcov(
  object,
  component = c("all", "location", "threshold", "common", "nonparallel"),
  type = c("model", "HC0", "HC1"),
  cluster = NULL,
  cadjust = TRUE,
  block.size = 8192L,
  ...
)

## S3 method for class 'fastglm_ordinal_partial'
summary(object, vcov. = NULL, ...)

## S3 method for class 'summary.fastglm_ordinal_partial'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'fastglm_ordinal_partial'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'fastglm_ordinal_partial'
confint(
  object,
  parm,
  level = 0.95,
  method = c("wald", "profile"),
  control = list(),
  ...
)

## S3 method for class 'fastglm_ordinal_partial'
predict(
  object,
  newdata = NULL,
  type = c("probability", "cumulative", "class", "link", "response"),
  se.fit = FALSE,
  offset = NULL,
  ...
)

tidy.fastglm_ordinal_partial(
  x,
  conf.int = FALSE,
  conf.level = 0.95,
  exponentiate = FALSE,
  component = c("all", "location", "threshold", "common", "nonparallel"),
  vcov. = NULL,
  ...
)

Arguments

Details

The coefficient vector contains thresholds, common slopes, and threshold-specific slopes, in that order. Within the last block, all selected predictors for the first threshold precede those for the second threshold. object$gamma also stores these slopes as a predictor-by-threshold matrix. Selecting component="location" includes both kinds of slope.

Thresholds describe the cumulative predictors when the covariates and offset equal zero. That covariate combination can lie outside the valid prediction domain, so the reported thresholds need not be ordered in isolation. The cumulative predictors are constrained to be ordered at positive-frequency training rows and checked again at prediction rows. Invalid predictions report the affected rows. Probability clipping or threshold reordering is not used to repair an invalid prediction.

Model-based covariance is inverse expected Fisher information, including covariance between all thresholds and slopes. Wald intervals and delta-method prediction standard errors are first-order approximations under an interior correctly specified model with independent observations or integer frequency counts and nonsingular limiting information. Numerical convergence alone does not establish those assumptions. With firth=TRUE, the full-model cumulative-link adjustment targets the leading mean bias of the reported thresholds, common slopes, and threshold-specific slopes. These adjusted equations do not generally maximize a scalar penalized likelihood. Empty-category reductions require additional parameter constraints and are not supported for these fits.

Value

Coefficients, covariance, predictions, or a summary as appropriate. Printing invisibly returns the supplied object.

Obtains predictions and optionally estimates standard errors of those predictions from a fitted generalized linear model object.

Description

Obtains predictions and optionally estimates standard errors of those predictions from a fitted generalized linear model object.

Usage

## S3 method for class 'fastglm'
predict(
  object,
  newdata = NULL,
  type = c("link", "response", "terms"),
  se.fit = FALSE,
  dispersion = NULL,
  offset = NULL,
  terms = NULL,
  ...
)

Arguments

Value

Predictions on the requested scale. With 'se.fit=TRUE', a 'fastglm_prediction' list containing 'fit', 'se.fit', and 'residual.scale'; see fastglm-prediction-methods for printing and table extraction.

Predictions and inference for two-component models

Description

The response prediction averages over both components. For hurdle models, 'count' is the conditional mean among positive counts and 'zero' is the probability of observing zero. For zero-inflated models, 'count' is the mean of the untruncated count distribution. Prediction standard errors use the delta method and include uncertainty in theta when it is estimated.

Usage

## S3 method for class 'fastglm_hurdle'
predict(
  object,
  newdata = NULL,
  type = c("response", "count", "zero", "prob", "link"),
  se.fit = FALSE,
  at = 0:10,
  sink = NULL,
  ...
)

## S3 method for class 'fastglm_zi'
predict(
  object,
  newdata = NULL,
  type = c("response", "count", "zero", "prob", "link"),
  se.fit = FALSE,
  at = 0:10,
  sink = NULL,
  ...
)

Arguments

Value

Predicted values, or a 'fastglm_prediction' list with 'fit' and 'se.fit'; see fastglm-prediction-methods for printing and table extraction.

Predict ordinal probabilities or outcome categories

Description

Predict ordinal probabilities or outcome categories

Usage

## S3 method for class 'fastglm_ordinal'
predict(
  object,
  newdata = NULL,
  type = c("probability", "cumulative", "class", "link", "response"),
  se.fit = FALSE,
  offset = NULL,
  ...
)

Arguments

Details

With collapsed empty interior categories, probability and cumulative uncertainty refer to the reduced category model. The removed category has fitted probability zero and unavailable probability standard error.

Value

A vector or matrix, or a fastglm_prediction list with fit, se.fit, and covariance.scope. See fastglm-prediction-methods for printing and table extraction.

Predictions and residuals from an ordinal source fit

Description

Predictions and residuals from an ordinal source fit

Usage

## S3 method for class 'fastglm_ordinal_stream'
predict(
  object,
  newdata = NULL,
  type = c("probability", "response", "cumulative", "class", "link"),
  se.fit = FALSE,
  offset = NULL,
  sink = NULL,
  ...
)

## S3 method for class 'fastglm_ordinal_stream'
residuals(
  object,
  type = c("probability", "cumulative", "pearson"),
  sink = NULL,
  ...
)

Arguments

Value

Collected predictions or residuals, or invisibly the number of blocks supplied to 'sink'.

Stream predictions or diagnostics to a callback

Description

Prediction for a large dataset can create a vector too large to retain in memory. A sink is a user-supplied function that receives one block of results at a time, for example to write predictions to a file or accumulate a summary. Without a sink, the result blocks are collected into a single R object, whose memory requirement grows with the number of observations.

Usage

## S3 method for class 'fastglm_stream'
predict(
  object,
  newdata = NULL,
  type = c("link", "response"),
  se.fit = FALSE,
  dispersion = NULL,
  sink = NULL,
  ...
)

Arguments

Value

Predictions (a 'fastglm_prediction' list when 'se.fit=TRUE'), or invisibly the number of chunks delivered to a sink. Each prediction block supplied to a sink uses the same structure. See fastglm-prediction-methods.

Probability-based residuals for ordinal models

Description

Measures discrepancies between observed categories and the fitted category distribution without assigning numerical distances to category labels.

Usage

## S3 method for class 'fastglm_ordinal'
residuals(object, type = c("probability", "cumulative", "pearson"), ...)

Arguments

Details

For an individual outcome y, the probability-scale residual is F(y-)+F(y)-1, where F is the fitted cumulative distribution. Positive values indicate outcomes above the fitted distribution. For a count response, this method returns the mean of the individual residuals within each row. Probability-scale residuals lie in [-1,1]; for ordinal data they generally have neither a normal nor a uniform distribution.

Cumulative residuals have one column per threshold. For row frequency m, the Pearson version divides cumulative observed count minus m F(j) by \sqrt{m F(j)(1-F(j))}. These columns are dependent and are descriptive diagnostics, not independent goodness-of-fit tests. Zero-frequency rows have missing residuals. Rows omitted under na.exclude are restored with missing values.

Value

A named vector for probability-scale residuals or a matrix for cumulative or Pearson residuals.

References

Li, C. and Shepherd, B. E. (2012). A new residual for ordinal outcomes. Biometrika, 99, 473–480. doi:10.1093/biomet/asr073.

Simulate categorical responses from an ordinal fit

Description

Generates independent categorical observations conditional on the retained predictors and offsets. Frequency rows are simulated as multinomial counts, so a frequency of ten generates ten observations at that predictor row.

Usage

## S3 method for class 'fastglm_ordinal'
simulate(object, nsim = 1L, seed = NULL, ...)

Arguments

Details

The original category labels and their order are preserved. A collapsed category has fitted probability zero and therefore simulated count zero. This is simulation from the fitted distribution; it does not quantify uncertainty about a category absent from that distribution. Row totals are effective frequencies, including counts times fitting weights, and must not exceed .Machine$integer.max. Rows omitted under na.exclude are restored with missing values.

Value

For an individual response with unit frequency in every retained row, a data frame of ordered factors, one column per simulation. Otherwise a named list of count matrices, one per simulation, with rows matching the fitted design and columns in the declared category order. Count matrices already include the fitting frequencies and should be refitted with unit weights. They are not counts to which the original weights are reapplied.

'summary()' method for 'fastglm' fitted objects

Description

'summary()' method for 'fastglm' fitted objects

Usage

## S3 method for class 'fastglm'
summary(
  object,
  dispersion = NULL,
  correlation = FALSE,
  symbolic.cor = FALSE,
  ...
)

Arguments

Value

A 'summary.fastglm' object

See Also

[summary.glm()]

Examples

x <- matrix(rnorm(10000 * 10), ncol = 10)
y <- 1 * (0.25 * x[,1] - 0.25 * x[,3] > rnorm(10000))

fit <- fastglm(x, y, family = binomial())

summary(fit)

Tidy summaries of ordinal fits

Description

Optional methods for the broom generics, with thresholds distinguished from location slopes and the estimation method retained in model summaries.

Usage

tidy.fastglm_ordinal(
  x,
  conf.int = FALSE,
  conf.level = 0.95,
  exponentiate = FALSE,
  component = c("all", "location", "threshold"),
  vcov. = NULL,
  ...
)

glance.fastglm_ordinal(x, ...)

augment.fastglm_ordinal(x, data = NULL, newdata = NULL, se_fit = FALSE, ...)

Arguments

Details

Thresholds have no reported zero-null statistic or p-value. Nonconverged fits retain point estimates with unavailable uncertainty. For bias-reduced fits, the likelihood, AIC, and BIC in glance() use the ordinary likelihood evaluated at those estimates; this is not a maximized likelihood or a bias-reduced information criterion.

augment() appends a matrix column .fitted of category probabilities, an ordered-factor column .class, and, for training rows, probability-scale residuals in .resid. The optional matrix column .se.fit contains model-based probability standard errors. When original data are not supplied, retained model-frame or design values are used; covariates in omitted rows are unavailable and restored as missing.

Value

A data frame with one row per parameter for tidy(), one row per model for glance(), or one row per observation for augment().

Methods for cumulative-link ordinal fits

Description

Methods for cumulative-link ordinal fits

Usage

## S3 method for class 'fastglm_ordinal'
update(object, formula., ..., evaluate = TRUE)

## S3 method for class 'fastglm_ordinal'
coef(object, component = c("all", "location", "threshold"), ...)

## S3 method for class 'fastglm_ordinal'
vcov(
  object,
  component = c("all", "location", "threshold"),
  type = c("model", "HC0", "HC1"),
  cluster = NULL,
  cadjust = TRUE,
  block.size = 8192L,
  ...
)

## S3 method for class 'fastglm_ordinal'
summary(object, vcov. = NULL, ...)

## S3 method for class 'fastglm_ordinal'
confint(
  object,
  parm,
  level = 0.95,
  method = c("wald", "profile"),
  control = list(),
  ...
)

Arguments

Value

Coefficients, covariance, or model summaries, as appropriate.

'vcov()' method for 'fastglm' fitted objects

Description

'vcov()' method for 'fastglm' fitted objects

Usage

## S3 method for class 'fastglm'
vcov(object, dispersion = object$dispersion, type = NULL, cluster = NULL, ...)

Arguments

Value

The estimated variance-covariance matrix of the fitted coefficients. For rank-deficient fits, rows and columns corresponding to aliased coefficients are filled with 'NA'.

Heteroskedasticity-consistent (HC) variance estimators for 'fastglm' objects

Description

Methods for 'sandwich::vcovHC()' on objects of class '"fastglm"' and '"fastglmFit"'. Load 'sandwich' ('library(sandwich)') before calling 'vcovHC(fit)'; otherwise no 'vcovHC' generic is in scope.

Usage

vcovHC.fastglm(object, type = c("HC3", "HC2", "HC1", "HC0"), ...)

vcovHC.fastglmFit(object, type = c("HC3", "HC2", "HC1", "HC0"), ...)

Arguments

Details

Computes the Eicker-Huber-White sandwich estimator 'bread ‘cov.unscaled') and 'meat = X’ diag(omega_i) X'. With 's_i = w_i^2 * r_i' the score contribution from observation 'i', the omegas are:

'HC0'

'omega_i = s_i^2'

'HC1'

'HC0' rescaled by 'n / (n - p)'

'HC2'

'omega_i = s_i^2 / (1 - h_i)'

'HC3'

'omega_i = s_i^2 / (1 - h_i)^2'

where 'r_i' is the working residual '(y - mu) / mu.eta(eta)', 'w_i^2 = prior.weight * mu.eta(eta)^2 / variance(mu)' is the IRLS working weight, and ‘h_i = w_i^2 * x_i’ (X' W X)^(-1) x_i' is the IRLS leverage. Equivalent to 'sandwich::vcovHC.glm'.

Requires the model matrix 'x' stored on the fitted object (set automatically by 'fastglm()', 'fastglmPure()', and 'fastglm_fit()' since version 0.0.6).

Value

A 'p x p' heteroskedasticity-consistent variance-covariance matrix.

Examples

if (requireNamespace("sandwich", quietly = TRUE)) {
  x <- cbind(1, matrix(rnorm(500 * 4), ncol = 4))
  y <- rbinom(500, 1, plogis(x %*% c(0.2, 0.3, -0.4, 0.1, 0.2)))
  fit <- fastglm(x, y, family = binomial())
  sandwich::vcovHC(fit)
  sandwich::vcovHC(fit, type = "HC0")
}