Cox Regression

Overview Cox regression, which implements the proportional hazards model or duration model, is designed for analysis of time until an event or time between events. One or more predictor variables, called covariates, are used to predict a status (event) variable. The classic univariate example is time from diagnosis with a terminal illness until the event of death (hence survival analysis). Cox regression is also used for policy adoption/diffusion studies (see Jones & Branton, 2005). The central statistical output is the hazard ratio. Note that unlike parametric models discussed in the section on event history models (EHA), Cox regression is semi-parametric and does not require the researcher to specify a baseline hazard rate or estimate absolute risk. For this reason, Cox regression may be preferred over parametric EHA models when there is no clear theoretical reason for positing a particular baseline hazard ratio. In most cases there is no such strong, clear reason and the more stringent data assumptions of parametric EHA models are not justified, making Cox models the better choice. Stata is the preferred software package for Cox regression and survival analysis. In addition to Stata, Limdep is another statistical package with extensive support for event history models, including Cox models. In Stata, declare the data with the stset command, then execute Cox regression with the stcox command. For ordinary Cox regression in SPSS, select Analyze, Survival, Cox Regression; enter the Time variable; enter the Covariates variable(s); enter the Status variable (the event variable) and Define Event to specify the value of the event occurring (ex., death = 1); in Options you may wish to check that you want 'Display baseline function' to get the time-only effect to compared to the covariate effect. Cox regression is a particular model within the broader category of event history analysis. For important related treatment, see the separate discussion of event history methods. See also the separate discussion of Kaplan-Meier survival analysis, a procedure for estimating survival and hazard functions but not covariate effects. See also the life tables procedure, used for descriptive, actuarial studies of duration where time is the only salient variable and censored and uncensored cases do not differ.

Contents Key concepts and terms Model fit Models Stepwise Cox regression Stratified Cox regression Cox regression in Stata Cox regression in SPSS Statistics Assumptions SPSS output Frequently asked questions Bibliography

Key Terms and Concepts

