pchreg: Piecewise Constant Hazards Models

Description Usage Arguments Details Value Author(s) References See Also Examples


This function estimates piecewise exponential models on right-censored, left-truncated data. The effect of covariates, and not just the baseline hazard, varies across intervals. Moreover, a special handling of zero-risk regions is implemented. Differently from the phreg function available in the eha package, this function is mainly intended to be used as a nonparametric maximum likelihood estimator.





an object of class “formula”: a symbolic description of the regression model. The response must be a Surv object as returned by Surv (see ‘Details’).


either a numeric vector of two or more unique cut points or a single number (greater than or equal to 1) giving the number of intervals into which the time variable is to be cut. If missing, the number of intervals is set to max(5, min(50, ceiling(n1/q/5))), where n1 is the number of events, and q is the number of predictors.


an optional data frame containing the variables in the model.


an optional vector of weights to be used in the fitting process.


either NULL, or an object created with splinex (see ‘Details’).


The left side of formula must be of the form Surv(time, event) if the data are right-censored, and Surv(time0, time, event) if the data are right-censored and left-truncated (time0 < time). Using Surv(time) is also allowed and indicates that the data are neither censored nor truncated. Note that the response variable (and thus the breaks) can be negative.

To fit the model, the time interval is first divided in sub-intervals as defined by breaks. When the location of breaks is not specified, the empirical quantiles of time[event == 1] are used as cut points. If there is a probability mass, this may result in two or more breaks being equal: in this case, an interval that only includes the mass point is created automatically.

A different costant hazard (exponential) model is then fitted in each sub-interval, using Poisson regression to model the log-hazard as a linear function of covariates. Within each interval, the risk of the event may be zero at some covariate values. For each covariate x, the algorithm will try to identify a threshold c such that all events (in any given interval) occur when x < c (x > c). A zero risk will be automatically fitted above (below) the threshold, using an offset of -100 on the log-hazard.

This type of model can be utilized to obtain a nonparametric maximum likelihood estimator of a conditional distribution, achieving the flexibility of nonparametric estimators while keeping the model parametric in practice. Users unfamiliar with this approach are recommended reading Geman and Hwang (1982) for an overview, and the paper by Ackerberg, Chen and Hahn (2012) describing how this approach can be applied to simplify inference in two-step semiparametric models.

The number of parameters is equal to the number of intervals, multiplied by the number of covariates. Both quantities are usually supposed to increase with the sample size. The special function splinex is a handy tool that facilitates implementing the linear predictor. When splinex is not NULL, each column of the original design matrix (as defined by formula) is automatically replaced with the corresponding spline basis. See the documentation of splinex for details.


An object of class “pch”, which is a list with the following items:


the matched call.


a matrix of regression coefficients. Rows correspond to covariates, while columns correspond to different time intervals.


the used cut points, with attributes 'h' indicating the length of each interval, and 'k' denoting the number of intervals.


the estimated asymptotic covariance matrix.


the value of the maximized log-likelihood, with attribute “df” indicating the number of free model parameters.


the fitted hazard values in each interval.


the fitted cumulative hazard values at the end of each interval.


the model frame used.


the model matrix.


a code indicating the convergence status. It takes value 0 if the algorithm has converged successfully; 1 if convergence has not been achieved; and 2 if, although convergence has been achieved, more than 1% of observations have an associated survival numerically equal to zero, indicating that the solution is not well-behaved or the model is misspecified.

The accessor functions summary, coef, predict, nobs, logLik, AIC, BIC can be used to extract information from the fitted model. This function is mainly intended for prediction and simulation: see predict.pch.


Paolo Frumento <paolo.frumento@ki.se>


Ackerberg, D., Chen, X., and Hahn, J. (2012). A Practical Asymptotic Variance Estimator for Two-Step Semiparametric Estimators. The Review of Economics and Statistics, 94(2), 481-498.

Friedman, M. (1982). Piecewise Exponential Models for Survival Data with Covariates. The Annals of Statistics, 10(1), pp. 101-113.

Geman, S., and Hwang, C.R. (1982). Nonparametric Maximum Likelihood Estimation by the Method of Sieves. The Annals of Statistics,10(2), 401-414.

See Also

predict.pch, splinex


  # using simulated data
  n <- 1000
  x <- runif(n)
  time <- rnorm(n, 1 + x, 1 + x)
  cens <- rnorm(n,2,2)
  y <- pmin(time,cens) # censored variable
  d <- (time <= cens) # indicator of the event

  model <- pchreg(Surv(y,d) ~ x, breaks = 20)

  # see the documentation of predict.pch

Search within the pch package
Search all R packages, documentation and source code

Questions? Problems? Suggestions? or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.