Improved Differentially Private Analysis of Variance

TL;DR

This paper introduces a new differentially private ANOVA test with higher statistical power, requiring significantly less data to detect effects compared to previous methods.

Contribution

It develops a novel test statistic $F_1$ for private ANOVA, providing rigorous reference distribution computation and improved data efficiency.

Findings

The $F_1$ statistic outperforms the traditional $F$-statistic in private settings.

The new test achieves an order of magnitude better data efficiency.

Experimental results show only 7% of the data needed compared to previous methods.

Abstract

Hypothesis testing is one of the most common types of data analysis and forms the backbone of scientific research in many disciplines. Analysis of variance (ANOVA) in particular is used to detect dependence between a categorical and a numerical variable. Here we show how one can carry out this hypothesis test under the restrictions of differential privacy. We show that the $F$-statistic, the optimal test statistic in the public setting, is no longer optimal in the private setting, and we develop a new test statistic $F_1$ with much higher statistical power. We show how to rigorously compute a reference distribution for the $F_1$ statistic and give an algorithm that outputs accurate $p$-values. We implement our test and experimentally optimize several parameters. We then compare our test to the only previous work on private ANOVA testing, using the same effect size as that work. We see an order of magnitude improvement, with our test requiring only 7% as much data to detect the effect.

Equations

$$\mathrm{Pr}[T\ge t\mid T=f(\mathbf{X})\text{ and }\mathbf{X}\leftarrow H_0].$$

$$\mathrm{Pr}[T\ge t^{\ast }\mid T=f(\mathbf{X})\text{ and }\mathbf{X}\leftarrow H_A]$$

$$\text{SSA}(\mathbf{x})=\sum _{j=1}^k{n}_j({\bar{y}}_j-\bar{y}{)}^2.$$

$$\text{SSE}(\mathbf{x})=\sum _{i=1}^N(y_i-{\bar{y}}_{c_i}{)}^2.$$

$$F(\mathbf{x})=\frac{\text{SSA}(\mathbf{x})/(k-1)}{\text{SSE}(\mathbf{x})/(N-k)}=\frac{\text{MSA}(\mathbf{x})}{\text{MSE}(\mathbf{x})}.$$

$$\text{Pr}[f(\mathbf{x})\in S]\le e^{\epsilon }\text{Pr}[f({\mathbf{x}}^{\prime })\in S].$$

$${\sf Lap}(z\mid b)=\frac{1}{2b}\text{\emph{exp}}\bigg{(}{-\frac{\lvert z\rvert}{b}}\bigg{)}.$$

$$\text{SA}(\mathbf{x})$$

$$\text{SE}(\mathbf{x})$$

$$F_1(\mathbf{x})=\frac{\text{SA}(\mathbf{x})/(k-1)}{\text{SE}(\mathbf{x})/(N-k)}.$$

$$\text{SE}(\mathbf{x})=\sum _{j=1}^k\sum _{i\in C_j}|y_i-{\bar{y}}_{c_i}|.$$

$$\frac{n_a}{n_a+1}+\frac{n_b}{n_b+1}\le 2$$

$$\text{SA}(\mathbf{x})=\sum _{j=1}^k\sum _{i\in C_j}|\overline{y}-{\overline{y}}_{c_i}|$$

$$\hat{\sigma }=\frac{\widehat{SE}}{\tilde{N}}\sqrt{\frac{\pi }{2}},\text{ where }\tilde{N}=\sum _{j=1}^k{n}_j\sqrt{\left(1-\frac{1}{n_j}\right)}$$

$$\mathrm{Pr}[p\text{-value}\le \alpha \mid \mathbf{X}\leftarrow H_0]\le \alpha .$$

$$\text{SQA}(\mathbf{x})=\sum _{j=1}^k{n}_j{\left|{\overline{y}}_j-\overline{y}\right|}^q$$

$$\text{SQE}(\mathbf{x})=\sum _{i=1}^N{\left|y_i-{\overline{y}}_{c_i}\right|}^q$$

$$F_q(\mathbf{x})=\frac{\text{SQA}/(k-1)}{\text{SQE}/(N-k)}$$

$$2\bigg{(}\frac{N}{2}\bigg{)}^{(1-q)}+1$$

$$N-N\bigg{(}1-\frac{2}{N}\bigg{)}^{q}+1$$

$$N\bigg{(}\frac{3}{N}\bigg{)}^{q}+1$$

$$N-N\bigg{(}1-\frac{3}{N}\bigg{)}^{q}+1$$

$$\widehat{SE}=\left(\sum _{j=1}^k\sum _{i\in C_j}|y_i-{\bar{y}}_{c_i}|\right)+L_2$$

$$y_i$$

$${\bar{y}}_{c_i}$$

$$y_i-{\bar{y}}_{c_i}$$

$${\tau }^2$$

$$=Var(y_i)+Var({\bar{y}}_{c_i})-2Cov(y_i,{\bar{y}}_{c_i})$$

$$={\sigma }^2+\frac{{\sigma }^2}{n_{c_i}}-2Cov(y_i,{\bar{y}}_{c_i}).$$

$$E(y_i{\bar{y}}_{c_i})$$

$$=\frac{1}{n_{c_i}}E(y_i^2+y_i{S}_{-i})$$

$$=\frac{1}{n_{c_i}}\left(E\left(y_i^2\right)+E\left(y_i{S}_{-i}\right)\right),$$

$$Z^2={\left(\frac{y_i-\mu }{\sigma }\right)}^2=\frac{1}{{\sigma }^2}\left(y^2-2y_i\mu +{\mu }^2\right).$$

$$E(y^2)$$

Full text

1 Introduction

A universal and recurring challenge in scientific research is determining whether a measured effect is real. That is, researchers wish to determine if the effect observed in a particular dataset indicates a similar effect in the broader world from which the sample was drawn. The most common statistical tool to make this determination is a hypothesis test. The particular form of the hypothesis test is driven by the scientific question and the data at hand.

Hypothesis testing is a common tool in population association studies, where the goal is to identify whether genetic variation is associated with disease risk balding2006tutorial . Consider a study looking at the effect of a particular gene’s mutation on some health outcome (e.g., blood pressure or weight). The data may include the mutation status of that gene (which may harbor one or more mutations on one or both copies of DNA). The researcher’s goal is to determine if the gene’s mutation status has an impact on that health outcome. When the health outcome is measured using a numerical variable, the natural hypothesis test to use is the one-way ANOVA (analysis of variance) test, treating the gene’s mutation status as a categorical variable.

The first step in conducting the ANOVA is to calculate a single number, the $F$-statistic, which measures the variation in group means compared to the variation in individual data points. The $F$-statistic is constructed so that if the expected value of the health outcome is the same in all groups, the expected value of $F$ is 1. If this is not the case, the value can be dramatically larger. Seeing a high value of $F$, a researcher will compare this value to the reference distribution, i.e., the distribution of $F$ that would occur if the gene had no impact on the health outcome. The result of this comparison is a $p$-value, the probability that the observed $F$ could occur by chance. If the $p$-value is low, the analyst can conclude that this gene must indeed affect the given health outcome. (For more detail on how ANOVA is used in this setting, see myers2003researchdesign .)

The analysis described above assumes that the researcher has full access to the database. However, there are many settings in medicine, psychology, education, and economics (not to mention private-sector data analysis) where the database is not available to the analyst due to privacy concerns. A well-established solution is to allow the researcher to issue queries to the data which are proven to satisfy differential privacy. Differential privacy requires the addition of random noise to statistical queries and guarantees that the results reveal very little about any individual’s data.

