An end-to-end real-data analysis with fbrglm

1. Overview

The other vignettes describe fbrglm’s API surface (Getting started) and walk each supported family on synthetic data (Families and model types). This vignette does one thing end-to-end on a real dataset:

  • read a real, small, publicly-available dataset;
  • split it into training and test sets;
  • fit a regularized logistic model with fbrglm();
  • predict on the test set and evaluate;
  • inspect and reuse the fit.

fbrglm does not replace glmnet’s optimizer — the numerical fit is still glmnet::glmnet() / glmnet::cv.glmnet(). What fbrglm adds is a glm()-shaped formula/data interface that stores the training-time model recipe (terms, xlevels, contrasts) on the fit object, restores glm()’s rank-deficient NA convention, and attaches complete-case bookkeeping (nobs_info) to the returned model. Those additions are what a working analyst notices in an end-to-end pipeline; the numeric coefficients are the same as raw glmnet given the same data and the same lambda.

2. Data

We use MASS::birthwt, a classic obstetric cohort collected at Baystate Medical Center, Springfield, MA in 1986. Each row is one mother; the response low indicates a low-birth-weight birth (< 2500 g). The dataset ships with the MASS package, so no download is needed.

library("fbrglm")
library("MASS")

data("birthwt", package = "MASS")

dim(birthwt)
#> [1] 189  10
head(birthwt, 3)
#>    low age lwt race smoke ptl ht ui ftv  bwt
#> 85   0  19 182    2     0   0  0  1   0 2523
#> 86   0  33 155    3     0   0  0  0   3 2551
#> 87   0  20 105    1     1   0  0  0   1 2557

The integer columns for race, smoking status, and first-trimester visits are more readable as factors, so we convert them up front. We also add a derived indicator any_visit. This one is deliberate: it is a linear function of the ftv dummies and the intercept, and would trip up a raw glmnet pipeline that silently accepts collinear columns. It appears often in real ETL — a downstream “cleaner” column bundled alongside the raw feature.

birthwt$race  <- factor(birthwt$race,  1:3, c("white", "black", "other"))
birthwt$smoke <- factor(birthwt$smoke, 0:1, c("no", "yes"))
birthwt$ftv   <- factor(birthwt$ftv)
birthwt$any_visit <- as.integer(birthwt$ftv != "0")

table(low = birthwt$low)                 # outcome distribution
#> low
#>   0   1 
#> 130  59
sapply(birthwt, function(v) sum(is.na(v)))  # no missing values here
#>       low       age       lwt      race     smoke       ptl        ht        ui 
#>         0         0         0         0         0         0         0         0 
#>       ftv       bwt any_visit 
#>         0         0         0

3. Train / test split

A 150 / 39 split with a pinned seed. The split is not stratified on race; that’s on purpose — it lets the test batch happen to miss one of the three race levels, which is exactly the kind of drift a real scoring pipeline sees.

set.seed(93)
idx   <- sample.int(nrow(birthwt), 150)
train <- birthwt[idx, ]
test  <- birthwt[-idx, ]

c(train_n = nrow(train), test_n = nrow(test))
#> train_n  test_n 
#>     150      39
table(train$low)
#> 
#>   0   1 
#> 100  50
table(test$low)
#> 
#>  0  1 
#> 30  9

## The test batch, by chance, contains no rows at race == "black".
## This is the kind of factor drift that crashes a naive glmnet caller.
table(test$race)
#> 
#> white black other 
#>    19     0    20

Nothing about the fit needs to be done differently to defend against this drift — fbrglm freezes the fit-time xlevels on the fit object and reuses them at predict time.

4. Preprocessing

fbrglm takes a formula + data.frame and builds the numeric design matrix internally. That means the preprocessing you actually have to write is limited to what you want on the data frame (factor recoding, column selection, complete-case filtering if desired), rather than everything that turns a data frame into a numeric matrix.

