Identifying Actionable Treatment Variation in Episodes of Low Back Pain

Objective

Identify actionable treatment variation within patients with low back pain.

---

Data Preparation

For this analysis, data was pulled from the following tables in the data warehouse.

• V_ETG_MED_CLAIMS
• V_MEDICAL_SERVICE
• V_STD_TOS_CATEGORY
• V_CLAIM_PROCEDURE
• V_ICD10_DIAG_GRP
• V_INCURRED_DATE

Patients were members of BCBS from 2016 who had episodes of low back pain without complications. Ages were between 18-65. ETG outliers were excluded. Members with missing episode numbers were dropped. There were 30,000 episodes included in this analysis out of an original 36,300.

#read in the SAS file 
df = read_sas("C:/Users/scasti01/Desktop/Cost and Use/R/Claims/episodes.sas7bdat") %>%
  #convert data types to proper format
  mutate(MEM_NUM = as.factor(MEM_NUM),
         EPI_TYPE = as.factor(EPI_TYPE),
         etg = as.character(ETG_FULL),
         ETG_OUTLIER = as.factor(ETG_OUTLIER),
         COMPLICATION = as.factor(COMPLICATION),
         COMORBIDITY = as.factor(COMORBIDITY),
         SEVERITY = as.factor(SEVERITY),
         TOS_CAT1_DESC = as.factor(TOS_CAT1_DESC),
         TOS_CAT2_DESC = as.factor(TOS_CAT2_DESC),
         TOS_CAT3_DESC = as.factor(TOS_CAT3_DESC)
         ) %>%
  #exclude members with complications
  #eclude episodes with episodes number 0 (missing values)
  filter(COMPLICATION == 0, EPISODE_NUM != 0, !is.na(EPISODE_NUM) ) %>%
  #find the time difference between each appointment and the start of th episode
  mutate(lag = difftime( INCURRED_DT_DAY, EPI_START_DT, units = "days"),
         surgery_indicator =  ifelse(test = (str_detect(CPT_PROC_DESC, "SURG"))|| any(TOS_CAT2_DESC == "SURGICAL"),
                                                   yes = T, 
                                                   no = F)) %>%
  #count the number of episodes for each patient
  group_by(MEM_NUM) %>%
  mutate( nth_episode = length(unique(EPISODE_NUM))) %>%
  ungroup() %>%
  #create clinical T/F variables for clinical categories
  group_by(EPISODE_NUM) %>%
  summarise( 
            #does the episode have surgery?  these are the ETGs that have surgery at some point in the episode    
             surgery = as.factor(any(ETG_FULL  %in% c(712208001, 712208011, 712208111, 712208101))),
             
             #When does the surgery occur?
             surgery_lag = ifelse( test = any(surgery_indicator),
                                   yes = min(difftime( INCURRED_DT_DAY, EPI_START_DT, units = "days")),
                                   no = 1000),
             #How many physical therapy claims are on record?
             PT_count = sum(ifelse( test = (TOS_CAT2_DESC == "PHYSICAL THERAPY"),
                                    yes = T, 
                                    no = F)),
             
             #did the patient go to more than 7 PT appointments?
             PT = as.factor(any(TOS_CAT2_DESC == "PHYSICAL THERAPY" & PT_count >= 7)),
             #Was there an MRI at any point in the episode?
             MRI_any = as.factor(any(TOS_CAT3_DESC == "MRI" )),
             #Was there an MRI in the first 6 weeks of the episode, and before surgery if there was a surgical outcome?
             MRI_6_weeks = as.factor(any(TOS_CAT3_DESC == "MRI" & (lag <= 6*7) & (lag <= surgery_lag))),
             #was there an MRI in the first 8 weeks and before surgery?
             MRI_8_weeks = as.factor(any(TOS_CAT3_DESC == "MRI" & (lag <= 8*7) & (lag <= surgery_lag))),
             #Did the patient have a spinal manipulation (chiropractic treatment) before surgery?
             spinal_manipulation = as.factor(any(TOS_CAT3_DESC == "SPINAL MANIPULATION" & (lag <= surgery_lag))),
             #did the patient have a diagnostic injection prior to surgery?
             diag_injection = as.factor(any(TOS_CAT3_DESC == "DIAGNOSTIC INJECTION" & (lag <= surgery_lag))),
             #where there xrays taken before surgery?  In the first 6, 8 weeks, and prior to surgery?
             xray_any = as.factor(any(str_detect( CPT_PROC_DESC, "X-RAY") & (lag <= surgery_lag))),
             xray_6_weeks = as.factor(any(str_detect( CPT_PROC_DESC, "X-RAY") & (lag <= 6*7) & (lag <= surgery_lag))),
             xray_8_weeks = as.factor(any(str_detect( CPT_PROC_DESC, "X-RAY") & (lag <= 8*7) & (lag <= surgery_lag))),
             #what was the ETG severity?
             severity1 = as.factor(any(SEVERITY == 1)),
             severity2 = as.factor(any(SEVERITY == 2)),
             severity3 = as.factor(any(SEVERITY == 3)),
             #how old was the member?
             member_age = as.numeric(member_age[1]),
             nth_episode = as.factor(nth_episode[1]),
             MEM_NUM = MEM_NUM[1],
             comorbidity = as.factor(COMORBIDITY[1])
             
  )#end summarise

