coxed: An R Package for Computing Duration-Based Quantities from the Cox Proportional Hazards Model

The Cox proportional hazards model is one of the most frequently used estimators in duration (survival) analysis. Because it is estimated using only the observed durations’ rank ordering, typical quantities of interest used to communicate results of the Cox model come from the hazard function (e.g., hazard ratios or percentage changes in the hazard rate). These quantities are substantively vague and difficult for many audiences of research to understand. We introduce a suite of methods in the R package coxed to address these problems. The package allows researchers to calculate duration-based quantities from Cox model results, such as the expected duration (or survival time) given covariate values and marginal changes in duration for a specified change in a covariate. These duration-based quantities often match better with researchers’ substantive interests and are easily understood by most readers. We describe the methods and illustrate use of the package. Introduction The Cox proportional hazards model (Cox, 1972) is frequently used for duration (survival) analysis in a myriad of disciplines including the health sciences, social sciences, operations research, and engineering. For many researchers who employ the Cox model, the chief concept of substantive interest is the duration of an event, such as the survival time of a patient or the duration of a civil war. However, the standard methods of reporting results from the Cox model—which are based in the hazard function—communicate no specific information about duration. As a result, standard interpretations of Cox model results are often substantively vague and difficult for many audiences of research to understand. Here we introduce an R package implementation of Cox proportional hazards model with expected durations, or COX ED (Kropko and Harden, 2020). The COX ED suite of methods available in the coxed package provides a more intuitive approach to communicating results from the Cox model. Specifically, it computes duration-based quantities of interest, such as the expected time until event occurrence according to the estimated model. These quantities have long been available with parametric duration models, but in some instances researchers may not wish to make the distributional assumptions required of those estimators. The COX ED methods allow researchers to stay within the Cox model framework, but communicate results in the language of time. This affords more conceptual precision when conversing with other researchers and makes the results of the analysis more intuitive and accessible for general audiences. The Methodology The goal of COX ED is to generate expected durations for individual observations and/or marginal changes in expected duration given a change in a covariate from the Cox model. Specifically, the methods can compute (1) the expected duration for each observation used to fit the Cox model, given the covariates, (2) the expected duration for a “new” observation with a covariate profile set by the analyst, or (3) the first difference, or change, in expected duration given two new observations. There are two different methods of generating duration-based quantities in the package. The first method employs a generalized additive model (GAM) to map the model’s exponentiated linear predictor values to duration times. The second method calculates expected durations by using nonparametric estimates of the baseline hazard and survivor functions. We present overviews of these methods here. See Kropko and Harden (2020) for additional details, including simulation results comparing the two methods. Importantly, both approaches use coefficient estimates from the Cox model, so researchers must first estimate the model just as they always have. COX ED is a postestimation procedure, not a new estimator. All of the choices required of applied researchers in estimating the Cox model must be made first, at the estimation stage, before proceeding to implement COX ED.1 1Additionally, because it is used after estimation, more extensive modeling features—such as non-linear effects or time-varying covariates—can be incorporated into the use of COX ED. The R Journal Vol. XX/YY, AAAA 20ZZ ISSN 2073-4859 CONTRIBUTED RESEARCH ARTICLE 2 Method 1: GAM The GAM approach to COX ED proceeds according to five steps. As is noted above, the first step is model estimation. Then the method computes expected values of risk for each observation by matrixmultiplying the covariates, X, by the estimated coefficients from the model, β̂, then exponentiating the result. This creates exp(X β̂), or the exponentiated linear predictor (ELP). Then the observations are ranked from smallest to largest according to their values of the ELP. This ranking is interpreted as the expected order of failure; the larger the value of the ELP, the sooner the model expects that observation to fail, relative to the other observations. The next step is to connect the model’s expected risk for each observation (ELP) to duration time (the observed durations). A GAM fits a model to data by using a series of locally-estimated polynomial splines set by the user (Hastie and Tibshirani, 1990). It is a flexible means of allowing for the possibility of nonlinear relationships between variables. COX ED uses a GAM to model the observed durations as a function of the linear predictor ranks generated in the previous step. More specifically, the method utilizes a cubic regression spline to draw a smoothed line summarizing the bivariate relationship between the observed durations and the ranks (for more details, see Wood, 2006, 2011).2 The GAM fit can be used directly to compute expected durations, given the covariates, for each observation in the data. However, for most researchers it is more useful to assess how a change to a particular covariate of interest corresponds to changes in expected duration. In order to examine such marginal changes, it is necessary to create two or more “new” observations corresponding to theoretically-interesting, hypothetical covariate profiles. For example, the analyst might set an indicator variable to 0 and 1 or a continuous variable to a “low” and a “high” value. COX ED allows the covariates in the model to vary naturally over the entire data, then averages over them in the computations.3 For instance, to estimate the effect of an increase in a covariate X1 from 0 to 1 on the expected duration, we use the following steps: (a) Set X1 to 1 for the entire data (all N observations) and calculate the ELP for every observation, then take an average value of those computations (the median is the default). (b) Repeat step (a) while setting X1 equal to 0. (c) Take the values obtained in steps (a) and (b) and append them to the list of ELP values from the original Cox model in which X1 is left as exogenous data. Then compute new rankings of the linear predictor values from this list, which is now length N + 2. (d) Pass the list of rankings from step (c) to the GAM as new data to generate expected values. Note that a new GAM is not estimated at this step. Rather, expected durations are generated for each observation—including the two new ones created in steps (a) and (b)—using the previously estimated GAM. This produces point estimates of the expected durations for those two new observations. (e) Compute the difference between the two estimates obtained in step (d): the expected duration for the data in which X1 is set to 1 and the expected duration for the data in which X1 is set to 0. This quantity is a point estimate for the marginal effect, or first difference, corresponding to the change in X1 from 0 to 1. Finally, to produce estimates of uncertainty, the GAM approach repeats this process via bootstrapping. The method generates bootstrap samples of the data and re-estimates the Cox model coefficients on each bootstrap sample.4 At each iteration, this produces a new vector of actual durations and a new ranking of ELP values, which are then used to fit a new GAM. This process results in a distribution of expected durations for each independent variable profile (e.g., step d) and a distribution of the marginal effect (step e). These distributions can be used to produce standard errors and confidence intervals for the estimates.5 Importantly, by bootstrapping the entire process, this step incorporates the uncertainty from the Cox model estimation and the uncertainty from the GAM. 2The GAM is fit with the uncensored observations only. If the sample contains a large proportion of censored observations, the NPSF method (see below) may be preferable to the GAM method. 3This default option can be changed at the discretion of the analyst. 4Standard bootstrapping at the observation level or bootstrapping at the group level (Cameron et al., 2008) are both available. 5By default, the method computes the standard errors of each quantity as the standard deviation of its bootstrap distribution. The halfwidth of the confidence interval is then computed by multiplying a tunable critical value based on the standard normal distribution by the standard error. The default critical value is 1.96 (i.e., a 95% confidence interval). Fully non-parametric confidence intervals based on quantiles or bias-corrected quantiles of the bootstrap distribution are also available (see below). The R Journal Vol. XX/YY, AAAA 20ZZ ISSN 2073-4859 CONTRIBUTED RESEARCH ARTICLE 3 Method 2: Nonparametric step-function One drawback to the GAM approach is that it uses two statistical models (Cox model and GAM), which yields two sources of estimation uncertainty. An alternative approach comes from the method proposed by Cox and Oakes (1984, 107–109) for estimating the cumulative baseline hazard function. This method is nonparametric and results in a step-function representation of the cumulative baseline hazard; we refer to it as the nonparametric step-function (NPSF) approach. Cox and Oakes (1984, 108) show that the cumulative baseline


