causalreg

Note: This is a fork of the causalreg CRAN package by Vinciotti and Wit, refactored for code clarity and efficiency. The statistical method is unchanged.

Causal discovery in generalized linear models (GLMs) and generalized additive models (GAMs) via Pearson risk invariance.

Given a response variable and a set of candidate covariates, causalreg identifies the subset of covariates that are causal parents of the response within a structural causal model. The key idea is that the Pearson risk (expectation of squared Pearson residuals) equals 1 if and only if the model is correctly specified with respect to its causal parents.

Installation

Install the development version from GitHub:

# install.packages("pak")
# pak::pak("franciscorichter/causalreg")
install.packages("causalreg", repos = c('https://franciscorichter.r-universe.dev', 'https://cloud.r-project.org'))

# or, equivalently:
# install.packages("remotes")
# remotes::install_github("franciscorichter/causalreg")

Or install from CRAN:

install.packages("causalreg")

Quick start

Causal Poisson GLM

library(causalreg)

# Simulate data: X1 is a cause of Y, X2 is an effect
n <- 1000
set.seed(123)
X1 <- rnorm(n)
Y  <- rpois(n, exp(X1))
X2 <- log(Y + 1) + rnorm(n, 0, 0.3)
data <- data.frame(X1, X2, Y)

# Exhaustive search with chi-square test
result <- cglm(Y ~ X1 + X2, "poisson", data, pval = "chi-square", search = "all")
result$model.opt
#> [1] "Y ~ X1"

Causal logistic GLM

n <- 2000
set.seed(123)
X1 <- rnorm(n)
Y  <- rbinom(n, 1, exp(X1) / (1 + exp(X1)))
flip <- rbinom(n, 1, 0.1)
X2 <- (1 - flip) * Y + rnorm(n, 0, 0.3)
data <- data.frame(X1, X2, Y)

# Stepwise search with bootstrap test
set.seed(1)
result <- cglm(Y ~ X1 + X2, "binomial", data, pval = "bootstrap", search = "stepwise")
result$model.opt
#> [1] "Y ~ X1"

Causal Poisson GAM (nonlinear effects)

n <- 1000
set.seed(123)
X1 <- rnorm(n)
Y  <- rpois(n, exp(sin(X1)))
X2 <- log(Y + 1) + rnorm(n, 0, 0.5)
data <- data.frame(X1, X2, Y)

result <- cgam(Y ~ s(X1) + s(X2), "poisson", data, pval = "chi-square", search = "all")
result$model.opt
#> [1] "Y ~ s(X1)"

Main functions

Function Description
cglm() Causal discovery within a generalized linear model
cgam() Causal discovery within a generalized additive model

Both functions share the same interface:

Return values

Exhaustive search (search = "all")

Element Description
$model.opt Formula string of the selected causal model
$models List of all evaluated model formula strings
$pv P-values for the Pearson risk test (one per model)
$bic BIC values (one per model)
$pearsonrisk Pearson risk values (one per model)

Stepwise search (search = "stepwise")

Element Description
$model.opt Formula string of the selected causal model
$models List of model formula strings visited during the search

How it works

The method rests on a population-level identity. For a response $Y$ and a set of candidate parents $\mathbf{X}$, the Pearson risk is the expected squared Pearson residual

$$R_P = \mathbb{E}\left[\frac{\bigl(Y - \mu(\mathbf{X})\bigr)^2}{V\bigl(\mu(\mathbf{X})\bigr)}\right],$$

where $\mu(\mathbf{X})$ is the conditional mean and $V(\cdot)$ is the family's variance function. When the model is correctly specified with respect to the true causal parents, the conditional variance of $Y$ equals $V(\mu(\mathbf{X}))$, so each squared Pearson residual has expectation 1 and therefore

$$R_P = 1 \qquad \text{(exactly, for the causal model).}$$

If the model includes non-causal variables or omits causal ones, the mean or variance structure is misspecified and $R_P \neq 1$.

Empirically, $R_P$ is estimated by the average of the squared Pearson residuals,

$$\widehat{R}P = \frac{1}{n}\sum{i=1}^{n} \widehat{r}_i^{,2}, \qquad \widehat{r}_i = \frac{Y_i - \widehat{\mu}(\mathbf{X}_i)}{\sqrt{V\bigl(\widehat{\mu}(\mathbf{X}_i)\bigr)}},$$