In this paper we propose a new statistic for ANOVA, called $F_{1}$, that is specifically tailored to the differentially private setting. This statistic measures the same variations as the $F$ statistic, but uses $|a-b|$ instead of $(a-b)^{2}$ to measure the distance between $a$ and $b$. In the public setting the $F_{1}$ is a worse test statistic than the traditional $F$-statistic, but we show that in the private setting it has much higher power than the previously published differentially private $F$-statistic. That is, we show that it can detect effects with a little as 7% of the data that was previously required. (In one example, an effect that took 5300 data points to detect 90% of the time with $\epsilon=1$ in the prior work takes only 350 data points to detect using our new hypothesis test.)

1.1 Contributions and organization

We first review differential privacy, hypothesis testing, and the body of work that lies at the intersection of the two fields (Section 2). In Section 3 we then present a new test statistic, $F_{1}$, for ANOVA in the private setting. While there is some work on differentially private hypothesis testing, designing a new test statistic explicitly tailored for compatibility with differential privacy has been done by few others rogers2017new .

In Section 3.2 we give a private algorithm for computing an approximation of $F_{1}$ by applying the Laplacian mechanism to the computation of several intermediate values. Section 3.3 then describes how to compute the correct reference distribution for $F_{1}$ to in order to compute accurate $p$-values, which are the end result used by practitioners. Computing the reference distribution is complicated by the fact that, unlike the traditional $F$-statistic, $F_{1}$ is not scale-free.

We implement the private $F_{1}$-statistic and apply the method to different simulated datasets in Section 4. The computational experiments allow us to optimize $\rho$, a parameter that determines the allocation of our privacy budget between the two important intermediate values. We also compare our method to prior work campbell2018diffprivanova , and show an order of magnitude improvement in statistical power.

Finally, in Section 5 we present a generalization of the $F_{1}$ statistic that allows for an arbitrary exponent in the distance measure (besides the absolute value from $F_{1}$ and the $L^{2}$-norm from the traditional $F$). We find that the $L^{1}$-norm used in the $F_{1}$ statistic is the most powerful across a wide range of scenarios.

2 Background

We begin by discussing hypothesis testing in general, the one-way analysis of variance (ANOVA) test in particular, and differential privacy. Readers familiar with one or more of these topics should be able to skip the relevant sections. We then discuss how these topics come together in private hypothesis testing and related work in this area.

2.1 Hypothesis Testing

Hypothesis tests are common tools for making statistical inferences from data. The end goal of a hypothesis test is to determine whether a data set is consistent with a proposed model. This model is called the null hypothesis, denoted $H_{0}$, and it suggests a mechanism by which the data could have been generated. The mechanism is chosen to be scientifically meaningful, for example: the variable of interest has the same distribution across all of the treatment groups.

The comparison between $H_{0}$ and an observed data set is made using a test statistic. A test statistic $f$ is simply a function from the data set to the real numbers. The goal is to design a test statistic with a known distribution when the data comes from $H_{0}$, but which will follow a markedly different distribution under other scenarios. The question then becomes, for a given database $\mathbf{x}$ with $f(\mathbf{x})=t$, how likely is a value at least as extreme as $t$ to occur if $\mathbf{x}$ was drawn from $H_{0}$. To compute this probability, we need to compare $t$ to the reference distribution.

Definition 1** (Reference Distribution).**

Suppose $f$ is a function that computes a test statistic. The reference distribution for $f$ is the probability distribution of the statistic $T$ when $T=f(\mathbf{X})$ and $\mathbf{X}$ is drawn from a distribution consistent with $H_{0}$.

This reference distribution is used to calculate a $p$-value. A $p$-value is the probability, under the reference distribution, of drawing a statistic at least as extreme as the observed statistic.

Definition 2** ($p$-value).**

For a given test statistic $t=f(\mathbf{x})$ and null hypothesis $H_{0}$, the $p$-value is defined as

[TABLE]

The $p$-value provides context for the observed statistic by positioning it in the range of statistics that could be observed under $H_{0}$.

Typically, researchers choose a significance threshold $\alpha$ and reject the null hypothesis when their calculated $p$-value is less than $\alpha$. The $\alpha$-level determines the probability of a type I error, which occurs when an analyst rejects the null hypothesis despite it being true. The value of the statistic that demarcates this rejection region is called the critical value, denoted by $t^{*}$. That is, $\Pr[T\geq t^{*}\mid T=f(\mathbf{X})\text{ and }\mathbf{X}\leftarrow H_{0}]=\alpha$.

When one develops a test statistic, a primary goal is to maximize statistical power. The power of a test quantifies how effectively it can detect a deviation from $H_{0}$. It is the probability of rejecting when $H_{0}$ is false. Generally, the power depends on both the amount of data and the effect (i.e., how different the true distribution $H_{A}$ is from $H_{0}$).

Definition 3** (Statistical Power).**

For a specific alternate hypothesis $H_{A}$, the statistical power of a hypothesis test is

[TABLE]

2.2 One-way ANOVA

Consider again our example from earlier, where a researcher has, for a set of individuals, both blood pressure measurements and the mutation status of a particular gene. This is a classic setting for a one-way analysis of variance (ANOVA) test.

Index each individual person or observation in a database $\mathbf{x}$ with $i\in\{1,\ldots,N\}$. Each $i$ is associated with a group or category $c_{i}$ (e.g., mutation status) and a numerical value $y_{i}\in\mathbb{R}$ (e.g., blood pressure). We use $\bar{y}$ to represent the mean of all $N$ numerical values. Index each group with $j\in\{1,\ldots,k\}$. Each $j$ is associated with $C_{j}=\{i\mid c_{i}=j\}$, the set of indices of observations in group $j$. Denote the size of the set $n_{j}$ and the mean of the values in that set $\bar{y}_{j}$. (Note that this means $\bar{y}_{c_{i}}$ is the mean of numerical values in the same group as observation $i$.)

The null hypothesis in a one-way ANOVA is that the $y_{i}$ follow identical normal distributions regardless of their group. This motivates the test statistic used in an ANOVA test, the $F$-statistic, which measures the ratio of the variation between group means (weighted by the size of the groups) and the variation of individuals within each group. If all groups had equal means, the variation between group means would be proportional to the variation between individual observations.111Our approach relies upon simulating normally distributed data in correspondence with the traditional normality assumption. The one-way ANOVA test is known to be robust to deviations from normality schmider2010 , so our approach should be applicable even in settings where the normality assumption is suspect. Readers interested in more about ANOVA generally are referred to cox1974theoretical .

Definition 4** ($F$-Statistic).**

Given a database $\mathbf{x}$ with $k$ groups and $N$ total entries, the $F$-statistic is the ratio of two values, traditionally called $\textit{SSA}(\mathbf{x})$ and $\textit{SSE}(\mathbf{x})$. The Sum of Squared errors of All category means (SSA) is a measure of variance between group means, weighted by group size:

[TABLE]

The Sum of Squared Errors of all observations (SSE) is a measure of variance within groups:

[TABLE]

The $F$-statistic is the ratio of SSA and SSE, each divided by their respective degrees of freedom. These adjusted values are called the Mean Sum of All category errors (MSA) and Mean of Sum of Squared Errors (MSE). We can now finish defining the $F$-statistic.

[TABLE]

Under the null hypothesis that the $y_{i}$ follow identical normal distributions regardless of their group, the reference distribution of the $F$-statistic is known exactly. This comes from recognizing that $\textit{SSA}(\mathbf{x})$ is drawn from $\sigma^{2}\chi_{k-1}^{2}$, the chi-squared distribution with $k-1$ degrees of freedom scaled by within-group variance, and $\textit{SSE}(\mathbf{x})$ is drawn from $\sigma^{2}\chi^{2}_{N-k}$. The ratio of these values, therefore, has a reference distribution that is scale-free (not dependent on $\sigma$) and can be calculated knowing only $N$ and $k$.

