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.

# Cite R Packages

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},
#>   }