Simultaneous Inference
in General Parametric Models
∗
Torsten Hothorn
Institut f
¨
ur Statistik
Ludwig-Maximilians-Universit
¨
at M
¨
unchen
Ludwigstraße 33, D–80539 M
¨
unchen, Germany
Frank Bretz
Statistical Methodology, Clinical Information Sciences
Novartis Pharma AG
CH-4002 Basel, Switzerland
Peter Westfall
Texas Tech University
Lubbock, TX 79409, U.S.A
March 15, 2013
Abstract
Simu ltaneous inference is a common problem in many areas of application. If
multiple null hypotheses are tested simultaneously, the probability of rejecting er-
roneously at least one of them increases beyond the pre-specified significance level.
Simu ltaneous inference procedures have to be used which adjust for multiplicity and
thu s control the overall type I error rate. In this paper we describe simultaneous infer-
ence procedures in general parametric models, where the experimental questions are
specified through a linear combination of elemental model parameters. The frame-
work described here is quite general and extends the canonical theory of multiple
comparison procedures in ANOVA models to linear regression problems, generalized
linear models, linear mixed effects models, the Cox model, robust linear models, etc.
Seve ral examples using a variety of different statistical models illustrate the breadth
∗
This is a preprint of an article published in Biometrical Journal, Volume 50, Number 3, 346–363.
Copyright
➞
2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim; available online
http://www.
biometrical-journal.com
.
1
of the results. For the analyses we use the R add-on package multcomp, which pro-
vides a convenient interface to the general approach adopted here.
Key words: multiple tests, multiple comparisons, simultaneous confidence intervals,
adjusted p-values, multivariate normal distribution, robust statistics.
1 Introduction
Multiplicity is an intrinsic problem of any simultaneous inference. If each of k, say, null
hypotheses is tested at nominal level α, the overall type I error rate can be substantially
larger than α. That is, the probability of at least one erroneous rejection is larger tha n
α for k ≥ 2. Common multiple comparison procedures adjust for multiplicity and thus
ensure that the overall type I error remains below the pre-sp ecifie d significance level α.
Examples of such multiple comparison procedures include Dunnett’s many-to-one compar-
isons, Tukey’s all-pairwise comparisons, sequential pairwise contrasts, comparisons with
the average, changep oint analyses, dose-response contrasts, etc. These procedures are all
well established for classical regression and ANOVA models allowing for cova riates and/or
factorial treatment structures with i.i.d.˜normal errors and constant variance, see Bretz
et˜al.
(2008) and the references therein. For a general reading on multiple co mparison
procedures we refer to
Hochberg and Tamhane (1987) and Hsu (1996).
In this paper we aim at a unified description of simultaneous inference procedures in para-
metric models with generally correlated parameter estimates. Each individual null hypothe-
sis is specified through a linear combination of elemental model parameters and we allow for
k of such null hypotheses to be tested simultaneously, regardless of the number of elemental
model parameters p. The general framework described here extends the current canoni-
cal theory with respect to the following aspects: (i) model assumptions such as normality
and homoscedasticity are relaxed, thus allowing for simultaneous inference in generalized
linear models, mixed effects models, survival models, etc.; (ii) arbitrary linear functions of
the elemental parameters are allowed, not just contrasts of means in AN(C)OVA models;
(iii) computing the r eference distribution is feasible for arbitrary designs, especially for
unbalanced designs; and (iv) a unified implementation is provided which allows for a fast
transition of the theoretical results to the desks of data analysts interested in simultaneous
inferences for multiple hypotheses.
Accordingly, the paper is organized as follows. Section˜
2 defines the general model and ob-
tains the asymptotic or exact distribution of linear functions of elemental model parameters
under rather weak conditions. In Section˜
3 we describe the framework for simultaneous
inference procedures in general parametric models. An overview about important applica-
tions of the methodology is given in Section˜4 followed by a s hort discussion of the software
implementation in Section˜
5. Most interesting from a practical point of view is Section˜6
where we analyze four rather challenging problems with the tools developed in this paper.
2 Model and Parameters
In this section we introduce the underlying model assumptions and derive some asymptotic
results necessary in the subsequent sections. The results from this section form the basis
for the simultaneous inference procedures described in Section˜3.
Let M((Z
1
, . . . , Z
n
), θ, η ) denote a semi-parametric statistical model. The set of n obser-
vations is described by (Z
1
, . . . , Z
n
). The model contains fixed but unknown elemental
parameters θ ∈ R
p
and other (random or nuisance) parameters η. We are primarily in-
terested in the linear functions ϑ := Kθ of the parameter vector θ as specified through
the consta nt matrix K ∈ R
k,p
. In what follows we describe the underlying model assump-
tions, the limiting distribution of estimates of our parameters of interest ϑ, as well as the
corresponding test statistics for hypotheses about ϑ and their limiting joint distribution.
Suppose
ˆ
θ
n
∈ R
p
is an estimate of θ and S
n
∈ R
p,p
is an estimate of cov(
ˆ
θ
n
) with
a
n
S
n
P
−→ Σ ∈ R
p,p
(1)
for some positive, nondecreasing sequence a
n
. Furthermore, we assume that a multivariate
central limit theorem holds, i.e.,
a
1/2
n
(
ˆ
θ
n
− θ)
d
−→ N
p
(0, Σ). (2)
If both (
1) and (2) are fulfilled we write
ˆ
θ
n
a
∼ N
p
(θ, S
n
). Then, by Theorem 3.3.A in Serfling
(1980), the linear function
ˆ
ϑ
n
= K
ˆ
θ
n
, i.e., an estimate of our parameters of interest, also
follows an approximate multivariate normal distribution
ˆ
ϑ
n
= K
ˆ
θ
n
a
∼ N
k
(ϑ, S
⋆
n
)
with covariance matrix S
⋆
n
:= KS
n
K
⊤
for any fixed matrix K ∈ R
k,p
. Thus we need not
to distinguish between elemental parameters θ or derived parameters ϑ = Kθ that are of
interest to the r esearcher. Instead we simply assume for the moment that we have (in
analogy to (
1) and (2))
ˆ
ϑ
n
a
∼ N
k
(ϑ, S
⋆
n
) with a
n
S
⋆
n
P
−→ Σ
⋆
:= KΣK
⊤
∈ R
k,k
(3)
and that the k parameters in ϑ are themselves the parameters of interest to the researcher.
It is assumed that the diagonal elements of the covariance matrix are positive, i.e., Σ
⋆
jj
> 0
for j = 1, . . . , k.
Then, the standardized estimator
ˆ
ϑ
n
is again asymptotically normally distributed
T
n
:= D
−1/2
n
(
ˆ
ϑ
n
− ϑ)
a
∼ N
k
(0, R
n
) (4)
where D
n
= diag(S
⋆
n
) is the diagonal matrix given by the diagonal elements of S
⋆
n
and
R
n
= D
−1/2
n
S
⋆
n
D
−1/2
n
∈ R
k,k
1
is the correlation matrix of the k-dimensional statistic T
n
. To demonstrate (4), note that
with (3) we have a
n
S
⋆
n
P
−→ Σ
⋆
and a
n
D
n
P
−→ diag(Σ
⋆
). Define the sequence ˜a
n
needed to
establish ˜a-convergence in (
4) by ˜a
n
≡ 1. Then we have
˜a
n
R
n
= D
−1/2
n
S
⋆
n
D
−1/2
n
= (a
n
D
n
)
−1/2
(a
n
S
⋆
n
)(a
n
D
n
)
−1/2
P
−→ diag(Σ
⋆
)
−1/2
Σ
⋆
diag(Σ
⋆
)
−1/2
=: R ∈ R
k,k
where the convergence in probability to a constant follows from Slutzky’s Theorem (The-
orem 1.5.4,
Serfling, 1980) and therefore (4) holds. To finish note that
T
n
= D
−1/2
n
(
ˆ
ϑ
n
− ϑ) = (a
n
D
n
)
−1/2
a
1/2
n
(
ˆ
ϑ
n
− ϑ)
d
−→ N
k
(0, R).
For the purposes of multiple comparisons, we need convergence of multivariate probabilities
calculated for the vector T
n
when T
n
is assumed normally distributed with R
n
treated
as if it were the true correlation matrix. However, such probabilities P(max(|T
n
| ≤ t)
are continuous functions of R
n
(and a critical value t) which converge by R
n
P
−→ R as
a consequence of Theorem 1.7 in
Serfling (1980). In cases where T
n
is assumed multi-
variate t distributed with R
n
treated as the estimated correlation matrix, we have similar
convergence as the degrees of freedom approach infinity.
Since we only assume that the parameter estimates are asymptotically normally distributed
with a consistent estimate of the associated covariance matrix being available, our frame-
work covers a large class of statistical models, including linear regression and ANOVA
models, generalized linear models, linear mixed effects models, the Cox model, robust lin-
ear models, etc. Standard software packages can be used to fit such models and obtain
the estimates
ˆ
θ
n
and S
n
which are essentially the only two quantities that are needed for
what follow s in Section˜
3. It should be noted that the elemental parameters θ are not
necessarily means or differences of means in AN(C)OVA models. Also, we do not restrict
our attention to contrasts of such means, but allow for any set of constants leading to the
linear functions ϑ = Kθ of interest. Specific examples for K and θ will be given later in
Sections˜
4 and 6.
3 Global and Simultaneous Inference
Based on the results from Section˜
2, we now focus on the derivation of suitable inference
procedures. We start considering the general linear hypothesis (
Searle, 1971) formulated
in terms of our parameters of interest ϑ
H
0
: ϑ := Kθ = m.
2
Under the conditions of H
0
it follows from Section˜2 that
T
n
= D
−1/2
n
(
ˆ
ϑ
n
− m)
a
∼ N
k
(0, R
n
).
This approximating distribution will now be used as the reference distribution when con-
structing the inference procedures. The global hypothesis H
0
can be tested using standard
global tests, such as the F - or the χ
2
-test. An alternative approach is to use maximum
tests, as explained in Subsection˜
3.1. Note that a small global p-value (obtained from one
of these procedures) leading to a rejection of H
0
does not give further indication about
the nature of the significant result. Therefore, one is often interested in the individual null
hypotheses
H
j
0
: ϑ
j
= m
j
.
Testing the hypotheses set {H
1
0
, . . . , H
k
0
} simultaneously thus requires the individual as-
sessments while maintaining the familywise error rate, as discussed in Subsection˜3.2
At this point it is worth considering two special cases. A stronger a ssumption than asymp-
totic normality of
ˆ
θ
n
in (
2) is exact normality, i.e.,
ˆ
θ
n
∼ N
p
(θ, Σ). If the covariance matrix
Σ is known, it follows by standard arguments that T
n
∼ N
k
(0, R), when T
n
is normalized
using fixed, known variances. Otherwise, in the typical situation of linear models with
normal i.i.d. errors, Σ = σ
2
A, where σ
2
is unknown but A is fixed and known, the exact
distribution of T
n
is a k-dimensional multivariate t
k
(ν, R) distribution with ν degrees of
freedom (ν = n − p − 1 for linear models), see
Tong (1990 ).
3.1 Global Inference
The F - and the χ
2
-test are classical approaches to assess the global null hypothesis H
0
.
Standard results (such as Theorem 3.5,
Serfling, 1980) ensure that
X
2
= T
⊤
n
R
+
n
T
n
d
−→ χ
2
(Rank(R)) when
ˆ
θ
n
a
∼ N
p
(θ, S
n
)
F =
T
⊤
n
R
+
T
n
Rank(R)
∼ F(Rank(R), ν) when
ˆ
θ
n
∼ N
p
(θ, σ
2
A),
where Rank(R) and ν are the corresponding degrees of freedom of the χ
2
and F distri-
bution, respectively. Furthermore, Rank(R
n
)
+
denotes the Moore-Penrose inverse of the
correlation matrix Rank(R).
Another suitable scalar test statistic for testing the global hypothesis H
0
is to consider
the maximum of the individual test statistics T
1,n
, . . . , T
k,n
of the multivariate statistic
T
n
= (T
1,n
, . . . , T
k,n
), leading to a max-t type test statistic max(|T
n
|). The distribution
of this statistic under the conditions of H
0
can be handled through the k-dimensional
distribution
P(max(|T
n
|) ≤ t)
∼
=
t
Z
−t
· · ·
t
Z
−t
ϕ
k
(x
1
, . . . , x
k
; R, ν) dx
1
· · · dx
k
=: g
ν
(R, t) (5)
3