Distributed Statistical Analyses: A Scoping Review and Examples of Operational Frameworks Adapted to Health Analytics

Introduction

Health Analytics at Scale

Learning health systems (LHSs) are coming of age and are being deployed to address important health challenges at different scales. The framework starts by leveraging health data created across various activities. It obviously includes data points from clinics and hospitals, but the perimeter of data required to meaningfully and optimally address important problems is much wider and includes research cohorts, biobanks, quantified self-data, environmental exposures, and social service delivery.

While some questions might be addressed at the scale of an individual organization, LHSs focus on system interactions and often require the analysis of processes and outcomes from various organizations. For example, to fully understand a cancer care trajectory, multiple data sources from multiple organizations will need to be examined to cover all relevant aspects (both within the traditional health system and in the community). This often implies organizations that are at least regional or of a wider scope (provinces, states, and countries), such as in the context of the Health Data Research Network Canada or Health Data Research United Kingdom. Similarly, comparing various approaches is often a fruitful way to identify the best approaches and understand what works, why, and how to scale the promising projects. It can also be a way to amass a critical number of observations in the context of rarer diseases. Nevertheless, working with data from multiple sources, from multiple organizations, and located in multiple jurisdictions poses significant challenges.

Traditionally, the analytical methods used by researchers in health-related and other domains have relied on data pooling (sometimes referred to as data centralization)—all required data are physically copied to a single location where analysis can take place. However, when working with data from multiple jurisdictions (even when part of the same country, such as the Canadian provinces and territories), data pooling is often very difficult, if not impossible, for ethical, legal, and social acceptability reasons.

Therefore, there is a pressing need to offer analytical methods allowing for the analysis of such data without the need to physically copy the data in a central location.

The primary intent of this study was to lay the foundations for future practical assessments of the feasibility of conducting distributed statistical analyses in health-related contexts. We achieved this by reviewing the literature on existing methods, evaluating which ones could be applied to a class of widely used regression models, and precisely identifying their operational and information sharing requirements in a notational framework that allows for straight comparisons. This unified framework enables the description of methods’ operational workflows, quantities exchanged, and algorithmic implementation and operates under assumptions commonly satisfied in health analytics.

This paper is structured as follows. We begin with a formal description of the distributed analytical framework considered in this study, followed by a discussion on the challenges associated with its implementation in health analytics. Next, we state our specific objectives and outline the methodology used to achieve them, including a scoping review and our approach to establishing the common notational framework for method description and comparison. After presenting the results pertaining to each objective, we discuss our findings and remaining challenges regarding distributed statistical analyses for health analytics.

Distributed Analysis

Overview

This study was concerned with frameworks in which the data needed for a statistical analysis consisted of the data about n individuals (referred to as the analytical data set), which are not all stored in a single source but are partitioned among K locations that will be called nodes hereafter. Therefore, the mereological sum of all the data held at each node forms the analytical data set. Data can be partitioned horizontally or vertically (or in a mixed way).

A horizontal partition implies that all data pertaining to a given individual can be found in a single node. If we assume that patients receive care only in 1 province, Canadian provincial health administrative data sets hosted by organizations such as Population Data BC, the Institute for Clinical Evaluative Sciences in Ontario, or the Manitoba Centre for Health Policy in Manitoba can be part of a horizontal partition. A clinical trial in which each recruiting site captures all data for a given participant is another example.

A vertical partition occurs when all data of a certain type are available in a single node for a group of individuals. A classic example is a hospital with its various information systems. All pathology results can be found in the pathology system, all billing information can be extracted from the finance system, and all x-rays are accessible in the picture archiving and communication system. However, to obtain the full picture of the care received by a patient, multiple systems need to be interrogated. Similarly, in the research setting, health administrative data may be in a provincial data center, and genomics data could be held in a research institute.

A mixed partition occurs when both principles partly apply—some individuals may have their data spread out across nodes, and different individuals may be present in different nodes.

Assumptions

The difficulties in conducting analyses on a large scale mentioned previously are often associated with horizontally partitioned data, and this work focused on this type of partition. Therefore, the methods presented in this paper might not be directly applicable to vertically partitioned data.

One group of approaches often labeled as distributed analysis involves calculations at each participating node and exchanges of the resulting aggregated statistics with a coordinating center (CC), which can itself also perform additional calculations based on the received aggregated statistics. The CC can be an organization not responsible for a data node or a data node taking on the additional role of CC for a given analysis.

It is important to note that, whether in the more traditional way of data pooling or using distributed approaches (in which the data are not copied centrally), data sources will be different on multiple levels. They will represent information using data models with significant variability in terms of structure and technology but also in terms of semantics. This situation also leads to heterogeneous data in which the presence of predictors and outcomes is likely to be different in different nodes. Different approaches (eg, data mediation or extract, transform, and load) have been developed to address these issues, and this work assumed that one of them was applied so that the data nodes mentioned hereafter are assumed to share the same structure, the same technological syntax, and the same semantics, as well as no missing data.

Horizontally Partitioned Statistical Analytics

In what follows, the field that pertains to the statistical analysis of horizontally partitioned and semantically homogeneous data that cannot be consolidated into a central location will be called horizontally partitioned statistical analytics (HPSA).

Methodological contributions to this field have arisen from several streams of literature. Meta-analysis and meta-regression methods [1] can be viewed as part of HPSA (eg, by considering that each node-specific data set belongs to a different pseudostudy). However, their scope is narrower compared to that of HPSA because they typically assume that only established study-level estimates are available as data. Conversely, HPSA allows for the sharing of additional summary statistics between the nodes and the CC, such as gradients and Hessians, to ensure the best possible performance at the global level. As meta-analysis does not leverage any supplementary information that could be obtained from studies with access to patient-level data, it can be susceptible to biased estimation, especially in settings with rare outcomes or in the presence of data nodes with limited sample sizes [2]. As meta-analysis and meta-regression methods have been extensively covered in the literature, approaches specifically designed for the analysis of already established study-level estimates will not be discussed hereafter.

An important research community that has generated a significant amount of analytical contributions is concerned with the massive data setting. There, a data set often cannot be processed by a single server and, therefore, is split across multiple machines, which are then considered as nodes able to perform computations and send aggregated results to a CC tasked with fitting a global model from them. The methodological avenues proposed in this setting share similarities with the ones designed for the multi–research facility setting involved in LHSs but also have important differences. For example, in a massive data setting, the experimenter has control over the distribution of individuals across nodes, which is typically not the case in multi–research facility studies. Thus, while these approaches share mechanistic similarities and have been suggested as options to consider in the health domain, some hypotheses may not hold. In regression settings, it is often reasonable to assume that the regression link between the response and covariate predictors is the same across nodes. However, assuming that the sampling distribution of covariates involved is equal across nodes is unrealistic in the health domain, particularly due to the presence of data centers that may systematically involve different types of patients. For example, certain clinics may predominantly serve older individuals. While this may not affect the estimation of parameter values, it can have implications for computing CIs to ensure the validity of inferences.

So far, 2 reviews discussing methods applicable to horizontally partitioned data have been published [3,4]. However, their focus is on the massive data setting, which works almost invariably under the assumption of even sampling distribution of covariates and equal sample sizes across nodes, and statistical inference tasks beyond parameter estimation are barely covered. This makes them less helpful for health analytics purposes as most studies involving data analyses rely on CIs or hypothesis testing in settings in which predictors’ distribution and sample sizes vary across nodes.

Contemporary Challenges in HPSA

Overview

The problem is 3-fold. First, there is a need to raise awareness regarding the existence of HPSA approaches among researchers aiming at undertaking statistical analyses from horizontally partitioned data, especially in health analytics. The reflex is often to request data pooling because it is perceived as the sole option. This has been the tendency of requests made by researchers to the Health Data Research Network Canada. Practitioners are usually concerned with finding the most appropriate statistical model that will take into account as many of the features of their specific context of application as possible. Consequently, a clear and unifying mapping of the state of the HPSA field is needed for them to be informed of the scope of existing methods available for their analyses to see whether alternatives to pooling exist.

Second, as underlined previously, methodological contributions came from research fields whose working assumptions can be fundamentally different from the ones researchers would be willing to make in health analytics. To ensure proper use of statistical inference techniques, it is necessary that the underlying assumptions of existing methods be adequately identified and understood. If necessary, these methods should be adapted to suit the specific requirements of health applications, thereby ensuring accurate and reliable results.