2.3 Differential Privacy

Differential privacy is a security definition for the release of information about a database of private records. Here we outline the foundational definitions and theorems in differential privacy; everything below first appeared in the seminal paper of Dwork et al. dwork2006calibrating .

Suppose we have a database $\mathbf{x}$ containing sensitive information that we want to study. In particular, we want to publish the output $f(\mathbf{x})$ of a function $f$ (also sometimes called a mechanism) on our database while protecting the privacy of the individuals whose data was collected. Differential privacy promises that an adversary will learn approximately nothing about an individual as a result of their presence in $\mathbf{x}$. Informally, this is done by requiring that the probability of seeing any particular output is roughly the same regardless of what information a given individual submitted to the database.

As above, we use $N$ to represent the number of rows in $\mathbf{x}$, where a “row” is simply the set of data associated with a single individual. We now define neighboring databases, which differ in only one row.

Definition 5** (Neighboring Databases).**

Two databases $\mathbf{x}$ and $\mathbf{x}^{\prime}$ are neighboring if $\mathbf{x}$ can be transformed to $\mathbf{x}^{\prime}$ by changing only one individual’s data (where a change is an in-place modification, not a full addition or removal).

To protect the privacy of individuals in a database, differentially private requires that the output of a query on any two neighboring databases should look nearly identical.

Definition 6** (Differential Privacy).**

A (randomized) mechanism $f$ with range R is $\varepsilon$-differentially private if for all $S\subseteq R$ and for all neighboring databases $\mathbf{x}$ and $\mathbf{x}^{\prime}$

[TABLE]

The parameter $\varepsilon$ is called the privacy parameter, and its choice is a policy decision. The lower the chosen $\varepsilon$, the stronger the privacy guarantee. Note that because neighboring databases are the same size, $N$ can always be released without compromising privacy.222Differential privacy can also be defined in terms of databases that differ by an addition/deletion, rather than by a change in a row. For most applications these definitions are equivalent except for a change in $\epsilon$ by a factor of two (with the version here being the more stringent interpretation of $\epsilon$). The one significant difference is that the other definition does not result in $N$ being public, which is important for our work here.

Differential privacy has several useful properties. One of the most useful is composition:

Theorem 1** (Composition).**

Suppose $f$ and $g$ are respectively $\varepsilon_{1}$- and $\varepsilon_{2}$-differentially private mechanisms. Then, a mechanism $h$ that returns the results of applying $f$ and $g$ to $\mathbf{x}$, $h(\mathbf{x})=(f(\mathbf{x}),g(\mathbf{x}))$, is $(\varepsilon_{1}+\varepsilon_{2})$-differentially private.

In other words, the privacy guarantee decreases, but does not disappear, when a database is queried multiple times. Composition allows database administrators to issue researchers a privacy budget, which researchers can then divide up as they wish between any number of different queries.

Another defining feature of differential privacy is its resistance to post-processing.

Theorem 2** (Post-Processing).**

Let $f$ be an $\varepsilon$-differentially private mechanism, and let $g$ be an arbitrary function. Then, $h(\mathbf{x})=g(f(\mathbf{x}))$ is also $\varepsilon$-differentially private.

This theorem allows us to do any computation we desire on the output of our differentially-private mechanism without diminishing the privacy guarantees. We will utilize this property to compute $p$-values of the private $F$-statistic. The $p$-values will be automatically private without additional argument.

Our algorithms are constructed by taking building blocks and combining them with composition and post-processing, but the fundamental building blocks are made private using the Laplace mechanism, the oldest and maybe simplest method for achieving differential privacy. The Laplace mechanism allows the conversion of any function $f$ into a private approximation $\hat{f}$. One must first compute (or bound) the sensitivity of the function, the maximum effect on the output that a single row can have.

Definition 7** (Sensitivity).**

The sensitivity of a (deterministic) real-valued function $f$ on databases is the maximum of $\lvert f(\mathbf{x})-f(\mathbf{x}^{\prime})\rvert$ taken over all pairs $(\mathbf{x},\mathbf{x}^{\prime})$ of neighboring databases.333Sensitivity and the Laplace mechanism can be defined on functions with output in $\mathbb{R}^{n}$, but we only need the one-dimensional version.

The Laplace mechanism will use random noise drawn from the Laplace distribution.

Definition 8** (Laplace Distribution).**

The Laplace Distribution (centered at 0) with scale $b$ is the distribution with probability density function

[TABLE]

We use ${\sf Lap}(b)$ to represent the Laplace distribution with scale $b$.

We can now present the Laplace mechanism.

Theorem 3** (Laplacian Mechanism).**

Let $f$ be a function with sensitivity bounded above by $s$. Let $L$ be a random variable drawn from ${\sf Lap}(s/\varepsilon)$. Then the function $\hat{f}(\mathbf{x})=f(\mathbf{x})+L$ is $\varepsilon$-differentially private.

2.4 Differentially Private Hypothesis Testing

In order to create a differentially private hypothesis test, we need a private function $f$ of a database to serve as our test statistic. This could be a differentially private estimate of an existing test statistic, or it could be a new test statistic altogether. Because randomization is essential to differential privacy, $f$ will be randomized. The same statistic on the same database may yield different outputs each time it is computed.

In addition to a test statistic $f$, we require a suitable reference distribution to calculate the corresponding $p$-value. While it may be tempting to compute the $p$-value using the reference distribution for the non-private statistic one is estimating, this may yield wildly inaccurate results campbell2018diffprivanova , because adding noise to the statistic increases the probability of outlier output values. Instead, we must compute the reference distribution for the noisy statistic. Only then can we calculate an accurate $p$-value.

The goal of differentially private hypothesis testing is to create a private test statistic and method of computing the $p$-value that maximizes statistical power, ideally approaching the power of the equivalent test in the classical non-private setting.

2.5 Related Work

There has been a moderate amount of work on differentially private hypothesis testing, but because there are many hypothesis tests most individual tests have received only a small amount of attention, and some very common tests have not seen a private analogue developed at all.

Several papers have addressed testing the value of a mean or the difference of means solea2014differentially ; d2015differential ; ding2018comparing . Hypothesis tests using coefficients of a linear regression to test for dependence between continuous variables is extremely common in many academic disciplines, but only recently has a method for carrying this analysis out privately been developed sheffet2015differentially ; barrientos2017differentially . and Nguyên and Hui propose a test for surival analysis data nguyen2017differentially . There is one prior work on private ANOVA testing, that of Campbell et al. campbell2018diffprivanova . We will discuss this result in greater depth in the next section.

The chi-squared test, which tests for the independence of two categorical variables, has received the most study. Vu and Slavković vu2009differential give an analogue to the test and also compute accurate $p$-values. Many private chi-squared tests have been specifically motivated by genome-wide association studies (GWAS) fienberg2011privacy ; uhlerop2013privacy ; johnson2013privacy . These give p-value calculations, but they are only accurate in the limit as $N$ grows large. Other work has used Monte Carlo simulations (as we do in this work) to give more accurate reference distributions for small $N$ gaboardi2016differentially ; wang2015revisiting . Rogers and Kifer rogers2017new instead propose a new statistic with an asymptotic distribution more similar to its non-private analogue. We note that this is one of few papers that, like the present work, proposes test statistics intended for the private setting, rather than simply approximating the accepted test statistic from the classic public setting. Very few of these papers carefully measure the power of the test they develop. Rogers and Kifer rogers2017new and Gaboardi et al. gaboardi2016differentially are notable exceptions, giving power curves for several different approaches.

