The Influence of Misspecified Covariance on False Discovery
Control when Using Posterior Probabilities
Ye Liang1 , Joshua D. Habiger1 and Xiaoyi Min2

arXiv:1709.05269v1 [] 15 Sep 2017

1

2

Department of Statistics, Oklahoma State University, Stillwater, OK, USA
Department of Mathematics and Statistics, Georgia State University, Atlanta, GA, USA

Abstract
This paper focuses on the influence of a misspecified covariance structure on false
discovery rate for the large scale multiple testing problem. Specifically, we evaluate
the influence on the marginal distribution of local fdr statistics, which are used in
many multiple testing procedures and related to Bayesian posterior probabilities. Explicit forms of the marginal distributions under both correctly specified and incorrectly
specified models are derived. The Kullback-Leibler divergence is used to quantify the
influence caused by a misspecification. Several numerical examples are provided to
illustrate the influence. A real spatio-temporal data on soil humidity is discussed.
Keywords: Multiple testing, Bayes, Dependent data, Divergence, Spatio-temporal.

1

Introduction

Large scale multiple testing arises from many practical problems, from genetic studies
to public health surveillance. Benjamini and Hochberg (1995) introduced the concept of
the false discovery rate (FDR) and proposed a powerful testing procedure, usually referred as the BH procedure. The BH procedure relies on a positive dependence assumption (Benjamini and Yekutieli, 2001), while adaptive BH procedures (Storey et al., 2004;
Liang and Nettleton, 2012) rely on an independence or weak dependence structure. Efron
(2007) noted that correlation may result in overly liberal or overly conservative testing
procedures. Though the BH procedure is valid under different dependence assumptions
(Farcomeni, 2007; Wu, 2008), Sun and Cai (2009) showed that failing to model dependence
can result in inefficiency. To address that problem, Sun and Cai (2009) and Sun et al. (2015)
propose a procedure using local significance index, which is a Bayesian posterior probability. Efron et al. (2001) described the connection between the FDR and Bayes procedures,
where a posterior probability is referred as the local false discovery rate (Lfdr). Sun and Cai
(2009)’s local significance index reduces to the Lfdr under independence. In fact, there is
a rich history of using Bayesian approaches for multiplicity adjustment. Scott and Berger
1

(2006) and Scott and Berger (2010) discussed Bayesian multiplicity adjustment in variable
selections. Muller et al. (2006) had a comprehensive discussion on the connection between
the FDR, Bayesian multiple testing, and procedures using posterior probabilities.
In Sun et al. (2015), the procedure depending on unknown parameters is called the “oracle” procedure, in the sense that we know all nuisance parameters as an oracle. Although the
oracle procedure is proved to control FDR at the nominal level and be optimal in terms of
false non-discovery rate (FNR), a data driven “adaptive” procedure relies on correctly specifying the model, including the prior specification, and/or consistent parameter estimation.
For large-scale data, dependence or covariance, if in Gaussian models, is often estimated
based on a structured model. The model choice or the structure choice itself may be debatable and parameter estimation remains challenging. For example, in spatial modeling,
the estimation of covariance relies on structured covariance specifications and in practice,
one may have multiple choices of specifications. Intuitively, the choice of specification will
influence the data driven procedure and may eventually lead to different decisions.
In this paper, we explore the influence of a misspecified covariance structure on the testing
procedure. Specifically, we study the sampling distributions of the local fdr statistics under
both correctly and incorrectly specified covariance structures. We derive explicit expressions
for those distributions under a general model setting. We propose to use the Kullback-Leibler
divergence as a quantitative measure for the influence. We show in both a simulation study
and a real application that the influence of a misspecification leads to unappealing results.
The paper is organized as follows. Section 2.1 gives a basic setup of this problem. Section
2.2 gives sampling distributions of the test statistics. Section 2.3 provides formulas for
computing the Kullback-Leibler divergence. Section 2.4 shows several numerical examples.
Section 3 provides real data that arose from the Oklahoma soil monitoring network. Section
4 is a discussion.

2
2.1

Main Results
The Setup

In this paper, we consider a general model for m observations y = (y1 , . . . , ym)′ :
yi = θi + ǫi ,

i = 1, . . . , m,

(1)

where θi is the latent state and ǫi is the noise term which independently follows N(0, σ 2 ). The
dependence of observations is introduced through their latent states θ = (θ1 , . . . , θm )′ . For
2

instance, θi can be a realized spatial process θ(s), for which a spatial dependence structure
can be specified. Consider a one-sided hypothesis
H0i : θi ≥ θ0i

versus H1i : θi < θ0i ,

(2)