Third, data custodians have to be properly informed on data sharing requirements entailed by the use of a specific HPSA method applicable to a given research setting. While HPSA avoids the complexities of pooling data, there are still flows of information that have to be acceptable to data stewards. However, even in basic statistical scenarios, available methods are often presented in a way that makes them challenging to compare in terms of information sharing requirements and operational complexity. Therefore, there is a need for clearer and more accessible presentations of these methods to facilitate decision-making regarding data sharing and operational implementation.

Although it would be ideal to offer managers a comprehensive operational workflow for each identified method to evaluate the information shared and execution complexity, with their accompanying underlying modeling assumptions, the abundance and diversity of available approaches make it unfeasible to accomplish this in a single paper. In fact, methods often differ in terms of their targeted application beyond their distributed aspect. For example, differences may exist in the studied model (eg, linear, logistic or Cox regression, and additive models), the dimensionality and sparsity of the predictor variable space, the use of regularization or shrinkage, the presence of missingness, confounders, imbalances, and heterogeneity.

Objectives

The objectives of this study were as follows:

1. To identify and map, from the literature, methodological approaches that make it possible to perform CI estimation and hypothesis testing from a horizontally partitioned data set
2. Among the approaches identified, to describe the ones that allow for the conduct of generalized linear model (GLM) analyses and identify their distributional assumptions
3. On the basis of the approaches identified for GLM-based inferences, to present methods adapted to the setting of uneven sampling distributions across nodes and compare them in terms of information sharing requirements and operational complexity

A scoping review methodology was chosen to achieve objective 1 of mapping the state of the field of HPSA that pertains to inference procedures. For our second objective (objective 2), we identified from the articles selected from the literature search the ones that presented a methodological framework for conducting statistical inference procedures from a GLM with horizontally partitioned data. We then used these frameworks to derive and describe GLM estimators that are applicable to horizontally partitioned data sets. For each identified method, we analyzed and reported its communication workflow and the distributional assumptions. For our third objective (objective 3), we first used statistical theory to adapt the identified procedures to the unequal sample size and uneven covariate distribution setting. Algorithms and mathematical expressions for the quantities involved are reported. For conciseness, we present mathematical formulas for estimation procedures of CIs only. Expressions involved for hypothesis testing are similar and can be deduced following the close connection between CIs and hypothesis tests in GLMs (eg, see Agresti [5]).

The mathematical description of the GLM setting considered for this analysis is described in the following section along with the mathematical notations to be used.

Mathematical Foundations and Notation

Notation

In the following, lowercase letters in bold will represent vector-valued quantities, whereas uppercase letters in bold will denote matrices. The jth element of any vector a ∈ ℝ^p will be denoted as [a]_j. Similarly, the entry at position (j, l) of any matrix A∈ ℝ^p×pwill be denoted as [A]_jl. If g is a real-valued and invertible function, we will use g^(–1) to represent its inverse. In addition, if f_θ is a real-valued function that depends on a parameter vector θ and is twice continuously differentiable, ∇_θf_θ and ∇²_θf_θ will, respectively, indicate the gradient and Hessian matrix of f_θ with respect to θ.

Model Mathematical Assumptions

A mathematical depiction of the horizontally partitioned data framework studied in this paper is as follows. There are n individuals horizontally partitioned across K data storage nodes. Each node’s data set is denoted by ${\mathcal{D}}^{\left(k\right)}=\{{\mathbf{z}}_i^{\left(k\right)}={\left(x_{1i}^{\left(k\right)},\dots ,x_{pi}^{\left(k\right)},y_i^{\left(k\right)}\right)}^{\mathrm{\top }}{\}}_{i=1}^{n^{\left(k\right)}}$, where 1≤ k≤K. Here, ${\mathbf{z}}_i^{\left(k\right)}$ represents the measurements on the ith individual at node k, where ${\displaystyle y_i^{\left(k\right)}\in \mathbb{R}}$ denotes their response variable and ${\displaystyle {\left[x_{1i}^{\left(k\right)},\dots ,x_{pi}^{\left(k\right)}\right]}^{\mathrm{\top }}\in {\mathbb{R}}^p}$ denotes their covariate vector. The total sample size at node k is denoted by n^(k). The combined data set ${\displaystyle {\mathfrak{D}}^{(1)},...,{\mathfrak{D}}^{(K)}}$ make up the whole data set without any duplicated individuals, indicating that $\sum _{k=1}^K{n}^{\left(k\right)}=n$.

Throughout the analysis, it is assumed that each ${\mathbf{z}}_i^{\left(k\right)}$ is independent across 1≤ i≤n^(k) and 1≤ k≤K and there are no missing data. In addition, the size of the covariate space (ie, the dimension of ${\left[x_{1i}^{\left(k\right)},\dots ,x_{pi}^{\left(k\right)}\right]}^{\mathrm{\top }}$, which is equal to p representing the number of features to include as predictors in the GLM) is assumed to be low, eliminating the need for regularization or variable selection. Finally, it is assumed that each node possesses a nonnegligible proportion of the whole data set. Specifically, for each k ∈ {1, ... , K}, the quantity n^(k)/n is bounded away from 0 and 1 as the sample size n tends to infinity, denoted as n^(k)/n → p^(k) ∈ (0,1).

Mathematical Description of the GLM Framework

The formulation of the GLM considered in this paper encompasses various commonly used regression models, such as linear regression, logistic regression, Poisson regression, and probit models. It assumes that the density or probability mass function of each response variable (known as the random components) belongs to the exponential family of distributions. Within this formulation, the mean of the response variable is expressed as a function of a linear combination of the corresponding covariate vector. Formally, it assumes that there exist unknown parameters ${\displaystyle {\mathit{\beta }}^{\star }\in {\mathbb{R}}^{p+1}}$ and ${\displaystyle {\varphi }^{\star }>0}$ and known model-specific functions b, c, g, and h such that, with ${\mathbf{x}}_i^{\left(k\right)}={\left[x_{0i}^{\left(k\right)},x_{1i}^{\left(k\right)},\dots ,x_{pi}^{\left(k\right)}\right]}^{\mathrm{\top }}$ and $x_{0i}^{\left(k\right)}=1$, ${\displaystyle y_i^{\left(k\right)}\mid {\mathbf{x}}_i^{\left(k\right)}\sim f\left(\cdot \hspace{0.25em};\hspace{0.25em}{\mathbf{x}}_i^{\left(k\right)},{\mathit{\beta }}^{\star },{\varphi }^{\star }\right)\hspace{0.17em},}$ where, for any ${\displaystyle \mathit{\beta }={\left[{\beta }_0,{\beta }_1,\dots ,{\beta }_p\right]}^{\mathrm{\top }}\in {\mathbb{R}}^{p+1}}$ and $\varphi$,

${\displaystyle f\left(y;{\mathbf{x}}_i^{\left(k\right)},\mathit{\beta },\varphi \right)=\mathrm{e}\mathrm{x}\mathrm{p}\left[\frac{yh\left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{\left(k\right)}\right)-b\left\{h\left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{\left(k\right)}\right)\right\}}{\varphi }+c\left(y,\varphi \right)\right].}$(1)

In formula 1, b is such that ${\displaystyle b\prime \left\{h\left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{\left(k\right)}\right)\right\}=E\left(y_i^{\left(k\right)}\mid {\mathbf{x}}_i^{\left(k\right)}\right)=g^{\left(-1\right)}\left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{(k)}\right)}$, with b’(x)=∂b(x)/∂x. In this framework, g is called link function, the term ${\displaystyle h\left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{(k)}\right)}$ is usually referred to as the natural parameter, and b is referred to as the cumulant function. $\varphi$ is often called the dispersion parameter and is either known (eg, with $\varphi =1$) or unknown. When h(x)=x (ie, h is the identity function), the link g is called canonical.

The logistic regression model is obtained upon taking $\varphi =1$, h(x)=x, b(x)=log(1+ e^x), $c\left(y,\varphi \right)=0$, and g(x)=log{x/(1 – x)}. The linear regression model with homoscedastic residual error variance $\varphi$ is derived upon setting h(x)=x, b(x)=x²/2, $c\left(y,\varphi \right)=-y^2/\left(2\varphi \right)-\mathrm{l}\mathrm{o}\mathrm{g}\left(2\pi \varphi \right)/2$, and g(x)=x. Hence, both the logistic and the linear regression models rely on a canonical link function in the exponential family distribution.

Methods

Methodology Related to Objective 1

Overview

Scoping reviews are well suited to efficiently map key concepts within a research area [6]. They are widely acknowledged for their ability to clarify working definitions and conceptual boundaries in a specific topic or field [7], facilitating a shared understanding among researchers regarding the status of the research area. These considerations make the scoping review methodology well designed to achieve objective 1.

