Cox regression with linked data

Thanh Huan Vo, Valérie Garès, Li-Chun Zhang, André Happe, Emmanuel Oger, Stéphane Paquelet, Guillaume Chauvet

<figure><img src="image_1.png" /></figure>

► To cite this version:

Thanh Huan Vo, Valérie Garès, Li-Chun Zhang, André Happe, Emmanuel Oger, et al.. Cox re-
gression with linked data. Statistics in Medicine, 2023, Statistics in Medicine, 43 (2), pp.296-314.
10.1002/sim.9960 . hal-04336258

HAL Id: hal-04336258

https://hal.science/hal-04336258v1

Submitted on 27 Jan 2026

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
RESEARCH ARTICLE

Statistics
in Medicine

WILEY

Cox regression with linked data

Thanh Huan Vo<fnref n="1" />,<fnref n="2" />,✉

Valérie Garès<fnref n="1" />

Li-Chun Zhang<fnref n="3" />,<fnref n="4" />

André Happe<fnref n="5" />

Emmanuel Oger<fnref n="5" />

Stéphane Paquelet<fnref n="2" />

Guillaume Chauvet<fnref n="6" />

**Correspondence**

Guillaume Chauvet, Univ Rennes, ENSAI, CNRS, IRMAR-UMR 6625, F-35000 Rennes, France.

Email: guillaume.chauvet@ensai.fr

Record linkage is increasingly used, especially in medical studies, to combine data from different databases that refer to the same entities. The linked data can bring analysts novel and valuable knowledge that is impossible to obtain from a single database. However, linkage errors are usually unavoidable, regardless of record linkage methods, and ignoring these errors may lead to biased estimates. While different methods have been developed to deal with the linkage errors in the generalized linear model, there is not much interest on Cox regression model, although this is one of the most important statistical models in clinical and epidemiological research. In this work, we propose an adjusted estimating equation for secondary Cox regression analysis, where linked data have been prepared by a third-party operator, and no information on matching variables is available to the analyst. Through a Monte Carlo simulation study, the proposed method is shown to lead to substantial bias reductions in the estimation of the parameters of the Cox model caused by false links. An asymptotically unbiased variance estimator for the adjusted estimators of Cox regression coefficients is also proposed. Finally, the proposed method is applied to a linked database from the Brest stroke registry in France.

**KEYWORDS**

adjusted estimating equation, Cox regression, linkage error, secondary analysis, variance estimation

# 1 | INTRODUCTION

Record linkage, also known as data matching, is a process of combining data from different sources that refer to the same individuals or entities. Nowadays, data are collected everywhere by different sectors, and the ability of combining information from several databases can lead to novel knowledge for analysts. For example, record linkage is widely used in epidemiology and medical studies to enrich data on clinical performance and other health-related information.<fnref n="1" />,<fnref n="2" /> In national censuses, population data files obtained at different times can be linked to create longitudinal data sets.<fnref n="3" /> Record linkage may also be applied early in a survey to link the sampling frame and administrative data.<fnref n="4" /> The linked data allows for statistical analysis (eg, Cox regression) which would not be possible with data collected solely by means of the survey.

The record linkage process is straightforward if unique identifiers (eg, Social Security Number) are available and free of error in both databases. However, this information is often not available, or sometimes cannot be used due to ethical reasons. In such cases, record linkage methods may only use partial identifying information shared between databases, such as name, address, and gender. The variables used for comparison are called matching variables. Over the last decades, several methods have been developed to link data efficiently,<fnref n="5" />,<fnref n="6" /> such as the frequentist approach⁷⁻⁹ and the

This is an open access article under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2023 The Authors. *Statistics in Medicine* published by John Wiley & Sons Ltd.

<footnote n="1">¹Univ Rennes, INSA, CNRS, IRMAR-UMR 6625, Rennes, France</footnote>

<footnote n="2">²IRT b<>com, Laboratoire IA, Rennes, France</footnote>

<footnote n="3">³Department of Social Statistics and Demography, University of Southampton, Southampton, UK</footnote>

<footnote n="4">⁴Statistics Norway, Oslo, Norway</footnote>

<footnote n="5">⁵Univ Rennes, EA 7449 REPERES, Rennes, France</footnote>

<footnote n="6">⁶Univ Rennes, ENSAI, CNRS, IRMAR-UMR 6625, Rennes, France</footnote>
Bayesian approach.<sup>10,11</sup> However, because the matching variables are not unique and are likely to contain inaccuracies, linkage errors are unavoidable. The two kinds of record linkage errors are false links (false positives, ie, a non-matched pair predicted as a link), and missed links (false negatives, ie, a matched pair failed to be predicted as a link). Ignoring these errors may cause substantial bias in the analysis model,<sup>12</sup> causing misleading inference. It is therefore important to account for linkage errors in statistical analysis.

In published literature, two positions are usually considered to account for linkage errors in statistical analysis. Under the primary analysis framework, the data analyst is supposed to be granted access to the full linkage process, including knowledge of matching variables. From this perspective, Scheuren and Winkler<sup>13</sup> made use of the two highest matching weights of each record pair to reduce the bias of ordinary least square estimators under a linear regression model. However, the proposed estimators are not unbiased in full generality. Lahiri and Larsen<sup>14</sup> discussed this problem and proposed unbiased estimators in the same context, using the posterior matching probabilities obtained from the Fellegi-Sunter record linkage model. Hof and Zwinderman<sup>15</sup> extended the method by Lahiri and Larsen<sup>14</sup> for multiple links, and also proposed alternative estimators based on weighted least square methods, both for linear and logistic regression models. Recently, Han and Lahiri<sup>16</sup> adapted the approach by Lahiri and Larsen<sup>14</sup> to provide a system of estimating equations, which may lead to unbiased estimators under a generalized linear model.

In some applications, the analysis step is separated from the record linkage, for example, when the matching variables contain confidential information. This is the secondary analysis framework, in which the data analyst is only provided access to the final linked data, whereas the (unknown) record linkage process has been performed by a third-party operator.<sup>17</sup> Starting from this perspective, Chambers<sup>18</sup> proposed the exchangeable linkage error (ELE) model, and bias-corrected estimating equations for both linear and logistic regression modeling. Under the ELE model, it is assumed that linked records may be split into distinct blocks inside which the probability of correct linkage and the probability of incorrect linkage are constant. Following this work, several authors<sup>19-22</sup> developed methods for secondary analysis of linked data. Recently, Zhang and Tuoto<sup>23</sup> proposed a pseudo ordinary least square method for secondary linkage-data linear regression analysis, which can accommodate heterogeneous linkage errors and incomplete match space problems.

Although the Cox proportional hazard model<sup>24</sup> is of routine use for survival analysis, comparatively very few papers have focused on accounting for record linkage errors in this context.<sup>25</sup> performed a simulation study emphasizing the impact of missing matches on the parameter estimation of the Cox model, but did not propose any solution to obtain unbiased estimators for the model parameters. Hof et al<sup>26</sup> proposed a joint modeling for survival analysis and probabilistic record linkage. However, this analysis model is developed under a primary analysis viewpoint, while in many applications, a secondary analysis is more likely. In this work, we reason from the secondary analysis position. We propose a model to account for record linkage errors, and an estimation method to correct for the bias caused by false link errors in the Cox regression model.

The article is organised as follows. In Section 2, we propose a new estimating equation, which leads to an approximately unbiased estimation of the parameters for the Cox model with linked data. A variance estimator is also proposed. In Section 3, we evaluate the proposed estimator and the associated variance estimator through simulation studies. In Section 4, an application on a real dataset is presented. Finally, possible further research is discussed in Section 5.

## 2 | COX REGRESSION ANALYSIS WITH LINKED DATA

### 2.1 | Cox regression model

The Cox proportional hazard model<sup>24</sup> is the most popular method to assess the effect of covariates X on a survival time. This is therefore one of the most important models in medical research. Suppose that a random sample of n units is available. For each unit i = 1, ..., n, we let $\tilde{T}_i$ be a non-negative random variable, which denotes the duration between a time origin and the time of occurrence of some event of interest. We suppose that $\tilde{T}_i$ is right censored, which means that the event is observed only if it occurs before censoring time $C_i$. For units i = 1, ..., n, we therefore observe $T_i = \min(\tilde{T}_i, C_i)$. We let $\delta_i = \mathbb{1}\{\tilde{T}_i \le C_i\}$ denote the variable indicating whether the duration time is observed prior to censoring. The vector of covariates is denoted as $\mathbf{X}_i = (X_i^1, ..., X_i^p)^T$. In this section, we first suppose that $\mathbf{X}_i$ is observed for any unit in the sample. According to the Cox model, the hazard function of an event at time t is given by

$$ \lambda(t|\mathbf{X}_i) = \lambda_0(t) \exp(\mathbf{X}_i^T \beta_0), \quad (1) $$

where $\beta_0 = (\beta_{01}, ..., \beta_{0p})^T$ is a p-vector of unknown parameters and $\lambda_0(t)$ is a common baseline hazard function. Assuming that the survival times are observed on a finite interval, and that C is independent of $\tilde{T}$ conditionally on X, a consistent
estimator $\hat{\beta}$ of $\beta_0$ may be obtained by solving the estimating equation:

$$H_0(\hat{\beta}) \equiv \frac{1}{n} \sum_{i=1}^{n} \delta_i \left\{ \mathbf{X}_i - \frac{\sum_{j=1}^{n} Y_j(T_i) \exp(\mathbf{X}_j^T \hat{\beta}) \mathbf{X}_j}{\sum_{j=1}^{n} Y_j(T_i) \exp(\mathbf{X}_j^T \hat{\beta})} \right\} = 0, \quad (2)$$

where $Y_j(t) = 1_{(T_i \ge t)}$ is an at-risk indicator.²⁷ We call (2) the theoretical estimating equation. This is also the maximum partial likelihood (mpl) estimation. Under some mild assumptions, a consistent estimator of the covariance matrix of $\hat{\beta}$ is given by Reference 27

$$\hat{\mathcal{V}}_{\text{mpl}}(\hat{\beta}) = \{-n\nabla H_0(\hat{\beta})\}^{-1}. \quad (3)$$

## 2.2 | Linkage error model

Suppose that we have a dataset A of $n_A$ time-to-event data. If the covariates $\mathbf{X}_i$ were known for any unit $i \in A$, the parameter of the Cox model would be estimated by solving the theoretical estimating Equation (2). However, if the covariates are not known in database A, Equation (2) may not be solved in practice.

In order to obtain the needed covariates, a linkage is performed with a dataset B of size $n_B \ge n_A$, containing in particular the auxiliary variables $\mathbf{X}_i$. For any unit $i$ in A, we note $\mathbf{Z}_i$ for the vector of auxiliary values resulting from the linkage process. Reasoning from the secondary analysis perspective, we do not have access to the matching variables and do not know the actual linkage process.

We assume that the linkage error is non-informative of the regression model, that is, may depend on the errors in the matching process, but not on the model covariates nor on the survival time.²⁸ This is the key assumption of most secondary analysis approaches in the literature, for which Zhang and Tuoto²³ have proposed a diagnostic test. Adopting the modelling approach of Copas and Hilton,²⁹ we suppose that both databases are partitioned into blocks $A_v$ and $B_v$, $v = 1, \dots, V$, and that the record linkage is performed independently in these blocks. Also, we suppose that for any entity $i \in A_v$, we have:

