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.

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)))
#> [1] 4.440892e-16

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.

Xb <- bigmemory::as.big.matrix(X)
mapped_fit <- fastglm(Xb,raw$y,poisson(),method="llt")
max(abs(coef(mapped_fit)-coef(dense_fit)))
#> [1] 2.63678e-16

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.

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)))
#> [1] 5.551115e-16
c(passes=stream_fit$passes, chunks=stream_fit$chunk.reads)
#> passes chunks 
#>      9     36
summary(stream_fit)
#> Streaming generalized linear model
#> 
#> Call:
#> fastglm_streaming(chunk_callback = source, family = poisson())
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.23670    0.05091   4.650   <1e-04 ***
#> x            0.30854    0.02626  11.747   <1e-04 ***
#> gb          -0.08615    0.07286  -1.182    0.237    
#> gc          -0.07566    0.07255  -1.043    0.297    
#> gd          -0.04880    0.07152  -0.682    0.495    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>  Dispersion Deviance Residual df  AIC
#>           1     1404        1195 3410
#> Observations: 1200
#> Convergence: converged

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

qr_stream <- fastglm_streaming(source,family=poisson(),method="qr")
max(abs(coef(qr_stream)-coef(stream_fit)))
#> [1] 3.885781e-16
head(predict(qr_stream,type="response",se.fit=TRUE)$se.fit)
#> [1] 0.05558614 0.06720906 0.05796655 0.05398368 0.06837212 0.05648853

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.

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)))
#> [1] 1.045563
c(passes=bias_stream$passes,blocks_read=bias_stream$chunk.reads)
#>      passes blocks_read 
#>           7          28

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.

chunk_means <- numeric(4)
predict(stream_fit,type="response",
        sink=function(value,k) chunk_means[k] <<- mean(value))
chunk_means
#> [1] 1.259740 1.254676 1.267344 1.254906
fastglm_average(stream_fit)
#>   estimate  std.error conf.low conf.high
#> 1 1.259167 0.03239299 1.195678  1.322656

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.

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))
#> (Intercept)           x          gb          gc          gd 
#>  0.04634219  0.02516485  0.06741032  0.06989078  0.07090223
fastglm_average(stream_fit,vcov.=V)
#>   estimate  std.error conf.low conf.high
#> 1 1.259167 0.03475857 1.191041  1.327292

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.

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)
#> Proportional-odds regression
#> 
#> Call:
#> fastglm::fastglm_ordinal_streaming(source = housing_source, firth = TRUE)
#> Estimation: mean bias reduction
#> 
#> Location coefficients (positive means higher outcomes):
#>               Estimate Std. Error z value Pr(>|z|)    
#> InflMedium     0.56442    0.10496   5.378  < 1e-04 ***
#> InflHigh       1.28379    0.12667  10.135  < 1e-04 ***
#> TypeApartment -0.57013    0.11873  -4.802  < 1e-04 ***
#> TypeAtrium    -0.36499    0.15675  -2.328 0.019887 *  
#> TypeTerrace   -1.08673    0.15149  -7.174  < 1e-04 ***
#> ContHigh       0.35891    0.09535   3.764 0.000167 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Thresholds:
#>                       Estimate Std. Error
#> threshold:Low|Medium   -0.4943     0.1245
#> threshold:Medium|High   0.6883     0.1252
#> 
#> Uncertainty: expected Fisher information
#> Sample: 1681 observations represented by 72 positive-frequency rows
#> Convergence: converged (7 iterations)
fastglm_diagnostics(ordinal_stream)
#> Numerical diagnostics
#> Convergence: converged (7 iterations)
#> Estimable coefficients (rank): 8
#> Boundary flag: no
#> 
#>                                             Check     Value
#>  Largest absolute adjusted score (reported units) 9.665e-08
#>                       Standardized adjusted score 9.693e-09
#>              Remaining change in linear predictor 1.584e-09
#>           Reciprocal information condition number    0.1449
#> Score scaling: Raw scores use reported predictor units; standardized scores use
#>   the centered/scaled fitting coordinates.

To allow contact with other residents to have a separate association at each satisfaction threshold, the same source can be fitted with nonparallel=~Cont:

ordinal_partial_stream <- fastglm_ordinal_streaming(
    housing_source,nonparallel=~Cont,firth=TRUE)
summary(ordinal_partial_stream)
#> Partial proportional-odds regression
#> 
#> Call:
#> fastglm::fastglm_ordinal_streaming(source = housing_source, firth = TRUE, 
#>     nonparallel = ~Cont)
#> Estimation: mean bias reduction
#> 
#> Common slopes:
#>               Estimate Std. Error z value Pr(>|z|)    
#> InflMedium      0.5674     0.1050   5.403   <1e-04 ***
#> InflHigh        1.2832     0.1267  10.129   <1e-04 ***
#> TypeApartment  -0.5683     0.1188  -4.783   <1e-04 ***
#> TypeAtrium     -0.3631     0.1568  -2.316   0.0206 *  
#> TypeTerrace    -1.0935     0.1516  -7.214   <1e-04 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Threshold-specific slopes:
#>  Predictor   Threshold Estimate Std. Error z value Pr(>|z|)
#>   ContHigh  Low|Medium   0.4424     0.1092   4.049  < 1e-04
#>   ContHigh Medium|High   0.2848     0.1065   2.675 0.007477
#> 
#> Thresholds at zero covariates and offset:
#>                       Estimate Std. Error
#> threshold:Low|Medium   -0.4476     0.1281
#> threshold:Medium|High   0.6456     0.1279
#> 
#> Uncertainty: expected Fisher information
#> Sample: 1681 observations represented by 72 positive-frequency rows
#> Convergence: converged (8 iterations)

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

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
#>              Low    Medium      High
#> block1 178.28195 143.13541 214.58264
#> block2  77.02767  46.15468  53.81765
#> block3 182.85196 166.19114 279.95690
#> block4 129.26579  89.31968 120.41453

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.