surgery date test

Do the tos_cat dates agree with the surgery dates?

testdata = 
read_sas("C:/Users/scasti01/Desktop/Cost and Use/R/Claims/episodes.sas7bdat") %>%
  #convert data types to proper format
  mutate(MEM_NUM = as.factor(MEM_NUM),
         EPI_TYPE = as.factor(EPI_TYPE),
         etg = as.character(ETG_FULL),
         ETG_OUTLIER = as.factor(ETG_OUTLIER),
         COMPLICATION = as.factor(COMPLICATION),
         COMORBIDITY = as.factor(COMORBIDITY),
         SEVERITY = as.factor(SEVERITY),
         TOS_CAT1_DESC = as.factor(TOS_CAT1_DESC),
         TOS_CAT2_DESC = as.factor(TOS_CAT2_DESC),
         TOS_CAT3_DESC = as.factor(TOS_CAT3_DESC)
         ) %>%
  #exclude members with complications
  #eclude episodes with episodes number 0 (missing values)
  filter(COMPLICATION == 0, EPISODE_NUM != 0, !is.na(EPISODE_NUM) ) %>%
  #find the time difference between each appointment and the start of th episode
  mutate(lag = difftime( INCURRED_DT_DAY, EPI_START_DT, units = "weeks"))



surgery_test = testdata %>%
  #count the number of episodes for each patient
  group_by(MEM_NUM) %>%
  mutate( nth_episode = length(unique(EPISODE_NUM))) %>%
  ungroup() %>%
  group_by(EPISODE_NUM) %>%
  summarise( 
    surgery_etg = as.factor(any(ETG_FULL  %in% c(712208001, 712208011, 712208111, 712208101))),
    surgery_toscat = as.factor(any((str_detect(CPT_PROC_DESC, "SURG"))|| any(TOS_CAT2_DESC == "SURGICAL"))),
    hit = ifelse(surgery_etg == surgery_toscat, 
                 yes = T, no = F))

confusionMatrix(surgery_test$surgery_etg, surgery_test$surgery_toscat)

testdata %>% 
  filter(ETG_FULL %in% c(712208001, 712208011, 712208111, 712208101)) %>%
  select(HCPCS_PROC_DESC) %>%
  unique()

Validation Set Split

A validation set was created to prevent overfitting the model to the data. A randomly-sampled 20% validation, 80% non-validation split was used. This 80% non-validation set was further split during a 10-fold cross validation which was repeated 3 times.