$$\mathbf{Z}_i = \begin{cases} \mathbf{X}_i & \text{with probability } \alpha_v, \\ \mathbf{X}_{(j)} & \text{with probability } 1 - \alpha_v, \end{cases} \quad (4)$$

where (j) stands for some unit randomly selected in database $B_v$. In other words, it is supposed that for any $i \in A_v$, the correct entity is linked to i with probability $\alpha_v$, otherwise the unit j linked to i is randomly selected in $B_v$. The correct linkage probability of record pairs within a block can be different in practice. However, most relevant approaches for secondary analysis seem to be robust to the failure of this assumption when the linkage errors are non-informative.²³ It should be noted that we implicitly assume that A is a subset from B, and that all entities in A can therefore have some matching records in B. Also, we assume that there is at most one link for each record of both databases. In practice, there will often be some entities of A which remain unlinked after the linkage process. This may be due to errors in the matching variables, or to the fact they are not sufficiently discriminant for identifying links. Such incomplete record linkage can be problematic for further analysis if the missed links are not at random.²⁵ There are some discussions on this incomplete matching space problem.¹⁹,²³,³⁰ This problem is out of the scope of our work. We therefore assume that the linkage is complete, or alternatively that any missing links are independent on the time of event and model covariates.

## 2.3 | Adjusted estimating equation

By naively treating the linked covariates $\mathbf{Z}_i$ as if they were the true covariates $\mathbf{X}_i$ for the units $i \in A$, an estimator of $\beta_0$ may be obtained by solving the following equation:

$$H_{\text{naive}}(\hat{\beta}) \equiv \frac{1}{n_A} \sum_{v=1}^{V} \sum_{i \in A_v} \delta_i \left\{ \mathbf{Z}_i - \frac{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) \exp(\mathbf{Z}_j^T \hat{\beta}) \mathbf{Z}_j}{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) \exp(\mathbf{Z}_j^T \hat{\beta})} \right\} = 0. \quad (5)$$
We call (5) the naive estimating equation. Since some units are incorrectly linked, it may lead to biased estimates, see the simulation results in Section 3.

We propose a bias-corrected estimating equation, accounting for the fact that from the hit-miss model (4), the covariates may be incorrectly linked. We first introduce some notations. Let us define

$$g(\beta, \mathbf{X}_i) = \exp(\mathbf{X}_i^\top \beta) \text{ and } h(\beta, \mathbf{X}_i) = \exp(\mathbf{X}_i^\top \beta) \mathbf{X}_i.$$

Also, let $\bar{\mathbf{X}}_{B_v}, \bar{\mathbf{g}}_{B_v}(\beta)$ and $\bar{h}_{B_v}(\beta)$ denote the means of $\mathbf{X}_i, g(\beta, \mathbf{X}_i)$ and $h(\beta, \mathbf{X}_i)$ over $B_v$, respectively. The linkage-error adjusted estimating equation (AEE) is given by

$$\bar{H}(\beta) = \frac{1}{n_A} \sum_{v=1}^{V} \sum_{i \in A_v} \delta_i \left\{ \mathbf{X}_i^*(\alpha_v) - \frac{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) h_j^*(\alpha_v, \beta)}{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) g_j^*(\alpha_v, \beta)} \right\} = 0 \quad (6)$$

where, for any $i \in A_v$,

$$\begin{align}
\mathbf{X}_i^*(\alpha_v) &= \alpha_v^{-1} \mathbf{Z}_i - (\alpha_v^{-1} - 1) \bar{\mathbf{X}}_{B_v}, \nonumber \\
g_j^*(\alpha_v, \beta) &= \alpha_v^{-1} g(\mathbf{Z}_j, \beta) - (\alpha_v^{-1} - 1) \bar{\mathbf{g}}_{B_v}(\beta), \nonumber \\
h_j^*(\alpha_v, \beta) &= \alpha_v^{-1} h(\mathbf{Z}_j, \beta) - (\alpha_v^{-1} - 1) \bar{h}_{B_v}(\beta). \tag{7}
\end{align}$$

We prove in Appendix A that $\bar{H}(\beta)$ is an (approximately) conditionally unbiased estimator for the function $H_0(\beta)$ involved in the theoretical estimating equation. Solving the proposed AEE therefore leads to a consistent estimator of $\beta$, see the simulation results in Section 3.

Since there is no closed-form solution for the estimating equations considered above, an iterative method like the Newton-Raphson algorithm is commonly used in practice. Also, the probabilities $\alpha_v$ may be (somewhat arbitrarily) specified by the record linkage practitioner, or estimated from a validation sample$^{18,23}$ if their true values are unknown.

## 2.4 | Variance estimator

In this section, we discuss variance estimation for the estimator of the parameter $\beta_0$ obtained by solving the AEE given in (6). We first note that several sources of variance need to be accounted for: (a) the (usual) variability associated to solving a sample-based estimating equation, (b) the variability associated to the linkage process, and (c) the variability associated to the estimation of the probabilities $\alpha_v, v = 1, \dots, V$. Using the variance estimator given in (3) fails to account for all these sources of variability, and therefore leads to an underestimation of the variance, see the simulation results in Section 3.

We propose a sandwich-like variance estimator, which reads as follows:

$$\hat{\mathcal{V}}_{\text{AEE}}(\hat{\beta}) = \{\nabla \bar{H}(\hat{\beta})\}^{-1} \times \hat{\mathcal{V}}\{\bar{H}(\beta_0)\} \times \{\nabla \bar{H}(\hat{\beta})\}^{-1}, \quad (8)$$

$$\text{with } \hat{\mathcal{V}}\{\bar{H}(\beta_0)\} = \hat{\mathcal{V}}_1\{\bar{H}(\beta_0)\} + \hat{\mathcal{V}}_2\{\bar{H}(\beta_0)\}. \quad (9)$$

The first component $\hat{\mathcal{V}}_1\{\bar{H}(\beta_0)\}$ in (9) accounts for the variability in (c). Under the assumption that the validation samples $S_v$ used for such estimation are selected in the datasets $A_v$ through simple random sampling without replacement, this variance estimator is

$$\hat{\mathcal{V}}_1\{\bar{H}(\beta_0)\} = \sum_{v=1}^{V} \bar{H}_{2,v}(\hat{\alpha}_v, \hat{\beta}) \{\bar{H}_{2,v}(\hat{\alpha}_v, \hat{\beta})\}^{\top} \times \left(\frac{1}{n_{S_v}} - \frac{1}{n_{A_v}}\right) \frac{n_{S_v}}{n_{S_v} - 1} \frac{1 - \hat{\alpha}_v}{\hat{\alpha}_v^3},$$

where $n_{S_v}$ is the sample size of the validation set $S_v$, and

$$\bar{H}_{2,v}(\alpha_v, \beta) = \frac{1}{n_A} \sum_{i \in A_v} \delta_i \left\{ (\mathbf{Z}_i - \bar{\mathbf{X}}_{B_v}) - \frac{\sum_{j \in A_v} Y_j(T_i) \left\{ \{h(\beta, \mathbf{Z}_j) - \bar{h}_{B_v}(\beta)\} - R_i^*(\alpha_v, \beta) \{g(\beta, \mathbf{Z}_j) - \bar{g}_{B_v}(\beta)\} \right\}}{\sum_{j \in A_v} Y_j(T_i) g_j^*(\alpha_v, \beta)} \right\}.$$
with

$$R_i^*(\alpha_v, \beta) = \frac{\sum_{j \in A_v} Y_j(T_i) h_j^*(\alpha_v, \beta)}{\sum_{j \in A_v} Y_j(T_i) g_j^*(\alpha_v, \beta)}$$

The second component $\hat{\mathcal{V}}_2\{\overline{H}(\beta_0)\}$ in (9) accounts for both the variability in (a) and (b). We have

$$\hat{\mathcal{V}}_2\{\overline{H}(\beta_0)\} = \frac{s_H^2(\hat{\beta})}{n_A}$$

where

$$s_H^2(\beta) = \frac{1}{n_A - 1} \sum_{v=1}^{V} \sum_{i \in A_v} \left\{ H_i(\beta) - \frac{1}{n_A} \sum_{v=1}^{V} \sum_{j \in A_v} H_j(\beta) \right\}^2$$

and

$$H_i(\beta) = \delta_i \left\{ \mathbf{X}_i^*(\hat{\alpha}_v) - \frac{\sum_{v=1}^V \sum_{j \in A_v} Y_j(T_i) h_j^*(\hat{\alpha}_v, \beta)}{\sum_{v=1}^V \sum_{j \in A_v} Y_j(T_i) g_j^*(\hat{\alpha}_v, \beta)} \right\}$$

The derivation of this variance estimator is explained in detail in Appendix B. It is evaluated empirically in the next section through a simulation study.

## 3 | A SIMULATION STUDY

In this section, we evaluate the performance of the proposed estimator for the parameter of the Cox model, and the associated variance estimator. The data generation process is first presented in Section 3.1. The estimation methods that we evaluate are presented in Section 3.2, along with the performance indicators. The simulation results are given in Section 3.3. To facilitate interpretation and to study the influence of different simulation parameters, we first consider in Section 3.3.1 scenarios with a single block. Scenarios with multiple blocks and different levels of linkage quality are considered in Section 3.3.2.

### 3.1 | Data generation

Assume that there are two datasets A with $n_A$ individuals, and B with $n_B \ge n_A$ individuals. We first generate the $n_B$ units in database B with $p = 2$ covariates, including a continuous variable $X_1 \sim N(0, 1)$ and a binary variable $X_2 \sim Bernoulli(0.7)$. Given the p-vector of coefficients $\beta = (\beta_1, \beta_2)^T = (0.5, -0.5)^T$, the true survival time $\tilde{T}^B$ is generated as

$$\tilde{T}^B = - \frac{\log(U)}{\lambda \exp(\mathbf{X}^T \beta)} \qquad (10)$$

where U follows a standard uniform distribution,³¹ and $\lambda$ is fixed as equal to 1 for simplicity. A constant censoring time is chosen (from 100 000 independent data generation runs) to yield a censoring rate of approximately 0.25 over all the simulation runs.

Without loss of generality, we suppose that the units in dataset A are the $n_A$ first ones in dataset B. In other words, a pair of individuals $(a_i, b_j)$ for $i \in A$ and $j \in B$ is a match if $i = j = 1, \dots, n_A$. The survival times $T_i^A$ for $i \in A$ are therefore obtained as $T_i^A = T_i^B$ for $i = 1, \dots, n$. Given the values of $\alpha$, the linked values Z for covariates in database A are obtained according to the linkage error model (4).
If there are multiple blocks, data for each block were generated independently as follows. Firstly, for each block $v$, we generate $n_{B_v}$ observations $(T, \delta, \mathbf{X})$ from the Cox model described in Equation (10). Note that the value of the true parameters $\beta$ and the distribution of $\mathbf{X}$ are the same over blocks $v$. Then, we choose randomly $n_{A_v} \le n_{B_v}$ survival times $(T, \delta)$ for block $A_v$. All generated $n_{B_v}$ values of $\mathbf{X}$ will be placed in block $B_v$. Secondly, given the value of $\alpha_v$ for block $v$, $n_{A_v}$ linked values $\mathbf{Z}$ for block $A_v$ are obtained by the linkage error model (4). Inside each block $A_v$, an audit sample of 10% of the units is selected by simple random sampling without replacement, and the proportion of correct links in the audit sample is used as the estimator $\hat{\alpha}_v$.