predictors <- c("age", "lwt", "race", "smoke", "ptl", "ftv", "any_visit")
str(train[, c("low", predictors)])
#> 'data.frame':    150 obs. of  8 variables:
#>  $ low      : int  1 1 1 0 1 0 0 0 1 0 ...
#>  $ age      : int  22 24 29 36 26 30 21 31 17 25 ...
#>  $ lwt      : int  130 105 130 202 190 112 124 100 120 130 ...
#>  $ race     : Factor w/ 3 levels "white","black",..: 1 2 1 1 1 1 3 1 1 1 ...
#>  $ smoke    : Factor w/ 2 levels "no","yes": 2 2 1 1 2 1 1 1 2 1 ...
#>  $ ptl      : int  0 0 0 0 0 0 0 0 0 0 ...
#>  $ ftv      : Factor w/ 6 levels "0","1","2","3",..: 2 1 3 2 1 2 1 4 4 3 ...
#>  $ any_visit: int  1 0 1 1 0 1 0 1 1 1 ...

No further transformation is needed. Concretely, we are not calling model.matrix() at any point in this script; that step happens inside fbrglm() at fit time and again at predict time (against the frozen xlevels), so the two design matrices are guaranteed to agree.

5. Fit with fbrglm

We select lambda by cross-validation (lambda = "cv_min") with 5 folds to keep the vignette runtime short. In an applied setting you would typically use the default 10 folds.

fit <- fbrglm(
    low ~ age + lwt + race + smoke + ptl + ftv + any_visit,
    data   = train,
    family = "binomial",
    lambda = "cv_min",
    nfolds = 5
)
#> fbrglm: dropped 1 linearly dependent column(s) (rank 11 / 12): any_visit. Their coefficients will be reported as NA, as in stats::glm().

print() shows the call, family, chosen lambda, and the complete-case bookkeeping in one glance:

print(fit)
#> Formula-based regularized GLM (fbrglm)
#> 
#> Call:
#> fbrglm(formula = low ~ age + lwt + race + smoke + ptl + ftv + 
#>     any_visit, data = train, family = "binomial", lambda = "cv_min", 
#>     nfolds = 5)
#> 
#> Family:        binomial
#> Alpha:         1
#> Lambda rule:   cv_min
#> Lambda value:  0.01612082
#> Inference:     none
#> Observations:  total = 150, dropped = 0, used = 150
#> Non-zero terms: 9 / 12
#> Rank deficiency: 1 / 12 columns dropped (NA in coef): any_visit

summary() prints the same header plus the (regularized) coefficient vector. any_visit is reported as NA: fbrglm caught the linear dependence at fit time (QR-pivot on the column-centered design matrix) and preserved stats::glm()’s convention. A numeric zero elsewhere in the coefficient table therefore unambiguously means “shrunk to zero by the penalty” rather than “not identifiable in the design”.

summary(fit)
#> fbrglm summary
#> ==============
#> 
#> Call:
#> fbrglm(formula = low ~ age + lwt + race + smoke + ptl + ftv + 
#>     any_visit, data = train, family = "binomial", lambda = "cv_min", 
#>     nfolds = 5)
#> 
#> Family:        binomial
#> Lambda method: cv_min
#> Lambda value:  0.01612082
#> Inference:     none
#> 
#> Observations:
#>   total = 150, dropped (missing) = 0, used = 150
#> 
#> Coefficients: (1 not defined because of singularities: any_visit)
#>             Estimate
#> (Intercept)   0.0411
#> age          -0.0180
#> lwt          -0.0082
#> raceblack     0.6548
#> raceother     0.5022
#> smokeyes      0.9401
#> ptl           0.5126
#> ftv1         -0.1931
#> ftv2          0.0000
#> ftv3          0.2935
#> ftv4          0.0000
#> ftv6         -0.3270
#> any_visit         NA
#> 
#> Non-zero predictors (9): age, lwt, raceblack, raceother, smokeyes, ptl, ftv1, ftv3, ftv6
#> 
#> Note: no standard errors, z-values, or p-values are reported under
#> infer = "none". Classical inference does not account for the
#> shrinkage bias from L1/L2 penalisation or for the data-driven
#> selection of lambda. Use infer = "split" / "selective" (planned)
#> for valid post-selection inference; see vignette("fbrglm").

Note that no standard errors, z-values, p-values, or CIs are printed: regularized point estimates do not have a valid classical Wald sampling distribution, so fbrglm suppresses them by default under infer = "none".

Two accessors expose the state a downstream analyst usually wants:

nobs(fit)
#> [1] 150
fit$nobs_info
#> $n_total
#> [1] 150
#> 
#> $n_dropped_missing
#> [1] 0
#> 
#> $n_used
#> [1] 150
fit$rank_info[c("rank", "ncol", "rank_deficient", "dropped_cols")]
#> $rank
#> [1] 11
#> 
#> $ncol
#> [1] 12
#> 
#> $rank_deficient
#> [1] TRUE
#> 
#> $dropped_cols
#> [1] "any_visit"

tidy() returns the same information as a tibble for downstream data-wrangling:

tidy(fit)
#> # A tibble: 13 × 4
#>    term        estimate nonzero lambda
#>    <chr>          <dbl> <lgl>    <dbl>
#>  1 (Intercept)  0.0411  TRUE    0.0161
#>  2 age         -0.0180  TRUE    0.0161
#>  3 lwt         -0.00824 TRUE    0.0161
#>  4 raceblack    0.655   TRUE    0.0161
#>  5 raceother    0.502   TRUE    0.0161
#>  6 smokeyes     0.940   TRUE    0.0161
#>  7 ptl          0.513   TRUE    0.0161
#>  8 ftv1        -0.193   TRUE    0.0161
#>  9 ftv2         0       FALSE   0.0161
#> 10 ftv3         0.293   TRUE    0.0161
#> 11 ftv4         0       FALSE   0.0161
#> 12 ftv6        -0.327   TRUE    0.0161
#> 13 any_visit   NA       FALSE   0.0161

6. Prediction

The test batch we set aside contains only two of the three training race levels (table(test$race) from Section 3). We predict on it directly — no manual model.matrix() alignment is needed.

p_test <- predict(fit, newdata = test, type = "response")
summary(p_test)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.1259  0.2682  0.3376  0.3464  0.4208  0.7202
head(cbind(p_test, low = test$low), 5)
#>         p_test low
#> [1,] 0.3091704   0
#> [2,] 0.4096995   0
#> [3,] 0.4585141   0
#> [4,] 0.4585141   0
#> [5,] 0.4183535   0

Under a hand-built glmnet call, the same input would raise "The number of variables in newx must be N" because model.matrix(~ ., test) would produce fewer columns than the training matrix. fbrglm reconstructs the design matrix against the stored xlevels, so the training-time layout is preserved regardless of which levels happen to appear at scoring time.

If a scoring batch contains a truly unseen factor level (rather than a subset of the training levels, as here), the default behavior is to raise a checked error. Production scoring pipelines that must not halt on a single unseen row can opt into predict(fit, newdata, on_new_levels = "na"), which returns NA for the affected rows and continues.

7. Evaluation

We report three metrics — accuracy at the 0.5 threshold, log-loss, and ROC AUC — all computed in base R to keep the vignette dependency footprint minimal.

y_test <- test$low
p_vec  <- as.numeric(p_test)

## Accuracy at threshold 0.5
acc <- mean((p_vec > 0.5) == (y_test == 1))

## Log-loss
log_loss <- -mean(y_test * log(pmax(p_vec, 1e-12)) +
                  (1 - y_test) * log(pmax(1 - p_vec, 1e-12)))

## ROC AUC via the trapezoidal rule
auc <- {
    ord  <- order(p_vec, decreasing = TRUE)
    y_o  <- y_test[ord]
    tp   <- cumsum(y_o == 1) / max(1, sum(y_test == 1))
    fp   <- cumsum(y_o == 0) / max(1, sum(y_test == 0))
    sum(diff(c(0, fp)) * (c(0, tp[-length(tp)]) + tp) / 2)
}

round(c(accuracy = acc, log_loss = log_loss, auc = auc), 3)
#> accuracy log_loss      auc 
#>    0.744    0.564    0.567

These numbers are modest; birthwt is a small cohort with weak signals and the point of the vignette is the workflow, not the fit quality. In an applied study you would compare several models on these same metrics; the pattern of fbrglm() calls stays identical.

8. Visualization

Two quick views. First, the distribution of predicted probabilities on the test set, split by the observed outcome — a sanity check that the model is separating the two classes at all:

op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
hist(p_vec[y_test == 0], breaks = seq(0, 1, by = 0.05),
     main = "low = 0 (normal weight)", xlab = "Predicted probability",
     col = "grey85", border = "white")
