Full matching

Load EM object and cleaned data

In our test set iteration, we created an EM object (resultsEM) and saved it to a file. Let’s load that now so we can use it for matching the full dataset. We’ll setup the environment as well.

load("EM_iteration_2.RData")
load("EMS_cleaned_and_crosswalk.RData")
load("ED_cleaned.RData")

library(Rcpp)
library(fastLink)
library(tidyverse)
library(lubridate)

We’ll also set up our thresholds for full and partial matches.

c1 <- 1; c2 <- 3
c3 <- 120; c4 <- 200
c5 <- 1; c6 <- 1
nc <- 24 ## number of cores to use (up to)
cut <- 0.0001 ## probability cutoff for returned matches before we do deuplication

Because of the amount of computer memory on the NSSP RStudio server, we’re not going to try to link all the records at once. Instead, we’re going to chunk the records into 12 time periods, each overlapping by one day. We do this to eliminate edge effects. You can set the number of chunks to be higher if you run into problems such as crashes, but you will very likely get crashes if you go below 12.

num_chunks=12 
edge_overlap=1440 # 1 day's worth of minutes

Next, we’ll set up the linkage run. Here, we’ll do the same pre-processing we did for the test iterations, just now it’s on the full dataset.

# make zip a comparable number
df.EMS$numeric_zip<-as.numeric(df.EMS$Patient_Postal_Code_E6_8)
all_ED_data_cleaned$numeric_zip<-as.numeric(all_ED_data_cleaned$ZipCode)

# calc minutes
df.EMS<-df.EMS %>% 
mutate(
  minutes_since_jan_1_2018 = as.numeric(time_length(as.period(interval("2018-01-01",Dates_Arrived_at_Destination_E5_10)),unit="mins"))
)
all_ED_data_cleaned<-all_ED_data_cleaned %>% 
mutate(
  minutes_since_jan_1_2018 = as.numeric(time_length(as.period(interval("2018-01-01",C_Visit_Date_Time)),unit="mins"))
)

# create ID
df.EMS$ID <- 1:nrow(df.EMS)
all_ED_data_cleaned$ID <- 1:nrow(all_ED_data_cleaned)

# first get overlapping dates for both datasets - set min and max dates
date_lower_bound<-max(
 min(df.EMS$minutes_since_jan_1_2018, na.rm=T), 
 min(all_ED_data_cleaned$minutes_since_jan_1_2018, na.rm=T)
)
date_upper_bound<-min(
 max(df.EMS$minutes_since_jan_1_2018, na.rm=T), 
 max(all_ED_data_cleaned$minutes_since_jan_1_2018, na.rm=T)
)

# next divide into chunks 
chunk_bounds<-seq(date_lower_bound, date_upper_bound, length.out = num_chunks+1)
complete_matches<-tibble()
full_match_checking<-tibble()

Now we’ll actually run through all the chunks. For each one, we’ll use our pre-computed EM object to score the matches. As each chunk is processed, it will be merged with the previous chunks and saved, in case there’s any problems. This will help you to figure out where the failure was.

