A fully objective Bayesian approach for the
Behrens-Fisher problem using historical studies
Antoine Barbieri1,2

Jean-Michel Marin1

Karine Florin3

(Antoine.Barbieri@umontpellier.fr)

(Jean-Michel.Marin@umontpellier.fr)

(Karine.Florin@sanofi.fr)

arXiv:1611.06873v1 [] 21 Nov 2016

1
2

Institut Montpelliérain Alexander Grothendieck (IMAG), Université de Montpellier
Institut régional du Cancer Montpellier - Val d’Aurelle (ICM), Unité de Biométrie
3 Unité de Biostatistiques, Sanofi, Montpellier, France

Abstract
For in vivo research experiments with small sample sizes and available historical data, we propose
a sequential Bayesian method for the Behrens-Fisher problem. We consider it as a model choice
question with two models in competition: one for which the two expectations are equal and one
for which they are different. The choice between the two models is performed through a Bayesian
analysis, based on a robust choice of combined objective and subjective priors, set on the parameters
space and on the models space. Three steps are necessary to evaluate the posterior probability of each
model using two historical datasets similar to the one of interest. Starting from the Jeffreys prior,
a posterior using a first historical dataset is deduced and allows to calibrate the Normal-Gamma
informative priors for the second historical dataset analysis, in addition to a uniform prior on the
model space. From this second step, a new posterior on the parameter space and the models space
can be used as the objective informative prior for the last Bayesian analysis. Bayesian and frequentist
methods have been compared on simulated and real data. In accordance with FDA recommendations,
control of type I and type II error rates has been evaluated. The proposed method controls them even
if the historical experiments are not completely similar to the one of interest.

Keywords: Behrens-Fisher problem; objective Bayesian approach; model choice; small sample sizes.

1

Introduction

In preclinical research, experiments are routinely performed using exactly the same protocol under the
same conditions. Historical data are often available with the same control, the same reference product or
similar products to the one tested. A specific feature of in vivo pharmacology experiments is the small
sample size per experiment, due to ethical considerations (Landis et al., 2012). Frequently, the objective
of such experiments is to evaluate a treatment effect in comparison to a control (placebo or reference
compound).
Before dealing with the experimental context used to evaluate the effects of a compound at different doses or to compare different treatments, we focused on a simple study design of an experimental
treatment group compared to a control group in parallel. The measurements are supposed to be normally distributed with expectation µt when the experimental treatment is administered and µc if not.

1

The variances of the two Gaussian distributions are unknown and assumed to be different. The aim is to
discriminate between the two hypotheses

H0 : µc = µt = µ
(1)
H1 : µc 6= µt
,
which is known as the Behrens-fisher problem (Kim and Cohen, 1998).
In the case of hypothesis testing, the frequentist methods are routinely applied using the p-value
interpretation. These methods do not take into account the historical data and are poorly adapted for small
sample sizes. The frequentist approach for solving this problem is the Student test with Satterthwaite
correction. However, the small sample size causes a lack of information in the statistical analysis, hence
the idea of including some prior information in the inferential process, to increase the discriminative
power of the procedure. With the frequentist method, a crude way to take into account historical datasets
is to pool the experiments. One can test if the treatment effect is significant after pooling the data. In
contrast, the Bayesian paradigm is considered to be a well-suited method for handling small samples, in
particular since it allows taking into account prior information. In the Bayesian paradigm, two directions
can be adopted:
• the first one is to use the historical datasets to determine informative prior distributions;
• the second one is to use a hierarchical model which assumes a distribution across studies with a
parameter that controls the variation.
There are numerous studies on the second case (Neuenschwander et al., 2010; Hobbs et al., 2011; Viele
et al., 2014; Schmidli et al., 2014), which is an interesting and natural modeling of the problem. However,
the results are very sensitive to the prior distribution on the between-experiments dispersion parameter
(Gelman, 2006). In this work, we investigated the first case with constraint on the sample size. In this
Bayesian situation, work has previously been performed on sample size calculation (Whitehead et al.,
2008, 2015), but that won’t be addressed in our work.
For a given prior distribution, a typical procedure to discriminate between H0 and H1 is to construct
a Bayesian credible interval of level 1 − α (α typically equal to 5%) on µc − µt and to decide in favour
of H0 if 0 belongs to the credible interval and H1 if not (Kim and Cohen, 1998; Ghosh and Kim, 2001).
This is not the fully Bayesian approach. In the Bayesian paradigm, a prior distribution on the space
of hypotheses must be defined and then their posterior distribution calculated (Marin and Robert, 2013;
Kass and Raftery, 1995). Indeed, each hypothesis is associated with a model: the model M0 under
which µc = µt and the model M1 under which µc 6= µt . Then, prior probability of model M0 , denoted
Pr (M0 ), directly related to M1 (with Pr (M1 ) = 1 − Pr (M0 )) is set and the posterior probability of
M0 is calculated. Using this last probability and a loss function, the experimenter can decide whether to
reject or not H0 .
There are very few works on a fully Bayesian approach for the Behrens-Fisher problem. To our
knowledge, there is only one proposal of Moreno et al. (1999) in which intrinsic and fractional Bayes
factors are calculated to avoid the improper prior difficulty. Here, we propose instead to use some
historical data. But, it is not easy to calibrate an objective prior distribution using historical datasets such
that the model posterior probabilities are correctly defined. In this work, we present a solution in three
steps involving two similar historical datasets used sequentially to calibrate the prior. In this work, we
propose a solution in three steps involving two similar historical datasets used sequentially to calibrate
2

the prior. The first allows building an informative prior on the model parameters from an improper and
non informative prior. The informative priors on the model and parameter spaces are calibrated at the
completion of the second step using the proper informative prior on parameter from the first step and the
second considered historical dataset. The last step returns the posterior probability of each model, given
the informative prior from the second step, allowing to make the statistical inference on the experiment of
interest. The plan of the paper is as follows. In section 2, the methodology is introduced and the proposed
Bayesian method is detailed. Section 3 then presents a simulation study to compare our method and the
frequentist methods. Finally, section 4 applies both approaches on real data from in vivo behavioral
pharmacology experiments.

2
2.1

Methods
Recap on Bayesian model choice

Let y denote the observed dataset and J the number of models in competition. In the parametric Bayesian
paradigm, model j ∈ {0, . . . , J − 1}, denoted Mj , is defined with
• a likelihood function `j (θj |y) with unknown parameter θj ∈ Θj ;
• a prior distribution denoted by πj (θj ) on the parameter space of Mj .
The Bayesian model choice is done according to the model posterior probabilities (Robert, 2007; Scott
and Berger, 2005). Given prior probabilities on the models Pr (Mj ), the posterior distribution on the
model space is deduced by
(
)
Z
Pr (Mj |y) ∝

`j (θj |y)πj (θj )dθj

Pr (Mj )

.

Θj

With Pr(y|Mj ) =

R