There is also a significant body of work looking at how quickly private approximations of test statistics converge to their limiting distributions (e.g., smith2008efficient ; wasserman2010statistical ; smith2011privacy ). These are important theoretical results, but they do not usually yield practical tests. Unless $N$ is very large (in which case the details of the test do not matter very much anyway) the distribution of the test statistic is not close enough to that of the standard public to allow accurate computation of $p$-values.

2.6 Prior work on private ANOVA

The only previous work on differentially private ANOVA testing that the authors are aware of is Campbell et al. campbell2018diffprivanova Using the ANOVA test as defined above, they analyze the sensitivity of the SSA and SSE with the assumption that all data was normalized to be between [math] and $1$ and add Laplacian noise proportional to these sensitivities to the public computation of the SSA and SSE. Their algorithm then uses post-processing to calculate the noisy $F$-statistic, and returns this in addition to the noisy SSA and SSE (Algorithm 1).

Normally, the $F$ distribution is used to calculate a $p$-value for the $F$-statistic. However, Campbell et al. find that the distribution of the private estimate $\widehat{F}$ differs too much from the $F$ distribution for this to be acceptable. Furthermore, they find that it is no longer scale-free, meaning that the distribution depends on the within-group variance $\sigma^{2}$.

Fortunately, the SSE is an estimate of $\sigma$, so using this estimate they computed an estimated distribution on $\widehat{F}$ through simulation. They could then compare a given value of $\widehat{F}$ to this distribution to obtain a $p$-value.

To assess the power of their private ANOVA algorithm, they simulate databases with three equal-sized groups with values drawn from $\mathcal{N}(0.35,0.15),$ $\mathcal{N}(0.5,0.15)$, and $\mathcal{N}(0.65,0.15)$ respectively. For several $(N,\varepsilon)$-pair choices, they generate many sets of data, apply the private ANOVA test, calculate the $p$-value, and record the percentage of simulations with $p$-values less than 0.05. They find that when $\varepsilon=1$, they need over five thousand data points to detect this effect (compared to two or three dozen data points in the public setting). Our goal in this paper is to reduce the gap between the public and private setting.

3 A New Test Statistic

In this section we describe our hypothesis test. This begins with the introduction in Section 3.1 of $F_{1}$, a new test statistic for the ANOVA setting. In Section 3.2 we then show how to privately compute a private approximation of $F_{1}$. Finally, in Section 3.3 we calculate $p$-values for the private $F_{1}$ statistic. This means simulating a correct reference distribution against which we can compare our output.

3.1 The $F_{1}$ Statistic

Our goal is to define a statistic that releases similar information as the $F$-statistic, but has higher power (Definition 3) for reasonable privacy guarantees. We focused on two approaches to improve the power of the ANOVA calculation: reducing the amount of Laplacian noise by decreasing sensitivity, and making SSA and SSE numerically larger so that the noise has less influence over the total value of the statistic. We achieved both goals by taking the absolute values of the summand terms in the SSA and SSE, rather than squaring them.

As before, let $k$ denote the number of categories in the database, $c_{i}$ be the category and $y_{i}$ be the numerical value associated with observation $i$, $n_{j}$ be the size of category $j$, and $N$ be the total size of the data set. Additionally, $\bar{y}$ is the grand mean of the entries in database $\mathbf{x}$ and $\bar{y}_{j}$ is the mean of the entries in group $j$.

Definition 9** ($F_{1}$-statistic).**

Given a database $\mathbf{x}$ with $k$ groups and $N$ total entries, the $\textit{SA}(\mathbf{x})$ and $\textit{SE}(\mathbf{x})$ calculations are defined as follows:

[TABLE]

The $F_{1}$-statistic is the ratio of SA and SE, each divided by their respective degrees of freedom.

[TABLE]

The $F_{1}$-statistic measures variation between group means compared to variation within groups in essentially the same way as the original $F$-statistic. The SA grows as the group means diverge. What constitutes a “large” variation between group means depends on the variation between individual items, so SE, which measures this individual-level variation, provides a sense of scale for the SA value.

In the next section, we show that the sensitivities of SA and SE in the $F_{1}$-statistic are less than half as large as the sensitivities of SSA and SSE in the original $F$-statistic. Further, because the summand terms are restricted to $[0,1]$, the SA and SE values are larger, meaning they can tolerate the addition of more noise before losing their usefulness.

3.2 A Private Approximation of $F_{1}$

The sensitivity of $F_{1}$ is very high. (In the worst case, $\textit{SE}(\mathbf{x})$ is almost zero and very small changes can have huge effects on $F_{1}(\mathbf{x})$.) As a result, we can’t simply apply the Laplace mechanism to $F_{1}$. Instead, we choose to apply it individually to the SA and SE functions, and then use composition and post-processing to compute an estimate of $F_{1}$. We must therefore bound the sensitivities of SA and SE.

We assume that the number of valid category values, $k$, is fixed and public, but the number of entries in each group is not. (This includes the possibility that one or more categories exist as valid entries but do not appear in the actual database.) We also assume that there are maximum and minimum possible values for the data, and that the computation first uses these to normalize the data, mapping it to the interval $[0,1]$.

Theorem 4** (SE Sensitivity).**

The sensitivity of the SE calculation in Definition 9 is bounded above by 3.

Proof.

Suppose neighboring databases $\mathbf{x}$ and $\mathbf{x^{\prime}}$ differ by some row $r$. Say that in $\mathbf{x}$, $c_{r}=a$, and in $\mathbf{x^{\prime}}$, $c_{r}=b$. Let $n_{a}$ be the size of category $a$ excluding $r$ and let $n_{b}$ be the size of category $b$ excluding $r$. We begin by expressing the SE calculation as nested summations indexing over group size and entries within each group.

[TABLE]

Call $t_{i}=\lvert y_{i}-\bar{y}_{c_{i}}\rvert$, and let $c_{i}\neq a,b$. Then, $t_{i}$ does not change between $\mathbf{x}$ and $\mathbf{x^{\prime}}$. Now, suppose $i\neq r$ but $c_{i}=a$. It follows that $\Delta t_{i}\leq 1/(n_{a}+1)$, since the only change comes from $\bar{y}_{a}$. There are $n_{a}$ such terms, so the total contribution from these terms is at most $n_{a}/(n_{a}+1)$. If $a\neq b$, we must also consider $t_{i}$ where $i\neq r$ but $c_{i}=b$, for which we have $\Delta t_{i}\leq 1/(n_{b}+1)$. Thus, the terms in groups $a$ and $b$ excluding row $r$ together contribute

[TABLE]

if $a\neq b$ and just $n_{a}/(n_{a}+1)<1$ otherwise.

Now, consider $\Delta t_{r}$. Since $y_{r}$, $\bar{y}_{a}$, and $\bar{y}_{b}$ are all in the interval $[0,1]$, the difference between $t_{r}$ in database $\mathbf{x}$ and in $\mathbf{x^{\prime}}$ is at most 1. Thus, the total sensitivity of SE is bounded above by 3.

∎

Theorem 5** (SA Sensitivity).**

The sensitivity of the SA calculation in Definition 9 is bounded above by 4.

Proof.

Again, suppose neighboring databases $\mathbf{x}$ and $\mathbf{x^{\prime}}$ differ by some row $r$ in both the categorical and numerical values, with $c_{r}=a$ in $\mathbf{x}$, and $c_{r}=b$ in $\mathbf{x^{\prime}}$. Denote the number of entries in groups $a$ and $b$ not including row $r$ by $n_{a}$ and $n_{b}$, respectively. We begin by expressing the SA calculation as two sums indexing over groups and entries within each group.