Scoping studies use systematic searches of relevant databases, using specific keywords to define the boundaries of the research field. However, identifying these keywords can be challenging, particularly when relevant papers are scattered across different research streams or in independent clusters that do not reference each other. To address the risk of overlooking significant methodological contributions due to a limited number of keywords, a snowballing literature search was initially conducted to generate a comprehensive list of keywords related to HPSA. The scoping review then proceeded with a systematic literature search using the identified keywords. It is worth noting that, as the planning of the scoping review is independent of the search approach, the guidelines presented in the work by Arksey and O’Malley [6] are still appropriate.

Methodology Pertaining to the Snowballing Keyword Search

Snowballing is generally used as a literature search method aimed at identifying papers belonging to a given field [8]. It typically consists of three steps: (1) initiate searches in prominent journals or conference proceedings to gather an initial set of papers, (2) conduct a backward review by examining the reference lists of the relevant articles discovered in step 1 (continue iterating until no new papers are found), and (3) perform a forward search by identifying articles that cite the papers identified in the previous steps.

To avoid selection bias, the initial set of papers for the snowballing approach in step 1 is sometimes generated through a search in Google Scholar [9]. The latter strategy was used in this study, too.

As mentioned previously, in this review, the snowballing search strategy was used in preparation for the application of the scoping review protocol with the goal of identifying relevant keywords. Specifically, the starting set of papers was assembled by screening titles and abstracts from the first 50 papers generated through a Google Scholar search using the strings distributed inference and federated inference. The main inclusion criterion was “presents, applies or discusses a statistical inference method to analyse horizontally partitioned data.” The backward and forward snowballing step approaches were then applied.

From the set of keywords found in the selected papers, a list of those relevant to HPSA but not directly associated with any specific method was retained for the scoping review step. It is worth noting that, as the objective of this scoping review was to identify statistical inference methods for horizontally partitioned data, keywords linked to method identifiers had to be excluded from the retained list to avoid preselection bias in the scoping review phase of this project.

The selected keywords that were identified from the snowballing literature search were distributed algorithms, distributed estimation, distributed inference, distributed learning, distributed regression, federated inference, federated estimation, federated learning, privacy-protecting algorithm, privacy-preserving algorithm, and aggregated inference.

Methodology Pertaining to the Scoping Review

The scoping review methodological framework by Levac et al [10] (see also the work by Arksey and O’Malley [6]) was followed. The steps are briefly described in this section. A detailed protocol is available in Multimedia Appendix 1 [2,4,6,10-19].

We conducted a comprehensive search across 4 bibliographic databases—MEDLINE, Scopus, MathSciNet, and zbMATH—to encompass the interdisciplinary nature of the topic and identify relevant research articles. Our search strategies were based on 2 key concepts: distributed data and statistical inference. In addition to the keywords obtained from the snowballing step, we incorporated terms such as confidence interval to target articles focusing specifically on statistical inference. To ensure the inclusion of recent advancements, our search was limited to papers published from 2000 onward. This cutoff date was chosen to account for the emergence of distributed data, the prevalence of massive data sets, and advancements in technology. It was set conservatively to capture any early developed methods and ensure comprehensive coverage of the topic.

After completing the primary search, a 2-stage selection process was used. Initially, 2 authors (MPD and FCL) collaborated to screen all articles identified through the search strategy based on their titles and abstracts. Subsequently, the full texts of the selected articles were independently reviewed by both authors to finalize the selection. This rigorous approach ensured a thorough evaluation of each article’s relevance and eligibility for inclusion.

The primary inclusion criterion for the selection process was as follows: presents a solution for conducting inferential statistics on horizontally partitioned data. This criterion was used to ensure that the chosen articles specifically addressed methods associated with performing statistical inference on horizontally partitioned data.

The following exclusion criteria were derived directly from objective 1: (1) does not address inferential statistics, including CIs, hypothesis testing, or asymptotic normality; (2) does not provide a methodological contribution; and (3) presents a solution for encryption or secret sharing.

To ensure the inclusion of validated approaches, the selection process only considered published papers that had full-text availability in English or French. Discussion papers were excluded as they do not present novel methods or approaches.

Exclusion was considered if any of the exclusion criteria were met or if any of the inclusion criteria were not met.

Finally, the references of each included article from the databases were assessed to identify any relevant articles that may not have been captured during the initial screening due to specific keywords. This additional step in the selection process was necessary given the broad range of vocabulary used to describe applicable approaches in our context.

Data extraction for the included articles was conducted by one author (MPD) and followed a collectively developed data-charting form. Model type (parametric regression, semiparametric regression, nonparametric regression, or not specific to regression) and number of communications from the CC to the nodes (0 or ≥1) were among the data extracted. All methods from the included articles were subsequently classified according to their specified characteristics, as outlined in the protocol. In addition, as part of the analysis, we conducted a screening of the general distributed approaches commonly used across all specific methods.

Methodology Related to Objective 2

Overview

To achieve objective 2, a total of 3 steps were taken. First, we identified methodological approaches from articles included in this scoping review that enable parameter and CI estimations from horizontally partitioned data within a standard GLM framework. Methods designed specifically for the particular cases of linear or logistic regression were also reported but were not analyzed in detail. Second, we extracted workflows for each approach to determine the information exchanged between data storage nodes and the CC. Third, we analyzed the mathematical assumptions necessary for parameter estimation and the consistency of CI procedures. We specifically reported the assumptions related to the distribution of node-specific covariates.

Identification of the Approaches

To identify approaches that enabled the fitting of any GLM using horizontally partitioned data, 2 authors (FCL and MPD) independently assessed all articles included in this scoping review. The reviewers specifically looked for articles that discussed approaches applicable to the GLM class described in the Mathematical Foundations and Notation section, including likelihood-based methods, M-estimation, and estimating equations. In addition, we identified and reported articles that specifically focused on regression settings for linear or logistic regression. However, unless the method described was considered easily adaptable to the GLM framework, these articles were not retained for detailed analysis.

A method was selected if it provided an algorithm for fitting GLMs using horizontally partitioned data, aligning with the characteristics outlined in the Mathematical Foundations and Notation section. In cases in which an article presented asymptotic normality results for the estimators but did not provide an estimator for the asymptotic variance-covariance matrix, the article was still retained, and an estimator for the asymptotic variance was derived using the available calculated quantities.

As our GLM framework assumes no missing values, low dimensionality, and a small number of nodes relative to the total sample size, any terms related to these specific conditions mentioned in an article’s methodology were disregarded. Consequently, the calculations for CIs were adjusted accordingly. If an article solely focused on one of these aspects without contributing to the overall methodology, it was not included in the final selection.

Methodological components regarding parameter estimation and CI procedures were extracted from the screened articles. Specifically, the focus was on understanding how parameters should be estimated within a horizontally partitioned framework and how CIs should be computed for these parameters. For each article, the formulas related to quantities shared among the nodes and quantities calculated by the CC were derived and analyzed. These formulas were examined within a workflow that indicated the necessary circulation of information for the procedure to be executed. The derived workflows constituted the first part of our defined unified framework for HPSA approach comparisons in GLM settings.

For the reported results, the rationale behind each method that was deemed suitable for fitting GLMs was documented, along with the corresponding reference to the paper included in this scoping review in which the method was introduced or discussed.

Articles that discussed approaches specifically applicable to the cases of linear or logistic regression were also mentioned but not elaborated on in detail.

Methodology Related to Objective 3

In most statistical settings with horizontally partitioned data, it is commonly assumed that the sample sizes of the data nodes are equal and that the distribution of covariates is the same across all nodes. However, when the number of nodes is fixed and relatively small compared to the sample sizes, it is possible to adapt a particular approach to handle situations in which the sample sizes and covariate distributions vary across nodes. This can be achieved by combining the theoretical arguments presented in the original article on the method in question with the principles of asymptotic statistics theory concerning maximum likelihood estimation.

To adapt a given approach for situations in which sample sizes and covariate distributions differ across nodes, the following steps were taken:

1. The formulas for the relevant quantities were modified to emphasize the changes caused by this scenario. It was ensured that the adapted quantities were equivalent to their counterparts presented in the original article for an equal sample size setting.
2. Using asymptotic theory, an asymptotic normality result was derived for the estimators of interest considering a set of assumptions that accommodated potential variations in sample sizes and covariate sampling distribution across nodes while still enabling meaningful theoretical arguments.
3. Formulas for the asymptotic variances were derived. Statistical theory on maximum likelihood estimation was used to obtain consistent estimators for asymptotic variances. The latter estimators were derived under the constraint that they had to be calculated without requiring any additional communication round between the CC and the nodes. Thus, throughout the adaptation process, the communication workflow remained unchanged compared to the original method.

