Package 'causalreg' reference manual

Title:	Causal Generalized Linear Models
Description:	An implementation of methods for causal discovery in a structural causal model where the conditional distribution of the target node is described by a generalized linear model conditional on its causal parents.
Authors:	Veronica Vinciotti [aut, cre], Ernst C. Wit [aut], Francisco Richter [aut]
Maintainer:	Veronica Vinciotti <[email protected]>
License:	GPL-3
Version:	0.3.0
Built:	2026-10-02 13:30:25 UTC
Source:	https://github.com/franciscorichter/causalreg

Bootstrap p-value for Pearson risk = 1 test (fully in C++)

Description

Bootstrap p-value for Pearson risk = 1 test (fully in C++)

Usage

boot_pval_cpp(X, y, family, B = 100L)
boot_pval_cpp(X, y, family, B = 100L)

Arguments

X

Design matrix

y

Response vector

family

"poisson" or "binomial"

B

Number of bootstrap replicates

Value

Two-sided bootstrap p-value

Causal generalized additive model

Description

This function does a search for a causal submodel within the generalized additive model provided.

Usage

cgam(
  formula,
  family,
  data,
  alpha = 0.05,
  pval = c("bootstrap", "chi-square"),
  B = 100,
  search = c("all", "stepwise"),
  ncores = 1L,
  fast_gam = TRUE,
  direction = c("forward", "backward"),
  ...
)
cgam(
  formula,
  family,
  data,
  alpha = 0.05,
  pval = c("bootstrap", "chi-square"),
  B = 100,
  search = c("all", "stepwise"),
  ncores = 1L,
  fast_gam = TRUE,
  direction = c("forward", "backward"),
  ...
)

Arguments

formula

A formula object.

family

The response family as a character string, "poisson" or "binomial". A family object such as poisson() is not accepted.

data

A data frame containing the variables in the model.

alpha

Significance level for statistical test.

pval

If pval="bootstrap", a bootstrap test is conducted to test whether Pearson risk is 1. When family="poisson" a chi-squared test can be conducted by setting pval="chi-square".

B

Number of bootstrap samples when pval="bootstrap". Default is 100.

search

"all" (the default) fits every non-empty submodel of the candidate terms; "stepwise" runs a stepwise search in the direction set by direction.

ncores

Number of cores for parallel computation. Default is 1 (sequential). When ncores > 1, model evaluations are distributed across cores using forking (parallel::mclapply) on Unix/macOS and a PSOCK cluster (parallel::parLapply) on Windows, so parallelization is available on all platforms. The backend can be overridden with options(causalreg.parallel = ) ("auto", "fork", "psock", or "sequential").

fast_gam

Logical; only relevant when pval="bootstrap". If TRUE, the smoothing parameters selected on the original data (once per candidate model) are held fixed across that model's bootstrap resamples, so each resample fit skips mgcv's REML/GCV smoothing-parameter selection. This typically gives a 3-4x speedup (composable with ncores) and the same model selection, but it is an approximation: individual bootstrap p-values can differ from the fully re-selected version, more so at small B. Default is TRUE; set fast_gam=FALSE to re-select the smoothing parameters on every resample (the exact bootstrap used in earlier versions, slower). Has no effect on cglm or on the chi-square test.

direction

Direction of the stepwise search, '"forward"' (default) or '"backward"'. Ignored when 'search = "all"'. Forward starts from the intercept model and adds the variable giving the largest p-value; backward starts from the full model and removes the variable whose removal gives the largest p-value, stopping once the model is no longer rejected. Both then prune by BIC. For a binomial response whose candidate set contains a categorical variable, the backward search is used and a message is emitted, because the Pearson risk of a binary regression on categorical covariates alone is mathematically equal to 1.

...

Further arguments to be passed to the gam function.

Value

A list containing the selected causal submodel and search diagnostics.

References

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

Examples

##############################
#causal Poisson gam##########
n<-1000
set.seed(123)
X1<-rnorm(n,0,1)
Y<-rpois(n,exp(sin(X1)))
X2<-log(Y+1)+rnorm(n,0,0.5)
data<-data.frame(X1, X2, Y)
cm_all<-cgam(Y ~ s(X1)+s(X2),"poisson",data,pval="chi-square",search="all")
cm_all$model.opt
cm_step<-cgam(Y ~ s(X1)+s(X2),"poisson",data,pval="chi-square",search="stepwise")
cm_step$model.opt