Introduction
The Cox proportional hazards model (Cox, 1972) is frequently used for duration (survival) analysis in a myriad of disciplines including the health sciences, social sciences, operations research, and engineering. For many researchers who employ the Cox model, the chief concept of substantive interest is the duration of an event, such as the survival time of a patient or the duration of a civil war. However, the standard methods of reporting results from the Cox model-which are based in the hazard function-communicate no specific information about duration. As a result, standard interpretations of Cox model results are often substantively vague and difficult for many audiences of research to understand.
Here we introduce an R package implementation of Cox proportional hazards model with expected durations, or COX ED (Kropko and Harden, 2020). The COX ED suite of methods available in the coxed package provides a more intuitive approach to communicating results from the Cox model. Specifically, it computes duration-based quantities of interest, such as the expected time until event occurrence according to the estimated model. These quantities have long been available with parametric duration models, but in some instances researchers may not wish to make the distributional assumptions required of those estimators. The COX ED methods allow researchers to stay within the Cox model framework, but communicate results in the language of time. This affords more conceptual precision when conversing with other researchers and makes the results of the analysis more intuitive and accessible for general audiences.

The methodology
The goal of COX ED is to generate expected durations for individual observations and/or marginal changes in expected duration given a change in a covariate from the Cox model. Specifically, the methods can compute (1) the expected duration for each observation used to fit the Cox model, given the covariates, (2) the expected duration for a "new" observation with a covariate profile set by the analyst, or (3) the first difference, or change, in expected duration given two new observations. There are two different methods of generating duration-based quantities in the package. The first method employs a generalized additive model (GAM) to map the model's exponentiated linear predictor values to duration times. The second method calculates expected durations by using nonparametric estimates of the baseline hazard and survivor functions. We present overviews of these methods here. See Kropko and Harden (2020) for additional details, including simulation results comparing the two methods. Importantly, both approaches use coefficient estimates from the Cox model, so researchers must first estimate the model just as they always have. COX ED is a postestimation procedure, not a new estimator. All of the choices required of applied researchers in estimating the Cox model must be made first, at the estimation stage, before proceeding to implement COX ED. 1