Θj `j (θj |y)πj (θj )dθj ,

the posterior probabilities become

Pr (Mi ) Pr (y | Mi )
Pr (Mi | y) = PJ−1
.
j=0 Pr (Mj ) Pr (y | Mj )

(2)

The key quantities in equation (2) are the Pr (y | Mj ) which are called marginal or integrated likelihoods. In absence of a loss function, we select the model one with the largest posterior probability.

2.2

The model

A set of two independent and identically distributed (iid) samples c and t are considered. The first one
associated with the control group is assumed to be normally distributed of size nc with expectation µc and
variance σc2 . The second one of size nt is associated with the treated group and assumed to be normally
distributed with expectation µt and variance σt2 . We focused on the Behrens-Fisher problem where the
aim is to test (1). In the fully Bayesian paradigm presented in the previous section, this hypothesis testing
is equivalent to a model choice problem between the two following models:
2 ) and t ∼
2
• M0 under which c ∼iid N (µ, σ0,c
iid N (µ, σ0,t ),

3

2 ) and t ∼
2
• M1 under which c ∼iid N (µc , σ1,c
iid N (µt , σ1,t ).
2 , σ 2 ), θ = (µ , µ , σ 2 , σ 2 ) and y = (c, t). According to model j, the posterior
Let θ0 = (µ, σ0,c
1
c t 1,c 1,t
0,t
distribution of θj is:
πj (θj |y) ∝ `j (θj |y) π (θj ) .

The aim of our work is to introduce a methodology that uses historical datasets to calibrate objective
prior distributions on θ0 and θ1 and also on the model index. We assume to have at disposal two previous
datasets from similar experiments.

2.3

Our objective Bayesian approach

Typically, to solve the Behrens-Fisher problem, the practitioners use the following Bayesian approach:
• to forget the model choice aspect of the problem;
• to construct a credible interval on the difference between the two expectation;
• to reject H0 if 0 does not belong to the credible interval.
In this case, the Bayesian paradigm is just used to construct a credible interval. The prior distribution on
the two hypotheses, i.e. the prior probabilities Pr(M0 ) and Pr(M1 ), do not exist. In the contrast, a fully
Bayesian approach considers prior probabilities and deduces posterior probabilities on both hypotheses.
The benefit of using such a procedure is to weight each hypothesis according to some knowledge such
as the data. Using these weights, the experimenter can then decide whether to reject H0 or not. The
setting of the prior distribution can be controversial, especially if it is based on expert opinion. However,
we show below how to calibrate the prior in presence of historical data. Thus, a robust choice taking
advantage of historical data to build an objective prior is proposed: starting from a non informative prior,
the historical data are used to deduce an informative proper prior.
Priors must be defined for θ0 , θ1 and also for the probabilities Pr(M0 ) and Pr(M1 ). Our proposal
is to use a posterior distribution based on the historical data as a prior for the experiment. An initial
prior distribution is needed and we propose to rely on Jeffreys prior (Box and Tiao, 1973). This prior
maximizes the distance between the prior and the posterior and is the standard non informative proposal,
on the two parameters θ0 and θ1 . It is defined as:
v
"
#
u
q
2 log ` (θ | y)
u
∂
j
j
t
πj (θj ) ∝ |I F (θj )| =
−Eθj
, j = 0, 1
(3)
∂θjt ∂θj
where I F (θj ) is Fisher information matrix and is interpreted as the amount of information provided by
the observation y on θj . Unfortunately, Jeffreys prior distributions on θ0 and θ1 are improper, since the
integrals are infinity. The corresponding posteriors are well-defined but, when an improper prior is used
on the parameters, the posterior in the model space is not properly defined (Robert, 2007; Gelman et al.,
2013). Indeed, the posterior probabilities of both models depend on the unknown normalizing constant.
To compute this posterior, we propose to proceed in three steps using two historical datasets similar to
the experiment we want to analyse. The experiment used in the first step is referred to as the experiment
1 or the first experiment. The experiment used in the second step is referred to as experiment 2 or the
second experiment. Finally, the third experiment (also named the experiment 3 or the experiment of
interest) used in the last step is the dataset on which the inference is made. The three steps are:
4

1. A first set of historical data, together with the Jeffreys priors on θ0 and θ1 are considered to deduce
a first proper posterior on θ0 and θ1 ;
2. A second set of historical data is considered with the use the first posterior as prior. We introduce
of a Laplace (flat) prior on the two hypotheses, namely Pr(M0 ) = Pr(M1 ) = 1/2, and deduce a
new posterior on θ0 and θ1 and also on the hypotheses;
3. In the last and third step of our procedure, we deduce the posterior probabilities of the models M0
and M1 which are computed using the posterior of the second step as an objective informative
prior.
In the next two subsections, the marginal likelihoods for the three step procedure are detailed. The prior
and the posterior distributions for each step of our proposal are indexed. Let us denote by:
• y1 = (c1 , t1 ) the first historical set where c1 is an iid sample of size nc1 and t1 an iid sample of
size nt1 ;
• y2 = (c2 , t2 ) the second one (of sizes nc2 and nt2 );
• y3 = (c3 , t3 ) the data to analyse (of sizes nc3 and nt3 );
• ci and ti denote the empirical expectations of ci and ti , for i = 1, 2, 3 and γci =
2
Pnti
ti − ti .
and γti = i=1

2.4

Pnci

i=1 (ci

− ci )2

Computation of integrated likelihoods under model M1

Let us first consider the case of model M1 for which all the computations are explicit. Under model
M1 , the likelihood of θ1 for dataset yi is:
)
(
i
 −nci
−1 h
2
2
2
exp
nci (ci − µc ) + γci
×
`1 (θ1 |yi ) = 2πσ1,c
2
2σ1,c
(
)
h
i
 −nti
2
−1
2
2πσ1,t 2 exp
nti ti − µt + γti
(4)
2
2σ1,t
Regarding the first step of the analysis, the Jeffreys’ prior on θ1 used is:
π11 (θ1 ) ∝ σc2 σt2

− 3

2

.

The posterior distribution of θ1 can be then deduced. Indeed,
π11 (θ1 |y1 ) ∝ `1 (θ1 | y1 ) π11 (θ1 )




2
2
2
2
|y1
∝ π11 µc |σ1,c
, y1 π11 σ1,c
|y1 π11 µt |σ1,t
, y1 π11 σ1,t

(5)

Equation (5) shows posterior independence between parameters of the two groups of the first step. The
conditional posterior distributions of parameters are deduced such as the location parameters µc and µt

5

2 and σ 2 belongs to the inverse-gamma
are Gaussian (N ) and the marginal posterior distributions of σ1,c
1,t
family (IG). Thus, the prior distributions that will be used at the second step of the analysis are:
!
!
2
2
σ
σ
1,c
1,t
2
2
, µt |y1 , σ1,t
∼ N t1 ,
,
µc |y1 , σ1,c
∼ N c1 ,
n c1
nt1
2
σ1,c
∼ IG

Using

π12 (θ1 )

=

π11 (θ1 |y1 ),

n

c1

2

,

γ c1 
,
2

`1 (θ1 |y2 )π12 (θ1 )dθ1