# now loop through each chunk
for (current_chunk in 1:num_chunks) {
    print(paste("merging chunk",current_chunk))
    temp_lower_bound<-chunk_bounds[current_chunk]
    temp_upper_bound<-chunk_bounds[current_chunk+1]
    if (current_chunk<num_chunks) {
        temp_upper_bound <- temp_upper_bound+edge_overlap
    }
    temp.EMS<-df.EMS %>%
     filter(minutes_since_jan_1_2018 >= temp_lower_bound & minutes_since_jan_1_2018 <= temp_upper_bound)
    temp.ED<-all_ED_data_cleaned %>%
     filter(minutes_since_jan_1_2018 >= temp_lower_bound & minutes_since_jan_1_2018 <= temp_upper_bound)

    if (nrow(temp.EMS) <= 1 | nrow(temp.ED) <= 1) { break }

    print(paste("calculating gammas"))

    print(paste("   g1"))
    g1 <- gammaNUMCKpar(temp.EMS$Age_In_Years, temp.ED$Age, cut.a = c1, cut.p = c2, n.cores = nc)
    print(paste("   g2"))
    g2 <- gammaNUMCKpar(temp.EMS$minutes_since_jan_1_2018, temp.ED$minutes_since_jan_1_2018, cut.a = c3, cut.p = c4, n.cores = nc)
    print(paste("   g3"))
    g3 <- gammaNUMCKpar(temp.EMS$numeric_zip, temp.ED$numeric_zip, cut.a = c5, cut.p = c6, n.cores = nc)
    print(paste("   g4"))
    g4 <- gammaKpar(temp.EMS$C_Biosense_Facility_ID, temp.ED$Hospital, n.cores = nc)
    print(paste("   g5"))
    g5 <- gammaKpar(temp.EMS$Patient_Gender_E6_11, temp.ED$Sex, gender=T, n.cores = nc)
    nr1 <- nrow(temp.EMS); nr2 <- nrow(temp.ED)
    list.obj <- list(g1, g2, g3, g4, g5)

    print(paste("Getting matches"))
    list.m <- matchesLink(list.obj, nobs.a = nr1, nobs.b = nr2, em = resultsEM, thresh = cut, n.cores = nc)
    matches.1 <- temp.EMS[list.m$inds.a, ]
    matches.2 <- temp.ED[list.m$inds.b, ]

    # quick count
    print(paste("Merging match results"))
    temp.0 <- abs(matches.1$Age_In_Years - matches.2$Age) 
    gamma.1 <- ifelse(temp.0 <= c1, 2, ifelse(temp.0 <= c2, 1, 0))
    rm(temp.0)

    temp.0 <- abs(matches.1$minutes_since_jan_1_2018 - matches.2$minutes_since_jan_1_2018) 
    gamma.2 <- ifelse(temp.0 <= c3, 2, ifelse(temp.0 <= c4, 1, 0))
    rm(temp.0)

    temp.0 <- abs(matches.1$numeric_zip - matches.2$numeric_zip) 
    gamma.3 <- ifelse(temp.0 <= c5, 2, ifelse(temp.0 <= c6, 1, 0))
    rm(temp.0)

    temp.0 <- (matches.1$C_Biosense_Facility_ID == matches.2$Hospital) 
    gamma.4 <- ifelse(temp.0 == T, 2, 0)
    rm(temp.0)

    temp.0 <- (matches.1$Patient_Gender_E6_11 == matches.2$Sex) 
    gamma.5 <- ifelse(temp.0 == T, 2, 0)
    rm(temp.0)

    gammas <- data.frame(cbind(gamma.1, gamma.2, gamma.3, gamma.4, gamma.5))
    # gammas <- data.frame(cbind(gamma.1, gamma.2, gamma.3, gamma.4))
    matches.1 <- cbind(matches.1, gammas)
    matches.2 <- cbind(matches.2, gammas)

    matches.1$ido <- 1:nrow(matches.1)
    matches.2$ido <- 1:nrow(matches.2)

    matches.1 <- merge(matches.1, EM, by = c('gamma.1', 'gamma.2', 'gamma.3', 'gamma.4', 'gamma.5'), all.x = T)
    matches.2 <- merge(matches.2, EM, by = c('gamma.1', 'gamma.2', 'gamma.3', 'gamma.4', 'gamma.5'), all.x = T)

    matches.1 <- matches.1[order(matches.1$ido), ]
    matches.2 <- matches.2[order(matches.2$ido), ]

    print(paste("deduplicating"))

    complete_matches<-bind_rows(
     complete_matches,
     as_tibble(bind_cols(matches.1,matches.2))
    ) 
    complete_matches<-complete_matches %>%
     select(-c(counts1, gamma.11, gamma.21, gamma.31, gamma.41, gamma.51, p.gamma.j.m1, p.gamma.j.u1, prob_m1, weights1, ido1))
    complete_matches$Match_ID = paste(complete_matches$Unique_Incident_Number, complete_matches$Visit_ID, sep="")

    # dedup
    complete_matches<-complete_matches %>%
     group_by(Unique_Incident_Number) %>%
     top_n(1, prob_m)
         
    # save matches
    save(complete_matches, full_match_checking, file=paste("matches_up_to_chunk_",current_chunk,".RData",sep=""))
}

