This study utilised data from the Canadian Primary Care Sentinel Surveillance Network (CPCSSN), a pan-Canadian collection of electronic medical records (EMR) encompassing over 2 million primary care patients spread across 14 regional networks [21]. CPCSSN has pan-Canadian coverage of primary care practices, with participating practices situated in the majority of provinces across Canada. Since its establishment in 2008, CPCSSN continuously receives EMR data extracted by contributing networks from participating primary care practitioners. CPCSSN contains information on consenting patients’ diagnoses, procedures, laboratory tests, medication prescriptions, and referrals. As CPCSSN data is pooled at the national level, it is one of the few data sources that allows for national level data analysis while also including meaningful differences across jurisdictional boundaries rather than creating artificial boundaries as has been done by others [22]. This allows for the assessment of this novel approach by providing a basis for comparison against the actual national values.

Adult patients (18 years of age and older) with diabetes (n = 42,341) were identified using validated case definitions [23, 24]. A unique baseline visit was defined for each patient as their first visit with a participating primary care practitioner between January 1, 2014, and December 31, 2014. In keeping with using CKD and hypertension as a case study, our analyses were restricted to patients without an existing diagnosis of CKD at any time prior to baseline and those who were not marked deceased at their baseline visit.

Diagnoses of CKD were identified using eGFR. Specifically, two or more laboratory values separated by at least 90 days but not more than 1 year apart reporting an eGFR less than 60 mL/min/1.73 m², regardless of age as well as two or more laboratory values separated by at least 90 days but not more than 1 year apart reporting an eGFR less than 75, 60, and 45 mL/min/1.73 m² for ages younger than 40, 40 to 64, and 65 years or older, respectively [25].

The variables examined in this study include hypertension, sex, and age. Hypertension is a dichotomous variable based on any diagnosis prior to baseline [23]. Sex is a dichotomous variable based on what was recorded in the EMR. Age is treated as a continuous variable ranging from 18 to 104, without any transformations. Age was calculated as baseline visit year minus birthyear.

All analyses were conducted using R version 4.4.1.

Five logistic regression models were fit separately for British Columbia (n = 1,700), Alberta (n = 8,002), Manitoba (n = 6,089), Ontario (n = 23,946), and the Maritimes (n = 2,604) which includes patients from Newfoundland and Labrador, New Brunswick, Nova Scotia, and Prince Edward Island. These regions were aggregated due to sample size. Additionally, a global logistic regression model was fit which encompasses all individual patient data across the provinces (n = 42,341). In each model, hypertension serves as the primary exposure variable, while sex and age are included as covariates to control for potential confounding effects. The province of Quebec has been excluded from the analyses as there were no recorded patients with the outcome of interest. The province of Saskatchewan and the Canadian territories are not included as CPCSSN does not currently collect data in those regions.

Common effect meta-analysis was employed to pool the coefficients from the five regional models using the R package ‘meta’. A separate meta-analysis was conducted for each model parameter to synthesise the provincial estimates of CKD across Canada. The inverse variance method was used for pooling [26]. The heterogeneity across regions was assessed using Cochran’s Q, which is a chi-square statistic comparing the effect size from each regional model to the overall pooled effect. Regional heterogeneity was measured using I², which is the proportion of total variation in effects due to true heterogeneity rather than chance.

The federated likelihood method is built upon a core set of assumptions that support both the validity and interpretability of inferences drawn from distributed data sources. These assumptions are critical for ensuring the reliability and applicability of the method, and are outlined as follows:

1. Common Effect Assumption:

A central assumption is that of a common effect. The fundamental objective of this method is to estimate a global parameter that is common across all sites. If true heterogeneity is present, more flexible models such as mixed-effects approaches may be warranted which is beyond the scope of this study.

2. Scientific Pooling Assumption (Justifiability of Pooled Analysis):

The primary assumption is that data from all participating jurisdictions or institutions can be justifiably pooled for unified analysis. This requires that the research question, study populations, and measured variables are sufficiently comparable to allow for meaningful inference. This assumption should be critically assessed a priori: if conceptual alignment across jurisdictions is lacking, then pooling – no matter how technically feasible – is inappropriate.

3. Consistency of Model Structure Across Sites:

To construct a valid joint likelihood, all participating sites must use an identical model structure. This includes consistent specification of the outcome variable, covariates, link function, and distributional assumptions. Structural inconsistencies, such as differences in variable inclusion, outcome definitions, missingness, or data transformations, would undermine compatibility and compromise the validity of joint inference.