Z

`1 (θ1 |y2 )π11 (θ1 |y1 )dθ1

=

= (2π)



n

t1

2

,

γt1 
.
2

the integrated likelihood of y2 becomes:

Z
Pr(y2 |M1 ) =

2
σ1,t
∼ IG

nc2 +nt
2
2

γt1 +γt2
2

+



 1 γt1  n2t1 γc1  n2c1
2
nc1 nt1
2
2
n 
n  ×
(nc1 + nc2 ) (nt1 + nt2 )
Γ 2t1 Γ 2c1
 


n +n
n +n
Γ t1 2 t2 Γ c1 2 c2

nt1 nt2 (t1 −t2 )

2

 nt1 +nt2 
2

2(nt1 +nt2 )

γc1 +γc2
2

+

nc1 nc2 (c1 −c2 )
2(nc1 +nc2 )

2

 nc1 +nc2

.

(6)

2

Similarly, using π12 (θ1 ) = π11 (θ1 |y1 ), the posterior distribution obtained in the second step is:
π12 (θ1 |y2 ) ∝ `1 (θ1 | y2 ) π12 (θ1 )




2
2
2
2
∝ π12 µc |σ1,c
, y2 π12 σ1,c
|y2 π12 µt |σ1,t
, y2 π12 σ1,t
|y2 .

(7)

Like for the first step, equation (7) shows posterior independence between parameters of the two groups
for the second step. Used as prior for the third step, the posterior distribution of each parameter is thus
deduced for the second step and corresponds to:
!
!
2
2
σ1,c
σ1,t
nc1 c1 + nc2 c2
nt1 t1 + nt2 t2
2
2
,
, µt |y2 , σ1,t ∼ N
,
,
µc |y2 , σ1,c ∼ N
nc1 + nc2
nc1 + nc2
nt1 + nt2
nt1 + nt2
!
2
γ
+
γ
n
n
(c
−
c
)
n
+
n
c
c
c
c
1
2
c
c
2
1
2
2
, 1
+ 1 2
,
σ1,c
∼ IG
2
2
2 (nc1 + nc2 )
2 !
n
n
t
−
t
n
+
n
γ
+
γ
t
t
1
2
t
t
t
t
1
2
1
2
2
2
σ1,t
∼ IG
, 1
+
.
2
2
2 (nt1 + nt2 )

6

Thus, using π13 (θ1 ) = π12 (θ1 |y2 ), the integrated likelihood of y3 becomes:
Z
Pr(y3 |M1 ) =
`1 (θ1 |y3 )π13 (θ1 )dθ1
Z
=
`1 (θ1 |y3 )π12 (θ1 |y2 )dθ1

= (2π)


nc3 +nt

3

2

γc1 +γc2
2

×



+

(nc1 + nc2 ) (nt1 + nt2 )
(nc1 + nc2 + nc3 ) (nt1 + nt2 + nt3 )
nc1 nc2 (c1 −c2 )2
2(nc1 +nc2 )

Γ



nc1 +nc2
2

 nc1 +nc2
2

Γ



1

2

2

+

nt1 nt2 (t1 −t2 )

Γ

nt1 +nt2 +nt3
2

βt

 nt1 +nt2

nt1 +nt2
2

nc1 +nc2 +nc3
2



,

(8)

βc

where
2
2

nt1 nt2 t1 − t2
γt1 + γt2 + γt3
(nt1 + nt2 ) nt3
nt1 t1 + nt2 t2
+
βt =
+
− t3 ,
2
2 (nt1 + nt2 )
2 (nt1 + nt2 + nt3 )
nt1 + nt2
γc + γc2 + γc3
nc nc (c1 − c2 )2
(nc1 + nc2 ) nc3
βc = 1
+ 1 2
+
2
2 (nc1 + nc2 )
2 (nc1 + nc2 + nc3 )

2.5



nc1 c1 + nc2 c2
− c3
n c1 + n c2

2
.

Calculation of integrated likelihoods under model M0