## 3.2 | Methods and performance indicators

For each scenario, we consider the following estimation methods. The Theoretical is obtained by solving the theoretical estimating equation (2) with the true values of covariates $\mathbf{X}$. This is a benchmark estimation strategy, since it cannot be applied on linked data in practice. The Naive is obtained by solving the naive estimating Equation (5) with linked data. The Validation is obtained by solving the theoretical estimating Equation (2) with only correct linked pairs in the validation set. Note that, contrarily to Theoretical, this method may be used in practice if an audit sample is available. For each of these three methods, the variance of the estimator of the parameter in the Cox model is estimated by using the variance estimator $\hat{\mathbf{V}}_{\text{mpl}}(\hat{\beta})$ in Equation (3), implemented by means of R SURVIVAL package.

For each scenario, we also consider estimation methods making use of the proposed approach. The TAE (theoretical adjusted estimating equation) is obtained by solving the proposed estimating Equation (6) with the theoretical value of $\alpha_v$. The AEE (adjusted estimating equation) is obtained by solving the proposed estimating Equation (6), where $\alpha_v$ is estimated by taking the proportion of correct links in the audit sample. For each method, the Newton-Raphson algorithm is applied with a maximum of 20 iterations and an initial parameter value $\beta = (0, 0)^T$. We also report the number of time (Fails) when the Newton-Raphson algorithm does not converge. For AEE, the variance is estimated by using $\hat{\mathbf{V}}(\hat{\beta})$ in Equation (B18). For TAE, the variance is estimated by setting $\hat{\mathbf{V}}_1\{\hat{\mathbf{H}}(\hat{\beta}_0)\} = 0$ in $\hat{\mathbf{V}}(\hat{\beta})$. For both TAE and AEE, we also compare to the variance estimator $\hat{\mathbf{V}}_{\text{mpl}}(\hat{\beta})$ in Equation (3).

The data generation and the estimation process are repeated $R = 1,000$ times. Over these simulations, we compare the estimation methods in terms of the Monte Carlo bias

$$B_{\text{MC}}(\hat{\beta}) = \frac{1}{R} \sum_{r=1}^{R} (\hat{\beta}^{(r)} - \beta),$$

with $\hat{\beta}^{(r)}$ the estimator computed on the $r$-th sample. We also compute the Monte Carlo standard deviation:

$$Sd_{\text{MC}}(\hat{\beta}) = \sqrt{\frac{1}{R-1} \sum_{r=1}^{R} (\hat{\beta}^{(r)} - \bar{\hat{\beta}})^2}.$$

For the variance estimation methods, we compute the Monte Carlo estimates of standard deviation

$$\hat{Sd} = \sqrt{\frac{1}{R} \sum_{r=1}^{R} \hat{\mathbf{V}}^{(r)} (\hat{\beta}^{(r)})^2}.$$

with $\hat{\mathbf{V}}^{(r)}$ a variance estimator computed on the $r$-th sample. The Monte Carlo estimate of standard deviation is compared to the true standard deviation $Sd(\hat{\beta})$, approximated by $Sd_{\text{MC}}(\hat{\beta})$.

## 3.3 | Simulation results

### 3.3.1 | One block situation

In this section, we consider the situation when the data sets are generated as presented in Section 3.1, with $V = 1$ block only. We consider two cases. In the first one, the sample sizes $n_A = 1,000$ and $n_B = 2,000$ are held fixed, and we let the
<table><caption>TABLE 1 Simulation results in case 1 with three different values for the probability of correct link $\alpha \in \{0.75, 0.85, 0.95\}$.</caption><thead><tr><th rowspan="2">α</th><th rowspan="2">Methods</th><th rowspan="2">Fails</th><th colspan="5">β̂<sub>1</sub></th><th colspan="5">β̂<sub>2</sub></th></tr><tr><th>B<sub>MC</sub></th><th>Sd<sub>MC</sub></th><th>Sd̂<sub>mpl</sub></th><th>Sd̂<sub>AEE</sub></th><th>CP</th><th>B<sub>MC</sub></th><th>Sd<sub>MC</sub></th><th>Sd̂<sub>mpl</sub></th><th>Sd̂<sub>AEE</sub></th><th>CP</th></tr></thead><tbody><tr><td>*</td><td>Theoretical</td><td>0</td><td>0.000</td><td>0.039</td><td>0.040</td><td></td><td>0.961</td><td>0.003</td><td>0.080</td><td>0.080</td><td></td><td>0.950</td></tr><tr><td rowspan="4">0.75</td><td>Naive</td><td>0</td><td>0.147</td><td>0.041</td><td>0.039</td><td></td><td>0.050</td><td>0.143</td><td>0.081</td><td>0.081</td><td></td><td>0.577</td></tr><tr><td>Validation</td><td>0</td><td>0.017</td><td>0.160</td><td>0.156</td><td></td><td>0.941</td><td>0.003</td><td>0.318</td><td>0.302</td><td></td><td>0.936</td></tr><tr><td>TAEE</td><td>0</td><td>0.007</td><td>0.072</td><td>0.041</td><td>0.069</td><td>0.945</td><td>0.013</td><td>0.124</td><td>0.081</td><td>0.129</td><td>0.957</td></tr><tr><td>AEE</td><td>4</td><td>0.009</td><td>0.082</td><td>0.041</td><td>0.085</td><td>0.962</td><td>0.015</td><td>0.131</td><td>0.081</td><td>0.138</td><td>0.962</td></tr><tr><td rowspan="4">0.85</td><td>Naive</td><td>0</td><td>0.092</td><td>0.040</td><td>0.039</td><td></td><td>0.347</td><td>0.088</td><td>0.081</td><td>0.080</td><td></td><td>0.799</td></tr><tr><td>Validation</td><td>0</td><td>0.016</td><td>0.149</td><td>0.146</td><td></td><td>0.955</td><td>0.000</td><td>0.296</td><td>0.283</td><td></td><td>0.931</td></tr><tr><td>TAEE</td><td>0</td><td>0.002</td><td>0.055</td><td>0.041</td><td>0.059</td><td>0.964</td><td>0.007</td><td>0.103</td><td>0.080</td><td>0.113</td><td>0.969</td></tr><tr><td>AEE</td><td>0</td><td>0.005</td><td>0.063</td><td>0.041</td><td>0.066</td><td>0.969</td><td>0.010</td><td>0.110</td><td>0.080</td><td>0.118</td><td>0.972</td></tr><tr><td rowspan="4">0.95</td><td>Naive</td><td>0</td><td>0.033</td><td>0.041</td><td>0.040</td><td></td><td>0.862</td><td>0.029</td><td>0.083</td><td>0.080</td><td></td><td>0.928</td></tr><tr><td>Validation</td><td>0</td><td>0.015</td><td>0.139</td><td>0.137</td><td></td><td>0.961</td><td>0.004</td><td>0.276</td><td>0.266</td><td></td><td>0.939</td></tr><tr><td>TAEE</td><td>0</td><td>0.001</td><td>0.045</td><td>0.040</td><td>0.051</td><td>0.965</td><td>0.003</td><td>0.089</td><td>0.080</td><td>0.101</td><td>0.977</td></tr><tr><td>AEE</td><td>0</td><td>0.000</td><td>0.048</td><td>0.040</td><td>0.054</td><td>0.973</td><td>0.004</td><td>0.090</td><td>0.080</td><td>0.103</td><td>0.981</td></tr></tbody></table>

probability of correct link $\alpha$ vary in $\{0.75, 0.85, 0.95\}$. In the second one, the probability of correct link is held fixed, equal to 0.85. We let $n_A$ vary in $\{500, 1000, 2000\}$, with $n_B = 2n_A$.

The simulation results obtained in Case 1 are presented in Table 1. As expected, the Theoretical method leads to an unbiased estimation of the parameters. The Naive method leads to severely biased estimators, especially with the smaller value $\alpha = 0.75$. The bias ranges from 0.029 to 0.147, corresponding to an absolute relative bias between 5.8% and 29.0%. This bias decreases as the probability of correct link increases, as expected. The proposed methods TAEE and AEE lead to approximately unbiased estimation of the parameters, with a larger variability for AEE as expected. The bias under AEE ranges from 0.000 to 0.015, corresponding to a reduction of the relative bias (as compared to Naive) ranging between 5.0% and 27.6%. We note that the variability under both TAEE and AEE is but only moderately increased, as compared to Theoretical. The Validation method also leads to unbiased estimators of the Cox regression coefficients, but with a larger variability than both TAEE and AEE.

We now turn to the variance estimators. The variance estimator $\hat{V}_{mpl}(\hat{\beta})$ (3) performs well for Theoretical, Naive and Validation, but underestimates the variability of the estimators obtained under TAEE and AEE. This is due to the fact that this variance estimator only accounts for the variability of the sample-based estimating equation. On the other hand, the proposed variance estimator performs well, except for $\hat{\beta}_1$ when $\alpha = 0.75$, in which case the variance is underestimated. We have also computed coverage probabilities (CP) for normality-based confidence intervals with a nominal coverage of 95%. We note that the coverage probability is very poorly respected in case of Naive, even in situations when the bias is moderate.

The simulation results obtained in Case 2 are presented in Table 2. We observe no qualitative difference compared to Case 1. The TAEE and AEE lead to almost unbiased estimations for the regression coefficients, and the proposed variance estimator performs well for both methods. The bias obtained under the Naive method does not decrease as the sample size increases. As could be expected, the variability obtained under any estimation method decreases as the sample size increases.

### 3.3.2 | Multiple blocks

