FlexRL-RL-vignette

FlexRL is a package proposing flexible probabilistic Record Linkage with diagnostic tools for downstream inference on linked data. The record linkage algorithm StEM() uses a Stochastic Expectation Maximisation approach to combine 2 data sources and outputs the set of records referring to the same entities.

This vignette uses example subsets from the SHIW data, from the NLTCS data, and synthetic data to showcase how FlexRL works, how StEM() adapts to dynamic PIVs, how RL_diagnostics may guide downstream inference.

We use 4 Partially Identifying Variables to link the records: 2 variables such that dynamics = "stable" and 2 variables that may change over time. Among these 2 dynamic PIVs, we will model the changes of one (dynamics = "structured") and we will let the algorithm flexible about the other (dynamics = "flexible").

PIVs_config <- list( V1 = list(dynamics = "stable",
                               bound_mistakes = c(0.10,0.10),
                               fix_mistakes = c(NA,NA)),
                     V2 = list(dynamics = "stable",
                               bound_mistakes = c(0.10,0.10),
                               fix_mistakes = c(NA,NA)),
                     V3 = list(dynamics = "flexible",
                               bound_mistakes = c(NA,NA),
                               fix_mistakes = c(NA,NA)),
                     V4 = list(dynamics = "structured",
                               bound_mistakes = c(NA,NA),
                               fix_mistakes = c(0.02,0.02),
                               cond_hazard_cov = list(cov1 = c("Xe", "Xf"),
                                                      cov2 = c()))
)
PIVs <- names(PIVs_config)
PIVs_stable <- sapply(PIVs_config, function(x) x$dynamics != "structured")
PIVs_type <- list( V1 = FALSE, V2 = FALSE, V3 = FALSE, V4 = TRUE )

With PIVs that evolve over time, when modelling is possible, it is natural to parametrise the probability that underlying true values of the registered data for a pair of linked records coincide using a survival function. It is possible to use an exponential model, whose baseline hazard is constant in time; a Weibull model, whose baseline hazard is monotone in time; a Gompertz model, whose baseline hazard grows or decays exponentially in time; a piecewise-constant model, whose baseline hazard takes one value on each interval defined by given cut points; a custom model, supplied as a survival function S(X, alpha, times) together with a number of parameters and their starting values. These models may assume proportional hazards, accelerated failure time, proportional odds, additive hazards, … If the covariates are categorical, make sure to encode them into dummies and drop one category (an intercept is included automatically). Since only the agreement of the true values at the observed time gap is used, the likelihood is the same for every model, -sum(Hequal \* log(S) + (1 - Hequal) \* log(1 - S)), and the parameters are estimated with stats::nlminb() at each M-step.

n_values  <- c(8, 9, 10, 15)
p_mistake <- list(V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                  V3 = c(0.02, 0.02), V4 = c(0.02, 0.02))
p_missing <- list(V1 = c(0.005, 0.005), V2 = c(0.005, 0.005),
                  V3 = c(0.005, 0.005), V4 = c(0.005, 0.005))
cond_hazard_params <- list(V1 = c(), V2 = c(), V3 = c(), V4 = c(-0.35, -0.50, -0.70))

gen_data <- simulate_data(
  PIVs_config = PIVs_config, n_values = n_values, n_records = c(450, 500), 
  n_links = 400, p_mistake = p_mistake, p_missing = p_missing, 
  cond_hazard_params = cond_hazard_params, enforce_estimability = TRUE,
  model_dynamics = survival_model("exponential")
)

We pre-process the data for the record linkage task: remove the records for which the linking variables are outside of the common support, check that PIVs are not strongly associated with each other (which may degrade record linkage performance), map PIV values to natural numbers and encode missing values as 0. The StEM algorithm denotes the smallest source as A and the largest as B.

prep_data <- prepare_data(
  data1 = gen_data$data1, data2 = gen_data$data2, label1 = "1", label2 = "2", 
  PIVs_config = PIVs_config, same_mistakes = TRUE, uniq_id = "entity_id",
  restrict_support_intersection = TRUE
)
#> '2' is the larger source, saved as B; '1' saved as A.

