Introduction to causalreg

Overview

The causalreg package implements causal discovery via Pearson risk invariance for generalized linear models (GLMs) and generalized additive models (GAMs). Given a response and a set of candidate covariates, it identifies which covariates are causal parents of the response within a structural causal model.

The key idea: if a GLM or GAM is correctly specified with respect to the true causal parents, the Pearson risk (expected squared Pearson residuals) is equal to 1. The package performs this test from observational data across subsets of covariates to find the causal model.

Causal Poisson regression

Consider a simple setting where X1 causes Y (Poisson), and X2 is a downstream effect of Y:

library(causalreg)

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)

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

# All models considered
unlist(result$models)
#> [1] "Y ~ X1"      "Y ~ X2"      "Y ~ X1 + X2"

# Their p-values (acceptance means no evidence to reject Pearson risk = 1)
result$pv
#> [1] 1.155966e-01 1.143346e-39 3.205287e-57

# Their BIC values
result$bic
#> [1] 2725.446 2245.723 2193.793

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

# Models visited during the search
unlist(result_step$models)
#> [1] "Y ~ 1"  "Y ~ X1"

Causal logistic regression

For binomial models, the chi-square approximation does not hold, so we use the bootstrap test:

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)

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

The bootstrap test correctly identifies the causal model.

Causal GAMs for nonlinear relationships

When the relationship between covariates and response is nonlinear, use cgam() with smooth terms:

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

The cgam() function uses mgcv::gam() internally and handles smooth terms (s()) in the formula. The Pearson risk invariance principle applies equally to GAMs.

Larger example: 5 covariates

In a more realistic setting with 5 candidate covariates and a binomial response:

set.seed(12)
n <- 3000
X1 <- rnorm(n)
X2 <- rnorm(n, X1, 0.5)
X3 <- rnorm(n, 0, 1)
X4 <- rnorm(n, X2, 0.5)
Y  <- rbinom(n, 1, exp(0.8 * X2 - 0.9 * X3) / (1 + exp(0.8 * X2 - 0.9 * X3)))
flip <- rbinom(n, 1, 0.1)
X5 <- (1 - flip) * Y + flip * (1 - Y) + rnorm(n, 0, 0.3)
dat <- data.frame(X1, X2, X3, X4, X5, Y)

# Exhaustive search (evaluates all 2^5 - 1 = 31 subsets)
set.seed(1)
mod_all <- cglm(Y ~ X1 + X2 + X3 + X4 + X5, "binomial", dat,
                pval = "bootstrap", search = "all")
mod_all$model.opt
#> [1] "Y ~ X2 + X3"

# Stepwise search (much faster)
set.seed(1)
mod_step <- cglm(Y ~ X1 + X2 + X3 + X4 + X5, "binomial", dat,
                 pval = "bootstrap", search = "stepwise")
mod_step$model.opt
#> [1] "Y ~ X2 + X3"

Both search strategies correctly identify X2 and X3 as the causal parents of Y.

Summary of recommendations

Reference

Polinelli, A., V. Vinciotti and E.C. Wit. (2026). “Causal generalized linear models via Pearson risk invariance.” Journal of Causal Inference.