#bigger simulation with 7 covariates
set.seed(123)
n<-1000
X1<-rnorm(n=n,sd=sqrt(0.04))
X2<-X1+rnorm(n=n,sd=sqrt(0.04))
X3<-X1+X2+rnorm(n=n,sd=sqrt(0.04))
m<-sin(X2*5)+X3^3
Z<-m+rnorm(n=n,sd=sqrt(0.04))
X4<-X2+rnorm(n=n,sd=sqrt(0.04))
X5<-Z+rnorm(n=n,sd=sqrt(0.04))
X6<-Z+rnorm(n=n,sd=sqrt(0.04))
X7<-X6+rnorm(n=n,sd=sqrt(0.04))
Y<-qpois(pnorm(Z, mean = m, sd = sqrt(0.04)), lambda=exp(m))
dat<-data.frame(X1, X2, X3, X4, X5, X6, X7,Y)
fml<- Y~s(X1)+s(X2)+s(X3)+s(X4)+s(X5)+s(X6)+s(X7)
mod.all <-cgam(fml,"poisson",dat,pval="chi-square",search="all")
mod.all$model.opt
mod.step <-cgam(fml,"poisson",dat,pval="chi-square",search="stepwise")
mod.step$model.opt
####################################
#causal logistic gam################
n<-1000
set.seed(123)
X1<-rnorm(n,0,1)
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)
cm_all<-cgam(Y ~ s(X1)+s(X2),"binomial",data,pval="bootstrap",search="all")
cm_all$model.opt
set.seed(1)
cm_step<-cgam(Y ~ s(X1)+s(X2),"binomial",data,pval="bootstrap",search="stepwise")
cm_step$model.opt

##############################
#causal Poisson gam##########
n<-1000
set.seed(123)
X1<-rnorm(n,0,1)
Y<-rpois(n,exp(sin(X1)))
X2<-log(Y+1)+rnorm(n,0,0.5)
data<-data.frame(X1, X2, Y)
cm_all<-cgam(Y ~ s(X1)+s(X2),"poisson",data,pval="chi-square",search="all")
cm_all$model.opt
cm_step<-cgam(Y ~ s(X1)+s(X2),"poisson",data,pval="chi-square",search="stepwise")
cm_step$model.opt

#bigger simulation with 7 covariates
set.seed(123)
n<-1000
X1<-rnorm(n=n,sd=sqrt(0.04))
X2<-X1+rnorm(n=n,sd=sqrt(0.04))
X3<-X1+X2+rnorm(n=n,sd=sqrt(0.04))
m<-sin(X2*5)+X3^3
Z<-m+rnorm(n=n,sd=sqrt(0.04))
X4<-X2+rnorm(n=n,sd=sqrt(0.04))
X5<-Z+rnorm(n=n,sd=sqrt(0.04))
X6<-Z+rnorm(n=n,sd=sqrt(0.04))
X7<-X6+rnorm(n=n,sd=sqrt(0.04))
Y<-qpois(pnorm(Z, mean = m, sd = sqrt(0.04)), lambda=exp(m))
dat<-data.frame(X1, X2, X3, X4, X5, X6, X7,Y)
fml<- Y~s(X1)+s(X2)+s(X3)+s(X4)+s(X5)+s(X6)+s(X7)
mod.all <-cgam(fml,"poisson",dat,pval="chi-square",search="all")
mod.all$model.opt
mod.step <-cgam(fml,"poisson",dat,pval="chi-square",search="stepwise")
mod.step$model.opt
####################################
#causal logistic gam################
n<-1000
set.seed(123)
X1<-rnorm(n,0,1)
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)
cm_all<-cgam(Y ~ s(X1)+s(X2),"binomial",data,pval="bootstrap",search="all")
cm_all$model.opt
set.seed(1)
cm_step<-cgam(Y ~ s(X1)+s(X2),"binomial",data,pval="bootstrap",search="stepwise")
cm_step$model.opt

Causal generalized linear model

Description

This function does a search for a causal submodel within the generalized linear model provided.

Usage

