The process proceeds (arbitrarily) from small to large values. The first bin is, initially, from a₀ to a₁, where a₀ is the minimal value in the data and a₁, is the smallest data value for which [a_o, a₁] has at least K observations from one group and either 0 or at least K observations from the other group. The next tentative bin limit, a₂, is found in the same way, looking at the interval (a₁, a₂]. This continues so long as a new bin limit can be found. When a limit cannot be found, the number of unbinned data in at least one group is between 0 and K. Tentatively extend the upper limit of the previously formed bin to the maximal value of this group as the next limit. The unbinned data from the other group might permit continuation of the process, blocking off new bins in which that group has counts of at least K, vs. counts of 0 for the first group. When that group has fewer than K unbinned data, replace the last bin limit by the maximal value in the second group (See Algorithm S1 in Supplementary material for details.).

The initial bin boundaries a₀, …, a_B produced by the algorithm above are actual data values and, unless many subjects share the same value, violate the privacy condition. There is a simple fix for a₁, …, a_B−1. All values in the jth bin are ≤ a_j and all values in the j+1st bin are >a_j. So we can replace a_j by c_j = wa_j+(1−w)v_j+1 where v_j+1 is the smallest value in the (j+1)st bin and w is a uniform random variable on (0, 1). The extreme boundaries a₀ and a_B are the minimum and maximum in the data, so a different approach is needed. One option is to take c₀ = −∞ and c_B = ∞. Another option is to impose natural limits; for example, if by definition a variable cannot assume negative values, we could choose c₀ = 0. A final option is to extend the bin limits by “privacy buffers”. To make these reasonably close to the data, we base them on the observed gaps between successive observations in the extreme bin. For example, compute c_B as aB+d̄B, where d̄B is the mean difference between consecutive data points in the last bin. (If d̄B=0, c_B = a_B, but this is now privacy preserving, as all observations in the last bin are equal to one another, with more than K in each group that has data.) Similarly, compute c₀ as a0-d̄1.

A new algorithm is needed to add the data from a new center, preserving all bin boundaries from the first center. The simple option of increasing the frequency counts in each current bin is not an option, as the incremental table from the new center will typically not be K-anonymous. Further, the incremental counts for some existing bin might be so large that data from the new center could actually be used to split it into two or more bins.

Algorithm S3 in Supplementary material is used to add the information from a new center to an existing summary table. We first iterate over the current bins, creating finer bins if possible. Then we remove any counts that are not K-anonymous by combining and redistributing data from adjacent cells. Pseudo-code for Algorithm S3 and for two algorithms called by it are given in the Supplementary material.

Splitting an existing bin into two bins forces us to reallocate the previous frequencies. We do so proportionally to the relative frequencies from the new center. For example, suppose a bin with a current count of 27 for one group is split into two new bins, which have equal counts at the new center. Then we split the 27 equally to the two new groups, adding 13.5 to each. Note that this procedure can result in counts that are not integers.

After creating new bins wherever possible, we iterate again and fix bins where the new center has frequencies between 0 and K. Proceeding from bin 1 to bin B, these non-private bins are combined with the next bin to the right until all counts from the new center are either 0 or at least K. Then the total counts are distributed among the original bins proportionally to the relative frequencies of the bins in the current table. Supplementary Table 1 shows an example that illustrates how the algorithm works.

The extreme bin limits c₀ and c_B must be compared with the minimum and maximum values, respectively, in the new center. If the new center has a more extreme data value, we need to revise these bin limits. We do so by applying the buffer method that was used to find c₀ and c_B in the largest center, but now adding buffers that depend only on the data in the extreme bin from the new center.

The MWU statistic can be defined as follows. Denoting the observations in the two groups by X₁, …, X_n and Y₁, …, Y_m,

with