Method 1: GAM
The GAM approach to COX ED proceeds according to five steps. As is noted above, the first step is model estimation. Then the method computes expected values of risk for each observation by matrixmultiplying the covariates, X, by the estimated coefficients from the model,β, then exponentiating the result. This creates exp(Xβ), or the exponentiated linear predictor (ELP). Then the observations are ranked from smallest to largest according to their values of the ELP. This ranking is interpreted as the expected order of failure; the larger the value of the ELP, the sooner the model expects that observation to fail, relative to the other observations. The next step is to connect the model's expected risk for each observation (ELP) to duration time (the observed durations). A GAM fits a model to data by using a series of locally-estimated polynomial splines set by the user (Hastie and Tibshirani, 1990). It is a flexible means of allowing for the possibility of nonlinear relationships between variables. COX ED uses a GAM to model the observed durations as a function of the linear predictor ranks generated in the previous step. More specifically, the method utilizes a cubic regression spline to draw a smoothed line summarizing the bivariate relationship between the observed durations and the ranks (for more details, see Wood, 2006Wood, , 2011 The GAM fit can be used directly to compute expected durations, given the covariates, for each observation in the data. However, for most researchers it is more useful to assess how a change to a particular covariate of interest corresponds to changes in expected duration. In order to examine such marginal changes, it is necessary to create two or more "new" observations corresponding to theoretically-interesting, hypothetical covariate profiles. For example, the analyst might set an indicator variable to 0 and 1 or a continuous variable to a "low" and a "high" value. COX ED allows the covariates in the model to vary naturally over the entire data, then averages over them in the computations. 3 For instance, to estimate the effect of an increase in a covariate X 1 from 0 to 1 on the expected duration, we use the following steps: (a) Set X 1 to 1 for the entire data (all N observations) and calculate the ELP for every observation, then take an average value of those computations (the median is the default).
(b) Repeat step (a) while setting X 1 equal to 0.
(c) Take the values obtained in steps (a) and (b) and append them to the list of ELP values from the original Cox model in which X 1 is left as exogenous data. Then compute new rankings of the linear predictor values from this list, which is now length N + 2.
(d) Pass the list of rankings from step (c) to the GAM as new data to generate expected values. Note that a new GAM is not estimated at this step. Rather, expected durations are generated for each observation-including the two new ones created in steps (a) and (b)-using the previously estimated GAM. This produces point estimates of the expected durations for those two new observations.
(e) Compute the difference between the two estimates obtained in step (d): the expected duration for the data in which X 1 is set to 1 and the expected duration for the data in which X 1 is set to 0. This quantity is a point estimate for the marginal effect, or first difference, corresponding to the change in X 1 from 0 to 1.
Finally, to produce estimates of uncertainty, the GAM approach repeats this process via bootstrapping. The method generates bootstrap samples of the data and re-estimates the Cox model coefficients on each bootstrap sample. 4 At each iteration, this produces a new vector of actual durations and a new ranking of ELP values, which are then used to fit a new GAM. This process results in a distribution of expected durations for each independent variable profile (e.g., step d) and a distribution of the marginal effect (step e). These distributions can be used to produce standard errors and confidence intervals for the estimates. 5 Importantly, by bootstrapping the entire process, this step incorporates the uncertainty from the Cox model estimation and the uncertainty from the GAM. 2 The GAM is fit with the uncensored observations only. If the sample contains a large proportion of censored observations, the NPSF method (see below) may be preferable to the GAM method.
3 This default option can be changed at the discretion of the analyst. 4 Standard bootstrapping at the observation level or bootstrapping at the group level (Cameron et al., 2008) are both available. 5 By default, the method computes the standard errors of each quantity as the standard deviation of its bootstrap distribution. The halfwidth of the confidence interval is then computed by multiplying a tunable critical value based on the standard normal distribution by the standard error. The default critical value is 1.96 (i.e., a 95% confidence interval). Fully non-parametric confidence intervals based on quantiles or bias-corrected quantiles of the bootstrap distribution are also available (see below).