head(prep_data$encodedA)
#>   V1 V2 V3 V4 change   date   Xe    Xf local_id entity_id source
#> 1  4  2  8  2  FALSE 0.0061 0.72  1.84        1         1      1
#> 2  5  1  2  2  FALSE 0.0071 2.48  0.91        2         2      1
#> 3  6  8  9 10  FALSE 0.0030 0.17  0.72        3         3      1
#> 4  3  8  7  6  FALSE 0.0084 0.98 -0.81        4         4      1
#> 5  4  8  1  5  FALSE 0.0063 0.14  2.14        5         5      1
#> 6  1  4  9 15  FALSE 0.0080 2.57  1.93        6         6      1

head(prep_data$encodedB)
#>   V1 V2 V3 V4 change    date local_id entity_id source
#> 1  4  2  8  2  FALSE 0.00899        1         1      2
#> 2  5  1  2  2  FALSE 0.00268        2         2      2
#> 3  6  8  9 10  FALSE 0.00752        3         3      2
#> 4  3  8  7  6  FALSE 0.00046        4         4      2
#> 5  4  8  1  5  FALSE 0.00455        5         5      2
#> 6  1  4  9 15  FALSE 0.00182        6         6      2

head(prep_data$true_pairs)
#>   1 2
#> 1 1 1
#> 2 2 2
#> 3 3 3
#> 4 4 4
#> 5 5 5
#> 6 6 6

The number of unique values per PIV gives information on their discriminating power:

prep_data$n_values
#> V1 V2 V3 V4 
#>  8  9 10 15

We gauge how hard the record linkage task is using a naive linkage approach: matching records on exact agreement of identifying information.

naive_fit <- naive_linkage(PIVs = PIVs, 
                           encodedA = prep_data$encodedA, 
                           encodedB = prep_data$encodedB)

df_results <- data.frame( matrix(NA, nrow = 11, ncol = 0) )
rownames(df_results) <- c("TP              ",
                          "FP              ",
                          "FN              ",
                          "sensitivity     ",
                          "FDP             ",
                          "hat FDP score   ",
                          "hat FDP score*  ",
                          "hat FDP synth*  ",
                          "max SMD         ",
                          "min support IoU ",
                          "max MMD         ")

true_pairs      <- do.call(paste, c(prep_data$true_pairs, list(sep="_")))
linked_pairs    <- do.call(paste, c(naive_fit, list(sep = "_")))
true_positive   <- length( intersect(linked_pairs, true_pairs) ) 
false_positive  <- length( setdiff(linked_pairs, true_pairs) ) 
false_negative  <- length( setdiff(true_pairs, linked_pairs) )
sensitivity     <- true_positive / (true_positive + false_negative) 
fdp             <- false_positive / (true_positive + false_positive) 
df_results[1:5,"Naive     "] <- c(true_positive,false_positive,false_negative,sensitivity,fdp)

df_results
#>                  Naive     
#> TP                   261.00
#> FP                   121.00
#> FN                   139.00
#> sensitivity            0.65
#> FDP                    0.32
#> hat FDP score            NA
#> hat FDP score*           NA
#> hat FDP synth*           NA
#> max SMD                  NA
#> min support IoU          NA
#> max MMD                  NA

We run record linkage with FlexRL. The data should contain encodedA, encodedB, n_values, PIVs_config, same_mistakes. StEM_iter, StEM_burnin and gibbs_iter, gibbs_burnin are total number of iterations (including burn-in) and iterations to be discarded as burn-in. If music_on the algorithm opens a short tune in the browser when it finishes. You can save_info_iter (save environment information at each iteration) to new_directory. You can set starting values of the algorithm for parameters gamma and phi with gamma0, phiA0, phiB0, and change the number of samples drawn a posteriori to provide a linkage estimate with n_post_sample (set to 1000 as default).

fit <- StEM(data = prep_data)

The missing values and potential mistakes in the registration, the size of the overlapping set of records between the files, the total number of entities in the population, the low discriminative power of PIVs, PIVs potential dynamics across time, PIVs distribution, dependencies among PIVs, are obstacles to the linkage.

delta_result <- fit$Delta
delta_result <- delta_result[delta_result$x>0.5, ]