If H₀ is true, the expected value of U is 0 and its variance is V=mn(N+1)3[1-∑r=1D(tr3-tr)N(N2-1)], where N = n+m, D is the number of distinct values in the data, and t_r is the number of observations that share the rth distinct value. The second term corrects the variance for the presence of ties in the data. If Y=dc+X, c∈ℝ, the distribution of U is stochastically increasing as a function of c. The power of the test depends on P(Y>X) and is high when this probability differs from 0.5.

The MWU test involves direct comparison of each data point in one group with each data point from the other group. As this includes comparisons of observations from different centers, it is impossible to compute the MWU statistic for a federated analysis. Two broad options are possible for federated analysis.

The next subsections present options for combining center-specific MWU statistics and the second analysis option, used in conjunction with our federated binning algorithm.

Denote by U_l the MWU from the lth center, based on n_l and m_l observations from the two groups, with N_l = n_l+m_l; and denote by V_l its variance under H₀. A simple way to form a federated test statistic is to sum the individual statistics over the centers and normalize them by their standard deviation, leading to

A simple generalization is to replace the sum of the statistics by a weighted sum, with an optimal choice of weights. It is convenient to do this using the normalized test statistics for each center, Zl=Ul/Vl0.5. The weighted test statistic is then

The choice of weights can be made to maximize the power of the test when the null hypothesis is not true, using the fact that

where δ_l is the standardized effect in center l. For the MWU statistic, the standardized effect can be expressed as

where Pl+=P(Y>X) and Pl-=P(Y<X) in center l. Although the formula permits the probability difference to vary over centers, the natural basis for defining the weighted sum statistic is to assume a constant difference, in which case the optimal weights depend on the sample sizes and, if present, the extent of tied data. See Equation S1 in the Supplementary material for derivation of the weights.

Fisher's method [23] combines the p-values from independent samples. The corresponding statistic is TF=-2∑l=1Llog(pl)→H0χ2L2 where p_l is the p-value from the MWU test result in the lth center.

We can compute the MWU statistic from the federated summary table generated by the algorithm described in Section 3. The table will have B bins whose frequencies are fx_i and fy_i. The frequencies sum to the total amount of data over all the centers, but need not be integers.

The MWU statistic for the federated table compares observations on the basis of their bins and is given by

where c₀<c₁ < ⋯ < c_B are the endpoints of the bins and

The variance of U_fed can be computed from the formula in Section 4.1, keeping in mind that all observations in the same bin are tied.

A simulation study was used to compare the different federated MWU tests to an analysis of the combined data. Our goals are to assess how the federated analysis affects the power of the tests, and to use the power analysis to compare the testing methods. We also vary the simulation settings to examine how the results and comparisons are affected by the number of centers in the study and by heterogeneity across centers.

We simulated situations with 1,500 observations in each group, divided over 3, 5, or 10 centers, with the number of observations unbalanced among the centers (see Table 1).

For the Mann-Whitney test, only the order of the observations is important, so any distribution can be used to simulate the data. Our model generates control group observations at center l as

and treatment group observations as

The possibility that centers may differ from one another is represented by αl~N(0,σα2). The difference between treatment and control at center l is βl~N(δ,σβ2) where δ is the overall difference, and σ_β represents heterogeneity of the treatment effect across centers. The terms ϵ_il, ϵ_jl~N(0, 1) are random errors. All random variables are independent of one another.

We simulated experiments with several different combinations of input parameters. We chose σ_α∈{0, 0.1, 0.2} and σ_β∈{0, 0.05, 0.06} to achieve between center variance, and δ∈{0, 0.05, 0.1}. Including δ = 0 allowed us to verify that the tests remain reliable when both groups have the same mean. Note, however, that the variance is slightly larger for the treatment group if σ_β>0, so that this setting does not fully match the null hypothesis of identical distributions.

