Risk Score Vignette

Introduction

Risk scores are sparse linear models that map an integer linear combination of covariates to the probability of an outcome occurring. Unlike regression models, risk score models consist of integer coefficients for often dichotomous variables. This allows risk score predictions to be easily computed by adding or subtracting a few small numbers.

Risk scores developed heuristically by altering logistic regression models have decreased performance, as there is a fundamental trade-off between the model’s simplicity and its predictive accuracy. In contrast, this package presents an optimization approach to learning risk scores, where the constraints unique to risk score models are integrated into the model-fitting process, rather than implemented afterward. This vignette demonstrates how to use the riskscores package to build a risk score model to predict breast cancer diagnosis.

library(riskscores)

Optimization Problem

The riskscores package uses a cyclical coordinate descent algorithm to solve the following optimization problem.

These constraints ensure that the model will be sparse and include only integer coefficients.

Loading Example Data

First we’ll load in an example dataset. In this example, we want to develop a risk score model that predicts whether a breast tissue sample is benign using features recorded during a biopsy. The breastcancer dataset was originally accessed from the UCI Repository and can be loaded into your environment from the riskscores package as so:

data("breastcancer")

This dataset contains 683 observations and 9 features. Our goal is to develop a risk score model that predicts whether a breast tissue sample is benign using 9 (or fewer) features recorded during a biopsy:

1. Clump thickness
   
2. Uniformity of cell size
   
3. Uniformity of cell shape
   
4. Marginal adhesion
   
5. Single epithelial cell size
   
6. Bare nuclei
   
7. Bland chromatin
   
8. Normal nucleoli
   
9. Mitoses
   

Data Preprocessing

Before building a risk score model, data often need to be preprocessed. Specifically, the dataset needs to have a binary outcome with all other variables containing either binary or integer values.

The breastcancer dataset is mostly ready to go. We’ll still need to split out our data into a matrix with all covariates (X) and a vector with the outcome data (y). In this case, the first column in our dataset contains the outcome variable.

y <- breastcancer[,1]
X <- as.matrix(breastcancer[,-1])

Cross Validation

The penalty coefficient λ₀ controls the sparsity of the model – a larger value of λ₀ will result in fewer non-zero coefficients. We can use cross validation to find the optimal λ₀ value that creates a sufficiently sparse model without sacrificing performance.

Ideally, each cross-validation fold should contain an approximately equal proportion of cases. The riskscores package contains the function stratify_folds() that creates fold IDs with an equal proportion of cases in each fold. These fold IDs can be entered into the cv_risk_mod() function under the foldids parameter. Otherwise, cv_risk_mod() will set random fold IDs.

foldids <- stratify_folds(y, nfolds = 5, seed = 1)

The cv_risk_mod() function runs cross validation for a grid of possible λ₀ values. If the user does not specify the vector of λ₀ values to test, the program constructs this λ₀ sequence. The maximum λ₀ in this sequence is the smallest value such that all coefficients in the logistic regression model are zero. The minimum λ₀ in the sequence is calculated using the user-defined lambda_ratio argument. The λ₀ grid is created by generating nlambda values linear on the log scale from the minimum λ₀ to the maximum λ₀. We’ve set nlambda to 25, so the program will construct an appropriate sequence of 25 λ₀ values to test using cross validation.

cv_results <- cv_risk_mod(X, y, foldids = foldids, nlambda = 25)

Running plot() on a cv_risk_mod object creates a plot of mean deviance for each λ₀ value in the grid. The number of nonzero coefficients that are produced by each λ₀ value when fit on the full data are listed at the top of the plot. The λ₀ value with the lowest mean deviance (“lambda_min”) is indicated in red, and its standard deviation is marked with a red dashed line. Its precise value can be accessed by calling cv_results$lambda_min. If we want a sparser model, we could increase λ₀ to “lambda_1se”, the largest value whose mean deviance is within one standard error of “lambda_min”. This value can be accessed by calling cv_results$lambda_1se. In our example, “lambda_min” creates a model with 8 non-zero coefficients and “lambda_1se” creates a model with 3 non-zero coefficients.

plot(cv_results)
#> Warning in ggplot2::geom_point(ggplot2::aes(x = log(x$lambda_min), y = min_mean), : All aesthetics have length 1, but the data has 25 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#>   a single row.
#> Warning in ggplot2::geom_linerange(ggplot2::aes(x = log(x$lambda_min), ymin = min_mean - : All aesthetics have length 1, but the data has 25 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#>   a single row.

cv_results$lambda_min
#> [1] 0.0005755162

cv_results$lambda_1se
#> [1] 0.03920946

To view a dataframe with the full cross-validation results (including both deviance and accuracy metrics), run cv_results$results.