4. Independence of Observations Across Jurisdictions (IID Assumption):

Finally, the method assumes that observations are independent and identically distributed (IID) across jurisdictions, with no overlapping or duplicated records. This condition is crucial to avoid double-counting, which can bias the likelihood.

These assumptions form the theoretical backbone of the federated likelihood framework. In practice, each must be carefully assessed within the specific context of the research question and available data. When any assumption is violated, alternative methods, such as mixed-effects modelling or meta-analytic approaches, may offer more suitable solutions.

As it is often impractical to share a closed-form solution to the likelihood function since such solutions typically don’t exist for complex models and must be approximated numerically, we investigate a brute force approach in which regional log-likelihood matrices are constructed and aggregated, eliminating the need to reconstruct full likelihood functions. This brute force approach utilises a grid search strategy to calculate the log-likelihood at every combination of plausible parameter values thereby constructing a matrix of log-likelihood values.

A predefined range of possible parameter values is created for each parameter in the model. These parameter values are drawn from the 95% confidence intervals associated with each of the estimates from the regional multivariable logistic regression models. Using confidence intervals to define this range is beneficial because they reflect the variability in our estimates, allowing us to incorporate regional variation into the grid search for selecting parameter values in our fully combined model. To construct the range of values for a given parameter, the smallest (most conservative) lower bound and the largest (most conservative) upper bound are identified across each of the regional confidence intervals for that parameter. These two values – representing the lowest and highest plausible values of the parameter across all models – form the endpoints of the parameter space for that vector. Once this inclusive range is established, it is segmented by generating a sequence of values at small, uniform increments. This study will explore the performance across two vector increments: 0.01 and 0.005. This results in a finely resolved vector that effectively captures the full extent of the parameter uncertainty across models. This approach ensures that no plausible values are excluded from the grid search and allows for the construction of a combined model that fully reflects the range of estimates derived from the ensemble of individual models.

Within each region, the log-likelihood was then calculated for each unique combination of values from the parameter vectors. The result of this process is a matrix of log-likelihood values, with each element of the matrix representing the likelihood of the data given a specific combination of parameter values. Technically, this will be a tensor for more than two parameters, but for simplicity we will refer to this as a matrix. This matrix effectively maps out the likelihood landscape for the regional model, allowing for a detailed examination of how the likelihood changes across different parameter combinations. By calculating the log-likelihood across all potential sets of parameter values, each matrix includes the "true" maximum log-likelihood, corresponding to the set of parameter values that best fit the regional data.

This process was then repeated independently for each regional model, generating a matrix of log-likelihood values for each region. The regional log-likelihood matrices were then summed. Summing the log-likelihoods at each intersection of parameter values leverages the additive property of log-likelihoods in independent models, allowing for the combination of information from the regions without directly accessing their raw data. The combined log-likelihood matrix represents the overall likelihood for the global model, effectively synthesising the regional information into a single estimate. An illustration of this process using simulated data is presented in Figure 1.

In addition, the 14.7% likelihood region was identified, representing the set of parameter values where the likelihood is at least 14.7% of the maximum likelihood. From this region, a 14.7% likelihood interval for each parameter can be derived, defining the range of parameter values around the maximum where the likelihood of the data is highest. The 14.7% cutoff is chosen because it approximately corresponds to the likelihood threshold associated with a 95% confidence interval under likelihood theory [27]. We report both likelihood and confidence intervals.

The parameter values that maximised the combined likelihood were identified and compared to the global model estimates to evaluate the performance of the likelihood-based approach. In addition, the estimates obtained through the likelihood method were directly compared to those generated by meta-analysis. The meta-analysis, which aggregates odds ratios across regions, and the likelihood-based method, which combines regional log-likelihoods, represent two distinct strategies for synthesising regional data into combined estimates. To assess the accuracy of each approach, the common effect meta-analysis estimates, along with the federated likelihood estimates, were compared to the global model by calculating the relative percent absolute bias (RPAB). This was defined as the absolute difference between each pooled estimate and the global estimate, divided by the global estimate. Comparing the RPABs across methods provides a measure of how closely each approach approximates the true global model and highlights potential limitations or biases inherent in either the meta-analysis or federated likelihood method.

A common effect meta-analysis was conducted to pool the coefficients from the five regional models. The regional coefficients for hypertension, and their standard errors (shown in Table 2) were used for the meta-analysis. The results are summarised in the forest plot presented in Figure 2.

