# Survival Modeling

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)
    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.

# References