tail(cv_results$results)
#> # A tibble: 6 × 6
#>   lambda0 mean_dev sd_dev mean_acc  sd_acc nonzero
#>     <dbl>    <dbl>  <dbl>    <dbl>   <dbl>   <int>
#> 1  0.0576     39.4   8.44    0.952 0.00662       3
#> 2  0.0845     48.6  12.9     0.941 0.0145        3
#> 3  0.124      50.5  15.7     0.937 0.0169        2
#> 4  0.182      60.9  15.7     0.917 0.0208        1
#> 5  0.267      86.5  25.9     0.874 0.0467        1
#> 6  0.392      69.1  10.8     0.900 0.0203        1

Fitting a Risk Score Model

We’ll fit a model on the full data using the function risk_mod(). We’ll use the “lambda_1se” value determined by cross-validation as our λ₀ parameter.

mod <- risk_mod(X, y, lambda0 = cv_results$lambda_1se)

The integer risk score model can be viewed by calling mod$model_card. An individual’s risk score can be calculated by multiplying each covariate response by its respective number of points and then adding all points together. In our example below, a patient with a ClumpThickness value of 1, a UniformityOfCellShape value of 10, and a BareNuclei value of 5 would receive a score of 10(1) + 8(10) + 7(5) = 125.

Each score can then be mapped to a risk probability. The mod$score_map dataframe maps an integer range of scores to their associated risk. For this example dataset, mod$score_map includes a range of integer scores from 25 to 200, which are the minimum and maximum scores predicted from the training data. The table below shows a sample of these scores mapped to their associated risk. We can see that a patient who received a score of 125 would have a 77.65% risk of their tissue sample being malignant.

The function get_risk() can be used to calculate the risk from a given score (or a vector of scores). Likewise, the function get_score() calculates the score associated with a given risk (or vector of risk probabilities).

get_risk(mod, score = 125)
#> [1] 0.9042011

get_score(mod, risk = 0.7765)
#> [1] 114.4883

We can evaluate the model’s performance under different classification thresholds using the get_metrics() function.

get_metrics(mod, threshold = seq(0.1, 0.9, 0.1))
#>   threshold_risk threshold_score  accuracy sensitivity specificity
#> 1            0.1            78.3 0.9531479   0.9832636   0.9369369
#> 2            0.2            86.8 0.9692533   0.9790795   0.9639640
#> 3            0.3            92.5 0.9663250   0.9665272   0.9662162
#> 4            0.4            97.1 0.9663250   0.9539749   0.9729730
#> 5            0.5           101.4 0.9633968   0.9456067   0.9729730
#> 6            0.6           105.7 0.9619327   0.9372385   0.9752252
#> 7            0.7           110.3 0.9546120   0.9037657   0.9819820
#> 8            0.8           116.0 0.9458272   0.8744770   0.9842342
#> 9            0.9           124.5 0.9355783   0.8326360   0.9909910

Generic Functions

summary(mod)
#> 
#> Intercept: -101.3896
#> 
#> Non-zero coefficients:               .
#> ClumpThickness 9
#> BareNuclei     6
#> BlandChromatin 8
#> 
#> Gamma (multiplier):  0.09507679 
#> Lambda (regularizer):  0.03920946 
#> 
#> Deviance:  143.4015 
#> AIC:  163.4015

coef(mod) # equivalently: mod$beta
#>                Intercept           ClumpThickness     UniformityOfCellSize 
#>                -101.3896                   9.0000                   0.0000 
#>    UniformityOfCellShape         MarginalAdhesion SingleEpithelialCellSize 
#>                   0.0000                   0.0000                   0.0000 
#>               BareNuclei           BlandChromatin           NormalNucleoli 
#>                   6.0000                   8.0000                   0.0000 
#>                  Mitoses 
#>                   0.0000

coef(mod) * mod$gamma
#>                Intercept           ClumpThickness     UniformityOfCellSize 
#>               -9.6397973                0.8556911                0.0000000 
#>    UniformityOfCellShape         MarginalAdhesion SingleEpithelialCellSize 
#>                0.0000000                0.0000000                0.0000000 
#>               BareNuclei           BlandChromatin           NormalNucleoli 
#>                0.5704607                0.7606143                0.0000000 
#>                  Mitoses 
#>                0.0000000

coef(mod$glm_mod)
#>                Intercept           ClumpThickness     UniformityOfCellSize 
#>               -9.6397973                0.8556911                0.0000000 
#>    UniformityOfCellShape         MarginalAdhesion SingleEpithelialCellSize 
#>                0.0000000                0.0000000                0.0000000 
#>               BareNuclei           BlandChromatin           NormalNucleoli 
#>                0.5704607                0.7606143                0.0000000 
#>                  Mitoses 
#>                0.0000000

plot(mod, score_min = 25, score_max = 200)