cglm(
  formula,
  family,
  data,
  alpha = 0.05,
  pval = c("bootstrap", "chi-square"),
  B = 100,
  search = c("all", "stepwise"),
  ncores = 1L,
  use_cpp = TRUE,
  direction = c("forward", "backward"),
  ...
)
cglm(
  formula,
  family,
  data,
  alpha = 0.05,
  pval = c("bootstrap", "chi-square"),
  B = 100,
  search = c("all", "stepwise"),
  ncores = 1L,
  use_cpp = TRUE,
  direction = c("forward", "backward"),
  ...
)

Arguments

formula

A formula object.

family

The response family as a character string, "poisson" or "binomial". A family object such as poisson() is not accepted.

data

A data frame containing the variables in the model.

alpha

Significance level for statistical test

pval

If pval="bootstrap", a bootstrap test is conducted to test whether Pearson risk is 1. When family="poisson" a chi-squared test can be conducted by setting pval="chi-square".

B

Number of bootstrap samples when pval="bootstrap". Default is 100.

search

"all" (the default) fits every non-empty submodel of the candidate terms; "stepwise" runs a stepwise search in the direction set by direction.

ncores

use_cpp

Logical; if TRUE (default), use fast C++ implementations for supported families.

direction

...

Further arguments to be passed to the glm function.

Value

A list containing the selected causal submodel and search diagnostics.

References

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

Examples

###################################
#causal Poisson glm#################
n<-1000
set.seed(123)
X1<-rnorm(n,0,1)
Y<-rpois(n,exp(X1))
X2<-log(Y+1)+rnorm(n,0,0.3)
data<-data.frame(X1, X2, Y)
cm_all<-cglm(Y ~ X1+X2,"poisson",data,pval="chi-square",search="all")
cm_all$model.opt
cm_step<-cglm(Y ~ X1+X2,"poisson",data,pval="chi-square",search="stepwise")
cm_step$model.opt

##########################
#causal logistic glm#######
n<-2000
set.seed(123)
X1<-rnorm(n,0,1)
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)
cm_all<-cglm(Y ~ X1+X2,"binomial",data,pval="bootstrap",search="all")
cm_all$model.opt
set.seed(1)
cm_step<-cglm(Y ~ X1+X2,"binomial",data,pval="bootstrap",search="stepwise")
cm_step$model.opt
#bigger simulation with 5 covariates
set.seed(12)
n<-3000
X1<-rnorm(n,0,1)
X2<-rnorm(n,X1,0.5)
X3<-rnorm(n,0,1)
X4<-rnorm(n,X2,.5)
Y<-rbinom(n,1,exp(.8*X2-.9*X3)/(1+exp(.8*X2-.9*X3)))
flip<-rbinom(n,1,0.1)
X5<-(1-flip)*Y+flip*(1-Y)+rnorm(n,0,.3)
dat<-data.frame(X1, X2, X3, X4, X5,Y)
set.seed(1)
mod.all <-cglm(Y~X1+X2+X3+X4+X5,"binomial",dat,pval="bootstrap",search="all")
mod.all$model.opt
set.seed(1)
mod.step <-cglm(Y~X1+X2+X3+X4+X5,"binomial",dat,pval="bootstrap",search="stepwise")
mod.step$model.opt

###################################
#causal Poisson glm#################
n<-1000
set.seed(123)
X1<-rnorm(n,0,1)
Y<-rpois(n,exp(X1))
X2<-log(Y+1)+rnorm(n,0,0.3)
data<-data.frame(X1, X2, Y)
cm_all<-cglm(Y ~ X1+X2,"poisson",data,pval="chi-square",search="all")
cm_all$model.opt
cm_step<-cglm(Y ~ X1+X2,"poisson",data,pval="chi-square",search="stepwise")
cm_step$model.opt

