Test sets

Let’s set up the environment.

library(tidyverse)
library(lubridate)
library(fuzzyjoin)

Before we create any test sets, we need to create 1-month subsets of our ED and EMS data. We’re going to use the month of October, 2018 because it’s in the middle of our date range, but more importantly, it has a manageable number of EMS runs for when we do test matching, as we’ll be manually evaluating at least one probabilistic linkage run.

We’re also going to convert the arrival date/time to an integer we can numerically compare: minutes since 1/1/2018.

# first load our data
load("~/Idaho_EMS/ED_cleaned.RData")
load("~/Idaho_EMS/EMS_cleaned_and_crosswalk.RData")

# subset EMS data to 10/18
df.EMS.10118_103118 <- df.EMS %>%
  filter(as.Date(Dates_Arrived_at_Destination_E5_10) > as.Date("2018-09-30") & as.Date(Dates_Arrived_at_Destination_E5_10) < as.Date("2018-11-01")) %>%
  mutate(minutes_since_jan_1_2018 = as.numeric(time_length(as.period(interval("2018-01-01",Dates_Arrived_at_Destination_E5_10)),unit="mins")))

# subset ED data to 10/18
all_ED_data_cleaned.10118_103118 <- all_ED_data_cleaned %>%
  filter(as.Date(C_Visit_Date_Time) > as.Date("2018-09-30") & as.Date(C_Visit_Date_Time) < as.Date("2018-11-01")) %>%
  mutate(minutes_since_jan_1_2018 = as.numeric(time_length(as.period(interval("2018-01-01",C_Visit_Date_Time)),unit="mins")))

Test set #1: fuzzy inner join

This set represents the minimal set which must be matched using probabilistic linkage. If probabilistic linkage does not find these matches, then performance is worse than naïve deterministic linkage.

A fuzzy inner join merges data sets based on a combination of exact and windowed matching. We’ll use exact matching for zip code, gender and hospital, and windowed matching for age and arrival time.

The purpose for this test set is as a sanity check that the probabilistic linkage will find at the very least least the lowest hanging fruit.

test_matches<-df.EMS.10118_103118 %>%
 mutate(age_upper_bound = Age_In_Years+1) %>%
 mutate(age_lower_bound = Age_In_Years-1) %>%
 mutate(minutes_upper_bound = minutes_since_jan_1_2018+120) %>%
 mutate(minutes_lower_bound = minutes_since_jan_1_2018-120) %>%
 fuzzy_inner_join(all_ED_data_cleaned.10118_103118, by = c("C_Biosense_Facility_ID" = "Hospital", "Patient_Gender_E6_11" = "Sex", "Patient_Postal_Code_E6_8" = "ZipCode", "age_upper_bound" = "Age", "age_lower_bound" = "Age", "minutes_upper_bound" = "minutes_since_jan_1_2018", "minutes_lower_bound" = "minutes_since_jan_1_2018"),   match_fun = list(`==`, `==`, `==`, `>=`, `<=`, `>=`, `<=`))

Let’s check the record count for this test set.

test_matches %>% select(Unique_Incident_Number) %>% distinct()
## # A tibble: 2,847 x 1
##    Unique_Incident_Number                                                  
##    <chr>                                                                   
##  1 Aberdeen_Springfield_Fire_District_AFD_18_100_AFD_18_100_80_Female_5_3_…
##  2 Aberdeen_Springfield_Fire_District_AFD_18_107_AFD_18_107_22_Female_12_6…
##  3 Ada_County_City_EMS_System_E16_18_0036105_87_Female_4_8_83714           
##  4 Ada_County_City_EMS_System_E16_18_0036845_85_Female_11_19_83703         
##  5 Ada_County_City_EMS_System_E42_18_0034450_63_Male_11_25_83714           
##  6 Ada_County_City_EMS_System_M13_18_0033245_21_Male_5_5_83706             
##  7 Ada_County_City_EMS_System_M13_18_0033260_56_Female_2_26_83705          
##  8 Ada_County_City_EMS_System_M13_18_0033296_64_Male_8_28_83705            
##  9 Ada_County_City_EMS_System_M13_18_0033365_49_Male_12_18_83706           
## 10 Ada_County_City_EMS_System_M13_18_0033382_40_Male_9_10_83705            
## # … with 2,837 more rows

Test set #2: Manual evaluation of matches

For this test set, we’re going to do a test linkage, then evaluate the matches by hand. We’ll use this test data for evaluating performance of subsequent iterations of our probabilistic linkage, as well as determination of correct cutoffs for posterior probability of matches, and for determining the correct priors for feeding into the expectation maximization algorithm (EM).