hist(p_vec[y_test == 1], breaks = seq(0, 1, by = 0.05),
     main = "low = 1 (low birth weight)", xlab = "Predicted probability",
     col = "grey60", border = "white")

par(op)

Second, the standard glm-style diagnostic panel produced by plot(fit). fbrglm’s deviance residuals are computed against the family’s dev.resids, so the two-arc pattern typical of a binomial GLM is visible. Raw glmnet does not expose residuals in the glm() sense, so this view is a genuine fbrglm add-on rather than a pass-through.

op <- par(mfrow = c(1, 3), mar = c(3.5, 3.5, 2, 0.5), mgp = c(2, 0.7, 0))
plot(fit)

par(op)

9. Inspecting and reusing the fit

as_glmnet() and as_cv_glmnet() return the underlying objects unchanged, so any function that consumes a raw glmnet handle keeps working:

class(as_glmnet(fit))
#> [1] "lognet" "glmnet"
class(as_cv_glmnet(fit))
#> [1] "cv.glmnet"

The fit object survives saveRDS() / readRDS() without losing the training-time model recipe, which is what makes hand-off between an exploration session and a scoring service safe:

tmp <- tempfile(fileext = ".rds")
saveRDS(fit, tmp)
fit_reloaded <- readRDS(tmp)
p_reloaded   <- predict(fit_reloaded, newdata = test, type = "response")

## Predictions on newdata are identical after round-trip.
isTRUE(all.equal(as.numeric(p_test), as.numeric(p_reloaded)))
#> [1] TRUE
unlink(tmp)

The frozen terms / xlevels / contrasts on the fit object are what makes this round-trip safe: predict() on fit_reloaded never consults the current session’s options("contrasts") or the levels present in test, and so cannot drift from the training-time encoding.

10. What fbrglm saved us from

A hand-rolled equivalent of the pipeline above would have to do, outside the model:

  • run model.matrix(~ ., train) and strip the intercept, in exactly the shape glmnet::glmnet() expects;
  • do the same at predict time, against a newdata whose factor levels may have narrowed — and remember to pad the missing dummy columns with zeros before calling predict();
  • keep training-time xlevels and contrasts in a separate side variable, otherwise a different R session with a different options("contrasts") will silently change the design matrix’s meaning;
  • run a rank check on the design matrix to catch the accidental linear dependency (any_visit in this example) — otherwise the glmnet coefficient for the dependent column comes back as a numeric 0, visually indistinguishable from a penalty result;
  • track how many rows were dropped by complete-case filtering and attach that count to the model somewhere the future analyst can find.

None of these steps involve any statistical decision. Moving them into the package leaves the analysis script to say formula + data + family and read back a fit object that behaves like stats::glm() for every downstream question.

11. Reproducibility

  • Data: MASS::birthwt, shipped with the MASS package. No external download.
  • Seed: set.seed(93) at the top of Section 3 pins the train / test split.
  • Software used to build this vignette:
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] MASS_7.3-66    glmnet_5.1     Matrix_1.7-6   fbrglm_0.1.0   rmarkdown_2.32
#> 
#> loaded via a namespace (and not attached):
#>  [1] jsonlite_2.0.0   compiler_4.6.1   Rcpp_1.1.2       jquerylib_0.1.4 
#>  [5] splines_4.6.1    yaml_2.3.12      fastmap_1.2.0    lattice_0.23-1  
#>  [9] R6_2.6.1         generics_0.1.4   shape_1.4.6.1    knitr_1.52      
#> [13] iterators_1.0.14 tibble_3.3.1     maketools_1.3.2  bslib_0.12.0    
#> [17] pillar_1.11.1    rlang_1.3.0      utf8_1.2.6       cachem_1.1.0    
#> [21] xfun_0.61        sass_0.4.10      sys_3.4.3        otel_0.2.0      
#> [25] cli_3.6.6        magrittr_2.0.5   digest_0.6.39    foreach_1.5.2   
#> [29] grid_4.6.1       lifecycle_1.0.5  vctrs_0.7.3      evaluate_1.0.5  
#> [33] glue_1.8.1       codetools_0.2-20 buildtools_1.0.0 survival_3.8-12 
#> [37] tools_4.6.1      pkgconfig_2.0.3  htmltools_0.5.9