Package {FlexRL}

Title:	Flexible Record Linkage and Linked Data Quality Assessment
Version:	1.0.0
Description:	Probabilistically link records that refer to the same entities across two data sources without a unique identifier, using partially identifying variables such as product code, brand, category, birth year, sex or postal code. 'FlexRL' implements a Stochastic Expectation Maximisation (StEM) approach to Record Linkage (Robach et al., 2025, <doi:10.1093/jrsssc/qlaf016>). The model accounts for registration errors (missing values and mistakes) and for variables that change over time, enforces one-to-one assignment, and has a low memory footprint. The package also provides tools for inference on linked data: two estimators of the false discovery proportion of a linkage (Robach et al., 2025, <doi:10.1002/sim.70292>), based on linkage scores and on synthetic data, and diagnostics comparing the linked sample with the source data. These tools also apply on the linkage output of other record linkage packages.
License:	GPL (≥ 3)
Encoding:	UTF-8
Depends:	R (≥ 4.1.0)
Imports:	Matrix (≥ 1.7), Rcpp (≥ 1.0.13), cli, graphics, grDevices, stats, utils, BRL, arf, diyar, fastLink, fedmatch, mice, multilink, reclin2, synthpop
LinkingTo:	Rcpp
Suggests:	knitr, rmarkdown
VignetteBuilder:	knitr
URL:	https://github.com/robachowyk/FlexRL, https://robachowyk.github.io/FlexRL/
BugReports:	https://github.com/robachowyk/FlexRL/issues
Config/roxygen2/version:	8.1.0
NeedsCompilation:	yes
Packaged:	2026-10-07 14:20:39 UTC; kayane
Author:	Kayané ROBACH [aut, cre, cph], Michel H. HOF [aut, cph]
Maintainer:	Kayané ROBACH <k.c.robach@amsterdamumc.nl>
Repository:	CRAN
Date/Publication:	2026-10-07 15:00:09 UTC

FlexRL: A Flexible Model for Record Linkage

Description

Links records that refer to the same entities across two data sources without a unique identifier, using partially identifying variables. These are stable (not changing over time, probability of mistakes could be bounded), flexible (dynamic but no information to model changes over time) or structured (dynamic and information to model changes over time, probability of mistakes could be fixed). The main function StEM() fits a latent-variable model by stochastic expectation-maximisation; it models registration errors (missing values and mistakes) and changes over time. prepare_data() prepares the data sources for record linkage, RL_diagnostics() gathers diagnostics (FDP estimation and discrepancy metrics for inference on the linked data.

Details

Methodological paper: doi:10.1093/jrsssc/qlaf016. False discovery proportion estimation: doi:10.1002/sim.70292. Experiments repository: https://github.com/robachowyk/FlexRL-experiments.

Author(s)

Kayané Robach

See Also

Useful links:

• https://github.com/robachowyk/FlexRL

• https://robachowyk.github.io/FlexRL/

• Report bugs at https://github.com/robachowyk/FlexRL/issues

Examples

# Link two simulated sources with 4 PIVs: two stable, one flexible, one 
# structured. The true links are known, so performance can be computed.
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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                      cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                  V3 = c(0.05, 0.05), 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 = log(c(0.7, 0.6, 0.5)))
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params, 
                           TRUE, survival_model("exponential") )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
fit <- StEM( data = prep_data, StEM_iter = 10, StEM_burnin = 5,
             gibbs_iter = 10, gibbs_burnin = 5, n_post_sample = 10 )

# linked pairs and performance against the true pairs
linked <- fit$Delta[fit$Delta$x > 0.5, ]
linked_pairs <- paste(linked$i, linked$j, sep = "_")
true_pairs   <- paste(prep_data$true_pairs[[1]], prep_data$true_pairs[[2]], sep = "_")
tp <- length(intersect(linked_pairs, true_pairs))
fp <- length(setdiff(linked_pairs, true_pairs))
fn <- length(setdiff(true_pairs, linked_pairs))
c(LinkageDecisionRule = 0.5, FDP = fp / (tp + fp), Sensitivity = tp / (tp + fn))

# diagnostics for inference on the linked data
diag <- RL_diagnostics(fit, prep_data$encodedA, prep_data$encodedB,
                       names(PIVs_config),
                       list(V1 = FALSE, V2 = FALSE, V3 = FALSE, V4 = TRUE),
                       names(PIVs_config),
                       true_pairs = prep_data$true_pairs, FDP_estimation = TRUE, 
                       RL_method = "FlexRL", data = prep_data,
                       StEM_iter = 5, StEM_burnin = 2, 
                       gibbs_iter = 5, gibbs_burnin = 2, n_post_sample = 5,
                       maxIter4CV = 1, n_repeats = 1)
diag # print(diag)
print(diag, threshold = 0.75)
plot(diag, "scores")
plot(diag, "distributions", threshold = 0.75)
plot(diag, "convergence")
plot(diag, "FDP")
plot(diag, "discrepancy")

Model specific, linkage score based, false discovery proportion at a given threshold

Description

Model specific, linkage score based, false discovery proportion at a given threshold

Usage

FDP_score(LinkScore, threshold)

Arguments

Value