Probabilistic matching setup

Next we need to set up the windows for full and partial matches. This is a weighting scheme that allows you to account for matches that may not be in the range you choose, but allow you to set how far you’re willing to go to have matches be considered.

# match and partial match windows for age (years)
c1 <- 1; c2 <- 3
# match and partial match windows for arrival time (minutes)
c3 <- 120; c4 <- 200
# match and partial match windows for zipcode (allowing last digit to be off by one)
c5 <- 1; c6 <- 1

Next, we’ll set up our cutoffs, number of cores, and do a little preprocessing to make zipcodes easier to work with.

nc <- 24 ## number of cores to use (up to)
cut <- 0.0001 ## probability cutoff for returned matches before we do deuplication

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

Next we’ll create a temporary unique IDs for sorting and merging operations later.

df.EMS.10118_103118$ID <- 1:nrow(df.EMS.10118_103118)
all_ED_data_cleaned.10118_103118$ID <- 1:nrow(all_ED_data_cleaned.10118_103118)

Now we’ll set up the types of matching for each variable. FastLink has different match functions depending on how you want to match your variables: exact matching, windowed matching, weighted matching using full and partial matches.

g1 <- gammaNUMCKpar(df.EMS.10118_103118$Age_In_Years, all_ED_data_cleaned.10118_103118$Age, cut.a = c1, cut.p = c2, n.cores = nc)
g2 <- gammaNUMCKpar(df.EMS.10118_103118$minutes_since_jan_1_2018, all_ED_data_cleaned.10118_103118$minutes_since_jan_1_2018, cut.a = c3, cut.p = c4, n.cores = nc)
g3 <- gammaNUMCKpar(df.EMS.10118_103118$numeric_zip, all_ED_data_cleaned.10118_103118$numeric_zip, cut.a = c5, cut.p = c6, n.cores = nc)
g4 <- gammaKpar(df.EMS.10118_103118$C_Biosense_Facility_ID, all_ED_data_cleaned.10118_103118$Hospital, n.cores = nc)
g5 <- gammaKpar(df.EMS.10118_103118$Patient_Gender_E6_11, all_ED_data_cleaned.10118_103118$Sex, gender=T, n.cores = nc)
list.obj <- list(g1, g2, g3, g4, g5)

Next we’ll set up our priors. For this first iteration we’re going to assume all records should match. For the prior weight, for this first iteration we really don’t know if we’re correct in our estimation, so we’ll set this weight to be low.

# calculate number of rows of each dataset
nr1 <- nrow(df.EMS.10118_103118); nr2 <- nrow(all_ED_data_cleaned.10118_103118)

# Prior: Number of matches we expect to get.
PL<-nr1/(nr1*nr2)

# the weight to give the prior. 
WL<-0.2

The next step in the probabilistic linkage is to calculate the EM object for scoring matches.

counts <- tableCounts(list.obj, nobs.a = nr1, nobs.b = nr2, n.cores = nc)
##     Parallelizing calculation using OpenMP. 7 threads out of 48 are used.
resultsEM <- emlinkMARmov(patterns = counts, tol = 1e-04, nobs.a = nr1, nobs.b = nr2, prior.lambda = PL, w.lambda=WL)
EM <- data.frame(resultsEM$patterns.w)
EM$prob_m <- resultsEM$zeta.j
EM <- EM[order(EM[, "weights"]), ] 

# save for subsequent use
save(EM, resultsEM,file="~/Idaho_EMS/EM_iteration_1.RData")

Now we’re going to pull out the actual matches.

list.m <- matchesLink(list.obj, nobs.a = nr1, nobs.b = nr2, em = resultsEM, thresh = cut, n.cores = nc)
##     Parallelizing calculation using OpenMP. 7 threads out of 48 are used.
matches.1 <- df.EMS.10118_103118[list.m$inds.a, ]
matches.2 <- all_ED_data_cleaned.10118_103118[list.m$inds.b, ]

At this point we’ll do a couple of little checks to:

  1. See how many of our EMS runs are in this 1st iteration matched set (either good or bad)
  2. Make sure we pulled in all the EMS runs that were in the fuzzy join set
# quick counts
matches.1 %>% select(Unique_Incident_Number) %>% distinct() %>% summary()
##  Unique_Incident_Number
##  Length:3997           
##  Class :character      
##  Mode  :character
# now do a quick sanity check vs. our fuzzy join set. If there's any rows in this result then we have a big problem
test_matches %>%
 filter(!Unique_Incident_Number %in% matches.1$Unique_Incident_Number) %>%
 select(Unique_Incident_Number) %>%
 distinct()