for every i simultaneously. This type of one-sided hypothesis is often of interest in many
practices. In spatial epidemiology, one may want to determine which regions or locations have
disease rates higher than some given threshold θ0 . Here, data will be spatially correlated
and each hypothesis will be one-sided. Similarly, in agricultural studies, one may want
to determine which locations or time periods have soil moisture levels lower than a given
threshold θ0 , indicating a risk of drought, and the data will be either spatially or temporally
correlated and each hypothesis will be one-sided. A more general hypothesis would be
H0i : θi ∈ Θ0i versus H1i : θi ∈ Θ1i . We do not consider a precise (or two-sided) hypothesis
in this paper, but some comments are given in the discussion section.
The dependence of latent states θ is usually specified through a prior model. For example,
consider a normal-inverse-gamma prior on (θ, σ 2 ),
θ | σ 2 ∼ Nm (θ0 , gσ 2Σ) and σ 2 ∼ IG(α, β).

(3)

For simplicity, we assume that Σ is a known covariance structure and g is a known scale
parameter. The use of g here is the same as that in Zellner’s g-prior for Bayesian variable
selection problems. The g value could be fixed, estimated or have a hyperprior (Liang et al.,
2008). The prior specification (3) in fact induces a marginal probability for each hypothesis:
P (H0i ) = P (θi ≥ θ0i ) = 0.5.

Sun and Cai (2009) and Sun et al. (2015) showed that, to control FDR when data are
dependent, the posterior probability hi = P (H0i | y) is useful. The posterior probability
hi is viewed as a test statistic, called local index of significance in their work. The oracle
procedure orders hm = (h(1) , . . . , h(m) ) and rejects all H(i) , i = 1, . . . , k such that
(
)
i
1X
k = max i :
(4)
h(j) ≤ α∗ ,
i j=1

where α∗ is the nominal level. The procedure mimics the Benjamini-Hochberg procedure,
in which p-value is the test statistic. Sun and Cai (2009) showed that this oracle procedure
controls FDR at level α∗ and has the smallest FNR among all FDR procedures at α∗ for
a hidden Markov model. Sun et al. (2015) further showed that in a spatial random field
model, this oracle procedure controls FDR at level α∗ and has the smallest missed discovery
rate (MDR). A data-driven procedure, however, depends on the estimation of other nuisance
3

parameters. The covariance Σ is especially important in this case as it describes the dependence. Our objective is to determine if the procedure is sensitive when the covariance is
incorrectly specified or estimated, and if so, to quantify the sensitivity.

2.2
2.2.1

Sampling Distribution of Test Statistics
Known variance of noise

We now focus on how the distribution of test statistics (h1 , . . . , hm ) is influenced by a misspecified covariance structure. Assume that the data are generated from the true underlying
process:
yi ∼ N(θi , σ02 ) and θ ∼ N(θ0 , Σ1 ),

(5)

for i = 1, . . . , m. Assume known σ02 and consider model (1) with priors θ ∼ N(θ0 , gΣ1) and
θ ∼ N(θ0 , gΣ2), the latter of which has a misspecified covariance structure. The scale g
determines the strength of the prior. The intuition is that both g and Σ2 will influence the
test statistics (h1 , . . . , hm ) and FDR control.
Lemma 1 Under the correct covariance Σ1 , hi marginally has the following CDF:

r
aii −1
Φ (hi ) ,
F (hi ) = Φ
bii

(6)

−1
where aii is the ith diagonal element in A = (1/σ02 I + 1/gΣ−1
and bii is the ith diagonal
1 )
−1 −1
−1 −1
2
2
2
element in B = (I + σ0 /gΣ1 ) (σ0 I + Σ1 )(I + σ0 /gΣ1 ) .
−1
−1
2
Under the misspecified covariance Σ2 , A = (1/σ02 I+1/gΣ−1
and B = (I+σ02 /gΣ−1
2 )
2 ) (σ0 I+
−1
Σ1 )(I + σ02 /gΣ−1
2 ) .

Lemma 1 shows explicitly how the sampling distribution is altered by a misspecification.
Observe that F (hi ) is completely determined by the ratio aii /bii . Figure 1 shows different
shapes of both CDF and pdf under different ratio values. Note that when aii = bii , the
sampling distribution is Uniform(0, 1). Consider g → ∞, for which the prior becomes non2
2
informative, then aii /bii → σ02 /(σ02 + σ1,ii
), where σ1,ii
is the ith diagonal element in Σ1 . The
ratio becomes irrelevant to the misspecified covariance Σ2 . In other words, the behaviors of
the correct specification and the incorrect specification will be similar when g is large, which
seems intuitive. Because FDR procedures based on hi reject H0i if hi ≤ C for some C, it
is important to note that P (Hi ≤ hi ) = F (hi ) can be substantially influenced by covariance
4

misspecification. For example, in Figure 1, we see F (0.2) ranges from 0.03 to 0.35. It should
be noted that some recommend rejecting H0i if hi ≤ 0.2 (Efron, 2012). In this setting, the
rejection probability ranges from 0.03 to 0.35.
Unlike the independent case, h1 , . . . , hm are now dependent when the data are dependent.
The change under a misspecified structure is revealed in their joint distribution.
Theorem 1 Using the same definition for A and B in Lemma 1, h1 , . . . , hm have a joint
CDF