Method 2: Nonparametric step-function
One drawback to the GAM approach is that it uses two statistical models (Cox model and GAM), which yields two sources of estimation uncertainty. An alternative approach comes from the method proposed by Cox and Oakes (1984, 107-109) for estimating the cumulative baseline hazard function. This method is nonparametric and results in a step-function representation of the cumulative baseline hazard; we refer to it as the nonparametric step-function (NPSF) approach. Cox and Oakes (1984, 108) show that the cumulative baseline hazard function can be estimated after fitting a Cox model byĤ where τ j represents time points earlier than t, d j is a count of the total number of failures at τ j , (τ j ) is the remaining risk set at τ j , andψ(l) represents the ELP from the Cox model for observations still in the risk set at τ j . The NPSF method uses equation (1) to calculate the cumulative baseline hazard at all time points in the range of observed durations with the following steps.
(a) Tied durations are handled by collapsing the dataset by unique duration. The method calculates d j , the numerator in equation (1), for all time points τ j by summing the indicator for a noncensored failure within each unique duration (d j = 0 only if all observed durations at τ j are right-censored). Additionally, it sums the ELPs for all observations with the same duration, because these observations leave the risk set at the same time.
(b) The NPSF approach calculates a running sum, in reverse, for the collapsed ELPs. That is, at the first time point this sum includes the ELP for observations at every time point. At the second time point, this sum includes the ELP for every observation except for those with the earliest observed duration. At the last time point, this sum is equal to the sum of only the ELPs of observations with the latest observed duration. These sums represent the denominator of equation (1).
(c) For each time point, the method divides the number of failures d j by the sum of ELPs for observations still in the risk set.
(d) Finally, the method calculates the running sum of the ratios we derived in the previous step. This running sum is the non-parametric estimate of the cumulative hazard function.
This procedure yields a stepwise function. Time points with no failures do not contribute to the cumulative hazard, so the function is flat until the next time point with observed failures.
The NPSF approach next obtains expected durations and marginal changes in expected duration by first calculating the baseline survivor function from the cumulative hazard function, usinĝ Each observation's survivor function is related to the baseline survivor function bŷ whereψ(i) is the ELP for observation i. These survivor functions can be used directly to calculate expected durations for each observation. The expected value of a non-negative random variable can be calculated by where F(.) is the cumulative distribution function for X. In the case of a duration variable t i , the expected duration is where T is the largest possible duration and S(t) is the individual's survivor function. The NPSF method approximates this integral with a right Riemann-sum by calculating the survivor functions at every discrete time point from the minimum to the maximum observed durations, and multiplying these values by the length of the interval between time points with observed failures: To calculate a marginal effect, the NPSF approach to COX ED follows the same strategy employed in the GAM approach. It creates two new covariate profiles, setting a variable of interest to two theoretically interesting values. It calculates expected values from each profile, then computes the difference in the two estimates. Finally, the method bootstraps to obtain a standard error and/or confidence intervals for this point estimate.

