--- title: "Large-data estimation and inference with fastglm" author: "Jared Huling" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Large-data estimation and inference with fastglm} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(collapse=TRUE, comment="#>", message=FALSE, eval=requireNamespace("MASS",quietly=TRUE)) ``` Large regression analyses can require more memory than is available on a single computer. Much of that memory is used by the design matrix, whose rows represent observations and whose columns represent the predictors used in the regression. Memory can be reduced by storing only nonzero entries, accessing a matrix from a file, or reading a limited number of rows at a time. We describe when each approach is useful, then consider how to obtain predictions and standard errors without loading the complete dataset. ```{r} library(fastglm) set.seed(82) n <- 1200 raw <- data.frame(x=rnorm(n),g=factor(rep(letters[1:4],length.out=n))) raw$y <- rpois(n,exp(.2+.3*raw$x)) X <- model.matrix(~x+g,raw) ``` # Sparse and mapped matrices When most design entries are zero, a sparse matrix saves memory by recording the nonzero entries and their locations. The `dgCMatrix` class from the *Matrix* package supplies this representation. It is useful, for example, when a model includes indicators for many categories. The example below converts a small design to sparse form to demonstrate the interface; a design with few zeros may take more memory in this representation. Sparse and file-backed fitting use Cholesky decompositions, which solve the weighted regressions underlying GLM fitting. The choices `"llt"` and `"ldlt"` are two forms of this decomposition. Both require a full-rank design among observations with positive weight: no predictor column may be an exact linear combination of the others. Redundant columns must therefore be removed before these fits. ```{r} Xs <- as(X,"dgCMatrix") sparse_fit <- fastglm(Xs,raw$y,poisson(),method="llt") dense_fit <- fastglm(X,raw$y,poisson(),method="llt") max(abs(coef(sparse_fit)-coef(dense_fit))) ``` A `bigmemory::big.matrix` can store an ordinary matrix in a file and make its entries available as they are needed. This is called memory mapping. The fitting routine reads blocks of rows and avoids making a second complete copy of the design. It still retains several vectors with one entry per observation, so the memory requirement continues to increase with the number of rows. The following example illustrates the same interface using a small matrix held in memory. A view selects rows or columns of a `big.matrix` without copying the complete matrix. Views are supported when the matrix uses contiguous storage. For matrices stored as separate columns, the block-reading approach in the next section provides an alternative. ```{r} Xb <- bigmemory::as.big.matrix(X) mapped_fit <- fastglm(Xb,raw$y,poisson(),method="llt") max(abs(coef(mapped_fit)-coef(dense_fit))) ``` # Preserving a model across chunks When even the observation-level working vectors are too large to retain, the data can be read in smaller blocks, called chunks. Fitting by repeatedly reading these blocks is called streaming. A pass is a visit to all blocks. GLM fitting generally requires several passes because each coefficient update depends on the current fitted means. The function `get_data(k)` below supplies block `k`. A function supplied for another function to call when needed is often called a callback. Here its role is simply to read the requested rows. `fastglm_source()` records this reader together with the number of blocks and the model specification. The same predictor column must have the same meaning in every block. For example, a factor level absent from one block should not disappear from its design matrix, and a spline should use the same knots throughout. The `prototype` argument supplies an example data frame from which the factor levels, coding of categories, and transformations are established. This stored specification is called the formula schema. In the small example we use the complete data frame as the prototype; for a large analysis, a smaller example can be used if it establishes all intended levels and transformations. ```{r} chunk_rows <- 300L get_data <- function(k) raw[((k-1L)*chunk_rows+1L):(k*chunk_rows),] source <- fastglm_source(get_data, n_chunks=4L, formula=y~x+g, prototype=raw) stream_fit <- fastglm_streaming(source, family=poisson()) max(abs(coef(stream_fit)-coef(dense_fit))) c(passes=stream_fit$passes, chunks=stream_fit$chunk.reads) summary(stream_fit) ``` The coefficient summary uses quantities saved during fitting, so `summary()` requires no additional reads. If a reader already constructs the model matrix, it can return `list(X=..., y=..., weights=..., offset=...)` directly. Here `X` contains predictor rows and `y` contains the corresponding responses; weights and offsets are optional. The fitting routine checks each block and applies the selected family's response initialization, including trial totals for grouped binomial outcomes. A file or database reader may need to return to its starting position before a new pass. The optional `reset` function performs that operation before the first block is requested. Each request must reproduce the same rows, values, and order. Reader errors stop the fit, and the data must remain unchanged during fitting and subsequent inference. The algorithm reuses weighted-regression summaries calculated while checking a proposed coefficient update. Additional passes can still be needed to obtain starting values, shorten an unsuccessful step, fit the intercept-only comparison model, or calculate final summaries. The recorded `passes` and `chunk.reads` values help assess this cost when reading the data is slower than the regression calculations. # A more stable solve for related predictors Streaming fits use the `"llt"` Cholesky method by default. At each coefficient update, this method combines weighted products of predictor columns into a small matrix. These products are inexpensive to accumulate across blocks, but they can lose precision when predictors are nearly linear combinations of one another. For example, two measurements may differ only slightly while still carrying distinct information. Rescaling their units can help with differences in magnitude, but it cannot remove their near redundancy. The option `method="qr"` provides an alternative. QR decomposition summarizes the weighted design using orthogonal transformations, which preserve lengths and angles up to numerical rounding. The streaming calculation updates a small triangular matrix as each block arrives. It does not retain the complete design or form the matrix of weighted crossproducts used by the Cholesky methods. This permits more reliable coefficient updates for difficult predictor configurations, at the cost of additional arithmetic within each pass. The blockwise approach follows the bounded-memory GLM calculations described by [Zietkiewicz and Kosmidis (2024)](https://arxiv.org/abs/2307.07342). ```{r} qr_stream <- fastglm_streaming(source,family=poisson(),method="qr") max(abs(coef(qr_stream)-coef(stream_fit))) head(predict(qr_stream,type="response",se.fit=TRUE)$se.fit) ``` For this well-conditioned example, the two methods give essentially the same answer. The QR option is most useful when a Cholesky fit reports that the weighted design is numerically unreliable. It is available for both ordinary and bias-reduced streaming GLMs. Predictions, leverage values, robust covariance, and default average-prediction standard errors use the retained triangular factor to preserve the benefit of the more stable fit. Supplying a covariance matrix explicitly to `fastglm_average()` instead uses that matrix, as requested. The complete weighted design must still have full column rank. QR does not make an exactly redundant predictor identifiable, and this streaming implementation reports a rank failure rather than choosing coefficients for redundant columns. A single block may contain too few rows or lack some factor levels: rank is assessed after combining all blocks. Storage consists of the current block, small records of block sizes, and summaries whose size grows with the square of the number of coefficients. Thus, this approach addresses a large number of observations; a very large number of predictors can still require substantial memory. The default solver and the existing `"auto"` choice are unchanged. Bias-reduced QR fits check the remaining change in the fitted linear predictors together with the joint adjusted score and, when estimated, dispersion. This is useful when nearly redundant predictors have large coefficients that mostly cancel: the size of an individual coefficient update can then give a poor description of the remaining change in the fitted model. The checks use the current working weights and leverage values and require additional source passes. A better-scaled basis is used for verification when the design requires additional numerical precision. Fits that still fail these checks retain their nonconvergence flag. The bias-reduction vignette illustrates the saved diagnostics. # Bias reduction with repeated block reads A large dataset can still contain small or imbalanced subgroups for which bias reduction is useful. The same `firth=TRUE` option is available for sparse matrices, mapped matrices, and repeatable sources. The bias-reduction vignette describes the estimating equations and their supported families. For Gamma-log and Gaussian-identity fits, some of the matrix calculations depend only on the predictors and prior weights. These calculations are retained after their initial construction, reducing the work performed on subsequent passes. Poisson-square-root fits can also reuse the information matrix while continuing to update the response-dependent adjustment. The reader still visits each block in order because fitted means and response contributions must be updated. ```{r} set.seed(149) positive <- raw positive$y <- rgamma(n,shape=2,scale=exp(.2+.3*positive$x)/2) positive_source <- fastglm_source( function(k) positive[((k-1L)*chunk_rows+1L):(k*chunk_rows),], 4L,formula=y~x+g,prototype=positive) bias_stream <- fastglm_streaming(positive_source,Gamma("log"), firth=TRUE,tol=1e-9) bias_dense <- fastglm(X,positive$y,Gamma("log"),firth=TRUE,tol=1e-9) max(abs(coef(bias_stream)-coef(bias_dense))) c(passes=bias_stream$passes,blocks_read=bias_stream$chunk.reads) ``` Gamma dispersion calculations also combine identical prior weights within each block, avoiding repeated distribution calculations for observations that contribute the same term. This combination uses only the current block and does not retain all weights from the source. The fitting state includes matrices whose dimensions depend on the number of coefficients and, for Gamma-log, a coefficient-length adjustment vector. In-memory sparse fits additionally retain one leverage value per observation, consistent with their existing observation-level working vectors. The optimization reduces arithmetic within a pass; it does not remove the need for repeatable data or imply fewer reader passes. # Prediction and robust inference A fitted streaming model retains the instructions for reading its data. Calling `predict()` rereads those data, calculates predictions block by block, and normally joins them into a single result. That result can itself be too large for memory, even though the model was fitted successfully. The `sink` argument avoids collecting the complete result. It takes a function that receives the predictions for one block and the block number. That function can write the predictions to a file or reduce them to a small summary before the next block is read. In the following example, only four block means are retained. The name `sink` refers to this destination for the calculated values. ```{r} chunk_means <- numeric(4) predict(stream_fit,type="response", sink=function(value,k) chunk_means[k] <<- mean(value)) chunk_means fastglm_average(stream_fit) ``` The same argument is available for `residuals()` and `hatvalues()`. Residuals compare observed outcomes with fitted means, while leverage values describe the relative influence of an observation's predictor pattern in the weighted-regression step. Processing their values one block at a time avoids retaining a full diagnostic vector. Standard errors can also be calculated by reading blocks. Robust covariance estimation uses the observed variation in contributions to the estimating equations, allowing the response variance to differ from the variance specified by the fitted family. With correlated observations, contributions are first summed within groups, such as patients or sites. This gives cluster-robust standard errors under the assumption that different groups are independent and that sufficiently many groups are available. The `get_cluster(k)` function below supplies the group identifiers in the same row order as data block `k`. The `HC1` option applies a finite-sample rescaling of this covariance calculation. ```{r} cluster_ids <- rep(seq_len(120),each=10) get_cluster <- function(k) cluster_ids[((k-1L)*chunk_rows+1L):(k*chunk_rows)] V <- fastglm_vcov(stream_fit,type="HC1",cluster=get_cluster) sqrt(diag(V)) fastglm_average(stream_fit,vcov.=V) ``` The square roots of the diagonal entries of `V` are the coefficient standard errors. Passing `V` to `fastglm_average()` uses the same covariance to quantify uncertainty in the average fitted response, conditional on the observed predictor values. These calculations retain one block of design rows, matrices whose size depends on the number of coefficients, and small records of the number of rows in each block. Clustered inference additionally retains a contribution vector for each group, so its memory requirement grows with the number of groups. Returning every prediction or residual in a single vector still requires memory proportional to the number of observations; a supplied `sink` function is useful when that output would be too large. # Ordered outcomes from repeatable blocks Ordinal models with common or selected threshold-specific slopes can also be fitted without retaining the complete design or response. `fastglm_ordinal_streaming()` uses the same cumulative-link likelihood and mean bias-reduction equations as `fastglm_ordinal()`. The source must preserve the category order and predictor definitions across all blocks. A category or factor level may be absent from an individual block; the category totals and predictor rank are assessed across the complete source. Selected threshold-specific effects use `nonparallel=~term` for formula sources, or predictor column names or indices for matrix sources. These partial models require positive total frequency in every declared category. The housing satisfaction data provide a small example with frequency weights. Here the reader selects rows from an in-memory table to make the construction visible. For a dataset larger than memory, the reader would obtain each block from a file or database instead. ```{r} housing <- MASS::housing housing_source <- fastglm_source( function(k) housing[((k-1L)*18L+1L):(k*18L),,drop=FALSE], n_chunks=4L,formula=Sat~Infl+Type+Cont, prototype=housing,weights="Freq") ordinal_stream <- fastglm_ordinal_streaming(housing_source,firth=TRUE) summary(ordinal_stream) fastglm_diagnostics(ordinal_stream) ``` To allow contact with other residents to have a separate association at each satisfaction threshold, the same source can be fitted with `nonparallel=~Cont`: ```{r} ordinal_partial_stream <- fastglm_ordinal_streaming( housing_source,nonparallel=~Cont,firth=TRUE) summary(ordinal_partial_stream) ``` The varying effects must produce ordered cumulative probabilities at the observed predictor values. Predictions check that condition again, since it need not hold when extrapolating beyond the observed predictor range. The [ordinal guide](ordinal-fastglm.html#ordinal-nonparallel) explains the model and its interpretation. Logit, probit, and complementary log-log links are supported with either maximum likelihood or mean bias reduction. For bias reduction, the information matrix is first accumulated over every block. A subsequent pass uses that global information to calculate the adjustment. Fitting separate bias-reduced models to the blocks would not produce the same estimator. Coefficient estimates, model-based covariance, Wald intervals, and numerical diagnostics use the ordinary ordinal interfaces. Predictions can be calculated at supplied predictor values, or the original source can be read again. When the complete prediction matrix would be too large, the `sink` argument receives one block at a time. As above, this function is the destination for the calculated values; it can write them to a file or retain only a summary. The following example retains the predicted number of residents in each satisfaction category for each block. ```{r} category_totals <- matrix(0,4L,nlevels(housing$Sat), dimnames=list(paste0("block",1:4),levels(housing$Sat))) predict(ordinal_stream,sink=function(value,k) { i <- ((k-1L)*18L+1L):(k*18L) category_totals[k,] <<- colSums(value*housing$Freq[i]) }) category_totals ``` Residuals and frequency weights can also be processed through a block destination. Returning all predictions, residuals, or weights in one object still requires memory proportional to the number of rows. Fitting storage depends on the largest supplied block, block-size records, and matrices whose dimensions depend on the number of coefficients. The source must itself read blocks without retaining the complete dataset for this memory benefit to apply. Peak process memory also includes temporary allocations and memory managed by R. Retaining a smaller design therefore does not guarantee a lower peak for every dataset size. Streaming can use more peak memory than dense fitting for moderate datasets and requires additional time to reread the observations. Its practical benefit should be assessed at the intended data size and with the intended reader, especially when the complete design would be too large to retain. The fitted source retains neither the complete design nor the response. Prediction scores, calibration, common-slope tests, and robust covariance can use the original source directly. Profile intervals, general constrained coefficient-score inference, and refitting bootstraps require retained observations. Average prediction contrasts can use explicitly supplied evaluation data. The source must continue to return the same rows and values whenever it is reread. ## Robust covariance for ordered outcomes `fastglm_vcov()` rereads an ordinal source once to accumulate the observed sensitivity of the ordinary scores and their variability. Its covariance agrees with the retained-data calculation, including for partial models and mean-bias-reduced estimates. The latter use the first-order ordinary-score convention explained in the [ordinal guide](ordinal-fastglm.html#ordinal-robust). ```{r ordinal-source-covariance} V_ordinal <- fastglm_vcov(ordinal_partial_stream,type="HC0") fastglm_average(ordinal_partial_stream,newdata=housing, weights=housing$Freq,vcov.=V_ordinal) ``` Cluster identifiers can describe observations distributed across several source blocks. The next example uses the synthetic predictor rows defined earlier, with an ordered outcome and 120 group identifiers. Repeating the identifiers across blocks illustrates how the score contributions are combined for a group. ```{r ordinal-source-clusters} set.seed(110305) ordered_data <- raw ordered_data$rating <- ordered(cut(.5*ordered_data$x+rlogis(n), c(-Inf,-.5,.75,Inf),labels=c("low","middle","high"))) ordinal_clusters <- rep(seq_len(120),length.out=n) ordinal_source <- fastglm_source( function(k) ordered_data[((k-1L)*chunk_rows+1L):(k*chunk_rows),], n_chunks=4L,formula=rating~x+g,prototype=ordered_data) ordinal_cluster_fit <- fastglm_ordinal_streaming( ordinal_source,nonparallel=~x,firth=TRUE) ordinal_cluster_ids <- function(k) ordinal_clusters[((k-1L)*chunk_rows+1L):(k*chunk_rows)] V_ordinal_cluster <- fastglm_vcov( ordinal_cluster_fit,type="HC1",cluster=ordinal_cluster_ids) summary(ordinal_cluster_fit,vcov.=V_ordinal_cluster) ``` The cluster function returns one identifier for each original row in the requested block. A label has the same meaning throughout the source; factor labels are matched by their values rather than their block-specific integer codes. All individuals represented by a frequency row must belong to its stated cluster. The uncertainty calculation requires sufficiently many independent clusters and an identifiable interior model. It does not add cluster-specific model parameters. Fitting, covariance calculation, and prediction must use unchanged source rows and values. Clustered covariance retains one score vector per distinct group, so its memory requirement grows with the number of groups and coefficients even though the complete design is not retained. ## Assessment and testing without collecting the observations {#ordinal-source-assessment} Prediction losses and calibration summarize how the fitted category probabilities relate to observed outcomes. These summaries require only accumulated totals, so `fastglm_ordinal_assess()` can compute them in one source pass. The default result contains the overall scores and calibration tables; it does not retain an entry for every input row. ```{r ordinal-source-assessment} ordinal_assessment <- fastglm_ordinal_assess(ordinal_partial_stream,bins=5) ordinal_assessment ``` These calculations use the fitting outcomes and therefore describe apparent performance. Separate evaluation data can still be supplied explicitly. To inspect individual losses, `keep.rows=TRUE` collects the row-score table in memory. Alternatively, a function supplied as `sink` receives each original block's row scores. This function is the output destination: it can write the result and release the block before the next is read. ```{r ordinal-source-assessment-output} score_directory <- tempfile("ordinal-scores-") dir.create(score_directory) save_scores <- function(value,k) saveRDS(value,file.path(score_directory,paste0("block-",k,".rds"))) fastglm_ordinal_assess(ordinal_partial_stream,sink=save_scores) ``` The saved tables include global row numbers, represented frequencies, prediction losses, and any exclusion reason. Retaining all row scores requires storage proportional to the number of observations; saving blocks to files moves that storage out of the fitted R object. The compact assessment summaries themselves depend on the number of outcome categories and calibration bins. ```{r ordinal-source-assessment-cleanup, include=FALSE} unlink(score_directory,recursive=TRUE) ``` Common-slope tests also operate on the original source. For the partial housing model, the following test asks whether the threshold-specific contact slopes can be replaced by one common contact effect: ```{r ordinal-source-po-test} fastglm_po_test(ordinal_partial_stream,terms="Cont") ``` The default score test refits the null by maximum likelihood. `null.estimation="meanBR"` requests the established ordinary-score calculation at a bias-reduced null estimate, and `method="lr"` compares full and restricted partial models fitted by maximum likelihood. The [ordinal guide](ordinal-fastglm.html#ordinal-partial-comparisons) explains their hypotheses and first-order interpretation. These tests need additional source passes because they refit models and calculate information across all blocks. No observation-level design or outcome array is collected. An optional diagnostic pass can locate rows with the smallest fitted probabilities or cumulative-predictor gaps: ```{r ordinal-source-boundary-diagnostics} fastglm_diagnostics(ordinal_partial_stream,boundary=TRUE,max.rows=3) ``` The printed rows include their source block and observed category frequency. Original predictor values remain available in `boundary.details`. Only a bounded number of extreme rows are retained, and small probabilities alone do not establish that a bias-reduced solution is absent. # Parquet and database sources Parquet is a file format that stores data by column and is commonly used for large tables. `fastglm_source_arrow()` reads one Parquet file as one block, so each file must be small enough to load into memory. A larger table can be divided into several files before fitting. For data in a database, `fastglm_source_dbi()` runs one supplied query per block. The queries should divide the observations without overlap and return them in a fixed order. The database contents used by the queries must remain unchanged across passes, for example by reading a fixed snapshot of the table. ```{r, eval=FALSE} parquet_source <- fastglm_source_arrow( files=c("part-1.parquet","part-2.parquet"), formula=y~x+g, prototype=prototype) fit <- fastglm_streaming(parquet_source, family=poisson()) queries <- c("SELECT y, x, g FROM outcomes WHERE id <= 10000 ORDER BY id", "SELECT y, x, g FROM outcomes WHERE id > 10000 ORDER BY id") database_source <- fastglm_source_dbi(connection,queries, formula=y~x+g,prototype=prototype) ``` These functions use [Arrow's Parquet reader](https://arrow.apache.org/docs/r/reference/read_parquet.html) and [DBI query retrieval](https://dbi.r-dbi.org/reference/dbGetQuery.html) to obtain each block. The analyst supplies the common model specification and a stable set of observations. This separation permits the same fitting routine to operate on files, databases, or other repeatable readers. # Count and two-part models from the same readers Negative-binomial, hurdle, and semicontinuous two-part models can also use repeatable sources. A two-component formula gives the count or positive-outcome predictors before `|` and the binary predictors after it. For example, the following hurdle model describes whether a count is positive and, separately, its distribution among positive counts. ```{r} hurdle_source <- fastglm_source(get_data,4, formula=y~x+g|x,prototype=raw) hurdle_stream <- fastglm_hurdle(hurdle_source,dist="poisson") fastglm_average(hurdle_stream) ``` `fastglm_nb(source)` estimates negative-binomial coefficients and theta from blocks, while `fastglm_twopart(source)` fits the binary and positive Gamma regressions. Their prediction functions accept `sink` functions with the same role described above: each output block is saved or summarized before the next is read. Robust covariance and average prediction contrasts also read blocks. Clustered inference retains a contribution vector for each group, and bootstrap resampling of rows or groups additionally requires row indices or counts proportional to the number of observations. The count and semicontinuous-outcome vignettes explain the quantities these models estimate and illustrate their uncertainty calculations.