The meta-analysis for hypertension produced pooled estimates that were compared to the known global model estimate of 0.25 (95% CI: 0.20 to 0.31). The common effect meta-analysis yielded an estimated log-odds ratio of 0.26 (95% CI: 0.20 to 0.32), which closely approximated the global estimate and included it within its confidence interval. However, substantial heterogeneity was detected across regions (I² = 80%, p = 0.0005), suggesting meaningful variation in the association between hypertension and CKD across regions. This level of heterogeneity suggests that the assumption of a common effect may not hold. In contexts where effect heterogeneity is substantial, estimates may be biased or oversimplified because they have been pooled across jurisdictions without adequately accounting for true between-site variation, therefore alternative approaches such as mixed-effects or random-effects models could be more appropriate. Nonetheless, the primary aim of this study was to evaluate the performance of the federated likelihood method, rather than to identify the most suitable model for inference. As such, the use of a common effect model remains justifiable within the scope of our objectives.

A wide parameter space was generated based on the confidence intervals from the regional multivariable logistic regression models estimating the association between hypertension and CKD. The development of the parameter vectors is illustrated in Figure 3.

Within each region, the log-likelihood was calculated for every unique combination of values from the parameter vectors, resulting in 4,438,496 evaluations per region given the threshold of 0.01, and 61,151,895 evaluations per region given the threshold of 0.005. This process produced a matrix of log-likelihood values for each regional model, effectively mapping the likelihood landscape and enabling a detailed examination of how the likelihood varies across different parameter combinations. Each matrix includes the "true" maximum log-likelihood, corresponding to the set of parameter values that best fit the regional data. The regional log-likelihood matrices were then summed to create a combined log-likelihood matrix, which was subsequently maximised to identify the parameter values that best fit the data across regions. The results of this process are presented in Table 3.

At the regional level, the estimates obtained using the federated likelihood approach closely approximated the estimates from the regional logistic regression models. Across all regions, the differences in maximum log-likelihood values between the regional and federated likelihood models (both the 0.01 and 0.005 vector increments) were negligible – typically less than one unit – indicating strong agreement in model fit. Similarly, coefficient estimates from the federated likelihood models remained consistent with those from the regional models, suggesting that the federated estimation procedure preserved site-specific inference.

For the primary exposure of interest (hypertension), the combined federated likelihood model using the 0.01 vector increment produced an estimated log-odds ratio of 0.2476, closely approximating the global estimate of 0.2542. Using a smaller increment of 0.005 further improved alignment with the global model: the log-odds ratio increased to 0.2576, and the maximum log-likelihood (–15,564.56) approached that of the global model (–15,563.98), outperforming the 0.01 model (–15,567.06). These results suggest that tighter increments enhance estimation accuracy.

Additionally, adjustment for potential confounding was maintained across all modelling strategies. By ensuring that the same covariates (sex and age) were included across jurisdictional models, the federated likelihood method was able to control for key sources of confounding.

Overall, the likelihood method effectively estimated the true global model parameters without requiring individual-level data transfer. Although slight deviations were present, the magnitude of the discrepancies in effect size was small and did not impact the overall interpretation of the findings. These results highlight the utility of the federated likelihood approach for interjurisdictional health data analysis. Table 4 summarises all of the results previously discussed.

The common effect meta-analysis estimates, along with the federated likelihood estimates, were compared to the global model to evaluate their accuracy. Relative percent absolute bias (RPAB), defined as the absolute difference between each pooled estimate and the global estimate divided by the global estimate, was calculated to quantify the bias in each method. Figure 4 summarises the RPAB for the coefficients for hypertension across the two approaches compared to the global model. Among all methods, the federated likelihood model with the 0.005 increment achieved the lowest RPAB (1.34%), indicating the closest approximation to the global model. The model with the 0.01 increment also performed well but showed slightly higher bias (2.60%). Both federated likelihood approaches outperformed the common effects meta-analysis, which produced moderately accurate estimates (4.09% RPAB).

Overall, the federated likelihood method produced estimates that were highly consistent with those obtained from the common effects meta-analysis approach. Notably, the federated likelihood method provided the most accurate estimate for the primary effect of interest, demonstrating a smaller RPAB compared to the common effect meta-analysis. The federated likelihood method therefore offers a robust alternative that maintains both accuracy and interpretability in federated data analysis.