linked_pairs    <- do.call(paste, c(delta_result[,c("i","j")], list(sep = "_")))
true_positive   <- length( intersect(linked_pairs, true_pairs) ) 
false_positive  <- length( setdiff(linked_pairs, true_pairs) ) 
false_negative  <- length( setdiff(true_pairs, linked_pairs) )
sensitivity     <- true_positive / (true_positive + false_negative) 
fdp             <- false_positive / (true_positive + false_positive)  

df_results[1:5,"FlexRL    "] <- c(true_positive,false_positive,false_negative,sensitivity,fdp)

df_results
#>                  Naive      FlexRL    
#> TP                   261.00     262.00
#> FP                   121.00      44.00
#> FN                   139.00     138.00
#> sensitivity            0.65       0.66
#> FDP                    0.32       0.14
#> hat FDP score            NA         NA
#> hat FDP score*           NA         NA
#> hat FDP synth*           NA         NA
#> max SMD                  NA         NA
#> min support IoU          NA         NA
#> max MMD                  NA         NA

diag <- RL_diagnostics(fit = fit, encodedA = prep_data$encodedA, encodedB = prep_data$encodedB,
                      compare_vars = PIVs, vars_type_cont = PIVs_type, PIVs = PIVs, 
                      true_pairs = prep_data$true_pairs, FDP_estimation = TRUE, 
                      RL_method = "FlexRL", data = prep_data,
                      maxIter4CV = 1, n_repeats = 2)

df_results[6,"FlexRL    "] <- diag$FDP_measures$curve$RL_FDP_score[1]
df_results[7,"FlexRL    "] <- diag$FDP_measures$curve$augmRL_FDP_score[1]
df_results[8,"FlexRL    "] <- diag$FDP_measures$curve$augmRL_FDP_synth[1]

smd_idx <- grep("smd", names(diag$discrepancy_measures))
iou_idx <- grep("iou", names(diag$discrepancy_measures))
mmd_idx <- grep("mmd", names(diag$discrepancy_measures))

df_results[9, "FlexRL    "] <- max(diag$discrepancy_measures[1,smd_idx])
df_results[10,"FlexRL    "] <- min(diag$discrepancy_measures[1,iou_idx])
df_results[11,"FlexRL    "] <- max(diag$discrepancy_measures[1,mmd_idx])

diag # default, equivalent to print(diag)
#> <RL_diagnostics>
#> 
#>   Linked pairs  (score > 0.50):  306
#> 
#>   FDP           (true pairs):    0.144
#>   Sensitivity   (true pairs):    0.655
#> 
#>   FDP score estimate (RL task):            min threshold 0.50: FDP ~ 0.113
#>                                                threshold 0.50: FDP ~ 0.113
#>                                            max threshold 0.99: FDP ~ 0.004
#>   FDP score estimate (augmented RL task):  min threshold 0.50: FDP ~ 0.129
#>                                                threshold 0.50: FDP ~ 0.129
#>                                            max threshold 0.99: FDP ~ 0.004
#>   FDP synth estimate (augmented RL task):  min threshold 0.50: FDP ~ 0.152
#>                                                threshold 0.50: FDP ~ 0.152
#>                                            max threshold 0.99: FDP ~ 0.000
#>   The score-based estimator is valid when the linkage model is well calibrated to the data.
#>   The synthetic-data estimator is valid when the augmented task is equivalent to the
#>   original one, which requires links and non-links to have similar distributions; similar
#>   score-based FDP estimates on the original and augmented tasks support this prerequisite.
#> 
#>   Multivariate MMD  (linked A vs. A):  0.0031
#>   Multivariate MMD  (linked B vs. B):  0.0037
#> 
#>   IoU support V1  (linked A vs. A):  1.0000
#>   IoU support V2  (linked A vs. A):  1.0000
#>   IoU support V3  (linked A vs. A):  0.9000
#>   IoU support V4  (linked A vs. A):  1.0000
#>   IoU support V1  (linked B vs. B):  1.0000
#>   IoU support V2  (linked B vs. B):  1.0000
#>   IoU support V3  (linked B vs. B):  1.0000
#>   IoU support V4  (linked B vs. B):  0.8750
#> 
#>   Agreement V1  (linked A vs. linked B):  0.9935
#>   Agreement V2  (linked A vs. linked B):  0.9902
#>   Agreement V3  (linked A vs. linked B):  0.9967
#>   Agreement V4  (linked A vs. linked B):  0.7974
#>   Agreement V1  (true pairs):  0.9450
#>   Agreement V2  (true pairs):  0.9250
#>   Agreement V3  (true pairs):  0.9525
#>   Agreement V4  (true pairs):  0.7600
#> 
#>   SMD V1 values: {3, 4, 8, ...}  (linked A vs. A):  0.0664, 0.0694, -0.1025, ...
#>   SMD V2 values: {2, 5, 9, ...}  (linked A vs. A):  0.0775, 0.0370, -0.0528, ...
#>   SMD V3 values: {1, 3, 6, ...}  (linked A vs. A):  0.0509, -0.0744, 0.0539, ...
#>   SMD V4                         (linked A vs. A):  -0.0225
#>   SMD V1 values: {1, 7, 8, ...}  (linked B vs. B):  0.0704, -0.0578, -0.0742, ...
#>   SMD V2 values: {2, 4, 7, ...}  (linked B vs. B):  0.0589, 0.0454, -0.0631, ...
#>   SMD V3 values: {0, 3, 6, ...}  (linked B vs. B):  -0.0353, -0.0899, 0.0563, ...
#>   SMD V4                         (linked B vs. B):  -0.0345