fraction_Validation <- 0.20
fraction_NonValidation <- 0.80
# Compute sample sizes.
sample_Size_Validation   <- floor(fraction_Validation   * nrow(df))
sample_Size_NonValidation <- floor(fraction_NonValidation * nrow(df))
# Create the randomly-sampled indices for the dataframe. Use setdiff() to
# avoid overlapping subsets of indices.
indices_Validation    <- sort(sample(seq_len(nrow(df)), size = sample_Size_Validation))
indices_Not_Validation <- setdiff(seq_len(nrow(df)), indices_Validation)
# Finally, output the three dataframes for training, validation and test.
dfValidation <- df[indices_Validation, ]
dfNotValidation <- df[indices_Not_Validation,]

Adjusting for Class Inbalance

The surgical rate is only ~6%, which leads to a class inbalance in surgical and non-surgical outcomes. To account for this, the minority population can be oversampled and the majority (surgical) population undersampled. A technique known as Sythetic Minority Oversampling was used from the library SMOTE. To test performance, two versions were used: the first with both over-sampling of the surgical class and undersampling of the non-surgical class, and a second with only over-sampling of the surgical class. For brevity this report only contains the first version, which performed better based on sensitivity and accuracy ratings.

pre_processed_data = dplyr::as_data_frame(dfNotValidation) %>% 
  dplyr::select(-MRI_6_weeks, 
                -xray_6_weeks, 
                -MRI_any, 
                -xray_any, 
                - MEM_NUM,
                - EPISODE_NUM,
                - surgery_lag) %>%
  as.data.frame()

Prior to the oversampling there is about a 95%-5% split between those episodes with surgery and without.

t1 = Sys.time()
rf.data.smote = SMOTE(surgery ~., 
                      data <- pre_processed_data, 
                      perc.over = 1500,
                      perc.under = 100,
                      k = 9)
t2 = Sys.time()
run_time = t1 - t2
run_time

Time difference of -1.45 secs

dim(rf.data.smote)

[1] 23436    13

summary(rf.data.smote$surgery)

 TRUE FALSE 
12096 11340 

After the bootstrap there are roughly an equal number of cases with and without surgery.

Fitting a Model

Because of the number of binary predictors and need to compensate for overfitting, a random forest was used with an automated tuning function from the library caret.

t1 = Sys.time()
# Create model with default paramters
control <- trainControl(method = "repeatedcv", number=10, repeats=3)
seed <- 7
metric <- "Kappa"
set.seed(seed)
mtry <- sqrt(ncol(rf.data.smote))
tunegrid <- expand.grid(.mtry=mtry)
rf_smote <- train(surgery ~., 
                    data= rf.data.smote, 
                    method="rf", 
                    metric=metric, 
                    tuneGrid=tunegrid, 
                    trControl=control)
t2 = Sys.time()
rf_runtime = t2 - t1
rf_runtime

Time difference of 5.044917 mins

Testing the Model

Testing was conducted both against the bootstrapped data and the un-inflated validation set, which contained the raw observations. The cutoff value for the random forest was optimized based on the sensitivity and accuracy of the model.

Test against bootstrapped data

predictor = rf_smote$finalModel$votes[,1]
response = rf.data.smote$surgery 
roc(response, predictor, percent=TRUE,
            # arguments for auc
            plot=TRUE, 
            auc.polygon=TRUE, 
            max.auc.polygon=TRUE, 
            grid=TRUE,
            print.auc=TRUE, 
            show.thres=TRUE)

Call:
roc.default(response = response, predictor = predictor, percent = TRUE,     plot = TRUE, auc.polygon = TRUE, max.auc.polygon = TRUE,     grid = TRUE, print.auc = TRUE, show.thres = TRUE)

Data: predictor in 12096 controls (response TRUE) > 11340 cases (response FALSE).
Area under the curve: 94.19%

The Area Under the Curve (AUC) of 0.942 indicates that the model fits the bootstrapped data well. The interpretation is the probability that a randomly selected subject who ends up with surgery has a test result indicating greater suspicion than that of a randomly chosen subject who does not end up having surgery. (Hanley & McNeil, 1982)