These steps ensured the mathematical correctness of adapting the approaches to handle different sample sizes and covariate distributions across nodes. Importantly, these adaptations maintained consistency with the original method’s communication workflows. To adapt and compare the methods based on their information sharing requirements, common assumptions and a unified mathematical and algorithmic notation were necessary. Ultimately, these assumptions and the notation, along with the workflow types derived for objective 2, enabled us to establish a unified notational framework for approach comparisons when performing HPSA based on GLMs in the context of health analytics.

We describe the statistical estimates of interest. A standard GLM typically includes 1 or 2 unknown parametric components. The first are the β parameters, which are commonly assumed to be unknown. The second component is the nuisance parameter $\varphi$, which can be either known (eg, in logistic models) or unknown (eg, in linear models). In practical applications, when $\varphi$ is unknown, its estimated value is often not the main focus, although the latter is necessary to estimate the asymptotic variance of the β parameter estimates.

In the upcoming analysis, we will assume that the parameter $\varphi$ is unknown and estimated using the recommended approach in the selected methods. However, in cases in which $\varphi$ is known, the process becomes simpler. This involves substituting the known value of $\varphi$ and disregarding the estimation step. It is important to highlight that estimating $\varphi$ requires additional information to be shared between the nodes and the CC but it does not necessitate any extra communication round between them.

The estimation processes for both the β and $\varphi$ parameters are discussed. In addition, we explain how to compute an estimator for the asymptotic variance specifically for the estimator of ${\displaystyle {\mathit{\beta }}^{\star }}$. It is important to note that the results presented in the following section can be modified and extended to develop a similar procedure for estimating ${\varphi }^{\star }$.

Using these results, based on an estimator of β (eg, ${\displaystyle \hat{\mathit{\beta }}}$) and a formula for the estimator of the asymptotic variance-covariance matrix involved in its associated asymptotic normality result (eg, ${\displaystyle \hat{\mathbf{\Sigma }}}$), Wald-type (1–α) CIs can be computed for each component of ${\displaystyle {\mathit{\beta }}^{\star }}$ using the following formula:

${\displaystyle {\left[\hat{\mathit{\beta }}\right]}_j\pm z_{1-\alpha /2}\sqrt{{\left[\hat{\mathbf{\Sigma }}\right]}_{jj}/n}\text{for }j\in \left\{1,...,p+1\right\}}$

Regarding the reported results, for each approach considered, we present the formulas necessary to compute the final estimates of the β parameters and their corresponding CIs. The presentation of these formulas was designed to emphasize the communication workflow. Furthermore, a comprehensive algorithm is provided outlining the step-by-step process.

In addition, the asymptotic normality of the β parameter estimators is stated accompanied by the formula for the asymptotic variance and its consistent estimator. Detailed proofs for these results can be found in Multimedia Appendix 2 [20-22].

Results

Results Related to Objective 1

Search Outcomes From the Scoping Review

As presented in Figure 1, a total of 1407 articles were initially identified across all 4 databases after removing duplicates. Subsequently, most of these articles (1274/1407, 90.55%) were excluded based on the evaluation of titles and abstracts, leaving 9.45% (133/1407) of the articles for eligibility assessment through full-text review. Following this assessment, 29 articles were included from the databases. In addition, by reviewing the references of the included articles, 12 more articles were identified and added to the study.

‎

Regarding the additional 12 articles obtained through the assessment of references of the included articles, it was observed that most of them did not mention statistical inference or related terms in their abstracts [2,11,23]. Consequently, these articles were not captured in the initial database search results. Furthermore, some articles directly referred to the specific method used without including any keywords related to horizontally partitioned data in their abstracts or titles [24,25], which greatly reduced the chance of initially identifying them. However, during the process of reviewing the references of the included articles, all the relevant papers that were initially identified through the snowballing strategy were eventually retrieved either through the search strategy or the selection process based on the references of the included articles.

Results of the Scoping Review

Each article included in this scoping review put forth one or multiple methodological approaches pertaining to objective 1. The similarities and differences regarding the communication schemes involved and their background of origin are summarized in this section.

First, all the selected articles discussed one or more statistical procedures that operate on horizontally partitioned data using one of the communication schemes depicted in Figures 2-5.

‎

‎

‎

‎

In workflow I, as shown in Figure 2, each node calculates summary statistics from its own samples, and the results are sent to the CC. The CC combines the information provided by each node to produce the final estimates. This communication approach is commonly referred to as “one-shot” or “noniterative” in the literature, although not always consistently.

In workflow II, as shown in Figure 3, multiple communication rounds are allowed between the CC and the data storage nodes. This allows for iterative interactions between the nodes and the CC to refine the estimates.

Some approaches fundamentally differ from the 2 previous workflows by assigning a different role to one of the nodes (eg, node 1) compared to the others. These approaches operate using workflow III, as illustrated in Figure 4, where node 1 follows a distinct communication pattern compared to the other nodes. In the papers included in this scoping review that discussed these approaches, node 1 was invariably designated as the CC. However, in the context of this paper, their roles were differentiated. The additional step performed by the CC, which involves data aggregation, can be particularly well suited for privacy protection purposes in practice.

The particular setting shown in workflow IV in Figure 5 requires 2 back-and-forth communication exchanges between each node and the CC at each iteration. This communication pattern distinguishes this workflow from the others.

In light of the preceding discussion, from an operational standpoint, 2 categories of workflows emerge. On the one hand, there are workflows that do not necessitate any communication from the CC to the nodes, which are captured in workflow I. On the other hand, there are workflows that involve one or more communication exchanges from the CC to the nodes, which are captured in workflows II, III, and IV.

To emphasize similarities among the methods presented in the articles in this review and facilitate the identification of methods suitable for specific purposes, a systematic classification is presented in Table 1. The articles are categorized based on the type of models used and the number of communications from the CC to the individual nodes.

Most of the methods were published the big or massive data and multi-machine methodological setting, whereas some were reported within the context of health analytics. Within the big or massive data and multi-machine methodological setting, many methods involved an initial step of random data partitioning among multiple machines. However, certain methods assumed a scenario in which data were already stored on separate machines, as observed in the studies by Fan et al [28] and Jordan et al [13]. Furthermore, it is worth noting that no articles published before 2010 were included, aligning with our initial hypothesis regarding the identification of contemporary methodological settings. Most of the included articles (30/41, 73%) were published after the year 2018.

Most articles (33/41, 81%) addressed a setting in which a CC existed external to the nodes, as exemplified by articles such as those by Lin and Xi [52], Volgushev et al [48], and Yue et al [37]. In contrast, as mentioned previously, some studies (8/41, 19%) designate one of the nodes to assume this central role, as demonstrated in the study by Chang et al [32].

The methods identified through our search strategy shared a common characteristic of using a global model that incorporated population-level parameters. In some cases, these parameters may also include node-specific components to accommodate node-specific statistical heterogeneity in the outcome-predictors relationship, which captures deviations from the population-level conditional probability distribution of the outcome given the predictors.

A few of the reported methods had the capability to yield results identical to those obtained if the individual line data were pooled from all nodes [33,43].

While most articles (31/41, 76%) featured methods related to regression models, including semiparametric and nonparametric designs, a few (10/41, 24%) reported results for other modeling frameworks. These included methods for M-estimation, U-statistics, symmetrical statistics, and natural parameter estimation, some of which encompassed regression models as a specific instance.

Results Related to Objective 2

In total, 6 approaches were selected as applicable to the standard GLM framework. They all assumed that nodes had equal sample sizes and identical distributions for the covariates.

Simple Averaging

One of the simplest methods for horizontally partitioned data analysis, often referred to as the “simple averaging method” or the “divide-and-conquer” approach, has been extensively studied in the literature, for example, the studies by Zhang et al [31] and Shamir et al [59], which were included in our scoping review. This method operates through workflow I in Figure 2. In this approach, node-level model estimates are gathered and averaged at the CC to generate the final estimates.

In the context of GLM, each node is initially tasked with calculating the maximum likelihood estimator (MLE) of the ${\displaystyle {\mathit{\beta }}^{\star }}$ and ${\varphi }^{\star }$ parameters using their respective data. In addition, the Hessian matrix of the log-likelihood function with respect to the β parameters must be computed for constructing Wald-type CIs. The estimated parameters and the computed Hessian matrix are then transmitted to the CC.