Under model M0 , the likelihood of θ0 for dataset yi is:
`0 (θ0 |yi ) =

 −nci
 −nti
2
2
2
2
2πσ0,c
2πσ0,t

!2 
2 c + n σ2 t
 n σ2 + n σ2 

n
σ
ti 0,c i
ci 0,t
ti 0,c
ci 0,t i
× exp −
µ
−
2 σ2
2 + n σ2


2σ0,c
nci σ0,t
ti 0,c
0,t

2 

 γ
c
−
t
n
n
γ
i
i
c t
t
c
 .
× exp − 2i − 2i −  i i
 2σ0,t 2σ0,c 2 n σ 2 + n σ 2 
c
t
i

0,t

i

0,c

The Jeffreys prior on θ0 used in the first step of the analysis is:
π01 (θ0 ) ∝

2 + n σ2
nc σ0,t
t 0,c
6 σ6
σ0,c
0,t

!1
2

.

For each step in the analysis, the posterior distribution of µ is a Gaussian distribution with expectations
2 and σ 2 , the posterior
and variances described in Table 1. However, due to the dependence between σ0,c
0,t

7

2

2(nt1 +nt2 )



 

n +n +n
Γ c1 2c2 c3

nt +nt +nt
1
2
3
2



γt1 +γt2
2


2 , σ2
distribution of the couple σ0,c
0,t is not explicit. For instance, for the first step of the analysis, the
2 , σ 2 ) is:
posterior density of (σ0,t
0,t

2 

 γ
−nt
−nc1
1



n
n
c
−
t
γ
−1
−1
c t
1
1
t
c
2
2
2
2
2
2
 .
σ0,t
exp − 21 − 21 −  1 1
π01 σ0,c
, σ0,t
|y1 ∝ σ0,c
 2σ0,c 2σ0,t 2 n σ 2 + n σ 2 
c1 0,t

t1 0,c

(9)
2
2
Table 1: Under model M0 , posterior expectation and variance of µ conditional on σ0,c
and σ0,t
for the three

steps of the analysis.
2 , σ2 , y
µ|σ0,c
0,t i

 2

2 ,y
E µ|σ0,c
, σ0,t
i
2 c + n σ2 t
nc1 σ0,t
1
t1 0,c 1

 2

2 ,y
V µ|σ0,c
, σ0,t
i
2 σ2
σ0,c
0,t

2 + n σ2
nc1 σ0,t
t1 0,c

2 (n c + n c ) + σ 2
σ0,t
c1 1
c2 2
0,c nt1 t1 + nt2 t2

2 + n σ2
nc1 σ0,t
t1 0,c
2 σ2
σ0,c
0,t

2 (n + n ) + σ 2 (n + n )
σ0,t
c1
c2
t1
t2
0,c

2 (n c + n c + n c ) + σ 2
σ0,t
c1 1
c2 2
c3 3
0,c nt1 t1 + nt2 t2 + nt3 t3

2 (n + n ) + σ 2 (n + n )
σ0,t
c1
c2
t1
t2
0,c
2 σ2
σ0,c
0,t

2 (n + n + n ) + σ 2 (n + n + n )
σ0,t
c1
c2
c3
t1
t2
c2
0,c

2 (n + n + n ) + σ 2 (n + n + n )
σ0,t
c1
c2
c3
t1
t2
c2
0,c

Step 1
(conditioned on y1 , first historical data)

Step 2
(conditioned on y2 , second historical data)

Step 3
(conditioned on y3 , interested experiment)

A Markov Chain Monte Carlo (MCMC) algorithm through the software WinBUGS (Cowles, 2004;
Kery, 2010) is used to obtain a posterior sample from (9). Then, as prior for the second and third steps
of the analysis, a product of Inverse-Gamma distributions with parameters calibrated using the posterior
sample given by WinBUGS is considered. According to the same process for the steps of the analysis
i = 1, 2, this corresponds to:





i σ 2 |y = π i σ 2 |b
bc ∼ IG α
bc
π
α
,
β
b
,
β

i
c
c
0
0,c
0
MCMC methods
2
2
 0,c i i 
 i i ,
π0i σ0,c
, σ0,t
|yi −−−−−−−−−−−−−→

i
2
2
2
b
Independence
hypothesis
|
{z
}
π0 σ0,t |yi = π0 σ0,t |b
αti , βti ∼ IG α
bti , βbti
unknown

where the parameters α
b and βb are estimated by the maximum likelihood method on the MCMC paths.
2 , σ 2 ), target of the
Note that, for the second step of theanalysis, we have the posterior density of (σ0,c
0,t
second MCMC scheme, such as:
π02


2
2
σ0,c
, σ0,t
|y2

∝

)
 12
2
2
σ0,t
nc1 + σ0,c
nt1
βbt1
exp − 2
2
2
σ0,t
σ0,t
(nc1 + nc2 ) + σ0,c
(nt1 + nt2 )
(


 )

2
2
2
2
2
2

 2
nc1 σ0,t + nt1 σ0,c nc2 σ0,t + nt2 σ0,c
nc2 σ0,t c2 + nt2 σ0,c t2
2
2

× exp − 2 2  2
− E µ|σ0,c , σ0,t , y1
2
2
2
nc2 σ0,t
+ nt2 σ0,c
2σ0,c σ0,t σ0,t (nc1 + nc2 ) + σ0,c
(nt1 + nt2 )
(
2 )
nc2 nt2 c2 − t2
γc2
γt2
 .
× exp − 2 − 2 −
2
2
2σ0,c
2σ0,t
2 nc2 σ0,t
+ nt2 σ0,c
2 
σ0,c

−nc2
2

−α
b c1 −1

2 
σ0,t

−nt
1
2

−α
b t1 −1

(

βbc
exp − 2 1
σ0,c

)

(

For steps 2 and 3, the prior distributions are:
 

 
π02 (θ0 ) = π01 µ|σc2 , σt2 , y1 π01 σc2 |b
αc1 , βbc1 π01 σt2 |b
αt1 , βbt1 ,
 

 
αc2 , βbc2 π02 σt2 |b
αt2 , βbt2 .
π03 (θ0 ) = π02 µ|σc2 , σt2 , y2 π02 σc2 |b
8

(10)

Thus, the integrated likelihoods are:
Z
Pr (y2 |M0 ) =
`0 (θ0 | y2 ) π02 (θ0 ) dθ0
Z
=
`0 (θ0 | y2 ) π01 (θ0 |y1 ) dθ0
Z
 

  2
2
2
2
2
2
dσ0,t
|b
αt1 , βbt1 dµdσ0,c
=
`0 (θ0 | y2 ) π01 µ|σ0,c
, σ0,t
, y1 π01 σ0,c
|b
αc1 , βbc1 π01 σ0,t

Z

`0 (θ0 | y3 ) π03 (θ0 ) dθ0

Z

`0 (θ0 | y3 ) π02 (θ0 |y2 ) dθ0

Z


 
  2