In this section, we consider the situation when the data sets are generated as presented in Section 3.1, with $V = 3$ blocks only. We take $(n_{A_1}, n_{A_2}, n_{A_3}) = (250, 500, 250)$ and $(n_{B_1}, n_{B_2}, n_{B_3}) = (500, 1000, 500)$. Also, we consider a first scenario where $(\alpha_1, \alpha_2, \alpha_3) = (0.6, 0.7, 0.8)$; a second scenario where $(\alpha_1, \alpha_2, \alpha_3) = (0.7, 0.8, 0.9)$; a third scenario where $(\alpha_1, \alpha_2, \alpha_3) = (0.8, 0.9, 1.0)$.
<table><caption>TABLE 2 Simulation results in case 2 with three different values for the sample size $n_A$.</caption><thead><tr><th rowspan="2">n<sub>A</sub></th><th rowspan="2">Methods</th><th rowspan="2">Fails</th><th colspan="4">β̂<sub>1</sub></th><th colspan="4">β̂<sub>2</sub></th></tr><tr><th>B<sub>MC</sub></th><th>Sd<sub>MC</sub></th><th>Sd̂<sub>mpl</sub></th><th>Sd̂<sub>AEE</sub></th><th>CP</th><th>B<sub>MC</sub></th><th>Sd<sub>MC</sub></th><th>Sd̂<sub>mpl</sub></th><th>Sd̂<sub>AEE</sub></th><th>CP</th></tr></thead><tbody><tr><td rowspan="5">500</td><td>Theoretical</td><td>0</td><td>0.002</td><td>0.056</td><td>0.057</td><td></td><td>0.954</td><td>0.005</td><td>0.113</td><td>0.114</td><td></td><td>0.955</td></tr><tr><td>Naive</td><td>0</td><td>0.089</td><td>0.057</td><td>0.056</td><td></td><td></td><td>0.087</td><td>0.113</td><td>0.114</td><td></td><td>0.876</td></tr><tr><td>Validation</td><td>0</td><td>0.033</td><td>0.222</td><td>0.215</td><td></td><td>0.951</td><td>0.024</td><td>0.435</td><td>0.419</td><td></td><td>0.949</td></tr><tr><td>TAEE</td><td>0</td><td>0.009</td><td>0.078</td><td>0.058</td><td>0.085</td><td>0.963</td><td>0.010</td><td>0.145</td><td>0.114</td><td>0.161</td><td>0.972</td></tr><tr><td>AEE</td><td>1</td><td>0.015</td><td>0.104</td><td>0.058</td><td>0.104</td><td>0.976</td><td>0.015</td><td>0.161</td><td>0.114</td><td>0.172</td><td>0.977</td></tr><tr><td rowspan="5">1000</td><td>Theoretical</td><td>0</td><td>0.000</td><td>0.039</td><td>0.040</td><td></td><td>0.961</td><td>0.003</td><td>0.080</td><td>0.080</td><td></td><td>0.950</td></tr><tr><td>Naive</td><td>0</td><td>0.092</td><td>0.040</td><td>0.039</td><td></td><td>0.347</td><td>0.088</td><td>0.081</td><td>0.080</td><td></td><td>0.799</td></tr><tr><td>Validation</td><td>0</td><td>0.016</td><td>0.149</td><td>0.146</td><td></td><td>0.955</td><td>0.000</td><td>0.296</td><td>0.283</td><td></td><td>0.931</td></tr><tr><td>TAEE</td><td>0</td><td>0.002</td><td>0.055</td><td>0.041</td><td>0.059</td><td>0.964</td><td>0.007</td><td>0.103</td><td>0.080</td><td>0.113</td><td>0.969</td></tr><tr><td>AEE</td><td>0</td><td>0.005</td><td>0.063</td><td>0.041</td><td>0.066</td><td>0.969</td><td>0.010</td><td>0.110</td><td>0.080</td><td>0.118</td><td>0.972</td></tr><tr><td rowspan="5">2000</td><td>Theoretical</td><td>0</td><td>0.000</td><td>0.028</td><td>0.028</td><td></td><td>0.945</td><td>0.000</td><td>0.056</td><td>0.057</td><td></td><td>0.960</td></tr><tr><td>Naive</td><td>0</td><td>0.092</td><td>0.029</td><td>0.028</td><td></td><td>0.111</td><td>0.092</td><td>0.056</td><td>0.057</td><td></td><td>0.640</td></tr><tr><td>Validation</td><td>0</td><td>0.006</td><td>0.103</td><td>0.100</td><td></td><td>0.932</td><td>0.003</td><td>0.197</td><td>0.197</td><td></td><td>0.948</td></tr><tr><td>TAEE</td><td>0</td><td>0.001</td><td>0.039</td><td>0.029</td><td>0.041</td><td>0.953</td><td>0.000</td><td>0.071</td><td>0.057</td><td>0.080</td><td>0.971</td></tr><tr><td>AEE</td><td>0</td><td>0.002</td><td>0.043</td><td>0.029</td><td>0.046</td><td>0.964</td><td>0.001</td><td>0.075</td><td>0.057</td><td>0.082</td><td>0.969</td></tr></tbody></table>

Let $\bar{\alpha}$ be the weighted average of $\alpha_1, \dots, \alpha_v$ defined as

$$ \bar{\alpha} = \frac{\sum_{i=1}^{V} n_{A_v} \alpha_v}{\sum_{i=1}^{V} n_{A_v}}. $$

This leads to a percentage of correct links approximately equal to $\bar{\alpha} = 70\%$ in Scenario 1, $\bar{\alpha} = 80\%$ in Scenario 2 and $\bar{\alpha} = 90\%$ in Scenario 3. In this context, we also consider two additional versions of our proposed methods, when we are unable to access to the value $\alpha_v$ of each block, but we have only access to their weighted average: TAEE-$\bar{\alpha}$ where the TAEE is used with $V = 1$ and true value of $\bar{\alpha}$, and AEE-$\bar{\alpha}$ where the AEE is used with $V = 1$ and estimated value of $\hat{\bar{\alpha}}$.

The simulation results are presented in Table 3, and confirm the good results of the proposed methods observed in the situation of one block. Scenarios 1 and 2 are the cases when the behaviour of the Naive method is particularly poor, with a very large bias due to a larger number of false links, and very poor coverage for the confidence intervals. On the other hand, AEE performs well in reducing the estimation bias even in this situation. The proposed variance estimator also performs well in these cases. The standard errors of TAEE and AEE estimators decrease as $\bar{\alpha}$ increase, that is, by going from Scenario 1 to Scenario 3 in Table 3. As explained in Section 2.4, there are three sources of variance in the estimation process: (a) the variability associated to solving a sample-based estimating equation, (b) the variability associated to the linkage process, and (c) the variability associated to the estimation of the probabilities $\alpha_v$. Since the sample size is kept constant, the term (a) is likely not affected by the value of $\alpha_v$. The term (b) decreases as $\alpha_v$ increases, as the variance in the hit-miss model (4) does so. The term (c) also decreases as $\alpha_v$ increases, as illustrated by the fact that $\hat{V}_1$ depends on $(1 - \alpha)/\alpha^3$, which is decreasing as $\alpha \to 1$. Concerning the coverage probability of normality-based confidence intervals, we note that they are well respected under the proposed methods, although the confidence intervals are slightly conservative when the variance estimators are so.

When the block-specific true link rate is not correlated with the block-specific distribution of $T$ and $X$, for example, this multiple blocks simulation set up, a single-$\bar{\alpha}$ adjustment (TAEE-$\bar{\alpha}$ and AEE-$\bar{\alpha}$) can still perform well. The main result in Table 3 concerning AEE and AEE-$\bar{\alpha}$ is that they both lead to virtually unbiased estimators. The bias is indeed always smaller with AEE-$\bar{\alpha}$, but the difference is no greater than 0.009, which is very small as compared to the value of the parameters ($\beta_1 = 0.5$ and $\beta_2 = -0.5$). Closeness between TAEE and TAEE-$\bar{\alpha}$ confirms the somewhat favourable simulation setup of non-informative linkage error. Reduced bias of AEE-$\bar{\alpha}$ compared to AEE-$\bar{\alpha}$ may be due to the non-linearity of
<table><caption>TABLE 3 Simulation results with three blocks with different linkage quality.</caption><thead><tr><th rowspan="2">Scenario</th><th rowspan="2">Methods</th><th rowspan="2">Fails</th><th colspan="5">β̂<sub>1</sub></th><th colspan="5">β̂<sub>2</sub></th></tr><tr><th>B<sub>MC</sub></th><th>Sd<sub>MC</sub></th><th>Sd<sub>mpl</sub></th><th>Sd<sub>AEE</sub></th><th>CP</th><th>B<sub>MC</sub></th><th>Sd<sub>MC</sub></th><th>Sd<sub>mpl</sub></th><th>Sd<sub>AEE</sub></th><th>CP</th></tr></thead><tbody><tr><td>*</td><td>Theoretical</td><td>0</td><td>0.002</td><td>0.040</td><td>0.039</td><td></td><td>0.953</td><td>0.002</td><td>0.078</td><td>0.078</td><td></td><td>0.944</td></tr><tr><td rowspan="6">1</td><td>Naive</td><td>0</td><td>0.171</td><td>0.041</td><td>0.039</td><td></td><td>0.010</td><td>0.171</td><td>0.082</td><td>0.081</td><td></td><td>0.440</td></tr><tr><td>Validation</td><td>0</td><td>0.021</td><td>0.161</td><td>0.161</td><td></td><td>0.961</td><td>0.002</td><td>0.322</td><td>0.315</td><td></td><td>0.945</td></tr><tr><td>TAEE</td><td>1</td><td>0.018</td><td>0.097</td><td>0.042</td><td>0.136</td><td>0.944</td><td>0.013</td><td>0.143</td><td>0.081</td><td>0.144</td><td>0.947</td></tr><tr><td>AEE</td><td>17</td><td>0.030</td><td>0.136</td><td>0.042</td><td>0.128</td><td>0.961</td><td>0.022</td><td>0.177</td><td>0.081</td><td>0.183</td><td>0.960</td></tr><tr><td>TAEE-ᾱ</td><td>0</td><td>0.017</td><td>0.086</td><td>0.043</td><td>0.139</td><td>0.943</td><td>0.013</td><td>0.139</td><td>0.081</td><td>0.141</td><td>0.951</td></tr><tr><td>AEE-ᾱ</td><td>9</td><td>0.021</td><td>0.108</td><td>0.042</td><td>0.144</td><td>0.964</td><td>0.016</td><td>0.153</td><td>0.081</td><td>0.168</td><td>0.963</td></tr><tr><td rowspan="6">2</td><td>Naive</td><td>0</td><td>0.118</td><td>0.041</td><td>0.039</td><td></td><td>0.167</td><td>0.120</td><td>0.084</td><td>0.080</td><td></td><td>0.660</td></tr><tr><td>Validation</td><td>0</td><td>0.018</td><td>0.151</td><td>0.150</td><td></td><td>0.948</td><td>0.002</td><td>0.294</td><td>0.293</td><td></td><td>0.946</td></tr><tr><td>TAEE</td><td>0</td><td>0.007</td><td>0.066</td><td>0.041</td><td>0.064</td><td>0.952</td><td>0.003</td><td>0.118</td><td>0.080</td><td>0.122</td><td>0.953</td></tr><tr><td>AEE</td><td>1</td><td>0.015</td><td>0.086</td><td>0.041</td><td>0.081</td><td>0.969</td><td>0.010</td><td>0.129</td><td>0.080</td><td>0.135</td><td>0.961</td></tr><tr><td>TAEE-ᾱ</td><td>0</td><td>0.007</td><td>0.064</td><td>0.041</td><td>0.063</td><td>0.955</td><td>0.003</td><td>0.116</td><td>0.080</td><td>0.120</td><td>0.959</td></tr><tr><td>AEE-ᾱ</td><td>0</td><td>0.009</td><td>0.073</td><td>0.041</td><td>0.075</td><td>0.966</td><td>0.006</td><td>0.123</td><td>0.080</td><td>0.127</td><td>0.963</td></tr><tr><td rowspan="6">3</td><td>Naive</td><td>0</td><td>0.060</td><td>0.041</td><td>0.040</td><td></td><td>0.662</td><td>0.061</td><td>0.082</td><td>0.080</td><td></td><td>0.882</td></tr><tr><td>Validation</td><td>0</td><td>0.016</td><td>0.143</td><td>0.140</td><td></td><td>0.945</td><td>0.006</td><td>0.272</td><td>0.275</td><td></td><td>0.950</td></tr><tr><td>TAEE</td><td>0</td><td>0.005</td><td>0.052</td><td>0.041</td><td>0.056</td><td>0.965</td><td>0.004</td><td>0.097</td><td>0.080</td><td>0.108</td><td>0.967</td></tr><tr><td>AEE</td><td>0</td><td>0.007</td><td>0.058</td><td>0.041</td><td>0.062</td><td>0.973</td><td>0.006</td><td>0.102</td><td>0.080</td><td>0.112</td><td>0.971</td></tr><tr><td>TAEE-ᾱ</td><td>0</td><td>0.004</td><td>0.051</td><td>0.041</td><td>0.055</td><td>0.965</td><td>0.004</td><td>0.096</td><td>0.080</td><td>0.107</td><td>0.972</td></tr><tr><td>AEE-ᾱ</td><td>0</td><td>0.005</td><td>0.055</td><td>0.041</td><td>0.060</td><td>0.970</td><td>0.005</td><td>0.099</td><td>0.080</td><td>0.109</td><td>0.973</td></tr></tbody></table>

