Prediction of diabetes, in an interpretable way

library(tidyverse)
library(Matrix)
library(jsonlite)
library(xgboost)
library(broom)
library(GGally)

md <- src_sqlite("MelbDatathon2017.sqlite", create=FALSE)

# Only use patients in the kaggle training set
max_patient_id <- 279200
patient_ids <- 1:max_patient_id

patients <- tbl(md,"Patients") %>% collect(n=Inf) %>%
    mutate(
        Patient_ID = as.integer(Patient_ID),
        age=ifelse(year_of_birth==1900,NA,2016-year_of_birth)) %>%
    filter(Patient_ID <= max_patient_id)

drug_illness <- tbl(md,"ChronicIllness") %>%
    transmute(illness=ChronicIllness, Drug_ID=as.integer(MasterProductID)) %>%
    collect(n=Inf) %>%
    mutate(illness=ifelse(illness=="Chronic Obstructive Pulmonary Disease (COPD)","COPD",illness))

drugs <- tbl(md,"Drugs") %>% collect(n=Inf) %>%
    rename(Drug_ID=MasterProductID) %>%
    mutate(Drug_ID=as.integer(Drug_ID)) %>%
    left_join(drug_illness, by="Drug_ID")

patient_drug_pre2016 <- tbl(md,"Transactions") %>%
    filter(
        Patient_ID <= max_patient_id,
        DISPENSE_WEEK < "2016-01-01") %>%
    distinct(Patient_ID, Drug_ID) %>%
    collect(n=Inf) %>% ungroup()

patient_illness_pre2016 <- inner_join(patient_drug_pre2016, drug_illness, by="Drug_ID") %>%
    distinct(Patient_ID, illness)


patient_drug_2016 <- tbl(md,"Transactions") %>%
    filter(
        Patient_ID <= max_patient_id,
        DISPENSE_WEEK >= "2016-01-01") %>%
    distinct(Patient_ID, Drug_ID) %>%
    collect(n=Inf) %>% ungroup()

patient_illness_2016 <- inner_join(patient_drug_2016, drug_illness, by="Drug_ID") %>%
    distinct(Patient_ID, illness)

Predicting diabetes in 2016 based on past chronic illness

Logistic regression, the predictor variables are indicator variables for whether a patient has taken a drug associated with a chronic illness in the past.

y <- patient_ids %in% filter(patient_illness_2016, illness=="Diabetes")$Patient_ID
table(y)

y
 FALSE   TRUE 
224854  54346 

illnesses <- sort(unique(drug_illness$illness))

get_illness_col <- function(illness_wanted) {
    result <- cbind(as.numeric(patient_ids %in% filter(patient_illness_pre2016, illness==illness_wanted)$Patient_ID))
    colnames(result) <- illness_wanted
    result
}
x <- data.frame(lapply(illnesses, get_illness_col))

fit <- glm(y ~ ., data=x, family="binomial")
summary(fit)

Call:
glm(formula = y ~ ., family = "binomial", data = x)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9654  -0.2030  -0.1711  -0.1483   3.1787  

Coefficients:
               Estimate Std. Error  z value Pr(>|z|)    
(Intercept)    -4.21676    0.02456 -171.677  < 2e-16 ***
Anti.Coagulant -0.09212    0.03524   -2.614  0.00894 ** 
COPD           -0.12032    0.02937   -4.097 4.18e-05 ***
Depression     -0.17879    0.01965   -9.100  < 2e-16 ***
Diabetes        5.55266    0.02037  272.576  < 2e-16 ***
Epilepsy       -0.12389    0.02612   -4.744 2.10e-06 ***
Heart.Failure  -0.07967    0.01924   -4.141 3.46e-05 ***
Hypertension    0.01431    0.01817    0.788  0.43092    
Immunology     -0.15469    0.12524   -1.235  0.21679    
Lipids          0.42454    0.02219   19.134  < 2e-16 ***
Osteoporosis   -0.36857    0.03929   -9.382  < 2e-16 ***
Urology        -0.07445    0.04452   -1.672  0.09445 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 275232  on 279199  degrees of freedom
Residual deviance:  93115  on 279188  degrees of freedom
AIC: 93139

Number of Fisher Scoring iterations: 7

df <- tidy(fit) %>%
    mutate(conf95_low=estimate-std.error*1.96, conf95_high=estimate+std.error*1.96) %>%
    arrange(term != "(Intercept)", -abs(estimate)) %>%
    select(term, estimate, conf95_low, conf95_high) 
df %>%
    mutate(estimate=round(estimate,3), conf95_low=round(conf95_low,3), conf95_high=round(conf95_high,3))

Viewing these coefficients in a plot, with confidence intervals:

#Confidence interval calculation is very slow
#ggcoef(fit)
ggplot(df,aes(y=factor(term,rev(term)),x=estimate,xmin=conf95_low,xmax=conf95_high)) + 
    geom_point() + geom_errorbarh(height=0) + geom_vline(xintercept=0) + labs(y="",x="contribution to log odds")

Drugs predictive of diabetes in 2016

Using L1 regularization to produce a sparse logistic regression model. glmnet was very slow for this size of model with L1 regularization, so I’ve used XGBoost’s linear mode. The predictor variables are indicator variables for whether a patient has taken a drug in the past.

drug_ids <- sort(unique(patient_drug_pre2016$Drug_ID))

x_drug <- sparseMatrix(
    i=patient_drug_pre2016$Patient_ID,
    j=match(patient_drug_pre2016$Drug_ID, drug_ids),
    x=rep(1, nrow(patient_drug_pre2016)),
    dims=c(max_patient_id, length(drug_ids)))
colnames(x_drug) <- as.character(drug_ids)

# Alpha chosen simply to produce a manageable list of nontrivial coefficients
param <- list(booster="gblinear", eta=0.1, alpha=100)
xfit_drug <- xgboost(data=x_drug, label=y, objective="binary:logistic", param=param, 
    nrounds=300, print_every_n=50)

[1] train-error:0.133983 
[51]    train-error:0.076100 
[101]   train-error:0.076411 
[151]   train-error:0.076468 
[201]   train-error:0.076468 
[251]   train-error:0.076472 
[300]   train-error:0.076479 

dump <- xgb.dump(xfit_drug, dump_format="json") %>% fromJSON
df <- data_frame(
        Drug_ID = as.integer(colnames(x_drug)),
        coef = dump$weight[[1]]) %>%
    left_join(drugs, by="Drug_ID") %>%
    arrange(-abs(coef)) %>%
    filter(abs(coef) >= 0.01) %>%
    select(coef, illness, MasterProductFullName, GenericIngredientName)

df %>% arrange(-coef) %>% mutate(coef=round(coef,digits=3))

Unsurprisingly, taking a diabetes drug prior to 2016 predicts continuing to take a diabetes drug in 2016.

Let us examine the non-diabetes drugs in this table:

df %>% filter(is.na(illness) | illness != "Diabetes") %>% arrange(-coef) %>% mutate(coef=round(coef,digits=3))

Proton-pump inhibitors may be interesting (the “-prazole”s):

https://www.ncbi.nlm.nih.gov/pubmed/22886351 etc