2
2
2
2
2
dσ0,t
|b
αt1 , βbt1 dµdσ0,c
`0 (θ0 | y3 ) π02 µ|σ0,c
, σ0,t
, y2 π02 σ0,c
|b
αc1 , βbc1 π02 σ0,t

Pr (y3 |M0 ) =
=
=

The integrations over µ are explicit, in contrast to the variance parameters. Numerical integration is thus
used to solve this difficulty through the R package cubature and the adaptInegrate function which implements an adaptive quadrature method (code by Steven G. Johnson and by Balasubramanian Narasimhan,
2013).

3

Simulation study

In this section, the relative efficiency of the proposed method is studied and compared with classical approaches concerning the Behrens-Fisher problem using simulation. First, we assess the frequentist properties of our Bayesian approach as recommended by the FDA (Food and Drug Administration, 2010).
Our methodology is compared with two currently used approaches: Student test with Satterthwaite correction applied only on the experiment of interest and Student test with Satterthwaite correction applied
on the three experiments pooled. Using the frequestist approaches, the null hypothesis is rejected if the
p-value is less than 5%. Then, our method’s behavior is also assessed in the similar cases when a weight
is affected to historical data.
The output of our Bayesian method is the posterior probability P r (M1 |y) which is directly interpretable and more intuitive than a p-value. Given this probability, the experimenter has to discriminate
between the two hypotheses. Thus, a decision rule has to be defined according to a chosen threshold.
In our case, we accept that the expectations are different if the posterior probability is greater than a
certain threshold p such as P r (M1 |y) > p. The first intuitive threshold is p = 0.5. Indeed, in the
case of two models, if the posterior probability of one model is greater than 0.5, this model is more
probable than the other one, given the data. We also propose another arbitrary threshold which is more
conservative for M0 : p = 0.8. If the posterior probability of the model is greater than 0.8, we assume
the probability is enough to decide with a minimal risk of concluding wrongly to a difference. Finally,
a threshold (tailor-made) was also calibrated to control the type I error rate given the studied situations.
From the simulations under the model M0 , the latter threshold corresponded to the posterior probability
P r (M1 |y) associated with the 5% rank.

9

We simulated three datasets (experiments) using Gaussian distributions in the context of the Complete Freund’s Adjuvant (CFA) protocol, detailed in the real data analysis section. The values of the
parameters in the simulations are chosen according to prior knowledge gained from experimental results
in the CFA protocol. The small sample sizes (nci )i=1,2,3 = 10 and (nti )i=1,2,3 = 10 are considered.
From all the simulated data, we defined µc = 2.94 and µt = µc + µc δ where δ is the percentage of
effect. Concerning the standard deviation of the control group, (σci )i=1,2,3 = 0.6 whatever the situation.
Regarding the treated group standard deviations (σti )i=1,2,3 , we considered the different values 0.6, 1.5
and 3, respectively associated with the minimum, median and maximum values deduced from historical
data. For all simulation studies, each experiment is simulated a thousand times (N = 1000).

3.1

Comparison between the Bayesian and frequentist approaches

The first situation presented in Table 2 is the ideal case where the three experiments are similar. The
expectations of the treated group were deduced from several effects relative the control group: 0% for
the verification of the type I error, then 30%, 40% and 50%. We began to study a difference of 30% because this was considered to be the minimum meaningful effect in the biological experts’ opinion. Both
historical datasets were similar with the same standard deviation for control groups (σci = 0.6)i=1,2
and treated groups (σti = 1.5)i=1,2 . Concerning the experiment of interest (experiment 3), the different
standard deviations were studied for the treated group. When there was no difference between the expectations of treated and control groups, both frequentist approaches controlled the type I error rate: the
difference between the two expectations was detected in maximum 5% of cases. Regarding the Bayesian
approach, the choice of the threshold to ensure a sufficient posterior probability of the model M1 was
important. Indeed, a threshold of 0.5 did not seem strict enough because the type I error rate was not
controlled: we concluded to a difference between the two expectations in 10% to 20% of cases while it
was not existant. However, a threshold of 0.8 allowed a better control of the type I error rate than for
the frequentist approaches, and even if it is conservative for the M0 . Whatever the considered standard
deviation for the treated group in the third experiment, the calibrated threshold verifying this error to 5%
did not exceed 0.7.
When the standard deviation of experiment 3 treated group varied from those of experiments 1 and
2 (Table 2), the power of the three approaches was affected. Indeed, it decreased when treated group
variability in step 3 increased. In all cases, the Bayesian approach was better than the Student test
performed only on the experiment 3 for each variation. The frequentist approach where data from the
three experiments were pooled was the most powerful. Indeed, this approach was very robust in this
case because to pool data from Gaussian distribution with the same expectation was like considering one
sample with an expectation and a variance from the three experiment’s variances. Thus, the Student test
was the best approach because it was in the ideal case: two Gaussian samples with the size of 30 unit.
However, the experiment of interest is the last and the two first are taken into account from an informative
point of view. We want to conclude only on the last experiment and not the all of them.
To pool datasets can influence wrongly the statistical test. The Bayesian approach uses information
from the two previous experiments to calculate the posterior probability of each model about the third
experiment. Table 3 presents the results from different situations which illustrate the historical data influence when the location parameters (µti )i=1,2,3 are different between experiments i. For each situation,
the difference between the control group expectation and the treated group expectation were distinct by
experiments. The Student test performed only on the experiment of interest was not influenced by the
10

Table 2: Efficiency (in percentage) of methods to detect a difference between the two sample expectations when
σt3 changes given a difference of δ percent. The thresholds p to respect the type I error rate to 5% for the values
of σt3 min, median and max are respectively 0.672, 0.698 and 0.660.

δ
σt3
T.test
T.test (pool)
Pr (M1 |y) > 0.5
Pr (M1 |y) > p
Pr (M1 |y) > 0.8

min
4.3
5.4
16.4
4.9
1.3

0%
med
5.0
4.3
17.1
5
1.7

max
5.0
5.3
12.7
5
2.3

min
87.1
91.4
94.1
88.3
76.7

30%
med
34.5
81.2
82.8
68.8
59.1

max
14.4
58.6
53.0
40.3
29.1

min
98.5
99.3
98.8
97.3
95.5

40%
med
59.6
97.2
94.4
89.2
84.8

max
18.3
79.7
63.8
53.5
42.0

min
100
99.9
100
100
100

50%
med
76.7
99.8
98.3
96.1
94.2

max
29.8
94.5
74.3
65.2
56.4

historical datasets because it did not take them into account. However, the Bayesian approach and the
application of the Student test on the pooled data were affected by historical data information, and concluded sometimes wrongly on the difference concerning the experiment of interest. In situations 6 to 9,
there was no distinction between treated and control group concerning the third experiment, in contradiction with historical experiments. In these cases, the approach with pooled data gave the same weight between the three experiments concluding wrongly with a high percentage. Thereby bringing no difference
between the three experiments which was not our aim in this work. The consequence of an insignificant
effect for the first two experiments (first two columns) was the decrease of the method’s power, compared to the previous results where the three experiments were similar. The Bayesian approach made a
difference between the data of each experiment, giving the conclusion on the last experiment. We noticed
a natural weighting for our method where the conclusion on the last experiment was less influenced by
the first historical experiment than the second one. Indeed, when the difference between the two groups
changed only on the experiment considered in the second step (situation 8 versus 9), the percentage to
conclude wrongly on the difference varied more strongly than when this change was performed only on
experiment of the first step (situation 6 versus 8). As mentioned before, considering a threshold of 0.5
did not allow to respect the good frequentist properties. Indeed, the type I error rate was not controlled
to 5.4%. The threshold p = 0.8 seemed more relevant to ensure the decision. The threshold p = 0.698
was also relevant to ensure the decision as the type I error rate was relatively controlled while the power
was higher than when using p = 0.8. The historical data has to be chosen appropriately so that the
information objectively helps the decision. Thus, the Bayesian approach is to be the golden mean of the
two frequentist approaches using the Student test. Indeed, the information from historical experiments
impacted on the posterior probability with a balanced influence where the experiment used in the second
step was more informative than the first one associated with the first step. The Bayesian method takes
into account the information provided by the historical data. It is important to pay attention to the experiments set a priori as they affect the results. The Student test on pooled data is not suitable for our aim,
it is not taken into account in the next studies.
Table 4 shows the influence of the variability of two prior experiments on the Bayesian method. Using
a median variability for the prior experiments, the results are similar to those in Table 2 where σt3 = 1.5
(median value). Also, the variability modification of these experiments did not change Student test
power because it was applied only on the third experiment. The power of the Bayesian method seemed
to increase when the variability increases for the experiment 1. In this case, the prior of the first step was
more vague and allowed to calibrate an objective prior more efficiently in the second step. However, the
11

Table 3: Efficiency (in percentage) of methods to detect a difference between the two sample expectations when
the location parameter (µti )i=1,2,3 changes from (µci )i=1,2,3 = 2.94 between experiments following the percentage δ. The variability of three experiments being equivalent, the threshold associated with the posterior probability
is p = 0.698.

Situations
δ of Experiment 1
δ of Experiment 2
δ of Experiment 3
T.test
T.test (pool)
Pr (M1 |y) > 0.5
Pr (M1 |y) > p
Pr (M1 |y) > 0.8

1
10%
10%
40%
53.7
44.1
67.4
44.7
32.4

2
10%
10%
50%
74.9
58.1
78.5
59.4
45.6

3
30%
30%
10%
8.3
61
57.3
38.7
29.0

4
50%
30%
10%
8.2
80.7
45.9
31.9
24.1

5
50%
10%
10%
8.9
56.6
21.8
13.1
8.9

6
50%
10%
0%
5.4
44.0
13.7
6.1
3.4

7
50%
30%
0%
5.5
66.2
30.8
19.9
15.2

8
10%
10%
0%
4.3
9.6
20.7
6.8
3.5

9
10%
20%
0%
6.1
16.2
31.9
13.1
7.9

Bayesian method power decreased when the variability increased for the second experiment. Indeed, the
weight of the second experiment was more important than the first, and a prior less diffused improved the
method power. Thereby, the variability of the second experiment impacted the power more importantly
than the first experiment variability.
Table 4: Efficiency (in percentage) of methods to detect a difference between the two sample expectations given
a difference of δ percent, when (σti )i=1,2 changes, and σt3 = (σtk )k6=i = 1.5.

i=1

i=2

δ
σti
Pr (M1 |y) > 0.5
Pr (M1 |y) > p
Pr (M1 |y) > 0.8
Pr (M1 |y) > 0.5
Pr (M1 |y) > p
Pr (M1 |y) > 0.8

0%
min
8.4
5 (p=0.59)
1.8
13.6
5 (p=0.67)
1.2

max
13.4
5 (p=0.618)
1.0
15.2
5 (p=0.703)
2.0

30%
min max
57.2 79.0
53.4 66.6
40.1 43.5
93.2 53.5
86.8 36.8
77.4 28.1

40%
min max
75.0 93.5
71.9 87.9
61.7 71.1
98.7 64.2
97.6 47.7
94.7 41.3

50%
min max
87.3 97.9
84.9 95.2
79.0 88.8
99.9 74.5
99.5 63.5
99.1 56.2

Our method was sensitive to prior experiments. It is important to choose similar experiments where
the variability and treatment effect are close. This proximity between experiments is all the more important between the second experiment and the experiment of interest.

3.2

Influence of weighting for historical experiments

Although, the three step procedure includes a natural weighting through the updating step of the posterior
probability calculation in the model space, an additional weighting could be applied to decrease the prior
information influence. The aim is to downweight the historical data to some degree (Ibrahim and Chen,
2000), by for example reducing the sample size using a multiplication factor between 0 and 1 given
the same empirical expectations and variances. We characterised this weighting (in percent) with the
couple (w1 , w2 ) being the two weights respectively associated with experiment 1 and 2. They allowed
to decrease the sample size of these prior experiments in the computation and caused more vague prior
12

distributions on the parameters. In prior distributions for the third step (π03 (θ0 ) for M1 with the equation
(7), and π03 (θ0 ) for M0 with the equation (10)), considering the quantities (w1 nc1 , w1 nt1 , w2 nc2 , w2 nt2 )
instead of (nc1 , nt1 , nc2 , nt2 ) impacted the deviation parameters of the prior distributions (Gaussian and
inverse gamma distributions) and caused diffused priors. The previous simulation results showed a natural weighting between the two historical experiments: the experiment used in the second step of the
Bayesian method influenced more than the one used in the first step. Thus, this natural weighting is
pertinent and suitable. In this subsection, we considered an additional weighting with the same weights
for both historical experiments (i = 1, 2) such as w = w1 = w2 . The previous simulation studies were
performed once more with the two following weights: w = 1/3 and w = 1/2.
Table 5 is equivalent to Table 2 with the weighting. The thresholds calibrated to control the type I
error rate to 5% increased with the incertitude. Indeed, when the three experiments were similar, these
thresholds were 0.8 and 0.9 respectively for the weights w = 1/2 and w = 1/3. Since the percentage
to conclude wrongly on the different expectations increased with the incertitude, the lower the weight
was, the higher the threshold to be considered. Therefore, regarding δ > 0, the results were close with
and without weights w. With weight, the percentage to detect a difference between expectation seemed
to increase for a higher value of σt3 and seemed to decrease for a lower value of σt3 in comparison with
results of Table 2.
Table 5: Efficiency (in percentage) of methods to detect a difference between the two sample expectations when
σt3 changes, given a difference of δ percent and the weighting w. The thresholds p to respect the type error to 5%
for the values of σt3 min, median and max are respectively 0.82, 0.89.6 and 0.73 for w = 31 and respectively 0.753,
0.798, 0.669 for w = 21 .

Pr
Pr
Pr
Pr
Pr
Pr

δ
σt3
M1 |y, w =
M1 |y, w =
M1 |y, w =
M1 |y, w =
M1 |y, w =
M1 |y, w =

1
3
1
3
1
3
1
2
1
2
1
2



> 0.5
>p
> 0.8
> 0.5
>p
> 0.8

min
16.0
5.0
5.7
16.2
5.0
3.8

0%
median
25.5
5.0
7.4
22.1
5.0
5.0

max
10.4
5.0
3.4
9.1
5.0
2.4

min
79.0
67.6
68.5
85.6
75.6
73.1

30%
med
74.5
57.9
63.9
78.5
64.3
64.2

max
47.7
39.3
36.2
47.7
38.0
30.6

min
90.3
84.7
85.2
95.4
92.4
91.0

40%
med
89.1
80.1
83.5
91.2
84.3
84.3

max
63.5
57.8
55.1
63.3
56.9
49.2

min
96.2
94.0
94.2
98.1
97.0
96.7

50%
med
93.7
89.0
91.6
96.3
93.5
93.3

max
74.0
69.2
67.5
75.6
70.7
66.5

Table 6 is equivalent to Table 4 of the previous subsection. The Bayesian method seemed conservative with the M0 when the weighting was introduced. Thus, we found the same conclusions than the
simulations without the weights, but the Bayesian method power was lower with the weight than without.
In light of these results, we advise to consider a higher threshold with the weighting than without. The
Bayesian method using p = 0.9 was more powerful than Student test and ensured the good frequentist
properties.
Table 7 shows the trend for the Bayesian approach to detect wrongly the difference when the weighting is taken into account (situations 6 to 9). Contrary to expectation, the proposed method was influenced
more by the prior experiments given the weighting than without. Indeed, for the situations 1 and 2, the
power to detect a difference between the two groups for the third experiment decreased with the weighting (Table 7 versus Table 3). In contrast, the situations 3 and 4 showed the power was more important
with weighting than without weighting, while the difference between the two groups was mainly observed for the prior experiments.
The different simulation studies showed the importance of the choice of experiments to calibrate the
13

Table 6: Efficiency (in percentage) of methods to detect a difference between the two sample expectations when
(σti )i=1,2 changes, given a difference of δ percent and the weighting w. N = 1000 is the simulation number
where : n = 10, µc = 2.94, σc = 0.6 et σt3 = (σtk )k6=i = 1.5. Min corresponds to σti = 0.6 and max to σti = 3.

i=1

i=2

i=1

i=2

Pr
Pr
Pr
Pr
Pr
Pr
Pr
Pr
Pr
Pr
Pr
Pr

Effect
σti
M1 |y, w
M1 |y, w
M1 |y, w
M1 |y, w
M1 |y, w
M1 |y, w
M1 |y, w
M1 |y, w
M1 |y, w
M1 |y, w
M1 |y, w
M1 |y, w

=
=
=
=
=
=
=
=
=
=
=
=

1
3
1
3
1
3
1
3
1
3
1
3
1
2
1
2
1
2
1
2
1
2
1
2



> 0.5
>p
> 0.8
> 0.5
>p
> 0.8
> 0.5
>p
> 0.8
> 0.5
>p
> 0.8

0%
min
max
2.5
5.3
5 (p=0.187) 5 (p=0.515)
1.1
1.7
1.7
8.6
5 (p=0.255) 5 (p=0.76)
0.5
3.4
3.2
7.8
5 (p=0.367) 5 (p=0.581)
0.9
1.3
3.3
8.1
5 (p=0.429) 5 (p=0.604)
0.9
2.4

30%
min max
30.5 45.9
37.8 45.7
25.7 34.2
85.5 28.1
89.6 22.5
77.6 19.5
35.6 57.3
39.4 52.6
26.8 35.7
92.0 33.5
93.1 29.9
82.6 23.0

40%
min max
45.9 70.7
52.5 69.5
40.0 58.7
97.4 32.9
98.5 28.4
95.6 26.5
54.3 78.6
58.0 75.5
46.2 63.1
98.4 41.9
98.9 39.5
95.9 31.9

50%
min max
63.0 87.2
67.6 86.7
57.1 80.7
99.6 43.3
99.9 38.7
99.4 35.7
72.6 93.0
75.0 91.7
64.9 83.7
99.8 52.8
99.8 48.4
99.7 40.7

Table 7: Efficiency (in percentage) of methods to detect a difference between the two sample expectations, given
the weighting w, when the location parameter (µti )i=1,2,3 changes from (µci )i=1,2,3 = 2.94 between experiments
following a percentage δ.

Pr
Pr
Pr
Pr
Pr
Pr

Situations
δ of Experiment 1
δ of Experiment 2
δ of Experiment
 3
M1 |y, w = 13  > 0.5
M1 |y, w = 13  > 0.8
M1 |y, w = 13  > 0.9
M1 |y, w = 12  > 0.5
M1 |y, w = 12  > 0.8
M1 |y, w = 12 > 0.9

1
10%
10%
40%
47.6
30.4
22.7
56.0
32.3
22.3

2
10%
10%
50%
51.2
32.8
25.8
60.8
36.8
27.4

3
30%
30%
10%
67.1
53.7
45.2
67.4
45.8
36.3

4
50%
30%
10%
54.0
44.7
39.7
53.4
40.5
32.7

5
50%
10%
10%
21.7
14.2
11.8
22.7
12.9
9.3

6
50%
10%
0%
14.5
9.9
8.1
14.0
7.1
4.8

7
50%
30%
0%
48.7
37.5
31.0
43.7
28.9
22.0

8
10%
10%
0%
30.9
13.9
9.6
29.6
10.5
5.3

9
10%
20%
0%
45.0
28.8
20.5
42.3
20.1
13.0

prior. The weighting did not seem to clarify or reduce the prior information in our case, and this choice is
more delicate with the weighting. The choice of these historical experiments is important and they have
to be similar (protocol, compound, ...) to the experiment of interest.

4

Real data analysis

In this section, the Bayesian approach was applied to experimental results obtained in a model of inflammation, for which historical behavioral data are available. The objective of the study is to evaluate
the effect of test compounds in mice after intra-plantar administration of a pro-inflammatory substance,
Complete Freund’s Adjuvant (CFA). CFA rapidly produces local inflammation at the side of administration, and responses measuring spontaneous behaviors (licking, protection of the administrated paw) or
14

behaviors such as paw withdrawal evoked by mechanical of thermal stimuli. The data from these experiments are measured as reaction time (latency) of mice to withdraw their paw upon stimulation. The aim
is to compare the effects of a compound versus those of a vehicle.
In the present in vivo experiments, a radiant heat source was focused on the inflamed hind-paw of
every mouse in each group, and the measured variable was a latency to withdraw the paw from the heat
source in seconds. Previous analyses on historical data showed that normality was quite satisfying, but
there was often a problem of heterogeneity of variance. The selection of historical experiments was made
by researchers following our recommendations. The three experiments were similar: the control groups
received saline solution (vehicle) and Ibuprofen for the treated groups. The sample size was eight for
each group. Figure 1 shows the latency distribution for the three used experiments. The latency was
more important in the treated group than in the control group. In addition to the same conditions, Figure
1 shows a similar trend for the three experiments. We are in an ideal case for the application of the
three-step Bayesian method.

Figure 1: Distribution of latency in seconds for each group of the three experiments. Each experiment is composed
of two groups: control (mice receiving the vehicle) and treated (mice receiving Ibuprofen). The third experiment
(right panel) is the experiment of interest, and the experiments 1 and 2 (left and middle panels) are prior datasets.

The p-value associated with the Student test with Satterthwhaite correction was 0.042 for experiment 3. The null hypothesis was then rejected but the p-value was near to the usual threshold of 5%.
The heterogeneity between variances would have been greater given the lower sample size, the Student

15

test would reject a significant effect of Ibuprofen in comparison with control. Our Bayesian method
allows to make a statistical decision for the last experiment given the information from the two first experiments. The Bayesian method advantage is the direct interpretation of the posterior probabilities of
models. Without assuming a weighting on historical experiments since this was not pertinent based on
the simulation results, the posterior probability associated with a similar expectation for the two groups
was Pr (M0 |y) = 0.006, and the posterior probability associated with the effect of Ibuprofen compared
to the vehicle was Pr (M1 |y) = 0.994%. The posterior probability of M1 reflected the belief in difference between the two groups given the historical data. Thus, whatever the threshold p previously
discussed to validate the hypothesis H1 associated with Pr (M1 |y), we easily concluded that the latency
was different between the two groups. Both approaches concluded identically for this example: there
was an effect of Ibuprofen on the paw withdrawal latency to heat in the CFA protocol. The treated group
latency was statistically greater than the control group latency. Moreover, the Bayesian approach seemed
to be a safer method to discriminate between the two hypotheses.

5

Discussion

We propose a fully Bayesian technique to solve the Behrens-Fisher problem as an alternative to the
usual frequentist method which is not appropriate and could lack of power with small sample size. The
corresponding hypothesis testing is considered as a model choice question in the Bayesian paradigm. The
proposed method takes into account information from historical data to calibrate an objective informative
prior on parameter distribution and also on the model space in two steps. A third and later step is
performed to provide the posterior probability of each model. This quantity reflects the confidence in
the model (or hypothesis) conditionally on the data. Its interpretation is easy. Despite the small sample
sizes, the three step Bayesian method has good frequentist properties as seen in the results of simulation,
thanks to the objective informative prior built with the historical data.
Usually, the chosen model is the one with the greatest posterior probability. Considering two models,
the natural threshold of 0.5 can be used to detect an expectation difference between the two groups (M1 )
or not (M0 ). The arbitrary choice of this threshold is a difficult task. The threshold of 0.8 corresponding
to our expectation on this type of study is interesting according to our simulation studies because of its
good frequentist properties. Optimally, the threshold could be calibrated for each specific protocol in
order to control the type I error at 5%.
In most cases, the proposed methodology is more powerful than the others. However, as any Bayesian
method, it depends on a pertinent choice of historical data. In the preclinical research context, where the
in vivo pharmacology experiments are routinely performed using the same protocol, several experiments,
which meet the Pocock’s criteria (Pocock, 1976), could be used to apply the three steps methodology proposed. The method naturally weights the two historical datasets and an overweight is not needed. Indeed,
the historical data used in the second step carries more information (the informative prior probability in
the model space is built on the second experiment only). The order of the historical data used is important. We recommend, considering the second step, the most recent historical experiment to decrease
potential biases. If only one historical experiment is available, the proposed Bayesian method is also
suitable in two steps and the results are the posterior probabilities of the second step. However this
probability would not be as accurate as it would be after three steps since the informative prior on model
space would be missing.

16

As mentioned in the introduction, the hierarchical models (Neuenschwander et al., 2010) for the
prior building were not studied in this work. However, it could be interesting to expand our work with
the comparison between this kind of MAP prior and also with power prior approach (Ibrahim et al.,
2015). In this paper, we focused on the comparison of two groups. But, in most cases, the objective of
this kind of preclinical research experiments is to compare more than one group with a control one, at
least three groups are considered: control, treated and reference. It would be interesting to develop the
proposed method to meet the objective in the case of several groups (Scott and Berger, 2005). With three
groups for example, we have to discriminate between four models such as:
• Model M0 : no difference between the means of the three groups or no effect of reference compound and tested compound;
• Model M1 : just reference compound has an effect;
• Model M2 : just tested compound has an effect;
• Model M3 : both compounds have an effect versus placebo.

Acknowledgements
We would like to thank Nicolas Bonnet from Sanofi CEP statistical department in Montpellier for his
support and review. Thanks also to Noelle Boussac from Sanofi research statistical department for her
support on this work, and thanks to the researchers from Aging Pain Sanofi department in Montpellier
for whom these analyses were performed; thanks in particular to Isabel A. Lefevre for her editorial
assistance and Veronique Menet, Eric Dufour for providing the experimental data and their availability
for exchanges.

References
Box, G. E. P. and Tiao, G. C. (1973). Bayesian inference in statistical analysis. Addison-Wesley Pub.
Co.
code by Steven G. Johnson, C. and by Balasubramanian Narasimhan, R. (2013). cubature: Adaptive
multivariate integration over hypercubes. R package version 1.1-2.
Cowles, M. K. (2004). Review of WinBUGS 1.4. The American Statistician, 58(4):330–336.
Food and Drug Administration, F. (2010). Guidance for the Use of Bayesian Statistics in Medical Device
Clinical Trials.
Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models (comment on
article by Browne and Draper). Bayesian Analysis, 1(3):515–534.
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013). Bayesian
Data Analysis, Third Edition. CRC Press.

17

Ghosh, M. and Kim, Y.-H. (2001). The behrens-fisher problem revisited: A bayes-frequentist synthesis.
Canadian Journal of Statistics, 29(1):5–17.
Hobbs, B. P., Carlin, B. P., Mandrekar, S. J., and Sargent, D. J. (2011). Hierarchical commensurate and
power prior models for adaptive incorporation of historical information in clinical trials. Biometrics,
67(3):1047–1056.
Ibrahim, J. G. and Chen, M.-H. (2000). Power prior distributions for regression models. Statistical
Science, 15(1):46–60.
Ibrahim, J. G., Chen, M.-H., Gwon, Y., and Chen, F. (2015). The power prior: theory and applications.
Statistics in Medicine.
Kass, R. E. and Raftery, A. E. (1995). Bayes Factors. Journal of the American Statistical Association,
90(430):773–795.
Kery, M. (2010). Introduction to WinBUGS for Ecologists: Bayesian approach to regression, ANOVA,
mixed models and related analyses. Academic Press, Amsterdam, 1 edition edition.
Kim, S.-H. and Cohen, A. S. (1998). On the Behrens-Fisher Problem: A Review. Journal of Educational
and Behavioral Statistics, 23(4):356–377.
Landis, S. C., Amara, S. G., Asadullah, K., Austin, C. P., Blumenstein, R., Bradley, E. W., Crystal, R. G.,
Darnell, R. B., Ferrante, R. J., Fillit, H., Finkelstein, R., Fisher, M., Gendelman, H. E., Golub, R. M.,
Goudreau, J. L., Gross, R. A., Gubitz, A. K., Hesterlee, S. E., Howells, D. W., Huguenard, J., Kelner,
K., Koroshetz, W., Krainc, D., Lazic, S. E., Levine, M. S., Macleod, M. R., McCall, J. M., Moxley Iii,
R. T., Narasimhan, K., Noble, L. J., Perrin, S., Porter, J. D., Steward, O., Unger, E., Utz, U., and
Silberberg, S. D. (2012). A call for transparent reporting to optimize the predictive value of preclinical
research. Nature, 490(7419):187–191.
Marin, J.-M. and Robert, C. (2013). Bayesian Essentials with R. Springer Science & Business Media.
Moreno, E., Bertolino, F., and Racugno, W. (1999). Default Bayesian analysis of the Behrens–Fisher
problem. Journal of Statistical Planning and Inference, 81(2):323–333.
Neuenschwander, B., Capkun-Niggli, G., Branson, M., and Spiegelhalter, D. J. (2010). Summarizing
historical information on controls in clinical trials. Clinical Trials (London, England), 7(1):5–18.
Pocock, S. J. (1976). The combination of randomized and historical controls in clinical trials. Journal of
Chronic Diseases, 29(3):175–188.
Robert, C. (2007). The Bayesian Choice: From Decision-Theoretic Foundations to Computational
Implementation. Springer Science & Business Media.
Schmidli, H., Gsteiger, S., Roychoudhury, S., O’Hagan, A., Spiegelhalter, D., and Neuenschwander, B.
(2014). Robust meta-analytic-predictive priors in clinical trials with historical control information.
Biometrics, 70(4):1023–1032.

18

Scott, J. G. and Berger, J. O. (2005). An exploration of aspects of Bayesian multiple testing. Journal of
Statistical Planning and Inference.
Viele, K., Berry, S., Neuenschwander, B., Amzal, B., Chen, F., Enas, N., Hobbs, B., Ibrahim, J. G., Kinnersley, N., Lindborg, S., Micallef, S., Roychoudhury, S., and Thompson, L. (2014). Use of historical
control data for assessing treatment effects in clinical trials. Pharmaceutical Statistics, 13(1):41–54.
Whitehead, J., Cleary, F., and Turner, A. (2015). Bayesian sample sizes for exploratory clinical trials
comparing multiple experimental treatments with a control. Statistics in Medicine, 34(12):2048–2061.
Whitehead, J., Valdés-Márquez, E., Johnson, P., and Graham, G. (2008). Bayesian sample size for
exploratory clinical trials incorporating historical data. Statistics in Medicine, 27(13):2307–2327.

19