adjustment, such that the additional variance of AEE-α adjustment is manifested in terms of the bias of adjustment. Moreover, a single-ᾱ adjustment can provide a smaller variance. In practice, this is very helpful when the analyst cannot conduct auditing, and when the linker can only provide a single overall estimate of α. However, the linkage error may be informative, such as when β and α vary across the blocks in a correlated manner. Block-specific adjustment would then be clearly more helpful at reducing the bias than adjustment by a single ᾱ.

Some additional simulation results are presented in the supplementary material. In particular, we have studied the situation when the non-informative assumption is not true. The simulation results in Tables S1 and S2 indicate that the Cox parameters estimated under TAEE and AEE will be more biased when α is dependent on variables from the Cox model. This is especially true when α is small and dependent on T̄ (see Table S2). This emphasizes the importance of the non-informative linkage error assumption. We note, however, that the proposed methods still perform better in this case than the Naive method. We have also performed a sensitivity analysis, evaluating the performance of TAEE with incorrect values for the parameter α. The results are presented in Table S3. As could be expected, the bias in the estimated parameters increases with the error in α, but the estimator remains less biased than with the Naive method if the error is moderate.

# 4 | APPLICATION

## 4.1 | Data description

The proposed model is fitted to a linked dataset between a registry of strokes, denoted by AVC ("Accident Vasculaire Cérébral"), and an extraction of the national health information system of France, denoted by SNDS ("Système national des données de santé"). The AVC recorded all stroke cases of patients aged 15 years and older, who have lived in the Brest
<table><caption>TABLE 4 Description of the linkage process.</caption><thead><tr><th>Steps</th><th>Number of agreements among nine matching variables</th><th>Number of record pairs</th></tr></thead><tbody><tr><td>1</td><td>9</td><td>1,792</td></tr><tr><td>2</td><td>8</td><td>170</td></tr><tr><td>3</td><td>7</td><td>11</td></tr><tr><td>4</td><td>6</td><td>1,500</td></tr><tr><td>5</td><td>5</td><td>58</td></tr><tr><td>6</td><td>4</td><td>4</td></tr><tr><td>Total</td><td></td><td>3535</td></tr></tbody></table>

<table><caption>TABLE 5 Description of the linked database.</caption><thead><tr><th>Variable</th><th>Description</th><th>Source</th></tr></thead><tbody><tr><td>Time</td><td>Time (in days) between the first stroke and death or end of follow-up (31/12/2018)</td><td>AVC</td></tr><tr><td>Censoring</td><td>If the patient died before 01/01/2019: 1 = Yes, 0 = No</td><td>AVC</td></tr><tr><td>Age</td><td>Age (in years) at the first stroke</td><td>SNDS</td></tr><tr><td>Gender</td><td>Sex: 0 = Male, 1 = Female</td><td>SNDS</td></tr><tr><td>Type AVC</td><td>Type of stroke (0 = Ischemic, 1 = Hemorrhagic)</td><td>SNDS</td></tr></tbody></table>

area from 2008 to the end of 2018. SNDS is an extraction from the French health information system, and contains patients for whom at least one medical service or hospitalization were recorded since 2008 while they were living in the Brest area. Due to the limited information in the registry, there is a demand of linking AVC and SNDS to enrich the registry for further analyses.

The linkage was performed by a separate team, and due to confidentiality restrictions, we were not allowed to access to the matching data and have limited knowledge about the linkage. A deterministic record linkage method was used. This is the simpler linkage approach, which ideally requires agreement on all matching variables, or otherwise on a (large) subset of these variables. In the linkage process, there are nine matching variables, and the linkage is implemented sequentially. In the first step, it is required that the nine matching variables agree for a pair to be viewed as a link. The corresponding pairs are then suppressed, and among the remaining ones it is asked that eight matching variables agree for a pair to be viewed as a link. The procedure continues on similarly. The process is summarized in Table 4.

After performing the linkage process, a dataset of 3,535 patients has been obtained. It contains the survival time, the censoring indicator and three covariates (age, gender, type of stroke). We suppose that these covariates were obtained from SNDS by the linkage process, and may therefore be affected by linkage errors. A description of the dataset is presented in Table 5. In this application, we are interested in comparing the risk of death after the first stroke between males and females, taking into account the age and the type of stroke.

## 4.2 | Cox regression analysis

In this application, we use the Cox regression model (1) to model the relationship between the survival time and three explanatory variables (age, gender, type of stroke). We consider AVC as database A and SNDS as database B in our proposed model. In the naive approach, we use the linked data as if it was directly observed. However, the simulation results in Section 3.3 show that linkage errors lead to biased estimators of the regression coefficients. Therefore, we also use the adjusted estimating Equation (6).

For the record pairs obtained at each step, the percentage of matching variables which are in agreement are seen as a proxy of the probability that the matching is correct. For example, for the 1,500 pairs obtained at step 4, the probability that the matching is correct is estimated as 6/9 = 0.667. We suppose that the linked dataset is comprised of two blocks, so
<table><caption>TABLE 6 Estimated coefficients (coef), estimated standard deviation of the estimated coefficients (SD), and the hazard ratio (hr = exp(coef)) of the naive method and the AEE method from linked data.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Naive method</th><th colspan="3">AEE</th></tr><tr><th>coef</th><th>SD</th><th>hr</th><th>coef</th><th>SD</th><th>hr</th></tr></thead><tbody><tr><td>Age</td><td>0.059</td><td>0.002</td><td>1.061</td><td>0.070</td><td>0.001</td><td>1.073</td></tr><tr><td>Sex</td><td>-0.120</td><td>0.047</td><td>0.887</td><td>-0.145</td><td>0.067</td><td>0.865</td></tr><tr><td>Type AVC</td><td>0.773</td><td>0.058</td><td>2.165</td><td>0.846</td><td>0.082</td><td>2.330</td></tr></tbody></table>

as to avoid the possibility of dependency between the linkage process performed into the different blocks. The estimates of $\alpha_v$ for each block v are obtained as follows:

* Block 1: 1,792 record pairs are obtained from Step 1, with $\hat{\alpha}_1 = 9/9 = 1$.

* Block 2: 1,743 remaining record pairs, with $\hat{\alpha}_2 = \frac{170 \times 8/9 + 11 \times 7/9 + 1500 \times 6/9 + 58 \times 5/9 + 4 \times 4/9}{1743} \approx 0.694$.

Besides, because the covariates are not available for any units in the SNDS, the adjustment terms in (7) cannot be computed since the proposed approach requires full access to the set of covariates in database B. We therefore use the proxy solution suggested in Equation (C1), which requires that the covariates are known on the linked dataset only. Simulations in Appendix C show that if the database A may be seen as a random sample from the database B, or when the sampling leading to A is independent of the covariates, this method leads to comparable results as the method proposed in Section 2.3.

In Table 6, we present the estimations arising from both the Naive and the AEE methods. The two methods decidedly lead to different estimations. If the Naive method is used, the hazard ratio of sex is 0.887, which means that given the same age and the same type of stroke, the female's risk of death after the first stroke is 0.887 times smaller than male's. On one hand, this ratio from the adjusted estimating equation approach is just 0.865.

## 5 | DISCUSSION

In this work, our simulations proved that the naive use of linked data may lead to substantial bias in a Cox regression model. Therefore, under the secondary analysis position where the analyst can access to linked data only, we have proposed an adjusted estimating equation for linked data, which can correct the bias from the naive estimating equation. A variance estimator, which can capture three sources of variability has also been proposed. However, proving the asymptotic normality of the resulting estimators remains challenging.

Through various simulation scenarios with one block and also multiple blocks, the proposed adjusted estimating equation is shown to lead to substantial bias reductions as compared to the naive estimating equation. Additional simulations study the non-information linkage assumption and the sensitivity analysis of $\hat{\alpha}$ are also presented in the Supplementary material.

In our modeling, it is assumed that the probability of correct linkage is identical within blocks. The assumption may not be completely realistic, since the linkage errors vary over different individuals. However, the detrimental effect may be limited as long as the linkage errors are non-informative, as shown both by simulation and with real-life linkage data in Reference 23. It would be desirable to study modelling approaches that allow for individual linkage errors, as well as to carry out empirical validations of adjustment methods derived under the simplified assumption of linkage errors. A key to both would be the possibility to work with the data linkers directly on such problems.

We have also proposed different variants of the approach for scenarios where information is limited. For example, when the block-specific linkage rate $\alpha_v$ is not available for each block, our method still works well by using the average true link rate $\bar{\alpha}$. If the analysts are not able to fully access the covariates in database B, we proposed to use the adjustments in (C1) in the Appendix, which still maintain the good performance of the AEE if A is a random sample from B. Detailed simulation results are presented in Table S4 of the Supplementary material. In addition, a linear approximated estimating equation (LAEE), which can provide better estimation than AEE with small sample size, is given in Table S5 of the Supplementary material.
Although the proposed method has improved on the naive estimation, there are perspectives that need to be developed. In this work, we assumed that observations on survival time are already available and all explanatory variables are obtained from another database. In practice, there are some cases when a part of the covariates is also available in A, and only a part of the covariates is acquired from B by linkage. In addition, the covariates can be obtained from several sources with different linkage processes. The proposed model should be developed to adapt to these cases.

We also supposed that the survival time and the censoring indicator are observed in database A, while the explanatory variables are obtained from database B by a linkage process. However, the opposite situation may occur in practice: the covariates may be available for the units in A, while the survival time needs to be obtained from another database B by a linkage process. The proposed adjustment in Equation (6) only accounts for the error associated to $Z_i$. If $T_i$ and $\delta_i$ are linked from dataset B, they are prone to linkage errors which need to be accounted for in modifying the estimating equation. This requires a different adjustment approach.

## DATA AVAILABILITY STATEMENT

