postlink: An R Package for Post-Linkage Data Analysis
2026-10-21
Source:vignettes/postlink.Rmd
postlink.RmdAbstract
The advent of the digital age has reshaped the processes by which
data are collected, shared, and analyzed. An important trend within this
landscape is the integration of disparate data sources, often collected
independently and stored in isolation. Consequently, methods for data
integration are needed to realize the value of all existing data sources
in combination. Among these methods, Record Linkage (RL) plays a pivotal
role by enabling the connection of individual records across multiple
data sets. RL spans a wide array of contemporary applications, including
the linkage of survey data, administrative and electronic health
records, as well as historical censuses and vital records. Open-source
software to facilitate RL when records are linked using
quasi-identifiers is widely available. However, software for analyzing
linked files that accounts for the uncertainty and error-prone nature of
RL is lacking. Specifically, there is lack of software in the secondary
analysis setting in which RL and downstream analysis are conducted
separately. We describe postlink, an R package dedicated to
analyzing data integrated using error-prone record linkage. A central
feature of the package is its extensible, formula-based architecture,
offering a suite of post-linkage counterparts to R’s standard modeling
toolkit.
1. Introduction
Information integration across multiple data sources is increasingly being applied in various scientific and policy contexts. Record linkage (RL) and privacy preserving record linkage (PPRL) refer to a set of tools that identify records on the same entity across data sources. In the absence of unique identifying information, RL and PPRL are prone to errors. Records from different entities may be incorrectly linked (false matches), and records that correspond to the same entity may be missed (false non-matches). These errors can result in biased and inaccurate estimates (Neter et al., 1965). To obtain statistically valid estimates, analyses of linked datasets should adjust for these errors (Chambers & Diniz da Silva, 2020; Kamat & Gutman, 2026; Wang et al., 2022). However, because of privacy considerations, the original datasets and the information pertinent to adjustment may be limited or absent. This challenge is referred to as secondary analysis, in which the linkage is performed by one organization and the analyses of the linked data-sets are implemented by other organizations. The National Institute on Aging (NIA) LINKAGE project, which enables access to NIA-funded study data linked to administrative records from the Centers for Medicare & Medicaid Services, is a prominent example of such a setting (National Institute of Aging, 2025). We use the term post-linkage data analysis (PLDA) to describe statistical analyses of data that have been generated through record linkage or privacy preserving record linkage.
Multiple R packages have been developed to link records in the
absence of unique identifiers. Some of these applications are based on
the Fellegi-Sunter model (Fellegi & Sunter, 1969),
RecordLinkage (Sariyar & Borg, 2025),
fastLink (Enamorado et al., 2023). More recently,
Bayesian record linkage packages were developed, bstrl
(Betancourt et al.,
2022), BRL (Sadinle, 2017) and blink
(Steorts, 2020).
In addition, privacy preserving record linkage package is available
through PPRL (Schmiedel, 2025). To evaluate the quality of
the linkage, the ER-Evaluation package (Binette & Reiter,
2023) has also been developed in Python to assess linkage
quality using benchmark data.
Similarly, Python record linkage tools such as Splink
(Linacre et al.,
2022), recordlinkage (De Bruin, 2019),
dedupe (Gregg
& Eder, 2022), gfs_sampler (Farley & Gutman,
2020), and anonlink (Data61, 2017) have been developed to
offer flexible and scalable linkage solutions.
Beyond R and Python, record linkage functionality is available across
the statistical and commercial software landscape, including Stata
packages such as reclink2 (Wasi & Flaaen, 2015) and
dtalink (Kranker, 2018), SAS-based
tools like The Link King (Campbell et al., 2008), and
standalone applications such as Link Plus (Campbell et al.,
2008) and FEBRL (Christen, 2008).
Despite software available to perform and evaluate record linkage
across platforms, dedicated tools to conduct downstream statistical
analyses that account for potential linkage errors remain absent. A
common statistical tool to estimate the relationship between an outcome
and a set of covariates are linear and generalized linear models. These
models are efficiently implemented with the lm() and
glm() commands in R (R Core Team, 2024). Applying these methods
directly to linked datasets assumes that all records are correctly
linked, with no erroneous matches. Although statistical methods to
adjust for linkage errors have been proposed (Kamat & Gutman, 2026), the
corresponding software is often not readily accessible. This poses a
barrier to the adoption and application of these methods in practice.
Chambers (2009) provides R functions
to adjust for linkage errors within the regression framework. These
methods assume that the probability of a mismatch is constant across
records within defined subgroups (i.e., exchangeable linkage errors).
However, the available R functions are largely limited to estimating
equations for linear and logistic regressions, and are geared toward
simulation and methodological illustration. As such, they are not well
suited for routine use by analysts or subject-matter experts and may be
computationally inefficient.
We have developed an R package, postlink (Bukke et al.,
2026), that provides a cohesive and user-friendly
implementation of PLDA methods across a range of settings and input
types. The package is available from the Comprehensive R Archive Network
(CRAN) at https://CRAN.R-project.org/package=postlink, and is
implemented using a common regression modeling interface while adjusting
for falsely linked records. The postlink package focuses on
methods for secondary analysis. In primary analysis settings, where
individual files are accessible, methods that jointly perform record
linkage and analysis with direct propagation of uncertainty are more
appropriate.
The postlink R package consolidates three foundational
frameworks for post-linkage data analysis. First, the software
implements expectation-maximization (EM) based mixture modeling, which
absorbs the precursor package pldamixture (Bukke et al.,
2024) and facilitates contingency table analysis based on
methods developed in Slawski et al. (2025). Second, the package
incorporates a Bayesian mixture modeling approach proposed in Gutman et al. (2016). This approach enables the
estimation of model parameters after adjustments for linkage errors, as
well as using multiple imputations of the latent match status (correctly
vs. incorrectly linked records). The imputed match status can then be
used to fit post-linkage models not originally specified, using standard
multiple imputation combining rules (Rubin, 1996). Third, the package
integrates bias-adjusted estimating equations following the work of
Chambers (2009) for generalized
linear models and Vo et al. (2024) for Cox proportional hazards
regression. Within postlink, the methods in Chambers (2009) are extended to
encompass the Poisson and Gamma regression models. In addition, our
implementation improves the scalability of these methods for larger
datasets and their accessibility for analysts.
Subsequent sections of the paper are structured as follows. Section 2 provides an overview of the
postlink package. Section 3
details the three adjustment methods implemented in the R package. Section 4 illustrates the usage of the
package with a case study. We conclude and discuss future directions in
Section 5. Further demonstrations and details
regarding the package functionality are available on the
postlink website at https://postlink-group.github.io/postlink/. The
long-term goal of postlink is to progressively expand its
functionality to support a broad range of linkage and post-linkage
analysis scenarios.
2. The postlink Architecture
The postlink package provides an extensible
computational framework for post-linkage data analysis (PLDA). Given the
diversity in record linkage pipelines—both in underlying data structures
and in the amount of linkage information available—the software is built
using a modular, object-oriented R architecture. As illustrated in
Figure 1, the workflow leverages S3 method dispatch to decouple the
specification of linkage uncertainty from the substantive statistical
modeling. This separation serves three primary software design
objectives: (1) providing a unified interface for diverse post-linkage
scenarios; (2) allowing practitioners to propagate linkage uncertainty
using standard regression paradigms; and (3) separating the internal
mechanics of the package from its public application programming
interface (API) to facilitate extensibility.
The user workflow is consequently partitioned into two phases: 1) Adjustment Specification and 2) Estimation & Inference. In the first phase, the focus is on formulating the nature of the data and the linkage errors. The second phase focuses on the scientific questions that need to be addressed using the data.
Figure 1: Schematic of the typical two-phase
postlink workflow, illustrating the decoupling of the
adjustment specification from the outcome model fitting.
Phase 1: Adjustment Specification
In the first phase, the analyst defines the linked dataset and the
available information about the linkage process using a constructor
function. These functions, prefixed with adj, correspond to
different frameworks (e.g., weighting-based methods or mixture modeling
approaches). Section 3 describes the specific
adjustment methods and linkage error scenarios currently supported. The
constructors perform an initial validation layer to check for coherence
between the data structure and the assumptions of the chosen adjustment
methodology. Upon successful validation, the constructor returns a
structured S3 Adjustment object. To accommodate large-scale
data, this object is engineered as a lightweight container. This
construction minimizes memory overhead by retaining essential metadata
and accessing the underlying data by reference to avoid duplication. The
Adjustment object also implements context-specific methods
for standard generics (e.g., print(), plot()),
enabling diagnostic inspection of the linkage assumptions before model
fitting.
Phase 2: Estimation and Inference
The second phase handles model fitting through a suite of wrapper
functions designed to mirror standard R regression interfaces. The
current core functions include plglm() for generalized
linear models, plcoxph() for Cox proportional hazards,
plsurvreg() for parametric survival models, and
plctable() for contingency table analyses.
Using S3 dispatch, these wrappers retrieve the encapsulated
specifications from the Adjustment object and perform the
corresponding internal inference routine. This design allows analysts to
use the familiar formula-based syntax of core R functions (e.g.,
stats::glm(), survival::coxph()), while the
package manages the underlying computational complexity.
The resulting fitted model objects integrate with the R ecosystem,
supporting standard post-modeling generics such as
summary(), confint(), vcov(), and
predict(). For Bayesian adjustments, multiple imputation
pooling is supported with a generic mi_with() function.
3. Adjustment Methods for Secondary Analysis
We describe the theoretical foundation for the secondary analysis
methods that are supported in postlink. We assume that the
linked dataset, denoted as LD,
comprises of n record pairs. For entity
i \in \{1, \dots, n\}, we observe a
single outcome variable Y_{Ai} and a
set of covariates \mathbf{X}_{Ai}
originally recorded in data source A, alongside a set of covariates
\mathbf{X}_{Bi} originally recorded in
data source B. Formally, the linked data is defined as LD = \left(\mathbf{X}_{B}, \mathbf{Y}_A,
\mathbf{X}_A \right), where \mathbf{X}_{B} =
\{\mathbf{X}_{B1},\ldots,\mathbf{X}_{Bn}\}, \mathbf{Y}_{A} = \{Y_{A1},\ldots,Y_{An}\} and
\mathbf{X}_{A} =
\{\mathbf{X}_{A1},\ldots,\mathbf{X}_{An}\}. Suppose that an
analyst utilizes dataset LD to estimate
the conditional distribution p(Y_{Ai} \mid
\mathbf{X}_{Ai},\mathbf{X}_{Bi};\boldsymbol{\beta}), where \boldsymbol{\beta} is a set of unknown
parameters governing this distribution. As LD is formed by merging two separate data
sources without unique identifiers, it contains the observed paired
outcome Y_{Ai}^*, which may differ from
the true outcome Y_{Ai} because of
linkage errors (Y_{Ai}^* \neq Y_{Ai}).
As a result, LD includes both true and
false links, and only limited information on overall linkage quality may
be available (e.g., estimated mismatch rates within blocks).
When the original data sources are of the same size and they both represent the same entities, the observed vector \mathbf{Y}_{A}^* is a permuted version of true vector \mathbf{Y}_{A}, such that \mathbf{Y}_{A}^* = \Delta\mathbf{Y}_{A}, where \Delta is an n\times n permutation matrix with exactly one entry equal to one in each row and column. When there are no errors in the linkage, \Delta is the identity matrix and \mathbf{Y}_{A}^*=\mathbf{Y}_{A}. This notation cannot be utilized when the dataset includes some records that represent the same entities, and some records that do not. To address this limitation, we rely on an indicator variable C_i with C_i=1 for correct matches, i.e., Y_{A_i}^*=Y_{A_i}, and C_i=0 otherwise, 1 \leq i \leq n.
A naive estimator \widehat{\boldsymbol{\beta}}_N of the
parameter \boldsymbol{\beta} is
obtained by regressing \mathbf{Y}_A^*
on \mathbf{X} = [\mathbf{X}_A \;
\mathbf{X}_B]. This estimator is typically biased when \mathbf{Y}^*_{A} \neq \mathbf{Y}_{A} (Kamat & Gutman,
2026; Neter et al., 1965; Wang et al., 2022). To mitigate this
bias, postlink provides three distinct adjustment
frameworks, which are accessible through the package’s two-phase
adj*() and pl*() architecture. Table 1
provides an overview of these frameworks, their constructor functions,
and the statistical modeling they currently support.
Table 1: An overview of three PLDA frameworks implemented in
postlink
| Framework | Constructor | Primary Use Case | Supported Models |
|---|---|---|---|
| Weighting | adjELE() |
When block-level mismatch rates are given (e.g., from audit samples). | GLMs, CoxPH |
| Mixture Modeling | adjMixture() |
When continuous or categorical paradata on linkage quality (e.g., string distances) are available. | GLMs, CoxPH, Contingency Tables |
| Bayesian Mixture | adjMixBayes() |
When exploring multiple scientific models or to incorporate prior knowledge. | GLM, Parametric Survival |
3.1 Weighting
The first approach relies on weighting the estimating equations to correct for mismatch bias. Suppose that \mathbf{X} = \mathbf{X}_B, so that there are no covariates residing with the outcome in data source A. In addition, dataset LD is partitioned into Q blocks, where block q \in \{1,\dots, Q\} is of size n_q, and n = \sum_{q=1}^{Q} n_q. Block partitioning is commonly used to reduce the computational complexity of linkage algorithms. Blocks are formed by requiring exact agreement on categorical linking variables across datasets, such that all individuals within a block share the same values on these variables. Let \mathbf{Y}^{*}_{Aq}, \mathbf{X}_{Bq} and \Delta_q denote the response variable, the matrix of covariates, and a binary matrix of order n_q in block q, respectively. Under the non-informative linkage assumption (Chambers, 2009; Kamat & Gutman, 2026), \Delta_q is independent of \mathbf{Y}_{Aq} given \mathbf{X}_{Bq}, which implies that E(\mathbf{Y}_{Aq}^* \mid \mathbf{X}_{Bq}) = E(\Delta_q \mid \mathbf{X}_{Bq}) E(\mathbf{Y}_{Aq} \mid \mathbf{X}_{Bq}) = E_q f_q(\boldsymbol{\beta}), where E_q = E(\Delta_q \mid \mathbf{X}_{Bq}), and f_q(\boldsymbol{\beta}) = E(\mathbf{Y}_{Aq} \mid \mathbf{X}_{Bq}) is the regression function under correctly linked data.
When the data is correctly linked, an estimating equation for \boldsymbol{\beta} is of the form \sum_{q=1}^{Q} \mathbf{G}_q(\mathbf{X}_{Bq}, \boldsymbol{\beta}) \left\{\mathbf{Y}_{Aq}-f_q(\boldsymbol{\beta}) \right\}=0, \tag{1} where \mathbf{G}_q(\mathbf{X}_{Bq}, \boldsymbol{\beta}) is a weighting matrix that is independent of \mathbf{Y}_{Aq}. However, we only have access to \mathbf{Y}^*_{Aq}. Replacing \mathbf{Y}_{Aq} with \mathbf{Y}^*_{Aq} in Equation (1) results in the naive estimator \widehat{\boldsymbol{\beta}}_N. This estimator is biased when some of the units are incorrectly linked (Chambers, 2009). To correct for the bias in \widehat{\boldsymbol{\beta}}_N, we modify Equation (1) to adjust for linkage errors, \sum_{q} \mathbf{G}_q(\mathbf{X}_{Bq}, \boldsymbol{\beta}) \left\{ \mathbf{Y}^*_{Aq} - E_q f_q(\boldsymbol{\beta}) \right\} = 0. \tag{2} Chambers (2009) constructed E_q under the exchangeable linkage error (ELE) model [E_q]_{ij} = \begin{cases} \lambda_q, & \text{if record pair } (i,j) \text{ is a link},\\ \frac{1 - \lambda_q}{n_q - 1}, & \text{otherwise}, \end{cases} \tag{3} where \lambda_q is the proportion of correct links in block q, and it is either known to the analyst or estimated from audit samples.
The choice of \mathbf{G}_q(\mathbf{X}_{Bq}, \boldsymbol{\beta}) in Equation (2) leads to different estimators. For example, in linear regression, \mathbf{G}_q = \mathbf{X}_{Bq}^\top yields a ratio-type adjustment estimator, which is similar to the one proposed by Scheuren & Winkler (1993). Setting \mathbf{G}_q = \mathbf{X}_{Bq}^\top E_q results in the estimator proposed by Lahiri & Larsen (2005). Setting \mathbf{G}_q = \mathbf{X}_{Bq}^\top E_q^\top (\boldsymbol{\Sigma}_q + \mathbf{V}_q)^{-1} leads to the best linear unbiased estimator (BLUE) of \boldsymbol{\beta}, where \mathbf{V}_q = \mathrm{Var}(\boldsymbol{\Delta}_q E(\mathbf{Y}_q \mid \mathbf{X}_{Bq}) \mid \mathbf{X}_{Bq}) and \boldsymbol{\Sigma}_q is a diagonal matrix containing the observation-specific error variances, \sigma_i^2 = \mathrm{Var}(Y_i \mid \mathbf{X}_{Bi}), 1 \leq i \leq n. Since the BLUE depends on unknown quantities, it can be approximated using an iterative procedure in which these quantities are replaced by successively refined estimates.
This suite of weighting methods has been extended to generalized linear models (Chambers, 2009, pp. 22–23) and time-to-event models (Vo et al., 2024). Corresponding sandwich-type variance estimators have also been derived (Chambers, 2009).
In postlink, this procedure is implemented using the
two-phase workflow. First, the analyst specifies the ELE assumptions
using the adjELE() constructor:
adjELE(
linked.data,
m.rate, # Block-specific or global mismatch rates
audit.size, # Optional: Sample sizes for clerical review audit
blocks, # Vector defining the blocking structure
weight.matrix # Specifies the G_q matrix formulation
)The m.rate argument supplies block-specific or global
probabilities of mismatch, 1-\lambda_q.
The blocks argument is an indicator for the blocking
structure, which can be omitted if no blocks are present. If the
mismatch rates are estimated from a clerical review, the
audit.size argument can be supplied to accordingly inflate
the variance estimation. The form of the weighting matrix can be
specified in weight.matrix, allowing the user to select
between the ratio-type adjustment ("ratio"), the
Lahiri-Larsen estimator ("LL"), the best linear unbiased
estimator ("BLUE"), or "all" to compute all of
them simultaneously.
Once the adjustment object is validated and constructed, it is passed
to a standard modeling wrapper, such as plglm() for GLMs or
plcoxph() for survival analysis using Cox proportional
hazard model:
plglm(formula, family = "gaussian", adjustment, control, ...)
plcoxph(formula, adjustment, control, ...)The formula describes the relationship between \mathbf{Y}_q^* and \mathbf{X}_{Bq}, while the
family describes the link function. These arguments
correspond to the regression mean function f_q(\boldsymbol{\beta}). The
family is set to "gaussian" by default, but
can take the values "binomial", "poisson", or
"gamma". The control argument allows the user
to specify starting values for the iterative algorithm used to estimate
\boldsymbol{\beta}. By default, the
naive estimator is used as the initial value for ratio-type and
Lahiri–Larsen weight matrix adjustments. For the BLUE weight matrix, the
Lahiri–Larsen estimator serves as the default starting value for
estimating \boldsymbol{\beta}.
The backend implementation differs from that of Chambers (2009) primarily in how
Equation (2) is handled. By relying on vectorized operations for
block-wise processing and simplifying E_q
M into a sum of column means and scalar constants, the internal
computational routine reduces the computational complexity from \mathcal{O}(n_q^2 d) down to \mathcal{O}(n_q d), where d represents the number of covariates in the
model (see Appendix A). The nleqslv package (Hasselman, 2023) is
used to solve the non-linear system of equations.
3.2 Inferring true links via latent mixture modeling (EM)
An alternative set of methods that adjusts for linking errors attempts to estimate the configuration of match status indicators \mathbf{C}= \{C_i\}. For record pair i, let \mathbf{Z}_i denote auxiliary variables available at the secondary analysis stage that may be associated with C_i. Following Gutman et al. (2016) and Slawski et al. (2025), a two-component mixture model for \mathbf{C} can be specified as \begin{aligned} Y_{Ai},\mathbf{X}_{Ai},\mathbf{X}_{Bi} \mid \mathbf{Z}_i,C_{i}=1 & \sim p_m(Y_{Ai},\mathbf{X}_{Ai},\mathbf{X}_{Bi} \mid \mathbf{Z}_i;\boldsymbol{\beta}_m), \\ Y_{Ai},\mathbf{X}_{Ai},\mathbf{X}_{Bi} \mid \mathbf{Z}_i,C_{i}=0 & \sim p_u(Y_{Ai},\mathbf{X}_{Ai},\mathbf{X}_{Bi} \mid \mathbf{Z}_i;\boldsymbol{\beta}_{u}), \\ C_i & \sim \text{ Bernoulli }(h(\mathbf{Z}_i;\boldsymbol{\xi})), \end{aligned} \tag{4} where \boldsymbol{\beta}_m and \boldsymbol{\beta}_u govern the joint distributions p_m and p_u of (Y_{Ai}, \mathbf{X}_{Ai},\mathbf{X}_{Bi}) among true links and false links, respectively. Moreover, h(\cdot\,;\boldsymbol{\xi}) denotes the conditional probability of a true link given \mathbf{Z}_i, with \boldsymbol{\xi} denoting the associated parameters. Slawski et al. (2025) do not consider \mathbf{X}_{Ai} co-occurring with the outcome in data source A and model p_m(Y_{Ai}, \mathbf{X}_{Bi} \mid \mathbf{Z}_i; \boldsymbol{\beta}_m) and p_u(Y_{Ai}, \mathbf{X}_{Bi} \mid \mathbf{Z}_i;\boldsymbol{\beta}_u). In addition, they assume that p_u(Y_{Ai}, \mathbf{X}_{Bi} \mid \mathbf{Z}_i;\boldsymbol{\beta}_u) = p_{uY}(Y_{Ai};\boldsymbol{\beta}_u) \cdot p_{u\mathbf{X}}(\mathbf{X}_{Bi} \mid \mathbf{Z}_i;\boldsymbol{\beta}_u), and estimate \boldsymbol{\beta}_m by maximizing the likelihood resulting from (4) using the EM algorithm (Dempster et al., 1977) while a plug-in estimator is used for p_u. The function h(\cdot\,;\boldsymbol{\xi}) in Equation (4) can be parameterized in different ways. For instance, h(\cdot\,;\boldsymbol{\xi}) can be modeled according to a logistic regression model, with specific constraints on the linear predictors to facilitate estimation. An intercept-only specification of this logistic regression model corresponds to assuming a constant prior mismatch probability. Standard errors for the estimate of \boldsymbol{\beta}_m are computed using a normal approximation for the sampling distribution of the point estimator.
The EM-based mixture model requires specification using
adjMixture():
adjMixture(
linked.data,
m.formula, # Formula for the mismatch indicator model
m.rate, # Optional: global mismatch rate constraint
safe.matches # Optional: indicator for known true links
)The parameter m.formula takes a one-sided formula object
to model the conditional probability of a false link given the auxiliary
variables \mathbf{Z}_i, 1-h(\cdot\,;\boldsymbol{\xi}). The default is
an intercept-only model (~ 1), corresponding to a constant
marginal mismatch rate. The parameter m.rate allows the
user to supply an estimate to constrain the global mismatch rate. The
parameter safe.matches is a logical vector that indicates
records that are known to be correct matches, fixing their correct
linkage probability to 1 during the estimation process.
As the adjustment specification is decoupled from the modeling
syntax, the same Adjustment object can be passed to
plglm(), plcoxph(), or plctable()
for analyses. For example, fitting an adjusted contingency table
requires only the formula and the adjustment object:
plctable(formula, adjustment, control, ...)The internal routine retrieves the S3 object, checks for coherence, and dispatches the customized EM algorithm, which iteratively updates the posterior probability that a record is a correct match in the E-step, and maximizes model parameters in the M-step.
3.3 Bayesian mixture modeling and multiple imputation
While the EM approach yields robust maximum likelihood estimates for
a specified p_m, researchers may also
wish to propagate linkage uncertainty into downstream analyses that are
not easily accommodated within an EM framework. In addition, analysts
may use complex models to identify erroneous links but report results
from a simpler primary analysis model. To support these settings,
postlink implements a Bayesian mixture model (Gutman et al.,
2016) for estimating the parameters in Equation (4). Note
that in the Bayesian formulation, \mathbf{X}_i includes both \mathbf{X}_{Ai} and \mathbf{X}_{Bi}, and parameter constraints
can be imposed through the specification of prior distributions.
Posterior inference can be performed using the data augmentation
procedure (Tanner
& Wong, 1987), which alternates between imputing \mathbf{C} and updating the model parameters.
When p_m coincides with the analysis
model of interest, inference based on posterior draws of \boldsymbol{\beta}_m incorporates the
uncertainty about \mathbf{C}. In
settings where p_m differs from the
analysis model, post-linkage inference can be obtained using multiple
imputation (MI, Rubin,
1987). Multiple imputed realizations of \mathbf{C} are drawn from the posterior
distribution. Each realization is used to fit the analysis model
separately, and resulting point and interval estimates are combined
across imputations using standard multiple imputation combining rules
(Little &
Rubin, 2002; Rubin, 1996). The
underlying software for sampling from the posterior distribution is
based on Stan (Stan
Development Team, 2024) and the rstan package
(Stan Development Team,
2025), which employs the No-U-Turn Sampler to draw samples
from the posterior distribution of the parameters (Hoffman & Gelman,
2014). Given the kth draw of
the parameters \boldsymbol{\beta}^{(k)}_{m}, \boldsymbol{\beta}^{(k)}_{u}, and \boldsymbol{\xi}^{(k)}, imputed values of
\mathbf{C}^{(k)} are generated from the
corresponding Bernoulli distribution. Currently, the Bayesian mixture
model assumes that h(\cdot\,;\boldsymbol{\xi}) is constant
across all units.
The Bayesian adjustment is initialized using
adjMixBayes(), which optionally accepts customized prior
distributions and the current specifications of the prior distribution
for the different models is described in Appendix B,
adj_bayes <- adjMixBayes(
linked.data,
priors # Optional custom hyperparameters
)The model is then fitted using a standard wrapper, where the MCMC
iterations and computational cores can be controlled using the
control argument. In addition to GLMs, the package also
provides plsurvreg() for time-to-event outcomes. This
function extends the Bayesian mixture framework to survival data by
specifying p_m and p_u as parametric survival models. The
function call is
plsurvreg(formula, dist = "weibull", adjustment, control, ...)The parameters formula and data specify the
survival regression model, where the response must be a
Surv() object with right-censoring. The dist
option defines the parametric form of the survival outcome models p_m and p_u,
with current choices including "weibull" and
"gamma". The parameter control is a list that
specifies the MCMC settings, including the number of iterations,
burn-in, random seed, and number of computational cores. These options
control posterior sampling and can be tuned to balance computational
cost with estimation accuracy.
We note that the plcoxph() function implements the
semi-parametric Cox proportional hazards model, and treats the baseline
hazard non-parametrically. In contrast, the plsurvreg()
function uses a fully parametric survival model when performing Bayesian
inference.
Finally, to support multiple imputation, the package provides
mi_with() function:
mi_with(object, data, formula, family)This function extracts samples of the latent match indicators from
the glmMixBayes or survMixBayes objects. Using
these samples, it repeatedly fits a user-specified regression model
(formula and family) on the subset of records
classified as correct matches for each MCMC draw. Estimates across
imputations are aggregated using standard combination rules for multiple
imputation (Little
& Rubin, 2002; Rubin,
1996).
4. Illustrations
We illustrate the implementation of the different PLDA adjustment
methods and their varying inclusion of linkage paradata. To mimic a
secondary analysis environment where the original source files are
inaccessible, we provide a linked dataset LD1000 within
postlink. The LD1000 was generated through
one-to-one matching of records from two separate synthetic datasets.
R> data("LD1000", package = "postlink")
R> head(subset(LD1000, select = -c(True_Disease_Status, is_match)), n = 5)
#> # A tibble: 5 × 8
#> Disease_Status BMI Age Treatment fl_prob safe_match bm block_mrate
#> <int> <dbl> <int> <int> <dbl> <lgl> <dbl> <dbl>
#> 1 0 26.4 33 0 0.980 FALSE 12 0.45
#> 2 1 23.5 39 0 0.175 FALSE 6 0.35
#> 3 1 25.1 78 0 0.981 FALSE 6 0.35
#> 4 0 25.4 52 1 0.979 FALSE 3 0.25
#> 5 1 27.1 40 0 0.965 FALSE 4 0.35Dataset A consisted of 1,000 unique records and contained the
demographic and clinical predictors: BMI, Age,
and a binary Treatment indicator. Dataset B contained the
binary outcome variable, Disease_Status. To construct these
datasets, the RLdata10000 data from the
RecordLinkage package (Sariyar & Borg, 2025) were used
as identifiers of entities that possess two distinct records, which
often contain natural typographical errors. Of these identifiers, 800
were randomly selected to allocate the first instance to Dataset A and
the corresponding error-prone second instance to Dataset B. To complete
the 1,000-record cohorts, the remaining 200 records in Dataset A and 200
records in Dataset B were sampled from mutually exclusive sets of
identifiers representing unique, unmatchable entities. This design lead
to a mismatch rate of approximately 30% when the files were linked as
described below. Both datasets shared a set of error-prone
quasi-identifiers used for the linkage process: first name, last name,
and date of birth. The linkage process was strongly non-informative, so
that the linkage errors were independent of the variables used in the
post-linkage analysis (Kamat & Gutman, 2026).
The two files were linked probabilistically with the
fastLink package that implements of the Fellegi-Sunter
model (Fellegi &
Sunter, 1969). To maintain a cohort of n=1000 units and establish one-to-one matches
between records, the posterior probabilities from fastLink
were used to construct a cost matrix that was resolved via the
solve_LSAP function in the R package clue
(Hornik & Böhm,
2023). Within this linked dataset, a subset of records was
designated as “safe matches” if both the first and last names exhibit
exact string agreement across the datasets and the birth year exactly
agreed.
The primary objective is to fit a logistic regression model to
estimate the probability of the outcome (Disease_Status)
given the covariates (BMI, Age, and
Treatment):
R> formula_sci <- Disease_Status ~ BMI + Age + TreatmentFormally, the analysis model is as below, where both BMI
and Age are standardized to have a mean of zero and a
standard deviation of one for comparability of effect sizes across the
continuous and binary predictors. The true, unknown parameters are \beta_{0}=-1, \beta_{\texttt{BMI}}=1.2, \beta_{\texttt{Age}}=0.8, and \beta_{\texttt{Treatment}}=-1.5:
\text{logit}(P(Y_i = 1 \mid \texttt{BMI}_i, \texttt{Age}_i, \texttt{Treatment}_i)) = \beta_0 + \beta_{\texttt{BMI}}\texttt{BMI}_i + \beta_{\texttt{Age}}\texttt{Age}_i + \beta_{\texttt{Treatment}}\texttt{Treatment}_i
LD1000$BMI <- as.numeric(scale(LD1000$BMI))
LD1000$Age <- as.numeric(scale(LD1000$Age))In the secondary analysis setting, the analyst has access to the
variables of interest merged from both files
(Disease_Status, BMI, Age,
Treatment), optionally along with linkage paradata (e.g.,
posterior match probabilities fl_prob from the underlying
fastLink RL algorithm, safe match indicator
safe_match, blocking variable bm, or
block-wise mismatch rates block_mrate), but does not
observe the latent true match indicator (is_match) or true
linked values (True_Disease_Status). Rather than the raw
identifiers used for linking the records, which are typically
unavailable because of privacy constraints, the models rely on these
paradata as available to assess linkage quality; for example, Figure 2
shows the distribution of the fl_prob across correct
matches (is_match = TRUE) and mismatches
(is_match = FALSE).
Figure 2: Horizontal boxplots of fastLink posterior
match probabilities (fl_prob) in LD1000 by
latent match status (C = 0 for
mismatches; C = 1 for correct matches),
with blue diamonds denoting group means and semi-transparent points
indicating outliers beyond 1.5 times
the interquartile range. As most mismatches concentrate at 0 (where the C =
0 median and first quartile coincide at the left box edge) and
most correct matches concentrate near 1
(compressing the C = 1 box), record
pairs with intermediate or misclassified probabilities appear as
discrete outlier points whose darker shading reflects overlapping
observations.
We demonstrate the overall workflow for different adjustment methods using the available paradata and compare them to cases where we know the truth or use naive methods.
4.1 The Benchmarks (Models 1–3)
The first three models represent standard approaches that do not
utilize the postlink package:
Model 1 is the “Oracle”. This is logistic regression fit on the unobservable, true linked data.
R> fit_oracle <- glm(True_Disease_Status ~ BMI + Age + Treatment,
+ data = LD1000, family = "binomial")Model 2 is a standard “Naive” logistic regression fit on the probabilistically linked data that ignores the presence of mismatches.
R> fit_naive <- glm(formula_sci, data = LD1000, family = "binomial")Model 3 is a “Naive (Safe Matches Only)” approach, in which the analyst restricts the regression only to links deemed “safe”, discarding the rest of the dataset.
4.2 Aggregate Information (Model 4)
If the analyst has access to block-level mismatch rates based on
clerical review and knowledge of the random audit sample sizes per
block, the adjELE() constructor can be used. The blocks are
defined as the birth month (bm) variable.
R> adj_ele_audit <- adjELE(linked.data = LD1000,
+ m.rate = LD1000$block_mrate,
+ blocks = bm,
+ audit.size = rep(20, length(unique(LD1000$bm))),
+ weight.matrix = "BLUE")
R> print(adj_ele_audit)
#>
#> Adjustment Object: Exchangeable Linkage Errors
#>
#> * Linked Data:
#> Observations: 1,000
#> * Specification:
#> Weight Matrix: BLUE
#> Blocks: 12 distinct blocks
#> Mismatch Rate: Variable (Mean: 0.271, Range: 0.100 - 0.450)
#> Audit Sample: Same size 20In practice, the choice of blocking variables for post-linkage adjustment generally aligns with the blocking algorithms used by the underlying record linkage algorithm, because mismatch rates tend to be more similar within these blocks.
R> fit_ele_audit <- plglm(formula_sci, family = "binomial",
+ adjustment = adj_ele_audit)
R> summary(fit_ele_audit)
#>
#> Call:
#> plglm(formula = formula_sci, family = "binomial", adjustment = adj_ele_audit)
#>
#> Family: binomial
#>
#> --- Weighting Method: BLUE ---
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -0.841 0.143 -5.89 3.8e-09 ***
#> BMI 1.278 0.177 7.24 4.5e-13 ***
#> Age 0.719 0.125 5.76 8.4e-09 ***
#> Treatment -1.526 0.264 -5.78 7.5e-09 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The fit_ele_audit object returns results for \boldsymbol{\beta} = (\beta_{0},
\beta_{\texttt{BMI}}, \beta_{\texttt{Age}},
\beta_{\texttt{Treatment}}) in the analysis model, similar to the
standard glm object in R.
4.3 Probabilistic Paradata (Models 5–8)
If continuous or categorical paradata from the linkage process is
retained, it can be passed into the adjMixture()
constructor. We demonstrate four progressively information-rich
applications of the EM mixture model.
Model 5 relies on a anticipated global mismatch rate (0.3), which may be given by the linked data
provider based on a mean(LD1000$block_mrate) of 0.271:
R> adj_em_rate <- adjMixture(linked.data = LD1000,
+ m.formula = ~ 1,
+ m.rate = 0.3)
R> print(adj_em_rate)
#>
#> Adjustment Object: Mixture Model
#>
#> * Linked Data:
#> Observations: 1,000
#> * Specification:
#> Mismatch Model: ~1
#> Global Mismatch Rate: 0.3 (Constrained)
#> Safe Matches: None specifiedModel 6 alternatively incorporates a logical indicator for safe matches, allowing the algorithm to fix the match probability to 1 for highly confident linked records:
R> adj_em_safe <- adjMixture(linked.data = LD1000,
+ m.formula = ~ 1,
+ safe.matches = safe_match)
R> print(adj_em_safe)
#>
#> Adjustment Object: Mixture Model
#>
#> * Linked Data:
#> Observations: 1,000
#> * Specification:
#> Mismatch Model: ~1
#> Global Mismatch Rate: Unconstrained (Estimated from data)
#> Safe Matches: 84 (8.4%)In Model 7, we drop the global mismatch rate constraint and instead
incorporate the posterior match probabilities obtained from the
underlying record linkage algorithm (i.e., fastLink). This
allows the mixture model to predict the match indicators (C_i) using the linkage paradata (\mathbf{Z}_i), which are the auxiliary
covariates introduced in Equation (4).
R> adj_em_prob <- adjMixture(linked.data = LD1000,
+ m.formula = ~ fl_prob)
R> print(adj_em_prob)
#>
#> Adjustment Object: Mixture Model
#>
#> * Linked Data:
#> Observations: 1,000
#> * Specification:
#> Mismatch Model: ~fl_prob
#> Global Mismatch Rate: Unconstrained (Estimated from data)
#> Safe Matches: None specified Model 8 shows how the linkage paradata can also be used in combination, if all are available.
R> adj_em_all <- adjMixture(linked.data = LD1000,
+ m.formula = ~ fl_prob,
+ safe.matches = safe_match,
+ m.rate = 0.3)
R> print(adj_em_all)
#>
#> Adjustment Object: Mixture Model
#>
#> * Linked Data:
#> Observations: 1,000
#> * Specification:
#> Mismatch Model: ~fl_prob
#> Global Mismatch Rate: 0.3 (Constrained)
#> Safe Matches: 84 (8.4%)The fit_em_all object is summarized below. The outcome
model coefficients correspond to \boldsymbol{\beta} = (\beta_{0},
\beta_{\texttt{BMI}}, \beta_{\texttt{Age}},
\beta_{\texttt{Treatment}}). The mismatch model coefficients are
-\boldsymbol{\xi}, based on Equation
(4) in Section 3.2. Below, fl_prob is taken as a predictor
in this model based the specified adj_em_all. While the
negative estimate supports that higher fl_prob
probabilities predict correct matches, the concentration of values near
0 and 1 can inflate the standard error, contributing to a
non-significant p-value (see Figure
2).
R> fit_em_all <- plglm(formula_sci, family = "binomial",
+ adjustment = adj_em_all)
R> summary(fit_em_all)
#>
#> Call:
#> plglm(formula = formula_sci, family = "binomial", adjustment = adj_em_all)
#>
#> Outcome Model Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -0.859 0.139 -6.19 6.2e-10 ***
#> BMI 1.130 0.131 8.64 < 2e-16 ***
#> Age 0.722 0.122 5.94 2.9e-09 ***
#> Treatment -1.511 0.271 -5.57 2.5e-08 ***
#>
#> Mismatch Model Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 0.989 0.600 1.65 0.099 .
#> fl_prob -9.416 20.482 -0.46 0.646
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Average Correct Match Probability: 0.798In earlier Models 5-7, the fit object returned is similarly
formatted. The average correct match probability (0.798) refers to the mean of the posterior
correct match probabilities for all of the records. The probability for
each record is accessible through the fit object in the
match.prob vector. For safe matches, the probability is set
by default to 1.
R> str(fit_em_all$match.prob)
#> Named num [1:1000] 1 0.447 1 1 1 ...
#> - attr(*, "names")= chr [1:1000] "1" "2" "3" "4" ...4.4 Prior Knowledge (Model 9)
If the analyst has historical data or domain expertise regarding the
scientific parameters, the adjMixBayes() constructor allows
for the integration of this prior knowledge. In Model 5, the match rate
is set to a scalar. The Bayesian mixture framework accommodates
probability distributions that represent the uncertainty of the analyst
regarding the true match rate. In settings where relationship between
Y_{A} and \mathbf{X}_B is weak, it may be hard to
distinguish between true and false links without large number of units.
However, incorporating paradata in the form of prior distributions, such
as knowing the effect of age is positive and bounded or the approximate
error rate, can assist in the estimation process and result in more
reliable parameter estimates. The prior argument is a list
of prior distributions for model parameters, and it depends on the
family of distributions (See Appendix B).
In many record linkage applications one would expect that Y_{Ai} is independent or approximately
independent from \mathbf{X}_{Bi} among
false links (Guha
et al., 2022; Slawski et al.,
2025). This assumption may not be valid when the variables
used to create blocks are associated with Y_{Ai} and \mathbf{X}_{Bi} (Kamat et al., 2023). For this
simulation, we will maintain this independence assumption to demonstrate
the specification of prior distributions. Implementing this assumption
in the form of prior distribution assumes that a-priori the intercept
among mismatched links, intercept2, is centered at -1.09 and standard deviation of 3. The center
of the prior distribution is derived from the empirical log-odds of the
outcome prevalence in the linked dataset. The standard deviation of 3
result in relatively diffused prior distribution for the intercept of a
logistic regression model. To enforce independence between Y_{Ai} and \mathbf{X}_{Bi} among mismatched links, we
assumed that slope parameters, beta2, follow \mathcal{N}(0, 0.01), which reflects that the
coefficient is centered at zero with very small standard deviation. This
prior distribution, requires significant input from the observed data to
deviate from 0. Lastly, to incorporate that approximately 70\% of the links are true links, we assumed
that the prior distribution for the marginal proportion of true links,
theta, follows a Beta distribution with shape parameters 70
and 30. This implies that the mean of the distribution is 0.70 and the
standard deviation is 0.05.
R> adj_bayes <- adjMixBayes(linked.data = LD1000,
+ priors = list(
+ intercept2 = "normal(-1.09, 3)",
+ beta2 = "normal(0,0.01)",
+ theta = "beta(70,30)"))
R> print(adj_bayes)
#>
#> Adjustment Object: Bayesian Mixture
#>
#> * Linked Data:
#> Observations: 1,000
#> * Priors:
#> intercept2 : normal(-1.09, 3)
#> beta2 : normal(0, 0.01)
#> theta : beta(70, 30)
#> (Unspecified parameters will use defaults below)
#>
#> Defaults applied during fitting (for any unspecified):
#> Intercept: intercept ~ normal(0,10)
#> GLM Slopes: beta ~ normal(0,5) [binomial: normal(0,2.5)]
#> Surv Slopes: beta ~ normal(0,5) [weibull: normal(0,2)]
#> Mix Weight: theta ~ beta(1,1)The Bayesian mixture model returns posterior summaries for two latent components. The coefficients for the correct-match component (Component 1) correspond to the analysis model of primary interest. For the mismatch component (Component 2), the posterior means for the slope coefficients are centered near zero, which is consistent with the prior specification of conditional independence between the outcome and covariates among falsely linked records.
R> fit_bayes <- plglm(formula_sci, family = "binomial",
+ adjustment = adj_bayes,
+ control = list(seed = 2026))
R> summary(fit_bayes)
#> Call:
#> plglm(formula = formula_sci, family = "binomial", adjustment = adj_bayes,
#> control = list(seed = 2026))
#>
#> Family:
#> binomial
#>
#> (Component 1 = Correct-match):
#> Outcome Model Coefficients:
#> Estimates Std. Error 2.5 % 97.5 %
#> (Intercept) -0.788 0.284 -1.374 -0.262
#> BMI 1.412 0.252 0.963 1.951
#> Age 0.805 0.169 0.504 1.170
#> Treatment -1.683 0.354 -2.448 -1.063
#>
#> (Component 2 = Incorrect-match):
#> Outcome Model Coefficients:
#> Estimates Std. Error 2.5 % 97.5 %
#> (Intercept) -1.309 0.583 -2.599 -0.394
#> BMI 0.000 0.010 -0.020 0.019
#> Age 0.000 0.010 -0.018 0.019
#> Treatment 0.000 0.010 -0.020 0.020Note that these Component 2 outcome coefficients (\boldsymbol{\beta}_u) should be distinguished from the mismatch model coefficients (\boldsymbol{\xi}); the former describe the outcome relationship among false links, while the latter correspond to the probability of being a mismatch.
4.5 Results and Discussion
The coefficient estimates and 95% confidence intervals for all scenarios are visualized in Figure 3 and detailed in Table 2. These results illustrate the potential biases introduced by linkage errors and the relative usage of different post-linkage adjustment strategies.
Ignoring the \sim 30% mismatch rate (Model 2) results in estimates that are biased towards zero. For example, the estimate for \beta_{BMI} is attenuated by approximately 30% compared to the Oracle benchmark (from 1.163 by the Oracle down to 0.797 by the Naive), and its 95% confidence interval fails to include the true parameter. Relying on this unadjusted analysis would lead to underestimation of the relationship between the outcome and the covariate.
An alternative heuristic approach is to restrict the analysis to only high-confidence matches (Model 3). While this successfully eliminates the false links, it discards a significant portion of the dataset. By requiring exact name matches, the effective sample size is reduced from N=1000 to N=84, causing the standard error for the Treatment effect to more than triple (from 0.190 to 0.643). The confidence intervals are wider and the statistical power of the analysis is reduced. Moreover, when units with lower-confidence matches have different relationship between \mathbf{Y}_A and \mathbf{X}_B than high-confidence matches, bias may be introduced.
The adjustment methods implemented in postlink mitigate
these limitations to varying degrees. The ELE approach (Model 4) shifts
the estimates toward the true parameters’ values. Since the
audit.size argument was supplied, the standard errors
increase to account for the sampling variance inherent in the clerical
review, providing more appropriately calibrated inference. The EM
Mixture model using posterior match probabilities (Model 8) demonstrates
favorable performance. By incorporating continuous paradata, the
algorithm assigns probabilistic weights to the reliability of individual
records, recovering point estimates comparable to the Oracle model
without the large variance inflation observed in the exact matches only
approach. We note, however, that the simulated separation between the
high string distances of true matches and the near-zero distances of
forced random mismatches creates a comparatively straightforward
classification task. In applied settings where false links share higher
string similarity (e.g., matching ‘John Smith’ to ‘Jon Smith’
erroneously), the resulting overlap of the distributions p_m and p_u
in (4) would naturally increase the standard errors of the estimated
model parameters. The point estimates of the Bayesian model are
comparable to Model 5. However, because the prior distributions allow
for more variability in the mismatch rate and the estimates of
regression parameters among non-links leads to larger standard errors,
but still smaller than the high confidence matches analysis (Model
3).
Table 2: Logistic regression estimates and standard errors (SE) from the benchmark models and post-linkage adjustment methods.
| Model | Est. (Intercept) | SE (Intercept) | Est. (BMI) | SE (BMI) | Est. (Age) | SE (Age) | Est. (Treatment) | SE (Treatment) |
|---|---|---|---|---|---|---|---|---|
| 1. Oracle | -0.870 | 0.112 | 1.163 | 0.108 | 0.799 | 0.096 | -1.604 | 0.190 |
| 2. Naive (All Links) | -0.835 | 0.104 | 0.797 | 0.090 | 0.518 | 0.084 | -1.015 | 0.165 |
| 3. Naive (Safe Matches Only) | -0.855 | 0.380 | 0.838 | 0.334 | 0.506 | 0.302 | -1.522 | 0.643 |
| 4. ELE (Audit Sample) | -0.841 | 0.143 | 1.278 | 0.177 | 0.719 | 0.125 | -1.526 | 0.264 |
| 5. EM (Global Rate) | -0.841 | 0.146 | 1.219 | 0.261 | 0.713 | 0.141 | -1.483 | 0.341 |
| 6. EM (Safe Matches Only) | -0.837 | 0.136 | 1.129 | 0.198 | 0.666 | 0.120 | -1.400 | 0.288 |
| 7. EM (Posterior Probabilities) | -0.858 | 0.138 | 1.131 | 0.131 | 0.722 | 0.121 | -1.512 | 0.270 |
| 8. EM (All) | -0.859 | 0.139 | 1.130 | 0.131 | 0.722 | 0.122 | -1.511 | 0.271 |
| 9. Bayesian (Informative Priors) | -0.788 | 0.284 | 1.412 | 0.252 | 0.805 | 0.169 | -1.683 | 0.354 |
4.6 Decoupling Linkage and Analysis via Multiple Imputation
While the adjMixBayes() constructor performs joint
estimation of the linkage and outcome models, MCMC sampling is
computationally intensive. In practice, analysts frequently conduct
exploratory model building, assessing various interactions, polynomials,
or covariate combinations. Repeatedly fitting the joint Bayesian model
for each exploratory iteration introduces a computational
bottleneck.
The postlink package addresses this computational
complexity by using multiple imputation (Gutman et al., 2016), which is
implemented with the mi_with() generic function. Since the
Bayesian mixture model returns posterior draws of the latent match
status for every record, each draw effectively generates one completed
dataset. The mi_with() function extracts these latent
assignments, fits a user-specified regression model on the subset of
records classified as correct matches for each MCMC draw, and aggregates
the results using common combination rules for multiple imputation (Rubin, 1987). This
allows the analyst to run the MCMC procedure once and subsequently
explore alternative scientific models more efficiently. For example, if
the analyst is only interested in the marginal effect of
Treatment, they can pass the updated formula to the
pre-computed adjustment object:
R> formula_alternative <- Disease_Status ~ Treatment
R> pooled_mi_fit <- mi_with(object = fit_bayes,
+ data = LD1000,
+ formula = formula_alternative,
+ family = binomial())
R> print(pooled_mi_fit, digits = 3)
#> Pooled regression results across posterior match classifications:
#> Retained imputations (m): 9000
#>
#> Estimate Std.Error CI.lwr CI.upr df
#> (Intercept) -0.541 0.196 -0.925 -0.156 18844
#> Treatment -1.182 0.238 -1.649 -0.715 61465As detailed in the output, the marginal term for
Treatment is estimated to be significantly different from
zero at the 5% nominal level (95% CI: -1.649, -0.715), which is
consistent with the underlying data-generating process. This model
exploration utilized 9,000 posterior draws as imputed datasets,
demonstrating how linkage error estimation can be decoupled from
downstream scientific modeling. When using multiple imputation one
should attempt to avoid uncongenial models (Meng, 1994). Thus, it is
advisable to include many relationships of substantive interest when
classifying between true and false links, so that the analysis model is
more restrictive than the imputation model (Kamat & Gutman, 2026; Meng, 1994).
Figure 3: Correction for Attenuation Bias from Linkage Errors in Logistic Regression. The dashed line represents the true parameter values.
5. Conclusion
Record linkage is widely used to combine information of interest
scattered across multiple datasets. Datasets that were created using
linkage algorithms that are based on quasi-identifiers, often suffer
from linkage errors. These datasets are frequently analyzed without
regard for these linkage errors. The lack of easy-to-use software
implementing existing methods adjusting for such errors is likely a
major contributing factor. This lack is in contrast to the availability
of software for performing linkage. In this paper, we introduce the
postlink package, which equips researchers with three
state-of-the-art methods for post-linkage regression analysis through an
interface that mimics standard R modeling functions (e.g.,
stats::glm(), survival::coxph()). Our
implementations offer a variety of ways to incorporate information about
linkage quality in secondary analysis settings, which we expect to be
useful for practitioners working with linked data sets. The package
leaves room for various extensions, including additional types of
post-linkage analysis models, support for missed matches (i.e., records
that could not be linked) in addition to mismatches, and integration
with record linkage software. Alongside such enhancements, developing
tailored goodness-of-fit metrics, residual diagnostics, and exploratory
tools to evaluate adjustment specifications remains an important
direction for future research.
Computational Details
The results and illustrations in this article were obtained using R
4.5.1 and the postlink package version 0.1.2. R itself and
all standard dependent packages, including survival for
time-to-event outcomes, and nleqslv for solving estimating
equations, are available from the Comprehensive R Archive Network (CRAN)
at https://CRAN.R-project.org/.
The Bayesian mixture modeling routines in postlink use
Hamiltonian Monte Carlo sampling and require a functional C++ toolchain.
These routines depend on the rstan and
rstantools packages, which are also available from CRAN.
Cross-platform floating point differences may cause local numerical
deviations. For reproducibility, specific random seeds were set prior to
model fitting, and all results were generated using a Ubuntu 24.04
Docker container.
The postlink package and datasets used for the
application are open-source and publicly available. Full replication
scripts and the container Dockerfile for all results presented in this
article are provided in the supplementary materials.
Acknowledgments
Priyanjali Bukke, Jiahao Cui, Roee Gutman, and Martin Slawski were supported by NSF grant # 2411270. Bukke and Slawski were additionally supported by NSF grant # 2120318.
Appendix A: Computational and variance estimation details for the ELE model
A.1 Computation of expected design matrices
To solve the adjusted estimating equations (2), we compute the product of the expected linkage matrix E_q and the block-level design matrix \mathbf{X}_{Bq} (or a function of it, denoted by M_q). Constructing the n_q \times n_q matrix E_q and performing standard matrix multiplication requires \mathcal{O}(n_q^2 d) floating-point operations (flops) per block, which is computationally prohibitive for large datasets.
The postlink package circumvents this by exploiting the
structural properties of the ELE model. Let \alpha_q = 1 - \lambda_q represent the
mismatch rate. The product E_q M_q can
be analytically rewritten as a linear combination of the observed values
and their block-specific column means:
E_q M_q = c_{1q} M_q + c_{0q} \bar{M}_q, \tag{5}
where \bar{M}_q represents the
column means of M_q, and the scalars
are defined as c_{1q} = 1 - \alpha_q -
\alpha_q / (n_q - 1) and c_{0q} =
\alpha_q \frac{n_q}{n_q - 1}. This formulation reduces the
computational complexity to \mathcal{O}(n_q
d). The postlink implementation applies this
optimization globally across the internal glmELE() and
coxphELE() routines, using vectorized aggregation functions
(e.g., rowsum()) for model fitting and risk-set
evaluation.
A.2 Variance estimation and audit sample uncertainty
The variance of the ELE estimators is computed using a sandwich estimator, \widehat{\text{Var}}(\widehat{\boldsymbol{\beta}}) = J^{-1} V_H (J^{-1})^\top, where J is the Jacobian matrix of the estimating equations evaluated at \widehat{\boldsymbol{\beta}}. The variance of the estimating functions, V_H, is decomposed as V_H = V_1 + V_2.
The term V_2 represents the
empirical covariance of the estimating functions assuming the correct
match rates \lambda_q are fixed and
known. When these rates are estimated via a simple random audit sample
of size m_q drawn from the n_q records in block q (supplied via the audit.size
argument), the package computes V_1 to
account for this additional source of uncertainty. A finite population
correction (FPC) is applied to reflect the sampling fraction:
V_1 = \sum_{q=1}^Q \left( \frac{1}{m_q} - \frac{1}{n_q} \right) \left(
\frac{m_q}{m_q - 1} \right) \frac{\alpha_q}{\lambda_q^3} U_q U_q^\top,
\tag{6}
where U_q represents the
partial derivative of the block-specific adjusted estimating equation
with respect to \lambda_q, evaluated at
\hat{\boldsymbol{\beta}} and \hat{\lambda}_q. This accounts for audit
uncertainty so that the standard errors returned by standard
post-estimation generics, such as confint() and
vcov(), provide valid statistical inference even when
linkage quality is estimated from small clerical reviews (Vo et al.,
2024).
Appendix B: Default Prior Distributions for Each Family of Distributions
Bayesian adjustment methods require the specification of the prior
distribution. The specifications of all the prior distributions follow
the Stan programming language syntax (Stan Development Team, 2024). The prior
distributions include independent prior distributions for the intercept
among true links, intercept1, and the intercept among
non-links intercept2. To reduce the burden on analysts, we
assume the same independent prior distribution for all the slope
coefficients among true links, beta1, and the same
independent prior distribution among all false links,
beta2. All of the models support a prior distribution on
the marginal probability of being a true link, theta. The
other prior distributions are specific for each of the models. The prior
for GLM regressions are:
"gaussian" = list(
intercept1 = "normal(0,10)",
intercept2 = "normal(0,10)",
beta1 = "normal(0,5)",
sigma1 = "cauchy(0,2.5)",
beta2 = "normal(0,5)",
sigma2 = "cauchy(0,2.5)",
theta = "beta(1,1)"
)
"poisson" = list(
intercept1 = "normal(0,10)",
intercept2 = "normal(0,10)",
beta1 = "normal(0,5)",
beta2 = "normal(0,5)",
theta = "beta(1,1)"
)
"binomial" = list(
intercept1 = "normal(0,10)",
intercept2 = "normal(0,10)",
beta1 = "normal(0,2.5)",
beta2 = "normal(0,5)",
theta = "beta(1,1)"
)
"gamma" = list(
intercept1 = "normal(0,10)",
intercept2 = "normal(0,10)",
beta1 = "normal(0,5)",
beta2 = "normal(0,5)",
phi1 = "gamma(2,0.1)",
phi2 = "gamma(2,0.1)",
theta = "beta(1,1)"
)The prior distributions for the parametric time-to-event models:
"gamma" = list(
intercept1 = "normal(0,10)",
intercept2 = "normal(0,10)",
beta1 = "normal(0,5)",
beta2 = "normal(0,5)",
phi1 = "exponential(1)",
phi2 = "exponential(1)",
theta = "beta(1,1)"
)
"weibull" = list(
intercept1 = "normal(0,10)",
intercept2 = "normal(0,10)",
beta1 = "normal(0,2)",
beta2 = "normal(0,2)",
shape1 = "gamma(2,1)",
shape2 = "gamma(2,1)",
scale1 = "gamma(2,1)",
scale2 = "gamma(2,1)",
theta = "beta(1,1)"
)