[TABLE]

Consider the change from $\textit{SA}(\mathbf{x})$ to $\textit{SA}(\mathbf{x^{\prime}})$ as if it occurred in two steps. In the first step, the grand mean $\overline{y}$ is updated. Note that in the worst case, $\bar{y}$ can change by at most $1/N$ between $\mathbf{x}$ and $\mathbf{x^{\prime}}$. Then, since there are $N$ summands, each including the grand mean, this step changes the value by a maximum of $N(1/N)=1$ to the overall sensitivity.

In the second step we change the group means for groups $a$ and $b$. There will be $n_{a}$ terms containing $\overline{y}_{a}$, each of which will change by at most $1/n_{a}$, changing the overall value by at most 1. Similarly updating $\overline{y}_{b}$ changes $n_{b}$ terms each by at most $1/n_{b}$ for a total contribution of 1.

Finally, the change of $r$ between $\mathbf{x}$ and $\mathbf{x^{\prime}}$ contributes 1 to the overall sensitivity, and thus the sensitivity of the SA is bounded above by 4.

∎

Having proven these sensitivities, we can now introduce our private algorithm (Algorithm 2) for approximating $F_{1}$. The algorithm first estimates SA and SE, allocating part of the privacy budget to each one. We introduce a parameter $\rho\in(0,1)$ that determines the relative amount of the privacy budget spent on each intermediate value. The optimal value of $\rho$ will be experimentally determined in Section 4.2. We note that in addition to using a different test statistic, the prior work by Campbell et al. did not consider $\rho$ values other than 0.5.

Theorem 6**.**

Algorithm 2 is $\varepsilon$-differentially private.

Proof.

By the sensitivity bounds of SE and SA in Theorems 4 and 5 and the Laplace mechanism (Theorem 3), $\widehat{\textit{SA}}$ is $\rho\varepsilon$-differentially private and $\widehat{\textit{SE}}$ is $(1-\rho)\varepsilon$-differentially private. By the composition theorem (Theorem 1), outputting both is $\varepsilon$-differentially private. Since $k$ and $N$ are both public information, computing $\widehat{F_{1}}$ is post-processing (Theorem 2). ∎

3.3 Reference Distribution and $p$-values

As discussed previously, the test statistic on its own is not useful; we need a $p$-value to provide a sense of scale in the context of the null hypothesis. Computing a $p$-value begins with an accurate reference distribution. We numerically approximate this distribution through simulation. The intermediate values of the $F$ statistic, SSA and SSE, are drawn from $\sigma^{2}\chi_{k-1}^{2}$ (for SSA) and $\sigma^{2}\chi_{n-k}^{2}$ (for SSE). Campbell et al. campbell2018diffprivanova used this to easily sample from the correct distributions for SSA and SSE.

The SA and SE values needed to compute $F_{1}$ follow no similarly tractable distribution, so instead we simulate full databases according to the null hypothesis and for each calculate $\widehat{F_{1}}$. The distribution of these $\widehat{F_{1}}$ values approximates the reference distribution. The goal, given a database $\mathbf{x}$, is to simulate databases with the same size $N$ and number of groups $k$, same standard deviation $\sigma$, and same expected value $\mu$. The values $N$ and $k$ are public, so we can use those values. The expected value $\mu$ is not, but as long as it is safely inside the $[0,1]$ interval, its value has no effect on the distribution, so we simply always use 0.5.

Unfortunately, using an accurate $\sigma$ is more difficult. Unlike in the non-private setting, the reference distribution depends on the choice of $\sigma$, so an inaccurate value can cause incorrect results. We had two choices: either we could use some of our $\varepsilon$-budget to directly estimate the standard deviation of the $y_{i}$ in database $\mathbf{x}$, or we could devise an indirect method of estimating the standard deviation given $\widehat{\textit{SA}}$ and $\widehat{\textit{SE}}$. See Appendix D for further discussion of the first option. (It has higher power for low $n$ but takes longer to approach full power.) Here we focus on the second option by deriving an unbiased estimator $\hat{\sigma}$ for $\sigma$ that can be computed from $\widehat{\textit{SE}}$.

[TABLE]

See Appendix A for the proof that this is indeed an unbiased estimator.

Computing this estimator requires knowledge of each of the group sizes, which are private. Fortunately, $\tilde{N}$ is closely approximated by $N-k$, which we know. At the smallest database sizes that we considered (around $N=100$ with three equal-sized groups), this approximation has < 1% error. As the size of the database grows, this error shrinks to zero. Figure 1 visualizes the unbiasedness and asymptotically shrinking variance of $\hat{\sigma}$. (using the $\tilde{N}=N-k$ approximation). We further confirm in Section 4 that it is precise enough to compute valid $p$-values.

Another issue presented by the private estimation of $\sigma$ is that the Laplacian noise can be large enough to make the estimate negative. Negative standard deviation estimates are more likely to occur when the SE is small, i.e. when the database is small or when the within-group standard deviation is small.

This problem is unique to the private setting and we do the most conservative possible thing — we never reject the null hypothesis when the estimated standard deviation was negative. Our reasoning was that a negative standard deviation has no statistical meaning and any calculations made from such an estimate would be uninformative. This method for dealing with negative standard deviation estimates results in a type I error rate lower than $\alpha$, which means that our results are more conservative. As was mentioned earlier, the SA increases as the database size and effect size between groups increases. Thus, negative standard deviation estimates tend to occur for database sizes and effect sizes that are so small that the effect is undetectable anyway. As a result, this conservative choice does little to reduce the power of our test.

A formal description of the private $F_{1}$ ANOVA test in presented in Algorithm 3, which returns a Boolean indicating whether the null hypothesis $H_{0}$ is rejected.

Theorem 7**.**

Algorithm 3 is $\varepsilon$-differentially private.

Proof.

This follows immediately from Theorem 6 (privacy of the test statistic) and Theorem 2 (post-processing). Using the notation of Theorem 2, private_F1 (Algorithm 2) is the function $f$, which is proven $\epsilon$-differentially private by Theorem 6. The rest of the computation is the function $g$, and uses only the output of $f$. The whole of Algorithm 3 is therefore the composition of these two functions and as a result is itself private. ∎

4 Experimental Results

In this section we assess the properties and performance of Algorithm 3 through simulation.

4.1 Properties of the test

A $p$-value is considered valid if  casellaberger2002

[TABLE]

In other words, the actual type I error rate must be less than or equal to $\alpha$. This can be assessed by conducting many tests on simulated data sampled according to $H_{0}$ at various $\alpha$ and checking if the proportion of rejections is less than or equal to $\alpha$.

Figure 2 presents these simulations for several choices of $\epsilon$ and demonstrates that the $F_{1}$ test produces valid $p$-values. The type I error rate is lowest at low $\epsilon$ values because the high privacy guarantee requires a large amount of noise be added to SE, which can lead to a negative estimate for $\sigma$, which in turns leads to an automatic decision to retain $H_{0}$ (see section 3.3).

We note also that when we simulate data for calculating a reference distribution, we always simulate data with equal-size groups. I.e., we must confirm that the critical value of the reference distribution is highest when groups are of equal size. Fortunately, this appears to be the case. Appendix B contains both experimental and theoretical arguments for this claim, though not a complete analytic proof.

4.2 Optimal $\rho$

The computation of the private $F_{1}$ statistic requires the specification of $\rho\in(0,1)$, the parameter that determines the proportion of the privacy budget that is allocated to SA relative to SE. We determined an optimal value for $\rho$ by constructing power plots comparing database size to power for different $\rho$ values. We began by exploring the full range from $0.1$ to $0.9$ by $0.1$-increments to get a sense of the range of variability in power. After that initial pass, we tuned in to the value with higher precision. We considered many effect sizes and found that in every case $\rho\approx 0.7$ was the most powerful.