Our R programs for simulation results are available at Github, https://github.com/thanhhuanVO/Cox-regression-with-linked-data. The data used in the application of the methods may be obtained from a third party and are not publicly available. For all interested researchers, data are available via SNDS, https://www.snds.gouv.fr/SNDS/Accueil, subject to the authorization of CNIL (National Commission on Informatics and Liberty of France). The data that support the findings of this study are available on request from the corresponding author. The data are not publicly available due to privacy or ethical restrictions.

## ORCID

Thanh Huan Vo https://orcid.org/0000-0001-8764-9256

## REFERENCES

1. Harron K, Gilbert R, Cromwell D, Meulen J. Linking data for mothers and babies in de-identified electronic health data. *PloS One*. 2016;11:1-18.

2. Padmanabhan S, Carty L, Cameron E, Ghosh RE, Williams R, Strongman H. Approach to record linkage of primary care data from clinical practice research datalink to other health-related patient data: overview and implications. *Eur J Epidemiol*. 2018;34:91-99.

3. Zhang G, Campbell P. Data survey: developing the statistical longitudinal census dataset and identifying its potential uses. *Aust Econ Rev*. 2012;45:125-133.

4. Winkler W, Thibaudeau Y. *An Application of the Fellegi-Sunter Model of Record Linkage to the 1990 U.S.* Decennial Census in Technical report. Washington, DC: US Bureau of the Census; 1987.

5. Herzog T, Scheuren F, Winkler W. *Data Quality and Record Linkage Techniques*. New York: Springer-Verlag; 2007.

6. Christen P. *Data Matching: Concepts and Techniques for Record Linkage, Entity Resolution, and Duplicate Detection*. Berlin, Heidelberg: Springer Publishing Company; 2012.

7. Ivan F, Alan S. A theory for record linkage. *J Am Stat Assoc*. 1969;64:1183-1210.

8. Winkler WE. Using the EM algorithm for weight computation in the Fellegi-Sunter model of record linkage. Proceedings of the Section on Survey Research Methods, American Statistical Association. 1988 667-671.

9. Huan VT, Guillaume C, André H, Emmanuel O, Stéphane P, Valérie G. Extending the Fellegi-Sunter record linkage model for mixed-type data with application to the French national health data system. *Comput Stat Data Anal*. 2023;179:107656.

10. Tancredi A, Liseo B. A hierarchical Bayesian approach to matching and size population problems. *Ann Appl Stat*. 2011;5(2B):1553-1585.

11. Sadinle M. Bayesian estimation of bipartite matchings for record linkage. *J Am Stat Assoc*. 2017;112:600-612.

12. Neter J, Maynes ES, Ramanathan R. The effect of mismatching on the measurement of response errors. *J Am Stat Assoc*. 1965;60:1005-1027.

13. Scheuren F, Winkler W. Regression analysis of data files that are computer matched. *Surv Methodol*. 1993;19:39-58.

14. Lahiri P, Larsen MD. Regression analysis with linked data. *J Am Stat Assoc*. 2005;100:222-230.

15. Hof M, Zwinderman A. Methods for analyzing data from probabilistic linkage strategies based on partially identifying variables. *Stat Med*. 2012;31:4231-4242.

16. Han Y, Lahiri P. Statistical analysis with linked data. *Int Stat Rev*. 2019;87:S139-S157.

17. Zhang L-C. On secondary analysis of datasets that cannot be linked without errors. *Analysis of Integrated Data*. London: CRC/Chapman and Hall; 2019.

18. Chambers R. Regression analysis of probability-linked data. *Off. Statist. Res. Ser.*, 4. Statistics New Zealand, Wellington; 2009. Available from http://www3.stats.govt.nz/statisphere/official Statistics Research Series/Regression Analysis of Probability-Linked Data.pdf.

19. Kim G, Chambers R. Regression analysis under probabilistic multi-linkage. *Stat Neerland*. 2012;66:64-79.

20. Kim G, Chambers R. Regression analysis under incomplete linkage. *Comput Stat Data Anal*. 2012;56:2756-2770.

21. Chambers RC, Kim G. Secondary analysis of linked data. In: Harron K, Goldstein H, Dibben C, eds. *Methodological developments in data linkage, Chapter 5*. Hoboken: John Wiley & Sons; 2015:83-108.

22. Chambers R, Diniz da Silva A. Improved secondary analysis of linked data: a framework and an illustration. *J R Stat Soc Ser A Stat Soc*. 2020;183(1):37-59.
23. Zhang L-C, Tuoto T. Linkage-data linear regression. *J R Stat Soc A Stat Soc.* 2020;100:222-230.

24. Cox DR. Regression models and life-tables. *J R Stat Soc Ser B (Methodol)*. 1972;34:187-220.

25. Baldi I, Ponti A, Zanetti R, Ciccone G, Merletti F, Gregori D. The impact of record inkage bias in the Cox model. *J Eval Clin Pract.* 2010;16:92-96.

26. Hof M, Ravelli A, Zwinderman A. A probabilistic record linkage model for survival data. *J Am Stat Assoc.* 2017;112:1504-1515.

27. Andersen PK, Gill RD. Cox's regression model for counting processes: a large sample study. *Ann Stat.* 1982;10:1100-1120.

28. Chambers R, Salvati N, Fabrizi E, Diniz SA. Domain estimation under informative linkage. *Stat Theory Related Fields*. 2019;3:90-102.

29. Copas JB, Hilton FJ. Record linkage: statistical models for matching computer records. *J R Stat Soc Ser A Stat Soc.* 1990;153:287-320.

30. Goldstein H, Harron K, Wade A. The analysis of record-linked data using multiple imputation with data value priors. *Stat Med.* 2012;31:3481-3493.

31. Bender R, Augustin T, Blettner M. Generating survival times to simulate Cox proportional hazards models. *Stat Med.* 2005;24:1713-1723.

## SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.

**How to cite this article:** Vo TH, Garès V, Zhang L-C, et al. Cox regression with linked data. *Statistics in Medicine.* 2024;43(2):296-314. doi: 10.1002/sim.9960

## APPENDIX A. EXPECTATION OF THE ADJUSTED ESTIMATING EQUATION

The proposed adjusted estimating equation is given by

$$ \overline{H}(\beta) = \frac{1}{n_A} \sum_{v=1}^{V} \sum_{i \in A_v} \delta_i \left\{ \mathbf{X}_i^*(\alpha_v) - \frac{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) h_j^*(\alpha_v, \beta)}{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) g_j^*(\alpha_v, \beta)} \right\} = 0. \quad (A1) $$

Let $\mathcal{F} = \{(T_i, \delta_i), i = 1, \dots, n_A \text{ and } \mathbf{X}_j, j = 1, \dots, n_B\}$ denote the information related to the duration times and censoring indicators for the units in A, and to the true values of covariates for all the units in B. We have

$$ \begin{aligned} \mathbb{E}\{\overline{H}(\beta)|\mathcal{F}\} &= \frac{1}{n_A} \sum_{v=1}^{V} \sum_{i \in A_v} \mathbb{E}\left\{\delta_i \left[ \mathbf{X}_i^* - \frac{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) h_j^*(\alpha_v, \beta)}{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) g_j^*(\alpha_v, \beta)} \right] \middle| \mathcal{F} \right\} \\ &= \frac{1}{n_A} \sum_{v=1}^{V} \sum_{i \in A_v} \delta_i \left\{ \underbrace{\mathbb{E}(\mathbf{X}_i^*|\mathcal{F})}_{E_1} - \underbrace{\mathbb{E}\left(\frac{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) h_j^*(\alpha_v, \beta)}{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) g_j^*(\alpha_v, \beta)} \middle| \mathcal{F} \right)}_{E_2} \right\} \end{aligned} \quad (A2) $$

For each $i \in A_v$ and $j \in B_v$, let $l_{ij}$ be an indicator equal to 1 if unit i and j are linked, and to 0 otherwise. Then for each $i \in A_v$, we have $\mathbf{Z}_i = \sum_{j \in B_v} l_{ij} \mathbf{X}_j$, and

$$ \mathbb{E}(\mathbf{Z}_i|\mathcal{F}) = \sum_{j \in B_v} \mathbf{X}_j \mathbb{E}(l_{ij}|\mathcal{F}). $$

Under the non-informative assumption for the linkage process, we obtain from the hit-miss model (4) that

$$ \mathbb{E}(l_{ii}|\mathcal{F}) = \alpha_v + (1 - \alpha_v)(n_B)^{-1}, $$

$$ \mathbb{E}(l_{ij}|\mathcal{F}) = (1 - \alpha_v)(n_B)^{-1} \text{ for } j \in B \setminus \{i\}, $$
which leads to

$$ \mathbb{E}(\mathbf{Z}_i | \mathcal{F}) = \alpha_v \mathbf{X}_i + (1 - \alpha_v) \bar{\mathbf{X}}_{B_v} $$

From Equation (7) and under the non-informative linkage assumption, we have

$$ 
\begin{aligned}
E_1 &= \mathbb{E}\{\alpha_v^{-1}\mathbf{Z}_i - (\alpha_v^{-1} - 1)\bar{\mathbf{X}}_{B_v} | \mathcal{F}\} \\
&= \alpha_v^{-1}\mathbb{E}(\mathbf{Z}_i | \mathcal{F}) - (\alpha_v^{-1} - 1)\bar{\mathbf{X}}_{B_v} \\
&= \alpha_v^{-1}[\alpha_v\mathbf{X}_i + (1 - \alpha_v)\bar{\mathbf{X}}_{B_v}] - (\alpha_v^{-1} - 1)\bar{\mathbf{X}}_{B_v} \\
&= \mathbf{X}_i.
\end{aligned}
\quad (A3) $$

By using a first order Taylor approximation, we have up to negligible factors of order $O_p(n_A^{-1})$:

$$ E_2 \approx \frac{\mathbb{E}\left\{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) h_j^*(\alpha_v, \boldsymbol{\beta}) \mid \mathcal{F}\right\}}{\mathbb{E}\left\{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) g_j^*(\alpha_v, \boldsymbol{\beta}) \mid \mathcal{F}\right\}} \quad (A4) $$

where

$$ 
\begin{aligned}
\mathbb{E}\left\{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) h_j^*(\alpha_v, \boldsymbol{\beta}) \mid \mathcal{F}\right\} &= \sum_{v=1}^{V} \sum_{j \in A_v} \mathbb{E}\left\{Y_j(T_i) h_j^*(\alpha_v, \boldsymbol{\beta}) \mid \mathcal{F}\right\} \\
&= \sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) \mathbb{E}\left\{h_j^*(\alpha_v, \boldsymbol{\beta}) \mid \mathcal{F}\right\} \\
&= \sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) h(\boldsymbol{\beta}, \mathbf{X}_j).
\end{aligned}
$$

Similarly:

$$ \mathbb{E}\left(\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) g_j^*(\alpha_v, \boldsymbol{\beta}) \mid \mathcal{F}\right) = \sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) g(\boldsymbol{\beta}, \mathbf{X}_j). $$

Therefore,

