lillardhaz (Python)

A Python package for simultaneous-equations hazard/probit models with a Gaussian-copula correlation, in the style of Lillard (1993): a pair of processes — each a binary probit outcome or a continuous-time hazard duration — linked through a single correlation parameter between their underlying error terms, estimated jointly by maximum likelihood.

Four model families are supported via eq1/eq2 in fit_lillardhaz():

pgompertz is a piecewise-linear-in-time log-hazard (piecewise Gompertz) with an arbitrary, user-specified number of segments (nodes1/nodes2). Every combination can also be fit independently with corr=False (rho fixed at 0), which — as derived in the manual — reduces algebraically to two separate univariate fits, giving a natural nested baseline for a likelihood-ratio test of correlation.

Companion packages implementing the identical models: a Stata command and an R package. See docs/manual.html for the full model derivation.

Author

Nobutaka Fukuda, Tohoku University — nobutaka.fukuda@tohoku.ac.jp

Installation

pip install lillardhaz

Usage

import numpy as np
from lillardhaz import fit_lillardhaz

# Probit & log-normal hazard, correlated
fit = fit_lillardhaz(
    "probit", "lognormal", data,
    y1="y1", y2="time2", d2="event2",
    x1=["x1"], x2=["x2"],
)
print(fit)
fit.rho, fit.rho_se

# Piecewise Gompertz & piecewise Gompertz, correlated, 2 interior nodes each side
fit = fit_lillardhaz(
    "pgompertz", "pgompertz", data,
    y1="time1", d1="event1", y2="time2", d2="event2",
    x1=["x1"], x2=["x2"], nodes1=[5], nodes2=[4],
)

# Same, but independent (nested test of rho = 0)
fit0 = fit_lillardhaz(
    "pgompertz", "pgompertz", data,
    y1="time1", d1="event1", y2="time2", d2="event2",
    x1=["x1"], x2=["x2"], nodes1=[5], nodes2=[4], corr=False,
)
lr_stat = 2 * (fit.loglik - fit0.loglik)   # LR test of rho = 0, ~chisq(1)

# Predictions
fit.predict(data, kind="surv2")   # fitted eq2 survival (default)
fit.predict(data, kind="dens2")   # fitted eq2 density
fit.predict(data, kind="xb1")     # eq1 linear index

data is any mapping of column name to array-like (a plain dict of numpy arrays, or a pandas DataFrame).

Verification

Every eq1×eq2 combination, correlated and independent, was verified by simulating 20,000–25,000 observations under known parameters and confirming fit_lillardhaz() recovers them; see tests/test_models.py. See the manual's verification section for a summary table.

License

MIT — see LICENSE.

References

Lillard, L. A. 1993. Simultaneous equations for hazards: Marriage duration and fertility timing. Journal of Econometrics 56(1–2): 189–217.

Waite, L. J., & Lillard, L. A. 1991. Children and marital disruption. American Journal of Sociology 96(4): 930–953.