print(diag, threshold = 0.75)
#> <RL_diagnostics>
#> 
#>   Linked pairs  (score > 0.75):  242
#> 
#>   FDP           (true pairs):    0.074
#>   Sensitivity   (true pairs):    0.560
#> 
#>   FDP score estimate (RL task):            min threshold 0.50: FDP ~ 0.113
#>                                                threshold 0.75: FDP ~ 0.049
#>                                            max threshold 0.99: FDP ~ 0.004
#>   FDP score estimate (augmented RL task):  min threshold 0.50: FDP ~ 0.129
#>                                                threshold 0.75: FDP ~ 0.057
#>                                            max threshold 0.99: FDP ~ 0.004
#>   FDP synth estimate (augmented RL task):  min threshold 0.50: FDP ~ 0.152
#>                                                threshold 0.75: FDP ~ 0.022
#>                                            max threshold 0.99: FDP ~ 0.000
#>   The score-based estimator is valid when the linkage model is well calibrated to the data.
#>   The synthetic-data estimator is valid when the augmented task is equivalent to the
#>   original one, which requires links and non-links to have similar distributions; similar
#>   score-based FDP estimates on the original and augmented tasks support this prerequisite.
#> 
#>   Multivariate MMD  (linked A vs. A):  0.0058
#>   Multivariate MMD  (linked B vs. B):  0.0065
#> 
#>   IoU support V1  (linked A vs. A):  1.0000
#>   IoU support V2  (linked A vs. A):  1.0000
#>   IoU support V3  (linked A vs. A):  0.9000
#>   IoU support V4  (linked A vs. A):  1.0000
#>   IoU support V1  (linked B vs. B):  1.0000
#>   IoU support V2  (linked B vs. B):  1.0000
#>   IoU support V3  (linked B vs. B):  1.0000
#>   IoU support V4  (linked B vs. B):  1.0000
#> 
#>   Agreement V1  (linked A vs. linked B):  0.9917
#>   Agreement V2  (linked A vs. linked B):  0.9876
#>   Agreement V3  (linked A vs. linked B):  0.9959
#>   Agreement V4  (linked A vs. linked B):  0.8636
#>   Agreement V1  (true pairs):  0.9450
#>   Agreement V2  (true pairs):  0.9250
#>   Agreement V3  (true pairs):  0.9525
#>   Agreement V4  (true pairs):  0.7600
#> 
#>   SMD V1 values: {1, 4, 8, ...}  (linked A vs. A):  0.1055, 0.0917, -0.0985, ...
#>   SMD V2 values: {5, 6, 7, ...}  (linked A vs. A):  0.0660, 0.0736, -0.0640, ...
#>   SMD V3 values: {1, 3, 6, ...}  (linked A vs. A):  0.1134, -0.1360, 0.1251, ...
#>   SMD V4                         (linked A vs. A):  -0.0523
#>   SMD V1 values: {1, 2, 8, ...}  (linked B vs. B):  0.1307, 0.1007, -0.0702, ...
#>   SMD V2 values: {1, 6, 7, ...}  (linked B vs. B):  0.0859, 0.0817, -0.0960, ...
#>   SMD V3 values: {1, 3, 6, ...}  (linked B vs. B):  0.0950, -0.1531, 0.1279, ...
#>   SMD V4                         (linked B vs. B):  -0.0491