Supplementary Figure 1 shows the distributions of p-values for all the tests in the null setting δ = 0. The left panel includes heterogeneity across centers (σ_α = 0.1), but no effect heterogeneity, and shows a uniform distribution for all the tests, as desired. The right panel adds a small amount of effect heterogeneity (σ_β = 0.05. This results in a slightly wider spread of p-values for all the tests, so that actual type 1 errors are inflated from their nominal values. The fraction of p-values below 0.05 (0.01) was approximately 0.08 (0.025). The inflation was slightly weaker when more centers were included and slightly larger only for Fisher's test. The additional bias of Fisher's test is not surprising, as it is sensitive to the existence of an effect within a center, but not to having a consistent direction of the effect.

Figure 1 compares methods when δ≠0 across different parameters and numbers of centers. See also Supplementary Table S1. The federated table and weighted tests have p-value distributions that are very similar to those from combining all the data, indicating almost no loss of power. The sum test has higher p-values, hence consistently lower power. The p-values with Fisher's method are a bit higher when the treatment effect is consistent across centers (σ_β = 0). When the effect is not consistent, they are lower. However, as already seen, Fisher's test in this case fails to preserve type 1 error, with a bias toward low values.

Figure 2 focuses on how closely the federated test results compare with those from the combined test (i.e., using the full data) by comparing the p-values of each method on the same simulated data set. The Y axis presents log(p_iv/p_is) where i represents the simulation number, s is the combined test and v is the federated test. A federated test that produces the same p-values as the combined test has no loss of power. The more tightly concentrated are these distributions around 0, the more nearly identical are the p-values of the federated method to those of the combined method.

Across all the settings, the weighted test most closely replicates the p-value of the combined test. The federated table is also similar, but more variable, especially when δ≠0. In the top left panel, where H₀ is true, all methods are similar to the combined test. However, adding treatment heterogeneity (top right panel) induces bias in the p-values from Fisher's test, which will reject the null hypothesis too often. It also increases the variance of the log ratio for that test and for the sum. In all the settings with center heterogeneity (σ_α>0), the sum test gave, typically, slightly higher p-values than the combined test, hence had lower power.

To assess the power of the tests as a function of the effect size, we simulated p-values over a set of 4 increasing values of δ, when σ_α = 0.1 and σ_β = 0.05. Figure 3 compares the methods to the unconstrained test using log(p_iv/p_is) (Y-axis) where i represents the simulation number, s is the unconstrained method and v is the other method. Again the weighted test is most similar to the combined test, followed by the federated table. Table 2 shows the median of the p-value distributions with 10 centers; smaller medians correspond to higher power for the test. The medians for the weighted test are consistently the lowest ones; with even the modest heterogeneity present here, they are lower even than those for the combined test. The test from our federated table has slightly higher medians throughout than does the combined test. Similar quantiles were found for 3 and for 5 centers, indicating that, for the settings we examined, the number of centers has little effect on power.

A quantile can be estimated as the solution to a minimization problem,

where the target function is the quantile loss function. The optimization can be carried out on federated data by returning function and gradient values from each center, proceeding iteratively to compute Q^p,Loss. The need for an iterative algorithm to minimize the loss, has the drawback of communication inefficiency.

A more serious concern is that the quantile loss compromises privacy. The loss function within each center is piecewise linear with a change in derivative at each data value in the center. Thus the information from a collection of calls can be used to recover the original data values at the federated node.

Despite the privacy violation, we will include Q^p,Loss in the subsequent comparisons as a benchmark.

It is possible to exploit the loss function to compute approximate quantile estimators that are differentially private [24].

The binning algorithm we introduced in Section 3 can be used to compute a federated estimate of Q_p that is K-anonymous. Here, we apply the single group version of the algorithm which gives a summary table that has B bins with endpoints b₀<b₁<b₂ < ⋯ < b_B and frequencies f_{x, k}. Let F^x,i denote the cumulative distribution for the federated table at b_i.

A naive estimate is the smallest bin limit with cumulative frequency greater than 100p% of the data. However, restricting Q_p to the set of bin limits is an obvious drawback, especially for quantiles in the tails of the distribution. A simple improvement is to interpolate the estimated CDF from one bin limit to the next. Linear interpolation corresponds to the assumption of a uniform distribution within each bin. That may be reasonable for bins in the center of the data. However, it is not likely to work well in the tails, especially in the most extreme bins. We did attempt to use linear interpolation, but the results were poor and are not reported here.

We propose here a more sophisticated interpolation method based on the Yeo-Johnson transformation (“YJ”) [25], a power transformation used to achieve a distribution that is closer to the normal. The approach extends the well-known Box-Cox [26] transformation to also handle variables that can take on negative values. The transformation is defined by

Federated quantile estimates can be used to generate an alternative summary table, which presents a collection of quantiles. See Table 3 for an example, with estimates from optimizing the quantile loss and the YJ likelihood.

We compared the three quantile estimators using a simulation configuration similar to that in the testing chapter. As the quantiles are univariate summaries, we generated data and estimated quantiles only in one group. Another important difference is that the form of the underlying distribution affects the estimation results. In particular, methods may vary when faced with long rather than short tails. To gain insight into this issue, we chose the Gamma as the base distribution for assessing the quality of quantile estimation.

Each simulated data set included 1,500 observations, spread across 3, 5, or 10 centers exactly as described in Table 1. The observations were generated from the following model: x_il = ϵ_ilexp(α_l) where x_il is observation i at center l with αl~N(0,σα2) and ϵ_il~Gamma(r, 1) r∈{4, 10}. The skewness of Gamma is 2r, so the smaller value for r has a longer right tail.

For the Gamma data, heterogeneity across centers was induced using scale rather than location shifts. The value of σ_α was chosen to achieve between center heterogeneity similar in extent to that in Section 4. There the key term was the ratio σ_α/σ_ϵ, which was taken to be 0, 0.1 or 0.2. With Gamma data, the standard deviation of the homogeneous data is proportional to the median, so the analogous choice is to set σα=log(Q0.5+ϕσϵQ0.5), with ϕ similar to the values chosen above. We used only ϕ = 0.1 in our simulations for quantile estimation.

For each combination of the parameters, 2,000 simulations were run. The true quantile Q_p for each simulation was computed from the mixture (over centers) distribution by solving the equation below with l as the center index.

where N is the number of observations from all centers, n_l the observations in center l and Γ_r is the standard Gamma CDF with shape parameter r. A dominant part of the quantile estimation errors is the natural variability of the underlying Gamma distribution. As the standard deviation for Gamma(r, 1) is r, we summarized results via the normalized estimation error (Q^p-Qp)/r where Q^p is the estimator of Q_p.

The simulation results for estimating Q_0.98 are shown in Figure 4. This quantile is presented separately, as it is the most challenging case, in the right tail of a right-skewed distribution. Results for Q_0.02,Q_0.25,Q_0.5, Q_0.75 are depicted in Figure 5. Further detail is provided in Supplementary Tables S2–S5, which give, respectively, the estimated bias and standard deviation, the mean squared error (MSE), and the ratio of squared bias to variance for all the methods and all the quantiles.

The YJ data estimator achieved lower MSE than the quantile loss estimator. For the extreme quantiles, the decrease in MSE ranged from 14% to 44%. The “one pass” YJ table estimator was very accurate for estimating the median and the quartiles, but lost efficiency for the extreme quantiles with the more skewed of the two Gamma distributions and when the number of centers was large. In that setting, the estimator for Q_0.02 suffered from negative bias and its MSE was almost 3 times as large as for the quantile loss estimator; the MSE for Q_0.98 was about 80% larger.

For the settings we studied, variance was the dominant component of MSE. Bias was a substantial problem only in a small number of cases. The YJ methods had large positive bias for Q_0.98 when r = 4; however, when r = 10, and the distribution is itself closer to normal, the bias was negligible.