## # A tibble: 0 x 1
## # … with 1 variable: Unique_Incident_Number <chr>

We see here that 3997 of our EMS runs are in this matched set, and that all of our fuzzy join EMS runs are in the matched set.

Now we’re going to add in the actual match probabilities for the Fellegi-Sunter matching patterns.

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))
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), ]

complete_matches_test<-as_tibble(bind_cols(matches.1,matches.2))
complete_matches_test<-complete_matches_test %>%
 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_test$Match_ID = paste(complete_matches_test$Unique_Incident_Number, complete_matches_test$Visit_ID, sep="")

Now we need to de-duplicate. At this point an EMS run has many ED visit matches, and we want to move forward with only the matches with the highest posterior probability per EMS run.

complete_matches_test<-complete_matches_test %>%
 group_by(Unique_Incident_Number) %>%
 top_n(1, prob_m)

Now we’ll do another sanity check vs. our 1st test set

test_matches$Match_ID = paste(test_matches$Unique_Incident_Number, test_matches$Visit_ID, sep="")

complete_matches_test %>% 
 filter(prob_m >= 0.9) %>%
 nrow()
## [1] 3604
complete_matches_test %>% 
 filter(prob_m >= 0.9) %>% 
 filter(Match_ID %in% test_matches$Match_ID) %>%
 select(Unique_Incident_Number) %>% 
 distinct() %>% 
 summary()
##  Unique_Incident_Number
##  Length:2847           
##  Class :character      
##  Mode  :character

We can see here that all of our fuzzy join matches are in the probabilistic linkage set, with posterior probability > 0.9. We can also see that there are more high probability matches then if we had used fuzzy join, so we can see that probabilistic matching is adding value over naive fuzzy joining.

Now we’ll create the de-identified data set for hand-evaluation.

complete_matches_test_for_matching<-complete_matches_test %>%
 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(numeric_zip-numeric_zip1)<=c6,"Y",""),
 H = ifelse((C_Biosense_Facility_ID==Hospital),"Y","")) %>%
 mutate_at(vars(time_diff), funs(as.integer)) %>%
 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_test_for_matching$ID <- 1:nrow(complete_matches_test_for_matching)

# first we'll filter to get rid of times > 16h - these are so far out of our windows that we can throw those out for evaluation purposes
matching_de_id_complete_with_match_ids <- complete_matches_test_for_matching %>%
 filter(time_diff <= 960)
 
# Next we'll save a key so we can later use the results of our manual evaluation for performance measures.
matching_de_id_complete_key <- matching_de_id_complete_with_match_ids %>%
 data.frame() %>%
 select(ID, Match_ID)

# save for later
save(matching_de_id_complete_key, file="matching_de_id_complete_key.RData")
 
# Next well create the export dataset without our match ids. This dataset is completely de-identified
matching_de_id_complete_no_match_ids <- matching_de_id_complete_with_match_ids %>%
 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)

# now let's save this 
write_tsv(matching_de_id_complete_no_match_ids, path="matching_de_id_complete_no_match_ids.txt")

Hand-evaluation

First, download the de-identified file matching_de_id_complete_no_match_ids.txt to your computer. I suggest opening it in MS Excel. Then, add a column at the end named “call”. Next, look at the combination of match attributes in each row to determine whether it is an actual match or should not be a match. Put a “Y” in the “call” column if you think it is a match. Put a “N” in the “call” column if it is not a match. Reading the chief complaints and diagnoses, in combination with the other fields, will help you to determine this. Field descriptions:

  • time_diff = absolute difference in minutes between arrival time recorded by EMS and the C_Visit_Date_Time field in BioSense
  • age_diff = absolute difference in years between ED and EMS matched records
  • A = “Y” if age_diff within the partial match window
  • T = “Y” if time_diff within the partial match window
  • G = “Y” if gender matches
  • Z = “Y” if zip codes match within partial match window
  • H = “Y” if hospital matches

I found that it took me less than a full day to do the hand matching. for the ~4,000 matches in this set.

After you perform the hand-evaluation, save your results with the “call” column as a text file. In Excel this is:

File->Save as…->File Format(pulldown)Tab delimited Text (.txt)

and use the file name: matching_de_id_complete_no_match_ids_with_call.txt

Next upload to your RStudio directory using the upload button.

We can now load this file and use it for evaluation.

matching_de_id_complete_no_match_ids_with_call<-read_tsv("~/matching_de_id_complete_no_match_ids_with_call.txt",col_names=T)
load(matching_de_id_complete_key)

called_matches<-left_join(matching_de_id_complete_key, matching_de_id_complete_no_match_ids_with_call[,c("ID","call")], by = c("ID" = "ID"), copy = FALSE)