A list with FDP_score Numeric, 1 - mean(score | score > threshold), the estimated FDP at threshold and n_linked Integer, number of linked records.

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                      cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                   V3 = c(0.05, 0.05), 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 = log(c(0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params, TRUE )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
fit <- StEM( data = prep_data, StEM_iter = 10, StEM_burnin = 5,
             gibbs_iter = 10, gibbs_burnin = 5, n_post_sample = 10 )
FDP_score(fit$Delta$x, 0.5)
FDP_score(fit$Delta$x, 0.75)

Model agnostic, based on synthetic data, false discovery proportion at a given threshold

Description

Model agnostic, based on synthetic data, false discovery proportion at a given threshold

Usage

FDP_synth(idxA, idxB, LinkScore, threshold, n_records_A, n_records_B, n_synth)

Arguments

Value

A list with FDP_synth Numeric, proportion of synthetic falsely linked records, n_linked_real Integer, number of real linked records and n_linked_all Integer, total number of linked records (real and synthetic).

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                      cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                   V3 = c(0.05, 0.05), 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 = log(c(0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params, TRUE )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
PIVs <- names(PIVs_config)
n_synth <-as.integer(0.10 * nrow(prep_data$encodedB))
new_data <- synthesise("arf", prep_data$encodedA[, c(PIVs,"local_id","source",
                      "date",prep_data$PIVs_config$V4$cond_hazard_cov$covA)],
                      prep_data$encodedB[, c(PIVs,"local_id","source","date",
                      prep_data$PIVs_config$V4$cond_hazard_cov$covB)], PIVs, 
                      n_synth, TRUE)
new_data$dataA[PIVs][is.na(new_data$dataA[PIVs])] <- 0
new_data$dataB[PIVs][is.na(new_data$dataB[PIVs])] <- 0
# cannot model dynamics for synthetic data, set dates to 0
new_data$dataA$date[is.na(new_data$dataA$date)] <- 0
new_data$dataB$date[is.na(new_data$dataB$date)] <- 0
arguments <- list(data = prep_data, StEM_iter = 5, StEM_burnin = 2, 
                  gibbs_iter = 5, gibbs_burnin = 2, n_post_sample = 5)
fit_flexrl <- link_with_FlexRL(new_data$dataA, new_data$dataB, arguments)  
FDP_synth(fit_flexrl$idxA, fit_flexrl$idxB, fit_flexrl$LinkScore, 0.5,
          nrow(prep_data$encodedA), nrow(prep_data$encodedB), 
          n_synth)

Agreement rate of pairs on shared variables

Description

For a set of pairs, computes how often the two records agree on each variable in common_vars.

Usage

RL_agreement(data1, data2, common_vars, pairs, na.rm = TRUE, na_match = NULL)

Arguments

Value

List with agreements (named numeric vector, one entry per variable in common_vars) and, if available, true_agreements.

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                      cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                   V3 = c(0.05, 0.05), 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 = log(c(0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params, TRUE )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
fit <- StEM( data = prep_data, StEM_iter = 10, StEM_burnin = 5, 
             gibbs_iter = 10, gibbs_burnin = 5, n_post_sample = 10 )                         
linked_pairs <- fit$Delta[fit$Delta$x > 0.5, ]
RL_agreement( prep_data$encodedA, prep_data$encodedB,
              names(PIVs_config), linked_pairs )
RL_agreement( prep_data$encodedA, prep_data$encodedB,
              names(PIVs_config), prep_data$true_pairs )

Post-linkage diagnostics

Description

Gathers, from a fitted linkage, the diagnostics needed to judge the linked data before using it for inference: false discovery proportion estimate, agreement among PIVs of linked pairs, comparison of the linked subset with each source file (standardised mean differences, support overlap, maximum mean discrepancy).

Usage

RL_diagnostics(
  fit,
  encodedA,
  encodedB,
  compare_vars,
  vars_type_cont,
  PIVs,
  true_pairs = NULL,
  FDP_estimation = TRUE,
  ...
)

Arguments

Details

All measures are computed for every threshold from 0.50 to 0.99 (by 0.01) on the linkage scores; for a method without scores (or with a single score value) they are computed once, for all returned pairs. print(x, threshold = ) summarises diagnostics at one threshold, plot() shows the curves.

Value

An object of class "RL_diagnostics", a list with:

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                      cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                   V3 = c(0.05, 0.05), 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 = log(c(0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params, TRUE )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
PIVs <- names(PIVs_config)  
PIVs_type <- list(V1=FALSE, V2=FALSE, V3=FALSE, V4=TRUE)
fit_flexrl <- StEM( data = prep_data, StEM_iter = 5, StEM_burnin = 2,
                    gibbs_iter = 5, gibbs_burnin = 2, n_post_sample = 5 )
diag_flexrl <- RL_diagnostics(fit_flexrl, prep_data$encodedA, prep_data$encodedB,
                              PIVs, PIVs_type, PIVs, true_pairs = prep_data$true_pairs,
                              FDP_estimation = TRUE, RL_method = "FlexRL", 
                              data = prep_data,
                              StEM_iter = 5, StEM_burnin = 2, 
                              gibbs_iter = 5, gibbs_burnin = 2,
                              n_post_sample = 5,
                              maxIter4CV = 1, n_repeats = 1)
diag_flexrl # print(diag_flexrl)
print(diag_flexrl, threshold = 0.75)
plot(diag_flexrl, "scores")
plot(diag_flexrl, "distributions", threshold = 0.75)
plot(diag_flexrl, "convergence")
plot(diag_flexrl, "FDP")
plot(diag_flexrl, "discrepancy")

fit_brl <- link_with_BRL( prep_data$encodedA, prep_data$encodedB, 
                          list( flds = PIVs, 
                                types = rep("bi",length(PIVs)) ) )
diag_brl <- RL_diagnostics(fit_brl, prep_data$encodedA, prep_data$encodedB,
                           PIVs, PIVs_type, PIVs, true_pairs = prep_data$true_pairs,
                           FDP_estimation = TRUE, RL_method = "BRL", 
                           flds = PIVs, types = rep("bi",length(PIVs)),
                           maxIter4CV = 1, n_repeats = 1)
diag_brl # print(diag_brl)
print(diag_brl, threshold = 0.75)
plot(diag_brl, "scores")
plot(diag_brl, "distributions", threshold = 0.75)
plot(diag_brl, "FDP")
plot(diag_brl, "discrepancy")

Stochastic Expectation-Maximisation for record linkage

Description

Fits the FlexRL model with a Stochastic EM algorithm: each iteration runs a Gibbs sampler (alternating between simulating the true PIV values and the linkage matrix D) and then updates the model parameters (gamma, eta, alpha, phi) from the post-burn-in Gibbs draws. See the methodology paper https://10.1093/jrsssc/qlaf016 for details.

Usage

StEM(
  data,
  StEM_iter = 30,
  StEM_burnin = 15,
  gibbs_iter = 20,
  gibbs_burnin = 10,
  music_on = FALSE,
  new_directory = NULL,
  save_info_iter = FALSE,
  gamma0 = NULL,
  phiA0 = NULL,
  phiB0 = NULL,
  n_post_sample = 1000,
  model_dynamics = survival_model("exponential")
)

Arguments

Value

A list with: Delta sparse-matrix summary (i, j, x) of posterior linkage probabilities; a pair is a valid link candidate once x > 0.5 (one-to-one constraint), gamma, eta, alpha, phi the StEM chains for each parameter.

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                      cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                   V3 = c(0.05, 0.05), 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_params1 <- list( V1 = c(), V2 = c(), 
                             V3 = c(), V4 = log(c(0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params1, TRUE )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
fit1 <- StEM( data = prep_data, StEM_iter = 5, StEM_burnin = 2, 
              gibbs_iter = 5, gibbs_burnin = 2, n_post_sample = 5,
              model_dynamics = survival_model("exponential") )
apply(fit1$alpha$V4[2:5,], 2, mean)
cond_hazard_params1$V4
head(fit1$Delta[fit1$Delta$x > 0.5, ])

cond_hazard_params2 <- list( V1 = c(), V2 = c(), 
                             V3 = c(), V4 = log(c(0.3, 0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params2, TRUE,
                           model_dynamics = survival_model("weibull") )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
fit2 <- StEM( data = prep_data, StEM_iter = 5, StEM_burnin = 2, 
              gibbs_iter = 5, gibbs_burnin = 2, n_post_sample = 5,
              model_dynamics = survival_model("weibull"))
apply(fit2$alpha$V4[2:5,], 2, mean)
cond_hazard_params2$V4
head(fit2$Delta[fit2$Delta$x > 0.5, ])

cond_hazard_params3 <- list( V1 = c(), V2 = c(), 
                             V3 = c(), V4 = log(c(0.3, 0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params3, TRUE,
                           model_dynamics = survival_model("gompertz") )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
fit3 <- StEM( data = prep_data, StEM_iter = 5, StEM_burnin = 2, 
              gibbs_iter = 5, gibbs_burnin = 2, n_post_sample = 5,
              model_dynamics = survival_model("gompertz"))
apply(fit3$alpha$V4[2:5,], 2, mean)
cond_hazard_params3$V4
head(fit3$Delta[fit3$Delta$x > 0.5, ])

cond_hazard_params4 <- list( V1 = c(), V2 = c(), 
                             V3 = c(), V4 = log(c(0.6, 0.5, 0.4, 0.3, 0.2, 0.1)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params4, TRUE,
                           model_dynamics = survival_model("piecewise", cuts = c(1,2,3)) )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
fit4 <- StEM( data = prep_data, StEM_iter = 2, StEM_burnin = 1, 
              gibbs_iter = 2, gibbs_burnin = 1, n_post_sample = 2,
              model_dynamics = survival_model("piecewise", cuts = c(1,2,3)))
apply(fit4$alpha$V4, 2, mean)
cond_hazard_params4$V4
head(fit4$Delta[fit4$Delta$x > 0.5, ])

Estimate the false discovery proportion of a record-linkage method via model specific linkage scores

Description

Estimate the false discovery proportion of a record-linkage method via model specific linkage scores

Usage

compute_RL_FDP_score(
  encodedA,
  encodedB,
  PIVs,
  maxIter4CV = 10,
  n_repeats = 10,
  RL_method,
  ...
)

Arguments

Value

List with FDP_score_estimator, Linked_pairs: data frames (n_repeats rows x 50 thresholds, 0.50 to 0.99).

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                      cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                   V3 = c(0.05, 0.05), 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 = log(c(0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params, TRUE )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
PIVs <- names(PIVs_config)                          
compute_RL_FDP_score( prep_data$encodedA, prep_data$encodedB, PIVs, 1, 2,
                      "BRL", flds = PIVs, types = rep("bi",length(PIVs)) )
compute_RL_FDP_score( prep_data$encodedA, prep_data$encodedB, PIVs, 1, 2,
                      "FlexRL", data = prep_data, StEM_iter = 5, 
                      StEM_burnin = 2, gibbs_iter = 5, gibbs_burnin = 2,
                      n_post_sample = 5 )

Estimate the false discovery proportion of a record-linkage method via synthetic augmentation

Description

Repeatedly augments file B with synthetic records (synthesise()), runs the chosen record-linkage method link_with_*, and compares the synthetic-vs-real proportion among linked pairs to estimate the false discovery proportion, for a range of score thresholds (0.50 to 0.99). See https://doi.org/10.1002/sim.70292 for the method.

Usage

compute_augmRL_FDP_synth(
  synth_method,
  encodedA,
  encodedB,
  PIVs,
  n_synth = NULL,
  restrict_support_intersection = TRUE,
  maxIter4CV = 10,
  n_repeats = 10,
  RL_method,
  ...
)

Arguments

Value

List with FDP_score_estimator, FDP_synth_estimator, Linked_pairs_augm, Linked_pairs: data frames (n_repeats rows x 50 thresholds, 0.50 to 0.99).

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                      cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                   V3 = c(0.05, 0.05), 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 = log(c(0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params, TRUE )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
PIVs <- names(PIVs_config)                         
compute_augmRL_FDP_synth( "arf", prep_data$encodedA, prep_data$encodedB, PIVs, 
                          NULL, TRUE, 1, 2, "BRL", flds = PIVs, 
                          types = rep("bi",length(PIVs)) )
compute_augmRL_FDP_synth( "arf", prep_data$encodedA, prep_data$encodedB, PIVs, 
                          NULL, TRUE, 1, 2, "FlexRL", data = prep_data, 
                          StEM_iter = 5, StEM_burnin = 2, 
                          gibbs_iter = 5, gibbs_burnin = 2, n_post_sample = 5 )

Book-keeping data frame for parameters of PIVs dynamics

Description

Internal helper used by StEM() to accumulate, across Gibbs iterations, the covariates, true-value agreement indicator, and time gaps needed to re-estimate the survival (hazard) parameters of a dynamic structured PIV.

Usage

create_data_alpha(n_coef_unstable, stable)

Arguments

Value

An empty data frame with n_coef_unstable + 2 columns (the covariates/intercept, plus Hequal and times) if stable is FALSE; NULL if the PIV is stable (nothing to accumulate).

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                    cov2=c())) )
PIVs_stable <- sapply(PIVs_config, function(x) x$dynamics != "structured")
n_coef_unstable = c(0,0,0,3)
Valpha <- mapply(create_data_alpha, n_coef_unstable = n_coef_unstable,
                 stable = PIVs_stable, SIMPLIFY = FALSE)

Run BRL through the uniform wrapper interface

Description

Links dataA and dataB with the package and returns the linked pairs in the common format used by compute_RL_FDP_score() and compute_augmRL_FDP_synth(). dataA and dataB can be the outputs of synthesise() or the original encoded data sources.

Usage

link_with_BRL(dataA, dataB, arguments, ...)

Arguments

Details

More details about BRL on https://cran.r-project.org/web/packages/BRL/index.html

Value

Named list with idxA, idxB and LinkScore (row indices in dataA, dataB and linkage scores of potential linked records).

Run FlexRL through the uniform wrapper interface

Description

Links dataA and dataB with the package and returns the linked pairs in the common format used by compute_RL_FDP_score() and compute_augmRL_FDP_synth(). dataA and dataB can be the outputs of synthesise() or the original encoded data sources.

Usage

link_with_FlexRL(dataA, dataB, arguments, ...)

Arguments

Details

More details about FlexRL on https://cran.r-project.org/web/packages/FlexRL/index.html

Value

Named list with idxA, idxB and LinkScore (row indices in dataA, dataB and linkage scores of potential linked records).

Run diyar through the uniform wrapper interface

Description

Links dataA and dataB with the package and returns the linked pairs in the common format used by compute_RL_FDP_score() and compute_augmRL_FDP_synth(). dataA and dataB can be the outputs of synthesise() or the original encoded data sources.

Usage

link_with_diyar(dataA, dataB, arguments, ...)

Arguments

Details

More details about diyar on https://cran.r-project.org/web/packages/diyar/index.html

Value

Named list with idxA, idxB and LinkScore = NULL (row indices of default potential linked records in dataA, dataB, this package does not return scores).

Run fastLink through the uniform wrapper interface

Description

Links dataA and dataB with the package and returns the linked pairs in the common format used by compute_RL_FDP_score() and compute_augmRL_FDP_synth(). dataA and dataB can be the outputs of synthesise() or the original encoded data sources.

Usage

link_with_fastLink(dataA, dataB, arguments, ...)

Arguments

Details

More details about fastLink on https://cran.r-project.org/web/packages/fastLink/index.html

Value

Named list with idxA, idxB and LinkScore (row indices in dataA, dataB and linkage scores of potential linked records).

Run fedmatch through the uniform wrapper interface

Description

Links dataA and dataB with the package and returns the linked pairs in the common format used by compute_RL_FDP_score() and compute_augmRL_FDP_synth(). dataA and dataB can be the outputs of synthesise() or the original encoded data sources.

Usage

link_with_fedmatch(dataA, dataB, arguments, ...)

Arguments

Details

More details about fedmatch on https://cran.r-project.org/web/packages/fedmatch/index.html

Value

Named list with idxA, idxB and LinkScore (row indices in dataA, dataB and linkage scores of potential linked records).

Run multilink through the uniform wrapper interface

Description

Links dataA and dataB with the package and returns the linked pairs in the common format used by compute_RL_FDP_score() and compute_augmRL_FDP_synth(). dataA and dataB can be the outputs of synthesise() or the original encoded data sources.

Usage

link_with_multilink(dataA, dataB, arguments, ...)

Arguments

Details

More details about multilink on https://cran.r-project.org/web/packages/multilink/index.html

Value

Named list with idxA, idxB and LinkScore = NULL (row indices of default potential linked records in dataA, dataB, this package does not return scores).

Run reclin2 through the uniform wrapper interface

Description

Links dataA and dataB with the package and returns the linked pairs in the common format used by compute_RL_FDP_score() and compute_augmRL_FDP_synth(). dataA and dataB can be the outputs of synthesise() or the original encoded data sources.

Usage

link_with_reclin2(dataA, dataB, arguments, ...)

Arguments

Details

More details about reclin2 on https://cran.r-project.org/web/packages/reclin2/index.html

Value

Named list with idxA, idxB and LinkScore (row indices in dataA, dataB and linkage scores of potential linked records for the select_threshold() argument).

Log-likelihood of the linkage matrix

Description

Log-likelihood of the linkage matrix

Usage

log_lik(LLL, LLA, LLB, links, sumRowD, sumColD, gamma)

Arguments

Value

Numeric, the log-likelihood of the linkage matrix.

Examples

LLL <- Matrix::Matrix(0, nrow = 13, ncol = 15, sparse = TRUE)
LLA <- stats::runif(13, 0, 2)
LLB <- stats::runif(15, 0, 2)
links <- as.matrix(data.frame(idxA = c(5, 9, 11, 12, 13),
                              idxB = c(5, 9, 11, 13, 15)))
LLL[links] <- 0.67
sumRowD <- (seq_len(13) %in% links[, 1])
sumColD <- (seq_len(15) %in% links[, 2])
gamma <- 0.5
log_lik(LLL, LLA, LLB, links, sumRowD, sumColD, gamma)

Log number of possible linkage configurations

Description

Computes ⁠log(nB! / (nB - sumD)!)⁠, i.e. the log number of ways to choose sumD ordered links among n_records_B records of the larger file; used in log_lik() to normalise the likelihood of the linkage matrix.

Usage

log_possible_config(n_records_B, sumD)

Arguments

Value

Numeric, sum(log((n_records_B - sumD + 1):n_records_B)), or 0 if sumD == 0.

Examples

log_possible_config(n_records_B = 15, sumD = 5)

Maximum Mean Discrepancy between two sets of variables

Description

Joint (multivariate) measure of distributional discrepancy between x and y, using a Gaussian (RBF) kernel with bandwidth set by the median pairwise distance (biased estimator, Gretton et al. 2012).

Usage

mmd(x, y)

Arguments

Value

Numeric, the (biased) MMD estimate.

Examples

x <- matrix(rnorm(50), ncol = 5)
y <- matrix(rnorm(30) + 0.5, ncol = 5)
mmd(x, y)

Naive (deterministic exact-match) record linkage

Description

Links records that agree exactly on every non-missing PIV. Does not enforce the one-to-one assignment constraint, so should be used only to gauge the difficulty of the linkage task (amount of duplication, and discriminative power of the PIVs together).

Usage

naive_linkage(PIVs, encodedA, encodedB, na_match = TRUE, na_is_zero = TRUE)

Arguments

Value

Data frame with columns idxA, idxB: pairs of record indices that match.

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                      cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                   V3 = c(0.05, 0.05), 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 = log(c(0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params, TRUE )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
naive_linkage( names(prep_data$PIVs_config), 
               prep_data$encodedA, prep_data$encodedB )

Plot FDP estimates against the linkage decision threshold

Description

Plot FDP estimates against the linkage decision threshold

Usage

## S3 method for class 'FDP_curves'
plot(x, ...)

Arguments

Value

x, invisibly; called for its plotting side effect.

Plot post-linkage diagnostics

Description

Plot post-linkage diagnostics

Usage

## S3 method for class 'RL_diagnostics'
plot(x, type, threshold = NULL, ...)

Arguments

Value

x, invisibly; called for its plotting side effect.

Plot post-linkage diagnostics over thresholds

Description

Plot post-linkage diagnostics over thresholds

Usage

## S3 method for class 'discrepancy_curves'
plot(x, ...)

Arguments

Value

x, invisibly; called for its plotting side effect.

Monte Carlo convergence (trace) plots for a fitted StEM model

Description

Trace-plots the raw StEM chains of gamma, eta, alpha and phi across iterations, to visually judge whether the chains have stabilised and pick an adequate StEM_burnin for StEM(). One plot per parameter: gamma (proportion linked), eta (PIVs distribution), alpha (hazard coefficients, for dynamic PIVs only), and phi (agreement/missing rates).

Usage

plot_StEM_convergence(fit)

Arguments

Value

NULL, invisibly; called for its plotting side effect.

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                      cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                   V3 = c(0.05, 0.05), 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 = log(c(0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params, TRUE )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
fit <- StEM( data = prep_data, StEM_iter = 10, StEM_burnin = 3,
             gibbs_iter = 10, gibbs_burnin = 3, n_post_sample = 10 )
plot_StEM_convergence(fit)

Compare distributions of shared variables across data sets (histograms/barplots)

Description

Overlays, for each variable in common_vars, the empirical distribution in every data set of data_list (e.g. a baseline file vs. the linked set).

Usage

plot_distributions(data_list, common_vars, threshold, colours = NULL)

Arguments

Value

NULL, invisibly; called for its plotting side effect.

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                      cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                   V3 = c(0.05, 0.05), 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 = log(c(0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params, TRUE )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
fit <- StEM( data = prep_data, StEM_iter = 10, StEM_burnin = 5,
             gibbs_iter = 10, gibbs_burnin = 5, n_post_sample = 10 )
threshold_strict <- stats::quantile(fit$Delta$x, 0.75)
data_list = list(data_baseline = prep_data$encodedA, 
data_select = prep_data$encodedA[fit$Delta[fit$Delta$x > threshold_strict, "i"],])
common_vars = names(PIVs_config)
plot_distributions(data_list, common_vars)
common_vars = c("Xe", "Xf")
plot_distributions(data_list, common_vars)

Plot the distribution of linkage scores

Description

Plot the distribution of linkage scores

Usage

plot_linkage_scores(n_pairs, LinkScore)

Arguments

Value

NULL, invisibly; called for its plotting side effect.

Examples

fit_Delta_x <- c(0,0,0,0,0,0,0,0,0,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,
0.2,0.2,0.2,0.2,0.4,0.4,0.4,0.4,0.4,0.5,0.6,0.7,0.7,0.7,0.7,0.7,0.7,0.8,0.8)
plot_linkage_scores(1000, fit_Delta_x)

Prepare two data sources for StEM()

Description

Wraps the data-preparation steps needed before calling StEM(): tags each source with a source column, labels the larger source as B, adjusting cond_hazard_cov/bound_mistakes/fix_mistakes accordingly if A and B are swapped), drops records whose PIV values fall outside the support shared by both files, warns if two PIVs are strongly associated (Cramer's V > 0.3), which may degrade record linkage performance, encodes every PIV to natural numbers using levels pooled across both sources, and encodes missing values to 0.

Usage

prepare_data(
  data1,
  data2,
  label1,
  label2,
  PIVs_config,
  same_mistakes = TRUE,
  uniq_id = NULL,
  restrict_support_intersection = TRUE
)

Arguments

Value

A list ready to use as the data argument of StEM(): encodedA, encodedB, n_values, PIVs_config, same_mistakes, and true_pairs (NULL if uniq_id was not supplied).

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                      cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                   V3 = c(0.05, 0.05), 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 = log(c(0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params, TRUE )
prep_data <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
str(prep_data, max.level = 1)

Print post-linkage diagnostics

Description

Print post-linkage diagnostics

Usage

## S3 method for class 'RL_diagnostics'
print(x, threshold = 0.5, ...)

Arguments

Value

x, invisibly; called for its printing side effect.

Proportion of a data set with a given variable at a given level

Description

Proportion of a data set with a given variable at a given level

Usage

prop_level(df, var, level)

Arguments

Value

Numeric, proportion of rows where df[[var]] == level (na.rm=TRUE).

Examples

df <- data.frame(colour = sample(c("orange", "purple"), 100, replace = TRUE))
prop_level(df, "colour", "purple")

Simulate the linkage matrix D

Description

Given the current draw of true PIV values, samples the linkage matrix D from its conditional distribution, and returns the updated log-likelihood.

Usage

simulateD(
  data,
  linksR,
  sumRowD,
  sumColD,
  truepivsA,
  truepivsB,
  gamma,
  eta,
  alpha,
  phi,
  model_dynamics = survival_model("exponential")
)

Arguments

Value

A list (Dsample, Rcpp sampleD() output) with: links updated set of links, sumRowD updated sumRowD, sumColD updated sumColD, loglik updated value of the complete log likelihood, nlinkrec updated number of linked records

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                      cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                   V3 = c(0.05, 0.05), 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 = log(c(0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params, TRUE )
data_StEM <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
PIVs_stable <- sapply(data_StEM$PIVs_config, function(x)
                        x$dynamics != "structured")
FlexRL:::initDeltaMap()
linksR = base::matrix(0,0,2)
linksCpp = linksR
sumRowD = rep(0, nrow(data_StEM$encodedA))
sumColD = rep(0, nrow(data_StEM$encodedB))
nlinkrec = 0
survivalpSameH = base::matrix(1, nrow(linksR), length(data_StEM$n_values))
gamma = 0.5
eta = lapply(data_StEM$n_values, function(x) rep(1/x,x))
phi = lapply(data_StEM$n_values, function(x)  c(0.9,0.9,0.1,0.1))
n_coef_unstable = lapply( seq_along(PIVs_stable), function(idx)
 if(PIVs_stable[idx]){ 0 }else{
   ncol(data_StEM$encodedA[, data_StEM$PIVs_config[[idx]]$cond_hazard_cov$covA,
                             drop=FALSE]) +
   ncol(data_StEM$encodedB[, data_StEM$PIVs_config[[idx]]$cond_hazard_cov$covB,
                             drop=FALSE]) + 1 } )
alpha = lapply( seq_along(PIVs_stable),
                function(idx) if(PIVs_stable[idx]){ c(-Inf) }
                else{ rep(log(0.05), n_coef_unstable[[idx]]) }
              )
newTruePivs = simulateH(data=data_StEM, links=linksCpp,
                        survivalpSameH=survivalpSameH,
                        sumRowD=sumRowD, sumColD=sumColD, eta=eta, phi=phi)
truepivsA = newTruePivs$truepivsA
truepivsB = newTruePivs$truepivsB
Dsample = simulateD(data=data_StEM, linksR=linksR, sumRowD=sumRowD,
                   sumColD=sumColD, truepivsA=truepivsA, truepivsB=truepivsB,
                   gamma=gamma, eta=eta, alpha=alpha, phi=phi)
linksCpp = Dsample$links
linksR = linksCpp + 1

Simulate the true PIV values underlying the registered records

Description

Draw the latent true values of each PIV given the currently registered (possibly mistaken or missing) values, the current linkage status, and the current parameters.

Usage

simulateH(data, links, survivalpSameH, sumRowD, sumColD, eta, phi)

Arguments

Value

List with truepivsA and truepivsB, matrices (same shape as the input data) of simulated true PIV values.

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                    cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                  V3 = c(0.05, 0.05), 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 = log(c(0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params, TRUE )
data_StEM <- prepare_data( gen_data$data1, gen_data$data2, "1", "2",
                           PIVs_config, TRUE, "entity_id", TRUE )
PIVs_stable <- sapply(data_StEM$PIVs_config, function(x)
                        x$dynamics != "structured")
FlexRL:::initDeltaMap()
linksR = base::matrix(0,0,2)
linksCpp = linksR
sumRowD = rep(0, nrow(data_StEM$encodedA))
sumColD = rep(0, nrow(data_StEM$encodedB))
nlinkrec = 0
survivalpSameH = base::matrix(1, nrow(linksR), length(data_StEM$n_values))
gamma = 0.5
eta = lapply(data_StEM$n_values, function(x) rep(1/x,x))
phi = lapply(data_StEM$n_values, function(x)  c(0.9,0.9,0.1,0.1))
n_coef_unstable = lapply( seq_along(PIVs_stable), function(idx)
 if(PIVs_stable[idx]){ 0 }else{
   ncol(data_StEM$encodedA[, data_StEM$PIVs_config[[idx]]$cond_hazard_cov$covA,
                             drop=FALSE]) +
   ncol(data_StEM$encodedB[, data_StEM$PIVs_config[[idx]]$cond_hazard_cov$covB,
                             drop=FALSE]) + 1 } )
alpha = lapply( seq_along(PIVs_stable),
                function(idx) if(PIVs_stable[idx]){ c(-Inf) }
                else{ rep(log(0.05), n_coef_unstable[[idx]]) }
              )
newTruePivs = simulateH(data=data_StEM, links=linksCpp,
                        survivalpSameH=survivalpSameH,
                        sumRowD=sumRowD, sumColD=sumColD, eta=eta, phi=phi)
truepivsA = newTruePivs$truepivsA
truepivsB = newTruePivs$truepivsB

Simulate two linked data sources for record linkage benchmarking

Description

Creates two synthetic data sources of given sizes sharing a given number of common entities ("links"), each described by a set of Partially Identifying Variables (PIVs). For every PIV, choose the number of possible values, the proportion of mistakes and of missing values, and whether it is stable over time, flexible (may change but change is not modelled), or structured (expected to change over time, with a survival model for the hazard of change). For structured PIVs, enforce_estimability forces half of the linked pairs to have a near-zero time gap, which helps separate "mistake" from "change over time" when fitting the model.

Usage

simulate_data(
  PIVs_config,
  n_values,
  n_records,
  n_links,
  p_mistake,
  p_missing,
  cond_hazard_params,
  enforce_estimability,
  model_dynamics = survival_model("exponential")
)

Arguments

Value

A list with: data1, data2 the two simulated (encoded) data frames, n_values number of unique values per PIV, time_difference time gap between linked records (NA if no PIV is structured), proba_same_H matrix (n links x n PIVs) of probabilities that the true values coincide, true_pairs data frame with the true 1/2 indices of the linked records

Examples

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.03,0.03),
                               cond_hazard_cov = list(cov1=c("Xe", "Xf"),
                                                    cov2=c())) )
n_values  <- c( 5, 6, 7, 12 )
p_mistake <- list( V1 = c(0.02, 0.02), V2 = c(0.02, 0.02),
                  V3 = c(0.05, 0.05), 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 = log(c(0.7, 0.6, 0.5)) )
gen_data <- simulate_data( PIVs_config, n_values, c(150, 200), 100, 
                           p_mistake, p_missing, cond_hazard_params, TRUE )
str(gen_data, max.level = 1)

Standardised mean difference between a selected set and a baseline

Description

When missing values are encoded as 0m the smd will report information on the missingness in the linked sample vs. the source.

Usage

smd(data_select, data_baseline, var, continuous = TRUE)

Arguments

Value

Named list, one SMD per variable (continuous = TRUE) or per level (continuous = FALSE).

Examples

base <- data.frame(age = rnorm(200, 40, 10), sex = sample(c("M", "F"),
                    200, replace = TRUE))
select <- data.frame(age = rnorm(80, 43, 10), sex = sample(c("M", "F"), 80,
                    replace = TRUE, prob = c(0.6, 0.4)))
smd(select, base, "age")
smd(select, base, "sex", continuous = FALSE)

Intersection-over-union of two histogram supports

Description

Intersection-over-union of two histogram supports

Usage

support_iou(h1, h2)

Arguments

Value

Numeric, IoU of the ranges over which h1 and h2 have non-zero counts.

Examples

h1 <- graphics::hist(rnorm(200), plot = FALSE)
h2 <- graphics::hist(rnorm(200) + 1, plot = FALSE)
support_iou(h1, h2)

Survival models for the dynamics of a PIV

Description

A structured PIV may change between the two registrations. The probability that the true value of a linked pair is unchanged after a time gap t is modelled by a survival function S(t | X, alpha), where X holds an intercept and the covariates given in cond_hazard_cov and alpha the parameters estimated by StEM(). This function builds the model object used by StEM(), simulateD() and simulate_data().

Usage

survival_model(
  type = c("exponential", "weibull", "gompertz", "piecewise", "custom"),
  cuts = NULL,
  S = NULL,
  n_par = NULL,
  init = NULL
)

Arguments

Details

Available models (h is the hazard, ⁠lambda = exp(X alpha_cov)⁠ the proportional-hazards term):

"exponential"

⁠S(t) = exp(-lambda t)⁠; alpha = coefficients of X (default, the model of the methodology paper).

"weibull"

⁠S(t) = exp(-(lambda t)^k)⁠ with shape k = exp(alpha[1]); alpha = log-shape, then coefficients of X.

"gompertz"

⁠h(t) = lambda exp(g t)⁠, so ⁠S(t) = exp(-lambda (exp(g t) - 1) / g)⁠ with g = alpha[1]; alpha = g, then coefficients of X.

"piecewise"

piecewise-constant baseline hazard on the intervals defined by cuts, times lambda; alpha = one log-hazard per interval, then coefficients of the covariates (no intercept).

"custom"

S, n_par and init supplied by the user.

For every model the negative log-likelihood used in the M-step is ⁠-sum(Hequal log S + (1 - Hequal) log(1 - S))⁠, where Hequal indicates that the true values of the linked pair agree. Parameters are estimated with stats::nlminb() (numerical gradient).

Value

An object of class "survival_model": a list with type, S, negloglik, n_par, init and par_names.

Examples

X <- cbind(intercept = rep(1, 5))
times <- c(0.001, 0.2, 1.3, 1.5, 2)
Hequal <- c(TRUE, TRUE, TRUE, FALSE, FALSE)

expo <- survival_model("exponential")
expo$S(X, alpha = log(0.3), times)
stats::nlminb(expo$init(ncol(X)), expo$negloglik, X = X, times = times,
              Hequal = Hequal)$par

weib <- survival_model("weibull")
weib$S(X, alpha = c(0.5, -1), times)
stats::nlminb(weib$init(ncol(X)), weib$negloglik, X = X, times = times,
              Hequal = Hequal)$par

gomp <- survival_model("gompertz")
gomp$S(X, alpha = c(0.5, -1), times)
stats::nlminb(gomp$init(ncol(X)), gomp$negloglik, X = X, times = times,
              Hequal = Hequal)$par
              
pw <- survival_model("piecewise", cuts = 1)
pw$S(X, alpha = c(-2, -1), times)
stats::nlminb(pw$init(ncol(X)), pw$negloglik, X = X, times = times,
              Hequal = Hequal)$par

# a custom model: log-logistic with scale exp(X alpha[-1]) 
# and shape exp(alpha[1])
loglogistic <- survival_model("custom",
  S = function(X, alpha, times) {
    lambda <- exp(as.matrix(X) %*% alpha[-1])
    1 / ( 1 + (lambda * times)^exp(alpha[1]) )
  },
  n_par = function(n_cov) n_cov + 1,
  init = function(n_cov) c(stats::rnorm(1, 0, 0.1), stats::runif(n_cov, log(0.01), log(1))))
loglogistic$S(X, alpha = c(0, -2), times)
stats::nlminb(loglogistic$init(ncol(X)), loglogistic$negloglik, X = X, times = times,
              Hequal = Hequal)$par

Augment file B with synthetic records

Description

Fits a generative model on encodedB[, PIVs] and draws n_synth new synthetic records from it, appended to encodedB with source = "synthetic"; used by compute_augmRL_FDP_synth() to estimate the false discovery proportion of a record-linkage method without ground truth.

Usage

synthesise(
  method,
  encodedA,
  encodedB,
  PIVs,
  n_synth,
  restrict_support_intersection = TRUE
)

Arguments

Value

List with dataA (unchanged encodedA, support-restricted) and dataB (encodedB plus the synthetic records).