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
carData
package 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) #> Please see the documentation on https://data.lesslikely.com/concurve/ or by typing `help(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 “ratio
”.
ggcurve(z[[1]], measure = "ratio", nullvalue = TRUE)
ggcurve(z[[2]], type = "cd", measure = "ratio", nullvalue = TRUE)
curve_table(z[[1]], format = "image")
Lower Limit |
Upper Limit |
Interval Width |
Interval Level (%) |
CDF |
P-value |
S-value (bits) |
1.086 |
1.106 |
0.020 |
25.0 |
0.625 |
0.750 |
0.415 |
1.075 |
1.117 |
0.042 |
50.0 |
0.750 |
0.500 |
1.000 |
1.060 |
1.133 |
0.072 |
75.0 |
0.875 |
0.250 |
2.000 |
1.056 |
1.137 |
0.080 |
80.0 |
0.900 |
0.200 |
2.322 |
1.052 |
1.142 |
0.090 |
85.0 |
0.925 |
0.150 |
2.737 |
1.045 |
1.149 |
0.103 |
90.0 |
0.950 |
0.100 |
3.322 |
1.036 |
1.159 |
0.123 |
95.0 |
0.975 |
0.050 |
4.322 |
1.028 |
1.168 |
0.141 |
97.5 |
0.988 |
0.025 |
5.322 |
1.018 |
1.180 |
0.162 |
99.0 |
0.995 |
0.010 |
6.644 |
We could also construct a function for another predictor such as age
x <- curve_surv(mod.allison, "age") ggcurve(x[[1]], measure = "ratio")
ggcurve(x[[2]], type = "cd", measure = "ratio")
curve_table(x[[1]], format = "image")
Lower Limit |
Upper Limit |
Interval Width |
Interval Level (%) |
CDF |
P-value |
S-value (bits) |
0.938 |
0.951 |
0.013 |
25.0 |
0.625 |
0.750 |
0.415 |
0.930 |
0.958 |
0.028 |
50.0 |
0.750 |
0.500 |
1.000 |
0.921 |
0.968 |
0.048 |
75.0 |
0.875 |
0.250 |
2.000 |
0.918 |
0.971 |
0.053 |
80.0 |
0.900 |
0.200 |
2.322 |
0.915 |
0.975 |
0.060 |
85.0 |
0.925 |
0.150 |
2.737 |
0.911 |
0.979 |
0.068 |
90.0 |
0.950 |
0.100 |
3.322 |
0.904 |
0.986 |
0.081 |
95.0 |
0.975 |
0.050 |
4.322 |
0.899 |
0.992 |
0.093 |
97.5 |
0.988 |
0.025 |
5.322 |
0.892 |
0.999 |
0.107 |
99.0 |
0.995 |
0.010 |
6.644 |
That’s a very quick look at creating functions from Cox regression models.
Please remember to cite the packages that you use.
citation("concurve") #> #> Rafi Z, Vigotsky A (2020). _concurve: Computes and Plots Compatibility #> (Confidence) Intervals, P-Values, S-Values, & Likelihood Intervals to #> Form Consonance, Surprisal, & Likelihood Functions_. R package version #> 2.7.7, <URL: https://CRAN.R-project.org/package=concurve>. #> #> Rafi Z, Greenland S (2020). "Semantic and Cognitive Tools to Aid #> Statistical Science: Replace Confidence and Significance by #> Compatibility and Surprise." _BMC Medical Research Methodology_, *20*, #> 244. ISSN 1471-2288, doi: 10.1186/s12874-020-01105-9 (URL: #> https://doi.org/10.1186/s12874-020-01105-9), <URL: #> https://doi.org/10.1186/s12874-020-01105-9>. #> #> To see these entries in BibTeX format, use 'print(<citation>, #> bibtex=TRUE)', 'toBibtex(.)', or set #> 'options(citation.bibtex.max=999)'. citation("survival") #> #> Therneau T (2020). _A Package for Survival Analysis in R_. R package #> version 3.2-7, <URL: https://CRAN.R-project.org/package=survival>. #> #> Terry M. Therneau, Patricia M. Grambsch (2000). _Modeling Survival #> Data: Extending the Cox Model_. Springer, New York. ISBN 0-387-98784-3. #> #> To see these entries in BibTeX format, use 'print(<citation>, #> bibtex=TRUE)', 'toBibtex(.)', or set #> 'options(citation.bibtex.max=999)'. citation("carData") #> #> To cite package 'carData' in publications use: #> #> John Fox, Sanford Weisberg and Brad Price (2020). carData: Companion #> to Applied Regression Data Sets. R package version 3.0-4. #> https://CRAN.R-project.org/package=carData #> #> A BibTeX entry for LaTeX users is #> #> @Manual{, #> title = {carData: Companion to Applied Regression Data Sets}, #> author = {John Fox and Sanford Weisberg and Brad Price}, #> year = {2020}, #> note = {R package version 3.0-4}, #> url = {https://CRAN.R-project.org/package=carData}, #> } citation("survminer") #> #> To cite package 'survminer' in publications use: #> #> Alboukadel Kassambara, Marcin Kosinski and Przemyslaw Biecek (2020). #> survminer: Drawing Survival Curves using 'ggplot2'. R package version #> 0.4.8. https://CRAN.R-project.org/package=survminer #> #> A BibTeX entry for LaTeX users is #> #> @Manual{, #> title = {survminer: Drawing Survival Curves using 'ggplot2'}, #> author = {Alboukadel Kassambara and Marcin Kosinski and Przemyslaw Biecek}, #> year = {2020}, #> note = {R package version 0.4.8}, #> url = {https://CRAN.R-project.org/package=survminer}, #> }