r
r
a11 −1
amm −1
m
Φ (h1 ), · · · ,
Φ (hm ) ,
F (h1 , . . . , hm ) = Φb
b11
bmm
where Φm
b is the CDF for a multivariate normal Nm (0, Pb ), and Pb is the correlation matrix
of B.
The joint distribution of h1 , . . . , hm represents a multivariate surface in the space [0, 1]m .
Notice that jointly not only the ratio aii /bii plays a role but also the correlation structure of
B. Under a misspecified covariance structure, Pb will be altered as well. However, still, as
g → ∞, B → σ02 I + Σ1 , and there is no misspecification effect.
2.2.2

Unknown variance of noise

Assume the underlying data generating process (5). Suppose σ02 is unknown and consider
specifying model (1) with prior (3). As before, a correct covariance structure is Σ1 and a
misspecified structure is Σ2 .
Theorem 2 Under the correct covariance Σ1 , the test statistics h1 , . . . , hm jointly have the
following CDF:
r

r
a11 −1
amm −1
m
F (h1 , . . . , hm ) = Ξa,b
Ψ
(h1 ), · · · ,
Ψ
(hm ) ,
(7)
b11 m+2α
bmm m+2α
−1 −1
−1 −1
−1
2
where A = (I + 1/gΣ−1
1 ) , B = (I + 1/gΣ1 ) (σ0 I + Σ1 )(I + 1/gΣ1 ) , Ψm+2α is the
CDF for a univariate t-distribution with degrees of freedom m + 2α, and Ξm
a,b is the CDF for
the following random vector
s
m + 2α
zb ,
′
zb Czb + 2β

5

where zb ∼ N(0, Pb ) and C = (diag B)−1/2 (A−2 −A−1 )(diag B)−1/2 , where diag here denotes
a diagonal matrix.
−1
−1
2
Under the misspecified covariance Σ2 , A = (I + 1/gΣ−1
and B = (I + 1/gΣ−1
2 )
2 ) (σ0 I +
−1 −1
Σ1 )(I + 1/gΣ2 ) .
The main difference between Theorem 2 and Theorem 1 is that Ξm
a,b has a more complim
cated form than Φb . The ratio aii /bii still plays a role in the joint distribution and Ξm
a,b will
be affected by a misspecified structure.

2.3

Kullback-Leibler Divergence

Denote the sampling distribution of h under the correct covariance specification as fcor (h)
and under the misspecified covariance as fmis (h). We evaluate the influence of the misspecification by the Kullback-Leibler (KL) divergence
Z
fcor (h)
DKL (fcor k fmis ) = fcor (h) log
dh.
(8)
fmis (h)
The KL divergence here can be interpreted as the information loss when using fmis to approximate fcor . The following two corollaries are useful to approximate the KL divergence
in our cases.
Corollary 1 As a consequence of Lemma 1, the marginal density function of hi is


√
1
2
f (hi ) = ri exp
(1 − ri )φi ,
2
where ri = aii /bii and φi = Φ−1 (hi ).
Corollary 2 As a consequence of Theorem 1, the joint density function of (h1 , . . . , hm )′ is


1
1
1
1 2
1
1 ′
−1
−1
f (h1 , . . . , hm ) = R 2 Pb R 2 exp
φ (I − R 2 Pb R 2 )φ ,
2
where φ = (φ1 , . . . , φm )′ = (Φ−1 (h1 ), . . . , Φ−1 (hm ))′ and R = diag(r1 , . . . , rm ).
With Corollary 2, the term log{fcor (h)/fmis (h)} in expression (8) is analytically available.
Thus, the KL divergence DKL can be easily evaluated using Monte Carlo approximation.
Notice that we can draw from fcor exactly given the underlying model because hi = P (H0i |
6

y) is analytically available. To draw a sample of h from fcor , draw y from the underlying
model and then calculate h. Suppose we obtain a sample h1 , . . . , hL , the approximation is
L

DKL (fcor

fcor (hl )
1X
log
.
k fmis ) ≈
L l=1
fmis (hl )

A misspecified covariance will change the defined matrices A and B, and consequently change
matrices R and Pb . Notice the relationship R1/2 Pb−1 R1/2 = (diag1/2 A)B −1 (diag1/2 A).
Since the KL divergence will generally increase as the dimension m increases (in the independent case it is simply a sum of individual dimensions), we may also consider a relative
measure of influence DKL /m.
The KL divergence is computable under the general model (1) with Σ1 and Σ2 provided.
This easy-to-compute measure can be used to quantify the influence of a misspecification.
In practice, when there are multiple candidate covariances, we may assess the KL divergence
between those candidates.

2.4

Numerical Examples