and for each candidate subset the algorithm performs a statistical test of the null hypothesis $H_0\colon R_P = 1$ using either:

Among all subsets that pass the test ($p > \alpha$), the one with the lowest BIC is selected.

Choosing between search strategies

Choosing between p-value methods

Citation

If you use this package, please cite the paper describing the method:

Polinelli A, Vinciotti V, Wit EC (2026). "Causal generalized linear models via Pearson risk invariance." Journal of Causal Inference, 14(1), 20240043. https://doi.org/10.1515/jci-2024-0043

@article{polinelli2026causal,
  title   = {Causal generalized linear models via {Pearson} risk invariance},
  author  = {Polinelli, A. and Vinciotti, Veronica and Wit, Ernst C.},
  journal = {Journal of Causal Inference},
  year    = {2026},
  volume  = {14},
  number  = {1},
  pages   = {20240043},
  doi     = {10.1515/jci-2024-0043},
}

Changelog

v0.2.1 — Faster GAM bootstrap (June 2026)

The v0.2.0 C++ acceleration covered only cglm(). cgam() with pval = "bootstrap" stayed on pure R/mgcv and could be very slow, because every bootstrap resample re-ran mgcv's smoothing-parameter (REML/GCV) selection. An exhaustive search over 5 covariates (31 models, B = 100) on n = 5000 binomial data took roughly 43 minutes.

Timings on the n = 5000, 5-covariate, B = 100 example (Apple M1 Max); all three configurations select the same model Y ~ s(X2) + s(X3) + s(X5):

Configuration Time Speedup
ncores = 1, fast_gam = FALSE (v0.2.0 behavior) ~43 min 1.0×
ncores = 8, fast_gam = FALSE (exact, parallel) 20.0 min 2.1×
ncores = 8, fast_gam = TRUE (parallel + fixed sp) 7.7 min 5.6× 
# fast_gam = TRUE is the default; add ncores for large searches
cgam(fml, "binomial", data, pval = "bootstrap", search = "all", ncores = 8)

# exact (re-selected) bootstrap, matching v0.2.0 — slower
cgam(fml, "binomial", data, pval = "bootstrap", search = "all",
     ncores = 8, fast_gam = FALSE)

v0.2.0 — C++ acceleration via Rcpp (March 2026)

C++ fast path for GLM fitting (Rcpp/RcppArmadillo)

Benchmark

v0.1.2-fork — Refactored fork vs CRAN v0.1.2

The original CRAN package by Vinciotti and Wit (source) contained four near-identical internal files (cglm_all, cgam_all, cglm_step, cgam_step) with substantial code duplication. This fork refactors the codebase for clarity and efficiency while preserving the statistical method exactly.

Code cleanup and deduplication

Bug fixes

Computational efficiency

Packaging

C++ Acceleration (Rcpp)

Last updated: March 2026

Core GLM fitting, Pearson statistic computation, and bootstrap resampling have been reimplemented in C++ via RcppArmadillo. The fast path activates automatically for cglm() with "poisson" or "binomial" families. Set use_cpp = FALSE to use the original pure-R implementation. The C++ path does not cover cgam() (GAMs are fit with R/mgcv); for fast GAM bootstrap searches use fast_gam = TRUE and ncores instead (see v0.2.1 in the changelog).

Benchmark results (median wall-clock time, Apple Silicon ARM64, R 4.4.2):

Scenario R (ms) C++ (ms) Speedup
Poisson, χ², all, p=2, n=1,000 5.7 2.5 2.3×
Poisson, χ², all, p=5, n=1,000 51.8 14.9 3.5×
Poisson, χ², step, p=5, n=1,000 49.2 34.1 1.4×
Binomial, boot B=50, all, p=2, n=2,000 399.3 45.7 8.7×
Binomial, boot B=50, all, p=5, n=3,000 6,538.0 916.2 7.1×
Binomial, boot B=50, step, p=5, n=3,000 3,007.7 441.8 6.8×
Poisson, χ², all, p=3, n=10,000 80.3 25.3 3.2×

License

GPL-3