called_TPs<-called_matches %>%
 filter(call=="Y") %>%
 select(Match_ID) %>% distinct() %>% pull()

called_FPs<-called_matches %>%
 filter(call=="N") %>%
 select(Match_ID) %>% distinct() %>% pull()

Next we’re going to check what our actual match rate is. We’ll use this for subsequent runs to give a much better prior to the EM algorithm.

length(called_TPs)
## [1] 3314

For my hand evaluation, I got the number of true positive matches in this dataset to be 3,314. So for subsequent rounds I’m going to use that number to calculate the prior. And now that we are more confident in how many actual matches there are, we can increase the weight of the prior as well.

Probabilistic linkage: second iteration

In reality, there were many iterations to get to this point. That is, the point where we have the prior, and have the match weights and windows properly sorted out. To do that, you basically perform more iterations while varying the windows, and see how well you matching compared to your hand evaluation dataset. We’ll just do one of these here, and it will also be the final iteration to create the EM object for use in the full linkage.

We’re going to start after the setting up of the windows and weights.

# Prior: Number of matches we expect to get.
PL<-(3314/nr1)*nr1/(nr1*nr2)
# the weight to give the prior. 
WL<-0.8

counts <- tableCounts(list.obj, nobs.a = nr1, nobs.b = nr2, n.cores = nc)
##     Parallelizing calculation using OpenMP. 7 threads out of 48 are used.
resultsEM <- emlinkMARmov(patterns = counts, tol = 1e-04, nobs.a = nr1, nobs.b = nr2, prior.lambda = PL, w.lambda=WL)
EM <- data.frame(resultsEM$patterns.w)
EM$prob_m <- resultsEM$zeta.j
EM <- EM[order(EM[, "weights"]), ] 

# save for subsequent use
save(EM, resultsEM,file="~/Idaho_EMS/EM_iteration_2.RData")

list.m <- matchesLink(list.obj, nobs.a = nr1, nobs.b = nr2, em = resultsEM, thresh = cut, n.cores = nc)
##     Parallelizing calculation using OpenMP. 7 threads out of 48 are used.
matches.1 <- df.EMS.10118_103118[list.m$inds.a, ]
matches.2 <- all_ED_data_cleaned.10118_103118[list.m$inds.b, ]

# quick counts
matches.1 %>% select(Unique_Incident_Number) %>% distinct() %>% summary()
##  Unique_Incident_Number
##  Length:3995           
##  Class :character      
##  Mode  :character
# now do a quick sanity check vs. our fuzzy join set. If there's any rows in this result then we have a big problem
test_matches %>%
 filter(!Unique_Incident_Number %in% matches.1$Unique_Incident_Number) %>%
 select(Unique_Incident_Number) %>%
 distinct()
## # A tibble: 0 x 1
## # … with 1 variable: Unique_Incident_Number <chr>
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))
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), ]

complete_matches_test<-as_tibble(bind_cols(matches.1,matches.2))
complete_matches_test<-complete_matches_test %>%
 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_test$Match_ID = paste(complete_matches_test$Unique_Incident_Number, complete_matches_test$Visit_ID, sep="")

complete_matches_test<-complete_matches_test %>%
 group_by(Unique_Incident_Number) %>%
 top_n(1, prob_m)

test_matches$Match_ID = paste(test_matches$Unique_Incident_Number, test_matches$Visit_ID, sep="")

complete_matches_test %>% 
 filter(prob_m >= 0.9) %>%
 nrow()
## [1] 3587
complete_matches_test %>% 
 filter(prob_m >= 0.9) %>% 
 filter(Match_ID %in% test_matches$Match_ID) %>%
 select(Unique_Incident_Number) %>% 
 distinct() %>% 
 summary()
##  Unique_Incident_Number
##  Length:2847           
##  Class :character      
##  Mode  :character

The number of matches with prob_m > 0.9 is slightly less than the first iteration. This is because we’ve properly set the prior, and the algorithm is likely doing a better job of distinguishing between true positives and false positives.

Evaluation for posterior probability cutoff

Let’s do a little evaluation now. First let’s plot a histogram of our prob_m distribution for both our called TPs and our FPs.

TP_counts_by_prob_m<-complete_matches_test %>%
 select(Match_ID, prob_m) %>% distinct() %>%
 filter(Match_ID %in% called_TPs) %>%
 group_by(prob_m) %>%
 summarize(n=n()) %>%
 mutate(call="TP")