We now consider two examples of misspecified covariance. In each example, without loss of
generality, we set m = 900 and σ02 = 0.25. We will numerically evaluate the KL divergence
and perform a simulation study, in which we estimate FDR and FNR with Monte Carlo
replications of 1000.
Example 1 Positive spatial covariance vs. Independence.
Consider a regular spatial grid with unit distance one for generating the latent states θ. The
2
true covariance Σ1 = {σ1,ij
} has a positive decaying structure determined by an exponential
2
covariance function σ1,ij = exp{−ksi − sj k/ρ} with ρ = 5, where s represents a location.
A misspecified covariance is Σ2 = I. To get a rough idea of ri = aii /bii , let g=1 and
compute A and B. Under the correct specification, ri ranges from 0.12 to 0.16, and under
the misspecification, ri = 0.25.
Note that this misspecification is essentially to ignore the dependence and treat data as
independent observations. This is quite common in practice, where domain scientists are
hesitant to model complex covariance structures, though evidences suggest that data may
be correlated.

7

Example 2 Negative AR(2) covariance vs. Positive AR(2) covariance.
Consider a time series for generating the latent states θ. The true covariance Σ1 is determined by an AR(2) process: θi = ρ1 θi−1 + ρ2 θi−2 + εi , εi ∼ N(0, 1), and ρ1 = 1.5 and
ρ2 = −0.9. The autocorrelation function of this specification has an oscillating pattern
(mixed positive and negative values in Σ1 ). A misspecified covariance Σ2 is chosen to be the
covariance for an AR(2) process with ρ1 = 0.6 and ρ2 = 0.3, whose autocorrelation is always
positive. To get a rough idea of ri = aii /bii , let g=1 and compute A and B. Under the
correct specification, ri ranges from 0.088 to 0.14, and under the misspecification, ri ranges
from 0.20 to 0.25.
Note that this example shows a scenario where model is correctly specified but parameter
estimates are wrong. This example also compares a covariance matrix containing negative
values with a covariance matrix containing all positive values.
Results of Example 1 and 2 are shown in Figure 2 and Figure 3. We choose different g
values representing different strengths of information brought in by the prior dependence.
Four plots are shown in each result: the estimated FDR, the estimated FNR, the difference
between the rejection rate of the correct specification and that of the misspecification, i.e.
(#discoverycor − #discoverymis)/m, and DKL /m. We can reach the following conclusions
from these plots. First, as g increases, both the KL divergence and the difference between
rejection rates decrease, and also both the FDR and the FNR become closer, all of which
are as expected, indicating a decreasing misspecification influence. Second, the FNR under
the misspecification is universally higher than under the correct specification. Hence, misspecification results in an inefficient procedure. Also notice that the procedure tends to give
less discoveries under the misspecification than under the correct specification. Last but not
least, the FDR change is not monotonic and the comparison between the two specifications is
profound. Notice that the nominal level is 0.05 and g = 1 (or log10 g = 0) represents a “true”
scale. When both the structure Σ1 and the scale g are correct, from the top left plot in both
results, we can see that the FDR is controlled at the nominal level. This seems to suggest
that a correctly estimated scale of prior dependence is desired, which should neither be too
strong nor too weak. And this correct scale will only work as expected if the covariance
structure is correct as well.
Example 3 Positive spatial covariance revisit.
Revisit Example 1. Consider to fix g = 1. The correct specification is exactly the same as
the underlying model with ρ = 5. For the misspecification, let the spatial range parameter ρ
change from 0.1 to 20, representing the strength of dependence.
8

Results of Example 3 are shown in Figure 4. The FDR is maintained at 0.05 only when
ρ has the correct value. The pattern of change in FDR is also reflected in the KL divergence
plot. In this example, using DKL /m as a measure of influence seems to be reasonable.
The FNR, on the other hand, monotonically decreases as the dependence goes stronger.
We shall note here that, when g is large, the KL divergence becomes sensitive to detect a
misspecification as the measure approaches zero. However, in that case, the prior is vague
and the influence on FDR is negligible.

3

Real Data: Soil Relative Humidity

Oklahoma Mesonet (Illston et al., 2008) is a comprehensive observatory network monitoring
environmental variables across the state. One of the focuses of the network is soil moisture.
Extreme weather conditions, especially drought, severely impact Oklahoma’s agriculture,
which is a leading economy of the state. Soil moisture is fundamentally important to many
hydrological, biological and biogeochemical processes. The information is valuable to a wide
range of government agencies and private companies. We take a small dataset from their
data warehouse as an example of real application. Consider only one variable here: the
relative humidity, ranging from 0 to 100 percent. The dataset consists of monthly averages
in 2014 for 108 monitoring stations, which is in total 1,296 measurements. Consider each
hypothesis being H0i : θi ≥ 50 versus H0i : θi < 50 for detecting low humidity times and
locations.
We consider two different specifications for the dependence structure. Consider a spatiotemporal process, for a spatial location s and a time point t, y(s; t) = θ(s; t) + ǫ(s; t), where
ǫ(s; t) is pure error process with N(0, σ 2 ), and θ(s; t) is a stationary Gaussian process with
a constant mean µ and a separable covariance function:
C(h; τ ) = δC (s) (h)C (t) (τ ),
where h = ks − s′ k and τ = |t − t′ | are both Euclidean distances. Specify C (s) (h) =
exp{−h/ρ} and C (t) (τ ) = ατ . Specify priors for parameters: σ 2 ∼ IG(1, 1), δ ∼ IG(1, 1),
ρ ∼ Uni(0, +∞) and α ∼ Uni(0, 1). Posterior distributions are obtained through standard
Markov chain Monte Carlo (MCMC). We ensure that the chain is long enough to converge
and take 10, 000 MCMC samples. Posterior probabilities P (H0i | y) are approximated with
posterior samples.
The second model is specified that the process is independent over time but the variance

9

is not stationary over time:
C(h; t, t′ ) =

(
δt C (s) (h), if t = t′ ;
0,

otherwise.

This is a different specification from Model 1 and neither model is a reduced case of the
other. We use the same specification for C (s) (h) as in Model 1 and also the same prior
distributions for σ 2 , δt and ρ. We take 10, 000 MCMC samples and posterior probabilities
are approximated with MCMC samples.
For both models, we follow the FDR control procedure (4) at a nominal level 0.05. Results
of both are shown in Figure 5. Note that drought impacts Oklahoma mostly in the western
areas. Both results seem reasonable and meaningful for practitioners and they overlap on
most decisions. However, we do observe that, at nine time/location points, they do not agree
with each other. Those nine points are all rejected in Model 1 but neither in Model 2. Table
1 shows the observed values and model inferred upper credible bounds at the given nominal
level for the nine disagreed points. We can see that all nine points are boundary cases and
Model 2 results in higher upper bounds than Model 1, causing the disagreed decisions. Such
disagreed decisions will likely cause confusions in practice. As one must assume normality in
the first place before performing a small sample t test, we believe that, in a good practice,
it is necessary to clearly assume and carefully check the model specification before using
posterior probabilities from the model for testing.

4

Discussion

In this paper, we explore the influence of a misspecified covariance structure on the multiple
testing procedure using Bayesian posterior probabilities. We explicitly show the influence
on the test statistics and discuss the KL divergence as a measure of that influence. We
see from a simulation study that both the correct strength of dependence and structure of
dependence are necessary to ensure control of the FDR at the nominal level. We also see
that misspecified covariance can significantly impact efficiency, in terms of FNR. From a real
application, we see that different covariance specifications can result in different decisions.
This paper does not cover any discussion on a precise (or two-sided) hypothesis: H0i :
θi = 0 versus H1i : θi 6= 0. In that scenario, a mixture model is often assumed: f (y) =
p0 f0 (y) + p1 f1 (y). The local fdr by Efron et al. (2001) is p0 f0 (y)/f (y) under independence.
When data are dependent, it is unclear how to properly incorporate the dependence into
the mixture model. One practical example given by Brown et al. (2014) specifies: yi ∼
10

N(γi µi , σ 2 ), γi ∼ Bern(1 − p) and µ ∼ CAR(ρ, τ 2 ), which mimics the independent model in
Scott and Berger (2006). In the Bayesian framework, to compute P (H0i | y), we would need
P (H0i ), P (H1i ) and priors πi0 (θi ) under H0i and πi1 (θi ) under H1i . Moreover, πi0 (θi ) and
πi1 (θi ) should have a dependence structure for i = 1, . . . , m, in some way. A misspecified
covariance (or model) would be worth further investigation in this setting.

A
A.1

Appendix
Proof of Lemma 1

Proof. First, given the underlying true model (5), it is straightforward to derive the true marginal
distribution for y:
Z
(9)
p(y | θ)p(θ)dθ = N(θ0 , σ02 I + Σ1 ).
p(y) =
θ

If we estimate the posterior using the correct covariance structure, we will have the following
−1
−1
2
posterior distribution: θ | y ∼ N(θ (y) , Σ(y) ), where θ (y) = (1/σ02 I + 1/gΣ−1
1 ) (1/σ0 y + 1/gΣ1 θ0 )
(y)
(y)
−1
and Σ(y) = (1/σ02 I + 1/gΣ−1
1 ) . Marginally, θi | y ∼ N(θi , Σii ). Then,


(y)
θ − θ0i 
.
Hi = P (H0i | y) = P (θi ≥ θ0i | y) = Φ  iq
(y)
Σii

−1 2
Using (9), we have the marginal distribution: θ (y) ∼ N(θ0 , B), where B = (I + σ02 /gΣ−1
1 ) (σ0 I +
(y)
−1 −1
(y)
(y)
2
Σ1 )(I + σ0 /gΣ1 ) . Note that Σ is free of y, so let A = Σ . Marginally, θi ∼ N(θ0i , bii )
which leads to
(y)

θi − θ0i
q
∼ N(0, bii /aii ).
(y)
Σii

Now, under the true covariance, the CDF for Hi is given by
 


(y)
θ − θ0i 
≤ hi 
F (hi ) = P Φ  iq
(y)
Σii


(y)
θ − θ0i
= P  iq
≤ Φ−1 (hi )
(y)
Σii
r

aii −1
= Φ
Φ (hi ) .
bii

11

If a misspecified covariance Σ2 is used to estimate the posterior, then in the posterior distribution
θ | y ∼ N(θ (y) , Σ(y) ), we will have Σ2 , instead of Σ1 in both θ (y) and Σ(y) . As a consequence,
−1 and B = (I + σ 2 /gΣ−1 )−1 (σ 2 I + Σ )(I + σ 2 /gΣ−1 )−1 . The rest remains
A = (1/σ02 I + 1/gΣ−1
1
0
0
0
2 )
2
2
the same.

A.2

Proof of Theorem 1

Proof. According to the proof in Lemma 1, θ (y) ∼ N(θ0 , B). Hence, diag−1/2 (B)(θ (y) − θ0 ) ∼
N(0, Pb ), where Pb is the correlation matrix of B. Equivalently,
(y)

(y)

θ1 − θ01
θm − θ0m
√
,··· , √
b11
bmm

!′

∼ N(0, Pb ).

Then, the joint CDF of (H1 , . . . , Hm )′ is
F (h1 , . . . , hm ) = P (H1 ≤ h1 , . . . , Hm ≤ hm )
#
!
!
"
(y)
(y)
θm − θ0m
θ1 − θ01
≤ h1 , · · · , Φ
≤ hm
= P Φ
√
√
a11
amm
"
#
r
r
(y)
(y)
θ1 − θ01
θm − θ0m
a11 −1
amm −1
√
= P
≤
Φ (h1 ), · · · , √
≤
Φ (hm )
b11
bmm
b11
bmm

r
r
a11 −1
amm −1
m
= Φb
Φ (h1 ), · · · ,
Φ (hm ) .
b11
bmm

A.3

Proof of Theorem 2

Proof. First, given the underlying true process (5), marginally y ∼ N(θ0 , σ02 I + Σ1 ).

If the correct covariance is used, we have the posterior distribution
Z ∞
p(θ | y) ∝
p(y | θ)p(θ | σ 2 )p(σ 2 )dσ 2
0

−m−α
∝ (y − θ)′ (y − θ) + 1/g(θ − θ0 )′ Σ−1 (θ − θ0 ) + 2β
h
i−m−α
∝ (θ − θ (y) )′ (I + 1/gΣ−1 )(θ − θ (y) ) + (y − θ0 )′ (I + gΣ)−1 (y − θ0 ) + 2β
∝



− m+(m+2α)
2
1
(y) ′
(y) −1
(y)
1+
,
(θ − θ ) (V ) (θ − θ )
m + 2α

12


−1
which is tm (m + 2α, θ (y) , V (y) ), with location θ (y) = I + 1/gΣ−1
y + 1/gΣ−1 θ0 and scale


−1
V (y) = (m + 2α)−1 2β + (y − θ0 )′ (I + gΣ)−1 (y − θ0 ) I + 1/gΣ−1
. Similarly as Lemma 1,
−1 and B = (I + 1/gΣ−1 )−1 (σ 2 I + Σ )(I + 1/gΣ−1 )−1 .
define A = (I + 1/gΣ−1
)
1
0
1
1
1
(y)

(y)

Marginally θi | y ∼ t1 (m + 2α, θi , Vii ), then each test statistic is


(y)
θ
−
θ
i0 
Hi = P (H0i | y) = P (θi ≥ θi0 | y) = Ψm+2α  iq
.
(y)
Vii

q

(y)
((θ1 − θ10 )/

(y)

(y)

V11 , . . . , (θm −
In order to find the joint CDF, we need the joint distribution of
q
(y)
θm0 )/ Vmm )′ , which can be re-written as
s
− 1

1
m + 2α
2
(θ (y) − θ0 ) =
diag V (y)
(diag A)− 2 (θ (y) − θ0 ), (10)
′
−1
(y − θ0 ) (I + gΣ1 ) (y − θ0 ) + 2β
Given the marginal distribution of y, we have θ (y) − θ0 = A(y − θ0 ) ∼ N(θ0 , B), or, θ (y) − θ0 =
(diag B)1/2 zb . The quadratic term in (10) is
(y − θ0 )′ (I + gΣ1 )−1 (y − θ0 ) = (y − θ0 )′ (I − A)(y − θ0 )

= zb′ (diag B)−1/2 (A−2 − A−1 )(diag B)−1/2 zb .

The equation (10) is then
− 1

2
(y)
(θ (y) − θ0 ) =
diag V

s

1
1
m + 2α
(diag A)− 2 (diag B) 2 zb .
zb′ Czb + 2β

Or, equivalently,

′ s
r
r
(y)
(y)
m + 2α
 a11 θ1q − θ10 , . . . , amm θmq − θm0  =
′ Cz + 2β zb .
b11
bmm
z
(y)
(y)
b
b
V11
Vmm
The joint CDF of H1 , . . . , Hm is given by

F (h1 , . . . , hm ) = P (H1 ≤ h1 , . . . , Hm ≤ hm )






(y)
(y)
θm − θm0 
θ − θ10 
≤ h1 , . . . , Ψm+2α  q
= P Ψm+2α  1q
≤ hm 
(y)
(y)
V11
Vmm


(y)
(y)
θm − θm0
θ − θ10

≤ Ψ−1
= P  1q
≤ Ψ−1
m+2α (h1 ), . . . , q
m+2α (hm )
(y)
(y)
V11
Vmm

r
r
a11 −1
amm −1
m
Ψ
(h1 ), · · · ,
Ψ
(hm ) .
= Ξa,b
b11 m+2α
bmm m+2α

13

If the misspecified covariance is used, follow the same argument in Lemma 1.

A.4

Proof of Corollary 1

√
2
Proof. Notice that dΦ−1 (x)/dx = 1/ϕ(Φ−1 (x)), where ϕ(x) = 1/ 2πe−x /2 . Then
√

ri
√ −1
f (hi ) = ϕ ri Φ (hi )
−1
ϕ(Φ (hi ))


√
1
(1 − ri )[Φ−1 (hi )]2 .
=
ri exp
2

A.5

Proof of Corollary 2

Proof. Let ui =

√ −1
ri Φ (hi ). Using the chain rule,
m

f (h1 , . . . , hm ) =
=
=
=
=

Y ∂ui
∂ m Φm
∂mF
b
=
∂h1 · · · ∂hm
∂u1 · · · ∂um
∂hi
i=1
√
m
Y
ri
ϕm
(u
,
.
.
.
,
u
)
1
m
b
−1
ϕ(Φ (hi ))
i=1


Y

m
√
1 ′ −1
1 −1
− 12
2
|Pb | exp − u Pb u
ri exp
[Φ (hi )]
2
2
i=1
)
(m


Y√
1
1
1
1
1
ri |Pb |− 2 exp − φ′ R 2 Pb−1 R 2 φ + φ′ φ
2
2
i=1


1
1
1
1 2
1
1 ′
R 2 Pb−1 R 2 exp
φ (I − R 2 Pb−1 R 2 )φ .
2

Acknowledgement
This research is partially supported by National Science Foundation under Grant No. OIA1301789.

14

References
Benjamini, Y. and Hochberg, Y. (1995), ‘Controlling the false discovery rate: a practical and
powerful approach to multiple testing’, Journal of the Royal Statistical Society, Series B
57, 289–300.
Benjamini, Y. and Yekutieli (2001), ‘The control of the false discovery rate in multiple testing
under dependency’, Annals of Statistics 29, 1165–1188.
Brown, A., Lazar, N., Datta, G., Jang, W. and McDowell, J. (2014), ‘Incorporating spatial
dependence into Bayesian multiple testing of statistical parametric maps in functional
neuroimaging’, Neuroimage 84, 97–112.
Efron, B. (2007), ‘Correlation and large-scale simultaneous significance testing’, Journal of
the American Statistical Association 102, 93–103.
Efron, B. (2012), Large-scale inference: Empirical Bayes methods for estimation, testing and
prediction, Cambridge University Press.
Efron, B., Tibshirani, R., Storey, J. and Tusher, V. (2001), ‘Empirical bayes analysis of a
microarray experiment’, Journal of the American Statistical Association 96, 1151–1160.
Farcomeni, A. (2007), ‘Some results on the control of the false discovery rate under dependence’, Scandinavian Journal of Statistics 34, 275–297.
Illston, B., Basara, J., Fisher, D., Elliott, R., Fiebrich, C., Crawford, K., Humes, K. and
Hunt, E. (2008), ‘Mesoscale monitoring of soil moisture across a statewide network’, Journal of Atmospheric and Oceanic Technology 25, 167–182.
Liang, F., Paulo, R., Molina, G., Clyde, M. and Berger, J. (2008), ‘Mixtures of g priors for
Bayesian variable selection’, Journal of the American Statistical Association 103, 410–423.
Liang, K. and Nettleton (2012), ‘Adaptive and dynamic adaptive procedures for false discovery rate control and estimation’, Journal of the Royal Statistical Society, Series B
74, 163–182.
Muller, P., Parmigiani, G. and Rice, K. (2006), FDR and Bayesian multiple comparisons
rules, in ‘Proceedings of Valencia / ISBA 8th World Meeting on Bayesian Statistics’,
Benidorm (Alicante, Spain).

15

Scott, J. and Berger, J. (2006), ‘An exploration of aspects of Bayesian multiple testing’,
Journal of Statistical Planning and Inference 136, 2144–2162.
Scott, J. and Berger, J. (2010), ‘Bayes and empirical-Bayes multiplicity adjustment in the
variable-selection problem’, The Annals of Statistics 38, 2587–2619.
Storey, J., Taylor, J. and Siegmund, D. (2004), ‘Strong control, conservative point estimation
and simultaneous conservative consistency of false discovery rates: a unified approach’,
Journal of the Royal Statistical Society, Series B 66, 187–205.
Sun, W. and Cai, T. (2009), ‘Large-scale multiple testing under dependence’, Journal of the
Royal Statistical Society, Series B 71, 393–424.
Sun, W., Reich, B., Cai, T., Guindani, M. and Schwartzman, A. (2015), ‘False discovery
control in large-scale spatial multiple testing’, Journal of the Royal Statistical Society,
Series B 77, 59–83.
Wu, W. (2008), ‘On false discovery control under dependence’, Annals of Statistics 36, 364–
380.

16

Table 1: Soil Data: Disagreed decisions under Model 1 and 2, for the hypothesis H0i : θi ≥ 50
versus H0i : θi < 50. † denotes a rejection under the model. The nominal level is 0.05.
Time & Location Observed Value
January, Site 6
49.49
January, Site 70
49.55
49.62
January, Site 72
January, Site 87
49.70
49.62
January, Site 101
March, Site 8
49.94
March, Site 47
49.58
49.53
April, Site 38
June, Site 53
49.63

Upper C.I. under Model 1 Upper C.I. under Model 2
50.42†
51.03
50.52†
51.30
50.37†
50.89
50.67†
51.07
50.70†
51.21
50.78†
51.10
50.36†
51.22
50.19†
50.90
50.85†
51.83

17

5

1.0

4

a/b=0.2
a/b=1
a/b=5

0

0.0

1

0.2

2

0.4

f(h)

F(h)

3

0.6

0.8

a/b=0.2
a/b=1
a/b=5

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

h

0.4

0.6
h

Figure 1: Marginal CDF and pdf for hi .

18

0.8

1.0

0.6

0.10

Correct
Misspecified

FNR
0.3
0.0

0.02

0.1

0.04

0.2

FDR
0.06

0.4

0.08

0.5

Correct
Misspecified

0.0

0.5

1.0
1.5
log10 (g)

2.0

2.5

−0.5

0.0

0.5

1.0
1.5
log10 (g)

2.0

2.5

−0.5

0.0

0.5

1.0
1.5
log10 (g)

2.0

2.5

−0.5

0.0

0.5

1.0
1.5
log10 (g)

2.0

2.5

0.4

DKL m

0.0

0.00

0.1

0.05

0.2

0.3

Diff Rej Rate
0.10
0.15

0.5

0.20

0.6

0.25

0.7

−0.5

Figure 2: Results for Example 1. FDR, FNR, Difference of rejection rate (ratecor − ratemis )
and DKL /m. The sequence of g is from 0.2 to 500. The nominal level is 0.05.

19

0.6

0.10

Correct
Misspecified

FNR
0.3
0.0

0.02

0.1

0.04

0.2

FDR
0.06

0.4

0.08

0.5

Correct
Misspecified

0.0

0.5

1.0
1.5
log10 (g)

2.0

2.5

−0.5

0.0

0.5

1.0
1.5
log10 (g)

2.0

2.5

−0.5

0.0

0.5

1.0
1.5
log10 (g)

2.0

2.5

−0.5

0.0

0.5

1.0
1.5
log10 (g)

2.0

2.5

0.0

0.00

0.05

0.5

DKL m

Diff Rej Rate
0.10
0.15

1.0

0.20

0.25

1.5

−0.5

Figure 3: Results for Example 2. FDR, FNR, Difference of rejection rate (ratecor − ratemis )
and DKL /m. The sequence of g is from 0.2 to 500. The nominal level is 0.05.

20

0.7
0.4

DKL m

0.18

0.1

0.2

0.20

0.3

FNR
0.22 0.24

0.5

0.26

0.6

0.28

0.09
0.08
FDR
0.07
0.06

0.0

0.05
−2

−1

0
1
log ρ

2

3

−2

−1

0
1
log ρ

2

3

−2

−1

0
1
log ρ

2

3

Figure 4: Results for Example 3. FDR, FNR and DKL /m. The dashed vertical line is the
correct value ρ = 5. The sequence of ρ is (0.1, 0.5, 1, 5, 10, 20). The nominal level is 0.05.

21

Jan

Feb

Mar

Apr

May

Jun

Jul

Aug

Sep

Oct

Nov

Dec

Figure 5: Rejected sites for the year of 2014 in Oklahoma, from January to December. Black
(·) are the observations, red (+) are rejections using Model 1, and green (△) are rejections
using Model 2.

22