• Variables
  • Status variable. Also called the event or censoring variable, the status variable is the dependent in Cox regression. The classic example is the binary variable death in medical studies, with death equal to 0 for surviving or 1 for dead. However, the researcher may assign a range of values or a list of values to the "event" condition, which does not have to be "1". The status variable is analyzed in relation to a time variable (see below) with hazard or survival rates being the central output of Cox regression. In a study of policy adoption as the status event, logistic regression would focus on analyzing the variance in (the logit of) adoption at the time of data collection. In contrast, Cox regression focuses on analyzing the probability of adoption in any time period. Because apart from historical datasets one does not know the final status for all observations or the time to reach final status, data contain censored and uncensored cases, which is compatible with assumptions of Cox regression but not with those of logistic regression.

  • Time variable. The time variable measures duration to the event defined by the status variable. Normally the "time" variable is a simple counter of units of time since the start of the series. If time is measured by a counter variable in units of time, the Cox model assumes that the hazard rate increases linearly with time, conditional on covariates in the model. It is possible, however, for the time variable to be a logarithmic or some other function of the counter. The significance or non-significance of covariates in the model may vary by the type of time variable used. Smaller time units provide more time intervals, which in turn increases the power of Cox models (less chance of a Type II error: thinking there is no relation of the modeled variables to the dependent when in fact there is).
    • Analysis time is a time variable where t = 0 is the time of onset of risk. Onset of risk means when failure (or the "event") first becomes possible. Analysis time is labeled "t", whereas "time" is used when 0 has other meanings, such as start of measurement. The "origin" is the "time" when analysis time, t = 0. Hence t = time - origin. It is possible for time = t, if the event (ex., adoption) is possible immediately from the start of measurement but no cases have yet incurred the event (none have adopted). Stata assumes analysis time and by default t = time. However, if the time variable in the dataset does not reflect 0 being the onset of risk, but instead there is a different origin, the Stata function origin() can adjust the time variable to be an analysis time variable. Also, the scale() function can adjust the time interval to be what is convenient for analysis (ex., convert days to years). This is done in Stata using the stset command, discussed further below in the section on "Survival time data (st data)".

  • Covariates are the predictor/independent variables in a Cox model.Covariates can be categorical (ex. gender, race or continuous (ex. income, age). Covariates also may be either time-fixed or time-dependent, a difference which affects how the covariate is modeled in Cox procedures. For instance, "contiguity," (ex., coded 0 or 1 to indicate if a state was contiguous to a given state or not) would be time-fixed. "Median family income," which changes by year, would be time-varying (time dependent).
    • Lagging time-dependent covariates. It is recommended (ex., Box-Steffensmeier & Jones, 2004: 111) that time-dependent covariates be entered in lagged form. This is to avoid simultaneity of cause and effect. The covariate is seen as a cause of the event, but if it changes value literally at the same time as the event occurs, the cause-effect logic is lost even though the new value of the covariate is incorporated in the hazard ratio. Also, when measurement of the time of the event is imprecise, lagging helps assure that changes in the time-dependent covariate precede events.

    • Centering. If the researcher will be analyzing baseline hazard rates, the covariate data should be centered (subtract the mean) so that they have a natural zero point. Otherwise the baseline hazard rates, which are the time-only rates when all covariates are zero, are estimated for points which do not exist in the data set, resulting in misleading baseline hazard functions (see Box-Steffensmeier & Jones, 2004: 65).

  • Categorical covariates are explanatory variables which may be used in Cox regression. SPSS will convert them automatically into a set of dummy variables, omitting one category as is usual (by default, the last category, though SPSS allows manually specifying the first instead). Each dummy variable will have its own regression coefficient. It is not necessary to specify as categorical dichotomous covariates that are already indicator coded (0,1) unless the researcher wishes to specify groups for purposes of plots. The interpretation of this coefficient depends on the type of coding scheme:
    1. Indicator, aka "dummy coding," is the default: the regression coefficient compares the effect of the dummy with the reference category (the omitted category of the categorical covariate, usually the last category - SPSS allows the user to specify First or Last as the reference category).
    2. Deviation: The effect of each category other than the reference category is compared to the average effect of all categories.
    3. Repeated: The effect of each category is compared with the next category, except for the last category.
    4. Difference: Each category other than the first is compared to the average effect of all previous categories.
    5. Helmert: Each category other than the last is compared to the average effect of all subsequent categories.
    6. Polynomial: Categories are treated as equally spaced and the covariate is transformed into linear, quadratic, and cubic components, etc.

    The "Categorical Variable Codings" table documents the actual codes applied and is helpful when there is a need to recall what the omitted reference category is.

  • Data setup and example.Data setup is discussed below. In the figure below, the example as implemented in SPSS is predicting time to ratification of the U. S. Constitution. The status variable is "Status" and equals 1 for all states since there are no censored cases (all thirteen states eventually ratified the Constitution). The time variable is "Days," measuring in days how long it took from adoption of the Constitution until the given state ratified it. Other variables are either categorical factors (ex., size of state: small, medium, or large) or continuous covariates (ex., percent by which the vote for ratification passed) or dichotomouse covariates (ex., whether or not the state was a center for Bill of Rights pressure).
    

• Observations
  • Right-censored observations. A case is right-censored when time of the failure/event is known only to have occurred after time t.That is, a right-censored case is one for which the censoring event (the one indicated by the status variable) had not yet occurred by the end of the measurement period. Unless otherwise noted, a "censored observation" is right-censored.
    

  • Left-censored observations. A case is left-censored when the time of the failure/event is known only to have occurred before time t.

  • Truncated observations, also called left-truncated cases, are those which are unmeasured for all time periods up to a given time period, then are measured, usually because they fail (the event occurs) for them in the given time period. (Note right-truncated cases are the same as right-censored observations).

• Functions
  • Survival functions. The cumulative survival function is the percentage of cases which survive to a given point of time (ex., to when the dataset was collected). The baseline survival function is the percentage which would survive based on time alone. The covariate survival function is the percentage that would survive given the covariate(s). The standardized regression coefficient(s) of the covariate(s) is/are a measure of the relative importance of the covariate(s) to survival, controlling for time. Because hazard ratios lend themselves to more intuitive discussion, however, reports of results in Cox regression usually focus on hazard functions.
    

    The example above graphs the cumulative expected probabilities of a state taking the number of days reflected on the X axis prior to voting for ratification of the Constitution, for a hypothetical state which is at the mean of the predictor variables, which are VotePct (percent favoring ratification, reflecting tightness of the vote) and Size (small states and medium states vs. the reference category of large states). Note significance testing may show a covariate is not significant, however.

  • Hazards. The "hazard" is the event of interest occurring. In medical studies the hazard might be death. In industrial studies the hazard might be engine breakdown. However, the hazard may have a positive meaning, as in diffusion studies of time-to-adoption, where the hazard is the adoption of the innovation.

  • Hazard rates and hazard ratios. A hazard rate at a given time is the probability of the given event (dependent = 1) occurring in that time period, given survival through all prior time intervals. A hazard ratio, also called the hazard function, is the estimate of the ratio of the hazard rate in one group (ex., the treatment group) to the hazard rate in another group (ex., the placebo group), for a coded covariate (ex., placebo=0, treatment=1). For a continuous covariate, the hazard ratio is the ratio of the hazard rate given a one unit increase in the covariate to the hazard rate without such an increase. Thus, for example, where units are Democratic and Republican states, the hazard rate for Democratic states might equal the hazard ratio times the hazard rate for Republican states, assuming proportional hazards.

    A key assumption of the Cox model is proportional hazards: the ratio will remain constant over time. The Cox model says nothing about the absolute shape of the curve formed by two hazard rates over time, only that their ratio will be constant.Note that proportional hazards means that hazards are proportional over time, not that they are the same over time. The slopes of the proportional hazard rates for two groups may be downward, for example, indicating decreased hazard over time. Note also that hazard rates are not hazard ratios, and their respective interpretations differ (this is a confusion in a portion of the extant literature using Cox regression).

    • Partial likelihood methods and why Cox models are semi-parametric. Cox models do not assume any particular distribution for the shape of the hazard function, focusing instead on predicting the hazard ratio. Parametric models such as exponential or log-linear event history analysis models, in contrast, do specify the assumed shape of the hazard function. Cox models do not use maximum likelihood estimation but rather a maximum partial likelihood method that requires only the order of failure times, not the intervals between failure times, be known when estimating the hazard ratio. Actual survival times are not used in partial likelihood estimation of the hazard function.
      • Handling tied failure times. Ideally, partial likelihood methods would have no tied data, but rather a simple ordering of failure times. To handle the real-world fact that tied failure times exist, partial likelihood algorithms have been adapted to handle ties. The method of handling ties can be set by the researcher in SPSS and other software. The default method (in SPSS, Stata, and SAS) is the Breslow method. However, in addition to the Breslow method, three other methods are available: the Efron method, the exact marginal-likelihood method, and the exact partial-likelihood method. The Breslow method is adequate when there are few ties. The Efron method is considered more accurate than Breslow when ties are few. When ties are numerous, one of the exact methods may be selected. The choice of tie-breaking method rarely affects substantive results. In Stata, the general Cox regression command syntax is: stcox [varlist] [if] [in] [, options]. The options for the four available methods of handling ties are breslow, efron, exactm, and exactp. For further description of the Breslow and other algorithms, see Box-Steffensmeier & Jones (2004: 54-58).

    • Baseline hazard ratio. The hazard ratio may be partitioned into the baseline hazard ratio (depending on time alone) and the covariate hazard ratio (depending on the covariate(s), controlling for time). The difference between the baseline model and the model with covariates shows the effect of the covariates in the model. Warning: Box-Steffensmeier & Jones (2004: 89) note, "because the Cox estimate of the baseline hazard is so closely tied to the observed data, it is difficult to generalize these estimates to other settings." That is, Cox estimates of the baseline hazard ratio may be considered to be overfitted. This is discussed further below. Since researchers using Cox models are primarily focused on hazard ratios of the covariates, not on the baseline hazard ratio, this may matter little. Where the researcher regards time dependency as important, not just the effect of independent variables (covariates), parametric event history analysis models may still be preferred,
      

      As illustrated above, the baseline cumulative hazard for the intercept only model and the cumulative hazard at the mean of covariates in the full model is shown in the "Survival Table" of SPSS output. This is discussed further below, in the section on statistical output.

    • Hazard ratio with covariates.. The baseline hazard ratio represents the effect of the time variable alone, when all covariate(s) = 0. The hazard ratio tells the odds of an event occurring faster or slower given some covariate(s), but it does not tell how much faster or slower (though it is not uncommon in the literature to find the hazard ratio misinterpreted this way). Interpreting hazard ratios is discussed below, in the section on statistics.

      Example. Hazard ratios below 1.0 indicate that the more the covariate, the less the hazard. Hazard ratios above 1.0 indicate that the greater the covariate, the greater the hazard. Thus, in a model of electric generator life given type or ball bearings and given electric load, if "bearings" = 0 for old style and bearings=1 for new style, and the hazard ratio for bearings is .06, this means that going from old to new style bearings reduced the risk of the generator failing, controlling for load. The hazard ratio of .06 is the proportional change in hazard when the variable bearings increases by 1 unit (i.e., goes from 0 old style to 1 new style). If, however, the low and high confidence intervals on the bearings hazard ratio included 1.0, we could not be sure at the 95% confidence level that bearings really made a difference. For the same model as illustrated above for the covariate survival function, the covariate hazard function looks like this:

      
      • For time-invariant continuous covariates. The hazard ratio for non-time-varying covariates is the amount of change in the hazard of the event occurring for each unit change in the covariate. For instance, a hazard ratio of 1.12 means that there is a 12% increase in the rate of the event occurring for a 1 unit increase in the covariate, controlling for other covariates in the model. A hazard ratio of 1.05 means that for a 1 unit increase in the covariate, there is a 5% increase in the hazard rate of the event variable being studied. A hazard ratio of 1.1 for the covariate age would mean that a one-year increase in age would be associated with a .1 (10%) increase in the hazard rate. Ten years increase in age would correspond to 1.1¹⁰ = 2.59 = increasing the hazard rate by a factor of 2.59 = 159% increase.

      • For time-varying continuous covariates, the hazard ratio is the amount the rate of the event occurring changes for a unit change in the time-dependent function of the covariate.

      • For binary covariates. Where the covariate is a dichotomous grouping variable, such as gender, the hazard ratio for gender=1 is in comparison with the gender=0 group (ex., if the ratio is more than 1.0, the gender=1 group is more likely to incur the event). For instance, given the event not healed=0/healed=1, the covariate placebo=0 /treatment=1, and a computed hazard ratio of 3.0, we can say that a person in the treatment group who has survived to a given time has three times the chances (meaning odds, not probability) as a person in the placebo group of being healed in the next time increment. Since the hazard ratio is a constant, we can also say that by an odds of 3:1 a person in the treatment group is more likely to reach the healed state than a person in a placebo group. This is also the same as saying there is a 3/4 = 75% chance the person in the treatment group will be healed first. We can not say the treated person will heal three times as fast, nor that healing time is cut in a third, nor that three times as many treated persons will be healed by a given time (see Spruance et al., 2004).

        As a second example, for the event "governor reelected = 0, not reelected = 1," for the covariate "Republican state = 0, Democratic state = 1," a hazard ratio of 1.5 would mean that a governor in a Democratic state who had been in office to time t has an odds of 1.5:1 (or 3:2) of not being reelected at time t +1, compared to a governor in a Republican state. This is equivalent to saying that there is a 60% (3/5) chance that the duration in office until the event of non-reelection will occur sooner for a governor in a Democratic state compared to one in a Republican state.

    • Covariate hazard coefficients When the unexponentiated regression coefficient for the covariate hazard function is greater than 0 for a given covariate, hazard (non-survival) is increasing as that covariate increases. If less than 0, hazard is decreasing and survival likelihood increasing as that covariate decreases. The exponentiated coefficient is the hazard ratio and is interpreted relative to 1.0, not 0, as discussed above. Since a case with shorter time to event is more likely to incur the event, hazard rates are often also interpreted as probabilities of the event occurring for a given case. By the same token, hazard rates do not estimate absolute effects, only relative risk.

    • The cumulative hazard function is the non-survival rate, computed as the integral of hazard functions for all time periods.

    • Confidence intervals can be computed around a hazard ratio.

  • Median ratios. Because the hazard ratio shows relative but not absolute effect, other time-based measures may be used to assess the magnitude of the effect on duration times. The ratio of median times is the obvious candidate for such a measure. For instance, in a study of the effect of zinc lozenges on duration of common colds, the ratio of median healing times between zinc lozenge and placebo lozenge groups would measure the effect of zinc lozenges on absolute duration times. Life tables, discussed below, are among the means to compute median times for use in median ratios.

  • Life tables. Related but not part of SPSS's Cox module are life tables, accessed in SPSS by selecting Analyze, Survival, Life Tables. The life table will give the number entering and exiting the risk pool in any time interval, number exposed to risk, number terminating, proportion terminating, proportion surviving, cumulative proportions, and the hazard rate and its standard error for each time interval. Unrau & Coleman (2006), for instance, use life tables to analyze hazard rates for child abuse in terms of time from discharge from a social services program, suggesting the policy application being deciding how long after program termination to schedule follow-up social worker visits.

  
  
• Model fit using the likelihood ratio, AIC, and residual analysis is discussed in the section on event history analysis.

• Models. Various Cox regression models exist to fit various sets of data assumptions/situations.
  • Time-constant Cox regression models. In these models, covariates are constant over time for any given subject/observation (ex., gender=1 or driverstatus = 1). In SPSS, this option supports types of plots and saving of diagnostic variables not available in the time-dependent model. In SPSS, select Analyze, Survival, Cox Regression. In Stata, the stvary command checks to see if covariates are time-constant or time-dependent, and the stcox command executes Cox regression.

  • Time-dependent Cox regression models. In these models, one covariate varies over time and there may be constant covariates as well. This means relative hazard (ratio of observed to baseline hazard) varies over time. Hazards are still proportional over time, but only within the time blocks formed by changes in the covariates. That is, each time a significant covariate changes in value, there is a "jump" up or down in the hazard, but between jumps hazards are proportional. In SPSS, select Analyze, Survival, Cox w/Time Dependent Covariate. In Stata, the stcox command is used in conjunction with the tvc(varlist) option to declare time-varying covariates. Time-varying covariates may be continuously varying (ex., age increases 1 time unit each time t increases by 1) or may be discretely varying (ex., income may go up, be the same, or go down from time period to time period in no set pattern). The regression coefficient, b, remains the same for time-varying covariates but the effect varies according to the magnitude of the variable.

  • Frailty models. Frailty models address the situation where the same individual may experience the hazard more than once, raising the possibility that due to some unmeasured and perhaps unknown cause (that is, a cause of "unobserved heterogeneity"), some subjects may be more likely than others to experience repeated hazards. This likelihood is the "frailty" of the subject and in standard Cox models is an unmeasured effect. Frailty models model the frailty effect as a random effect. Thus frailty models are analogous to regression with random effects. By estimating frailty as a cause of unobserved heterogeneity as a random effect, coefficients for the measured variables are less biased. In addition, of course, the frailty effect (nu) is estimated and may be plotted on the y axis against caseid on the x axis, displaying which cases are the most frail. Frailty is assumed to be constant over time, independent of the covariates, and to be drawn from a given distribution (usually gamma) which the researcher must specify. Frailty models may be badly biased if frailty is correlated with the covariates (Hausman, 1978) or the wrong distribution is assumed (Blossfeld & Rohwer, 1995). Frailty models are supported by Stata but not by SPSS.

  • Conditional frailty models. Conditional frailty models modify frailty models to adjust for event dependence. Simulation studies by Box-Steffensmeier & DeBoef (2006) have demonstrated the superiority of conditional frailty models over standard frailty models under conditions of event dependence. Conditional frailty models stratify cases by event number (1 for the first experience of the event, 2 for the second, etc.). If the frailty variance estimate is significant in a conditional frailty model, then unobserved heerogeneity does affect the model when explaining event dependence. See also Box-Steffensmeier, DeBoef, & Joyce (2007), where conditional frailty models were custom programmed in the R language..
    • Event dependence exists when experiencing the event an earlier time affects the likelihood of experiencing the event a subsequent time. To check for event dependence, plot the cumulative hazard on y by time on x, stratifying by event number. When event depdendence is present, the different strata will display clearly different cumulative hazard function curves in standard Cox models.

  • Repeated events models, also called "multiple episode" models, are directed to situations where events repeat, such as multiple bouts with a disease and remission for patients; ot multiple periods of peace and war for nations. Repeated events are discussed further in the section on event history analysis, but may be implemenedt as Cox models also. Stata but not SPSS supports repeated events models.

  • Competing risks models, also called "multiple destinations" models, are directed at situations where the terminal event may occur for more than one reason. For example, termination of wars may occur through negotiation or defeat. Competing risks models treat different reasons as different events, allowing comparison of hazard functions across competing risks. Competing risks models are discussed further in the section on event history analysis, but may be implemented as Cox models also. Stata but not SPSS supports competing risks models.

  
  

• Stepwise Cox Regression. As in other forms of regression, Cox regression supports "stepwise" as well as "enter" (all model variables entered in one step) and "block" (variables entered in user-specified blocks) methods for entering the independent variables (the covariates). At each step, stepwise methods add the variable with the highest significant score. Also at each step, the residual chi-square is computed and displayed in the "Variables Not in the Equation" table. If the residual chi-square is significant, then at least one of the variables yet to be added to the model is significant.
  • Entry criterion. The score statistic is used by SPSS as the entry criterion. At each step, the variable with the highest score statistic in the "Variables not in the Equation" table is the next to be entered in the following step.

  • Removal criteria In the "Regression if Term Removed" table, a "loss chi-square" statistic is computed at each step, reflecting the contribution of the variable to the model. For any given variable, the variable is removed if the significance of its loss chi-square is greater than .10.

  • "Omnibus Tests of Model Coefficients" table uses -2LL to test change from the prior step, or change from the previous block (if block entry is used; otherwise this will be the same as for previous step). If the overall significance is .05 or less at any step, then at least one of the variables in the model at that point is significant. If the change from previous step significance is .05 or less at any step, then the variable added in that step is significant. If the change from previous block significance is .05 or less at any step, then the variable(s) added in that block is/are significant. In backwards stepwise, where one is removing a variable each step, if the significance of the change > .10, it is conventional to conclude that excluding that variable is justified. This removal criterion is usually based on the likelihood ratio based on maximum partial likelihood estimates, but the user may instead select the likelihood ratio based on conditional parameter estimates (similar but faster computationally), or the Wald statistic.

  
  
• Stratified Cox Regression. By entering a categorical covariate in the "Strata" box of the Cox regression dialog in SPSS, one will obtain separate baseline hazard functions for each value of the categorical variable. One would do this, of course, if one thought that different categories had different baseline functions which were not proportional (if they were proportional, one could use the would-be stratification variable as a covariate; proportionality may be checked by Log-Minus-Log survival plots, discussed below in the section on "Plots"). Because one must assume the same effects in all categories, only one set of pooled coefficients are computed for the covariates (predictors). The stratification variable is not treated as a predictor and no coefficients are computed for it.

  
  

• Cox Regression With or Without Time-Dependent Covariates in Stata. In Stata, Cox regression is executed with the stcox command after declaring a survival time dataset format with the stset command, discussed above. Time-dependent covariates, if any, are declared in the stcox command using the tvc(varlist) option. The sts, stir, and ltable commands generate statistical output related to survival analysis. The stcurve command can be used with stcox or streg to produce survival, hazard, and cumulative hazard function plots, which allow comparison of these functions across levels of covariates. The sts generate command adds new variables to the dataset based on previously-modeled hazard and related functions. In the variants on stcox discussed below, it is assumed that stset has already declared/defined the survival time dataset as described above. Examples are from the Stata manual.
  • Simple Cox regression with uncensored data. Command: stcox load bearings. For a dataset on how long electric generators last till failure, load and bearings are covariates which do not vary by time. The general command syntax is stcox (varlist). The time-to-failure variable is failtime, declared by stset earlier and so not mentioned in the stcox command. All cases (generators) have failed; there are no generators still working (no censored cases). The log likelihood and its probability are also printed. If the probability of the log likelihood is .05 or less, the model as a whole is significant. The main output table will display the hazard ratio, its standard error, its probability level, and its confidence intervals.

  • Cox regression with censored data. Command: stcox drug age. For a dataset on cancer treatment, where drug=1 means the patient received a cancer drug rather than a placebo. Earlier the stset command defined studytime as the time-to-event variable and defined the variable died as the event variable. If died=0, these patients were still alive at the end of the study and constituted censored cases. The more the computed hazard ratio for drug is below 1.0, the more the drug reduced the likelihood of death from cancer, controlling for age. The greater the hazard ratio for age is above 1.0, the more age increased the likelihood of death from cancer, controlling for drug treatment..

  • Cox regression with discrete time-varying covariates. Command: stcox age posttran surg year. For the Stanford heart transplant dataset. There are 1 or 2 records per patient, depending on whether they received a transplant. Earlier, stset defined t1 as the time variable and died as the even variable, and id as the id variable. The variable surg = 1 when the patient had had a previous heart surgery. The variable year was the year the patient was accepted into the transplant program. If posttran = 1, the patient received a transplant and thus posttran is a discrete time-varying covariate. The model is specified the same, however, as earlier examples.

  • Cox regression with continuous time-varying covariates. Command: stcox age, tvc(drug1, drug2) texp(exp(-.035*_t)) nolog. For a dataset on pneumonia, where data for drug1 and drug2 are dosage levels for two drugs respectively, and age is a covariate. Earlier, the stset command defined time as the time variable and cured as the event variable. Had the command been stcox age drug1 drug2, the computed hazard ratios would show the effect of age, drug1, or drug2, each controlling for the other two, assuming the dosage levels of drug1 and drug2 remained constant in any patient's body over time. However, the more complex command with the tvc() and texp() functions can handle the more realistic model assuming that the level of these drugs is time varying, specifically that the amount of the drug in the body diminishes over time. The tvc(drug1, drug2) option declares drug1 and drug2 to be time-varying covariates. The texp(exp(-.035*_t)) option specifies the function defining how the declared time-varying covariates change over time - in this case diminishing exponentially by the function exp(-,35t), where _t = t = analysis time. The nolog parameter specifies the function should not be logarithmic, but that is the default so is not needed (log would specify logarithmic). Output is similar to ordinary Cox regression but the hazard functions are computed differently and the output groups the covariates into non-time-dependent and time-dependent sets. The simple model without the tvc() and texp() functions gives the hazard ratio for, say, drug1 controlling for drug2 and age, and a hazard ratio is the proportionate change in hazard when the dosage level of drug1 increases by 1 unit. The more complex model with the tvc() and texp() functions gives the hazard ratio for drug 1 as a diminishing function of time, controlling for drug2 as a diminishing function of time and for age, and the hazard ratio for drug1 is the proportional change in hazard when blood concentration level (i.e., drug1*exp(-.35t)) increases by one unit.

  • Cox regression with shared frailty . An example given in the Stata manual is an experiment with catheter insertion and possible subsequent infection, with each individual experiencing two insertions (at different times) and thus the possibility of two separate infections. Here the unit of analysis is the insertion, grouped by subject, with presumed shared frailty. Command: stcox age female, shared(patient). In this dataset, patient is the patient id, but it is not used as a conventional id variable but rather as a shared frailty variable. Age and female are continuous and dichotomous covariates respectively. Hazard ratios will be computed differently but interpreted as before. Below the main hazard ratio table, Stata will print a value of theta, its standard error, and a log-likelihood test of theta. If the log-likelihood test of theta is significant (ex., < .05), then there is a significant frailty effect (in this case, a significant patient-level effect in addition to insertion-level effects). If we wish, in a second step we may test to see which patients are the least or most frail (i.e., contribute least or most to patient-level frailty). This would be done with the command stcox age female, shared(patient) effects(nu) followed by the command sort(nu) and list patient nu. This creates a table of patient by nu, which is a measure of frailty. The larger the nu, the more frail the patient, meaning the more likely to experience the hazard.

  • Cox regression with multiple-failure data. We may wish to analyze data in which the event of interest may occur more than once for the same case. Stata supports this and it involves taking each case with multiple failures and creating new cases with new id's, one for each failure. This is done with the stgen, egen, sort, replace, gen, and stset commands as described in StataCorp (2005: 136-138) but not discussed here.

  • Types of variance estimates. Stata supports conventional estimation of variance-covariance matrices by default, and three other user-specifiable alternatives, described in the Stata manual. To get the alternatives, add the options vce(robust), vce(bootstrap), or vce(jackknife) to the stcox command. Alternatively, robust is a synonym for vce(robust). Simply add a comma followed by robust after the stcox varlist.

  
  

• Cox Regression With and Without Time-Dependent Covariates in SPSS. Ordinary Cox regression assumes that any given observation's values on each covariate do not vary over time (ex., "churchattendancerate" of person with CaseID = 437 is the same in each time period; ex., gender=1 does not vary over time for an individual). Cox models can, however, be adapted for covariates which do vary over time for the same individuals. This requires different computation but the output tables are largely the same and interpreted the same. Such models are non-proportional hazard models.
  • Without time-dependent covariates. In SPSS, select Analyze, Survival, Cox Regression.

  • With time-dependent covariates

    • In SPSS, select Analyze, Survival, Cox w/ Time-Dep Cov...; define (optionally transform) in the 'Compute Time-Dependent Cov" dialog box; click the Model button to enter the Time variable, the Status variable (and define its event, ex. death = 1); and enter the Covariates. For time-dependent covariates of type (1) above, select T_COV_ and the covariate and click the ">=a*b>" button, to get a time-covariate interaction term entered into the Covariate list. For time-dependent covariates of type (2) above, click the Paste button to open the Syntax Editor, where you can enter a complex logical expression.

    • Selecting the Cox time-dependent option automatically inserts a time variable, T_, at the top of the variable list. In the 'Compute Time-Dependent Cov' dialog box, you can transform it (ex., for weekly data, T_/52 would transform it to annual) or leave it as T_. Either way, a new variable called T_COV_ is created for use in the analysis.
      1. Ordinary time-dependent covariates. If the variable is systematically related to time (the T_ variable), then an interaction term is created with time (ex., T_COV_*churchattendancerate).
      2. Segmented time-dependent covariates. If the variable is not systematically related to time, then one must create a logical expression which relates the variable to time at each time period (ex., let CA1 be church attendance in Time 1, CA2 in Time 2, etc., to Time 4: (T_ < 1) * CA1 + (T_ >= 1 & T_ < 2) * CA2 + (T_ >= 2 & T_ < 3) * CA3 + (T_ >= 3 & T_ < 4) * CA4). In this example, whichever logical sub-expression (ec., (T_ < 1) is true is evaluated as 1 and multiplied by the corresponding CA variable, and the others are zeroed out.

  
  

• Statistics
  • Hazard ratio, also called the odds ratio or Exp(B). The hazard ratio is the probability of the event occurring in time t + 1, given survival to time t. A hazard ratio of 1.0 indicates the variables in the model have no effect on time to event for the status variable. The more the hazard ratio is below 1.0, then the greater the covariate, the less the odds of the event occurring (increasing predicted survival times). The more above 1.0, the more the variables increase the odds of the event occurring (ex., death=1: diminishing predicted survival times).

  • Relative risk is the hazard ratio for the case where the covariate is a dichotomy, such that when coded 0,1, a 1 indicates presence of a characteristic. In such a case, the characteristic has no influence on the event when its relative risk is 1.0, and increases the probability of the event when its relative risk exceeds 1.0, etc. For instance, if a covariate is smoking (0, 1) with 1 being heavy smoking, and if the hazard ratio (Exp(B)) is 1.1, and if the event is death=1, then the risk of death is 1.1 times greater for heavy smokers than for light and non-smokers (smoking=0), controlling for any other covariates in the model.

  • Confidence intervals on Exp(B) are output by SPSS and other packages, giving the lower and upper 95% confidence limits around the value of Exp(B). If the value 1.0 lies within these confidence limits, one cannot be 95% confident that the covariate has any effect and must report it as non-significant.

  • SPSS. The hazard ratio appears as Exp(B) in the "Variables in the Equation" table in SPSS output.

  • Stata:. As discussed above in the section on simple Cox regression, the stcox command generates the hazard ratio by default, along with its standard error, p value, and confidence intervals.

  • Interpretation of odds ratios is discussed further in the sections on loglinear analysis and logistic regression.

  • Likelihood ratio test of the model, also called the omnibus test or -2LL or -2 log likelihood. If -2LL is significant, the model as a whole is significant. That is, if Sig.> .05 (the usual social science standard), then the covariate effect(s) cannot be assumed to be different from zero. This means that at least one of the covariates contributes significantly to the explanation of duration to event. It also means the model is significantly better than the null model, which is the time-only model when all covariates are 0.
    

    In the figure above, the null (intercept-only) model had -2LL = 45.104. The full model had -2LL = 32.224, a model chi-square difference of 12.88, which is significant at the .012 level. That is, the covariates significantly contribute to explaining days duration of states until ratification of the Constitution, which is the simple example used here..

    • Compared to other tests. The likelihood ratio test is preferred over the Wald test or the score test as ways to assess overall model significance of logistic models.

    • The score statistic. SPSS and some other packages also output a score statistic (aka, global chi-square or overall chi-square) as an alternative significance criterion for the model, but the likelihood ratio test is the standard test. However, score is used in stepwise Cox regression in SPSS, which at each step adds the variable with the highest significant score. For the illustration above, the significance level is the same (.012) whether by likelihood ratio chi-square or global (score) chi-square.

    • SPSS:. The -2 log likelihood statistic appears in the "-2 Log Likelihood" table in SPSS. Stepwise likelihood ratio tests appear in the 'Model Coefficients' table in SPSS.

    • Stata: The likelihood ratio test is generated by the stcox command.

    • Likelihood ratio tests are discussed further in the section on logistic regression

  • Regression coefficients. Most statistics packages display the regression coefficient (B) for each covariate, the standard error of B (SE), its Wald test significance value, the degrees of freedom (df), and the significance value of the coefficient, all similar to and interpreted as in logistic regression. SPSS does this in the "Variables in the Equation" table illustrated below. If Sig.> .05 (the usual social science standard), then the covariate effect cannot be assumed to be different from zero. That is, if sig(Wald) < .05, then the researcher concludes the variable is useful to the model. Positive regression coefficients mean the covariate increases hazard (an increased probability that death=1, for ex.), while negative coefficients correspond to reduced hazard. Note that the -2LL test is preferred over the Wald test when testing the model overall.
    
    • Example. In the illustration above, days until state ratification of the Constitution is predicted from various covariates. In Model 1, closeness of the vote (VotePct) is a significant predictor. But in Model 2, when a binary predictor is added, reflecting whether or not the state was heavily involved in the struggle to include a Bill of Rights in the Constitution (Rights), the Rights variable becomes significant and controlling for Rights, VotePct becomes non-significant. (Of course, with all 13 original states in the data, the data are not a sample and any effect, no matter how small, is not due to chance of sampling).

      Rights was coded 1=heavily involved, 0 = not, so the negative sign of the Rights coefficient in Model 2 signifies reduced hazard (of the event ratification), which translates into more days until ratification. Likewise, the positive sign of the VotePct covariate in Model 1 signifies increased hazard of ratification, which equates to fewer days until ratification. That is, the closer the vote (lower VotePct) or the being in the higher Bill of Rights category (1=heavily) both tend to increase days until ratification when considered separately, but when considered together, Rights controls for VotePct.

    • Odds ratio. Exp(B) in the output above is the odds ratio, which is also the hazard ratio for a given covariate. Exp(B) is the predicted change in the hazard for a unit increase in the predictor. Odds ratios of 1.0 mean the covariate is of no effect on the odds associated with the dependent. Odds ratios above 1.0 are associated with increased hazard of the event (in this case, sooner ratification and hence fewer days). Odds ratios below 1.0 are associated with decreased hazard of the event (in this case, later ratification and more days duration). Thus odds ratios above 1.0 will correspond to positive b coefficients, and odds ratios below 1.0 will correspond to negative b coefficients. For categorical covariates, one must interpret with respect to the reference category. For instance, in the illustration above, Size is categorical (1=small states, 2=medium states, and 3 (the reference category) = large states. The odds ratio for small (1) and medium (2) states are both less than 1.0, indicating that compared to large states (3, the reference category), being small or medium increased the hazard, meaning small and medium states ratified sooner. Thus large states took most days to ratify the Constitution. Medium states ratified in fewest days, indicated by the odds ratio for Size=2 being the lowest (furthest from 1.0).

    • Correlation matrix of regression coefficients is one of the tables in SPSS output and output from other packages. It is used to check for multicollinearity. Ideally, no pair of predictors is highly correlated. The illustration below shows multicollinearity is not indicated for the predictors in Model 2, discussed above.
      

    • Categorical variables. In SPSS and other packages, if a variable is categorical, then there will be an overall row (ex., "religion") as well as a row for each non-omitted value (ex., "religion(1)", "religion(2)", etc.). The overall row will not have entries for B, SE, Exp(b), or confidence intervals, but will have a Wald value and the corresponding significance (see below). This overall Wald significance tests the null hypothesis that all the effect coefficients for that categorical variable are zero. If the overall Wald significance is .05 or less, the researcher may conclude that at least one of the effect coefficients differs from zero.

    • SPSS: The unstandardized regression coefficient, "B" appears in the "Variables in the Equation" table in SPSS.

    • Stata: After fitting a Cox regression model with stcox, it can be reinvoked to display regression coefficients rather than hazard ratios: stcox. nohr. Output will be very similar to the default, but with coefficients substituted for hazard ratios.

  • Baseline hazard, survival, and cumulative hazard rates.

    • Baseline cumulative hazard is the hazard rate for the time-only model when all covariates = 0. As illustrated below, rather than being a single rate, the baseline cumulative hazard function is displayed for various times represented as rows, each with a corresponding baseline cumulative hazard rate. Generally this rate will increase with time. Baseline cumulative hazard rates are easier to interpret when numerical covariate data have been standardized (not the case illustrated below), causing covariate means to be zero, with baseline hazard rates interpreted as time-only rates for people at the mean of the covariate(s). When the covariates are all categorical, then the rate will be for people in category "0" for each covariate.
      

    • Survival rates. The "Survival" column in SPSS output gives the estimated survival rate for the row-specified time after the onset of risk of the event, for persons or other units of observation at the mean of the covariate(s). In the example above, onset of risk is the adoption of the Constitution, states are the units of observation, and the event is ratification of the Constitution (the "risk" event). The survival rate is the estimated percentage of cases which have not experienced the event of interest by the row-specified time. In the illustration above, by day 224 after adoption of the Constitution in Philadelphia, 36% of states had not ratified it.

    • Cumulative survival plots. In this plot, illustrated previously above, supported by SPSS and other packages, the X axis is still survival time. The Y axis, however, is cumulative survival. The curves represent a hypothetical individual(or other unit of analysis) with mean values on the covariate(s) at any given time as represented on the X axis. The curve(s) show how cumulative survival decreases over time for such hypothetical individuals. If a categorical variable has two values (ex., 0=non- or light smoking, 1=heavy smoking), then there will be a cumulative survival plot for each value, allowing comparison. Cumulative survival plots slope down left to right since cumulative survival decreases as survival time increases.

    • Cumulative hazard rates. The "Cumulative hazard" column in SPSS output is similar, but for the model where both time and the covariate (ex., age) are predictors. Mathematically, the cumulative hazard is the negative of the log of survival.

    • Cumulative hazard plots. In this plot, also illustrated previously above, supported by SPSS and other packages, the X axis is survival time. The Y axis is cumulative hazard. The curves represent a hypothetical individual (or other unit of analysis) with mean values on the covariate(s) at any given time as represented on the X axis. The curve(s) show how cumulative hazard increases over time for such hypothetical individuals. If a categorical variable has two values (ex., 0=non- or light smoking, 1=heavy smoking), then there will be a cumulative hazard plot for each value, allowing comparison. Cumulative hazard plots slope up left to right from the origin on the left.

    • SPSS: As discussed above, the "Survival Table" in SPSS output contains the baseline and predicted hazard rates. Rates are presented for a hypothetical person who scores at the mean on the covariate(s). Rows in the table are intervals of time (ex, .5 years, 1.0 years, 1.5 years, etc.). SPSS outputs a "Mean of covariates" table, giving covariate means (ex., mean age in a study of deaths from disease, or mean productivity score in a study of promotions).

    • Stata: In Stata, the sts graph command generates graphs of the Kaplan-Meier survival functions (st graph), the Nelson-Aalen cumulative hazard function (st graph, na), the estimated hazard function (st graph, hazard), and more, based on an analysis previously computed with the stcox command. These are plots with analysis time on the x axis and the survival estimate (1 to 0) or the cumulative hazard estimate (0 to 1) on the y axis. If the "by" parameter is added, the graph will display two or more function curves, one for each value of the by variable (ex., st graph, by(drug) will give a function plot for drug=0 and drug=1 if drug has two values). That is, the "by" parameter enables comparison of survival or hazard functions across levels of a discrete covariate. If one wants survival or hazard functions as a table rather than a graph, this can be accomplished with the sts list command, which also can take a "by" parameter.

  • Test of equality of survivor functions. Stata, using the sts test command, performs either the log-rank test (the default) or alternative tests (the Wilcoxon (Breslow) test of differences in survivor functions by group, the Cox test of equality, the Tarone-Ware test of equality, the Peto-Peto-Prentice test of equality), the generalized Fleming-Harrington test of equality). The maximum number of groups is 800. For instance, assuming an already-executed stcox command, the command sts test gender, wilcoxon gives the Wilcoxon test for equality of survivor functions by gender. Whatever test is chosen, one gets a chi-square value and the significance for that value. If the significance (p value) is =< .05, then the groups differ significantly by survivor function. The default log-rank test is appropriate when each failure time is to be given equal weight,as when the researcher thinks hazard functions are proportional across groups. The Wilcoxon test weights by number of subjects in the risk pool at the time of failure and is appropriate when the researcher thinks hazard functions vary non-proportionally across groups but censoring patterns are similar across groups. The Tarone-Ware test weights by the square root of the number of subjects remaining in the risk pool at the time of failure and is similar in assumptions to the Wilcoxon test; both weight earlier failure times more heavily, but Wilcoxon more so. The Peto-Peto-Prentice test is appropriate when hazard functions are assumed to be non-proportional across groups, but this test is not affected by similarity or dissimilarity of censoring patterns across groups.

  • Parametric survival models In Stata, the streg command supports survival models other than the Cox proportional hazard model. These are the exponential, Weibull, log-normal, log-logistic, Gompertz, and gamma models. Output for these models still contains a log likelihood test of the model as a whole, hazard ratios for the covariates, and also a test of the parametric shape parameter assumed for the baseline hazard function (this is reported as the p value for "/ln_p" following the hazard ratios.

  • Pattern plots. The Plots button in the SPSS dialog allows the user to specify categorical predictors to create pattern plots, where the patterns are the hazard or survival functions plotted separately for each level of the categorical variable, as illustrated below.
    

    In the example above, the categorical predictor is Size, referring to whether a state is small, medium, or large. With respect to the example of ratification of the Constitution, a major debate at the Constitutional Conventions concerned compromises between large and small states, so it may be of interest to compare states with regard to the hazard function (where the "hazard" is ratification of the Constitution, and time is days until ratification). In the illustration above, we see that the predicted hazard functions were indeed different for small, medium, and large states. (Note Size was not a significant predictor of duration once other variables were controlled, but then since the data are an enumeration of all 13 original states and not a random sample, significance does not have its normal meaning and relevance).

    

  • Outlier analysis with DfBeta. The Save option of the Cox regression dialog in SPSS, illustrated above allows computation of the DfBeta statistic for each case for each variable. DfBeta is a measure of how much a given case will affect the regression coefficient for a given covariate. The higher the DfBeta for a given case for a given covariate, the more removing that case from the dataset will change the b coefficient for that covariate. That is, the higher the dfBeta, the more the case is "influential" for that covariate. Influential cases can be spotted visually by plotting DfBeta for a variable against the Case ID variable. High DfBeta values may flag coding errors or sampling errors, or may call attention to clusters of cases requiring a different model.
    

    The saved DfBeta output above has four DfBetas to represent the four terms in Model 2 of the example used in this module, in order as entered in the model: VotePct, Size(1), Size(2), and Rights. In Model 2, Rights was the only significant predictor. In the "Variables in the Equation" table above, Rights had a parameter coefficient of -4.227. DFB4_1 estimates the change in this coefficient if that case is removed. Removing SC would have the most effect in a negative direction. Removing NC would have the most effect in a positive direction. A positive direction corresponds to increase risk (of ratification, the status event) and thus fewer days duration to ratification. NC had among the highest durations, so removing it would leave a dataset with fewer days duration on the average. However, the highest positive DfBeta does not necessarily correspond the the unit with the highest time score (that would be RI, which is not an outlier by the DfBeta criterion). Rather, DfBeta reflects effects on duration to event of a particular variable after other variables in the model are controlled. .

  • Plots. In addition to statistical output, the Plots button in Cox regression in SPSS supports cumulative hazard, cumulative survival, log-minus-log, and partial residual plots. Use of these plots is discussed above in the "Baseline hazard, survival, and cumulative hazard rates" section and below in the "Assumptions" section. The Plots button dialog for SPSS is shown below.
    

    
    

Assumptions

• Assumption of proportional hazards. Cox regression with time-invariant covariates assumes that the ratio of hazards for any two observations is the same across time periods. For instance, in a time-invariant Cox model the ratio of hazards for persons a and b should be the same this year as in the period 10 years from now. This can be a false assumption, as when 10 years from now person B is in their 70's, when mortality spikes, considering age as the covariate. This is a critical assumption of Cox regression and must be checked for each covariate. Gray (1996; quoted in Box-Steffensmeier & Zorn, 2001: 974) has reported as much as a 90% reduction in the power of significance tests (power = chance of false negatives, rejecting the existence of true covariate effects) when rates cross rather than are proportionate.

  If a covariate fails this assumption, then for hazard ratios that increase over time for that covariate, relative risk is overestimated (that is, for diverging hazards, coefficient estimates are inflated). For ratios that decrease over time, relative risk is often underestimated (that is, for converging hazards, coefficient estimates are deflated and biased toward zero). ["Converging" means that the hazard rates for two groups formed by a covariate factor are tending toward the same rate over time]. Correspondingly, standard errors are incorrect and significance tests are decreased in power (Box-Steffensmeier & Zorn, 2001: 972). It is common for a covariate to fail the assumption of proportional hazards, and the implication for estimation should be reported. There are alternative ways to check:

  • Partial residual plots, Graphical methods are used to examine covariates. This is a plot scaled Schoenfeld residuals on the y axis against time on the x axis, with one such plot per covariate. A lowess smoothing line summarizing the residuals should be close to the horizontal 0 reference line for the y axis, since the average value of residuals at an tiime should be zero if the effects of the covariate being plotted are proportional (see Box-Steffensmeier & Zorn, 2001: 978-981). Partial residual methods are the most common and preferred methods for testing for non-proportionality in Cox models.
    • In SPSS select "Partial residual plots" under the Plots button after first having saved partial residuals by checking "Partial residuals" in the "Save New Variables" dialog box under the Save button in the Cox regression dialog. The X axis is survival time. The Y axis is the partial residual for a given covariate. In a well-fitting model, distribution of residuals over time is random. This can be checked further in the Chart Editor by adding a loess smoothing line or linear regression line to show non-random trends. If random, fit lines should not diverge much from the Y-axis 0 reference line.

    • The Grambsch-Therneau global test of non-proportionality uses partial residuals. If the Grambsch-Therneau coefficient is non-significant, at least one of the covariates is non-proportional.

  • Martingale residual plots. If martingale residuals on Y are plotted against the linear predictor (the right-hand side of the model equation) on X, there should be no pattern of correlation if the proportional hazards assumption is met. SPSS does not save martingale residuals directly but they may be computed as mresid = event-haz_1, where event is the event variable and haz_1 is the variable saved under the Save option for the cumulative hazard function. The Save option saves the linear predictor values under the default variable name of X'Beta where XBeta is linear combination of mean corrected covariates times regression coefficients from the final model.

  • Survival probability plots. A plot of cumulative survival on the y axis and analysis time on the x axis may be generated for two or more groups of a covariate. If the lines cross, the covariate violates the proportional hazards assumption. This is an indication that it is a time-dependent variable.
    • SPSS: The covariate in question is entered as a "Strata" variable, not in the Covariates box. Inder the Plots button, select "Survival".

    • Stata: In Stata, the stcoxkm command may be run after defining data with stset but before running stcox (stcoxkm runs Cox itself, for comparision purposes). One can add a "by" parameter (ex., stcoxkm, by(gender) ) to get multiple pairs of predicted/observed curves, one pair for each value of a discrete covariate. There will be one curve for observed and one for predicted values. If the two lines are close together, the proportional hazards assumption is not violated. This tests if the proportional hazards assumption is valid for all groups. The sthplot command is a similar test, also supporting a "by" parameter: if the proportional hazards assumption is valid, the lines for the "by" variable should be parallel and not cross.

  • Log minus log plots (log-log plots or LML plots). Alternatively, use the log minus log test of proportionality. In this test, requested under the Plots button in SPSS, have specified the categorical covariate as the Strata variable. When entered as a Strata variable rather than as a covariate, proportional hazard functions are not enforced for each level of the categorical variable. The cumulative survival estimate after the ln(-ln) transformation is applied to the estimates. The X axis is survival time. The Y axis is log minus log. If the survival plots for the groups of a single categorical covariate are oughly parallel (and certainly should not cross), then the baseline survival functions are parallel and the researcher rejects the need to conduct stratified Cox regression. That is, if there is no violation, the hazard function lines or LML for each category should be parallel. Intersecting survival, hazard function, or LML lines indicate clear violation of the assumption of proportional hazards. The LML method is not recommended for multiple covariates or when a covariate is continuous (see Box-Steffensmeier & Zorn, 2001: 975-976).

  • Time interaction test. If the assumption of proportional hazards is not violated for a given numerical (continuous or categorical) covariate, then the interaction term between that covariate and time (ex., age*time or more commonly, age*log(time)) can be added to the model as in regression, and should have a regression coefficient not significantly different from zero. If the time interaction effect is significant for a covariate, then the proportional hazards assumption is violated and the covariate should be modeled as time-dependent. In the SPSS option for "Cox Regression with a Time-Dependent Covariate" one may add time-covariate interactions to the model and if the interaction is not significant, then the covariate in question is not time-dependent and would not violate the proportional hazards assumption in Cox regression.

  • Categorical covariates. For a given categorical covariate, one may compute the baseline hazard function for each category of that covariate. The shape of the baseline hazard functions should be similar if the assumption of proportional hazards is not violated. In SPSS, click Plots, check the "Hazard" checkbox, and enter the categorical covariate in the "Separate lines for" textbox.

  • Piecewise regression method. Though considered an imprecise "rule of thumb" method, one may divide the sample into observations above and below the median survival times, then model each sample separately to see if the estimated hazard ratio for each covariate coefficient is the same, thus supporting the proportional hazards assumption.

  • Harrell's rho. A rho coefficient may be computed for each covariate. A significant rho means that covariate violates the proportional hazards assumption.

  • Relation to Cox regression with time-dependent covariates. If a covariate fails the test of proportional hazards, this is evidence that it is time-dependent and one needs to abandon ordinary Cox regression in favor of Cox regression with time-dependent covariates. Alternatively, one may include a time-covariate interaction term in the model. As a third alternative, the covariate may be entered as a Strata variable, but then regression coefficients are not computed so this is only an option when the covariate is not of research interest.

• True starting time. The ideal model for survival analysis would be manufacturing of a motor or light bulb, where there is a true zero time (the time of manufacture, before which failure is logically impossible). In medicine, the true zero point is often birth, before which death from a disease is impossible. However, the true zero time in other analyses may be less clear. This is the case in most time-to-adoption studies, where what is adopted is an innovation or piece of legislation. If the zero point is arbitrary or ambiguous, this means that the data series will be different depending on starting point and hence the computed hazard rate coeffiicients will differ, perhaps markedly. (Except if there are no data on predictor variables for, say, years 1900-1980 in a study of 1900-2000, the coefficients will be the same as for a study of 1981-2000). If there is ambiguity about the true starting time, at a minimum the researcher should conduct a sensitivity analysis to see how coefficients may change according to different starting points for data on the predictor variables. Sensitivity analysis may or may not lead the researcher to conclude that Cox modeling is inappropriate.

• Clearly defined events. The status variable must be unambiguously defined, so that any subject for any time period, that subject may be clearly assigned (usually to status=0 or status=1, depending on whether the event of interest had occurred). "Ordinary" Cox modeling deals with situations where the unit of analysis has a risk of a particular event, not a series of different types of events. If the actual data have events representing multiple states, the researcher should use more complex multiple event models for such data. For example, the model for "diplomatic resolution" may be very different from the model for "military resolution," and subsuming both under "conflict resolution" will yield estimates unsatisfactory for explaining either state when using usual single-event Cox models. Rather, multiple states should be modeled explicitly.

• Absence of outliers. As in nearly all forms of analysis, outlier cases can bias estimates. See above for a discussion of statistical output used in analysis of possible outliers.

• No small samples. Precision of parameter estimates using the partial likelihood methods employed in Cox models can be much less than for maximum likelihood methods employed in parametric event history models. Therefore, according to Box-Steffensmeier & Jones (1997: 1434), "this [Cox] method should not be used with small samples."
• Proper model specification. As in other forms of regression, the Cox (and parametric EHA models) regression coefficients may change substantially and even reverse direction if previously omitted relevant variables are added to the model, or if irrelevant but correlated variables are removed from the model. In an event history analysis (EHA) context, this is discussed as the problem of unobserved heterogeneity, meaning bias introduced by omitting important explanatory variables. Unobserved heterogeneity (the effects of variance in important but unobserved variables) is associated with downward bias in duration dependence (Vermunt and Moors, 2005: 10).
  • Unobserved heterogeneity also biases estimates of covariate effects. Unobserved heterogeneity may be addressed by including random effects in the EHA model, using a time-constant latent covariate. See Heckman and Singer (1982) on random effects procedures. See Box-Steffensmeier & Jones (2004: Chapter 9) for an extended discussion of unobserved heterogeneity and ways to deal with it, including frailty models (including a random parameter in the hazard rate to represent unmeasured risk factors) and split-population models (dividing observations into a sample that will never experience the event and a sample at true risk of experiencing the event, so as to avoid the bias that arises when large numbers of cases are not "truly" at risk).

  • Model-trimming strategies apply to Cox and other EHA models. Since coefficient size is a function of other covariates in the model, including inappropriate covariates, when the hazard ratio for some covariates is found to be non-significant, the researcher should drop the most non-significant covariate from the model, re-run the analysis, and proceed stepwise until there are no more non-significant covariates in the final model. It is possible that substantive interpretation of the final covariates may differ from inferences that might have been made for the model including non-significant covariates.

• Few ties. Because Cox methods rely on order of events, handling ties poses a computational problem. Although there are methods for handling ties (ex., the Breslow method is most often used), as a rule of thumb, there should be 5% or fewer tied observations in the dataset to still assume insignificant bias (Prentice & Farewell, 1986: 14).

• Independent observations. EHA models assume the observation's event status at one point in time does not predict the observation's event status at a subsequent point in time. Lack of independence leads to error terms being correlated, which in turn leads to biased estimates of standard error and significance. This is the problem of autocorrelation in time series analysis. Since in diffusion studies and many other social science areas independence is not a sound assumption to make, time dependence must be modeled. For instance, a time variable may be included as a covariate, and/or a surrogate variable added such as number of neighbors adopting the innovation in a diffusion study.
  • Robust estimators. When data independence is an issue, the problem can be mitigated by the use of robust variance estimation, which relaxes assumptions about the distribution of error terms. Less commonly, clustered standard errors may be computed, which relaxes assumptions of independence even further. Clustering accounts for serial time or spatial dependence by dividing the data into groups defined by a grouping variable, then computes standard error across clusters.

    With time series data, one should assume that it is quite possible that data will be temporally dependent (the value at time t+1 is partly a function of its value at time t). This is related to the autocorrelation problem in time series analysis. In Cox and other EHA models, however, the researcher need not de-trend the data but only use "robust variance estimation", which refers to algorithms by Lin & Wei (1989) and Huber (1967) to adjust standard errors for time dependency. Robust estimation is usually the default in Cox and EHA software. It results in the same parameter estimates as standard variance (algorithms assuming independence) but higher standard errors. That is, robust estimation increases the possibility that parameters will be found to be non-significant. Robust estimation is recommended for parameter estimation for time-dependent covariates unless the researcher can demonstrate lack of time dependency in the data (that is, robust estimation should be used most of the time).

• Not applying single-event models to multiple event data. "Ordinary" Cox modeling deals with situations where the unit of analysis has a risk of an event, and after the event occurs the unit drops out of the risk pool. If the actual data are multiple event in nature, meaning the unit is still at risk even after the event occurs the first time, the use of single-event models in serial fashion assumes that experiencing an event in the past has no influence on experiencing the event in the future. Since this assumption may well not be met, hazard ratios will be biased. Instead, the researcher should use more complex multiple event models for such data.

• Exogenous covariates. Interpretation of hazard ratios in models with time-varying covariates assumes those covariates are exogenous (covariate values may affect duration to event, but duration does not cause the values of the covariate). Box-Steffensmeier & Jones (2004: 112) give the example of casualties as a time-dependent covariate in a model of war duration: the Cox model assesses if casualties affect war duration, but if war duration also causes casualties (as it would), casualties is an endogenous covariate and Interpretation of hazard ratios will be biased for that covariate. Unfortunately, there is no accepted method of dealing with endogenous covariates, yet leaving out a causally important endogenous covariate may be a form of model misspecification, with equally problematic implications. The researcher is forced to choose the path of lesser evil. This problem is not unique to Cox modeling.

• Factor invariance. It is assumed that the causal factor structure is the same at the end as at the beginning of the study period.

• Baseline distribution of survival times. Cox regression does not assume any particular baseline distribution of survival times, unlike parametric survival analysis models such as Weibull models, exponential models, and other models as found in Stata's streg procedure (see above). Put another way, Cox regression does not assume any particular distribution shape for the duration times of events. This is because the "dependent variable" in Cox regression is not the event or time to event, but rather the hazard rate. As such Cox regression is more robust than parametric models if the other assumptions of Cox regression are met.

• Hazard rate linearity. If time is measured by a counter variable in units of time, the Cox model assumes that the hazard rate increases linearly with time, conditional on covariates in the model.

• Log linearity. Covariates are assumed to be linearly related to the log of the hazard function. This is tested by running the model without the given covariate, then computing martingale residuals and plotting them on the y axis against the omitted covariate on the x axis. If the loess smoothing line through the scatterplot is close to linear, there is log linearity for that covariate. If log linearity is not present, one may have to transform the covariate (ex., use the square). SPSS does not save martingale residuals, but it does save values for the cumulative hazard; then use Transform, Compute Variable to compute martingale = event - Haz_1, where event is the event variable and Haz_1 is the saved cumulative hazard variable.

• No high multicollinearity. Not having high multicollinearity is an assumption of Cox regression, as in other forms of regression. The "Correlation Matrix of Regression Coefficients" table in SPSS output checks this. If there are multiple highly correlated covariates, one strategy is to include in the model only one variable from the set of intercorrelated variables.

• Random sampling of data is assumed. This is discussed further in the FAQ section below.

• No censoring patterns. Censored cases must be independent of the survival distribution. Censored data are the cases where the event never occurs (where the status variable remains equal to 0) for all time periods. There should be no pattern to these cases, which instead should be missing at random. For example, it could be all censored cases in a public policy study are cases which were ineligible for policy benefits, thereby affecting their status on the status variable. When there is no patterning, subjects entering in different time cohorts should be similar on the average.

Example of SPSS Cox Regression Output

• Cox Regression output from SPSS 14

Frequently Asked Questions

• Why can't we just use OLS or logistic regression to analyze time until event data?
  There are four main reasons:
  1. Censored data are not handled by traditional methods. Any given dataset on, say, disease and death, will contain data on people who have the disease and died (the uncensored observations) and people who have the disease but who have not yet died (the censored observations, meaning that the data on how long they will live is not yet known). Censored observations occur in all time to event data unless the data are historical, with all data present for all observations. Traditional regression methods would require either dropping censored cases, thereby risking sample selection bias, or treating censored cases the same as those for whom the event (ex., death) occurred in the final time period. thereby also biasing computed coefficients. Whereas the usual regression model uses ordinary least squares or maximum likelihood estimation of parameters, Cox regression uses partial likelihood methods, which do not assume uncensored data. In Cox regression, the computation of the regression coefficients is based only on the uncensored cases, but all cases are used when estimating the baseline hazard. Thus Cox regression uses all available information and is considered a full information method, whereas OLS and logistic regression are partial information methods when censored data are present.
  2. Time varying independent variables can be handled in Cox regression but not in traditional regression.
  3. Event distribution. Traditional regression requires events be well distributed over time. However, event analysis frequently centers on the analysis of rare events which are not well distributed. With such data, most time periods will have a value of zero. The large numbers of zeros may inflate correlations and parameter estimates in traditional regression, which is why techniques such as Poisson regression are more appropriate.
  4. Full effect vs. net effect. Also, OLS and logistic regression on cross-sectional data yields effect sizes which show net effect. When the variable in question "cuts both ways," effects could even cancel out, yielding an effect size of zero and the erroneous conclusion that the variable did not matter (ex., education in the short term decreases the likelihood of being employed but in the long term increases the likelihood of being employed; for a population ages 18-22, it is possible that these cross-cutting effects would cancel out). Even when there is not complete self-canceling, traditional regression on cross-sectional data lacks the ability to establish causal direction of the net effects. Moreover, if the process being studied is not stable over time, even the net effects found in cross-sectional regression will be misleading (if, for instance, the proportion of unemployed to employed varies greatly over time in a study of the effects of eduction on employment). To study such processes, one needs a time series method which traces the individual through various states (ex., unemployment, employment) over a long period of time during which various relevant events (ex., degree completion) may occur.
  5. For further discussion of why event history methods like Cox regression or survival analysis are preferred to using logit or probit regression in diffusion studies, click here.

• When would one use a parametric event history analysis model rather than a Cox model?
  Only when one has strong theoretical reasons for positing a particular distribution (shape) for the baseline hazard function, which is very rarely. Cox models can generate the same information as EHA parametric models without having to make as strong data assumptions. Box-Steffensmeier & Jones (2004: 66) write, "there are few instances we can think of where one would naturally prefer a parametric duration model over a Cox-type event history model for most kinds of social science applications." Also, as noted by Buckley & Westerland (2004), in the past there has been over-reliance on parametric EHA models such as logit and probit which, unlike Cox models, require the researcher to specify duration dependence (the shape of the hazard rate over time), misspecification of which can lead to computed standard errors which are too large or too small, leading to errors of inference.

• Describe data setup for Cox regression
  • Event history data setup. If there are no time-varying covariates, event history data setup has a code for the unit of analysis (ex., states), a 0-1 code for whether the risk event (the "censoring variable") ever occurred for that unit, a duration (time elapsed) count variable (ex., 12 for 12 time units since onset of risk), and one or more columns for the time-fixed covariate(s). If there is a time-varying covariate, there must also be a duration (time periods to event) variable and each unit of analysis must have a separate row for each time the covariate changes in value. More often, however, the counting process data setup is used for time dependent data. Of course, data setup must be compatible with the software used. Data setup in Stata is discussed below.

  • Counting process data setup. Time dependent data and multiple event data (where the risk event may occur more than once for the same unit of analysis) are usually entered in "counting process" format. Each unit of analysis (ex., states) has an id code and is represented by as many rows of data as there are time periods (ex., years). There is also a column for the time periods (ex., years from 1980, 1981, etc.). There are also columns for the start or stop interval number, with the start representing the start of risk and coded 0 (thus start/stop 3,4 would be the start of the 4th year of risk and the stop would be the end of the 4th year). The start of "0" would not necessarily be the same actual year for every unit if for some reason the start of risk varied by unit. There would also be a column for the risk event (the "censoring indicator"), coded 0 or 1, with 1 being the occurrence of the event. And there would be additional columns for the time-varying and time-fixed covariates.
    • Discontinuous risk intervals. Under unusual circumstances, some units of analysis might not have rows for some time periods (ex., some years) if risk did not exist for those periods. For example, in a study of peace=0, war=1, a country already at war would not be at risk of war, so that country would have a row for the year war started but would have no further rows until it was again at risk of another war (that is, the first year of peace).

• Cox data setup in Stata:
  • Survival time data (st data) is a data format in which each observation is a time span for a given observation (subject). There is an id variable which indicates the subject associated with the given observation. If the id variable is omitted, it is assumed all observations pertain to a single subject ("single-record st data"). There are variables t and t0, where the span is (t0, t), meaning the period from just after analysis time t0 up to and including time t. If there is no t0 variable, t0 is assumed to equal 0. For instance, the t0 variable might be labeled "Begin" and the t variable labeled "End." There may also be an event (a.k.a. status or failure) variable, d, for each observation, where d=1 if the observation failed (the event occurred) during that particular span but if not, is 0. For instance, the event variable d might be labeled "Died" or "Adopted." If there is no event variable, it is assumed that all observations fail (the event occurs) at time t. There may also be covariate variables (ex., Age, Income). In Stata, Cox regression and related survival analysis is performed on st data.

  • Declaring st data in Stata is the first step in Cox regression in Stata (unless you have ct data, discussed below), accomplished with the stset command. In the menu system, select Statistics, Survival Analysis, Setup & Utilities, Declare data to be survival-time data. Enter the time variable and the failure variable. The stset command declares the researcher's data format. At a bare minimum, the time variable must be declared. There are various variants depending on whether the st data also have an id variable, a t0 variable, or an event variable. Covariates will be utilized if in the dataset and do not have to be declared. Examples are from the Stata manual:
    1. Single record data. Command: stset failure. For a dataset on electric generator lifetimes, with three variables: failtime (the t or analysis time variable), load (a covariate), and bearings (another covariate). Failure is assumed to occur at time = failtime. The t0 variable is assumed to be 0.
    2. Single-record data with censoring. Command: stset failtime, failure(failed). For a similar dataset but with an event variable called 'failed' as well as failtime, load, and bearings. If the event variable is zero (failed = 0) for an observation, failtime is the time, t, at the point of observation and no failure occurred (ex., measurement was stopped at time t and the case is censored, meaning we know the generator will fail at some point after t but we do not know when).
    3. Multiple-record data. Command: stset t, id(patid) failure(died). For a dataset on patients where patid is the patient id; t is the analysis time at the point of measurement; died is the event variable (0 or 1, with 1 = died); and there are covariates also.
    4. Multiple-record data with multiple events. Command: stset day, id(patid) fail(code==402). For dataset with patid = patient id; day = the time variable (t); code = hospital patient status codes, where 402 = death; and there are covariates. So this command is saying day is the time variable, patid is the id variable, and code=402 is equivalent to the event variable = 1. (Note the actual command must have two equal signs).
    5. Multiple record data recording time rather than t (analysis time). Command: stset curday, id(patid) fail(code==402) origin(time adday). For a similar dataset, but with adday containing the day of admission (entered in time units) and curday containing the time variable (number of days since the ward opened in this case). The origin for a given patient is (curday - adday). Analysis is done on curday adjusted for adday. That is, the origin() function converts the time variable into an analysis time variable. Note: in Stata, dates may be displayed in date format but are, in fact, integers, so in this example, curday could be a date variable without any change to the command syntax.
    6. Multiple record data with time from event. Command: stset curday, id(patid) fail(code==402) origin(code==286). Let hospital code 286 mean "patient undergoes operation." Then for the same dataset as above, this command analyzes time from operation until death (code==402). That is, having an operation is considered the onset of risk in this command syntax. In the one above, admission to hospital was considered onset of risk.
    7. Multiple record dataset with delayed entry of observations. Command: stset curday, id(patid) fail(code--402) origina(time adday) enter(code=152). Let hospital code 152 indicate a patient is given a test. The enter() function adds a patient to the analysis only once a record indicates the patient has had the test (code==152). That is, observations for this patient after the one with code=152 will be in the sample.
    8. Scaling time data. Sometimes the original time variable must be adjusted not only for a different origin (ex., the original time variable has 0 as the start of measurement, not the onset of risk) using the origin() function, but also there is a need to convert time units using the scale() function. For instance, the original time variable may be in days, but the researcher wants years. Since there are 365.25 days in a year, adding the function scale(365.25) to the stset command line will rescale the timevar variable as desired. The default is scale(1), which accepts the time units as originally entered.
    9. Other functions. Not discussed here, Stata also provides a number of additional functions such as enter() and exit() for specifying when an observation is in the sample; and if(), ever(), never(), and after() for conditional inclusion,
    10. Temporary re-declaration of st data. Not discussed here, Stata supports the streset command to temporarily re-declare a previously declared st dataset, but with different options (ex., a different definition of time origin).
    11. Other data setup commands. The Stata command stfill is used to fill in missing values, as by carrying forward covariate values from the first observation. The stbase command resets all variables in a multiple-record dataset to the base values in the earliest record for the subject. The command stgen will create new covariates as functions of time variables and covariates (ex., create the variable evervoted if for any time period the variable voted=1). The command stsplit will create multiple records out of one record, as to add a time-varying covariate which will be the same for existing variables but have different values of the covariate for each of the new records.
    12. Error messages. Stata does a certain amount of data format error-checking, with error messages for "event time missing," "entry time missing," "multiple records at the same instant," and "overlapping records," among others.

  • Count time data (ct data) is an alternative data format in which each observation is a time. For each time, there are variables for the number known to fail during that time; the number of right-censored cases (see below); and the number of new cases added during the given time. In Stata, ct data is first converted to st data using the ctset (declare data to be count-time data) and cttost (convert ct data to st data) commands.

  • Snapshot data is a common real-world data format that must be converted to st data format. Rather than have separate t0 (begin) and t (end) variables, there is just a time-of-observation variable. In Stata one converts snapshot data to st data using the snapspan command (syntax: snapspan idvar timevar varlist). Stata will take the earliest timevar and make that t0 and use time units since t0 to create the analysis time, t, variable entries. The varlist variables in the converted st dataset will have the values from the corresponding observations of the snapshot dataset. Any variables not in varlist will have the values carried forward from the t0 record/observation.

  
  

• Why is no intercept coefficient reported for Cox models?
  In Cox models, the intercept is incorporated in the baseline hazard function.

  
  

• Since the Cox model does not posit any particular baseline hazard ratio, how can the baseline hazard function be retrieved?
  The survivor function can be estimated from the order of failure times, the risk at any given failure time, and the assumption of a constant hazard rate between failure times. From the survival function, the hazard function can be derived. See the summary by Box-Steffensmeier & Jones, 2004: 64-65; or the original articles by Kalbfleisch & Prentice, 1973, 1980. Note that the estimate of the baseline hazard function in Cox models is data-driven, whereas in parametric event history analysis models, the baseline hazard function is selected based on theory, or possibly based on comparisons of model fit among several alternative parametric models, each positing a shape of the baseline hazard function.

• Can I use Cox regression with non-random samples?
  Boehmke, Morey, & Shannon (2006) conducted Monte Carlo simulations on this question, coming to the conclusion that "sample selection issues can lead to biased parameter estimates, including the appearance of (nonexistent) duration dependence" (p.192). The authors further found that nonrandom selection of samples could lead to "inaccurate predicted hazard and survival functions" and "erroneous conclusions about what factors influence (or do not influence) the duration process of interest" (p. 205).

Bibliography

• Boehmke, Frederick J., Morey, Daniel S., & Shannon, Megan (2006). Selection bias and continuous-time duration models: Consequences and a proposed solution. American Journal of Political Science 50(1): 192–207.
• Blossfeld, Hans-Peter; Golsch, Katrin; & Rohwer, Götz (2007). Event history analysis with Stata. Mahwah, NJ: Lawrence Erlbaum Associates.
• Blossfeld, Hans-Peter; Rohwer, Götz (1995). Techniques of event history modeling. Mahwah, NJ: Lawrence Erlbaum.
• Box-Steffensmeier, Janet M. & DeBoef, Suzanna (2006). Repeated events survival models: The conditional frailty model. Statistics in Medecine 25: 3518-3533.
• Box-Steffensmeier, Janet M.; DeBoef, Suzanna; & Joyce, Kyle A. (2007). Event dependence and heterogeneity in duration models: The conditional frailty model. Political Analysis 15(3): 237-256.
• Box-Steffensmeier, Janet M. & Jones, Bradford S. (1997). Time is of the essence: Event history models in political science. American Journal of Political Science. 41(4): 1414-1461.
• Box-Steffensmeier, Janet M. & Jones, Bradford S. (2004). Event history modeling: A guide for social scientists. NY: Cambridge University Press.
• Box-Steffensmeier, Janet M. & Zorn, Christopher J. W. (2001). Duration models and proportional hazards in political science. American Journal of Political Science 45(4): 972-988.
• Buckley, Jack, & Chad Westerland. (2004). Duration dependence, functional form, and correct standard errors: Improving EHA models of state policy diffusion. State Politics and Policy Quarterly 4( 1): 94–113.
• Cleves, M. A.; Gould, W. W.; and Gutierrez, R. G. (2004). An introduction to survival analysis using Stata, rev. ed.. College Station, TX: Stata Press. Stata is a preferred software package for cox regression and survival analysis, and this is perhaps the most-recommended Stata text for it.
• Gray, Todd Edward (1996). .The impact of deviations from the proportional hazards assumption on power in the analysis of survival data.. Unpublished M.P.H. Thesis, Department of Biostatistics, Rollins School of Public Health, Emory University.
• Hausman, Jerry A. (1978). Specification tests in econometrics. Econometrica 46: 1251-1271.
• Hosmer, D. W. & Lemeshow, S. (1999). Applied survival analysis. NY: John Wiley & Sons.
• Huber, P. J. (1967) The behavior of maximum likelihood estimates under non-standard conditions. In Proceedings of the Fifth Berkely Symposium on Mathematical Statistics and Probability. Berkeley: University of California Press.
• Jones, Bradford S, & Branton, Regina P. (2005). Beyond logit and probit: Cox duration models of single, repeating, and competing events for state policy adoption. State Politics and Policy Quarterly 5(4): 420-443.
• Kalbfleisch, J. D. & Prentice, R. L. (1973). Marginal likelihoods based on X\Cox's regression and life model. Biometrika 60, 267-278.
• Kalbfleisch, J. D. & Prentice, R. L. (1980). The statistical analysis of failure rate data. NY: John Wiley.
• Klein, J. P. & Moeschberger, M. L. (1997). Survival analysis: Techniques for censored and truncated data. NY: Springer.
• Lin, D. Y. & L. J. Wei (1989). The robust inference for the Cox proportional hazards model. Journal of the American Statistical Association 84: 1074-1078.
• Prentice, R. L. & Farewell, B. T. (1986). Relative risk and odds ratio regression. Annual Review of Public Health 7: 335-338.
• Spruance, S. I.; Reid, J. E.; Grace, M.; & Samore, M. (2004). Hazard ratio in clinical trials. Antimicrobial Agents and Chemotherapy 48(8), 2787-2792.
• StataCorp (2005). Stata survival analysis and epidemiological tables reference manual release 9. College Station, TX: StataCorp LP. Lists all Stata commands and options pertaining to Cox regression and survival analysis, with examples.
• Unrau, Y. A. & Coleman, H. (2006). Evaluating program outcomes as event histories. Administration in Social Work 30(1): 45-65.