These summaries provide information on the linked data. We can also look at: the convergence of the StEM algorithm parameters, the linkage scores distribution, the data distributions between linked subset and original data sources, the divergence metrics of the linked data from the data sources, the FDP estimation.

plot(diag, "convergence")

plot(diag, "scores")

plot(diag, "distributions", threshold = 0.75)

plot(diag, "discrepancy")

Several other algorithms for Record Linkage have been developed so far in the literature. Each has its own qualities and flaws. StEM() usually outperforms in scenarios where the PIVs have low discriminative power with few expected registration errors; it is particularly efficient when there is enough information to model the dynamics of some PIV changing across time (such as postal code). In 2026, it is the only package taking account of dynamic PIVs. More can be found on this topic in the simulation setting of the methodology paper, code is provided in FlexRL-experiments. More can be found on the FDP estimation methods in this paper.

With link_with_<pkg>() wrappers we can link records with “multilink”, “fastLink”, “BRL”, “reclin2”, “diyar”, “fedmatch”.

brl_args <- list(flds = PIVs, types = rep("bi",length(PIVs)))
fit_brl <- link_with_BRL(prep_data$encodedA, prep_data$encodedB, brl_args)

delta_result    <- as.data.frame(fit_brl)
delta_result    <- delta_result[delta_result$LinkScore>0.5, ]
linked_pairs    <- do.call(paste, c(delta_result[,c("idxA","idxB")], list(sep = "_")))
true_positive   <- length( intersect(linked_pairs, true_pairs) )
false_positive  <- length( setdiff(linked_pairs, true_pairs) )
false_negative  <- length( setdiff(true_pairs, linked_pairs) )
sensitivity     <- true_positive / (true_positive + false_negative) 
fdp             <- false_positive / (true_positive + false_positive) 

df_results[1:5,"BRL       "] <- c(true_positive,false_positive,false_negative,sensitivity,fdp)

diag_brl <- do.call(RL_diagnostics, c(list(fit_brl, prep_data$encodedA, prep_data$encodedB,
                           PIVs, PIVs_type, PIVs, true_pairs = prep_data$true_pairs,
                           FDP_estimation = TRUE, RL_method = "BRL", maxIter4CV = 1, n_repeats = 2),
                           brl_args))

df_results[6,"BRL       "] <- diag_brl$FDP_measures$curve$RL_FDP_score[1]
df_results[7,"BRL       "] <- diag_brl$FDP_measures$curve$augmRL_FDP_score[1]
df_results[8,"BRL       "] <- diag_brl$FDP_measures$curve$augmRL_FDP_synth[1]

smd_idx <- grep("smd", names(diag_brl$discrepancy_measures))
iou_idx <- grep("iou", names(diag_brl$discrepancy_measures))
mmd_idx <- grep("mmd", names(diag_brl$discrepancy_measures))

df_results[9, "BRL       "] <- max(diag_brl$discrepancy_measures[1,smd_idx])
df_results[10,"BRL       "] <- min(diag_brl$discrepancy_measures[1,iou_idx])
df_results[11,"BRL       "] <- max(diag_brl$discrepancy_measures[1,mmd_idx])

multilink_args <- list(types = rep("bi",length(PIVs)), breaks = list(V1 = NA, V2 = NA, V3 = NA, V4 = NA),
                       duplicates = c(0,0), n_iter = 500, max_cc_size = 50)
fit_multilink <- link_with_multilink( prep_data$encodedA, prep_data$encodedB, multilink_args)

delta_result    <- as.data.frame(fit_multilink[c("idxA", "idxB")])
linked_pairs    <- do.call(paste, c(delta_result[,c("idxA","idxB")], list(sep = "_")))
true_positive   <- length( intersect(linked_pairs, true_pairs) )
false_positive  <- length( setdiff(linked_pairs, true_pairs) )
false_negative  <- length( setdiff(true_pairs, linked_pairs) )
sensitivity     <- true_positive / (true_positive + false_negative) 
fdp             <- false_positive / (true_positive + false_positive) 

