A fast and flexible machine learning library for nonparametric high-dimensional regression and classification with guarantees.
Code written for hapc ≤ 2.6.0 that uses family = "binomial" must be checked. R and Python behave identically.
2.7.0 was tagged but never reached PyPI (the wheel build failed). 2.7.1 is the first published release with these changes. It also replaces the scikit-learn dependency by scipy: the package no longer uses scikit-learn.
| What | ≤ 2.6.0 | 2.7.0 |
|---|---|---|
Binomial predictions from hapc() |
log-odds for norm = "1" and "2", probabilities for "sv" |
always probabilities; the linear predictor is in log_odds (cv.hapc / cv_hapc already returned probabilities and still do) |
probabilities, predicted_classes fields |
present in some paths | removed: use predictions (and predictions > 0.5 if classes are needed) |
| Intercept of binomial fits | fitted with the intercept fixed at 0, then calibrated afterwards (biased for unbalanced outcomes) | always an unpenalised intercept estimated jointly, in every norm |
intercept argument |
intercept = FALSE allowed (R, norm = "1") |
default NULL/None; intercept = FALSE is an error for binomial / logit-hazard; for gaussian it must agree with center (center = TRUE ⇔ intercept) |
norm = "2" logistic ridge λ |
effective penalty λ‖β‖² | glmnet's scale, -(1/n)·loglik + (λ/2)‖β‖² (same as glmnet(alpha = 0)): λ values from 2.6.0 are not comparable |
norm = "1" binomial basis |
built from the full HAL design (memory O(n²·Σ C(p,k))) and rebuilt for every fold × λ | built from the n×n kernel once per fold; λ along a warm-started path |
Python norm = "1" solver |
sklearn liblinear |
glmnet's algorithm (IRLS + coordinate descent, no sklearn); max_iter now counts IRLS steps (default 100) |
Grid edges (norm = "1"; also "sv"/"2" in R) |
— | extend_grid = k extends the λ grid at either edge; results report lambda_position (interior, lower_edge, upper_edge, null_model), n_extend, lambda_max |
| Kernel / cross-kernel | single-threaded | multi-threaded, bit-identical results; set the thread count with the env var HAPC_NUM_THREADS |
.Call("logistic_call") (R, unexported) |
present | removed (dead code) |
Migration: replace fit$probabilities / fit.probabilities by predictions, replace predictions used as log-odds by log_odds, and drop any intercept = FALSE.
The binomial norm = "sv" and norm = "2" paths are still pending a full review: they still materialise the HAL design and do not scale to n ≈ 2000 with p ≈ 30. Prefer norm = "1".
- Python API (rendered from docstrings): https://hapc.readthedocs.io —
configured via
.readthedocs.yamlanddocs/(Sphinx + autodoc). Build locally withpip install -e ".[docs]" && sphinx-build -b html docs docs/_build/html. - R API (rendered from roxygen): a pkgdown
site built by
.github/workflows/pkgdown.yaml(config in_pkgdown.yml). Build locally withRscript -e 'pkgdown::build_site()'.
- Python 3.8+
- C++ compiler (g++, clang, or MSVC)
- CMake 3.15+
- Eigen3
pip install hapcPrebuilt wheels are published for Linux (manylinux2014, x86_64), macOS (Intel + Apple Silicon) and Windows, for CPython 3.8–3.12. No compiler, CMake or Eigen is needed when a wheel is available.
The Linux wheels use the manylinux2014 baseline (glibc 2.17), so
pip install hapc works out of the box on HPC login/compute nodes —
no conda toolchain, devtoolset, or sysroot setup required:
pip install hapcIf you must build from the source distribution (niche architecture, very
old Python, or an air-gapped node), provide a C++17 compiler and either
let CMake fetch Eigen automatically (needs network) or install Eigen and
let find_package(Eigen3) find it:
# with conda compilers (recommended on HPC)
conda install -c conda-forge cxx-compiler cmake eigen
pip install hapc --no-binary hapcpip install git+https://github.com/meixide/hapc.gitOr with editable install for development:
git clone https://github.com/meixide/hapc.git
cd hapc
pip install -e .If installation fails, you may need to install build dependencies:
macOS:
brew install cmake eigenUbuntu/Debian:
sudo apt-get install cmake libeigen3-dev build-essentialWindows:
pip install cmake
# Install Visual Studio Build Tools or use conda
conda install -c conda-forge eigenimport numpy as np
from hapc.single import single_pcghal
from hapc.cv import pcghal_cv
# Generate sample data
X = np.random.randn(100, 5)
Y = X[:, 0] + 0.5 * X[:, 1] + np.random.randn(100) * 0.1
# Single fit with fixed lambda
result = single_pcghal(X, Y, maxdeg=2, npc=5, single_lambda=0.01)
print(f"Risk: {result.optimizer_output.risk:.6f}")
# Cross-validation to select lambda
lambdas = np.logspace(-4, 0, 10)
cv_result = pcghal_cv(X, Y, maxdeg=2, npc=5, lambdas=lambdas, nfolds=5)
print(f"Best lambda: {cv_result.best_lambda:.6f}")
# Make predictions
X_test = np.random.randn(20, 5)
result = single_pcghal(X, Y, maxdeg=2, npc=5, single_lambda=0.01, predict=X_test)
print(f"Predictions: {result.predictions}")from hapc.single import single_pcghal
result = single_pcghal(
X, Y,
maxdeg=2, # Maximum degree of interactions
npc=10, # Number of principal components
single_lambda=0.01,
predict=X_test # Optional: test data for predictions
)from hapc.single import single_pcghal
result = single_pcghal(
X, Y_binary,
maxdeg=2,
npc=10,
single_lambda=0.01,
predict=X_test
)from hapc.cv import pcghal_cv
cv_result = pcghal_cv(
X, Y,
maxdeg=2,
npc=10,
lambdas=np.logspace(-4, 0, 20),
nfolds=5
)
print(cv_result.best_lambda)Estimate the ATE E[Y(1)] − E[Y(0)] with HAPC nuisance models and a
doubly-robust (AIPW) efficient influence function. ate_hapc returns a point
estimate and a (1 − alpha) Wald confidence interval.
from hapc import ate_hapc
# W: covariates (n, p); A: binary treatment in {0,1} or {-1,+1}; Y: outcome
res = ate_hapc(W, Y, A, alpha=0.05, method="undersmooth")
print(res.estimate, res.lower, res.upper)Two bias-control strategies are available through method:
method="undersmooth"(default) — single-sample estimator. The outcome model is undersmoothed (λ pushed below the CV-optimal value) until the empirical influence function is withinσ / (√n · log n). This requires the full PC basis (npcs = n, the default) and a λ grid that reaches small λ (defaultslog_lambda_out_min = -10); otherwise the gate never reaches the low-bias regime andate_hapcemits a warning. Passreport_undersmoothing=Trueto print the|mean(EIF)|-vs-λ path.method="crossfit"— DML-style K-fold cross-fitting (cf_folds, default 5, stratified by treatment). Both nuisances are fit on the training folds and the influence function is evaluated out-of-fold, giving honest point estimates and coverage without undersmoothing. Recommended under good overlap.
Fit a discrete-time logistic hazard model with HAPC. You supply only the
observed right-censored data — baseline covariates X, the observed time
T = min(T_event, C), and the event indicator Delta = 1(T_event <= C) — and
the wrapper performs the person-period expansion (one row per
subject-per-interval-at-risk, hazard label = 1 at the event interval), prepends
the visit time as the first HAL covariate, and cross-validates the binomial fit.
Model. The discrete hazard is the conditional event probability in interval
t given survival up to t, modelled on the logit scale by a HAPC fit f of
the augmented covariate (t, x):
lambda(t | x) = P(T_event = t | T_event >= t, X = x)
logit lambda(t | x) = f(t, x)
Person-period likelihood. Under independent right-censoring the observed-data likelihood factorises over the at-risk intervals,
prod_i prod_{t <= T_i} lambda(t|x_i)^Y_it * (1 - lambda(t|x_i))^(1 - Y_it),
with Y_it = 1(T_event_i = t),
which is exactly the Bernoulli (logistic) likelihood of the expanded
person-period table — so a binomial HAPC fit of Y_it on (t, x_i) estimates
the discrete hazard (Cox 1972; Brown 1975; Allison 1982).
Survival. The conditional survival function follows by the product-limit
relation S(t | x) = prod_{s <= t} (1 - lambda(s | x)), returned for new
subjects when predict= is supplied.
from hapc import hazard_hapc
import numpy as np
# X: baseline covariates (n, p); T: observed times; Delta: 0/1 event indicator
fit = hazard_hapc(X, T, Delta, norm="1", max_degree=2, time_grid=np.arange(1, 7))
fit.hazard # estimated hazard per person-period row (CV predictions)
fit.best_lambda, fit.interior # CV-selected lambda; is it interior to the grid?
# survival curves S(t|x) for new subjects
fit = hazard_hapc(X, T, Delta, norm="1", predict=X_new)
fit.predict_survival # (m, K) survival probabilities over the gridlibrary(hapc)
# equivalent to cv.hapc(X, T, family = "logit-hazard", Delta = Delta, norm = "1")
fit <- hazard.hapc(X, T, Delta, norm = "1", max_degree = 2, time_grid = 1:6)
fit$hazard; fit$best_lambda; fit$interiornorm must be "1" (logistic LASSO) or "2" (logistic ridge); norm = "sv"
is not implemented for this family and is flagged.
Returns (Python HazardResult / R hapc_hazard):
hazard— cross-validated discrete hazard for each person-period rowlambdas,risk,best_lambda— CV grid, mean logistic deviance, selected λinterior— whetherbest_lambdais strictly inside the grid (sanity check)time_grid,ids/id,Y— the discrete grid and person-period bookkeepingpredict_hazard,predict_survival— hazard surface and survival curves for new subjects (only whenpredict=is given)cv— the underlying cross-validation result
Worked end-to-end examples (five hazard data-generating processes, with
true-vs-estimated hazard scatters and CV risk-vs-λ curves verifying an interior
optimum) are in
examples/hazard_logit_hazard_examples.R
and
examples/hazard_logit_hazard_examples.py.
References. Cox (1972, JRSS B); Brown (1975, Biometrics); Allison (1982, Sociological Methodology); Singer & Willett (2003, Applied Longitudinal Data Analysis); Benkeser & van der Laan (2016, IEEE DSAA).
Fit PC-GHAL with a single lambda value.
Parameters:
X(ndarray, shape (n, p)): Input featuresY(ndarray, shape (n,)): Response variablemaxdeg(int): Maximum degree of interactionsnpc(int): Number of principal componentssingle_lambda(float): Regularization parametermax_iter(int, default=100): Maximum iterationstol(float, default=1e-6): Convergence toleranceverbose(bool, default=False): Print progresspredict(ndarray, optional): Test data for predictionscenter(bool, default=True): Center the design matrix
Returns:
result.optimizer_output.alpha: Coefficientsresult.optimizer_output.risk: Final riskresult.optimizer_output.iter: Iterations until convergenceresult.predictions: Predictions on test data (if provided)
Cross-validation to select lambda.
Parameters:
lambdas(ndarray): Grid of lambda values to testnfolds(int, default=5): Number of CV folds- ...other parameters same as
single_pcghal
Returns:
cv_result.best_lambda: Optimal lambdacv_result.mses: CV errors for each lambdacv_result.best_model: Fitted model with best lambdacv_result.predictions: Predictions on test data (if provided)
Contributions welcome! The C++ core is shared between R and Python packages.
git clone https://github.com/meixide/hapc.git
cd hapc
pip install -e .
pytestMIT License - see LICENSE file
hapc.ate.tmle(W, Y, A) fits the initial outcome and propensity models with
HAPC and targets their cross-fitted predictions. It returns the substitution
ATE, EIF standard error, Wald interval and convergence/clipping diagnostics.
Binary and continuous complete-data outcomes are supported; no SuperLearner
dependency. See implementation and first simulations.
fit <- hapc.ate.tmle(W, Y, A, norm = "1", cf_folds = 5)
c(estimate = fit$estimate, lower = fit$lower, upper = fit$upper)