Implementation in R and empirical example
The methods described above are mostly automated in the package; analysts generally need only a coxph model object from the survival package (Therneau, 2015) or a cph model object from the rms package (Harrell, 2018), and, if covariate effects are desired, the name of the variable of interest and the two values of that variable they wish to input. 6 However, the functions also allow for several changes to default settings, such as the formulation of the GAM in the first approach or the computation of confidence intervals.
We illustrate the main features of the package with an empirical example. Martin and Vanberg (2003) examine the determinants of negotiation time among political parties forming coalition governments in Western Europe. The outcome variable in this analysis is the number of days between the beginning and end of the bargaining period. The covariates include the range of government-the ideological distance between the extreme members of the coalition-the number of parties in the coalition, as well as several others. Their main hypotheses predict negative coefficients on the range of government and number of parties variables. They expect that increases in the ideological distance between the parties and the size of the coalition correspond with decreases in the risk of government formation, or longer negotiation times.
The authors demonstrate support for their hypotheses with a sample of data on bargaining in Western European democracies between 1950 and 1995. They estimate a Cox model, then interpret the covariate effects with quantities based in the hazard rate. As an alternative, we employ COX ED with these data. We use the coxed() function to predict bargaining duration for every case in the data. Then test the first of their hypotheses by computing estimates of bargaining duration at different values of ideological range of government.
Next we use the GAM version of coxed() to examine expected durations and marginal changes in duration. 7 We can calculate standard errors and confidence intervals for any of these quantities with the bootstrap = TRUE option. By default the bootstrapping procedure uses 200 iterations (to set this value to a different number, use the B argument). 8 ed <-coxed(mv.cox, method = "gam", bootstrap = TRUE, B = 30) 6 Future versions of the software may accept Cox models estimated from other packages, such as timereg (Scheike and Zhang, 2011). 7 For an example of this analysis using the NPSF method, see the vignette for the coxed package. 8 Here we use 30 iterations simply to ease the computational burden of compiling this example. For more reliable results, set B to a higher value. There are different methods for calculating a bootstrapped confidence interval. The default method used by coxed() (setting the argument confidence = "studentized") adds and subtracts qnorm(level -(1 -level)/2) times the bootstrapped standard error to the point estimate, where level is the analyst's chosen threshold for evaluating statistical significance. The alternative approach is to take the (1 -level)/2 and level + (1 -level)/2 quantiles of the bootstrapped draws, which can be done by specifying confidence = "empirical". We recommend a higher number of bootstrap iterations for empirical confidence intervals. Additionally, the nonparametric bias corrected and accelerated (BC a ) method can be computed with confidence = "bca", which implements the bias correction and acceleration procedure in DiCiccio and Efron (1996) using code modified from the mediation package (Tingley et al., 2014).  Martin and Vanberg (2003). Entries report coefficients with standard errors in parentheses. These results represent a common approach to presenting Cox model output, but the coefficients themselves are not immediately intuitive. new.coalitions <-data.frame(postel = c(1, 1, 1, 0, 1), prevdef = c(0, 0, 1, 1, 0), cont = c(1, 0, 1, 0, 1), ident = c(1, 2, 2, 3, 3), rgovm = c(.3, .8, 1.1, .2, .35), pgovno = c(2, 3, 3, 2, 4), tpgovno = c(3.2, 0, 5, 0, 2.6), minority = c(0, 0, 1, 0, 0)) To forecast durations for these cases along with standard errors and confidence intervals, we use the coxed() function and place new.coalitions into the newdata argument: forecast The data used by coxed() to map rankings to durations are stored in the gam.data attribute, and can be used to visualize the fit of the GAM, as in Figure 1. Cox model LP rank (smallest to largest)

Duration
We use coxed() to provide an answer to the key question, "how much longer will negotiations take for an ideologically polarized coalition as compared to an ideologically homogeneous one?" Specifically, we call coxed() and specify two new datasets, one in which rgovm = 0 indicating that all political parties in the governing coalition have the same ideological position (i.e., a coalition of one party), and one in which rgovm = 1.24, indicating that the parties have very different ideological positions. 9 We use mutate() from the dplyr package (Wickham et al., 2018) to quickly create new data frames in which rgovm equals 0 or 1.24 for all cases, and set these two data frames as newdata and newdata2 inside coxed( These results demonstrate that a coalition in which the parties have average ideological differences will take 3.6 more days on average (with a median of 2.9 days) to conclude negotiations than a coalition in which all parties have the same position (i.e., a single-party government).
The NPSF method can be used to compute estimates of these same quantities simply by specifying method = "npsf" in the coxed() function. Additionally, the package includes a function called sim.survdata() designed for simple simulations of duration data that do not assume a distributional form for the baseline hazard. This method, which is fully described in Harden and Kropko (2019), can be useful in several applied and computational settings that involve the Cox model.

Conclusions
The Cox model is popular among applied researchers in a wide range of disciplines due to its inherent flexibility. However, this flexibility makes conveying the substantive meaning of results challenging. By using only the rank ordering of the observed duration times, the Cox model limits researchers to interpreting results in the language of hazard and changes in risk. This yields two key problems. First, it is substantively vague because hazard does not have a meaningful scale. This hinders researchers' capacity to determine whether an estimated effect is substantively "large" or "small." Furthermore, hazard-based interpretations require specialized knowledge to understand. This makes the research less accessible to general audiences, who may be able to learn from the work but cannot due to the means by which results are communicated.
The COX ED methods provide a solution to these problems by allowing researchers to compute duration-based quantities from the Cox model. Communicating results in the language of time allows for more substantive precision and is intuitive to a broad audience of readers. We demonstrate above that COX ED is straightforward to implement in R. The coxed package contains functions that allow researchers to use the methods even with minimal knowledge of R. Additionally, the functions are flexible; users can make several changes to many of their features to suit the problem at hand. Finally, the output from the functions provide point estimates, standard errors, and confidence intervals, so researchers can report their results with appropriate measures of uncertainty.
In sum, the coxed package provides a useful alternative for researchers to communicate results from the Cox model. It gives them the benefits of the intuitive quantities available in parametric models while retaining the desirable estimation properties of the Cox model. Thus, the analysis can be guided by appropriate modeling choices, but reported in an intuitive, accessible manner.