dev.new()
pred = prediction( predictor, response)
perf = performance(pred, measure="tpr",x.measure="cutoff")
plot(perf, 
     main = "True Positive Rate vs Cutoff",
     type = "s")

The true positive rate (sensitivity) is highest when the cutoff is set to close to zero. Because the goal is to correctly identify those patients who have surgery, the cutoff was set close to zero in order to optimize the sensitivity of the model, even at the cost of reduced accuracy. This does not effect the variable importance ranking, however.

dev.new()
pred = prediction( predictor, response)
perf = performance(pred, measure="fpr",x.measure="cutoff")
plot(perf, 
     main = "False Positive Rate vs Cutoff",
     type = "s")

Test Against Validation Set

rf_smote$finalModel$forest$cutoff = c(0.0001, 0.9999)
prediction = predict( object = rf_smote,
                      newdata = dfValidation)
confusionMatrix(prediction, dfValidation$surgery)

Confusion Matrix and Statistics

          Reference
Prediction TRUE FALSE
     TRUE   131  2521
     FALSE   60  3239
                                          
               Accuracy : 0.5663          
                 95% CI : (0.5536, 0.5789)
    No Information Rate : 0.9679          
    P-Value [Acc > NIR] : 1               
                                          
                  Kappa : 0.0343          
 Mcnemar's Test P-Value : <2e-16          
                                          
            Sensitivity : 0.68586         
            Specificity : 0.56233         
         Pos Pred Value : 0.04940         
         Neg Pred Value : 0.98181         
             Prevalence : 0.03210         
         Detection Rate : 0.02201         
   Detection Prevalence : 0.44564         
      Balanced Accuracy : 0.62410         
                                          
       'Positive' Class : TRUE            
                                          

Interpretation

Variable Importance

Logistic Regression of Top Predictors

A logistic regression allows for inferrence on whether or not each predictor increases or decreases the probability of having surgery. A variable with a positive coefficient indicates that this variable increases the likelyhood of the patient having surgery.

The sensitivity of the logistic model is only 42%. The sign of the coefficient is all that matters in this case.

control <- trainControl(method = "repeatedcv", number=10, repeats=3)
logit = train( surgery ~ .,
               data = rf.data.smote, 
               method = "glm",
               family = "binomial", 
               trControl = control)
prediction = predict( object = logit,
                      newdata = dfValidation)
confusionMatrix(prediction, dfValidation$surgery)

Confusion Matrix and Statistics

          Reference
Prediction TRUE FALSE
     TRUE    81  1296
     FALSE  110  4464
                                          
               Accuracy : 0.7637          
                 95% CI : (0.7527, 0.7745)
    No Information Rate : 0.9679          
    P-Value [Acc > NIR] : 1               
                                          
                  Kappa : 0.0497          
 Mcnemar's Test P-Value : <2e-16          
                                          
            Sensitivity : 0.42408         
            Specificity : 0.77500         
         Pos Pred Value : 0.05882         
         Neg Pred Value : 0.97595         
             Prevalence : 0.03210         
         Detection Rate : 0.01361         
   Detection Prevalence : 0.23139         
      Balanced Accuracy : 0.59954         
                                          
       'Positive' Class : TRUE            
                                          

Coefficients

exp(coef(logit$finalModel))

            (Intercept)             surgery_lag                PT_count                  PTTRUE         MRI_6_weeksTRUE spinal_manipulationTRUE      diag_injectionTRUE       xray_6_weeksFALSE 
             0.01892649              1.00467641              0.94601127              0.43084952              0.24220599              1.72237896              0.21415402              7.79866040 
          severity1TRUE          severity2FALSE           severity3TRUE              member_age            nth_episode2            comorbidity1 
             2.77189418              2.70571302              0.36142152              1.01904203              0.44337266              0.82503731