The final parameter estimates of are obtained by averaging the node-specific estimates while the local Hessians and estimates of ${\varphi }^{\star }$ are used to compute an estimator for the asymptotic variance.

Single Distributed Newton-Raphson Updating

The single distributed Newton-Raphson updating method is an iterative procedure that includes an additional communication round between the CC and the nodes compared to the simple averaging method. It was originally proposed as the “distributed one-step” method in the study by Beyan et al [60], but in this study, it is referred to by a different term to avoid any confusion regarding communication complexity. This method operates using workflow II, as depicted in Figure 3, with t=1 (where $T$ represents the number of cycles in the iteration scheme). It enhances the simple averaging estimators by incorporating a single distributed Newton-Raphson updating step.

In the context of GLM, each node first calculates the MLE of ${\displaystyle {\mathit{\beta }}^{\star }}$ and ${\varphi }^{\star }$ and transmits them to the CC. The CC aggregates these estimates using averaging and sends the result back to the nodes. The nodes then compute the gradient and the Hessian matrix of the log-likelihood function, evaluated at the received ${\displaystyle {\mathit{\beta }}^{\star }}$ and ${\varphi }^{\star }$ estimates. Subsequently, the gradient and the Hessian matrix are sent back to the CC, which averages them and computes a Newton-Raphson updating step based on the simple averaging estimates. An estimator for the asymptotic variance can be calculated by using the received Hessian matrices and the updated estimate of ${\varphi }^{\star }$.

Multiple Distributed Newton-Raphson Updatings

The multiple distributed Newton-Raphson updating method leverages the fact that, for standard GLMs, the algorithm typically used to calculate the MLE of ${\displaystyle {\mathit{\beta }}^{\star }}$ and ${\varphi }^{\star }$ in a centralized pooled setting can be executed in a distributed manner without any loss of information. This is possible because the algorithm relies on Newton-Raphson updatings (or, sometimes, Fisher scoring updatings) that are expressed using 2 sums of node-specific summary statistics, namely, local gradients and local Hessian matrices of the log-likelihood function, evaluated at the ${\displaystyle {\mathit{\beta }}^{\star }}$ and ${\varphi }^{\star }$ estimates from the previous iteration. A version of this method is proposed in the study by Wu et al [33] under the logistic regression framework. It operates through workflow II in Figure 3 for a general T≥1.

Distributed Estimating Equation

The class of estimating equation estimators is vast and encompasses a broad range of statistical estimation techniques, including likelihood-based approaches that rely on searching for critical points. The fundamental idea behind estimating equation methods is to establish a system of equations that involve both the sample data and the unknown model parameters. These equations are then solved to determine the parameter estimates. MLEs, which are obtained by setting the gradient of the log-likelihood function with respect to the unknown parameters equal to 0, belong to the class of estimating equation estimators.

The distributed estimating equation approach involves gathering summary statistics from nodes at the CC level, enabling the reconstruction of the estimating equations, or, more commonly, an approximation of them that would have been obtained in a pooled centralized setting. This method is discussed in the study by Lin and Xi [30] and operates using workflow I, as depicted in Figure 2.

In the context of GLMs, the distributed estimating equation approach involves initially assigning each node the task of computing and sending their local MLEs and the Hessian matrix of their local log-likelihood, evaluated at those MLEs, to the CC. The CC uses these received quantities to reconstruct the global estimating equations or an approximation thereof. This reconstruction ultimately leads to an analytical solution for obtaining the resulting estimates. CIs are computed using a combination of the Hessian matrices and the final estimator of ${\varphi }^{\star }$.

It is important to note that, when this approach is applied in the context of linear regression, it enables the acquisition of ${\displaystyle {\mathit{\beta }}^{\star }}$parameter estimates that are identical to those obtained in a pooled centralized setting.

Distributed Estimation Using a Single Gradient-Enhanced Log-Likelihood

This method differs fundamentally from the ones discussed thus far as it involves a distinct role for one particular node in obtaining model parameter estimates. It operates using workflow III, as depicted in Figure 4, and was proposed in the study by Jordan et al [13] under the name “Surrogate likelihood.” This approach relies on an approximation of the global likelihood by viewing it as an analytic function. It expands the global likelihood into an infinite series around an initial guess ${\displaystyle {\hat{\mathit{\beta }}}_{\text{SGE},\text{0}}}$ and replaces the higher-order derivatives (order of ≥2) of the global likelihood with those of a Taylor expansion of a node’s (eg, node k=1) local likelihood around the same value. By following this procedure, the so-called surrogate likelihood can be solved using data from node k=1and aggregated gradients from nodes k ∈ {2, ..., K}.

In the context of GLM, the CC first collects the necessary information to compute initial estimates for the parameters ${\displaystyle {\mathit{\beta }}^{\star }}$ and ${\varphi }^{\star }$. These initial estimates can be obtained through various approaches, such as a simple averaging estimator or the MLEs computed using data from node 1. These initial estimates are then transmitted to nodes k ∈ {2, ..., K}. Each of these nodes calculates the gradient of the log-likelihood function, evaluated at the received estimates, and sends it back to the CC. The CC averages these gradients and sends the result to node 1. Node 1 solves a gradient-enhanced log-likelihood using its own data and the received average gradient. The resulting estimate is sent back to the CC as the final estimate. To compute CIs, each node must send the Hessian matrix of its local log-likelihood function, evaluated at the initial received estimate.

The steps related to estimation can be repeated multiple times.

Distributed Estimation Using Multiple Gradient-Enhanced Log-Likelihoods

This method is in the spirit of the distributed estimation using a single gradient-enhanced log-likelihood approach described previously except that all nodes have to solve a gradient-enhanced log-likelihood instead of only one of them. Results pertaining to statistical inference are discussed in the study by Fan et al [28] under a penalized setting. A nonpenalized version of this method was introduced in the study by Shamir et al [59], although the latter did not discuss CIs or hypothesis testing and, hence, was not included in our scoping review. It operates through workflow IV depicted in Figure 5.

The following subsection presents the 6 approaches described in objective 2 within a unified notational framework that accounts for the peculiarities commonly encountered in health analytics. Algorithms were derived using this common notation to rigorously describe the methods and enable their comparison. Theoretical results regarding the estimators involved are detailed and proven in Multimedia Appendix 2. While this section is necessary to increase trust in HPSA methods by transparently showing precisely what information is shared with the CC through comprehensive mathematical formulas, the summary provided in Table 2 suffices for a high-level understanding of the overall picture.

Results Related to Objective 3

Notation for Shared Quantities

In what follows, let the log-likelihood of the data stored in node k (using D^(k)) be denoted by

${\displaystyle {\displaystyle {\mathfrak{l}}^{\left(k\right)}\left(\mathit{\beta },\varphi \right)=\frac{1}{n^{\left(k\right)}}\sum _{i=1}^{n^{\left(k\right)}}\left\{\frac{y_i^{\left(k\right)}h\left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{\left(k\right)}\right)-b\left\{h\left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{\left(k\right)}\right)\right\}}{\varphi }+c\left(y_i^{\left(k\right)},\varphi \right)\right\}.}}$

In addition, let D^(k) (β) ∈ ${\displaystyle \mathbb{R}}$^p+1 be such that

${\displaystyle {\displaystyle {\mathbf{D}}^{\left(k\right)}\left(\mathit{\beta }\right)=\frac{1}{n^{\left(k\right)}}\sum _{i=1}^{n^{\left(k\right)}}{\mathbf{x}}_i^{\left(k\right)}h\prime \left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{\left(k\right)}\right)\left[y_i^{\left(k\right)}-b\prime \left\{h\left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{\left(k\right)}\right)\right\}\right]}}$(2)

and define the (p+1)×(p+1) matrix V^(k) (β) as

$$\begin{gathered} {\displaystyle {\displaystyle \begin{array}{l}{\mathbf{V}}^{(k)}(\mathit{\beta })=\frac{1}{n^{(k)}}\sum _{i=1}^{n^{(k)}}{\mathbf{x}}_i^{(k)}{\left({\mathbf{x}}_i^{(k)}\right)}^{\mathrm{\top }}\left(h^{\mathrm{\prime }}{\left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{(k)}\right)}^2{b}^{\mathrm{\prime }\mathrm{\prime }}\{h\left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{(k)}\right)\} \\ -h^{\mathrm{\prime }\mathrm{\prime }}\left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{(k)}\right)\left[y_i^{(k)}-b^{\mathrm{\prime }}\{h\left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{(k)}\right)\}\right]\right).\end{array}}} \end{gathered}$$(3)

