Here, we’ll look at how to create consonance functions from the coefficients of predictors of interest in a Cox regression model.
We’ll use the
carData package for this. Fox & Weisberg, 2018 describe the dataset elegantly in their paper,
The Rossi data set in the
carDatapackage contains data from an experimental study of recidivism of 432 male prisoners, who were observed for a year after being released from prison (Rossi et al., 1980). The following variables are included in the data; the variable names are those used by Allison (1995), from whom this example and variable descriptions are adapted:
week: week of first arrest after release, or censoring time.
arrest: the event indicator, equal to 1 for those arrested during the period of the study and 0 for those who were not arrested.
fin: a factor, with levels “yes” if the individual received financial aid after release from prison, and “no” if he did not; financial aid was a randomly assigned factor manipulated by the researchers.
age: in years at the time of release.
race: a factor with levels “black” and “other”.
wexp: a factor with levels “yes” if the individual had full-time work experience prior to incarceration and “no” if he did not.
mar: a factor with levels “married” if the individual was married at the time of release and “not married” if he was not.
paro: a factor coded “yes” if the individual was released on parole and “no” if he was not.
prio: number of prior convictions.
educ: education, a categorical variable coded numerically, with codes 2 (grade 6 or less), 3 (grades 6 through 9), 4 (grades 10 and 11), 5 (grade 12), or 6 (some post-secondary).
emp1–emp52: factors coded “yes” if the individual was employed in the corresponding week of the study and “no” otherwise.
We read the data file into a data frame, and print the first few cases (omitting the variables emp1 – emp52, which are in columns 11–62 of the data frame):
library(concurve) library(carData) Rossi[1:5, 1:10] #> week arrest fin age race wexp mar paro prio educ #> 1 20 1 no 27 black no not married yes 3 3 #> 2 17 1 no 18 black no not married yes 8 4 #> 3 25 1 no 19 other yes not married yes 13 3 #> 4 52 0 yes 23 black yes married yes 1 5 #> 5 52 0 no 19 other yes not married yes 3 3
Thus, for example, the first individual was arrested in week 20 of the study, while the fourth individual was never rearrested, and hence has a censoring time of 52. Following Allison, a Cox regression of time to rearrest on the time-constant covariates is specified as follows:
library(survival) mod.allison <- coxph(Surv(week, arrest) ~ fin + age + race + wexp + mar + paro + prio, data = Rossi ) mod.allison #> Call: #> coxph(formula = Surv(week, arrest) ~ fin + age + race + wexp + #> mar + paro + prio, data = Rossi) #> #> coef exp(coef) se(coef) z p #> finyes -0.37942 0.68426 0.19138 -1.983 0.04742 #> age -0.05744 0.94418 0.02200 -2.611 0.00903 #> raceother -0.31390 0.73059 0.30799 -1.019 0.30812 #> wexpyes -0.14980 0.86088 0.21222 -0.706 0.48029 #> marnot married 0.43370 1.54296 0.38187 1.136 0.25606 #> paroyes -0.08487 0.91863 0.19576 -0.434 0.66461 #> prio 0.09150 1.09581 0.02865 3.194 0.00140 #> #> Likelihood ratio test=33.27 on 7 df, p=2.362e-05 #> n= 432, number of events= 114
Now that we have our Cox model object, we can use the
curve_surv() function to create the function.
If we wanted to create a function for the coefficient of prior convictions, then we’d do so like this:
z <- curve_surv(mod.allison, "prio")
Then we could plot our consonance curve and density and also produce a table of relevant statistics. Because we’re working with ratios, we’ll set the
measure argument in
ggcurve() to “
ggcurve(z[], measure = "ratio", nullvalue = TRUE)
ggcurve(z[], type = "cd", measure = "ratio", nullvalue = TRUE)
curve_table(z[], format = "image")
We could also construct a function for another predictor such as age
ggcurve(x[], type = "cd", measure = "ratio")
curve_table(x[], format = "image")
That’s a very quick look at creating functions from Cox regression models.