Figure 3 is an example of one of the many scenarios that were considered, which identifies .7 as the optimal value. The scale of the effect of $\rho$ on power is not dramatic but it was persistent across scenarios.

4.3 Performance of the test: power

There can be many tests for a given scientific setting that generate valid $p$-values and have identical type I error rates. What distinguishes them is their statistical power, or the probability that they reject $H_{0}$ when the database comes from a distribution under $H_{A}$.

The most common way to assess the power of a test is to generate a plot of power as a function either of database size or of effect size. Figure 4 fixes the effect size and shows power curves as a function of database size for four choices of $\epsilon$. An ideal test would very quickly develop power near 1 with very little data. In our private setting, it is clear that the cost of high privacy ($\epsilon=.1$) is roughly an order of magnitude more data than modest privacy ($\epsilon=1$) to achieve high power. Our private test (for reasonable values of $\epsilon$) still requires much more data than the public version, hundreds of data points as opposed to dozens.

Each point in a line of a power curve was computed from 10,000 simulations, each based on a synthetic database of a set sample size and a set effect size (distance between the group means). We considered the effect $H_{A}$ where $k=3$ groups, each distributed $\mathcal{N}(0.35,0.15),\mathcal{N}(0.5,0.15),\mathcal{N}(0.65,0.15)$ and equal group sizes, which was the same scenario as Campbell et al. campbell2018diffprivanova . Each simulation starts with a draw from $H_{A}$ and computes a $p$-value as described in Algorithm 3.

The most relevant comparison, however, is between the private $F_{1}$ statistic and the privatized version of the classical $F_{2}$ statistic proposed in Campbell et al., the only prior private version of ANOVA. This comparison can be seen in Figure 5, and the improvement is substantial. For example, at $\epsilon=1$, if one wanted to collect enough data to detect this effect with 80% probability, one would need  4500 observations if using the prior best test, whereas with the $F_{1}$ test presented here one would need only  300 data points, a 15-fold reduction in the necessary amount of data.

Figure 5 also demonstrates the degree to which the $F_{1}$ is well-suited to the private setting; the greatest improvement in power occurs under high privacy ($\epsilon=0.1$). As $\epsilon$ grows large, the difference between the two statistics shrinks. In the public setting, they are nearly indistinguishable, though the $F_{2}$ is narrowly more powerful at every database size.444This empirical result is consistent with theoretical results in the classical statistics literature that discuss conditions in which the $F_{2}$ is a most powerful test cox1974theoretical .

The improvement in power is the greatest practical contribution of our work: the ability to conduct a private ANOVA with an order of magnitude less data than the existing approach. This improvement can be attributed to two key characteristics of our new test. The first and most important is the notion of measuring distance using the $L^{1}$ norm. The second is the unequal apportionment of the privacy between the SA and SE terms.

5 Other considerations

In this section we present a straightforward generalization of the classic $F$ and $F_{1}$ to allow other exponents. In Section 5.1 we show that in this generalized class of statistics $F_{1}$ is indeed optimal. In Section 5.2 we discuss the experimental exploration of the full parameter space that guides the formulation of the $F_{1}$ statistic.

5.1 Varying the exponent

As seen is the previous sections, the change from squaring the differences in the original $F$-test to taking their absolute value with the $F_{1}$-statistic improved power significantly in the differentially private setting. It is not obvious that switching the exponent from 2 to 1 is optimal — perhaps some other exponent is superior. An exponent of 0 is clearly horrible, so there must be a local maximum in the power of the statistic for some exponent between 0 and 2.

In order to determine which exponent is in fact optimal, we further generalize the notion of an $F$-test. We define SQA and SQE, which are equivalent to SSA and SSE except that the summand is raised to the $q^{\text{th}}$ exponent, and we call the resulting statistic $F_{q}$. Note that $F_{1}$ (as defined earlier) is a special case of $F_{q}$ for $q=1$, and $F_{2}$ is the standard $F$-statistic.

Definition 10** ($F_{q}$).**

Given a database $\mathbf{x}$ with $k$ groups and $N$ total entries, define SQA and SQE as follows:

[TABLE]

[TABLE]

Then, $F_{q}$ is defined as

[TABLE]

We must now create a private approximation of $F_{q}$ for arbitrary $q$. To do this, we first bound the sensitivity of the SQA and SQE with the following two theorems, the proofs of which can be found in Appendix C.

Theorem 8** (SQE Sensitivity).**

The sensitivity of SQE is bounded above by

[TABLE]

when $q\in(0,1)$ and

[TABLE]

when $q\geq 1$. Note that both give an upper bound of 3 when $q=1$.

Theorem 9** (SQA Sensitivity).**

The sensitivity of SQA is bounded above by

[TABLE]

when $q\in(0,1)$ and

[TABLE]

when $q\geq 1$. Note that both give an upper bound of 4 when $q=1$.

Given these sensitivity bounds, we can calculate a private approximation of $F_{q}$ for simulated data using the same algorithm as for $F_{1}$, but with the sensitivities altered according to the choice of $q$. We can also simulate a reference distribution by adapting Algorithm 3, which was done to construct Figure 6. As is clear from these results, in terms of power, the optimal value of $q$ is 1.

We note that in the computation shown in Figure 6, we don’t estimate $\sigma$ using SQE. This is because we have not developed an estimator for $\sigma$ that can be computed from SQE (for $q\neq 1,2$). If another value of $q$ had indeed been optimal, the next step would have been to find such an estimator and confirm it was accurate enough to produce acceptable $p$-values. But since the power cannot possibly improve when switching to an estimated $\sigma$ value, this result is sufficient to show that other values of $q$ need not be considered.

5.2 Parameter tuning

With the generalization of the $F$-statistic, we add $q$ to the list of parameters that determine the power of a testing procedure. The parameters can be organized as follows:

Data Generation: $N$, $k$, $\sigma$, effect size

Private Algorithm: $\epsilon$, $q$, $\rho$

The analyst gets to select the parameters corresponding to the private algorithm. While $\epsilon$ is set based on privacy concerns, $q$ and $\rho$ should be set to maximize power, which our work suggest occurs at roughly 1 and 0.7, respectively. This conclusions is based upon an extensive exploration of the parameter space, a selection of which can be found in Appendix E.

The salient feature of these plots is that the choice of $q$ is much more consequential than the choice of $\rho$. In the setting where $\epsilon=.1$, we found that the result seen in Figure 5 – a greater than 10-fold reduction in database size to get equivalent power – holds across a range of difference data generation parameters.

By contrast, the effect of $\rho$ on power is much smaller; the database size reduction is closer to 1.1- or 1.2-fold when moving from $\rho=.5$ to $\rho=.7$ when $\epsilon=.1$.

Acknowledgements

This material is based upon work supported by the National Science Foundation under Grant No. SaTC-1817245 and by the Gillespie Family Student Research Fund.

Appendix A Appendix: An unbiased estimator for $\sigma$

We are concerned with finding an unbiased estimator for $\sigma$ given $\widehat{SE}$.

[TABLE]

Our first goal is to find $E(\widehat{SE})$. We begin by noting several relevant distributions.

[TABLE]

The term $\tau^{2}$ will be useful in determing the expected value of each term in $\widehat{SE}$, so we seek to express it in terms of known quantities.

[TABLE]

To solve for the covariance, noting that $Cov(X,Y)=E(XY)-\mu^{2}$, we start by finding

[TABLE]