As ${\displaystyle {\mathbf{D}}^{\left(k\right)}\left(\mathit{\beta }\right)=\varphi {\mathrm{\nabla }}_{\mathit{\beta }}{\mathfrak{l}}^{\left(k\right)}\left(\mathit{\beta },\varphi \right)}$, solving the equation ${\displaystyle {\mathbf{D}}^{\left(k\right)}\left(\mathit{\beta }\right)=0}$ yields the node-specific MLE of β, denoted hereafter using ${\displaystyle {\hat{\mathit{\beta }}}_{\text{MLE}}^{\left(k\right)}}$. The matrix ${\displaystyle {\mathbf{V}}^{\left(k\right)}\left(\mathit{\beta }\right)}$ is equal to ${\displaystyle -{\mathrm{\nabla }}_{\mathit{\beta }}{\mathbf{D}}^{\left(k\right)}\left(\mathit{\beta }\right)}$ and relates to the Fisher information matrix through the equation ${\displaystyle {\mathbf{V}}^{\left(k\right)}\left(\mathit{\beta }\right)=-\varphi {\mathrm{\nabla }}_{\mathit{\beta }}^2{\mathfrak{l}}^{\left(k\right)}\left(\mathit{\beta },\varphi \right)}$.

Finally, set

${\displaystyle {\displaystyle E^{\left(k\right)}\left(\varphi ,\mathit{\beta }\right)=\frac{1}{n^{\left(k\right)}}\sum _{i=1}^{n^{\left(k\right)}}\left[y_i^{\left(k\right)}h\left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{\left(k\right)}\right)-b\left\{h\left({\mathit{\beta }}^{\mathrm{\top }}{\mathbf{x}}_i^{\left(k\right)}\right)\right\}\right]-\frac{{\varphi }^2}{n^{\left(k\right)}}\sum _{i=1}^{n^{\left(k\right)}}\frac{\mathrm{\partial }}{\mathrm{\partial }\varphi }c\left(y_i^{\left(k\right)},\varphi \right)}}$(4)

and

${\displaystyle {\displaystyle F^{\left(k\right)}\left(\varphi \right)=\frac{2\varphi }{n^{\left(k\right)}}\sum _{i=1}^{n^{\left(k\right)}}\frac{\mathrm{\partial }}{\mathrm{\partial }\varphi }c\left(y_i^{\left(k\right)},\varphi \right)+\frac{{\varphi }^2}{n^{\left(k\right)}}\sum _{i=1}^{n^{\left(k\right)}}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\varphi }^2}c\left(y_i^{\left(k\right)},\varphi \right).}}$(5)

Because ${\displaystyle E^{\left(k\right)}\left(\varphi ,\mathit{\beta }\right)=-{\varphi }^2\left(\mathrm{\partial }/\mathrm{\partial }\varphi \right){\mathfrak{l}}^{\left(k\right)}\left(\mathit{\beta },\varphi \right)}$, when $\varphi$ is unknown, solving the equation ${\displaystyle E^{\left(k\right)}\left(\varphi ,{\hat{\mathit{\beta }}}_{\text{MLE}}^{\left(k\right)}\right)=0}$ for $\varphi$ yields its node-specific MLE of ${\varphi }^{\star }$. We have ${\displaystyle F^{\left(k\right)}\left(\varphi \right)=-(\mathrm{\partial }/\mathrm{\partial }\varphi )E^{\left(k\right)}\left(\varphi ,\mathit{\beta }\right)}$.

Simple Averaging

The simple averaging method follows upon execution of algorithm 1 (Textbox 1). First, each data node computes their local maximum by solving successively ${\displaystyle {\mathbf{D}}^{\left(k\right)}\left(\mathit{\beta }\right)=0}$ and ${\displaystyle E^{\left(k\right)}\left(\varphi ,{\hat{\mathit{\beta }}}_{\text{MLE}}^{\left(k\right)}\right)=0}$. To compute the CIs at the CC level, the entries of the (p+1)×(p+1) matrix ${\displaystyle {\mathbf{V}}_{\text{MLE}}^{\left(k\right)}={\mathbf{V}}^{\left(k\right)}\left({\hat{\mathit{\beta }}}_{\text{MLE}}^{\left(k\right)}\right)}$ have to be computed from formula 3 with ${\displaystyle \mathit{\beta }={\hat{\mathit{\beta }}}_{\text{MLE}}^{\left(k\right)}}$. Then, the set

${\displaystyle S_0^{\left(k\right)}=\left\{{\hat{\mathit{\beta }}}_{MLE}^{\left(k\right)},{\hat{\varphi }}_{\text{MLE}}^{\left(k\right)},{\mathbf{V}}_{\text{MLE}}^{\left(k\right)}\right\}}$(6)

is sent to the CC. The parameter estimates are then aggregated by the CC through averaging. Specifically, the CC computes

${\displaystyle {\displaystyle {\hat{\mathit{\beta }}}_{\text{SA}}=\sum _{k=1}^K{w}^{\left(k\right)}{\hat{\mathit{\beta }}}_{\text{MLE}}^{\left(k\right)}\text{and}{\hat{\varphi }}_{\text{SA}}=\sum _{k=1}^K{w}^{\left(k\right)}{\hat{\varphi }}_{\text{MLE}}^{\left(k\right)},}}$(7)

where w⁽¹⁾, ..., w^(K) are weights (ie, w^(k)≥0 and ${\displaystyle \sum _{k=1}^K{w}^{\left(k\right)}=1}$) used to combine each node’s contribution. Often, weights can be taken proportional to local sample sizes, leading to the choice w^(k)=n^(k)/n.

Wald-type CIs for ${\displaystyle {\mathit{\beta }}^{\star }}$ can be constructed based on the fact that the sequence ${\displaystyle \sqrt{n}\left({\hat{\mathit{\beta }}}_{\text{SA}}-{\mathit{\beta }}^{\star }\right)}$ converges in distribution to a centered normal random variable with covariance matrix

${\displaystyle {\displaystyle {\mathbf{\Sigma }}_{\text{SA}}={\varphi }^{\star }\sum _{k=1}^K\frac{w^{(k{)}^2}}{p^{(k)}}{\mathfrak{T}}_{{\mathit{\beta }}^{\star }}^{(k)},\text{where}{\mathfrak{T}}_{{\mathit{\beta }}^{\star }}^{(k)}=E\left\{{\mathbf{V}}^k\left({\mathit{\beta }}^{\star }\right)\right\}.}}$

As ${\displaystyle {\mathfrak{T}}_{{\mathit{\beta }}^{\star }}^{\left(k\right)}}$ is consistently estimated using ${\displaystyle {\mathbf{V}}_{\text{MLE}}^{\left(k\right)}}$ and ${\displaystyle {\varphi }^{\star }}$ by ${\displaystyle {\hat{\varphi }}_{\text{SA}}}$, and as p^(k) can be estimated using n^(k)/n, it follows that a consistent estimator for Σ_SA is given by

${\displaystyle {\displaystyle {\hat{\mathbf{\Sigma }}}_{\text{SA}}={\hat{\varphi }}_{\text{SA}}\sum _{k=1}^K\frac{n{w}^{(k{)}^2}}{n^{(k)}}{\mathbf{V}}_{\text{MLE}}^{(k)}.}}$

The simple averaging final estimates are then given by

${\displaystyle R=\left\{{\hat{\mathit{\beta }}}_{\text{SA}},{\hat{\mathbf{\Sigma }}}_{\text{SA}}\right\}.}$

Single Distributed Newton-Raphson Updating

The single distributed Newton-Raphson updating method follows upon execution of algorithm 2 (Textbox 2) with T=1. First, the CC gathers summary statistics to compute the simple averaging estimators of ${\displaystyle {\mathit{\beta }}^{\star }}$ and ${\displaystyle {\varphi }^{\star }}$ without their accompanying CI. Hence, for k∈{1, ..., K}, and with ${\displaystyle {\hat{\mathit{\beta }}}_{\text{MLE}}^{\left(k\right)}}$ as mentioned previously (equation 6), node k sends to the CC the quantities

${\displaystyle S_0^{\left(k\right)}=\left\{{\hat{\mathit{\beta }}}_{\text{MLE}}^{\left(k\right)},{\hat{\varphi }}_{\text{MLE}}^{\left(k\right)}\right\},}$(8)