FP_counts_by_prob_m<-complete_matches_test %>%
 select(Match_ID, prob_m) %>% distinct() %>%
 filter(Match_ID %in% called_FPs) %>%
 group_by(prob_m) %>%
 summarize(n=n()) %>%
 mutate(call="FP")

eval_counts_by_prob_m<-full_join(TP_counts_by_prob_m, FP_counts_by_prob_m, by="prob_m") %>%
 arrange(prob_m)

ggplot(eval_counts_by_prob_m, aes(x=substr(as.character(prob_m),1,5))) +
 geom_bar(aes(y=n.x, fill="TP"), stat="identity", alpha = 0.5) +
 geom_bar(aes(y=n.y, fill="FP"), stat="identity", alpha = 0.5) +
 labs(fill="call",x="posterior probability of match",y="count")

You can see that the TPs (green) overlap with the FPs (red) slightly, and that below 0.985 there are no TPs. Overall, 0.985 looks like a good candidate for a cutoff to maximize recovery of TPs, while minimizing FPs.

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_test %>% 
 data.frame() %>% select(prob_m) %>% distinct() %>% pull()

match_rates<-tibble()

for (prob in distinct_probs) {
 test_matches_for_match_rate<-complete_matches_test %>%
  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.10118_103118)
  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))

What this plot shows is that as we move the cutoff from 0 to the highest possible prob_m, we filter out more and more matches, presumably of lower quality, leading to a lower overall match rate.

So at our cutoff, 0.985, we’re getting around an 81% match rate.

Evaluation using ROC curve and MCC

We’re going to evaluate our chosen cutoff by how well it balances sensitivity (true positive rate, or TPR) and specificity (related to the false positive rate, or FPR). This type of plot is called a receiver operating curve (ROC).

Following that, we’ll calculate a correlation coefficient for how well our predictions of matching match the actual matching, called Matthews correlation coefficient (MCC).

First we’ll merge our calls with this iteration’s results, then get calculate our test statistics for each prob_m. The stats we’ll calculate:

  • sensitivity, recall, hit rate, or true positive rate (TPR) = TP / (TP+FN)
  • fall-out or false positive rate (FPR) = FP / (FP+TN)
  • precision or positive predictive value (PPV) = TP / (TP+FP)
  • specificity, selectivity or true negative rate (TNR) = TN / (TN+FP)
  • negative predictive value (NPV) = TN / (TN+FN)
  • false discovery rate (FDR) = FP / (FP+TP)
  • miss rate or false negative rate (FNR) = FN / (FN+TP)
  • false omission rate (FOR) = FN / (FN+TN)
if(!require(ggrepel)) install.packages("ggrepel")
if(!require(pracma)) install.packages("pracma")
library(ggrepel)


complete_matches_test <- complete_matches_test %>%
 mutate(Call = ifelse(Match_ID %in% called_TPs,"Y",ifelse(Match_ID %in% called_FPs,"N","")))

ROC_data<-tibble()

for (prob in distinct_probs) {
  temp_row<-complete_matches_test %>% 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)
}

Now we’ll merge in the match rates from before, and plot the ROC including match rates at each cutoff.

ROC_data<-inner_join(ROC_data, match_rates, by = c("prob" = "prob"))

# 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=paste(substr(as.character(prob),1,5),", ",substr(as.character(MR*100),1,2),"%",sep="")),
        size=3, box.padding = unit(1.0, "lines")
    ) + 
    scale_y_continuous(labels=scaleFUN)

The ROC curve shows that at our 0.985 cutoff (with 81% of our input EMS runs matched to an ED record), we’re maximizing true positive matches. If we decreased the cutoff, we would only be adding false positives.

Area under the curve (AUC)

Next we’ll calculate the area under the curve (AUC). This statistic is usually interpreted as the accuracy of a test. These are commonly accepted ranges for classification of a test best on AUC:

  • .90-1 = excellent
  • .80-.90 = good
  • .70-.80 = fair
  • .60-.70 = poor
  • .50-.60 = fail
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.9945443264739587

This result shows that our probabilistic linkage algorithm is a nearly perfect test for distinguishing matches.

Matthews bar graph

Matthews correlation coefficient (MCC) is interpreted exactly the same as a Pearson correlation coefficient:

  • -1 when there is perfect disagreement between actuals and prediction
  • 1 when there is a perfect agreement between actuals and predictions
  • 0 when the prediction may as well be random with respect to the actuals.

We’ve already calculated the stat - we just need to plot it.

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

You can see in this plot that we get the highest correlation between estimated and actual matches at the 0.985 cutoff.

In conclusion - we now have an EM object that performs very well at distinguishing true matches from false matches, and we also have determined what cutoff we should probably use going forward. Now we’re ready to do the full matches.