##########################
#causal logistic glm#######
n<-2000
set.seed(123)
X1<-rnorm(n,0,1)
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)
cm_all<-cglm(Y ~ X1+X2,"binomial",data,pval="bootstrap",search="all")
cm_all$model.opt
set.seed(1)
cm_step<-cglm(Y ~ X1+X2,"binomial",data,pval="bootstrap",search="stepwise")
cm_step$model.opt
#bigger simulation with 5 covariates
set.seed(12)
n<-3000
X1<-rnorm(n,0,1)
X2<-rnorm(n,X1,0.5)
X3<-rnorm(n,0,1)
X4<-rnorm(n,X2,.5)
Y<-rbinom(n,1,exp(.8*X2-.9*X3)/(1+exp(.8*X2-.9*X3)))
flip<-rbinom(n,1,0.1)
X5<-(1-flip)*Y+flip*(1-Y)+rnorm(n,0,.3)
dat<-data.frame(X1, X2, X3, X4, X5,Y)
set.seed(1)
mod.all <-cglm(Y~X1+X2+X3+X4+X5,"binomial",dat,pval="bootstrap",search="all")
mod.all$model.opt
set.seed(1)
mod.step <-cglm(Y~X1+X2+X3+X4+X5,"binomial",dat,pval="bootstrap",search="stepwise")
mod.step$model.opt

Evaluate multiple submodels in C++ (for exhaustive search)

Description

Given a full design matrix and a list of column index vectors, fit each submodel and return Pearson risk, p-value, BIC.

Usage

eval_submodels_cpp(X_full, y, family, col_indices, pval_method, B = 100L)
eval_submodels_cpp(X_full, y, family, col_indices, pval_method, B = 100L)

Arguments

X_full

Full design matrix (with intercept)

y

Response vector

family

"poisson" or "binomial"

col_indices

List of integer vectors, each specifying columns of X_full

pval_method

"chi-square" or "bootstrap"

B

Number of bootstrap replicates

Value

List of vectors: pearson_risk, pval, bic

Fit GLM and compute Pearson risk, BIC, and optionally bootstrap p-value

Description

Fit GLM and compute Pearson risk, BIC, and optionally bootstrap p-value

Usage

fast_fit_and_stat(X, y, family)
fast_fit_and_stat(X, y, family)

Arguments

X

Design matrix

y

Response vector

family

"poisson" or "binomial"

Value

List with pearson_stat, pearson_risk, bic

Compute BIC for a GLM fit

Description

Compute BIC for a GLM fit

Usage

fast_glm_bic(y, mu, family, p, n)
fast_glm_bic(y, mu, family, p, n)

Arguments

y

Response vector

mu

Fitted values

family

"poisson" or "binomial"

p

Number of parameters

n

Number of observations

Value

BIC value

Fast GLM fit via IRLS

Description

Fast GLM fit via IRLS

Usage

fast_glm_fit(X, y, family, max_iter = 25L, tol = 1e-08)
fast_glm_fit(X, y, family, max_iter = 25L, tol = 1e-08)

Arguments

X

Design matrix (n x p), including intercept column if needed

y

Response vector (n)

family

"poisson" or "binomial"

max_iter

Maximum IRLS iterations (default 25)

tol

Convergence tolerance (default 1e-8)

Value

List with fitted_values (mu) and converged flag

Compute log-likelihood for Poisson or Binomial GLM

Description

Compute log-likelihood for Poisson or Binomial GLM

Usage

glm_loglik_cpp(y, mu, family)
glm_loglik_cpp(y, mu, family)

Arguments

y

Response vector

mu

Fitted values

family

"poisson" or "binomial"

Value

Log-likelihood value

Compute Pearson chi-square statistic from y and mu

Description

Compute Pearson chi-square statistic from y and mu

Usage

pearson_stat_cpp(y, mu, family)
pearson_stat_cpp(y, mu, family)

Arguments

y

Response vector

mu

Fitted values

family

"poisson" or "binomial"

Value

Pearson chi-square statistic (sum of squared Pearson residuals)

Package 'causalreg'

Help Index

Bootstrap p-value for Pearson risk = 1 test (fully in C++)

Description

Usage

Arguments

Value

Causal generalized additive model

Description

Usage

Arguments

Value

References

Examples

Causal generalized linear model

Description

Usage

Arguments

Value

References

Examples

Evaluate multiple submodels in C++ (for exhaustive search)

Description

Usage

Arguments

Value

Fit GLM and compute Pearson risk, BIC, and optionally bootstrap p-value

Description

Usage

Arguments

Value

Compute BIC for a GLM fit

Description

Usage

Arguments

Value

Fast GLM fit via IRLS

Description

Usage

Arguments

Value

Compute log-likelihood for Poisson or Binomial GLM

Description

Usage

Arguments

Value

Compute Pearson chi-square statistic from y and mu

Description

Usage

Arguments

Value