Differentially private methods for managing model uncertainty in linear regression models

Journal of Machine Learning Research 25 (2024) 1-44 Submitted 12/21; Revised 2/24; Published 3/24 Differentially Private Methods for Managing Model Uncertainty in Linear Regression Models V´ıctor Pe˜na victor.pena.pizarr[email protected] Department d’Estad´ıstica i Investigaci´o Operativa Universitat Polit`ecnica de Catalunya Barcelona, Spain Andr´es F. Barrientos [email protected] Department of Statistics Florida State University Tallahassee, FL 32306, USA Editor: Moritz Hardt Abstract In this work, we propose differentially private methods for hypothesis testing, model averaging, and model selection for normal linear models. We propose Bayesian methods based on mixtures of g-priors and non-Bayesian methods based on likelihood-ratio statistics and information criteria. The procedures are asymptotically consistent and straightforward to implement with existing software. We focus on practical issues such as adjusting critical values so that hypothesis tests have adequate type I error rates and quantifying the uncertainty introduced by the privacy-ensuring mechanisms. Keywords: Confidential Data, Regression, Bayesian Methods, Information Criteria 1. Introduction Differential privacy (Dwork et al., 2006) is a formal framework for quantifying the privacy of randomized algorithms. Its theoretical properties have been studied extensively (see e.g. Dwork et al. (2014)) and it has been adopted by companies such as Google, Apple, and Microsoft (Garfinkel et al., 2018), as well as institutions like the U.S. Census Bureau (Abowd, 2018). In this article, we develop differentially private methods for normal linear models. We propose differentially private hypothesis tests for comparing nested models (in Section 4) as well as methods for model averaging and selection (in Section 5). We consider Bayesian methods based on mixtures of g-priors (Liang et al., 2008) and non-Bayesian methods that are built upon likelihood-ratio statistics and information criteria. In Bayesian hypothesis testing and model selection, prior distributions must be chosen carefully because their effect does not vanish as the sample size grows (Bayarri et al., 2012). Our work is based on mixtures of g-priors because, when combined with right-Haar priors on the common parameters, they satisfy a list of appealing criteria proposed in Bayarri et al. (2012). They are also conveniently implemented in the Rpackage library(BAS) (Clyde, 2020). c 2024 V´ıctor Pe˜na, Andr´es F. Barrientos. License: CC-BY 4.0, see https://creativecommons.org/licenses/by/4.0/. Attribution requirements are provided at http://jmlr.org/papers/v25/21-1536.html. Pe˜ na and Barrientos From a non-Bayesian perspective, we work with likelihood-ratio tests and information criteria. We make this choice because of their theoretical properties, intuitive appeal, and ease of use. The class of information criteria we consider includes the Akaike Information Criterion (AIC; Akaike (1974)) and the Bayesian Information Criterion (BIC; Schwarz (1978)), among others. Although we focus on normal linear models, information criteria are useful in model averaging and selection problems for more complex models. We recommend the monograph Claeskens et al. (2008) for an extensive overview of the approach. We enforce differential privacy with well-established techniques. For hypothesis testing, we use the subsample and aggregate technique (Nissim et al., 2007; Smith, 2011), which consists in splitting the data into subgroups and releasing perturbed averages. For model averaging and selection, we use sufficient-statistic perturbation, which consists in releasing a noisy version of a sufficient statistic (see, for example, McSherry and Mironov (2009); Vu and Slavkovic (2009) and Bernstein and Sheldon (2019)). 1.1 Related Work There is a growing literature on differentially private methods for linear regression models. For example, Amitai and Reiter (2018) estimate quantiles and posterior tail probabilities of coefficients, Barrientos et al. (2019) test the significance of individual regression coefficients, and Ferrando et al. (2022) propose methods for point and interval estimation that can be applied to the normal linear model. Lei et al. (2018) consider the problem of model selection based on information criteria, but do not consider model averaging or Bayesian approaches, and Bernstein and Sheldon (2019) propose a method for sampling from posterior distributions on regression coefficients, but do not consider hypothesis testing, model averaging or selection. Our methods for hypothesis testing involve data-splitting and censoring. Both operations can induce bias in the outputs. Evans et al. (2020) considers the effects of censoring in estimates that are asymptotically normal and proposes strategies to correct the bias induced by censoring. Covington et al. (2021) uses bags of little bootstraps (Kleiner et al., 2012) to find unbiased estimates and valid confidence intervals. Alternatively, Ferrando et al. (2022) use the parametric bootstrap for bias correction. In hypothesis testing, the bias induced by data-splitting and censoring leads to inappropriate critical values for rejecting null hypotheses. We address this issue by simulating the distribution of the differentially private test statistics under the null hypothesis and then find corrected critical values that are adequately calibrated. In general, our Bayesian methodology draws from the objective Bayesian literature for hypothesis testing, averaging, and selection, especially from Liang et al. (2008) and Bayarri et al. (2012). In these references, the classes of priors we consider here are proposed after showing that they satisfy a list of conceptually appealing criteria. 1.2 Main Contributions In Section 4, we argue that working on a logarithmic scale is natural for defining differentially private Bayes factors. We show that the methods are asymptotically consistent under regularity conditions that are similar to the ones needed for consistency when there are no privacy constraints. In Section 4.2, we describe a simple procedure to quantify the effects 2 Differentially private methods for managing model uncertainty in linear regression of the privacy-ensuring mechanisms. Then, in Section 4.3, we study the effects of censoring and data-splitting. In Section 5, we use sufficient-statistic perturbation to define methods for model averaging and selection. If we took a naive approach to the problem, the variance of the perturbation term would grow exponentially in the number of predictors. With sufficientstatistic perturbation, the variance of the perturbation term grows quadratically in the number of predictors. The methods for model selection are consistent under conditions that are similar to those needed for consistency without privacy constraints. In Section 5.1, we propose a strategy that quantifies the uncertainty introduced by the mechanisms that is analogous to the one pursued in Section 4.2. From a practical point of view, we give guidelines for maximizing the statistical utility of the methods in finite samples. We illustrate the performance of our methods in Sections 4.4 and 5.2. We include additional results from the simulation studies in the Appendix. The proofs of the Propositions stated in the main text can be found in the Appendix A. 1.3 Notation Extrema appear often in Section 4 because the methods are based on censored statistics. Our notation for them is a∨b= max(a, b) and a∧b= min(a, b). The censored statistics are of the form Tc= (T∨a)∧b: that is, Tc=awhenever T≤aand Tc=bwhenever T≥b(with the understanding that a<b). We use the following notation for probability distributions. The p-variate normal distribution with mean µand covariance matrix Σ is Np(µ, Σ), the Laplace distribution with location parameter µand scale parameter bis L1(µ, b), and the (p×p)-dimensional Wishart distribution with degree of freedom n>p−1 and positive-definite scale matrix Vis Wp(n, V ). In the case of matrices, all of them are assumed to be full-rank unless otherwise stated. Our notation for basic matrix operations and special matrices is as follows. The matrix transpose of Ais A0, the perpendicular projection operator onto the column space of Ais PA=A(A0A)−1A0, and the upper Cholesky factor of Ais A1/2. The (n×n)−dimensional zero matrix is 0n×nand the (n×n)-dimensional identity matrix is In. For vectors, the usual p-norm on R dis k·kp. The zero vector is 0n= (0,0, ... n),0)0and the vector of ones is 1n= (1,1, ... n),1)0. The expectation of a random variable Xis E(X) and the variance is Var(X). 2. Brief Review of Differential Privacy Differential privacy (Dwork et al., 2006, 2014) is a probabilistic property of randomized algorithms. Intuitively, differential privacy limits how much we can learn about individual entries in a data set. In the literature, randomized algorithms that ensure differential privacy are referred to as mechanisms. Conceptually, mechanisms are functions Mthat take data Das inputs and output a random M(D) that is, in some sense, private. In the definition of differential privacy, a key notion is that of neighboring data sets. Two data sets Dand ˜ Dare neighbors, which is denoted by D∼˜ D, if they only differ in one row. 3 Pe˜ na and Barrientos Definition 1 ((ε, δ)-differential privacy) Let ε > 0and 0≤δ≤1. A mechanism M satisfies (ε, δ)-differential privacy if, for all D∼˜ Dany M-measurable set S, we have P[M(D)∈S]≤exp(ε)P[M(˜ D)∈S] + δ. When δ= 0, we retrieve ε-differential privacy, which is the most popular formalization of privacy in the literature. When εis small, the privacy of Mincreases; in such cases, the probability distributions of M(D) and M(˜ D) are forced to be similar, so the output of Mis not very sensitive to small changes in D. On the other hand, as εincreases, the distributions of M(D) and M(˜ D) are allowed to be more different, which decreases the privacy of the output. When 0 < δ ≤1, Definition 1 allows M(D) to release outputs that lead to a “high privacy loss” with probability δ(see Section 2 in Dwork et al. (2014) for details). In this article, M(D) are perturbed versions of confidential statistics T(D). The scale of the perturbation depends on the global sensitivity of T(D), a concept we define below. Definition 2 (Global sensitivity) Let Tbe a statistic mapping data to R d. The global sensitivity of Tis defined as ∆p= supD∼˜ DkT(D)−T(˜ D)kp. A key property of differential privacy is that transformations of (ε, δ)-differentially private statistics are (ε, δ)-differentially private. In the literature, this is known as the postprocessing property of differential privacy. We use this property frequently; for example, we use it to find posterior probabilities of hypotheses given Bayes factors. 3. Brief Review of Bayes Factors and Information Criteria In this section, we review basic facts about Bayes factors and information criteria that are relevant for our purposes. This is not meant to be a comprehensive review; we refer the reader to Berger and Pericchi (2001), Liang et al. (2008), and Bayarri et al. (2012) for further background on Bayesian methods and Claeskens et al. (2008) for information criteria. 3.1 Hypothesis Testing In Section 4, we introduce differentially private methods for hypothesis testing. We cover Bayesian approaches that are based on Bayes factors and non-Bayesian approaches that are based on likelihood-ratio tests and information criteria. Now, we review some core concepts in Bayesian and non-Bayesian testing that are helpful for contextualizing our work. Let y= (y1, ... , yn) be a vector collecting independent and identically distributed observations from a statistical model with sampling density f(y|θ) = Qn i=1 f(yi|θ), for θ∈Θ. Our goal is testing H0:θ∈Θ0against H1:θ∈Θ1, where Θ0and Θ1are disjoint subsets of Θ. In the Bayesian paradigm, all unknowns have probability distributions associated to them, including the hypotheses H0and H1and the parameter θ. Before observing any data, the uncertainty in the hypotheses is reflected in the prior probabilities P(H0) and P(H1) = 1 −P(H0). The uncertainty in θis usually expressed conditionally through the prior distributions π(θ|H0) and π(θ|H1). Upon observing data, the uncertainty about θ 4 Differentially private methods for managing model uncertainty in linear regression and the hypotheses is updated in the posterior distribution, which is simply the conditional distribution of these unknowns given the data. A key quantity in Bayesian hypothesis testing is the Bayes factor, which is defined as B10 =Zθ∈Θ1 π(θ|H1)f(y|θ)ν(dθ)/Zθ∈Θ0 π(θ|H0)f(y|θ)ν(dθ) for a dominating measure ν(·). Bayes factors are important in Bayesian testing for several reasons. For example, if o10 =P(H1)/P(H0) are the prior odds of the hypotheses, then the posterior probability of H1is P(H1|y) = o10B10/(1 + o10B10), which depends on the data only through the Bayes factor B10. Bayes factors can be motivated as ratios of integrated likelihoods that quantify the evidence in favor of H1relative to H0(Berger et al., 1999). There have been efforts to categorize the strength of evidence in favor or against H1based on the magnitude of B10 alone. Perhaps the most popular categorization is “Jeffreys’ scale of evidence” (Jeffreys (1939), see Table 1). Table 1: Jeffreys’ scale of evidence for Bayes factors (Jeffreys, 1939) Bayes factor Interpretation B10 <1H0supported 1< B10 <101/2Evidence against H0, but not worth more than a bare mention 101/2< B10 <10 Evidence against H0substantial 10 < B10 <103/2Evidence against H0strong 103/2< B10 <102Evidence against H0very strong B10 >102Evidence against H0decisive Logarithms of Bayes factors are naturally symmetric about zero. To see this, assume that H0and H1are equally likely a priori. Then, the posterior probability of H1is the standard logistic function in log B10: that is, P(H1|D) = 1/[1+exp(−log B10)]. Changing the sign of log B10 leads to the complement 1 −P(H1|D), and P(H1|D)=1/2 if and only if log B10 = 0. In Bayesian hypothesis testing, the priors π(θ|H0) and π(θ|H1) must be chosen carefully. Vague priors on θ, which are commonplace in estimation problems, can lead to posterior probabilities that overwhelmingly support H0no matter what the data are. This phenomenon is known in the literature as Lindley’s paradox (see Lindley (1957) and Robert (2014)), and it can occur when H1represents a larger set than H0. In normal linear regression problems, mixtures of g-priors have been studied carefully in Liang et al. (2008) and Bayarri et al. (2012). They avoid Lindley’s paradox by having a fixed scale matrix and they have other desirable properties such as invariance of Bayes factors with respect to changes of measurement units and large-sample consistency. In Sections 4 and 5, we work with priors within this class. From a non-Bayesian perspective, we propose working with likelihood-ratio tests and information criteria. The likelihood ratio of H1to H0is defined as 5 Pe˜ na and Barrientos Λ10 = max θ∈Θ1 f(y|θ)/max θ∈Θ0 f(y|θ). The likelihood ratio Λ10 is similar to the Bayes factor B10, but instead of averaging the likelihoods with weights given by π(θ|H1) and π(θ|H0), the likelihoods are maximized under H1and H0. Likelihood ratio tests are standard within the field of statistics and their properties are well-characterized (see, for example, Lehmann and Romano (2005)). Under mild conditions, the transformed log-likelihood ratio 2 log Λ10 is asymptotically distributed as chi-squared under H0and consistent under H1. Finally, we define the class of information criteria I10 =n−ρ/2Λ10, which encompasses AIC for ρ= 2p/ log n, BIC for ρ=p, and the likelihood ratio statistic Λ10 for ρ= 0. For a fixed ρ, the information criterion I10 can be interpreted as a penalized likelihood ratio, where ρacts as a penalty for model complexity. For the normal linear model, likelihood-ratio tests and AIC fail to be consistent under H0; BIC, on the other hand, is consistent under H0and H1. We refer the reader to Chapter 4 in the monograph Claeskens et al. (2008) for a more general version of this result and a discussion on related issues. 3.2 Model Averaging and Selection In Section 5, we focus on model averaging and selection. The context is a regression problem where there is an outcome variable Y∈Rnand ppredictors collected in a design matrix X∈Rn×p. We do not know which variables in X, if any, should be included in our model for Ygiven X. This type of uncertainty is often referred to as model uncertainty. For an introduction to the topic with a strong Bayesian flavor, we recommend Draper (1995). Model uncertainty can be parameterized through a binary vector γ∈ {0,1}pthat indicates active predictors: γi= 0 if the ith predictor is not active and γi= 1 if it is. Conceptually, we assume that there is a true model generating the data identified by T∈ {0,1}p. Given finite data D, we are do not know what the truth (T) is. From a Bayesian perspective, we can put a prior on γand find posterior probabilities to quantify this uncertainty; from a non-Bayesian perspective, we can find a point estimate of γor average our uncertainty over it with rules inspired by Bayesian procedures. For each model, which we identify by its active predictors in γ∈ {0,1}p, we can compute a Bayes factor or information criterion relative to the null model, which does not include any active predictors. We denote these null-based Bayes factors and information criteria Bγ0and Iγ0, respectively. With these, we can perform model selection after maximizing over γor we can find model-averaged estimates with weights proportional to Bγ0or Iγ0. 4. Hypothesis Testing In this section, we describe differentially private methods for testing a null hypothesis H0 against an alternative H1. The methodology described here can be applied in general, but in Section 4.1 we focus on nested linear regression models, for which we have found theoretical guarantees. 6 Differentially private methods for managing model uncertainty in linear regression Our methods are based on the subsample and aggregate technique (Nissim et al., 2007). It consists in splitting the data into Mdisjoint subgroups, computing statistics within the subgroups, and averaging the results. The output is made differentially private by adding a perturbation term η. The variance of ηis increasing in the global sensitivity (i.e., the range) of the statistics involved. We apply the subsample and aggregate technique as follows. First, we split the data into Mdisjoint subgroups of sample size b1, b2, ... , bM(PM i=1 bi=n) and compute censored statistics Tc i= (Ti∨L)∧U∈[L, U] for i∈ {1,2, ... , M}, where Tiare raw statistics computed from confidential data. After censoring the Ti, we know that they have global sensitivity ∆ = U−L. Finally, we release the noisy average PM i=1 Tc i/M +η, where ηis a random perturbation term that ensures (ε, δ)-differential privacy. If δ= 0, then η∼ L1(0,∆/(Mε)); if 0 < δ ≤1, then η∼N1(0, a∆2/(2Mε)), where ais a constant that can be computed with Algorithm 1 in Balle and Wang (2018). The confidential statistics Tiare logarithms of Bayes factors or logarithms of information criteria. We justify working on a logarithmic scale for Bayes factors in the subsequent paragraphs. A similar argument can be used to justify this choice for information criteria. Let Bc 10,i be the censored Bayes factor in the ith subgroup. Naively, one could release the noisy average PM i=1 Bc 10,i/M +η. Unfortunately, that approach has undesirable properties. Under both the Laplace and analytic Gaussian mechanisms, ηis supported on R and symmetric around zero, whereas B10 is always non-negative and shows equal support to H0 and H1when B10 = 1. If there are no privacy constraints, Bayes factors satisfy B01 =B−1 10 , but in general PM i=1 Bc 01,i/M +η6= (PM i=1 Bc 10,i/M +η)−1. Alternatively, we propose working on a logarithmic scale, defining log ˜ B10 = M X i=1 log Bc 10,i/M +η, log Bc 10,i = (log B10,i ∨L)∧U. Logarithms of Bayes factors are supported on R and, as we argued in Section 3, they have a natural symmetry around zero, so it is sensible to add a zero mean perturbation term on that scale. After exponentiating, we obtain a geometric mean of censored Bayes factors with a multiplicative perturbation: ˜ B10 = exp(η) M Y i=1 Bc 10,i!1/M . Since the distribution of ηis symmetric, ( ˜ B10)−1is equal in distribution to ˜ B01, and it is exactly equal to ˜ B01 when η= 0 (i.e., when there are no privacy constraints). Geometric means of Bayes factors have appeared in the objective Bayesian literature in geometric intrinsic Bayes factors (Berger and Pericchi, 1996). The statistic ˜ B10 is based on censored Bayes factors that are supported on [L, U]. However, the support of ˜ B10 is not [L, U] after introducing the perturbation term η. This issue can be solved by censoring ˜ B10, defining B∗ 10 = ( ˜ B10 ∨L∗)∧U∗,(1) 7 Pe˜ na and Barrientos where L∗= exp(L) and U∗= exp(U). With B∗ 10, we can define the differentially private posterior probability of H1given Das P∗(H1|D) = [1 −P(H0)]B∗ 10/{P(H0) + [1 − P(H0)]B∗ 10}. Following the same reasoning, we can define a differentially private information criterion I∗ 10 = (˜ I10 ∨L∗)∧U∗(2) ˜ I10 = exp(η) M Y i=1 Ic 10,i!1/M Ic 10,i = (b−ρ/2 iΛ10,i ∨L)∧U. 4.1 Nested Linear Regression Models Consider the normal linear model Y=X0β0+Xβ +σW, W ∼Nn(0n, In), where X0∈ R n×p0and X∈ R n×pare full-rank and n>p+p0. In this section, we present differentially private methods for testing H0:β= 0pagainst H1:β6= 0p. The set of predictors in X0is common to H0and H1and it can be, for instance, an intercept 1n. For the Bayesian methods, we define our priors after reparameterizing the model. We rewrite the model as Y=X0ψ+V β +σW for V= (In−PX0)Xand ψ=β0+ (X0 0X0)−1X0 0Xβ. In this parameterization, the common predictors X0are orthogonal to V, which is specific to H1. If X0is an intercept 1n, the reparameterization simply centers the predictors in X. Our prior specification is π(ψ, σ2)∝1/σ2, π(β|σ2, H1) = Z∞ 0 Np(β|0p, gσ2(V0V)−1)π(g)ν(dg). where νis an appropriate dominating measure. We allow the prior measure on gto depend on nand p, but not on Y. The prior on (ψ, σ2) is the right-Haar prior for this problem, which is improper, whereas β|σ2, H1is a mixture of g-priors. This class of priors has strong theoretical support (see Liang et al. (2008) and Bayarri et al. (2012) for details). For example, it leads to Bayes factors and posterior probabilities of hypotheses that are invariant with respect to invertible linear transformations of the design matrix V, such as changes of units. This property would not be satisfied if we had chosen a diagonal covariance matrix for β|g, σ2, H1. The prior distribution for β0and σ2is improper, so the marginal distributions can only be defined up to arbitrary multiplicative constants. We use the same constants for both H0 and H1, which is justified by the principle of invariance (Berger et al., 1998) and predictive matching arguments (Bayarri et al., 2012). When the data are not confidential, we can report the Bayes factor (Liang et al., 2008): B10 =Z∞ 0 (g+ 1)(n−p−p0)/2[1 + g(1 −R2)]−(n−p0)/2π(g)ν(dg),(3) 8 Differentially private methods for managing model uncertainty in linear regression where R2=Y0PVY/Y 0(In−PX0)Yor, from a non-Bayesian perspective, we can report the information criterion I10 =n−ρ/2Λ10 =n−ρ/2(1 −R2)−n/2.(4) If there are privacy constraints, B10 or I10 cannot be released directly. We propose applying the subsample and aggregate technique to release differentially private versions of B10 and I10. That is, we split up the data set into Mdisjoint subgroups of sample sizes b1, b2, ... , bM(PM i=1 bi=n), define priors πi(gi) and penalties ρifor i∈ {1,2, ... , M}, and release B∗ 10 and I∗ 10 as defined in Equations (1) and (2), respectively. Proposition 1 below states that B∗ 10 and I∗ 10 are consistent under some regularity conditions. Consistency does not follow from Liang et al. (2008) for various reasons. First, there is a growing number of subgroups and the sample sizes within the subgroups grow to infinity. Another difference is that the Bayes factors are censored and there is a perturbation term η. Our result does not follow directly from Smith (2011), either. Smith (2011) studies the asymptotic behavior of averages PM i=1 Ti/M released with the subsample and aggregate method. If T∗ iis the “true value” of Ti, Lemma 6 in Smith (2011) assumes that √bi(Ti−T∗ i) is asymptotically normal, E(Ti)−T∗ i∈O(1/bi), and E{[√bi(Ti−T∗ i)]3} ∈ O(1).Furthermore, the results in Smith (2011) are for independent and identically distributed data. In our case, it is unclear whether the asymptotic conditions are satisfied (if they are, it would require proof), and we would need additional assumptions on, at least, the design matrices, the priors πi(gi), and the penalties ρi. Instead of verifying the conditions in Smith (2011), our proof uses a union bound and tail inequalities for R2. Proposition 1 Under the regularity conditions listed below, B∗ 10 and I∗ 10 are consisten: under H0,B∗ 10 →P0and I∗ 10 →P0; under H1,(B∗ 10)−1→P0and (I∗ 10)−1→P0. 1. Well-specified model: The raw confidential data are generated from the normal linear model described in this section. 2. Growth of Mand bi:limn→∞ M=∞and limn→∞ infi∈1:Mbi=∞. If p≥2, limn→∞ supi∈1:MM/bi= 0 and, if p= 1,limn→∞ supi∈1:MMplog bi/bi= 0. 3. Censoring limits: limn→∞ L=−∞,limn→∞ U=∞. 4. Privacy parameters: The privacy parameters εand δare such that limn→∞(U− L)/(Mε) = limn→∞ a(U−L)2/(Mε) = 0. 5. Design matrices: Under H1,limn→∞ infi∈1:Mβ0 TX0 iXiβT/(σ2 Tbi)>0, where βTand σ2 Tare the fixed true values of βand σ2, respectively. 6. Priors on giand penalties ρi: lim n→∞ sup i∈1:MZ∞ 0 bp/2 i(gi+ 1)−p/2πi(gi)ν(dgi)<∞, ρi≤max(p, log bi) lim n→∞ inf i∈1:MZ∞ bi bp/2 i(gi+ 1)−p/2πi(gi)ν(dgi)>0, ρi≥p. 9 Pe˜ na and Barrientos where βγ∈ R |γ|×1is a vector including the βjsuch that γj= 1 and Vγ∈ R n×|γ|is a matrix with the active variables in γ. Just as we did in Section 4, we parameterize the model so that X0and Vγare orthogonal. For convenience, we denote the null model where none of the variables are active as γ= 0. The matrix X0represents a set of predictors we are sure to include in our model. Usually, X0contains an intercept 1n, but it can be empty as well. If X0is empty, the methods in this section would be defined in analogous manner after replacing In−PX0by In. From a Bayesian perspective, our prior specification on the regression coefficients βγis the same we had in Section 4: we put mixtures of g-priors on βγ|σ2, γ and the right-Haar prior π(ψ, σ2)∝1/σ2on the common parameters. If γis not the null model, the expression for the non-private Bayes factor of model γto the null model, denoted Bγ0, is identical to the one in Equation (3) after substituting pby |γ|and R2by R2 γ=Y0PVγY/Y 0(I−PX0)Y. If γis the null model, we have Bγ0= 1. Given a prior distribution π(γ) on γ, the posterior probability of the model identified by γgiven the data Dis defined as P(γ|D) = P(γ)Bγ0/X ˜γ∈{0,1}p P(˜γ)B˜γ0, which depends on the data only through the Bayes factors Bγ0. We assume that the prior P(γ) can depend on p, but not on nor the design matrix. Most common choices for P(γ), like a uniform prior P(γ) = 2−por the hierarchical uniform prior recommended by Scott and Berger (2010) satisfy the condition. From a non-Bayesian perspective, the information criteria I∗ γ0are as in Equation 4 after substituting R2by R2 γand ρby ρ|γ|, where ρ|γ|is an increasing function in |γ|such as ρ|γ|=|γ|. We could use the method in Section 4.1 to release differentially private versions of all the Bγ0or Iγ0. However, if we released those 2pstatistics, the variance of the private statistics would increase exponentially in p. Instead, we propose working with a perturbed version of a sufficient statistic whose dimension increases quadratically in p. Let Z= (I−PX0)Y∈ R n×1be the null-centered outcome variable and V∈ R n×pbe the design matrix of the full model that includes all ppredictors, which we collect in a data matrix D= [V;Z]. Assuming Z0Z > 0 almost surely, we define G=D0D=V0V V 0Z Z0V Z0Z=V0V V 0Y Y0V Y 0(I−PX0)Y. The Gram matrix Gis a sufficient statistic for the normal linear model (see, for example, Seber and Lee (2012)). As a consequence, all of the R2 γ,Bγ0, and Iγ0can be constructed by taking appropriate subsets of G. We propose releasing a noisy version of the sufficient statistic G, a technique known in the differential privacy literature as sufficient-statistic perturbation (see, for example, McSherry and Mironov (2009); Vu and Slavkovic (2009) and Bernstein and Sheldon (2019)). We construct a differentially private version of Gby adding a random perturbation term, defining G∗=G+E, where Eis a random perturbation matrix that ensures differential privacy. To ensure ε-differential privacy, we use the Laplace mechanism. For (ε, δ)-differential 16 Differentially private methods for managing model uncertainty in linear regression privacy with 0 <ε<1 and 0 < δ ≤1, we use Algorithm 2 in Sheffet (2019), which we refer to as the Wishart mechanism. To establish the parameters of the distribution of E, we assume that there are lower and upper bounds for the data: that is, there are land usuch that each entry dij in Dis within the interval [l, u]. Since the response and predictors are centered, l < 0 and u > 0. The entries of G=D0Dare of the form Pn i=1 dijdij0. Here, d1jis the first row and jth column of D. If we replace this entry by another one, say ˜ d1j0, then the maximum absolute difference in the entries of Gis |d1jd1j0−˜ d1j˜ d1j0|. When computing the sensitivity ∆1, we only consider the entries where j≥j0because Gis symmetric. Therefore, the sensitivity ∆1can be upper-bounded as follows: ∆1= sup D∼ e DkD0D−˜ D0˜ Dk1 = sup D∼ e DX j≥j0|d1jd1j0−˜ d1j˜ d1j0| =X j≥j0 sup D∼ e D|d1jd1j0−˜ d1j˜ d1j0| ≤(p+ 1)(p+ 2) sup d11, d110∈(u,l)|d11d110| ≤(p+ 1)(p+ 2)(l2∨u2). For the Laplace mechanism, we define the perturbation term Eas a symmetric random matrix of the form E=     e11 e12 ... e1,p+1 e12 e22 ... e2,p+1 . . .. . ..... . . e1,p+1 e2,p+1 ... ep+1,p+1      , where ejj0iid ∼L1(0,∆1/ε), for j≥j0, and ej0j=ejj0. For the Wishart mechanism (i.e., Algorithm 2 in Sheffet, 2019), we need a uniform bound on the Euclidean norm of the rows of D. Given the assumption l < dij < u, we can use p(p+ 1)(l2∨u2) as a uniform bound. In this case, given the bound, we define the random perturbation as E=M−E(M), where M∼Wp+1(k, (p+ 1)(l2∨u2)Ip+1) and k=bp+ 1 + 28 log(4/δ)/ε2cwith 0 < ε < 1 and 0< δ ≤1. For both mechanisms, the variances of the entries in Ecan be high, especially for small values of εor δ. This, in turn, can lead to outputs that overestimate the number of active predictors. To avoid this issue, we propose post-processing G∗in two ways: hardthresholding off-diagonal elements as in Bickel et al. (2008) and adding a constant rto the diagonal elements. We propose thresholding the off-diagonal entries of G∗at eλ, the λ-th percentile of Eij. More formally, we define G∗∗ to be the matrix with typical element G∗ ij 1 (i=jor |G∗ ij| ≥ eλ). This choice can be justified as follows: if an off-diagonal entry of G∗ ij =Gij +Eij is not an extreme value in the distribution of Eij, it is likely that G∗ ij is essentially Eij and Gij is nearly zero. Adding a constant rto the diagonal elements of G∗can be seen as a ridge-type of regularization that can reduce the variability of the output. It can also be used to guarantee 17 Pe˜ na and Barrientos that the output is positive-definite because, after adding the Laplace perturbation terms, G∗ may not be positive-definite. In fact, the hard-thresholded matrix G∗∗ need not be positivedefinite even if G∗is positive-definite (Bickel et al., 2008). Given a symmetric matrix A, which here can be either G∗or G∗∗, the matrix Ar=A+rIp+1 is positive-definite as long as r > −eigmin(A), where eigmin(·) is a function returning the minimum eigenvalue of a matrix. In essence, we propose releasing a differentially private Gram matrix G∗, which we then transform to obtain a better estimate of G. Then, we use the formulas we would use if we had access to Gafter plugging in its differentially private estimate. In our applications (in Section 5.2), we consider methods that hard-threshold G∗and methods that do not. That is, we compare methods that are based on G∗∗ r=G∗∗+rIp+1 to methods based on G∗ r=G∗+rIp+1, where G∗is not hard-thresholded. In our experience, hard-thresholding is helpful when the ground truth is sparse, but it can be detrimental when most predictors are active. We justify this argument in Section 5.2. Proposition 2 establishes model-selection consistency under some assumptions. More precisely, we show that the differentially private Bayes factor of any model γto the true model, which can be expressed as B∗ γT =B∗ γ0/B∗ T0, converges to zero in probability for any γ6=T. The differentially private information criteria I∗ γT are also consistent under the assumptions listed below. We proved the result by characterizing the asymptotic behavior of G∗ r,G∗∗ r, and R2,∗ γ, and then bounding the Bayes factors and information criteria above and below. Just as we had in Section 4, the proof covers εand (ε, δ) differentially private methods. The proof does not follow from Liang et al. (2008) for several reasons, one of them being that we are not assuming that the response given the covariates is normal, since we assume that the data are bounded. Proposition 2 Let T∈ {0,1}pbe a vector indexing the truly active predictors. As n→ ∞, and under the regularity conditions listed below, B∗ γT →P0and I∗ γT →P0for any γ6=T. 1. Boundedness: The data Dare within the interval [l, u]for finite land u. 2. Regression mean and variance: E(Y|X0, V ) = X0ψ+VTβTand Var(Y| X0, V ) = σ2 TIn, where VT∈ R n×pTis a matrix that contains the truly active predictors. 3. Privacy parameters: The privacy parameters εand δare such that the matrix perturbation term E/n →P0. 4. Regularization parameters: λis fixed and ris so that limn→∞ r/n = 0. 5. Design matrices: limn→∞ V0V/n =S1, where S1is symmetric and positive-definite. 6. Priors on gand penalties ρ:for all 1≤ |γ| ≤ p, the prior π(g)satisfies lim n→∞Z∞ 0 np/2(g+ 1)−|γ|/2π(g)ν(dg)<∞ lim n→∞Z∞ n np/2(g+ 1)−|γ|/2π(g)ν(dg)>0, 18 Differentially private methods for managing model uncertainty in linear regression or from a non-Bayesian perspective, ρ|γ|is increasing in |γ|and satisfies |γ| ≤ ρ|γ|≤ |γ|∨log n. The framework here is different to the one in Section 4. We assume that the confidential data are bounded (Assumption 1) and do not assume that the distribution of the outcome given the predictors is normal. In other words, Bayes factors and information criteria are misspecified beyond the addition of the perturbation term E. Nonetheless, Proposition 2 shows that the methods are consistent. Our setup is also distinct to the one adopted in Lei et al. (2018), where it is simultaneously assumed that the response is normal (Assumption 1 in Lei et al. (2018)) and bounded (Assumption 4). Assumption 2 requires that the regression mean be well-specified and that the covariance of the response given the predictors be spherical. Assumption 3 forces the perturbation matrix Eto be so that E/n →P0. As we had in Section 4, it suffices to let εand δbe fixed for it to hold. Assumption 4 imposes conditions on the regularization parameters. The case of a non-thresholded matrix G∗ ris included as λ= 0. Assumption 5 is similar to the regularity condition on design matrices in Proposition 1. Finally, Assumption 6 is essentially the same as the assumptions on the priors and penalties in 2. Just as we had in Section 4, the differentially private Bayes factors with Zellner’s g-prior with g=n, the robust prior, Zellner-Siow, and BIC are all consistent, and so is BIC. In the common scenario where X0is an intercept 1n, the methods described in this section can be conveniently implemented in Rwith the bas.lm function in library(BAS) (Clyde, 2020). Given a set of predictors and an outcome variable, the bas.lm function enumerates Bayes factors for small to moderate pand samples from the model space for large p. The function outputs other statistics of interest such as posterior inclusion probabilities and model-averaged estimates. To use bas.lm for our problem, we need to generate a synthetic data set D(containing both centered predictors and outcome) whose sufficient statistic D0Dis equal to a fixed Gram matrix G,which can be G∗ ror G∗∗ r. Proposition 3 below shows how to obtain such a synthetic data set D= [V;Z] given a matrix G. Proposition 3 Let U ∈ R n×(p+1) be a full-rank matrix and define M= (In−1n10 n/n)U. Given a matrix G, we can generate a synthetic data set D= [V;Z]with the formula D= M(M0M)−1/2G1/2. The synthetic data set Dsatisfies the identities D0D=G,V01n= 0p, and Z01n= 0. Proposition 3 guarantees that the outputs we obtain from running bas.lm on the synthetic data Dare identical to what we would find by taking subsets of Gdirectly. In the proposition, the matrix Uis arbitrary; in practice, its entries can be simulated by sampling independently from the uniform distribution. 5.1 Quantifying the Uncertainty Introduced by the Mechanism With the non-thresholded methods, the private statistics are of the form G∗ r=G+E+rIp+1. Since rand the distribution of Eare both known, we can define a confidence set for the nonprivate Ggiven G∗ r. With such a set, it is possible to find confidence regions for summaries of interest T(G) like least-squares estimates or inclusion probabilities. To define a 1−αconfidence set for G, we first find E1−αsuch that P(E∈ E1−α)=1−α. For a matrix norm k·k, define E1−α={E:kE−E(E)k ≤ q1−α}where P[kE−E(E)k ≤ 19 Pe˜ na and Barrientos q1−α] = 1 −α. Then, ˜ C1−α={G∗ r−rIp+1 −E:E∈ E1−α}is a 1 −αconfidence set for G. Since Gis symmetric and positive-definite, we can intersect ˜ C1−αwith the set of symmetric positive-definite matrices S++ to define a 1−αconfidence region C1−α=˜ C1−α∩S++ whose volume is at most that of ˜ C1−α. The confidence set C1−αcan be transformed into T(C1−α) to produce confidence sets for summaries of interest T(G). We can approximate T(C1−α) with a rejection sampler. First, simulate E1, E2, ... , Ensim from the appropriate mechanism and compute kEi−E(E)kfor i∈ {1,2, ... , nsim}. Then, approximate the ((1 −α)×100)th percentile q1−αwith its empirical version ˆq1−α= inf{q: Pnsim i=1 1 (kEi−E(E)k ≤ q)/nsim ≥1−α}and define ˆ E1−α={Ei:kEi−E(E)k ≤ ˆq1−α, i = 1, . . . , nsim}. After that, we can find T(ˆ C1−α); that is, compute T(G∗ r−rIp+1 −Ei) for Ei∈ˆ E1−αand only keep those such that G∗ r−rIp+1 −Ei∈ S++. This computational strategy to construct T(C1−α) is an application of the parametric bootstrap: the quantity to be inferred is the unknown confidential summary T(G), and the only source of randomness is the noise injected into Gto make it differentially private. In general, the confidence set need not be an interval, but we can summarize the confidence set with a histogram. To do so, we define the bins of the histogram as Bk= [tk−1, tk) with min{T(ˆ C1−α)}=t0< t1< . . . < tK−1< tK= max{T(ˆ C1−α)}and their corresponding relative frequencies #Bkusing the number of elements of T(ˆ C1−α) that fall in Bk,k= 1, . . . , K. We denote the histogram summarizing T(ˆ C1−α) as Hist(T, ˆ C1−α) = {(Bk,#Bk)}K k=1. We can also report histograms in a density scale; that is, Hist(T, ˆ C1−α) = {(Bk, dk)}K k=1, where dk=#Bk |T(ˆ C1−α)|(tk−tk−1) and |T(ˆ C1−α)|denotes the cardinality of T(ˆ C1−α). While Hist(T, ˆ C1−α) displays the distribution of the estimator of T(G) constrained to T(ˆ C1−α), its support, denoted by [t0, tK], corresponds to an 1 −αconfidence interval. This is because t0= min T(ˆ C1−α) and tK= max T(ˆ C1−α). Therefore, analysts can directly report [t0, tK] as the confidence interval, if desired. Additionally, analysts can use the histogram to assess whether [t0, tK] provides a good representation of the set T(ˆ C1−α). They can also consider if T(ˆ C1−α) would be better represented by the union of non-adjacent intervals. 5.2 Empirical Evaluations We evaluate the performance of the methods described in this section in a simulation study and an application. The simulation study is similar to the one in Liang et al. (2008), whereas the real data set is a subset of the March 2000 Current Population Survey that was analyzed in Barrientos et al. (2019). We implement the methods with the Rpackage BAS (Clyde, 2020). In the simulation study, we found Bayes factors with the Zellner-Siow prior (ZS) and information criteria with BIC. The prior distribution on the model space π(γ) is the hierarchical uniform prior proposed in Scott and Berger (2010). From a non-Bayesian perspective, π(γ) acts as a function that weighs the information criteria. The results with Zellner-Siow and BIC are almost identical. We report the outputs based on the former here and show the results with the latter in Appendix C. 20 Differentially private methods for managing model uncertainty in linear regression We compare the results we obtain by hard-thresholding and not thresholding the Gram matrix G∗. In all cases, we add a regularization parameter rto the diagonal entries of G∗. For the Laplace mechanism, we set rto be the 99-th percentile of eigmin(E), which we find via simulation. For the Wishart mechanism, we use the analytical expression in Remark 2 of Sheffet (2019). 5.2.1 Simulation Study We simulate data from a normal linear model with ppredictors, where pis set to 2, 6, or 9. The sample size n(in thousands) varies from 5 to 10,000. The number of active predictors in the true model |T|depends on the value of pand ranges from 0 (null model is true) to p(full model is true). Specifically, if p= 2, we set |T| ∈ {0,1,2}; if p= 6, we set |T|∈{0,3,6}; and if p= 9, we set |T|∈{0,4,9}. The predictors are independently drawn from the uniform distribution on (−2,2). Following Hastie et al. (2017), we define the signal-to-noise ratio (SNR) as the variance of the regression mean (which is random, since we are simulating predictors and β) divided by σ2. In our simulations, we assume that the intercept is zero and βis a p-dimensional vector equal to b[1,...,1]0. We use optimization to find σ2and bsuch that SNR = 0.5 and the response falls within (−2,2) with high probability. For each combination of |T|and n, we simulate 1,000 data sets. All the data sets we simulated are such that the response falls in (−2,2). We consider ε∈ {0.5,0.9}and, in the case of the Wishart mechanism, we set δ= 1/n. We assess the performance of the methods by tracking Monte Carlo averages of predictive mean squared errors and the posterior probability of the true model. We define the predictive mean squared error as PMSE = n−1kVTβT−V β∗k2 2,where VTis a design matrix containing truly active predictors, βTis the true value of β, and β∗is the differentially private model-averaged posterior expectation. Figure 3 displays PMSEs for different values of n,p,ε, and |T|. As expected, the PMSEs for both private and non-private approaches decrease as the sample size increases. We observe that the PMSEs of the differentially private methods are smaller when εis 0.9 compared to when εis 0.5, and they are always higher than the non-private PMSEs. This is expected, as larger values of εshould lead to greater statistical utility. In most cases, the methods based on the Laplace mechanism have a lower PMSE than those based on the Wishart mechanism. However, the Wishart mechanism seems to perform better in the case when pis either 6 or 9, |T|is zero, and the sample size is small. Although this is slightly less evident in the figure, upon closer inspection, we can see that methods based on hard-thresholding tend to have a slightly lower PMSE when |T|is zero but cease to be advantageous when |T|is large and the sample size is small. Figure 4 displays the posterior probabilities of the true model, for the same values of n, p,ε, and |T|that we used in Figure 3. The results are consistent with what we observed for PMSEs. We also observe that, although all probabilities increase as the sample size increases, the rate at which they increase depends on p. For higher values of p, the rate of increase is lower because the computational complexity of the problem increases with the dimension of G. The number of entries in Gincreases quadratically in p, and so does the variance of the perturbation term added to ensure differential privacy. This fact affects the convergence rate of the posterior probabilities. 21 Pe˜ na and Barrientos p = 9 |T| = 0 p = 9 |T| = 4 p = 9 |T| = 9 p = 6 |T| = 0 p = 6 |T| = 3 p = 6 |T| = 6 p = 2 |T| = 0 p = 2 |T| = 1 p = 2 |T| = 2 10 100 1000 10000 10 100 1000 10000 10 100 1000 10000 1e−12 1e−08 1e−04 1e−12 1e−08 1e−04 1e−12 1e−08 1e−04 n (thousands) PMSE with Bayesian model averaging Method/Mechanism G (not DP) Lap Non−thresh (DP) Wish Non−thresh (DP) Lap Thresh (DP) Wish Thresh (DP) ε 0.5 0.9 Figure 3: Simulation study: Sample size (x-axis) against log(PMSE) (y-axis) with ZellnerSiow prior. 5.2.2 Application: Current Population Survey The data set includes n= 49,436 heads of households with non-negative incomes. We consider 6 predictors: age in years (β1), age squared (β2), marital status (β3), sex (β4), education (β5), and race (β6). All predictors are numeric or binary except for education, which is an ordinal variable. To reduce the number of coefficients in the model, we treat education as numeric, ranging from 1 (for less than 1st grade) to 16 (for doctoral degree). The binary predictors are: marital status (1: civilian spouse present; 0: otherwise), sex (1: male; 0: female), and race (1: white; 0: otherwise). The response variable is income. In this application, the non-private inclusion probabilities are all close to one. To provide a more challenging benchmark for our methods, we permute the rows for marital status and education in the design matrix to artificially make the inclusion probabilities for β3and β5close to zero. The predictors and the response are centered and rescaled to the interval (−0.5,0.5). Figure 5 displays the posterior expected values of β1,β3, and β4with the Zellner-Siow prior and ε= 0.9. We use the histograms described in Section 5.1 to define approximate 95% confidence sets for T(G) = E(βj|G). Our choice of matrix norm is the Frobenius norm. Specifically, we run our procedure 250 times and, for each run and a fixed collection of bins B1,...,BK, we summarize each T(ˆ C1−α) with its corresponding histogram Hist(T, ˆ C0.95) = {(Bk, dk)}K k=1. 22 Differentially private methods for managing model uncertainty in linear regression p = 9 |T| = 0 p = 9 |T| = 4 p = 9 |T| = 9 p = 6 |T| = 0 p = 6 |T| = 3 p = 6 |T| = 6 p = 2 |T| = 0 p = 2 |T| = 1 p = 2 |T| = 2 10 100 1000 10000 10 100 1000 10000 10 100 1000 10000 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 n (thousands) Posterior probability true model Method/Mechanism G (not DP) Lap Non−thresh (DP) Wish Non−thresh (DP) Lap Thresh (DP) Wish Thresh (DP) ε 0.5 0.9 Figure 4: Simulation study: Sample size (x-axis) against posterior probability of the true model (y-axis) with Zellner-Siow prior. If we let d(l) 1, . . . , d(l) Kbe the densities values of the histogram associated with the l-th run, l= 1,...,250, we define average histograms as Hist T(·) = E(βj| ·),ˆ C0.95=( Bk, dk=1 250 250 X l=1 d(l) k!)K k=1 , j ∈ {1,3,4}. The results are displayed in Figure 5. In all cases, the differentially private methods are close to the non-private answers. We can also see that the histograms are useful for quantifying the uncertainty introduced by the mechanism, since their spread increases when the thresholded and non-thresholded methods do not agree in their estimates. 5.3 Guidelines Thresholded methods tend to perform best when the true number of predictors is small. On the other hand, when most predictors are active, non-thresholded methods tend to outperform thresholded methods. The root cause behind this phenomenon is that thresholding shrinks the elements in G∗, which promotes sparsity. In practice, we recommend that users run analyses with both thresholded and nonthresholded methods. This can be done without affecting the privacy budget of the analyst 23 Pe˜ na and Barrientos posterior expected value of β1 density −2 −1 0 1 2 0 1 2 3 4 Laplace mechanism posterior expected value of β3 density −0.010 0.000 0.010 0 500 1000 1500 2000 Laplace mechanism posterior expected value of β4 density 0.00 0.01 0.02 0.03 0.04 0 20 40 60 80 Laplace mechanism posterior expected value of β1 density −2 −1 0 1 2 0 1 2 3 4 Wishart mechanism posterior expected value of β3 density −0.010 0.000 0.010 0 500 1000 1500 2000 Wishart mechanism posterior expected value of β4 density 0.00 0.01 0.02 0.03 0.04 0 20 40 60 80 Wishart mechanism Figure 5: Current population survey: Posterior expectations of coefficients β1, β3, and β4 with Zellner-Siow prior and ε= 0.9. The vertical solid lines are the non-private posterior expectations, whereas the dashed lines and dotted lines are averaged posterior expectations estimated with non-thresholded methods and thresholded methods, respectively. because both G∗ rand the thresholded matrix G∗∗ rare post-processed versions of the same differentially private matrix G∗. Finally, we recommend reporting confidence sets whenever possible, since we find them to be a valuable tool for quantifying the uncertainty introduced by the mechanisms. 6. Conclusions and Future Work In this article, we proposed differentially private methods for hypothesis testing, model averaging, and model selection for normal linear models. Under regularity conditions, the methods are consistent. The regularity conditions we have imposed are similar to the conditions used in the literature for establishing consistency of non-differentially private methods. Our methods for hypothesis testing are based on data-splitting and censoring statistics. We have studied the effects of these operations on the performance of the methods. In the case of data-splitting, increasing the number of subsets reduces the variance of the differentially private statistics, but it adds bias. In the case of censoring, more stringent 24 Differentially private methods for managing model uncertainty in linear regression censoring reduces the variance, but it can lead to substantial bias if the true, confidential statistic lies outside the uncensored range. The methods we proposed for model averaging and selection are based on a perturbed sufficient statistic. If we suspect that the ground truth is sparse, we recommend hardthresholding the perturbed sufficient statistic; however, if most predictors are active, hardthresholding can lead to underfitting. The methodology proposed here could be extended in a number of ways. It would be useful to extend the methods to generalized linear models through the framework proposed in Li and Clyde (2018). The implementation for hypothesis testing is straightforward, but our methods for model uncertainty, which are based on sufficient statistics, cannot be applied directly. This obstacle can be overcome using approximate sufficient statistics, as proposed in Huggins et al. (2017). This approach has been used successfully in estimation problems under differential privacy constraints in Kulkarni et al. (2021). It would also be interesting to extend the methods to survival models because health records are confidential. In this case, the extension could adapt the framework proposed in Castellanos et al. (2021). In our work, we have used off-the-shelf techniques for establishing differential privacy. While we observe that our proposals can be useful in practice, it might be possible to design more efficient mechanisms that are specifically tailored to the tasks we considered. This is another interesting avenue for future research. Acknowledgments The authors would like to thank the feedback from two anonymous reviewers that greatly improved the presentation and contents of the article. The research of the second author was supported by the National Science Foundation National Center for Science and Engineering Statistics [49100420C0002 and 49100422C0008] and the Test Resource Management Center (TRMC) within the Office of the Secretary of Defense (OSD), contract #FA807518D0002. 25 Pe˜ na and Barrientos We can show that Q/(n−p0)→Pσ2 Twith the Hanson-Wright inequality, which in turn implies that Q/n →Pσ2 T. We have E[Q/(n−p0)] = σ2 T. Define W=Y−E(Y). Then, E(W) = 0nand W2 iare uniformly bounded since both Yand E(Y) are by Assumption 1. Therefore, we can pick a finite constant Ksatisfying E[exp(W2 i/K2)] <2. Applying the Hanson-Wright inequality, for any given  > 0, P[|Q/(n−p0)−σ2 T|> ]≤2e−2c min[(K2k(I−PX0)/(n−p0)kF)−1,(k(I−PX0)/(n−p0)kop)−1]→n→∞ 0, because k(I−PX0)/(n−p0)kop =k(I−PX0)/(n−p0)kF= 1/√n−p0→0 as n→ ∞. This implies Q/n →Pσ2 Tand Z0Z/n →Pσ2 T+β0 TS3βT. We have shown that G/n →PG/n =S1S2βT β0 TS0 2σ2 T+β0 TS3βT=G∞. Therefore, we have G∗/n →PG∞. In the case of the non-thresholded matrix G∗ r, we have established that G∗ r/n =G∗/n + r/nIp+1 →PG∞since r/n →0 by Assumption 4. In the case of G∗∗ r, there is an indicator that can hard-threshold off-diagonal elements. Let G∗∗ ij /n =G∗ ij/n 1 (i=jor |G∗ ij|/n ≥ eλ/n) be the (i, j)-th entry of G∗∗. On the one hand, limn→∞ eλ/n = 0 because the variance of Eij is finite and does not depend on n. For i=jthe indicator is equal to 1. For i6=j: E[ 1 (|G∗ ij|/n ≥eλ/n)] = P[|G∗ ij|/n ≥eλ/n]→1 Var[ 1 (|G∗ ij|/n ≥eλ/n] = P[|G∗ ij|/n ≥eλ/n]−P[|G∗ ij|/n ≥eλ/n]2→0 1 (|G∗ ij|/n ≥eλ/n)→P1. By Slutzky’s lemma, we have that G∗∗ ij /n →G∞,ij for all i, j. This is enough to show that G∗/n →G∞. Finally, since G∗∗ r=G∗∗ +rIp+1 and r/n →0 by Assumption 4, we have G∗∗ r/n →PG∞, as required. Proposition 6 (Convergence of noisy R2) Given G∗ ror G∗∗ rand a model-indexing vector γ∈ {0,1}p, we can construct a differentially private version of R2 γ, denoted R2,∗ γ. Let T∈ {0,1}pbe the index of the true model. Under the regularity conditions stated in Proposition 5 and assuming that Z0Zis not equal to zero almost surely, the following are true: 1. If γdoes not nest T(i.e. if there exists isuch that Ti= 1 and γi= 0), then R2,∗ γ→R2 γ,∞and R2,∗ T→R2 T,∞with R2 γ,∞< R2 T,∞. 2. If the T= 0pis the null model, then nR2,∗ γis in Op(1). 3. If γnests T(i.e. if Ti= 1 implies γi= 1), then [(1 −R2,∗ T)/(1 −R2,∗ γ)](n−p0)/2is in Op(1). 32 Differentially private methods for managing model uncertainty in linear regression Proof We prove the three statements separately. Proof of statement 1. The proofs of R2 γ→PR2 γ,∞and R2,∗ γ→PR2 γ,∞are a direct consequence of Proposition 5. Note that R2 γ=Z0Vγ(V0 γVγ)−1V0 γZ Z0Z=Z0PVγZ/n Z0Z/n . On the one hand, n/Z0Z→P1/(σ2 T+β0 TS3βT). Then, V0 γZare subvectors of V0Z, so V0 γZ/n converges to a subvector of S2βT. Similarly, V0 γVγ/n converges to a submatrix of S1 and, since we assume that Vis full-rank, V0 γVγis invertible for any γ. This implies that R2 γconverges in probability to some constant R2 γ,∞, which has to be between zero and one because Vγ(V0 γVγ)−1Vγis a projection matrix and Z0PXZ≤Z0Zfor any projection matrix PX. The convergence of the noisy R2,∗ γto R2 γ,∞can be established after noting that R2,∗ γcan be constructed by taking submatrices of G∗ ror G∗∗ rand multiplying and dividing terms as needed. We can invoke Proposition 5 and Slutzky’s lemma and conclude that R2,∗ γ→PR2 γ,∞, as required. It remains to show that if γdoes not nest the true model T,R2 γ,∞< R2 T,∞. To see this, note that lim n→∞ E(Z0PVγZ/n) = lim n→∞β0 TX0 TPVγXTβT/n, lim n→∞Var(Z0PVγZ/n) = 0, and β0 TX0 TPVγXTβT< β0 TX0 TXTβT=β0 TX0 TPVTXTβT, which implies R2 γ,∞< R2 T,∞. Proof of statement 2. If Tis the null model, we show that nR2,∗ γis in Op(1). It is useful to write nR2,∗ γ= 1 √n(Z0V)∗ γh(V0V)∗ γ ni−1(V0Z)∗ γ1 √n (Z0Z)∗/n . By Proposition 5, we know that the denominator converges in probability to a constant. It is enough to show that the numerator is in Op(1). The matrix (V0V)∗ γ/n−1is in Op(1) by Proposition 5. It remains to show that (V0Z)∗ γ/√nis in Op(1). Let w∗= (V0Z) + E2and kλthe appropriate threshold for the off-diagonal elements of G∗∗. Then, k1 √n(V0Z)∗ γk2≤1 nkw∗ 1 (|w∗|> kλ)k2 ≤1 nkw∗k2 ≤1 nkV0Zk2+1 nkE2k2. It suffices to show that both the error E2/√nand the non-private V0Z/√nare in Op(1). Each of the Eij in Ehas a finite variance that does not depend on n. Therefore, each Eij/√nhas a variance that goes to zero and, since E2has pelements, kE2/√nkis in Op(1). It only remains to show that V0Z/√nis in Op(1), which we show using the Hanson-Wright inequality. 33 Pe˜ na and Barrientos First, note that kV0Z/√nk2=Z0V V 0Z/n. Then, since we are assuming that the true model is the null model: E(Z0V V 0Z/n) = σ2 ntr(V0V)→n→∞ c > 0. The limit is a positive constant because V0V/n converges to a symmetric positive-definite matrix by Assumption 5, which also implies limn→∞kV V 0/nkF<∞and limn→∞kV V 0/nkop < ∞.From here, we can apply the Hanson-Wright inequality to establish that |Z0V V 0Z/n − E(Z0V V 0Z/n)|is in Op(1) and, since E(Z0V V 0Z/n) converges to a constant, we conclude that Z0V V 0Z/n is in Op(1), as required. Proof of statement 3. Finally, we show that if γnests T, [(1−R2,∗ T)/(1−R2,∗ γ)](n−p0)/2 is in Op(1). Consider the non-private [(1 −R2 T)/(1 −R2 γ)](n−p0)/2. We can write 1−R2 T 1−R2 γ(n−p0)/2 =1 + 2 n−p0 Z0(PVγ−PVT)Z 2Z0(In−PVγ)Z/(n−p0)(n−p0)/2 ≤exp Z0(PVγ−PVT)Z 2Z0(In−PVγ)Z/(n−p0). The denominator Z0(In−PVγ)Z/(n−p0) converges to a constant. This fact follows directly given the asymptotic behavior of Z0PVγZ/n and Z0Z/n we just described in the proof of the first statement. The same is true for the private version of the statistic, using the argument we used for R2,∗ γ. The numerator Z0(PVγ−PVT)Zis in Op(1), which can be shown using the HansonWright inequality. The expectation is E[Z0(PVγ−PVT)Z] = σ2(|γ| − |T|). Let Qγ= Z0(PVγ−PVT)Z−σ2(|γ|−|T|). Applying the Hanson-Wright inequality P(Qγ> M)≤2 exp −cM K2min M K2(|γ|−|T|),1, where Kis a constant that can be chosen in a similar way as we did in Proposition 5. The right-hand side can be made arbitrarily close to zero by increasing M, which establishes that Qγand Z0(PVγ−PVT)Zare in Op(1). The same is true for the private version of the statistic, using an argument which is similar to the one we used for showing that nR2,∗ γis in Op(1) when the true model is the null model. Both (V0Z)∗ γ/√nand (V0Z)∗ T/√nconverge to their private versions (V0Z)γ/√nand (V0Z)T/√nbecause the error goes to zero and the indicator converges (see Proposition 5 and the earlier proof for nR2,∗ γfor more detailed versions of these arguments). The private (V0V)∗ γalso converge to their non-private versions (see e.g. Proposition 5). Therefore, the private (Z0(PVγ−PVT)Z)∗has the same asymptotic behavior as the private Z0(PVγ−PVT)Z, which we have shown to be in Op(1). This completes the proof. A.4 Proof of Proposition 2 in main text We consider two cases: one where the true model is the null model and another one where the true model is not the null model. 34 Differentially private methods for managing model uncertainty in linear regression True model is the null model: Let γbe a model that is not the null model. Then, we have B∗ γ0≤exp n−p0 2 R2,∗ γ 1−R2,∗ γ!Z∞ 0 (g+ 1)−|γ|/2π(g)ν(dg). The integral converges to zero by Assumption 6 and the exponential term is in Op(1) because, by Proposition 6, we know that nR2,∗ γis in Op(1). Therefore, for any model γ which is not the null model, B∗ γ0→P0. A similar argument works for information criteria. In such case, I∗ γ0≤n−ρ|γ|/2exp n 2 R2,∗ γ 1−R2,∗ γ!→P0. True model is not the null model: We study the asymptotic behavior of B∗ γT , where Tis the true model and γis a model that is not the true model. We split this case into two subcases: one where γnests the true model, and another one where γdoes not nest the true model. First, note that, since we are working with null-based Bayes factors, B∗ γT =B∗ γ0(1 −R2,∗ T)(n−p0)/2 B∗ T0(1 −R2,∗ T)(n−p0)/2 we will bound the numerator and denominator separately and put our bounds together. First, we bound the numerator: B∗ γ0(1 −R2,∗ T)(n−p0)/2=Z∞ 0 (g+ 1)−|γ|/2"1 + g(1 −R2,∗ T)−R2,∗ T 1 + g(1 −R2,∗ γ)#(n−p0)/2 π(g)ν(dg) ≤ 1−R2,∗ T 1−R2,∗ γ!(n−p0)/2Z∞ 0 (g+ 1)−|γ|/2π(g)ν(dg). Then, we bound the denominator: B∗ T0(1 −R2,∗ T)(n−p0)/2≥Z∞ n (g+ 1)−|T|/2"1−R2,∗ T 1 + g(1 −R2,∗ T)#(n−p0)/2 π(g)ν(dg) ≥"1−R2,∗ T 1 + n(1 −R2,∗ T)#(n−p0)/2Z∞ n (g+ 1)−|T|/2π(g)ν(dg). Putting the bounds together and using Assumption 6: B∗ γT ≤ 1−R2,∗ T 1−R2,∗ γ!(n−p0)/2 1−R2,∗ T 1 + n(1 −R2,∗ T)!−(n−p0)/2R∞ 0(g+ 1)−|γ|/2π(g)ν(dg) R∞ n(g+ 1)−|T|/2π(g)ν(dg) .n(|T|−|γ|)/2 1−R2,∗ T 1−R2,∗ γ!(n−p0)/2 exp(R2,∗ γ/(1 −R2,∗ γ)). 35 Pe˜ na and Barrientos When γnests T,n(|T|−|γ|)/2goes to zero and the remaining terms are in Op(1) by Proposition 6, so B∗ γT converges to zero in probability. When γdoes not nest T, we show that An=n(|T|−|γ|)/2[(1 −R2,∗ T)/(1 −R2,∗ γ)](n−p0)/2→P0 Let Rn= (1 −R2,∗ T)/(1 −R2,∗ γ), which converges in probability to a constant less than one. Taking logarithms log nn(|T|−|γ|)/2R(n−p0)/2 no=n−p0 2[(|T|−|γ|) log n/(n−p0) + log Rn], which diverges to −∞ in probability, so An→P0 The remaining term in the upper bound for B∗ γT is in Op(1). Therefore, we have shown that B∗ γT →P0. The proof for information criteria is essentially the same, but we do not need to bound integrals. Indeed, I∗ γT =n(ρ|T|−ρ|γ|)/2 1−R2,∗ T 1−R2,∗ γ!n/2 , and we can use the same arguments we used for B∗ γT to show that I∗ γT is consistent. A.5 Proof of Proposition 3 in main text Let U ∈ R n×(p+1) be a full-rank matrix and define M= (In−PX0)U. Given a Gram matrix G, we can generate a synthetic data set D= [V;Z] with the formula D=M(M0M)−1/2G1/2. In other words, we have D0D=G: D0D=G1/2M(M0M)−1/2M0M(M0M)−1/2G1/2=G1/2G1/2=G. The synthetic data are also centered (the same way that Vand Zare centered in our construction of D). This is true because Dis pre-multiplied by In−PX0, so it is orthogonal to the span of X0= 1n. Appendix B. Effects of censoring and data-splitting: High School and Beyond Survey In this section, we revisit the High School and Beyond Survey data set (Section 4.4) to study the effects of setting different censoring limits. For concreteness, we restrict our attention to the test of H02 against H12; that is, the hypothesis test where we wish to know if read scores are predictive of math scores when science scores are already a covariate in the model. The results for the test of H01 against H11 are similar. The simulation setup is the same as described in Section B. The only changes are the censoring limits. In Figure 6, we display the results for the differentially private posterior probability of H12 for different values of εand M. We censor the posterior probability of H12 at [0.35,0.65], [0.25,0.75], [0.01,0.99], and [0.001,0.999]. The true, non-private posterior probability of H12 is near 1. For all censoring limits, increasing the number of subgroups decreases the variance of the output, but it induces a bias that shrinks the probability to 0.5. More stringent censoring limits reduce the variance as well, but come at the cost of potential 36 Differentially private methods for managing model uncertainty in linear regression bias: for instance, censoring the posterior probability at [0.35,0.65] is clearly too stringent, since the true posterior probability is much higher than the upper limit. In Figure 7, we display the analogous result for likelihood ratios. In this case, the lower censoring limit is set to L= 0, which is a natural lower bound for likelihood ratios. The gray + and Mare calibrated critical values for rejection of H02 at significance at levels 0.01 and 0.05, respectively. The upper censoring limits Uare set to 1, 2, 7, and 10. We arrive at the same conclusions we reached with the Bayesian analysis. The true, non-private likelihood ratio is above 10, and if we censor at a much lower value (such as 1 or 2) the test fails to be powerful. The test performs best if the lower bound is 10, in which case the test is quite powerful, especially for M≥5 and δ≥0.25. P(H12| D) in [0.001, 0.999] M = 2 P(H12| D) in [0.001, 0.999] M = 5 P(H12| D) in [0.001, 0.999] M = 10 P(H12| D) in [0.01, 0.99] M = 2 P(H12| D) in [0.01, 0.99] M = 5 P(H12| D) in [0.01, 0.99] M = 10 P(H12| D) in [0.25, 0.75] M = 2 P(H12| D) in [0.25, 0.75] M = 5 P(H12| D) in [0.25, 0.75] M = 10 P(H12| D) in [0.35, 0.65] M = 2 P(H12| D) in [0.35, 0.65] M = 5 P(H12| D) in [0.35, 0.65] M = 10 0.25 0.50 1.00 2.00 4.00 8.00 0.25 0.50 1.00 2.00 4.00 8.00 0.25 0.50 1.00 2.00 4.00 8.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 .. P(H12| D) Figure 6: Distribution of P∗(H12 |D) as a function of ε,M, and censoring limits. The lower endpoint of the error bars is the first quartile of the distribution, the midpoint is the median, and the upper endpoint is the third quartile. The dashed lines are the non-private posterior probabilities. 37 Pe˜ na and Barrientos U = 10 M = 2 U = 10 M = 5 U = 10 M = 10 U = 7 M = 2 U = 7 M = 5 U = 7 M = 10 U = 2 M = 2 U = 2 M = 5 U = 2 M = 10 U = 1 M = 2 U = 1 M = 5 U = 1 M = 10 0.01 0.25 0.50 0.75 0.01 0.25 0.50 0.75 0.01 0.25 0.50 0.75 0.0 2.5 5.0 7.5 10.0 0.0 2.5 5.0 7.5 10.0 0.0 2.5 5.0 7.5 10.0 0.0 2.5 5.0 7.5 10.0 δ 2log*Λ10,2 Figure 7: Distribution of 2 log Λ∗ 10,1and 2 log Λ∗ 10,2as a function of δ,M, and censoring upper limit U. The gray + and Mare corrected critical values at the 0.01 and 0.05 significance levels, respectively. Appendix C. Additional Plots for Simulation Study In this section, we include additional plots for the simulation study in Section 5.2. We show results with BIC combined with least-squares estimates, as well as results for additional values of εthat were not included in the main text. The interpretation of the plots is the same: thresholded methods perform best when |T|is small, and non-thresholded methods perform best when |T|is large. 38 Differentially private methods for managing model uncertainty in linear regression p = 9 |T| = 0 p = 9 |T| = 4 p = 9 |T| = 9 p = 6 |T| = 0 p = 6 |T| = 3 p = 6 |T| = 6 p = 2 |T| = 0 p = 2 |T| = 1 p = 2 |T| = 2 10 100 1000 10000 10 100 1000 10000 10 100 1000 10000 1e−11 1e−07 1e−03 1e−11 1e−07 1e−03 1e−11 1e−07 1e−03 n (thousands) PMSE with Bayesian model averaging Method/Mechanism G (not DP) Lap Non−thresh (DP) Wish Non−thresh (DP) Lap Thresh (DP) Wish Thresh (DP) ε 0.5 0.9 Figure 8: Simulation study: Sample size (x-axis) against log(PMSE) (y-axis) with BIC prior. 39 Pe˜ na and Barrientos p = 9 |T| = 0 p = 9 |T| = 4 p = 9 |T| = 9 p = 6 |T| = 0 p = 6 |T| = 3 p = 6 |T| = 6 p = 2 |T| = 0 p = 2 |T| = 1 p = 2 |T| = 2 10 100 1000 10000 10 100 1000 10000 10 100 1000 10000 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 n (thousands) Posterior probability true model Method/Mechanism G (not DP) Lap Non−thresh (DP) Wish Non−thresh (DP) Lap Thresh (DP) Wish Thresh (DP) ε 0.5 0.9 Figure 9: Simulation study: Sample size (x-axis) against posterior probability of the true model (y-axis) with BIC prior. 40 Differentially private methods for managing model uncertainty in linear regression References John M Abowd. The US census bureau adopts differential privacy. In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 2867–2867, 2018. Milton Abramowitz, Irene A Stegun, et al. Handbook of mathematical functions, volume 55. Dover New York, 1964. Hirotugu Akaike. A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6):716–723, 1974. Gilad Amitai and Jerome Reiter. Differentially private posterior summaries for linear regression coefficients. Journal of Privacy and Confidentiality, 8(1), 2018. Borja Balle and Yu-Xiang Wang. Improving the Gaussian mechanism for differential privacy: Analytical calibration and optimal denoising. In International Conference on Machine Learning, pages 394–403. PMLR, 2018. Andr´es F Barrientos, Jerome P Reiter, Ashwin Machanavajjhala, and Yan Chen. Differentially private significance tests for regression coefficients. Journal of Computational and Graphical Statistics, 28(2):440–453, 2019. Maria J Bayarri, James O Berger, Anabel Forte, Gonzalo Garc´ıa-Donato, et al. Criteria for Bayesian model choice with application to variable selection. The Annals of Statistics, 40(3):1550–1577, 2012. James O Berger and Luis R Pericchi. The intrinsic Bayes factor for model selection and prediction. Journal of the American Statistical Association, 91(433):109–122, 1996. James O Berger and Luis R Pericchi. Objective Bayesian methods for model selection: Introduction and comparison. Lecture Notes-Monograph Series, pages 135–207, 2001. James O Berger, Luis R Pericchi, and Julia A Varshavsky. Bayes factors and marginal distributions in invariant situations. Sankhy¯a: The Indian Journal of Statistics, Series A, pages 307–321, 1998. James O Berger, Brunero Liseo, and Robert L Wolpert. Integrated likelihood methods for eliminating nuisance parameters. Statistical science, pages 1–22, 1999. Garrett Bernstein and Daniel R Sheldon. Differentially private Bayesian linear regression. Advances in Neural Information Processing Systems, 32:525–535, 2019. Peter J Bickel, Elizaveta Levina, et al. Covariance regularization by thresholding. The Annals of Statistics, 36(6):2577–2604, 2008. Lucien Birg´e. An alternative point of view on Lepski’s method. Lecture Notes-Monograph Series, pages 113–133, 2001. 41