$$ E_2 \approx \frac{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) h(\boldsymbol{\beta}, \mathbf{X}_j)}{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) g(\boldsymbol{\beta}, \mathbf{X}_j)} \quad (A5) $$

By plugging (A3) and (A5) into (A2), we obtain

$$ 
\begin{aligned}
\mathbb{E}\{\bar{H}(\boldsymbol{\beta}) | \mathcal{F}\} &\approx \frac{1}{n_A} \sum_{v=1}^{V} \sum_{i \in A_v} \delta_i \left\{ \mathbf{X}_i - \frac{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) h(\boldsymbol{\beta}, \mathbf{X}_j)}{\sum_{v=1}^{V} \sum_{j \in A_v} Y_j(T_i) g(\boldsymbol{\beta}, \mathbf{X}_j)} \right\} \\
&= H_0(\boldsymbol{\beta}).
\end{aligned}
\quad (A6) $$

## APPENDIX B. VARIANCE ESTIMATION FOR THE PROPOSED ADJUSTED ESTIMATOR

In this appendix, the derivation of the variance estimator is explained. For simplicity, we focus on the case $V = 1$ when a single block is used. The extension to multiple blocks is straightforward.
We first recall the main notations. A database B of size n<sub>B</sub> is first obtained, and the covariates X<sub>i</sub> are observed for all the units in B. We use the notations

$$
\bar{\mathbf{X}}_B = \frac{1}{n_B} \sum_{i=1}^{n_B} \mathbf{X}_i,
$$

$$
\bar{g}_B(\boldsymbol{\beta}) = \frac{1}{n_B} \sum_{i=1}^{n_B} g(\boldsymbol{\beta}, \mathbf{X}_i),
$$

$$
\bar{h}_B(\boldsymbol{\beta}) = \frac{1}{n_B} \sum_{i=1}^{n_B} h(\boldsymbol{\beta}, \mathbf{X}_i).
$$

We also note $\mathbf{X}_B \equiv \{\mathbf{X}_i\}_{i \in B}$ for the set of auxiliary variables in B.

A subsample A of size n_A is then selected in B, and the variable T_i is obtained for any unit i ∈ A. We note T_A ≡ {T_i}i∈A for the set of outcome values in A. The auxiliary variables are obtained in A by using record linkage, leading to the pseudo auxiliary variables Z_i for any unit i ∈ A. We note Z_A ≡ {Z_i}i∈A for the set of pseudo values in A.

Finally, a validation sample V of size n<sub>V</sub> is selected in A by simple random sampling, and the true auxiliary variables
X<sub>i</sub> are obtained for the units i ∈ V. By comparing the pseudo values Z<sub>i</sub> and the true values X<sub>i</sub> in V, we obtain an unbiased
estimator ˆα for the parameter α.

**B.1 Global estimating equation**

Using the unbiased estimator $\hat{\alpha}$ for the parameter $\alpha$ (see Equation 4), the global estimating equation for the parameter $\beta$ is

$$
\bar{H}(\boldsymbol{\beta}) \equiv \frac{1}{n_A} \sum_{i=1}^{n_A} \underbrace{\delta_i \left\{ \mathbf{X}_i^*(\hat{\alpha}) - \frac{\sum_{j=1}^{n_A} Y_j(T_i) h_j^*(\hat{\alpha}, \boldsymbol{\beta})}{\sum_{j=1}^{n_A} Y_j(T_i) g_j^*(\hat{\alpha}, \boldsymbol{\beta})} \right\}}_{H_i(\boldsymbol{\beta})} = 0, \quad (B1)
$$

where

$$
\mathbf{X}_i^*(\hat{\alpha}) = \frac{\mathbf{Z}_i}{\hat{\alpha}} - \frac{1 - \hat{\alpha}}{\hat{\alpha}} \bar{\mathbf{X}}_B,
$$

$$
g_j^*(\hat{\alpha}, \beta) = \frac{g(\beta, \mathbf{Z}_j)}{\hat{\alpha}} - \frac{1 - \hat{\alpha}}{\hat{\alpha}} \bar{g}_B(\beta),
$$

$$
h_j^*(\hat{\alpha}, \beta) = \frac{h(\beta, \mathbf{Z}_i)}{\hat{\alpha}} - \frac{1 - \hat{\alpha}}{\hat{\alpha}} \bar{h}_B(\beta). \quad (B2)
$$

Let us denote by $\boldsymbol{\beta}_0$ the true value of the parameter. Then we have

$$
\overline{H}(\hat{\boldsymbol{\beta}}) - \overline{H}(\boldsymbol{\beta}_0) = -\overline{H}(\boldsymbol{\beta}_0) \approx \{\mathbb{E}\nabla\overline{H}(\boldsymbol{\beta}_0)\}\{\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}_0\},
$$

with $\nabla \overline{H}(\boldsymbol{\beta})$ the differential of $\overline{H}(\boldsymbol{\beta})$. We obtain

$$
\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}_0 \approx -\{\mathbb{E}\nabla\overline{H}(\boldsymbol{\beta}_0)\}^{-1} \times \overline{H}(\boldsymbol{\beta}_0).
$$

It is thus sufficient to obtain a variance estimator for $\overline{H}(\beta_0)$, from which we can use the sandwich variance estimator

$$
\hat{\mathcal{V}}(\hat{\boldsymbol{\beta}}) = \{\nabla \overline{H}(\hat{\boldsymbol{\beta}})\}^{-1} \times \hat{\mathcal{V}}\{\overline{H}(\boldsymbol{\beta}_0)\} \times \{\nabla \overline{H}(\hat{\boldsymbol{\beta}})\}^{-1}. \quad (B3)
$$

The derivation of $\hat{\mathcal{V}}\{\overline{H}(\boldsymbol{\beta}_0)\}$ is explained in the next sections.
## B.2 Accounting for the estimation of $\alpha$

Since we have

$$ \frac{1}{\hat{\alpha}} = \frac{1}{\alpha} \times \frac{1}{1 + \frac{\hat{\alpha} - \alpha}{\alpha}} = \frac{1}{\alpha} \left[ 1 - \frac{\hat{\alpha} - \alpha}{\alpha} + o_p(n_V^{-0.5}) \right] \\ = \frac{1}{\alpha} - \frac{\hat{\alpha} - \alpha}{\alpha^2} + o_p(n_V^{-0.5}), $$

we may rewrite the quantities in (B2) as

$$ \mathbf{X}_i^*(\hat{\alpha}) = \underbrace{\frac{1}{\alpha}(\mathbf{Z}_i - \bar{\mathbf{X}}_B) + \bar{\mathbf{X}}_B}_{\mathbf{X}_i^*(\alpha)} - \frac{\hat{\alpha} - \alpha}{\alpha^2}(\mathbf{Z}_i - \bar{\mathbf{X}}_B) + o_p(n_V^{-0.5}), \\ g_j^*(\hat{\alpha}, \beta_0) = \underbrace{\frac{1}{\alpha}\{g(\beta_0, \mathbf{Z}_j) - \bar{g}_B(\beta_0)\} + \bar{g}_B(\beta_0)}_{g_j^*(\alpha, \beta_0)} - \frac{\hat{\alpha} - \alpha}{\alpha^2}\{g(\beta_0, \mathbf{Z}_j) - \bar{g}_B(\beta_0)\} + o_p(n_V^{-0.5}), \quad (B4) $$

$$ h_j^*(\hat{\alpha}, \beta_0) = \underbrace{\frac{1}{\alpha}\{h(\beta_0, \mathbf{Z}_j) - \bar{h}_B(\beta_0)\} + \bar{h}_B(\beta_0)}_{h_j^*(\alpha, \beta_0)} - \frac{\hat{\alpha} - \alpha}{\alpha^2}\{h(\beta_0, \mathbf{Z}_j) - \bar{h}_B(\beta_0)\} + o_p(n_V^{-0.5}). \quad (B5) $$

Let us denote $\epsilon = \frac{\hat{\alpha} - \alpha}{\alpha^2}$. By plugging (B4) and (B5) into Equation (B1), we have

$$ \frac{\sum_{j=1}^{n_A} Y_j(T_i) h_j^*(\hat{\alpha}, \beta)}{\sum_{j=1}^{n_A} Y_j(T_i) g_j^*(\hat{\alpha}, \beta)} = \frac{\sum_{j=1}^{n_A} Y_j(T_i) h_j^*(\alpha, \beta) - \epsilon \sum_{j=1}^{n_A} Y_j(T_i) [h(\beta_0, \mathbf{Z}_j) - \bar{h}_B(\beta_0)]}{\sum_{j=1}^{n_A} Y_j(T_i) g_j^*(\alpha, \beta) - \epsilon \sum_{j=1}^{n_A} Y_j(T_i) [g(\beta_0, \mathbf{Z}_j) - \bar{g}_B(\beta_0)]} + o_p(n_V^{-0.5}). $$

After some algebra, this leads to:

$$ \bar{H}(\beta_0) = \bar{H}_1(\beta_0) - \left(\frac{\hat{\alpha} - \alpha}{\alpha^2}\right) \bar{H}_2(\alpha, \beta_0) + o_p(n_V^{-0.5}), \quad (B6) $$

where

$$ \bar{H}_1(\beta_0) = \frac{1}{n_A} \sum_{i=1}^{n_A} \underbrace{\delta_i \{ \mathbf{X}_i^*(\alpha) - R_i^*(\alpha, \beta_0) \}}_{H_{1i}(\beta_0)} \quad (B7) $$

with $R_i^*(\alpha, \beta_0) = \frac{\sum_{j=1}^{n_A} Y_j(T_i) h_j^*(\alpha, \beta)}{\sum_{j=1}^{n_A} Y_j(T_i) g_j^*(\alpha, \beta)}$, and with

$$ \bar{H}_2(\alpha, \beta_0) = \frac{1}{n_A} \sum_{i=1}^{n_A} \delta_i \left[ (\mathbf{Z}_i - \bar{\mathbf{X}}_B) - \frac{\sum_{j=1}^{n_A} Y_j(T_i) \{[h(\beta_0, \mathbf{Z}_j) - \bar{h}_B(\beta_0)] - R_i^*(\alpha, \beta_0) [g(\beta_0, \mathbf{Z}_j) - \bar{g}_B(\beta_0)]\}}{\sum_{j=1}^{n_A} Y_j(T_i) g_j^*(\alpha, \beta_0)} \right]. \quad (B8) $$
By neglecting the terms which are $o_p(n_V^{-0.5})$, we obtain from (B6) that

$$ 
\begin{aligned}
\mathbb{V}[\bar{H}(\boldsymbol{\beta}_0)] &= \mathbb{V}[\mathbb{E}\{\bar{H}(\boldsymbol{\beta}_0)|\mathbf{X}_B, T_A, \mathbf{Z}_A\}] + \mathbb{E}[\mathbb{V}\{\bar{H}(\boldsymbol{\beta}_0)|\mathbf{X}_B, T_A, \mathbf{Z}_A\}] \\
&\approx \mathbb{V}[\bar{H}_1(\boldsymbol{\beta}_0)] + \mathbb{E}[\bar{H}_2(\boldsymbol{\beta}_0)\mathbb{V}\{\frac{\hat{\alpha} - \alpha}{\alpha^2}|\mathbf{X}_B, T_A, \mathbf{Z}_A\}\{\bar{H}_2(\boldsymbol{\beta}_0)\}^\top].
\end{aligned} 
\quad (B9) $$