Now we’ll get rid of duplicate rows from the overlaps on the edges of the chunks, and then save our matches.

# dedup for distinct rows
complete_matches<-distinct(complete_matches)

save(complete_matches, file="complete_matches.RData")

Random sample for hand-checking

We can check our matching results using the same methods we used for the test iterations. First we’ll create a random sample of 12,000 records for hand-checking.

complete_matches_for_matching<-complete_matches %>%
 mutate(time_diff = abs(minutes_since_jan_1_2018 - minutes_since_jan_1_20181),
 age_diff = abs(Age_In_Years - Age),
 A = ifelse(abs(Age_In_Years-Age)<=c2,"Y",""),
 T = ifelse(abs(minutes_since_jan_1_2018-minutes_since_jan_1_20181)<=c4,"Y",""),
 G = ifelse((Patient_Gender_E6_11==Sex),"Y",""),
 Z = ifelse(abs(Patient_Postal_Code_E6_8-ZipCode)<=c6,"Y",""),
 H = ifelse((C_Biosense_Facility_ID==Hospital),"Y","")) %>%
 select(Match_ID, Unique_Incident_Number, time_diff, age_diff, A, T, G, Z, H, Primary_Symptom_E9_13, Drug_Use_Indicators_Code_E12_19, Complaint_Organ_System_E9_12, Chief_Complaint_E9_5, Diagnosis_Combo, Category_flat, SubCategory_flat, Chief_Complaint_Combo, prob_m) %>%
 distinct()
complete_matches_for_matching$ID <- 1:nrow(complete_matches_for_matching)

# filter to get rid of times > 16h and convert to dataframe
matching_de_id_complete <- complete_matches_for_matching %>%
 filter(time_diff <= 960) %>%
 data.frame() %>%
 select(ID, time_diff, age_diff, A, T, G, Z, H, Primary_Symptom_E9_13, Drug_Use_Indicators_Code_E12_19, Complaint_Organ_System_E9_12, Chief_Complaint_E9_5, Diagnosis_Combo, Category_flat, SubCategory_flat, Chief_Complaint_Combo, prob_m)

matching_de_id_sample_12k <- matching_de_id_complete %>%
 sample_n(12000, replace=F)

# now save
save(matching_de_id_complete, complete_matches_for_matching, matching_de_id_sample_12k, file="matching_de_id_complete_full_matches.RData")
write_tsv(matching_de_id_sample_12k, path="matching_de_id_sample_12k.txt")

You can perform the hand-evaluation the same way as for the test iteration. To get good statistics, I would classify at least 4,000 of the random matches. Just to give an idea of the labor involved, C3 informatics classified the full 12,000 in one work day with 3 people doing the job.

Again, open the matching_de_id_sample_12k.txt in MS Excel, create a “call” column, and classify as good matches using Y or N. Save as a TSV text file named “consolidated_matches.txt” and upload to the RStudio server.

Evaluation

Evaluation for match rate

Match rate is the percentage of our input EMS runs that end up in a match that is above our cutoff.

# calculate number of unique EMS runs in matches for each cutoff
distinct_probs<-complete_matches %>% 
 data.frame() %>% select(prob_m) %>% distinct() %>% pull()

match_rates<-tibble()

for (prob in distinct_probs) {
 test_matches_for_match_rate<-complete_matches %>%
  filter(prob_m >= prob) %>%
  select(Unique_Incident_Number) %>%
  distinct() %>%
  nrow()
  temp_row=list()
  temp_row$prob=prob
  temp_row$MR=test_matches_for_match_rate/nrow(df.EMS)
  match_rates<-bind_rows(match_rates, temp_row)
}
 