df_results[1:5,"multilink "] <- c(true_positive,false_positive,false_negative,sensitivity,fdp)

diag_multilink <- do.call(RL_diagnostics, c(list(fit_multilink, prep_data$encodedA, prep_data$encodedB,
                                PIVs, PIVs_type, PIVs, true_pairs = prep_data$true_pairs,
                                FDP_estimation = TRUE, RL_method = "multilink", maxIter4CV = 1, n_repeats = 2),
                                multilink_args))

df_results[6,"multilink "] <- diag_multilink$FDP_measures$curve$RL_FDP_score[1]
df_results[7,"multilink "] <- diag_multilink$FDP_measures$curve$augmRL_FDP_score[1]
df_results[8,"multilink "] <- diag_multilink$FDP_measures$curve$augmRL_FDP_synth[1]

smd_idx <- grep("smd", names(diag_multilink$discrepancy_measures))
iou_idx <- grep("iou", names(diag_multilink$discrepancy_measures))
mmd_idx <- grep("mmd", names(diag_multilink$discrepancy_measures))

df_results[9, "multilink "] <- max(diag_multilink$discrepancy_measures[1,smd_idx])
df_results[10,"multilink "] <- min(diag_multilink$discrepancy_measures[1,iou_idx])
df_results[11,"multilink "] <- max(diag_multilink$discrepancy_measures[1,mmd_idx])

fastLink_args <- list(varnames = PIVs, tol.em = 1e-06, threshold.match = 0.5, return.all = TRUE, n.cores = 1 )
fit_fastLink <- link_with_fastLink( prep_data$encodedA, prep_data$encodedB, fastLink_args )

delta_result    <- as.data.frame(fit_fastLink)
delta_result    <- delta_result[delta_result$LinkScore>0.5, ]
linked_pairs    <- do.call(paste, c(delta_result[,c("idxA","idxB")], list(sep = "_")))
true_positive   <- length( intersect(linked_pairs, true_pairs) )
false_positive  <- length( setdiff(linked_pairs, true_pairs) )
false_negative  <- length( setdiff(true_pairs, linked_pairs) )
sensitivity     <- true_positive / (true_positive + false_negative) 
fdp             <- false_positive / (true_positive + false_positive) 

df_results[1:5,"fastLink  "] <- c(true_positive,false_positive,false_negative,sensitivity,fdp)

diag_fastLink <- do.call(RL_diagnostics, c(list(fit_fastLink, prep_data$encodedA, prep_data$encodedB,
                                PIVs, PIVs_type, PIVs, true_pairs = prep_data$true_pairs,
                                FDP_estimation = TRUE, RL_method = "fastLink", maxIter4CV = 1, n_repeats = 2),
                                fastLink_args))

df_results[6,"fastLink  "] <- diag_fastLink$FDP_measures$curve$RL_FDP_score[1]
df_results[7,"fastLink  "] <- diag_fastLink$FDP_measures$curve$augmRL_FDP_score[1]
df_results[8,"fastLink  "] <- diag_fastLink$FDP_measures$curve$augmRL_FDP_synth[1]

smd_idx <- grep("smd", names(diag_fastLink$discrepancy_measures))
iou_idx <- grep("iou", names(diag_fastLink$discrepancy_measures))
mmd_idx <- grep("mmd", names(diag_fastLink$discrepancy_measures))

df_results[9, "fastLink  "] <- max(diag_fastLink$discrepancy_measures[1,smd_idx])
df_results[10,"fastLink  "] <- min(diag_fastLink$discrepancy_measures[1,iou_idx])
df_results[11,"fastLink  "] <- max(diag_fastLink$discrepancy_measures[1,mmd_idx])

reclin2_args <- list(on = PIVs, formula = stats::as.formula(paste("~", paste(PIVs, collapse = "+"))),
                     type = "mpost", add = TRUE, variable = "selected", score = "mpost", threshold = 0.5)
fit_reclin2 <- link_with_reclin2( prep_data$encodedA, prep_data$encodedB, reclin2_args )