Under the assumption that the validation sample $S_V$ is selected in A by simple random sampling without replacement, we have

$$ \hat{\alpha} = \frac{1}{n_V} \sum_{i \in S_V} \mu_i \quad \text{where} \quad \mu_i = \begin{cases} 1 & \text{if linkage is correct,} \\ 0 & \text{otherwise.} \end{cases} $$

Since $\mu_i$ is a binary variable, it follows from standard results in survey sampling theory that an unbiased estimator for $\mathbb{V}\{\hat{\alpha}|\mathbf{X}_B, T_A, \mathbf{Z}_A\}$ is

$$ \hat{V}(\hat{\alpha}) = \left(\frac{1}{n_V} - \frac{1}{n_A}\right) \frac{n_V}{n_V - 1} \hat{\alpha}(1 - \hat{\alpha}). $$

Hence the second term in the right-hand side of (B9) may be estimated by

$$ \hat{V}_1[\bar{H}(\boldsymbol{\beta}_0)] = \bar{H}_2(\hat{\alpha}, \hat{\beta})\{\bar{H}_2(\hat{\alpha}, \hat{\beta})\}^\top \times \left(\frac{1}{n_V} - \frac{1}{n_A}\right) \frac{n_V}{n_V - 1} \frac{1 - \hat{\alpha}}{\hat{\alpha}^3}, \quad (B10) $$

where $\bar{H}_2(\hat{\alpha}, \hat{\beta})$ is obtained from (B8) by replacing $\boldsymbol{\beta}_0$ with $\hat{\beta}$ and $\alpha$ with $\hat{\alpha}$. This is the component of the variance estimator which accounts for the estimation of $\alpha$.

## B.3 Accounting for the linkage and estimation error

In this section, we focus on the first term in the right-hand side of (B9). We have

$$ \mathbb{V}[\bar{H}_1(\boldsymbol{\beta}_0)] = \mathbb{V}[\mathbb{E}\{\bar{H}_1(\boldsymbol{\beta}_0)|\mathbf{X}_B, T_A\}] + \mathbb{E}[\mathbb{V}\{\bar{H}_1(\boldsymbol{\beta}_0)|\mathbf{X}_B, T_A\}]. \quad (B11) $$

It follows from Equation (A6) in Appendix A that

$$ \mathbb{E}\{\bar{H}_1(\boldsymbol{\beta}_0)|\mathbf{X}_B, T_A\} \approx \frac{1}{n_A} \sum_{i=1}^{n_A} \underbrace{\delta_i\left\{\mathbf{X}_i - \frac{\sum_{j=1}^{n_A} Y_j(T_i)h(\boldsymbol{\beta}_0, \mathbf{X}_j)}{\sum_{j=1}^{n_A} Y_j(T_i)g(\boldsymbol{\beta}_0, \mathbf{X}_j)}\right\}}_{H_i(\boldsymbol{\beta}_0)}, \quad (B12) $$

which is the function associated to the theoretical estimating equation that we would solve if the covariates $\mathbf{X}_i$ were known without linkage error for the units $i \in A$. Secondly, note that conditionally on $\mathbf{X}_B$ and $T_A$, the terms $H_{1i}(\boldsymbol{\beta}_0)$ are approximately uncorrelated for $i = 1, \dots, n_A$. More precisely, it can be proved after some algebra that for any $i \neq j = 1, \dots, n_A$, we have

$$ \mathrm{Cov}\left(\delta_i\{\mathbf{X}_i^*(\alpha) - R_i^*(\alpha, \boldsymbol{\beta}_0)\}, \delta_j\{\mathbf{X}_j^*(\alpha) - R_j^*(\alpha, \boldsymbol{\beta}_0)\} \middle| \mathbf{X}_B, T_A\right) = O_p(n_A^{-1}). $$

Therefore, we obtain that

$$ \mathbb{V}\{\bar{H}_1(\boldsymbol{\beta}_0)|\mathbf{X}_B, T_A\} \approx \frac{1}{(n_A)^2} \sum_{i=1}^{n_A} \mathbb{V}\{H_{1i}(\boldsymbol{\beta}_0)|\mathbf{X}_B, T_A\}. \quad (B13) $$
where $H_{1i}(\cdot)$ is defined in (B7). From (B11), (B12) and (B13), we obtain that

$$ \mathbb{V}[\bar{H}_1(\boldsymbol{\beta}_0)] \approx \mathbb{V}\left(\frac{1}{n_A} \sum_{i=1}^{n_A} H_{ti}(\boldsymbol{\beta}_0)\right) + \mathbb{E}\left[\frac{1}{(n_A)^2} \sum_{i=1}^{n_A} \mathbb{V}\{H_{1i}(\boldsymbol{\beta}_0) | \mathbf{X}_B, T_A\}\right]. \quad (B14) $$

Now, we consider the sample dispersion term given by

$$ \begin{aligned} s_H^2(\boldsymbol{\beta}_0) &= \frac{1}{n_A - 1} \sum_{i=1}^{n_A} \left\{ H_i(\boldsymbol{\beta}_0) - \frac{1}{n_A} \sum_{j=1}^{n_A} H_j(\boldsymbol{\beta}_0) \right\}^2 \\ &= \frac{1}{2n_A(n_A - 1)} \sum_{i \neq j=1}^{n_A} \{H_i(\boldsymbol{\beta}_0) - H_j(\boldsymbol{\beta}_0)\}^2. \end{aligned} \quad (B15) $$

where $H_i(\cdot)$ is defined in (B1).

We have

$$ \begin{aligned} \mathbb{E}\left\{\frac{s_H^2(\boldsymbol{\beta}_0)}{n_A}\right\} &= \mathbb{E}\mathbb{E}\left\{\frac{s_H^2(\boldsymbol{\beta}_0)}{n_A} \middle| \mathbf{X}_B, T_A\right\} \\ &= \mathbb{E}\left[\frac{1}{2n_A^2(n_A - 1)} \sum_{i \neq j=1}^{n_A} \mathbb{E}\{H_i(\boldsymbol{\beta}_0) - H_j(\boldsymbol{\beta}_0) | \mathbf{X}_B, T_A\}^2\right] \\ &\quad + \mathbb{E}\left[\frac{1}{2n_A^2(n_A - 1)} \sum_{i \neq j=1}^{n_A} \mathbb{V}\{H_i(\boldsymbol{\beta}_0) - H_j(\boldsymbol{\beta}_0) | \mathbf{X}_B, T_A\}\right] \\ &\approx \mathbb{E}\left[\frac{1}{2n_A^2(n_A - 1)} \sum_{i \neq j=1}^{n_A} \{H_{ti}(\boldsymbol{\beta}_0) - H_{tj}(\boldsymbol{\beta}_0)\}^2\right] \\ &\quad \text{(where } H_{ti}(\cdot) \text{ is defined in (28)}) \\ &\quad + \mathbb{E}\left[\frac{1}{2n_A^2(n_A - 1)} \sum_{i \neq j=1}^{n_A} \mathbb{V}\{H_i(\boldsymbol{\beta}_0) | \mathbf{X}_B, T_A\} + \mathbb{V}\{H_j(\boldsymbol{\beta}_0) | \mathbf{X}_B, T_A\}\right] \\ &= \mathbb{E}\left[\frac{1}{n_A(n_A - 1)} \sum_{i=1}^{n_A} \left\{H_{ti}(\boldsymbol{\beta}_0) - \frac{1}{n_A} \sum_{j=1}^{n_A} H_{tj}(\boldsymbol{\beta}_0)\right\}^2\right] \\ &\quad + \frac{1}{n_A^2} \sum_{i=1}^{n_A} \mathbb{V}\{H_i(\boldsymbol{\beta}_0) | \mathbf{X}_B, T_A\} \\ &\approx \mathbb{V}[\bar{H}_1(\boldsymbol{\beta}_0)], \end{aligned} \quad (B16) $$

where the last line in (B16) follows from a comparison with Equation (B14). Therefore, $\mathbb{V}[\bar{H}_1(\boldsymbol{\beta}_0)]$ may be approximately unbiasedly estimated by replacing in (B15) the unknown parameter $\boldsymbol{\beta}_0$ with $\hat{\boldsymbol{\beta}}$, which leads to

$$ \hat{\mathbb{V}}_2[\bar{H}(\boldsymbol{\beta}_0)] = \frac{s_H^2(\hat{\boldsymbol{\beta}})}{n_A}. \quad (B17) $$

This is the component of the variance estimator, which accounts for both the linkage and estimation errors.

## B.4 Global variance estimator

By plugging (B10) and (B17) into (B9), we obtain:

$$ \hat{\mathbb{V}}\{\bar{H}(\boldsymbol{\beta}_0)\} = \hat{\mathbb{V}}_1\{\bar{H}(\boldsymbol{\beta}_0)\} + \hat{\mathbb{V}}_2\{\bar{H}(\boldsymbol{\beta}_0)\}. $$
The global variance estimator is therefore obtained from (B3) as:

$$ \hat{\mathcal{V}}(\hat{\beta}) = \{\nabla \overline{H}(\hat{\beta})\}^{-1} \times \{\hat{\mathcal{V}}_1\{\overline{H}(\beta_0)\} + \hat{\mathcal{V}}_2\{\overline{H}(\beta_0)\}\} \times \{\nabla \overline{H}(\hat{\beta})\}^{-1} \quad (B18) $$

## APPENDIX C. SAMPLING AFFECTATIONS

To compute $\overline{\mathbf{X}}_{B_v}$, $\overline{\mathbf{g}}_{B_v}(\boldsymbol{\beta})$ and $\overline{\mathbf{h}}_{B_v}(\boldsymbol{\beta})$ in (7), the AEE requires access to all the $\mathbf{X}$-vectors in $B$. In some cases, this may not be possible due to confidentiality reasons. In that case, we have access to only the linked dataset $A$. In this situation, we propose to approximate

$$ \overline{\mathbf{X}}_{B_v} \quad \text{with} \quad \overline{\mathbf{Z}}_{A_v} = \frac{1}{n_{A_v}} \sum_{i \in A_v} \mathbf{Z}_i, $$

$$ \overline{\mathbf{g}}_{B_v}(\boldsymbol{\beta}) \quad \text{with} \quad \overline{\mathbf{g}}_{A_v}(\boldsymbol{\beta}) = \frac{1}{n_{A_v}} \sum_{i \in A_v} \exp(\mathbf{Z}_i^\top \boldsymbol{\beta}), $$

$$ \overline{\mathbf{h}}_{B_v}(\boldsymbol{\beta}) \quad \text{with} \quad \overline{\mathbf{h}}_{A_v}(\boldsymbol{\beta}) = \frac{1}{n_{A_v}} \sum_{i \in A_v} \exp(\mathbf{Z}_i^\top \boldsymbol{\beta}) \mathbf{Z}_i. \quad (C1) $$

If $A_v$ is a random sample of $B_v$, (C1) can be a good approximation. A simulation study is presented in Table S4 of the Supplementary material.