where $S_{-i}=y_{1}+y_{2}+\ldots+y_{i-1}+y_{i+1}+\ldots+y_{n_{c_{i}}}$. The term $y_{i}^{2}$ will be chi-square distributed when standardized as follows:

[TABLE]

Therefore we can write

[TABLE]

Continuing from Eq. (3),

[TABLE]

Now we can return to the covariance term in Eq. (2),

[TABLE]

Now we can finish the calculation of $\tau^{2}$.

[TABLE]

The distribution of the absolute value of a normal random variable with mean 0 and variance $\tau^{2}$ is half normal with a single parameter $\tau$, which should properly be indexed by the observation. For ease of notation, call this random variable $W_{i}$.

[TABLE]

where $E(W_{i})=\tau_{i}\sqrt{\frac{2}{\pi}}$. To find $E(\widehat{SE})$, we can write it using the double sum notation as in Eq. (1). Since $L_{2}$ is a Laplace distribution centered at zero, we can use linearity of expectation to further simplify the expectation:

[TABLE]

The change in indices is justified by realizing that $W_{i}$ is the same for all $n_{c_{i}}$ elements in $C_{i}$ (in expectation). Continuing,

[TABLE]

Denote the sum $\tilde{N}$, which is a number that approaches $N$ as the group sizes get large. This allows us to express the expected value more concisely as

[TABLE]

The final step is to correct for this bias in our final estimator:

[TABLE]

Using this exact estimator requires knowledge of each of the group sizes, which are private. Instead of dedicating part of the $\epsilon$ budget to this estimation, we used $N-k$ in place of $\tilde{N}$. At the smallest database sizes that we considered (around $N=100$), this approximation accounts for $<1\%$ error. As the size of the database grows, this error shrinks to zero.

Appendix B Appendix: Validity under unequal $n_{j}$

Consider a specific setting in which $N=800$, $\sigma=0.1$, and $k=4$. Figure 7 shows the reference distributions of $F_{1}$ in four scenarios, each with a different group allocation of the 800 observations. The distribution in $s_{0}$, the equal group size scenario, generally takes the highest values; indeed it exceeds the other scenarios at every quantile. This represents the distribution that we use to calculate $p$-values and reject $H_{0}$ whenever an observed statistic is greater than the vertical dotted line (when $\alpha=.05$). This means that in the other three scenarios, featuring unequal allocation, the actual type I error rate (the proportion of the distributions beyond the dotted line) will be less than $\alpha$ and we meet the condition for valid $p$-values.

In the following, we show that under any scenario, the expected value of the $F_{1}$ statistic is maximized when the group sizes are equal. Recall that the form of the statistic is

[TABLE]

where

[TABLE]

and

[TABLE]

Only $SA$ is a function of the $n_{j}$, so we restrict our attention to that term as $E(F_{1})$ for different allocations will scale with $E(SA)$ by a multiplicative constant.

Denote each term in the summation $SA_{j}$. By using the approach used in Appendix A to find the distribution of $y_{i}-y_{c_{i}}$, it can be shown that $\bar{y}_{j}-\bar{y}$ is distributed normal mean zero, and variance $\omega_{j}^{2}=\frac{\sigma^{2}}{n_{j}}-\frac{\sigma^{2}}{N}$.

[TABLE]

Since each $|\frac{(\bar{y}_{j}-\bar{y})}{\omega_{j}}|$ is distributed standard half normal, the expectation of each evaluates to the constant $\sqrt{\frac{2}{\pi}}$. Now we may also factor out $\sigma$ from $\omega_{j}$, providing the final answer:

[TABLE]

Proposition 1**.**

The expectation of $SA$ is maximized when group sizes are all equal, i.e., when $n_{j}=n_{j^{\prime}}$, for all $j\neq j^{\prime}$

Proof.

It suffices to show that the sum between any two summands of the expectation is maximized when the other $k-2$ group sizes are fixed, so we begin by fixing each $n_{j}$ except for $n_{1}$ and $n_{2}$, without loss of generality. Define $c=n_{1}+n_{2}$ to be the sum between the two unfixed group sizes, so then the sum between the first two summands becomes:

[TABLE]

since $n_{1}$ and $c-n_{1}>0$. Now factoring constants and differentiating with respect to $n_{1}$:

[TABLE]

which has only one solution in $n_{1}$, when $n_{1}=\frac{c}{2}=\frac{n_{1}+n_{2}}{2}$, which implies $n_{1}=n_{2}$. This critical point is associated with a maximum expectation on the first two summands, as desired. Applying this result across all $n_{j},n_{j^{\prime}}$ pairs such that $j\neq j^{\prime}$, will result in maximizing each summand of the expectation when all $n_{j}=n_{j^{\prime}}$, hence maximizing $\text{E}(SA)$. ∎

Appendix C Appendix: Sensitivity Proofs for SQA and SQE

Recall that, from Definition 10, we have the following functions SQA and SQE:

[TABLE]

[TABLE]

Our goal here is to bound the sensitivities for SQA and SQE.

Theorem 7** (SQE Sensitivity).**

The sensitivity of SQE is bounded above by

[TABLE]

when $q\in(0,1)$ and

[TABLE]

when $q\geq 1$. Note that both give an upper bound of 3 when $q=1$.

Proof.

As in the previous proofs of the sensitivity of SE and SA, suppose neighboring databases $\mathbf{x}$ and $\mathbf{x^{\prime}}$ differ by some row $r$, with $c_{r}=a$ in $\mathbf{x}$ and $c_{r}=b$ in $\mathbf{x^{\prime}}$. Rewrite the SQE as a sum that indexes over group size and entries within each group.

[TABLE]

Let $t_{i}=\left|y_{i}-\overline{y}_{c_{i}}\right|^{q}$ for any entry $i$. Note that if $c_{i}\neq a,b$, then $t_{i}$ will not change between databases $\mathbf{x}$ and $\mathbf{x^{\prime}}$, as the group means of the other groups are not altered. Thus, unless $i$ is in group $a$ or $b$, $t_{i}$ will contribute nothing to the overall sensitivity of the SQE. For notational ease, let $z=y_{i}-\overline{y}_{c_{i}}$ and parameterize $t_{i}$ as a function of $z$. The maximum $z$ can change by is $1/n_{c_{i}}$. We now bound the sensitivity by individually bounding the sensitivity $\Delta t_{i}$ of each term with $c_{i}=a,b$.

Case 1: When $q$ is less than $1$, $t_{i}(z)$ is a concave function with minimum at $z=0$ that is symmetric about the y-axis and which monotonically increases for positive $z$. Because the slope of $t_{i}(z)$ is highest near $z=0$, the worst case sensitivity is between $z=0$ and $z=1/n_{c_{i}}$, and hence

[TABLE]

Note that when $c_{i}=a$ and $i\neq r$, $\Delta t_{i}\leq(1/n_{a})^{q}$. The analogous statement holds for group $b$. By multiplying these bounds by the number of terms in each group, we get

[TABLE]

where the first and second terms are the total change possible to terms in groups $a$ and $b$ respectively, and the final term is for the contribution of row $r$ itself, which we cannot bound other than by noting that its value both in $\mathbf{x}$ and $\mathbf{x^{\prime}}$ falls inside of $[0,1]$.

We must now take the worst-case value of this bound over all possible database sizes $n_{a}$ and $n_{b}$. Since $n_{a}$ is a positive integer and $q<1$, $n_{a}^{1-q}$ increases for a given $q$ as $n_{a}$ increases (and the same is true for $n_{b}$). So, the worst-case sensitivity will occur when as much of the total database is in groups $a$ and $b$ as possible. Write $n_{b}=N-n_{a}$. Then