V_ordinal <- fastglm_vcov(ordinal_partial_stream,type="HC0")
fastglm_average(ordinal_partial_stream,newdata=housing,
                weights=housing$Freq,vcov.=V_ordinal)
#>   category  estimate  std.error  conf.low conf.high
#> 1      Low 0.3373250 0.01103870 0.3156896 0.3589605
#> 2   Medium 0.2646364 0.01071245 0.2436404 0.2856324
#> 3     High 0.3980385 0.01139588 0.3757030 0.4203740

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)
#> Partial proportional-odds regression
#> 
#> Call:
#> fastglm::fastglm_ordinal_streaming(source = ordinal_source, firth = TRUE, 
#>     nonparallel = ~x)
#> Estimation: mean bias reduction
#> 
#> Common slopes:
#>    Estimate Std. Error z value Pr(>|z|)
#> gb  0.10356    0.16740   0.619    0.536
#> gc  0.07135    0.16289   0.438    0.661
#> gd  0.03407    0.16055   0.212    0.832
#> 
#> Threshold-specific slopes:
#>  Predictor   Threshold Estimate Std. Error z value Pr(>|z|)
#>          x  low|middle   0.4389     0.0560   7.837  < 1e-04
#>          x middle|high   0.4130     0.0664   6.219  < 1e-04
#> 
#> Thresholds at zero covariates and offset:
#>                       Estimate Std. Error
#> threshold:low|middle   -0.3709     0.1235
#> threshold:middle|high   0.8197     0.1236
#> 
#> Uncertainty: supplied covariance
#> Sample: 1200 observations represented by 1200 positive-frequency rows
#> Convergence: converged (7 iterations)

ordinal_assessment <- fastglm_ordinal_assess(ordinal_partial_stream,bins=5)
ordinal_assessment
#> Ordinal assessment (training source)
#> 
#> Prediction scores:
#>  Frequency Log loss Ranked probability score Normalized RPS
#>       1681    1.034                   0.4273         0.2136
#> Source blocks: 4
#> Rows assessed: 72
#> Row scores: Not retained; use keep.rows or sink to obtain them

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)
#> Ordinal assessment (training source)
#> 
#> Prediction scores:
#>  Frequency Log loss Ranked probability score Normalized RPS
#>       1681    1.034                   0.4273         0.2136
#> Source blocks: 4
#> Rows assessed: 72
#> Row scores: Supplied to callback

fastglm_po_test(ordinal_partial_stream,terms="Cont")
#> Partial proportional-odds expected-information score test
#> Terms: Cont
#> Null fit: maximum likelihood
#> Sample: 1681 observations represented by 72 rows
#> 
#>  Statistic df p-value
#>      2.435  1  0.1186
#> 
#> Note: Ordinary efficient score at the ML null; first-order chi-square
#>   approximation.

fastglm_diagnostics(ordinal_partial_stream,boundary=TRUE,max.rows=3)
#> Numerical diagnostics
#> Convergence: converged (8 iterations)
#> Estimable coefficients (rank): 9
#> Boundary flag: no
#> 
#>                                             Check     Value
#>  Largest absolute adjusted score (reported units) 3.178e-08
#>                       Standardized adjusted score 2.126e-09
#>              Remaining change in linear predictor 2.266e-10
#>           Reciprocal information condition number    0.1375
#>          Smallest fitted cumulative-predictor gap     1.093
#> Score scaling: Raw scores use reported predictor units; standardized scores use
#>   the centered/scaled fitting coordinates.
#> 
#> Rows near a probability boundary:
#> Positive-frequency rows examined: 72
#> Additional source passes: 1
#> 
#> Smallest fitted category probabilities:
#>  Row Block Category Probability Log probability Observed frequency
#>   43     3      Low      0.1022          -2.281                  3
#>   44     3      Low      0.1022          -2.281                  0
#>   45     3      Low      0.1022          -2.281                  0
#> 
#> Smallest gaps between adjacent cumulative predictors:
#>  Row Block Category Predictor gap Observed frequency
#>    1     1   Medium         1.093                  0
#>    2     1   Medium         1.093                 21
#>    3     1   Medium         1.093                  0
#> Interpretation: Small probabilities or gaps do not establish that an interior
#>   solution is absent. A displayed zero can be numerical underflow; inspect its
#>   log probability.

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.

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 and DBI query retrieval 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.

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)
#>   estimate  std.error conf.low conf.high
#> 1 1.258903 0.03303574 1.194155  1.323652

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.