delta_result    <- as.data.frame(fit_reclin2)
delta_result    <- delta_result[delta_result$LinkScore>0.5, ]
linked_pairs    <- do.call(paste, c(delta_result[,c("idxA","idxB")], list(sep = "_")))
true_positive   <- length( intersect(linked_pairs, true_pairs) )
false_positive  <- length( setdiff(linked_pairs, true_pairs) )
false_negative  <- length( setdiff(true_pairs, linked_pairs) )
sensitivity     <- true_positive / (true_positive + false_negative) 
fdp             <- false_positive / (true_positive + false_positive) 

df_results[1:5,"reclin2   "] <- c(true_positive,false_positive,false_negative,sensitivity,fdp)

diag_reclin2 <- do.call(RL_diagnostics, c(list(fit_reclin2, prep_data$encodedA, prep_data$encodedB, 
                               PIVs, PIVs_type, PIVs, true_pairs = prep_data$true_pairs, 
                               FDP_estimation = TRUE, RL_method = "reclin2", maxIter4CV = 1, n_repeats = 2),
                               reclin2_args))

df_results[6,"reclin2   "] <- diag_reclin2$FDP_measures$curve$RL_FDP_score[1]
df_results[7,"reclin2   "] <- diag_reclin2$FDP_measures$curve$augmRL_FDP_score[1]
df_results[8,"reclin2   "] <- diag_reclin2$FDP_measures$curve$augmRL_FDP_synth[1]

smd_idx <- grep("smd", names(diag_reclin2$discrepancy_measures))
iou_idx <- grep("iou", names(diag_reclin2$discrepancy_measures))
mmd_idx <- grep("mmd", names(diag_reclin2$discrepancy_measures))

df_results[9, "reclin2   "] <- max(diag_reclin2$discrepancy_measures[1,smd_idx])
df_results[10,"reclin2   "] <- min(diag_reclin2$discrepancy_measures[1,iou_idx])
df_results[11,"reclin2   "] <- max(diag_reclin2$discrepancy_measures[1,mmd_idx])

diyar_args <- list(probabilistic = TRUE, return_weights = TRUE)
fit_diyar <- link_with_diyar(prep_data$encodedA[, PIVs], prep_data$encodedB[, PIVs], diyar_args)

delta_result    <- as.data.frame(fit_diyar[c("idxA", "idxB")])
linked_pairs    <- do.call(paste, c(delta_result[,c("idxA","idxB")], list(sep = "_")))
true_positive   <- length( intersect(linked_pairs, true_pairs) )
false_positive  <- length( setdiff(linked_pairs, true_pairs) )
false_negative  <- length( setdiff(true_pairs, linked_pairs) )
sensitivity     <- true_positive / (true_positive + false_negative) 
fdp             <- false_positive / (true_positive + false_positive) 

df_results[1:5,"diyar     "] <- c(true_positive,false_positive,false_negative,sensitivity,fdp)

names4diag <- c(PIVs, "local_id", "source")
diag_diyar <- do.call(RL_diagnostics, c(list(fit_diyar, 
                             prep_data$encodedA[, names4diag], prep_data$encodedB[, names4diag], 
                             PIVs, PIVs_type, PIVs, true_pairs = prep_data$true_pairs, 
                             FDP_estimation = TRUE, RL_method = "diyar", maxIter4CV = 1, n_repeats = 2),
                             diyar_args))

df_results[6,"diyar     "] <- diag_diyar$FDP_measures$curve$RL_FDP_score[1]
df_results[7,"diyar     "] <- diag_diyar$FDP_measures$curve$augmRL_FDP_score[1]
df_results[8,"diyar     "] <- diag_diyar$FDP_measures$curve$augmRL_FDP_synth[1]

smd_idx <- grep("smd", names(diag_diyar$discrepancy_measures))
iou_idx <- grep("iou", names(diag_diyar$discrepancy_measures))
mmd_idx <- grep("mmd", names(diag_diyar$discrepancy_measures))

df_results[9, "diyar     "] <- max(diag_diyar$discrepancy_measures[1,smd_idx])
df_results[10,"diyar     "] <- min(diag_diyar$discrepancy_measures[1,iou_idx])
df_results[11,"diyar     "] <- max(diag_diyar$discrepancy_measures[1,mmd_idx])

fedmatch_args <- list(by = PIVs, match_type = "multivar", 
                      unique_key_1 = "unique_key_A", unique_key_2 = "unique_key_B",
                      multivar_settings = fedmatch::build_multivar_settings(
                        compare_type = rep("indicator", length(PIVs)),
                        wgts = rep(1/length(PIVs), length(PIVs))))