[TABLE]

is a downward facing parabola-like function with maximum value when $n_{a}=N/2$. So, the sensitivity of the SQE is bounded above by

[TABLE]

Case 2: When $q$ is greater than $1$, $t_{i}$ is a convex function with maximum at $z=1$, symmetric about the $y$-axis, and monotonically increasing for positive $z$. So, the worst case sensitivity is between $z=1$ and $z=1-1/n_{c_{i}}$, and hence

[TABLE]

Summing these bounds over all terms, we have

[TABLE]

Again, the sensitivity will be maximized when as much of the database is distributed between $n_{a}$ and $n_{b}$ as possible. To determine what the worst case allocation is, let

[TABLE]

i.e., $f$ is an expression for the upper bound of $\Delta\textit{SQE}$ with $n_{b}$ replaced by $N-n_{a}$ to maximize sensitivity. Then, we can maximize $f$ in terms of $n_{a}$:

[TABLE]

has a critical point at $n_{a}=N/2$, and

[TABLE]

is always negative. So, $f$ is concave down and $n_{a}=N/2$ is a global maximum. Hence, the worst case sensitivity occurs when the database is distributed equally between groups $a$ and $b$, i.e.,

[TABLE]

∎

Theorem 8** (SQA Sensitivity).**

The sensitivity of SQA is bounded above by

[TABLE]

when $q\in(0,1)$ and

[TABLE]

when $q\geq 1$. Note that both give an upper bound of 4 when $q=1$.

Proof.

Let $s_{j}=\left|\bar{y}_{j}-\bar{y}\right|^{q}$ for any group $j$. I.e., $s_{j}$ is the (unweighted) term in the calculation of the SQA that corresponds to group $j$. Note that as the grand mean changes between databases $\mathbf{x}$ and $\mathbf{x^{\prime}}$ in addition to the group means, all terms, not just those for groups $a$ and $b$, will contribute to the sensitivity of the SQA. Recall that the sensitivity of the grand mean is $1/N$, while the sensitivity of the group mean for groups $a$ and $b$ are $1/n_{a}$ and $1/n_{b}$ respectively.

Case 1: When $q<1$, $s_{j}(z)$ is a concave function with minimum at $z=0$, symmetric about the $y$-axis, and monotonically increasing for positive $x$. So, the worst case sensitivity of $\Delta s_{j}$ for $j\neq a,b$ is between $x=0$ and $x=1/N$, and the worst case sensitivity of $\Delta s_{j}$ for $j=a,b$ is between $z=0$ and $z=1/n_{j}+1/N$. Then, the total sensitivity of the SQA for $q<1$ is

[TABLE]

The addition of the $1$ comes from the fact that our data point $r$ that switches between groups contributes $|\bar{y}_{a}-\bar{y}|^{q}$ to the calculation of the SQA in database $\mathbf{x}$, and contributes $|\bar{y}_{b}-\bar{y}|^{q}$ to the calculation of the SQA in database $\mathbf{x^{\prime}}$; the difference between these two terms is bounded above by 1. Note that since $q<1$, $(1/z)^{q}>1/z$. Hence, $(1/N+1/n_{a})^{q}>(1/N)^{q}$ and thus the worst-case sensitivity occurs when all of $N$ is allocated to groups $a$ and $b$. I.e.,

[TABLE]

Then, as in the proof of the SQE’s sensitivity, to determine the worst-case sensitivity in terms of $N$, let

[TABLE]

and maximize this expression in terms of $n_{a}$.

[TABLE]

This is a symmetric expression between $N$ and $N-n_{a}$. So, there must be a critical point at $n_{a}=N/2$. Note also that

[TABLE]

since $N>n_{a}$ and $q<1$. So, $n_{a}=N/2$ is a global maximum, and hence

[TABLE]

Case 2: When $q$ is greater than $1$, $s_{j}$ is a convex function with minimum at $z=0$, symmetric about $z=0$, and monotonically increasing for positive $z$. So, the worst case sensitivity of $\Delta s_{j}$ for $j\neq a,b$ is between $z=1$ and $z=1-1/N$, and the worst-case sensitivity of $\Delta s_{j}$ for $i=a,b$ is between $z=1$ and $z=1-1/N-1/n_{i}$. Then, the total sensitivity of the SPA for $q<1$ is

[TABLE]

Note that

[TABLE]

Thus, $\Delta\textit{SQA}$ is maximized when $N$ is maximally allocated to groups $a$ and $b$. As in the proof for $q<1$, this occurs when $n_{a}=N/2=n_{b}$. Then,

[TABLE]

∎

Appendix D Appendix: Direct Calculation of $\sigma$

In our work, we used $\widehat{SE}$ to form an estimator for $\sigma$ when calculating the null distribution. We developed another version of a differentially-private ANOVA that calculates $\sigma$ directly using a portion of the privacy budget. We first define a generalized variance and prove the sensitivity bounds on this quantity.

Definition 11** ($\textit{VAR}_{q}$).**

Given a database $\mathbf{x}$ with $k$ groups and $n_{j}$ entries in the $j$-th group, the $\textit{VAR}_{q}$ calculation is defined as

[TABLE]

where $q$ is a positive real number.

Theorem 9** ($\textit{VAR}_{q}$-Sensitivity).**

The sensitivity of $\textit{VAR}_{q}$ is bounded above by

[TABLE]

when $q<1$, and is bounded above by

[TABLE]

when $q>1$. Note that these both give a bound of $2-1/N$ when $q=1$.

Proof.

Consider databases $\mathbf{x}$ and $\mathbf{x^{\prime}}$ which differ in entry $r$. Recall that the sensitivity of the grand mean is $1/N$. Let $t_{i}=\left|y_{i}-\bar{y}\right|^{q}$. When $q<1,t_{i}$ is a concave function with positive range, so the worst case sensitivity is between $y_{i}-\bar{y}=0$ and $y_{i}-\bar{y}=1/N$. That is,

[TABLE]

Every single term can be affected by at most $\left(\frac{1}{N}\right)^{q}$, except the term $r$, which can change by $1$. So,

[TABLE]

When $q>1,t_{i}$ is convex. The worst case sensitivity is between $y_{i}-\bar{y}=1$ and $y_{i}-\bar{y}=1-1/N$ Then,

[TABLE]

∎

We again use the Laplace mechanism, and the sensitivity of a database’s variance follows directly from Thm. 9:

Corollary 1**.**

The sensitivity of the variance of a database is bounded above by

[TABLE]

Algorithm 4 divides the privacy budget into $\rho_{1}$ for the SA, $\rho_{2}$ for the SE, and $\rho_{3}$ for the VAR calculations respectively. The values of $\rho$ are provided as input along with the database $\mathbf{x}$ and the $\varepsilon$ value. When $\widehat{SE}$ is used as the $\sigma$ estimate, we found a 70-30 split of the privacy budget between SA and SE was optimal (Figure 3). In the power analysis of Algorithm 4, we vary the proportion of the privacy budget dedicated to the VAR calculation. We fixed the proportion of the privacy budget for $\rho_{1}$ and $\rho_{2}$ to be a 70-30 split budget not used by $\rho_{3}$ (the VAR calcuation).

Unsurprisingly, allocating part of the budget to calculating the standard deviation negatively impacts the statistical power of the test (Fig. 8). However, it is surprising that this allocation does not impact the power of the test more strongly, since VAR has a larger sensitivity and requires adding noise to smaller values than the SA and SE.

Appendix E Appendix: Parameter Selection

Empirical results guiding parameter selection are presented on the following pages in Figure 9 ($\epsilon$ = .1) and Figure 10 ($\epsilon$ = 1).