# now plot
ggplot(data=match_rates, aes(x=substr(as.character(prob),1,7)),y=MR) +
 geom_bar(aes(y=MR, fill="MR"), stat="identity", alpha = 0.6) +
 geom_text(aes(y=MR, label=paste(substr(as.character(MR*100),1,2),"%",sep="")), hjust=1.2, color="white", size=4.0, angle=90) +
 labs(fill="MR",x="posterior probability of match",y="Match rate") + 
 theme(axis.text.x = element_text(angle = 90, hjust = 1))

So we’re still getting an 81% match rate in the full match set.

ROC curve

Now we’ll calculate the ROC curve for the full linkage results.

consolidated_matching<-read_tsv("consolidated_matches.txt",col_names=T)

distinct_probs<-consolidated_matching %>%
 select(prob_m) %>% distinct() %>% pull()

ROC_data<-tibble()
for (prob in distinct_probs) {
  temp_row<-consolidated_matching %>% data.frame() %>%
  summarize(
   prob=prob,
   TP = sum(prob_m>=prob & Call=="Y"),
   FP = sum(prob_m>=prob & Call=="N"),
   TN = sum(prob_m<prob & Call=="N"),
   FN = sum(prob_m<prob & Call=="Y"),
   P_exclusive = sum(prob_m==prob & Call=="Y"),
   N_exclusive = sum(prob_m==prob & Call=="N"),
  ) %>%
  mutate(
   TPR = TP / (TP+FN),
   FPR = FP / (FP+TN),
   PPV = TP / (TP+FP),
   TNR = TN / (TN+FP),
   NPV = TN / (TN+FN),
   FDR = FP / (FP+TP),
   FNR = FN / (FN+TP),
   FOR = FN / (FN+TN),
   MCC = sqrt(PPV*TPR*TNR*NPV) - sqrt(FDR*FNR*FPR*FOR)
  )
  ROC_data<-bind_rows(ROC_data, temp_row)
}

if(!require(ggrepel)) install.packages("ggrepel")
if(!require(zoo)) install.packages("zoo")
library(ggrepel)

# need to round probabilities to merge
# this is just to get 2 significant digits for the y-axis labels
scaleFUN <- function(x) sprintf("%.2f", x)
# this is the plot
ggplot(ROC_data,aes(FPR,TPR)) +
    geom_line(size = 2, alpha = 0.5, color="red3") +
    labs(
        title= "ROC curve", 
        x = "False Positive Rate (1-Specificity)", 
        y = "True Positive Rate (Sensitivity)"
    ) +
    geom_text_repel(
        data = ROC_data,
        mapping=aes(x=FPR, y=TPR,label=substr(as.character(prob),1,5)),
        size=3, box.padding = unit(1.0, "lines")
    ) + 
    scale_y_continuous(labels=scaleFUN)

The ROC curve shows that below the 0.985 prob_m cutoff, we will have a sharp increase in false positives.

library(zoo)

AUC_test<-ROC_data %>%
 select(FPR,TPR) %>%
 distinct()

# calc AUC
id <- order(AUC_test$FPR)
AUC <- sum(diff(AUC_test$FPR[id])*rollmean(AUC_test$TPR[id],2))
AUC

## [1] 0.9952

The AUC shows that this is an almost perfect test for distinguishing true positive matches from false matches.

Matthews correlation coefficient (MCC)

We’ll calculate the MCC just as we did on the test iteration.

# matthews bar graph
ggplot(ROC_data[!is.nan(ROC_data$MCC),],aes(substr(as.character(prob),1,5),MCC)) +
 geom_bar(stat="identity", width = 0.6, fill="tomato2") + 
 geom_text(aes(label=substr(as.character(MCC),1,5)), hjust=1.1, color="white", size=3.0, angle=90)+
 labs(title="Matthews correlation coefficient (MCC)", 
       x = "Posterior probability cutoff", 
       y = "MCC") +
 ylim(0,1) +
 theme(axis.text.x = element_text(angle=65, vjust=0.6))

The plot shows that the 0.985 cutoff still gives an extremely high correlation between estimated and actual matches. The 0.994 cutoff gives a microscopic advantage, but at this cutoff we’d be excluding many true positives.