fit_fedmatch <- link_with_fedmatch( prep_data$encodedA, prep_data$encodedB, fedmatch_args )

delta_result    <- as.data.frame(fit_fedmatch)
delta_result    <- delta_result[delta_result$LinkScore>0.5, ]
linked_pairs    <- do.call(paste, c(delta_result[,c("idxA","idxB")], list(sep = "_")))
true_positive   <- length( intersect(linked_pairs, true_pairs) )
false_positive  <- length( setdiff(linked_pairs, true_pairs) )
false_negative  <- length( setdiff(true_pairs, linked_pairs) )
sensitivity     <- true_positive / (true_positive + false_negative) 
fdp             <- false_positive / (true_positive + false_positive) 

df_results[1:5,"fedmatch  "] <- c(true_positive,false_positive,false_negative,sensitivity,fdp)

diag_fedmatch <- do.call(RL_diagnostics, c(list(fit_fedmatch, prep_data$encodedA, prep_data$encodedB, 
                               PIVs, PIVs_type, PIVs, true_pairs = prep_data$true_pairs, 
                               FDP_estimation = TRUE, RL_method = "fedmatch", maxIter4CV = 1, n_repeats = 2),
                               fedmatch_args))

df_results[6,"fedmatch  "] <- diag_fedmatch$FDP_measures$curve$RL_FDP_score[1]
df_results[7,"fedmatch  "] <- diag_fedmatch$FDP_measures$curve$augmRL_FDP_score[1]
df_results[8,"fedmatch  "] <- diag_fedmatch$FDP_measures$curve$augmRL_FDP_synth[1]

smd_idx <- grep("smd", names(diag_fedmatch$discrepancy_measures))
iou_idx <- grep("iou", names(diag_fedmatch$discrepancy_measures))
mmd_idx <- grep("mmd", names(diag_fedmatch$discrepancy_measures))

df_results[9, "fedmatch  "] <- max(diag_fedmatch$discrepancy_measures[1,smd_idx])
df_results[10,"fedmatch  "] <- min(diag_fedmatch$discrepancy_measures[1,iou_idx])
df_results[11,"fedmatch  "] <- max(diag_fedmatch$discrepancy_measures[1,mmd_idx])

cat("Record Linkage ")
#> Record Linkage
round(df_results[1:5,],2)
#>                  Naive      FlexRL     BRL        multilink  fastLink  
#> TP                   261.00     262.00     208.00     216.00     216.00
#> FP                   121.00      44.00      21.00      40.00      65.00
#> FN                   139.00     138.00     192.00     184.00     184.00
#> sensitivity            0.65       0.66       0.52       0.54       0.54
#> FDP                    0.32       0.14       0.09       0.16       0.23
#>                  reclin2    diyar      fedmatch  
#> TP                   250.00       82.0     268.00
#> FP                   102.00        0.0     173.00
#> FN                   150.00      318.0     132.00
#> sensitivity            0.62        0.2       0.67
#> FDP                    0.29        0.0       0.39

cat("FDP estimation ")
#> FDP estimation
round(df_results[6:8,],2)
#>                  Naive      FlexRL     BRL        multilink  fastLink  
#> hat FDP score            NA       0.11       0.11        NaN       0.23
#> hat FDP score*           NA       0.13       0.11        NaN       0.24
#> hat FDP synth*           NA       0.15       0.10       0.27       0.24
#>                  reclin2    diyar      fedmatch  
#> hat FDP score           0.3        NaN       0.09
#> hat FDP score*          0.3        NaN       0.09
#> hat FDP synth*          0.3       0.07        NaN

cat("Linked divergence ")
#> Linked divergence
round(df_results[9:11,],2)
#>                  Naive      FlexRL     BRL        multilink  fastLink  
#> max SMD                  NA       0.08       0.12       0.14       0.20
#> min support IoU          NA       0.88       0.88       0.88       0.88
#> max MMD                  NA       0.00       0.01       0.00       0.01
#>                  reclin2    diyar      fedmatch  
#> max SMD                0.28       0.38       0.17
#> min support IoU        0.88       0.88       0.88
#> max MMD                0.02       0.06       0.00