and the CC uses them to compute ${\displaystyle {\hat{\mathit{\beta }}}_{\text{SA}}}$ and ${\displaystyle {\hat{\varphi }}_{\text{SA}}}$ using the formulas in equation 7.

For reasons of convenience that will become clear later, the notation ${\displaystyle {\hat{\mathit{\beta }}}_{\text{NR},\text{0}}}$ and ${\displaystyle {\hat{\varphi }}_{\text{NR},\text{0}}}$ will be used instead of ${\displaystyle {\hat{\mathit{\beta }}}_{\text{SA}}}$ and ${\displaystyle {\hat{\varphi }}_{\text{SA}}}$, respectively. In this notation, the set of values

${\displaystyle C_1=\left\{{\hat{\mathit{\beta }}}_{\text{NR},\text{0}},{\hat{\varphi }}_{\text{NR},\text{0}}\right\}}$

is broadcasted to data nodes, which are then tasked with computing and sending back the quantities

${\displaystyle \begin{array}{c}{S}_1^{\left(k\right)}=\left\{{\mathbf{D}}_{\text{NR},\text{1}}^{\left(k\right)},{\mathbf{V}}_{\text{NR},\text{1}}^{\left(k\right)},E_{\text{NR},\text{1}}^{\left(k\right)},F_{\text{NR},\text{1}}^{\left(k\right)}\right\}\hspace{0.17em},\end{array}}$

where, for any integer t≥1, one defines

${\displaystyle \begin{array}{l}{\mathbf{D}}_{\mathrm{N}\mathrm{R},\mathrm{t}}^{(k)}={\mathbf{D}}^{(k)}({\hat{\mathit{\beta }}}_{\mathrm{N}\mathrm{R},\mathrm{t}-1})\\ {\mathbf{V}}_{\mathrm{N}\mathrm{R},\mathrm{t}}^{(k)}={\mathbf{V}}^{(k)}({\hat{\mathit{\beta }}}_{\mathrm{N}\mathrm{R},\mathrm{t}-1})\\ E_{\mathrm{N}\mathrm{R},\mathrm{t}}^{(k)}=E^{(k)}({\hat{\varphi }}_{\mathrm{N}\mathrm{R},\mathrm{t}-1},{\hat{\mathit{\beta }}}_{\mathrm{N}\mathrm{R},\mathrm{t}-1})\\ F_{\mathrm{N}\mathrm{R},\mathrm{t}}^{(k)}=F^{(k)}({\hat{\varphi }}_{\mathrm{N}\mathrm{R},\mathrm{t}-1},{\hat{\mathit{\beta }}}_{\mathrm{N}\mathrm{R},\mathrm{t}-1})\end{array}}$(9)

Upon receiving the ${\displaystyle S_1^{\left(k\right)}}$s from each node, the CC calculates the following weighted averages:

${\displaystyle {\displaystyle \begin{array}{cc}& {\overline{\mathbf{D}}}_{\text{NR},\text{1}}=\sum _{k=1}^K{w}^{\left(k\right)}{\mathbf{D}}_{\text{NR},\text{1}}^{\left(k\right)},\text{and}{\overline{\mathbf{V}}}_{\text{NR},\text{1}}=\sum _{k=1}^K{w}^{\left(k\right)}{\mathbf{V}}_{\text{NR},\text{1}}^{\left(k\right)},\\ & {\overline{E}}_{\text{NR},\text{1}}=\sum _{k=1}^K{w}^{\left(k\right)}{E}_{\text{NR},\text{1}}^{\left(k\right)},\text{and}{\overline{F}}_{\text{NR},\text{1}}=\sum _{k=1}^K{w}^{\left(k\right)}{F}_{\text{NR},\text{1}}^{\left(k\right)}.\end{array}}}$

This enables the CC to execute Newton-Raphson updates from ${\displaystyle {\hat{\mathit{\beta }}}_{\text{NR},\text{0}}}$ and ${\displaystyle {\hat{\varphi }}_{\text{NR},\text{0}}}$, respectively:

${\displaystyle \begin{array}{l}{\hat{\mathit{\beta }}}_{\mathrm{N}\mathrm{R},1}={\hat{\mathit{\beta }}}_{\mathrm{N}\mathrm{R},0}+{\overline{\mathbf{V}}}_{\mathrm{N}\mathrm{R},1}^{-1}{\overline{\mathbf{D}}}_{\mathrm{N}\mathrm{R},1}\\ {\hat{\varphi }}_{\mathrm{N}\mathrm{R},1}={\hat{\varphi }}_{\mathrm{N}\mathrm{R},0}+{\overline{F}}_{\mathrm{N}\mathrm{R},1}^{-1}{\overline{E}}_{\mathrm{N}\mathrm{R},1}.\end{array}}$(10)

It is shown in the Multimedia Appendix 2 that

${\displaystyle {\displaystyle \begin{array}{c}\sqrt{n}\left({\hat{\mathit{\beta }}}_{\text{NR},1}-{\mathit{\beta }}^{\star }\right)\to \mathcal{N}\left(0,{\mathbf{\Sigma }}_{\text{NR}}\right)\\ \text{where}{\mathbf{\Sigma }}_{\text{NR}}={\left(\sum _{k=1}^K{w}^{(k)}{\mathfrak{T}}_{{\mathit{\beta }}^{\star }}^{(k)}\right)}^{-1}\left\{{\varphi }^{\star }\sum _{k=1}^K\frac{w^{(k{)}^2}}{p^{(k)}}{\mathfrak{T}}_{{\mathit{\beta }}^{\star }}^{(k)}\right\}{\left(\sum _{k=1}^K{w}^{(k)}{\mathfrak{T}}_{{\mathit{\beta }}^{\star }}^{(k)}\right)}^{-1}.\end{array}}}$

As ${\displaystyle {\mathfrak{T}}_{{\mathit{\beta }}^{\star }}^{\left(k\right)}}$ is consistently estimated using ${\displaystyle {\mathbf{V}}_{\text{NR},\text{1}}^{\left(k\right)}}$ and ${\varphi }^{\star }$ by ${\displaystyle {\hat{\varphi }}_{\text{NR},\text{1}}}$, and as ${\displaystyle p^{\left(k\right)}}$ can be estimated using ${\displaystyle n^{\left(k\right)}/n}$, it follows that a consistent estimator for ${\displaystyle {\mathbf{\Sigma }}_{\text{NR}}}$ is given by

${\displaystyle {\displaystyle \begin{array}{c}{\hat{\mathbf{\Sigma }}}_{\text{NR}}={\left(\sum _{k=1}^K{w}^{(k)}{\mathbf{V}}_{\text{NR},1}^{(k)}\right)}^{-1}\left\{{\hat{\varphi }}_{\text{NR},1}\sum _{k=1}^K\frac{n{w}^{(k{)}^2}}{n^{(k)}}{\mathbf{V}}_{\text{NR},1}^{(k)}\right\}{\left(\sum _{k=1}^K{w}^{(k)}{\mathbf{V}}_{\text{NR},1}^{(k)}\right)}^{-1}\\ ={\left({\overline{\mathbf{V}}}_{\text{NR},1}\right)}^{-1}\left\{{\hat{\varphi }}_{\text{NR},1}\sum _{k=1}^K\frac{n{w}^{(k{)}^2}}{n^{(k)}}{\mathbf{V}}_{\text{NR},1}^{(k)}\right\}{\left({\overline{\mathbf{V}}}_{\text{NR},1}\right)}^{-1}.\end{array}}}$

The method’s final estimates are then given by

${\displaystyle R=\left\{{\hat{\mathit{\beta }}}_{\text{NR},\text{1}},{\hat{\mathbf{\Sigma }}}_{\text{NR}}\right\}.}$

Multiple Distributed Newton-Raphson Updatings

The multiple distributed Newton-Raphson updatings method follows upon execution of algorithm 2 with T>1.

The first communication cycle follows the same procedure described previously for the single distributed Newton-Raphson updating method. It involves distributively computing a simple averaging estimator and then performing a Newton-Raphson iteration starting from this estimator. The Newton descent is calculated as described in equation 10.

Formally, the algorithm begins with each data node k sending the set of quantities ${\displaystyle S_0^{\left(k\right)}}$, as described in equation 8, to the CC. Next, the CC calculates the simple averaging estimators using formula 7 and uses them to initialize ${\displaystyle {\hat{\mathit{\beta }}}_{\text{NR,Step=0}}={\hat{\mathit{\beta }}}_{\text{SA}}}$ and ${\displaystyle {\hat{\varphi }}_{\text{NR,Step=0}}={\hat{\varphi }}_{\text{SA}}}$.

The following steps are then repeated for a certain number of iterations. At iteration $t$, starting from $t=1$, the CC broadcasts the values ${\displaystyle C_t=\left({\hat{\mathit{\beta }}}_{\text{NR},\text{t-1}},{\hat{\varphi }}_{\text{NR},\text{t-1}}\right)}$ to the data nodes. The data nodes compute the quantities ${\displaystyle {\mathbf{D}}_{\text{NR},\text{t}}^{\left(k\right)}}$, ${\displaystyle {\mathbf{V}}_{\text{NR},\text{t}}^{\left(k\right)}}$, ${\displaystyle E_{\text{NR},\text{t}}^{\left(k\right)}}$, and ${\displaystyle F_{\text{NR},\text{t}}^{\left(k\right)}}$ as defined in equation 9 and send them back to the CC.

The CC then uses these quantities to perform a Newton update. Specifically, it calculates ${\displaystyle {\hat{\mathit{\beta }}}_{\text{NR},\text{t}}={\hat{\mathit{\beta }}}_{\text{NR},\text{t-1}}+{\mathbf{V}}_{\text{NR},\text{t}}^{-1}{\mathbf{D}}_{\text{NR},\text{t}}}$ and ${\displaystyle {\hat{\varphi }}_{\text{NR},\text{t}}={\hat{\varphi }}_{\text{NR},\text{t-1}}+E_{\text{NR},\text{t}}/F_{\text{NR},\text{t}}}$.

If the iterative cycle is repeated until convergence, the resulting estimates of ${\mathbf{\beta }}^{\star }$ are equivalent to the MLEs derived from pooled data. This is because, in GLMs, for MLEs, if both the pooled and distributed algorithms are initialized using the same values for ${\displaystyle {\hat{\mathit{\beta }}}_{\text{NR,Step=0}}}$ and ${\displaystyle {\hat{\varphi }}_{\text{NR,Step=0}}}$, then, at each subsequent iteration, the distributed Newton update computed by the CC will be identical to the update obtained in a pooled setting.

For the method to yield consistent estimates, it is not necessary to initialize it using simple averaging estimators. However, using simple averaging estimators as initialization may speed up convergence, particularly in large sample sizes, as these estimators are ${\displaystyle \sqrt{n}}$-consistent.

Let ${\displaystyle {\hat{\mathit{\beta }}}_{\text{NR},\text{T}}}$ denote the estimator obtained at convergence. As it is (nearly) equal to the pooled MLE of ${\displaystyle {\mathit{\beta }}^{\star }}$, we can deduce from the Multimedia Appendix 2 that

${\displaystyle {\displaystyle \begin{array}{c}\sqrt{n}\left({\hat{\mathit{\beta }}}_{\text{NR},T}-{\mathit{\beta }}^{\star }\right)\to \mathcal{N}\left(0,{\mathbf{\Sigma }}_{\text{NR}}\right)\\ \text{where}{\mathbf{\Sigma }}_{\text{NR}}={\left(\sum _{k=1}^K{w}^{(k)}{\mathfrak{T}}_{{\mathit{\beta }}^{\star }}^{(k)}\right)}^{-1}\left\{{\varphi }^{\star }\sum _{k=1}^K\frac{w^{(k{)}^2}}{p^{(k)}}{\mathfrak{T}}_{{\mathit{\beta }}^{\star }}^{(k)}\right\}{\left(\sum _{k=1}^K{w}^{(k)}{\mathfrak{T}}_{{\mathit{\beta }}^{\star }}^{(k)}\right)}^{-1}.\end{array}}}$(11)

Following the same reasoning used previously for the single distributed Newton-Raphson updating method, we can consistently estimate the variance-covariance matrix as

${\displaystyle {\displaystyle {\hat{\mathbf{\Sigma }}}_{\text{NR}}={\left({\overline{\mathbf{V}}}_{\text{NR},T}\right)}^{-1}\left\{{\hat{\varphi }}_{\text{NR},T}\sum _{k=1}^K\frac{n{w}^{(k{)}^2}}{n^{(k)}}{\mathbf{V}}_{\text{NR},T}^{(k)}\right\}{\left({\overline{\mathbf{V}}}_{\text{NR},T}\right)}^{-1}.}}$

Distributed Estimating Equation

The distributed estimating equation method follows upon execution of algorithm 3 (Textbox 3). First, each node is responsible for computing the MLEs ${\displaystyle {\hat{\mathit{\beta }}}_{\text{MLE}}^{\left(k\right)}}$ and ${\displaystyle {\hat{\varphi }}_{\text{MLE}}^{\left(k\right)}}$ of ${\displaystyle {\mathit{\beta }}^{\star }}$ and ${\varphi }^{\star }$, respectively, using its own data. These estimators, along with the Hessian matrix ${\displaystyle {\mathbf{V}}_{\text{MLE}}^{\left(k\right)}={\mathbf{V}}^{\left(k\right)}\left({\hat{\mathit{\beta }}}_{\text{MLE}}^{\left(k\right)}\right)}$ and ${\displaystyle F_{\text{MLE}}^{(k)}=F^{(k)}\left({\hat{\varphi }}_{\text{MLE}}^{(k)}\right)}$, are then sent to the CC. The set

${\displaystyle \begin{array}{c}{S}_0^{\left(k\right)}=\left\{{\hat{\mathit{\beta }}}_{\text{MLE}}^{\left(k\right)},{\hat{\varphi }}_{\text{MLE}}^{\left(k\right)},{\mathbf{V}}_{\text{MLE}}^{\left(k\right)},F_{\text{MLE}}^{\left(k\right)}\right\}\end{array}}$(12)

is transmitted to the CC. The CC calculates the weighted average of the Hessians and the ${\displaystyle F_{\text{MLE}}^{\left(k\right)}}$ values as follows:

${\displaystyle {\displaystyle \begin{array}{c}{\overline{\mathbf{V}}}_{\text{EE}}=\sum _{k=1}^K{w}^{\left(k\right)}{\mathbf{V}}_{\text{MLE}}^{\left(k\right)}\text{and}{\overline{F}}_{\text{EE}}=\sum _{k=1}^K{w}^{\left(k\right)}{F}_{\text{MLE}}^{\left(k\right)}.\end{array}}}$(13)

The parameter estimates can then be calculated as

${\displaystyle {\displaystyle \begin{array}{c}{\hat{\mathit{\beta }}}_{\text{EE}}={\overline{\mathbf{V}}}_{\text{EE}}^{-1}\sum _{k=1}^K{w}^{\left(k\right)}{\mathbf{V}}_{\text{MLE}}^{\left(k\right)}{\hat{\mathit{\beta }}}_{\text{MLE}}^{\left(k\right)}{\hat{\varphi }}_{\text{EE}}={\overline{F}}_{\text{EE}}^{-1}\sum _{k=1}^K{w}^{\left(k\right)}{F}_{\text{MLE}}^{\left(k\right)}{\hat{\varphi }}_{\text{MLE}}^{\left(k\right)}.\end{array}}}$(14)

It is shown in the Multimedia Appendix 2 that ${\displaystyle \sqrt{n}\left({\hat{\mathit{\beta }}}_{\text{EE}}-{\mathit{\beta }}^{\star }\right)}$ converges in distribution to a centered normal random variable with the variance-covariance matrix given by

${\displaystyle {\displaystyle {\mathbf{\Sigma }}_{\text{EE}}={\left(\sum _{k=1}^K{w}^{(k)}{\mathfrak{T}}_{{\mathit{\beta }}^{\star }}^{(k)}\right)}^{-1}\left\{{\varphi }^{\star }\sum _{k=1}^K\frac{w^{(k)2}}{p^{(k)}}{\mathfrak{T}}_{{\mathit{\beta }}^{\star }}^{(k)}\right\}{\left(\sum _{k=1}^K{w}^{(k)}{\mathfrak{T}}_{{\mathit{\beta }}^{\star }}^{(k)}\right)}^{-1}}}$

It can be consistently estimated by

${\displaystyle {\displaystyle {\hat{\mathbf{\Sigma }}}_{\text{EE}}={\left({\overline{\mathbf{V}}}_{\text{EE}}\right)}^{-1}\left\{{\hat{\varphi }}_{\text{EE}}\sum _{k=1}^k\frac{n{w}^{(k{)}^2}}{n^{(k)}}{\mathbf{V}}_{\text{MLE}}^{(k)}\right\}{\left({\overline{\mathbf{V}}}_{\text{EE}}\right)}^{-1}}}$