Fixed effects testing in high-dimensional linear mixed
models
Jelena Bradic1 and Gerda Claeskens2 and Thomas Gueuning∗2
1

arXiv:1708.04887v1 [stat.ME] 14 Aug 2017

2

Department of Mathematics, University of California San Diego
ORStat and Leuven Statistics Research Center, KU Leuven, Belgium

Abstract
Many scientific and engineering challenges – ranging from pharmacokinetic drug dosage allocation and personalized medicine to marketing mix (4Ps) recommendations – require an understanding of the unobserved heterogeneity in order to develop the best decision making-processes.
In this paper, we develop a hypothesis test and the corresponding p-value for testing for the
significance of the homogeneous structure in linear mixed models. A robust matching moment
construction is used for creating a test that adapts to the size of the model sparsity. When unobserved heterogeneity at a cluster level is constant, we show that our test is both consistent and
unbiased even when the dimension of the model is extremely high. Our theoretical results rely
on a new family of adaptive sparse estimators of the fixed effects that do not require consistent
estimation of the random effects. Moreover, our inference results do not require consistent model
selection. We showcase that moment matching can be extended to nonlinear mixed effects models
and to generalized linear mixed effects models. In numerical and real data experiments, we find
that the developed method is extremely accurate, that it adapts to the size of the underlying
model and is decidedly powerful in the presence of irrelevant covariates.

Keywords: p-values, random effects, misspecification, penalization, robustness

1

Introduction

In many applications, we want to use data to draw inferences about the common underlying effect
of a treatment. Examples include medical studies about the computation of effect sizes for assessing
the clinical importance of covariates on pharmacokinetic or pharmacodynamic responses, and to the
study of drug-drug interactions or studies that quantify the effects of different advertising mediums
that take into account other variables such as pricing, distribution points and competitor tactics
for example, commonly used by technology firms to optimize budget over these different mediums.
Historically, most datasets have been too small to meaningfully explore unobserved heterogeneity
beyond dividing the sample into independent subgroups. Recently, however, there has been an
explosion of empirical settings where it is potentially feasible to gather large-scale observations and
therefore better customize estimates across both population and individuals.
∗
Support of the Research Foundation Flanders, KU Leuven grant GOA/12/14 and of the IAP Research Network
P7/06 of the Belgian Science Policy is acknowledged. Jelena Bradic would like to acknowledge the support of the
National Science Foundation grant DMS - 1712481. The computational resources and services used in this work
were provided by the VSC (Flemish Supercomputer Center), funded by the Hercules Foundation and the Flemish
Government - department EWI.

1

An impediment to exploring unobserved heterogeneous effects is the fear that researchers will
iteratively search for subgroups with high treatment levels, and then report only the results for subgroups with extreme effects, thus exploring and utilizing heterogeneity that may be purely spurious.
Moreover, procedural restrictions have often been used to control for the unobserved randomness.
However, such procedural restrictions can make it difficult to encompass strong but unexpected heterogeneity which naturally occurs in practice. In this paper, we seek to address this challenge by
developing a method for hypothesis testing of fixed effects, while allowing models to have heterogeneous and random components, that yields valid asymptotic confidence intervals and p-values for the
true underlying fixed effect.
Classical approaches to test fixed effects in the presence of random effects include Breslow and
Clayton (1993), Kenward and Roger (1997) and Crainiceanu and Ruppert (2004); see, e.g., Verbeke
and Molenberghs (2009). These methods perform well in applications with a small number of covariates, but quickly break down as the number of covariates increases. In this paper, we explore the use
of doubly robust ideas from the literature to improve the performance of these classical methods in
the presence of an exploding number of covariates. We develop a family of moment matching tests,
which allows for flexible modeling of interactions in high dimensions by allowing models that are not
entirely sparse (not even approximately). The developed moments are related to Neyman orthogonalization principles in that they are robust to model misspecifications and estimation of nuisance
parameters; however, the tests differ in that they are also robust to the misspecification (and/or
mis-estimation) of random effects in the case of linear mixed models – this is especially important in
environments with many covariates with possible complex interactions and the presence of unobserved
model heterogeneity.
Despite their widespread success in estimation and inference in high-dimensional linear models,
there are important hurdles that need to be cleared before penalized estimators are directly useful
in linear mixed models. Ideally, an estimator of complete variability in the model needs to be
constructed, so that a researcher can use it to test hypothesis and establish sampling distribution.
However, in the case of mixed models, estimation of random effects is typically more difficult than
that of fixed effects one is interested in. Developing inferential tools that do not necessarily rely on
good quality estimates of random effects would therefore be extremely useful. Yet, such procedures
have been largely left undeveloped, even in the standard contexts.
This paper addresses these limitations, develops a robust method for fixed effect testing that
allows for a tractable asymptotic theory and valid statistical inference even when the estimation of
random effects is not consistent. Our proposed test is composed of an estimator of the correlation
between the model error under the null hypothesis and the error of carefully constructed feature
projections.
In the interest of generality, we begin our theoretical analysis by developing the desired consistency
and asymptotic normality results in the context of high-dimensional linear mixed models. We prove
these results for a particular variant of mixed models that uses Gaussian random effects with unknown
cluster variances while allowing the number of fixed effects to be much larger than the sample size.
However, the results do not rely on Gaussianity or quality of estimators of the unknown variance
components. This property is achieved by the doubly robust construction of the test statistic, as well
as new high-dimensional estimators of the fixed effects. We also show that such robustly motivated
estimator is consistent in l1 norm whenever the actual model is sparse. Our proof is built on within
cluster analysis and a martingale representation of the test statistic between the clusters. We show
that the consistency of the test adapts to the quality of estimation of the fixed effects by successfully
leveraging correlation properties among the features. Given these general results, we next show that
our ideas extend from the linear setting to the nonlinear setting as well as settings with unknown
cluster correlation. We also illustrate that the proposed fixed effect estimator can be utilized for
estimation of the unknown variance components.
2

Although our main focus in this paper is the inference on fixed or common effects, we note that
there are a variety of important applications of the p-value construction in a pure variable selection
context. Ryzhov et al. (2016) seek to improve the construction of marketing campaigns for non-profit
organizations by detecting which marketing strategy yielded a better donor retention (in terms of
continuous influx of small to medium donations). Donors are grouped naturally by the observed
frequency of donations thus yielding a natural longitudinal structure in the observations. Here we
need rigorous predictions for the probability that a donor would continue its activity if a specific
marketing strategy was selected as the most beneficial for the observed data. Our results would be
one of the first to develop rigorous variable selection that enables the use of large-scale data for this
purpose.

1.1

Related Work

There has been a longstanding understanding in the high-dimensional statistics literature that prediction methods based on regularization (e.g. Tibshirani (1996) and Fan and Li (2001)) perform
well outside of the class of linear models: If the goal is prediction, then we should define a proper
likelihood function and the method will be considered as good as long as the likelihood function is
convex (see Bühlmann and Van De Geer, 2011, for more details). However, good performance in prediction does not necessarily translate into good performance for estimation or inference about model
parameters. In fact, it can often be quite poor. In a Neyman orthogonalization framework we use to
formalize our inferential results, we show that the testing of significance of fixed effect can be done
asymptotically exactly regardless of the error in estimation of the random effects or the additional
nuisance parameters of the model. Thus, when evaluating estimators of the fixed effects, asymptotic
theory plays a much more important role than in the standard prediction context.
From a technical point of view, the main contribution of this paper is an asymptotic normality
theory enabling statistical inference in the context of high-dimensional linear mixed models while
simultaneously allowing misspecification in both the model and the random effects. Recent results
by Zhang and Zhang (2014); Van de Geer et al. (2014); Javanmard and Montanari (2014a); Cai et al.
(2017); Belloni et al. (2014), and others have established asymptotic properties of tests in a particular
variant of the sparse high-dimensional linear models Ren et al. (2015); Jankova et al. (2015); Ning
et al. (2017); Athey et al. (2016). To our knowledge, however, we provide the first set of conditions
under which tests are both asymptotically unbiased and Gaussian for the linear mixed models, thus
allowing for classical statistical inference; the estimator of the fixed effects used to achieve asymptotic
normality proposed in this paper is also new. We review the existing theoretical literature on linear
mixed models in more detail in Section 3.1.
A small but growing literature, including Belloni et al. (2015); Chernozhukov et al. (2015), has
considered the use of Neyman orthogonalization for the purposes of high-dimensional inference. These
papers use the Neyman orthogonalization method of Neyman (1959) together with sample splitting
and de-biasing (e.g. Chernozhukov et al. (2016)), and report confidence intervals or p-values resulting
from significance testing, obtained by Belloni et al. (2017). A limitation of the existing work is in that
it cannot allow for the unobserved heterogeneity in the model (like the presence of random effects for
example).
We view our contribution as complementary to this literature, by showing that Neyman orthogonalization need not only be viewed as successful in linear models (or generalized linear models), and
can instead be modified and used for rigorous asymptotic analysis in high-dimensional model with
random effects. We succeed in showing that our construction allows one more degree of robustness –
the random effects can be estimated rather poorly without changing the asymptotic null distribution
of the test. Even in low dimensions advancements of this type would be of considerable interests as
there are no fully reliable ways to identify the best covariance model; as many studies have revealed
biased inference for the fixed effects with an underspecified covariance structure; the problem is only
3

multiplied when there is a growing number of fixed effects in the model. The methodological and
theoretical tools developed here are useful beyond the specific class of algorithms studied in our paper.
In particular, our tools allow for a fairly direct analysis of variants of the Bayesian hierarchical linear
regression models (e.g. Quintana et al. (2016)).
Finally, we note a growing literature on estimating fixed effects in the presence of unobserved
heterogeneity using different regularization methods. Schelldorfer et al. (2011); Groll and Tutz (2014);
Hui et al. (2016) develop lasso-like methods for estimation of fixed effects in a sparse high-dimensional
linear mixed model setting. Non-convex methods were investigated by Wang et al. (2012) and Ghosh
and Thoresen (2016) among others. Fan and Li (2012) and others discuss variable selection procedures
for the fixed or random effects that enable off-the-shelf methods to be used for optimal estimation
(see a review article Müller et al. (2013) for more details). Additionally, Bonnet et al. (2015) estimate
heritability in sparse linear mixed models, relying on a special asymptotic scenario where the sizes of
the clusters are proportional to the number of observations. However, none of the approaches above
provides guarantees about honest p-values or confidence interval constructions.

2

Inference in Linear Mixed Models

We consider a linear mixed model in which observations are grouped. Suppose that we have a sample
of N subjects. Let i = 1, . . . , N be the grouping indexPand ni the number of observations in group i.
The total number of observations is denoted by n = N
i=1 ni . For the i-th subject we represent the
model as follows
yi = X i γ ∗ + Z i β ∗ + Wi bi + i
(1)
where yi = (yi1 , . . . , yini )> ∈ Rni collects the response variable of the i-th subject, X i ∈ Rni ×(p−1)
p−1 for j ∈
>
> > the design matrix of the fixed effects where x
with X i = (x>
ij ∈ R
i1 , xi2 , . . . , xini )
{1, . . . , ni }. Similarly we denote with Z i = (zi1 , zi2 , . . . , zini )> the design vector of the fixed effects
which we are interested in, with zij ∈ R. The vector of population specific fixed effect coefficients is
split into a (p − 1)-dimensional vector γ ∗ and a univariate β ∗ . The subject specific random effects are
defined through a q dimensional vector bi = (bi1 , . . . , biq )> ∈ Rq for which we assume bi ∼ Nq (0, ψ).
Here, ψ is the unknown covariance matrix of the random effects which can be correlated with each
other. The corresponding deterministic design matrix for group i is denoted with W i ∈ Rni ×q . In
addition, i is the random error vector with components independent and identically distributed with
mean zero and an unknown variance 0 < σ2 < ∞. We would like to remark that the independence
assumption can be generalized to a positive definite variance covariance structure. This generalization
still fits into the theoretical framework presented in Section 3. Nonetheless, for the sake of notational
simplicity, we restrict to the case of group independent errors.
The goal of this paper is to develop a test which is able to detect wether β ∗ is equal to β0 or not,
that is,
H0 : β ∗ = β0 versus H1 : β ∗ 6= β0 ,
(2)

for some given value β0 ∈ R. We study the asymptotics with N → ∞ and assume that ni is bounded in
√
N , while allowing a large number of fixed effects in that p  n, i.e. log p = o( n). The main difficulty
is that we can only ever test the fixed effects if the random effects are either known or estimated well.
However, in practice this is never achievable and we cannot directly test for the fixed effects using
existing tools and techniques. In general, we cannot estimate the variance components consistently
simply from the observed data without further restrictions on the data generating distribution. A
standard way to make progress is to assume models selection consistency of the fixed effects, i.e.,
that estimated fixed effects are correctly selected. The motivation behind this assumption is that
it enables direct dimensionality reduction as it effectively implies correct selection of the true fixed
effects; thus, imposing restrictive minimal signal strength assumptions together with a irrepresentable
4

condition. In this paper, we take a more indirect approach: We show that, under simple assumptions,
our approach can use moment conditions to achieve consistency of a test without needing to explicitly
estimate the variance component.
Let vectors Y , b and , and matrices X, Z be obtained by stacking vectors yi , bi and i and
matrices X i , Z i respectively, underneath each other, and let Ψ = diag(ψ, ..., ψ) ∈ RqN ×qN so that
W b ∼ Nn (0, W ΨW > ) with a block-diagonal matrix W = diag(W 1 , . . . , W N ) ∈ Rn×qN . In partic> >
qN is the random effect vector and Y |(X, Z) has mean Xγ ∗ + Zβ ∗ and
ular, b = (b>
1 , . . . , bN ) ∈ R
>
2
variance σ I n + W ΨW . We further standardize the design matrix X such that each column has
the same norm. The linear mixed effects model (1) can be rewritten as
Y = Xγ ∗ + Zβ ∗ + W b + .

(3)

The n components of the noise vector  are independent and identically distributed with mean zero
and variance 0 < σ2 < ∞. We assume that , b and W are mutually independent and that  is
independent of X, Z. Observe that we do not require the error to have Gaussian distribution (see
Schelldorfer et al., 2011; Fan and Li, 2012, where Gaussianity was explicitly assumed); if we restrict
our attention to Gaussian error then the independence above can be replaced with uncorrelatedness
with no significant changes in the proofs. Note that γ ∗ and β ∗ can be seen as being derived from an
original p-dimensional vector. The notation (3) emphasizes the fact that β ∗ is the component of this
original vector on which we want to perform hypothesis testing.

2.1

From the parametric null to the matching moments condition

At a high level, Wald and score tests can be thought of as two contrasting methods with likelihood
based constructions. Given a particular score function (typically a first derivative of the likelihood
function), a classical method such as Rao’s test, performs estimation only under the null hypothesis
and forms a test based on the first moment condition, that is the expectation of the score should
be close to zero. In contrast, orthogonalization methods seek to find appropriate directions that are
close to the score, where the closeness is defined with respect to a projection of certain kind. The
advantage of orthogonalization is that hypothesis testing of a parameter of primary interest can be
done rather efficiently in the presence of the nuisance parameters. This is achieved by constructing a
class of functions (dependent on the null hypothesis) that is orthogonal to the scores of the nuisance
parameters. Neyman then considers a particular construction of regressing the scores of the parameter
of interest onto the scores of the nuisance parameters. Then the test, formulated by exploring this
orthogonality, can potentially lead a to a substantial increase in power.
In this section, we seek orthogonalization principles that are adapted to the presence of random
effects in a linear regression model. Suppose first that we observe independent samples (yi , X i , Z i )
and want to build a orthogonalization criterion suitable for testing (2) in a model (3). We start by
splitting the feature space until we have partitioned it into Z i , a set of covariates in the null, and
X i , associated with the nuisance parameters. Then, given observations (X i , Z i ) only, we evaluate
the regression of Z i onto the remaining covariates X i by setting
Z = Xθ ∗ + U

(4)

where θ ∗ ∈ Rp−1 and with components of U are i.i.d. with mean zero and unknown variance
0 < σu2 < ∞. Heuristically, this strategy is well-motivated if we believe the covariates of the fixed
design have a Gaussian distribution and that the rows are roughly identically distributed. There are
possibly several other procedures on how to best design this regression; but for simplicity we consider
a linear case – other more elaborate cases follow from our methodology with easy extensions of the
proofs.
5

We assume that U is independent of X,  and b. However, uncorrelatedness can replace the independence with little changes in the proof whenever the errors in the model take Gaussian structure.
Note that the coefficients θ ∗ are functions of the variances and covariances of the design matrices X
and Z, i.e. θ ∗ and of σu2 depend on X = (X, Z). They can be explicitly computed if the rows of X are
i.i.d. with normal distribution Np (0, Σ). In that case, denoting by j the column of X corresponding
to Z, it follows from Anderson (1984) that the conditional density f (Zi |X i ) is a (p−1)-variate normal
−1
−1
2
(p−1)×(p−1)
with mean X >
i Σ−j,−j Σ−j,j and variance σu = Σj,j − Σj,−j Σ−j,−j Σ−j,j , where Σ−j,−j ∈ R
is obtained by removing the jth row and column of Σ, and Σ−j,j ∈ Rp−1 is the jth column of Σ
∗
without its jth component. It follows that θ ∗ = Σ−1
−j,−j Σ−j,j . Table 1 gives the exact values of θ
and σu2 for the setting considered in our numerical work. More details are provided in Appendix.
Design

θ∗

Toeplitz

Sparse

Equi-correlated

Dense

Particular case
 ∗
θ = (0, . . . , 0, −0.4, −0.4, 0, . . . , 0)
ρ = −0.5 ⇒
2
 ∗σu = 0.6
θ = (0.002, . . . , 0.002)
ρ = 0.8 ⇒
σu2 = 0.2

Table 1: Exact values of θ ∗ and σu2 for the settings considered in our simulation study. Note that
for the Toeplitz case, if Z the first (resp. the last) column of X is tested then θ ∗ is slightly different;
only its first (resp. last) component is non-zero, with value ρ. Furthermore in that case, σu2 = 1 − ρ2 .
Finally, given regression (3) together with (4) our procedure for testing (2), generates the pseudo
response
V = Y − Zβ0
(5)

and observes that under the null hypothesis V − Xγ ∗ = W b +  and Z − Xθ ∗ = U are uncorrelated
and are such that
E[n−1 (V − Xγ ∗ )> A(Z − Xθ ∗ )] = 0
(6)

for a n × n positive definite matrix A. This moment matching equation can be seen as a specific
orthogonality condition where Z i are treated as confounders in the model on V. Conversely, we
observe that under the alternative hypothesis the above two terms are correlated since V − Xγ ∗ =
Z(β ∗ − β0 ) + W b + ; we compute that
E[n−1 (V − Xγ ∗ )> A(Z − Xθ ∗ )] = σu2 n−1 trace(A)(β ∗ − β0 ).
The advantage of the above orthogonalization, (6), is that the form of a moment above resembles
double-robust constructions where only one of the two residuals,  or U , needs to be estimated well
enough. Moreover, we will see that for a particular choice of the matrix A the above observation leads
to optimal inference for the fixed effects without requiring any knowledge or even consistent estimation of the random effects. Such approach is of interest on its own right (even for low-dimensional
problems) and would be extremely beneficial for practical purposes in high-dimensional longitudinal
studies where estimation of random effects is particularly difficult.

2.2

Estimation of the unknowns

In our discussion so far, we have emphasized the flexible nature of our methods: for a wide variety
of the structure of γ ∗ , distributions of  and of the random effects, Ψ, our methods can be tailored,
we achieve both consistency and asymptotic normality, provided the sample size n scales at an
appropriate rate. Our results do, however, require the features corresponding to the fixed-effects are
not correlated extremely highly: features are sparsely correlated if the precision matrix is row-sparse,
6

i.e, each row has small number of non-zero elements. We discuss a new class of estimators for both
γ ∗ and θ ∗ that satisfy this condition and are adaptive to the structure of the random effects.
Our first algorithm, which we call a doubly robust estimator, achieves adaptivity in estimation
of γ ∗ regardless of its structure and the consistent estimation of the random effects. It achieves
honest estimation by incorporating a pseudo response vector V , (5), and a proxy correlation matrix
e directly intro its construction. To motivate our estimation procedure, let us observe that if  would
P
be normally distributed, the underlying joint density of (V , b) would be
f (V , b) = f (V |b)f (b)

(−n+qN )/2

= (2πσ )

−1/2

|Ψ|




1
1 > −1
>
exp − 2 (V − Xγ − W b) (V − Xγ − W b) − b Ψ b .
2σ
2
(7)

We treat f (V , b) as a quasi-likelihood function; our method does not require normality of ε. For a
given γ, the maximum likelihood estimator of b is then

−1
b
b(γ) = W > W + σ2 Ψ−1
W > (V − Xγ).

We can plug-in this estimator into (7) to construct a profile likelihood for γ. We define E = W > W +
σ2 Ψ−1 and P = (I n − W E −1 W > )(I n − W E −1 W > ) + σ2 W E −1 Ψ−1 E −1 W > and, dropping the
constants, we obtain the following profile log-likelihood:
L(γ, b
b(γ)) = −

1
(V − Xγ)> P (V − Xγ).
2σ2

(8)

By Lemma 1, P = (I n + σ−2 W ΨW > )−1 . Note that if  is normally distributed then W b +  ∼
Nn (0, σ2 P −1 ). However, observe that the above profile log-likelihood is a non-convex function of all
of the unknown parameters, γ, Ψ, σ . Moreover, in the presence of high-dimensional fixed effects, the
profile log-likelihood has too many stationary points.
Following the idea of Fan and Li (2012), we replace the unknown variance matrix σ−2 Ψ by a
known fixed symmetric matrix M which serves as a proxy. We discuss in the next section how to
choose M ; our theory allows for many choices of M . We define proxy versions of E, P and L by
e = W > W + M −1 and Pe = (I n + W M W > )−1 and obtain the
replacing σ−2 Ψ by M ; we define E
following proxy log-likelihood:
1
e b
L(γ,
b(γ)) = − 2 (V − Xγ)> Pe (V − Xγ).
2σ

(9)

Observe that the proxy log-likelihood does not have a good quality approximation of the true L.
Instead it serves as a proxy to give inspiration for a new estimator of the unknown fixed effects γ ∗ .
Provided that the l1 penalization of the above proxy log-likelihood heavily depends on the correct
model specification and the form of the likelihood function, we take on a different perspective and
b as the solution to the
design a new estimator that is inspired by the Dantzig selector. We define γ
following optimization problem
b = argminγ∈Rp−1
γ
subject to

kγk1

n−1 X Pe (V − Xγ)
≤ ηγ
∞
> e
−1
n V P (V − Xγ) ≥ η̄γ
>

Pe (V − Xγ)

7

∞

≤ µγ

(10)

p
√
with ηγ  (log n) n−1 log p, µγ  log n and 0 < η̄γ < n−1 trace(σ2 Pe P −1 ) suitably chosen tuning
parameters. The l1 norm induces sparse solutions whereas the constraints ensure good theoretical
e
properties. The above constraints arise from the score vector of the profile-log likelihood function L;
the first ensures that the score vector is minimized and the second ensures that the variance of the
residuals is properly guessed. However, observe that we do not assume that the profile log likelihood
is correctly specified, and for that matter a last constraint is needed to guarantee small size of the
estimated residuals in the model
1/2
1/2
Pe V = Pe Xγ + e,
(11)

e P −1 . The tuning parameter η̄γ provides a proxy for the unknown variance of the
where var(e) = σ2 P
error in the model above. Hence, the size of the residuals is constrained in a similar way as robust
loss functions are truncated to down-weight the effect of overly large outliers.
b defined in (10) is new and of potential interest in its own right. We show
Remark 1. The estimator γ
that under suitable conditions, kb
γ − γ ∗ k1 = OP (ηγ kγ ∗ k0 ) meaning that this estimator is consistent if
√
∗
kγ k0 = o( n/ log(p)/ log(n)). It is worth pointing that this result does not impose any restrictions
e correctly specifies variance parameters of the
on the minimum signal strength or that the matrix P
random effects and in that sense is honest and robust. Beside providing a reliable estimator, the
optimization problem (10) is fast since it can be reformulated as a linear program. Existing regularized
schemes, such as the one introduced by Schelldorfer et al. (2011) for example, do not provide such
quick implementation.
The estimation of θ ∗ is based on model (4). Unlike linear models, where node-wise lasso is sufficient
to estimate feature correlation, for the case of linear mixed models, we observe that the matrix Pe
should contribute to the estimation procedure – we see from (11) that the covariates are premultiplied
1/2
with Pe . Namely, the presence of random effects induces a larger than usual dependence in the
design matrix of the fixed effects. We incorporate such requirement in a new estimator defined below.
θb = argminθ∈Rp−1
subject to

kθk1
− Xθ)

n−1 X > (Z

n−1 X > Pe (Z
n−1 Z > (Z

− Xθ)

∞

≤ ηθ

∞

≤ ηθ0

(12)

− Xθ) ≥ η̄θ
kZ − Xθk∞ ≤ µθ
p
√
with tuning parameters ηθ , ηθ0  (log n) n−1 log p, µθ  log n and 0 < η̄θ < σu2 .
The first and the second constraint enable adaptive and automatic estimation of the correlation
of the features of the fixed effects without prior knowledge (or correct estimation) of the correlation
of the random effects. In particular, the first ensures that the gradient of the square loss is close to
zero, whereas the second one looks at the reweighed gradient where the weights are defined through
e . The presence of the second constraint allows us additional flexibility with the choice of
the matrix P
e . We can remove the second constraint if we impose that λmax (W > M W ) is bounded;
the matrix P
however, this would restrict our choices of the matrix M and in particular it would not allow for
M = log(n)In which we found to be beneficial whenever the variance in the random effects is
particularly large. Lastly, the last constraint excludes features X that are highly correlated with Z;
as our task was to remove heterogeneity (by de-correlating features), such reasoning is needed for a
successful test statistic.

2.3

Asymptotic inference for linear mixed models

Our results are achieved under the asymptotics where N → ∞ and ni /N → 0, a regime different
from that used in simple linear models, and require very mild conditions on the choice of the matrix
8

e , i.e., the matrix M , which needs to have the same diagonal structure as W W > . However, given
P
these high level conditions, we obtain a widely applicable result that applies to several different linear
mixed models.
In order to define a test statistic, we focus on the moment condition (6) and in particular consider
e . That is, the moment condition that we wish to test is now taking the form E[n−1 (V −
A = P
Xγ ∗ )> Pe (Z − Xθ ∗ )] = 0 whenever the null hypothesis H0 hold. For the alternative hypothesis, this
moment takes the form E[n−1 (V − Xγ ∗ )> Pe (Z − Xθ ∗ )] = σu2 n−1 trace(Pe )(β ∗ − β0 ).
We then proceed to form a test statistic. The moment above produces a doubly robust test
statistic as long as one of two component regression models is correctly specified and assuming that
there are no unmeasured confounders, giving the analyst two chances to correctly specify at least
one of the regression models, (3) and (4), respectively. In order to standardize the test statistic
b 2 . Whenever the models
b )k22 and σ
appropriately, we define σ
b2 = n−1 kPe (V − X γ
bu2 = n−1 k(Z − X θ)k
2
are correctly specified, they produce residual variance estimates. However, the quality of testing does
not rely on them being consistent estimators; in fact, as the proxy log-likelihood is not exact, this
will not be achievable; we show in Lemma 7 that σ
b2 = σ2 n−1 trace(Pe P −1 Pe ) + oP (1).
Now, we proceed to define the test statistic as
Tn =

b > Pe (V − X γ
b)
n−1/2 (Z − X θ)
.
σ
bu σ
b

(13)

The doubly robust test statistic combines the initial model under the null hypothesis with a regression model of the relationship between covariates and the feature related to each of the parameters
of interest in such a way that, as long as either the initial model or the feature regression model is
correctly specified, the effect of the initial null hypothesis on the transformed moment is correctly
estimated. In the case of linear mixed models, the initial model is only correctly specified when the
variance parameters Ψ are known, making the test statistic above particular useful for cases when
there is no such knowledge.
Our first result is that the test statistic has asymptotically a standard normal distribution whenever the null hypothesis holds (see Theorem 3). For this result we do not require sparsity of the
fixed effects, but rather sparsity of the precision matrix of the design matrix of the fixed effects.
However, we do not require any consistent estimation of the covariance parameters Ψ. To obtain
the asymptotic distribution under the alternative hypothesis, H1 : β ∗ = β0 + n−1/2 h, we need to
assume that the vector of nuisance parameters γ ∗ is sparse (or approximately) with a small number
of non-zero elements, although our simulations show that the power is preserved even if the fixed
effects have as many as n non-zero elements (see Theorem 4). To our knowledge, all existing results
regarding inference of linear mixed models require consistent estimation of the covariance parameters
and low-dimensionality of the fixed effects vector. Although some empirical work has demonstrated
that inconsistent variance estimation does not effect type I error control (in low-dimensional settings),
theoretical guarantees were never established.

3

Asymptotic Theory for Linear Mixed Models

In order to use the test statistics to provide formally valid statistical inference, we need an asymptotic
normality theory in the desired asymptotics. We begin by precisely describing the asymptotics under
the null hypothesis. We then proceed to describe the asymptotics under the alternative hypothesis
H1 : β ∗ = β0 + hn−1/2 .
Before stating our theoretical results, we give some definitions.

9

(14)

Definition 1. We say that a symmetric semi-definite positive matrix A satisfies the P -condition if
it has the same block diagonal structure as W > W and if each element of A is O(log(n)).
In Lemma 4, we show that Pe = (I n + W M W > )−1 satisfies the P -condition for any matrix M
having the same block diagonal structure as W > W . Furthermore, it holds that P and P −1 satisfy
the P -condition.
Next, we require that the errors  in our model have sub-Gaussian tails. A test statistic, as defined
in Section 2.3, can be used to give honest p-values or confidence intervals; in Section 4 we illustrate
its good properties even when the error is not sub-Gaussian but exhibits heavy-tails.
Definition 2. A random variable X is said to have an exponential-type tail with parameters (b, γ) if
∀x > 0, P(|X| > x) ≤ exp[1 − (x/b)γ )]. Furthermore, a random variable X is sub-Gaussian if it has
an exponential-type tail with parameters (b, 2).
In order to guarantee consistency, we also need to enforce that the design matrix of the fixed effects
satisfies some regularity conditions. The conditions are standard in the high-dimensional literature
and ensure that the population feature covariance matrix is a well conditioned matrix. Here, we
follow Rudelson and Zhou (2013), and achieve this effect by enforcing sub-Gaussianity in the design
in the following condition. Let Cmin , Cmax and κ denote positive constants.
Condition 1. The matrix Σ = E[X > X]/n ∈ R(p−1)×(p−1) is such that Cmin ≤ σmin (Σ) ≤ σmax (Σ) ≤
Cmax with σmin (Σ) and σmax (Σ) the minimal and maximal singular values of Σ. The vectors Σ−1/2 xi ,
 and u are centered with sub-Gaussian norms upper bounded by κ. Moreover, σu , σ ∈ [Cmin , Cmax ].
The elements of the symmetric and invertible matrix Ψ and of the deterministic design matrix W
are bounded.
The remaining definition is more technical. We use a regularity condition to control the shape
of the correlation in the features and in the variance covariance matrix Ψ. Such condition is used
regularly in high-dimensional estimation and inference (see Bickel et al., 2009).
Definition 3. We say that the restricted eigenvalue condition holds for a triplet (s, κ, X) if
min

J0 ⊆{1,...,p−1}
|J0 |≤s

min
δ J0c

δ6=0
≤kδ J0 k1

kXδk2
√
≥ κ.
n kδ J0 k2

1

Rudelson and Zhou (2013) show that the restricted eigenvalue condition holds with large probability if the sample size n is large enough and the rows of the design have sub-Gaussian distributions.
We note that our test statistic (13) does not depend on the inverse of Ψ. So although we started
this section with the nonsingularity assumption of Ψ, in practice our method can be directly applied
even when noise random effects exist.
Given these preliminaries, we state our main result on the asymptotic normality of the test statistic
√
Tn . As discussed previously, we require that the sample size scales as n/ log(p) = o(1) as N → ∞.
If the subsample size grows slower than this, the test statistic will still be asymptotically normal, but
the it may be asymptotically biased. Although we treat n, s and p as functions of N , for clarity, we
state the following result without this explicit dependence.
Theorem 1. Suppose that we have N independent observations (yi , X i , Z i , Wi ) ∈ Rni × Rni ×(p−1) ×
e be such that
Rni × Rn1 ×q . Suppose moreover that the features satisfy Condition 1. Let the matrix P
b
b and θ be estimators of the nuisance
Definition 1 holds. Given this data-generating process, let γ
parameters and the feature dependence as defined in (10) and (12). Finally, suppose that the sample
10

size scales as n1/2 / log(p)/ log(n) = o(1) and maxi ni /N = o(1) and that the models are normalized
so that kγ ∗ k2 = O(1) and kθ ∗ k2 = O(1). Then,
d

Tn = N (0, 1) +

√

n−1 trace(Pe )
n(β ∗ − β0 ) q
+ oP (rn ),
−1 e
−1
e
n trace(P P P )

where rn = n1/2 ηθ ηγ kθ ∗ k0 κ−2 η̄γ−1 σ /σu whenever β ∗ = β0 and
p
√
rn = n1/2 ηγ kγ ∗ k0 κ−2 ηθ0 σ /σu + ηθ kθ ∗ k0 σ /σu + n(β ∗ − β0 )N/nσ /σu whenever β ∗ 6= β0 .

The construction of p-values is based on the asymptotic distribution of the test statistic Tn . For
the null hypothesis H0 : β ∗ = β0 , we define the p-value for the two-sided alternative as P0 = 2(1 −
Φ(|Tn (β0 |))). Of course, we could also consider one-sided alternatives with an obvious modification.
Whenever the conditions of the Theorem 1 hold, then for any 0 < α < 1,
lim sup P [P0 ≤ α] = α if H0 holds .
n→∞

Furthermore, for any sequence an → 0 ( n → ∞) which converges sufficiently slowly, the statements
also hold when replacing α by an . A discussion about detection power of the method is given in
Section 3.
The proof of Theorem 1 is organized as follows. In Section 3.2, we provide bounds for the l1 norm
error in estimation of both γ ∗ and θ ∗ , while Section 3.3.1 studies the sampling distributions of the test
statistic under the null hypothesis and establishes asymptotic Gaussianity. Given a subsampling rate
and sparsity requirements in the two models, (3) and (4), we showcase the asymptotic distribution
under the alternative hypothesis in Section 3.3.2. Before beginning the proof, however, we relate our
result to existing results about linear mixed models in Section 3.1.

3.1

Theoretical background

There has been considerable work in understanding the theoretical properties of linear mixed models. The convergence and consistency properties of estimation have been studied by, among others,
McCulloch (1997); McGilchrist (1994); Lindstrom and Bates (1988) in low-dimensional setting, Goeman et al. (2011) with growing dimensions and Schelldorfer et al. (2011) in high-dimensional setting.
Meanwhile, the inference has been analyzed by Breslow and Clayton (1993) in low-dimensional setting. However, to our knowledge, our Theorem 1 is the first result establishing conditions under
which inference in linear mixed models can be performed despite the presence of a large number of
fixed effects and taking into account model selection. Moreover, Theorem 1 establishes this result
without even resorting to a proper (let alone consistent) estimation of the random effects therefore
enabling inference in linear models with misspecified random effects (the structure or distribution);
we believe that such result is unique in inferential theory of linear mixed models.
Probably the closest existing result is that of Fan and Li (2012), who showed that selection of
fixed effects can be achieved without resorting to proper estimation of random effects. They establish
model selection consistency results under conditions similar to ours; however, we note that they
e whereas we allow more general choices – for
impose a slightly more restrictive choice of the matrix P
example choice of M as an identity matrix is not allowed in their work, but it is in ours. Moreover,
the authors therein do not discuss asymptotic distributions and hypothesis testing problems. Thus,
their results cannot be used for valid asymptotic statistical inference about fixed effects.
Inference for linear mixed models that is robust to the normality assumption in random effects
or the model error has been studied exclusively in low-dimensional setting in Zhang and Davidian
(2001), who showed that normality can be substituted with a broader class of distributions and an
11

EM algorithm can be used for estimation. However, their results still heavily depend on a correct
likelihood specification. Identical conclusions hold for a class of semi-parametric linear mixed models
introduced in Lachos et al. (2010) or a linear mixed model with random mixture of normals (Verbeke
and Lesaffre, 1996; Song et al., 2007). See Heagerty and Kurland (2001) for an illustration of the
effects of inconsistent variance estimation for the likelihood based inference on fixed effects. Additional
work on this topic Litière et al. (2007), Hobert and Casella (1996) had only highlighted connections
to statistical tests and the effects of misspecification of the random effects, but did not establish any
formal justification for it.

3.2

Estimation properties

We start by bounding the bias of the newly introduced high-dimensional estimators. Our approach
relies on showing that the true parameter vector lies in the constraint set of the problems (10) and
(12) for appropriate choices of the tuning parameters defined therein. A sparsity assumption then
allows us to bound the bias. We show that the error of estimation is optimal and is not effected by
the misspecification of the random effects. In order to state the result, define kγ ∗ k0 as the number
of non-zero elements of a vector γ ∗ .
Theorem 2. Suppose Pe ∈ Rn×n satisfies the P -condition. Suppose that Condition 1 holds and that
∗
e 1/2 X). Then
there exists κ > 0 such that the restricted eigenvalue condition holds
p for (kγ k0 , κ, P √
there exist constants C1 , C2 , c0 > 0 such that for any ηγ ≥ C1 log(n) n−1 log p, µγ ≥ C2 log n and
η̄γ , ∈ (0, n−1 trace(σ2 Pe P −1 ) − c0 ), for a large enough constant C3 it holds with probability 1 − p−C3
kb
γ − γ ∗ k1 ≤ 8ηγ kγ ∗ k0 κ−2 , and

1/2
n−1 Pe X(b
γ − γ ∗)

2
2

≤ 16ηγ2 kγ ∗ k0 κ−2

(15)

.

(16)

This theorem
then directly
on the bias in estimation. Namely, with
p
√ translates into a bound
−1
−1
b of a s sparse
ηγ  log(n) n log p, µγ  log n and 0 < η̄γ < n trace(σ2 Pe P −1 ) the estimator γ
vector γ ∗ ,
p
kb
γ − γ ∗ k1 = OP (s log(n) log p/n).

The conditions required for the above result seem very mild. Finite sample oracle risk properties
have been established in Schelldorfer et al. (2011); however, the result therein crucially depends on
correct estimation of variance components. In contrast, our results hold for a wide variety of choices
of variance estimators; in particular, it remains valid for even non-consistent estimators – the only
assumption is on the structure of the matrix P illustrated through the P-condition. Furthermore, the
rate above matches the rate obtained by the mle estimator of Schelldorfer et al. (2011) highlighting
b and indicating optimality at esitmation.
excellent robustness properties of the proposed estimator γ

3.3

Testing properties

Our analysis proceeds in two steps: the first discusses asymptotic normality of the test statistic
whenever the null hypothesis holds, whereas the second discusses the case under the alternative
hypothesis.
3.3.1

Size

Analyzing specific linear mixed models can be challenging especially if the model (1) is not fully or
correctly specified. Below we establish asymptotic normality of the proposed test statistic. The result
is independent of the size of the sparsity of the nuisance parameters and in particular holds for fully
12

dense and ultra high-dimensional vectors. Moreover, it does not require consistent estimation of the
variance components.
Theorem 3. Suppose Pe ∈ Rn×n satisfies the P -condition. Suppose that Condition 1 holds and that
1/2
there exists κ > 0 such that the restricted eigenvalue condition holds for (kθ ∗ k0 , κ, Pe X). Further√
more, assume that kθ ∗ k0 = o( n/ log(p)/ log(n)) and a normalization of the form kγ ∗ k2 = O(1) and
d

kθ ∗ k2 = O(1). Then, under the null hypothesis H0 : β∗ = β0 , it holds that Tn → N (0, 1).

Theorem 3 establishes the asymptotic distribution of our test statistics under the null hypothesis;
it does not require minimal signal strength in the model (1) or an irrepresentable condition. Notably,
our statistic does not require γ ∗ to be sparse. This is a strong and unusual result in high-dimensions.
This finding is further illustrated numerically in Section 4 where we observe stable control over the
type I error rate regardless of the size of kγ ∗ k0 .
Furthermore, using some arguments of the proof of Theorem 4, interchangeability of sparsity
conditions can be shown: if γ ∗ is sparse then θ ∗ does not need to be sparse. Note that the
√
o( n/ log(p)/ log(n)) sparsity rate is up to a log(n) term matching those of simple linear models, see
for instance van de Geer et al. (2014) and Javanmard and Montanari (2014b). The additional log(n)
term is needed in controlling the random effects and can be though of as a price to pay for being able
to provide optimal estimation despite the correct specification of the random effects.
Besides providing the size property of the Tn statistic, Theorem 3 can also be used to construct
confidence intervals. The 1 − α confidence interval for β ∗ can be defined as

β : 1 − Φ−1 (1 − α/2) ≤ Tn (β) ≤ Φ−1 (1 − α/2)
where Φ is the standard Gaussian cumulative distribution function and where Tn (β) is the test
statistic derived under the null hypothesis β ∗ = β.

3.3.2

Power

In the previous section, we showed that the type I error is asymptotically equal to α for a given level
α ∈ (0, 1). In this section, we show that the test statistic furthermore has tight control of the type
II error and preserves power asymptotically while allowing inconsistent variance estimation. In this
regard we believe that our result stands out in the existing literature. In addition, we quantify the
asymptotic efficiency loss due to misspecification of the random effects.
Theorem 4. Suppose Pe ∈ Rn×n satisfies the P -condition, Condition 1 holds and that there exists

1/2
κ > 0 such that the restricted eigenvalue condition holds for (kθ ∗ k0 ∨ kγ ∗ k0 , κ, Pe X). Furthermore,
√
√
assume that both kθ ∗ k0 = o( n/ log(p)/ log(n)) and kγ ∗ k0 = o( n/ log(p)/ log(n)) and that the
normalization kγ ∗ k2 = O(1) and kθ ∗ k2 = O(1) holds. Then, under the alternative hypothesis H1 in
(14), it holds that with n → ∞
d

Tn = N (0, 1) + hσu σ−1 q

n−1 trace(Pe )

n−1 trace(Pe P −1 Pe )

+ oP (1).

Theorem 4 establishes the power properties of our test and requires both γ ∗ and θ ∗ to be sparse
√
with a o( n/ log(p)/ log(n)) sparsity rate. The deviation term depends on the choice of the matrix
Pe and by the Cauchy-Schwarz inequality it can be shown that
p
n−1 trace(Pe )
q
≤ n−1 trace(P ).
n−1 trace(Pe P −1 Pe )
13

This implies that using the proxy matrix Pe instead of the true unknown matrix P leads to an
asymptotic loss of power.
If we consider the particular, naive choice of the proxy matrix M = 0N q×N q then Pe reduces to
the identity matrix I N q . In that case the method reduces to the method entirely based on simple
linear model and leads to a significant loss of power. For a simple choice of the matrix M = cI N q
where c is a constant positive number, the matrix Pe takes the form (I n + cW W > )−1 .
Using Theorem 4 we can explicitly derive the expression for power of the Tn statistics under
nominal level α as




α
α
Power = 2 − Φ Φ−1 (1 − ) − Dn (h) + Φ Φ−1 (1 − ) + Dn (h)
2
2
where Dn (h) = hσu σ−1 q

e)
n−1 trace(P

.
Pe )
Furthermore, the length of the 1−α confidence interval of β ∗ constructed from Theorem 3 is given
by 2n−1/2 |hα | for hα satisfying 2Φ Φ−1 (1 − α/2) − Dn (hα ) = 1 − α which can easily be estimated
numerically.
We also note that the result can be extend to include fully dense models where kγ ∗ k∞ is not
allowed to grow rapidly with n but the l0 norm can be as large as p; or mixture settings with many
small elements and only a small number of strong elements (here the sum of the small elements would
be allowed to grow with p). We conjecture that the power function will look the same as it does in
Theorem 4.

4

eP
n−1 trace(P

−1

Numerical Experiments

In observational studies, accurate construction of p-values requires to overcome three potential sources
of bias. First, there is the initial linear mixed model and in particular the size of the sparsity of the
nuisance parameters γ ∗ . We need to make sure that our test is reasonably stable whatever the
structure of γ ∗ is. Then, second, there is the underlying precision matrix of the design matrix. Here
a bias can be introduced with presence of large feature correlations and/or dependencies. Third, we
need to make sure that the test is reasonably stable regardless of the exact structure of the random
effects; we need to make sure that the test is not biased by changing matrix Ψ, i.e. ψ. The simulations
here aim to test the ability of the introduced test to respond to all three of these factors.

4.1

Models and designs

We consider the mixed model

Y = Xβ ∗ + W b + 

with n = 200 observations coming from N = 50 different groups of size ni = 4. The vector β ∗ of
fixed effects is of length p = 500. The components of  are independent and come from a standard
> >
qN with b ∼
Gaussian distribution. The random effect vector is defined as b = (b>
i
1 , . . . , bN ) ∈ R
n×qN
Nq (0, ψ) for i = 1, . . . , N . We define W = diag(W 1 , . . . , W N ) ∈ R
with the stack matrix
> >
n×q consisting of the first q columns of X. We consider three models:
(W >
1 , . . . , WN) ∈ R
• Model 1 : the rows of X come from a multivariate normal distribution with zero mean and
Toeplitz covariance matrix Σ with Σij = (−0.5)|i−j| . The components of β ∗ are generated as
follows: β ∗ = 5a/ kak2 with a = (a1 , . . . , ap ), where aj is generated from a uniform distribution
on (0, 1) if j ≤ 3s/2 and j/3 is not an integer; otherwise, aj = 0. We then have kβk0 = s. The
model contains q = 2 random effects with covariance matrix ψ = diag(0.56, 0.56).
14

• Model 2 : the setting is the same as in Model 1 except that q = 3 and ψ = diag(3, 3, 2).
• Model 3 : the setting is the same as in Model 1 except for Σ which is such that Σij = 1 if i = j,
Σij = −ρ/(1 + ρ2 ) if |i − j| = 1, and Σij = 0 otherwise, where ρ varies.
As noted in the appendix, the Toeplitz design implies that the underlying correlation model (4) is
sparse with two non-zero components being equal (for ρ = −0.5) to −0.4. Thus, for models 1 and 2,
the vector θ ∗ is sparse. Conversely, model 3 is such that that the rows of Σ are sparse, while the row
of its inverse are not sparse sparse. This in turn implies that θ ∗ in (4) is a dense vector. Regarding
the random effect part, similar choices of q and ψ are made by Schelldorfer et al. (2011).

4.2

Tuning parameters

Tuning parameters need to be chosen in optimization problems (10) and (12). Based on our investigation, we suggest the following choices. The sensitivity of our procedure to the choice of the tuning
parameters is given in Table 4. Similarly to Zhu and Bradic (2016), we choose η̄γ = 0.05 V > Pe V /n.
p
√
We define ηγ = 0.5n−1 log p σ
b and µγ = 4 log n σ
b with σ
b = Pe (V − Xγ init ) where γ init is
2

obtained by a linear estimation (without the random effects) with the scaled lasso. Such an initial
estimator of γ is also used by Rohart et al. (2014). The choice of the tuning parameters η̄θ , ηθ , µη
and ηθ0 is done in a similar way.
Regarding the choice of the matrix Pe = (I n + W M W > )−1 serving as a proxy of P = (I n +
−2
σ W ΨW > )−1 , we show in Lemma 4 that there is a large pool of possible choices guaranteeing good
asymptotic results. Fan and Li (2012) suggest to use M = (log n)I N q . In the present simulation
−2
study we use the slightly more elaborate M = σ,init
diag(ψ init , . . . , ψ init ) where initial estimators are
derived similarly as in Rohart et al. (2014): using γ init obtained by a linear model estimation with
2[−1]
2[−1]
2[−1]
1
2
kV − Xγ init k22 , ψinit = 0.4
I q and σ,init
= 0.6σ
.
the scaled lasso, we compute σ
= n−df
q σ

4.3

Numerical results

We test H0,j : βj∗ = β0,j for j = 3, 4, 100 and study the behaviour of our procedure. These values of
j correspond to different types of coefficients; components 3 and 100 are not active while component
4 is active. Furthermore, in the Toeplitz case, component 100 is almost not correlated with any
of the active variables. We set H1,j : βj∗ = β0,j + hn−1/2 where the situation h = 0 corresponds
to size properties while h 6= 0 corresponds to power properties. All our results are based on 1000
replications and on a 5% nominal level. To the best of our knowledge, our procedure is the first one
to do hypothesis testing for high-dimensional mixed models which makes it hard to find competitors.
We compare our mixed model method to the linear model method of Zhu and Bradic (2016).
4.3.1

Gaussian mixed model

In Figure 1, we report size and power properties of our procedure for Model 1 and for different
components of β ∗ . On Figure 1 (a) we observe that the false discovery rate is close to the 5%
nominal level whatever the sparsity level. Our procedure can handle dense β ∗ and this is one of its
major features. Sparsity is required only for either θ ∗ or β ∗ to obtain reliable results under the null
hypothesis. Note that under the alternative hypothesis h 6= 0, our theoretical results require both of
them to be sparse, see Theorem 4. On Figure 1 (b), we provide the power of our test for Model 1 in
a sparse model (s = 5). Interestingly, the probability to reject the null hypothesis is not symmetric
in h. This can be explained by the fact that we use a penalized estimator which shrinks coefficients
towards zero. The correlation between the components compensates (here for h negative) or amplifies
(here for h positive) the bias of the estimator, which leads to an asymmetry of the results.
15

(b) Probability to reject H0

0.4
0.0

0.00

0.05

0.8

0.10

(a) Probability to reject the true H0

20

40

60

80

−6

100

−4

−2

0

2

4

6

deviation from the null

sparsity

Figure 1: Empirical rate of rejection of the null hypothesis under (a) the null hypothesis, and (b)
alternative hypothesis, at a 5% nominal level for Model 1. As described in the text, hypothesis
tests are performed on different components of β ∗ : the third component (black crosses), the fourth
component (red circles), and the 100-th component (blue triangles).

0.00

0.10

0.20

0.30

In Figure 2 we give size and power properties of our procedure for Model 3 in which θ ∗ is dense, i.e.
sθ = p in this particular case. We vary the size of the correlation ρ with larger values corresponding
to larger size of the coefficients of θ ∗ . We observe that the proposed test is remarkably resilient
and controls Type I error in finite samples. We observe that when γ ∗ is not sparse, which here
√
corresponds to sparsity larger than 5 (note that n/ log(p) ≈ 5 for n = 200, p = 500), our method
is not guaranteed to control Type I error; however, we see that our method strikingly controls errors
even in such cases whenever ρ < 0.5. It is worth pointing that no method is expected to perform
well when all elements are dense and correlation is high. Hence, we take this behavior to be close to
optimal.

0.10

0.20

0.30

0.40

rho

0.50

0.60

0.70

0.80

Figure 2: Empirical rate of rejection of the true null hypothesis at a 5% nominal level for Model 3
with dense vector θ ∗ . We vary the level of sparsity of γ ∗ and correlation level ρ. Hypothesis testing is
performed on the second component of β ∗ . Empirical Type I errors are plotted for different sparsity
levels: black crosses 2, red circles 3, blue triangles 4, orange diamonds 5, green upside-down triangles
10.
In Table 2, we compare our method, for model 1 and 2, to the method of Zhu and Bradic (2016)
designed for linear models and which thus does not take into account the group structure of the mixed
model. For both models, the probability to reject the null hypothesis is closer to the 5% nominal
level when the null hypothesis holds (h = 0) and is larger under the alternative hypothesis (h 6= 0),
compared to the linear model procedure.
16

h = −6
h = −4
h = −2
h=0
h=2
h=4
h=6

Probability to reject the null hypothesis at a 5% level
Model 1
Model 2
LM procedure MM procedure LM procedure MM procedure
0.85
0.96
0.46
0.65
0.60
0.77
0.30
0.39
0.29
0.30
0.19
0.16
0.11
0.06
0.12
0.05
0.11
0.16
0.07
0.05
0.28
0.56
0.07
0.18
0.55
0.89
0.13
0.42

Table 2: Probability of rejecting the null hypothesis as a function of the deviation from the null.
Nominal level is taken to be 5%. LM procedure consists in applying the linear model procedure of
Zhu and Bradic (2016) to the mixed model while MM procedure consists in applying our mixed model
procedure. Parameters are n = 200, p = 500, N = 50, s = 5. The 4th component of β ∗ is tested.
Model 1 contain 2 random effects and Model 2 contains 3 random effects with large magnitude.
4.3.2

Model with heavy-tailed design

In the next example we consider a model that departs from normality assumptions. Parameters
choices are done in the same way as in Model 1: n = 200, N = 50, p = 500, and q = 2.
• Model 4 : the setting is the same as in Model 1 except that the entries of Σ−1/2 xi and of  are
generated from a student t distribution with 10 degrees of freedom.
• Model 5 : the setting is the same as in Model 1 except that the entries of Σ−1/2 xi and of  are
generated from a student t distribution with 3 degrees of freedom.
We perform hypothesis test for the fourth component of β ∗ for Models 1, 4, and 5 and report
empirical sizes and powers in Figure 3. We observe empirical coverages close to the 5% nominal level
for the three models under the null hypothesis whatever the sparsity level is. We also observe that
the presence of heavier tailed errors results in a decrease of power.
(b) Probability to reject H0

0.0

0.03

0.4

0.05

0.8

0.07

(a) Probability to reject the true H0

20

40

60

80

−6

100

−4

−2

0

2

4

6

deviation from the null

sparsity

Figure 3: Empirical rate of rejection of the null hypothesis under (a) the null hypothesis, and (b)
alternative hypothesis, at a 5% nominal level for different distributions of the error: Gaussian (model
1, black crosses), Student with 10 degrees of freedom (model 4, red circles), and Student with 3
degrees of freedom (model 5, blue triangles). The hypothesis test is the H0,4 : β4∗ = β0,4 versus
H1,4 : β4∗ = β0,4 + hn−1/2 .
17

4.3.3

Choice of the proxy matrix

The choice of the proxy matrix Pe = (I n + W M W > )−1 has an influence on the power of the test.
Our default choice, presented in section 4.2 and inspired by Rohart et al. (2014) consists in working
−2
with M = σ,init
diag(ψ init , . . . , ψ init ). With these initial estimators, this is equivalent to work with
2
M = 3q
I N q . A second choice, proposed by Fan and Li (2012), is to work with M = log nI N q .
Simulation results for those choices of the proxy matrix are reported in Table 3 for models 1 and
2. We first observe that the two choices provide the 5% nominal level under the null hypothesis.
Secondly, we observe slight differences in terms of power which can be explained by how well the true
variance structure is approximated by the proxy matrix M .

h = −6
h = −4
h = −2
h=0
h=2
h=4
h=6

Probability to reject the null hypothesis at a 5% level
Model 1
Model 2
−2
−2
M = σ,init ψ init M = log nI N q M = σ,init ψ init M = log nI N q
0.96
0.91
0.65
0.79
0.77
0.66
0.39
0.49
0.30
0.25
0.16
0.19
0.06
0.05
0.05
0.05
0.16
0.15
0.05
0.09
0.56
0.49
0.18
0.30
0.89
0.85
0.42
0.64

Table 3: Probability to reject the null hypothesis in function of the deviation from the null for a 5%
nominal level, and in function of the choice of the proxy matrix. Parameters are n = 200, p = 500,
N = 50, s = 5. The 4th component of β ∗ is tested. Model 1 contains 2 random effects and Model 2
contains 3 random effects with large magnitude.

4.3.4

Sensitivity to tuning parameters choice

We conclude this simulation study by an analysis of the sensitivity of our procedure to the choice of
the tuning parameters. We consider Model 1 and two scenarios: (h, s) = (0, 40) and (h, s) = (4, 5)
corresponding to a non-sparse model under H0 and an sparse model under the true alternative H1 .
We perform hypothesis tests for the fourth component of β and provide results for different choices
of the tuning parameters. Results of Table 4 illustrate that the method is reasonably weakly sensitive
to the tuning parameters choice.

0.5µ
µ
2µ

Probability to reject H0 at a 5% level
h = 0 and s = 40
h = 4 and s = 5
1
1
η
2η
η
2η
2η
2η
0.06 0.05 0.05 0.48 0.48 0.53
0.04 0.05 0.05 0.59 0.59 0.59
0.06 0.05 0.06 0.35 0.34 0.34

1
2 η̄

η̄
2η̄

Probability to reject H0 at a 5% level
h = 0 and s = 40
h = 4 and s = 5
0.04
0.60
0.05
0.59
0.04
0.42

Table 4: Probability to reject the null hypothesis for different choices of the tuning parameters and
for Model 1 with j = 4 and (h, s) = (0, 40) or (h, s) = (4, 5), compared with the default choice
presented in Section 4.2. Notation η refers to ηγ , ηθ and ηθ0 ; µ refers to µγ and µθ ; η̄ refers to η̄γ and
η̄θ . If it is not specified, the default choice is used.

18

4.4

Hypothesis testing for Riboflavin data

We study the Riboflavin data which contains 4088 gene expressions from 111 observations divided
in 28 groups from size 2 to 6, and which has also been considered by Schelldorfer et al. (2011). The
response variable is the logarithm of the riboflavin production rate of Bacillus subtilis. In the context
of linear models a related dataset was considered in van de Geer et al. (2014) and Javanmard and
Montanari (2014b). The data of our interest has a clear group structure, and therefore we include a
random intercept in the linear model, resulting into the following mixed model
Y = Xβ ∗ + W b + 
where W contains 1 in the appropriates entries, Y is of length 111 and X is the 111 × 4089 design matrix with its first column containing only ones and its other columns corresponding to the
4088 standardized covariates. The vector b contains the 28 realizations of the random effects. Its
components are independent and identically distributed with mean zero and unknown variance ψ.
Riboflavin (vitamin B2) is an essential component of the basic metabolism. The riboflavin biosynthesis in bacteria was analyzed from a biological perspective using comparative analysis of genes,
operons and regulatory elements. However, little is known about the mechanisms of regulation of the
bacterial riboflavin genes.
Gene YpaA has been verified experimentally to be a transporter of riboflavin or related compounds, co-regulated with other riboflavin genes (Vitreschak et al., 2002). The RFN element (specific locus related to RNA folding) was only encoded upstream of the YpaA gene. Moreover, in B.
subtilis, riboflavin uptake was increased when YpaA was over-expressed and abolished when YpaA
was deleted (Vogl et al., 2007). Hence it makes sense to test a one sided hypothesis with the null
being that the corresponding β ∗ is smaller or equal to zero and alternative that it is larger than zero.
In this context, we have applied our test above while including all of the 4087 remaining genes in the
study. We have obtained a test statistic value of 1.60. Therefore our method is able to identify YpaA
gene at a 6% significance level; we should note that our sample size is extremely small and deviations
from 5% should be expected due to an extremely small sample size. It is also worth noting, that no
previous study of this dataset was able to identify YpaA as statistically significant gene. We also
have observed large correlation among the genes present in this study, which may additionally explain
why previous methods failed; observe that our method is doubly robust and is able to overcome high
correlations in the data.
The best studied system of the riboflavin biosynthesis in bacteria is the rib operon of Bacillus
subtilis. Our data contains information about the Bacillus subtilis and there are a number of genes
(but not all) related to rib operon. We are particularly interested in the ribB gene which was reported
as over-expressed in a number of chemical and biological studies (Mörtl et al., 1996) but has not been
yet studied from a statistical point of view. We have performed a single hypothesis test related to the
β ∗ corresponding to ribB gene; while including all of the remaining genes. We have found that our
test statistic is able to detect ribB with a 1% nominal value; the observed value of the test statistics
is 2.69. This result confirms the biological evidence and reconfirms wide-applicability of the proposed
test.
Another part of the rib operon is a ribC gene. The ribC gene was cloned and sequenced and it
was determined that it plays an essential role in the flavin metabolism of B. subtilis (Mack et al.,
1998). In particular it was observed that it is over-expressed and that it suppresses the riboflavin
overproduction. We have performed an one-sided hypothesis test to observe statistical significance
of this gene in the current dataset. We have observed a test statistic value of 1.01 leading to a 16%
significance level. Although larger than 10%, this study is one of the first statistical studies that was
able to detect ribC gene at any level. For completeness, we also report that an averaged size of the
selected sets in the original and the feature model was 17.90 and 741.10.
19

5

Further explorations

The methodology presented in the previous sections is extremely general and provides easy extensions
to a number of interesting problems and settings.

5.1

Power improvements

Our test statistic uses a proxy matrix Pe which has some influence on the power through the deviation
p
n−1 trace(Pe )
term hσu σ−1 √ −1
in
Theorem
4.
This
deviation
term
is
smaller
than
n−1 trace(P ) and
e −1 e
n

trace(P P

P)

it is expected that the “closer” Pe is to P , the bigger this deviation term is and thus the better the
b
power of the test statistics
p is. By using a good estimator P of P , we could expect the deviation term
−1
−1
to be close to hσu σ
n trace(P ).
This leads to the following scheme:
b by solving (10) with the proxy matrix Pe . We know by Lemma
Step 1: Compute an estimator γ
∗
6 that kγ − γ k1 is small.
Step 2: Obtain a better estimator Pb of P by solving (17) .
Step 3: Compute a new estimator of γ ∗ obtained by solving (10) with Pb instead of Pe and
construct a Tn statistics based on this new estimator and on Pb .
As maximum likelihood estimators are often biased for the estimation of the variance components,
we utilized the marginal log likelihood and propose to consider the following restricted maximum
likelihood estimator


1
1
1
b )> P (V − X γ
b ) + log(det(P −1 )) + log(det(X > P X)) .
Pb = arg min
(V − X γ
(17)
P ≥0 2
2
2

Theorems 3 and 4 are expected to hold with this new Tn statistic.
With this approach, not only do we obtain a better test statistic (more efficient) but we are
able to estimate the error variance of the initial model by a new estimator which can be defined as
σ
b2 = nb
σ 2 /trace(Pb ).

5.2

Multivariate Testing

For testing general multivariate hypotheses of the kind H0 : β ∗ = β 0 versus H1 : β ∗ 6= β 0 where
β ∗ ∈ Rd and d → ∞ (p → ∞) we can easily adapt the procedure of Section 2.
With a little abuse of notation let V = Y − Zβ 0 denote the pseudo-response vector (similar to
previous sections) where now a univariate parameter β0 is replaced with a d-dimensional counterpart
b as defined before, we consider d estimators θb(1) , θb(2) , · · · , θb(d) defined as follows
β 0 . Then with γ
b(j)
θ

such that

∈

arg min kθ (j) k1
θ (j) ∈Rp−d
kn−1 X > (Z (j) − Xθ (j) )k∞ ≤ ηθ,j
n−1 Z >
(j) (Z (j) − Xθ (j) ) ≥ η̄θ,j
kZ (j) − Xθ (j) k∞ ≤ µθ,j
kn−1 X > Pe (Z (j) − Xθ (j) )k∞ ≤ ηθ,j .

(18)

Then, we consider to reject the null hypothesis whenever Tn = maxj Tn,j is larger than the
bootstrap quantile q1−α to be defined below. The test statistics Tn,j are defined in the same spirit of

20

the previous section and take the form of
Tn,j = n−1/2

n
X

Tij ,

Tij =

i=1
2

b)
(Zi(j) − X i θb(j) )> Pe (vi − X i γ
σ
bu,j σ
b

2 = n−1 kZ
2
b )) , σ
where σ
bε2 = n−1 Pe (V − X γ
bu,j
(j) − Xθ (j) k2 .
2
The quantile q1−α is defined as a 1−α quantile of the distribution of a bootstrapped test statistics
e
Tn = maxj Ten,j with
n
X
Ten,j = n−1/2
ξi (Tij − n−1/2 Tn,j )
i=1

for a class of multipliers {ξi }ni=1 that are drawn independently from the data and that follow a
standard Gaussian distribution.
For a more general hypothesis H0,j : βj∗ = β0,j for j ∈ J , J ⊂ {1, 2, · · · , p} we can consider
the maximum as the test statistic and denote by Tn (β0,j ) the test statistics as in (13) where one
∗
element
of β ∗ is hold-out
h
i and the remaining ones are stacked in the vector γ . Then, we use GJ (c) =
P maxj∈J (|Ten,j |) ≤ c , Then, the p-value for H0,J : β ∗ = β0 , against the alternative being the
complement, is defined as PJ = 1 − GJ (maxj∈J |Tn (β0,j )|)).

5.3

Generalized Linear Mixed Models

We note that inference in generalized linear models is an extremely difficult problem, even in lowdimensional setting. Typical approaches are hindered by a difficult numerical integration and a
often non-analytic expression of a likelihood or profile likelihood function. In this subsection we
illustrate how the methodology introduced in Section 2 can be easily extended for this difficult
setting. Robustness of our approach provides flexibility regarding specifying the likelihood exactly.
Let Y be the observed data vector and, conditional on the random effects, b, assume that the
elements of Y are independent and drawn from a distribution in the exponential family (which, for
simplicity of exposition, we take with the canonical link). To complete the specification, assume a
distribution for b to be dependent on variance parameters, D:


yη i − c(η i )
∗
∗
(19)
+ d(y, φ)
fyi |b (y|b, γ , β , φ) = exp
a(φ)
where η i = X i γ ∗ + Z i β ∗ + Wi bi . With g denoting the link function, we have b0 = g −1
E[yi |b] = b0 (X i γ ∗ + Z i β ∗ + Wi bi ).
We propose the following test statistic for use in generalized linear mixed models
Tn,j =

b + Zβ 0 ))
n−1/2 (Z (j) − X θb(j) )> Pe (Y − b0 (X γ
σ
bu,j σ
b
2

2 = n−1 b00 (X γ
b + Zβ 0 )) , σ
b + Zβ 0 )kZ (j) − X θb(j) k22 .
bu,j
with σ
bε2 = n−1 Pe (Y − b0 (X γ
2

For the procedure to be adaptive to generalized linear models, the estimators of γ ∗ and θ ∗(j) need
to be carefully developed. Regarding the estimation of γ ∗ we adapt the estimator of Section 2.2 and

21

propose the following estimator
b∈
γ

arg minkγk1
γ ∈Rp−d

such that



n−1 X > Pe Y − b0 (Xγ + Zβ 0 )
≤ ηγ

∞
n−1 Y > Pe Y − b0 (Xγ + Zβ 0 ) ≥ η̄γ
kPe (Y − b0 (Xγ + Zβ 0 ))k∞ ≤ µγ

for suitable choices of tuning parameters ηγ 
estimator for θ ∗(j) can now be defined with
b(j)
θ

such that

∈

(20)

p
p
n−1 log(p), 0 < η̄γ < σε2 and µγ  log(n). The
arg min kθ (j) k1

θ (j) ∈Rp−d

b + Zβ 0 )X > (Z (j) − Xθ (j) )k∞ ≤ ηθ,j
kn−1 b00 (X γ
b + Zβ 0 )X > Pe (Z (j) − Xθ (j) )k∞ ≤ η̄θ,j
kn−1 b00 (X γ
−1
00
b + Zβ 0 )Z >
n b (X γ
(j) (Z (j) − Xθ (j) ) ≥ η̄θ,j
kZ (j) − Xθ (j) k∞ ≤ µθ,j ,

(21)

p
p
2 , and µ
log(n).
for suitable choices of tuning parameters ηθ,j  n−1 log(p), 0 < η̄θ,j < σu,j
θ,j 
>
n
Here, Z (j) = (z1,(j) , · · · , zn,(j) ) ∈ R . Observe that differently from the linear mixed models, in
b(j) are dependent. The
b and θ
the case of the generalized linear mixed models, the two estimators γ
above procedure can be solved using iterations of linear programs, much in the spirit of weighted
least squares methods. The final test statistic is now defined in a similar manner as before with
Tn = maxj Tn,j where the multiplier bootstrap of the previous subsection can be successfully applied.

6

Discussion

This paper proposed a class of test statistics for performing inference on fixed effects that allow for
high-dimensional and misspecified linear mixed models all while maintaining the benefits of classical
methods, i.e., an asymptotically normal and unbiased test statistic with valid and honest confidence
intervals. Our test statistics can be thought of as a doubly robust approach; it adapts to the sparsity
of the nuisance parameters (fixed) as well as the unknown structure of the random effects (both
variance and distributions). Such adaptivity seems essential for modern large-scale applications with
many features, various sparsity assumptions as well as distributional assumptions that cannot be
checked.
In general, the challenge in using adaptive methods as the basis for valid statistical inference
is that selection bias can be difficult to quantify. In this paper, pairing the initial model with a
complementary feature model enabled us to accomplish this goal in a simple yet principled way.
In our simulation experiments, our method provides better error control while achieving nominal
coverage rates in moderate sample sizes.
A number of important extensions and refinements are left open. Our current results only provide
point-wise confidence intervals; extending our theory to the setting of multivariate testing or the setting of a generalized linear mixed models, seems like a promising avenue for further work. Another
challenge is the selection of the proxy matrix M towards better efficiency and power. A systematic approach to design optimization and theory for such setting, would improve the finite sample
performance. In general, work can be done to identify methods that furthermore allow for inference
on the variance components and tests of heterogeneity even in more challenging circumstances, e.g.,
with small samples or a large number of covariates are likely to be bring impactful work to a broader
22

scientific audience.
In this document we provide details to all of the technical results of the main document. In
particular we provide detailed proofs of Theorems 1-4 and establish a sequence of useful Lemmas and
their respective proofs.

A

Proofs of the main results

A.1

Proof of Theorem 1

We observe that Theorem 1 is a direct consequence of the proofs of Theorems 3 and 4.

A.2

Proof of Theorem 2

Theorem 2 is a consequence of Lemmas 5, 6 and 7 presented in Section B.

A.3

Proof of Theorem 3 (under H0 )

Let us first notice that
b > X > Pe (V − X γ
b ) + n−1/2 (θ ∗ − θ)
b) .
σ
bu σ
bTn = n−1/2 U > Pe (V − X γ
|
{z
} |
{z
}
I1

(22)

I2

By Theorem 1.6 of Rudelson and Zhou (2013), the restricted eigenvalue condition holds for
(kθ ∗ k0 , κ, X). Furthermore, by Lemma 5 we have P(A) → 1 where
A = {θ ∗ and γ ∗ are in the feasible regions of (12) and (10) } .
This, together with the restricted eigenvalue condition and Lemma 6 implies that with probability
tending to 1
θb − θ ∗

1

≤ 8ηθ kθ ∗ k0 κ−2 and n−1 X(θb − θ ∗ )

2
2

≤ 16ηθ2 kθ ∗ k0 κ−2 .

(23)

Let us first bound I2 . We can write

|I2 | ≤ n1/2 θb − θ ∗

1

b)
n−1 X > Pe (V − X γ

∞

.

Using (23) and the first constraint of the optimization problem (10), it holds on M that
|I2 | ≤ 8n1/2 ηθ ηγ kθ ∗ k0 κ−2 .
−1 = (n−1 (V − X γ
b )> Pe Pe (V − X γ
b ))−1/2 . By Lemma 3 used with A = Pe ,
Let us now bound |b
σ |√
>
−1
−1
it holds that |b
σ | ≤ η̄γ
n−1 V V which is OP (1) by Lemma 2. We thus have

p
I2
≤ 8n1/2 ηθ ηγ kθ ∗ k0 κ−2 η̄γ−1 n−1 V > V = oP (1)
(24)
σ
b
√
where we used the sparsity condition kθ ∗ k0 = o( n/ log(p)/ log(n)) together with the rate of the
tuning parameters.
Let us now consider Iσb1 . We write
I1
= U >D
σ
b
23

b ) and observe that by construction, D > D = 1. Let us check that
with D = n−1/2 σ
b−1 Pe (V − X γ
max1≤i≤n |di ui | and max1≤i≤n d2i are both oP (1) in order to apply Lemma 9. First let us observe that
!
r
log n
.
max |di | = OP
1≤i≤n
n
Indeed, we have
p
(ii)
≤ n−1/2 η̄γ−1 n−1 V > V µγ = OP

(i)

b)
max |di | = n−1/2 σ
b−1 Pe (V − X γ

1≤i≤n

∞

r

log n
n

!

where (i) follows
from Lemma 3 and from the third constraint of (10) and (ii) follows from Lemma
√
2 and µγ  log n.
√
We also observe that max1≤i≤n |ui | = OP ( log n) by the union bound and the sub-Gaussian
property. We thus have
max |di ui | ≤ max |di | max |ui | = OP (n−1/2 log n) = oP (1)

1≤i≤n

1≤i≤n

1≤i≤n

and max1≤i≤n wi2 = OP (n−1 log n) = oP (1). It follows by Lemma 9 applied on {ui }ni=1 and {di }ni=1
that
I1 d
→ N (0, σu2 )
(25)
σ
b
d

By combining (22), (24) and (25) we have σ
bu Tn → N (0, σu2 ).
b
Since n−1/2 X(θ ∗ − θ)
= oP (1), we have
2

σ
bu = n−1/2 Z − X θb

2

= n−1/2 kU k2 + oP (1) = σu + oP (1).

d

Using Slutzky’s lemma, we obtain that Tn → N (0, 1) which ends the proof of Theorem 3.

A.4

Proof of Theorem 4 (under H1 )

First, note that by Lemma 5 we have P(A) → 1 where
A = {θ∗ and γ ∗ are in the feasible regions of (12) and (10) } .
Furthermore, by Theorem 1.6 of Rudelson and Zhou (2013) and since Pe satisfies the P -condition,
the restricted eigenvalue condition holds for (kθ ∗ k0 , κ, X) and for (kγ ∗ k0 , κ, Pe 1/2 X) with probability
tending to 1. It follows by Lemma 6 applied to problems (12) and (10) that we probability tending
to 1
θb − θ ∗

≤ 8ηθ kθ ∗ k0 κ−2

and

kb
γ − γ ∗ k1 ≤ 8ηγ kγ ∗ k0 κ−2

and

1

n−1 X(θb − θ ∗ )

2

1/2
n−1 Pe X(b
γ − γ ∗)

24

≤ 16ηθ2 kθ ∗ k0 κ−2

2
2
2

≤ 16ηγ2 kγ ∗ k0 κ−2 .

(26)

Since V = Xγ ∗ + Zhn−1/2 + W b + , we can write the following decomposition:
b
σ
bu σ
bTn = n−1/2 (W b + )> Pe (Z − X θ)
{z
}
|
T1

+n
|

−1/2

b + n−1 hZ > Pe (Z − X θ)
b .
b )> X > Pe (Z − X θ)
(γ − γ
{z
} |
{z
}
∗

T2

(27)

T3

Let us first consider T2 . We observe that

b
|T2 | ≤ n1/2 kb
γ − γ ∗ k1 n−1 X > Pe (Z − X θ)

∞

≤ n1/2 8ηγ kγ ∗ k0 κ−2 ηθ0 = oP (1)

(28)

where we used (26) and the fourth constraint of (12).
Regarding T3 , we write

b := T31 + T32 .
T3 = n−1 hZ > Pe U + n−1 hZ Pe X(θ ∗ − θ)

For T32 , we observe that Pe Z is sub-Gaussian since Z is sub-Gaussian and Pe satisfies the P condition. This together with (26) implies on M that
q
−1/2 e
−1/2
∗
b
|T32 | ≤ h n
PZ
n
X(θ − θ) = OP (ηθ kθ ∗ k0 ) = oP (1).
2

2

For T31 , we use the Markov law of large numbers for independent but non-identical random
variables. Working on a group level, we write
T31 = hn−1

N
X
j=1

Z j Pe j U j = hn−1

N
X

dj

j=1

with Pe = diag(Pe 1 , . . . , Pe N ) and Pe j ∈ Rnj ×nj and we observe that dj ⊥ dj 0 for g 6= j 0 . Let us show
that the mean and the variance of dj are bounded. We compute
∗> > e
2
e
e
E[dj ] = E[U >
j P j U j ] + E[θ X j P j U j ] = σu trace(P j )

where we make use of the independence of X and U . For the variance term, we make use of the
variance formula of a quadratic form (see appendix) and obtain
4
e
e e
Var[U >
j P j U j ] = 2σu trace(P j P j ).

This implies that the variance of dj is bounded. By using the Markov’s law of large numbers on
{dj }N
j=1 we obtain that
N
X
−1
N
(dj − σu2 trace(Pe j )) = oP (1).
j=1

−1 2
e
Multiplying this result by the bounded quantity h N
n , we get T31 − hn σu trace(P ) = oP (1). We thus
have
T3 = hσu2 n−1 trace(Pe ) + oP (1).
(29)

Let us now prove that T1 converges to a normal distribution under the alternative hypothesis.
We write
T1 /b
σu = r > Pe (W b + )
25

b Let us define ∗ = Pe (W b + ) and the sequences of random vectors
with r = n−1/2 σ
b−1 (Z − X θ).
n uoN
{r j }N
∗j
induced by the group structure of the mixed model. Let us show that the
j=1 and
j=1

conditions of Lemma 10 are satisfied for these sequences of random vectors.
First, we see that ∗j ∈ Rnj are mean zero, sub-Gaussian random vectors (as linear combination of
sub-Gaussian random variables) and such that var(∗j ) = K j with K j being a block diagonal element
of σ2 Pe P −1 Pe . By construction, the elements of K j are bounded which implies that λmax (K j ) is
P
>
>
bounded in N . Second, we notice that by definition of σ
bu , we have that N
j=1 r j r j = r r = 1.
As proven in Theorem 3, it holds that σ
bu = σu + oP (1). By the continuous mapping theorem, we
deduce that σ
bu−1 = σu + oP (1) = OP (1). Thus, by the third constraint of (12) we have that
krk∞ ≤ n−1/2 µθ /b
σu = OP (n−1/2 µθ ) = oP (1).

2
This implies that maxj r >
j r j ≤ maxj nj krk∞ = oP (1) since the nj are bounded. Furthermore,
√
we observe that since ∗ is sub-Gaussian, it holds by the union bound that k∗ k∞ = OP ( log n). We
thus have
p
∗
∗
−1/2
max r >
µθ log n) = oP (1).
j j ≤ k rk∞ k k∞ OP (n
1≤j≤n

n oN
We can now apply Lemma 10 on {r j }N
and
∗j
j=1

j=1

and derive that

σ
b−1 T
d
qP u 1
→ N (0, 1).
N
>
j=1 r j K j r j

By Lemmas 7 and 8, it holds that σ
b=
Slutsky’s Theorem implies that

qP
N

>
j=1 r j K j r j

(30)

+ oP (1). This, together with (30) and

d

σ
bu−1 σ
b−1 T1 → N (0, 1).

Now, using σ
bu = σu + oP (1) we can combine (28), (29) and (31) and obtain that

n−1 trace(Pe )
d
Tn = (b
σ −1 σ
bu−1 )(T1 + T2 + T3 ) = N (0, 1) + hσu σ−1 q
+ oP (1).
trace(n−1 Pe P −1 Pe )

(31)

(32)

This ends the proof of Theorem 4.

B

Auxiliary proofs and statements

Lemma 1. The matrix P = (I n − W E −1 W > )(I n − W E −1 W > ) + σ2 W E −1 Ψ−1 E −1 W > , with
E = W > W + σ2 Ψ−1 , can be rewritten as P = (I n + σ−2 W ΨW > )−1 .
Proof of Lemma 1. A proof can be found in the supplementary material of Fan and Li (2012).
Lemma 2. Under H0 : β ∗ = β0 and under conditions of Theorem 3, it holds that n−1 V > V = OP (1).
Proof of Lemma 2. Under H0 we have V = Xγ ∗ + ∗ with ∗ =  + W b. Recall that var(∗ ) =

N
P −1 . Let us work at a group level and consider the sequence of random variables V >
j V j j=1 . We
26

have
∗> ∗
∗> >
∗
E[V >
X j X j γ ∗ ] + E[∗>
j V j ] = E[j j ] + E[γ
j Xjγ ]
∗>
= trace(P −1
Σγ ∗ + 0
j ) + nj γ

= O(1)
by condition 1 and by the boundedness of the elements of P −1 . With similar arguments, we can
show that Var[V >
j V j ] is bounded. It follows from the Markov’s law of large number that
N −1

N
X
p
>
(V >
j V j − E[V j V j ]) → 0.
j=1

which implies that n−1 V > V = OP (1).
b ) ≥ η̄, then n−1 (V −
Lemma 3. Let A be a n × n symmetric matrix. If n−1 V > A(V − X γ
>
>
2
−1
b ) AA(V − X γ
b ) ≥ η̄ /(n V V ).
Xγ
Proof of Lemma 3.
For any a > 0, we have

b )> AA(V − X γ
b)
n−1 (V − X γ

b )> AA(V − X γ
b ) + a(η̄ − n−1 V > A(V − X γ
b ))
≥ n−1 (V − X γ
n
o
≥ min n−1 (V − Xγ)> AA(V − Xγ) + a(η̄ − n−1 V > A(V − Xγ))
γ

1
= aη̄ − a2 n−1 V > V
4

where the last equality follows from the first-order condition of the quadratic problem which is


1
−1 >
aV − A(V − Xγ) = 0.
n X A
2
Then by choosing a = 2η̄/(n−1 V > V ) we obtain
b )> AA(V − X γ
b ) ≥ η̄ 2 /(n−1 V > V ).
n−1 (V − X γ

Lemma 4. The matrices P and P −1 satisfy the P -condition where P = (I n + σ−2 W ΨW > )−1 .
Furthermore, for any N q × N q matrix M having the same block diagonal structure as W > W , it
holds that the matrix Pe = (I n + W M W > )−1 , if it exists, satisfies the P -condition.

Proof of Lemma 4. First, note that by construction P , P −1 and Pe have the correct block
diagonal structure. Furthermore, by condition 1 it is clear that the elements of P −1 are bounded.
Let us now show that each element of P and Pe is bounded. Recall that the `2 norm of a matrix
k A xk
A is defined as kAk2 = supx6=0 kxk 2 = σmax (A) where σmax (A) is the largest singular value of A.
2
By consequence, it suffices to show that σmax (Pe ) and σmax (P ) are bounded to guarantee that each

27

element is bounded.
n
o
σmax (Pe ) = max σ1 (Pe ), . . . , σn (Pe )
n
o
−1
−1
= max σ1 (Pe )−1 , . . . , σn (Pe )−1
n
o
= max (1 + σ1 (W M W > ))−1 , . . . , (1 + σn (W M W > ))−1
≤1

since by definition a singular value is always non-negative. The proof for P and P −1 follow from
similar arguments.
The following Lemma, combined with Lemma 4 guarantees the feasibility of the optimization
problems (12) and (10), which is required to derive the properties of our estimators.
Lemma 5. Let A ∈ Rn×n satisfying the P -condition. Suppose that condition 1 holds and that
β ∗ = β0 + hn−1/2 for a constant
h ∈ R. Then there
p
√ exist constants−1C1 , C2 , 2 c0 >−10 such that
0
−1
for any ηγ , ηθ , ηθ ≥ C1 log n n log p, µγ , µθ ≥ C2 log n, η̄γ , ∈ 0, n trace(σ AP ) − c0 and
η̄θ , ∈ 0, σu2 − c0 , it holds with probability tending to 1 that
n−1 X > A(V − Xγ ∗ ) ∞ ≤ ηγ
n−1 V > A(V − Xγ ∗ ) ≥ η̄γ
kA(V − Xγ ∗ )k∞ ≤ µγ

and

n−1 X > (Z − Xθ) ∞
n−1 Z > (Z − Xθ)
kZ − Xθk∞
>
−1
n X A(Z − Xθ) ∞

≤
≥
≤
≤

(33)

ηθ
η̄θ
µθ
ηθ0

(34)

Proof of Lemma 5. We prove that (33) holds with probability tending to 1. The proof for (34)
follows from similar arguments and is thus omitted.
Let us first show that there exists a constant
C1 such that with probability converging to 1 it
p
holds that n−1 X > A(V − Xγ ∗ ) ∞ ≤ C1 n−1 log p. Using the block diagonal structure of A =
(A1 , . . . , AN ) associated with the group structure, we observe that
n−1 X > A(V − Xγ ∗ ) =

N

N

g=1

g=1

1X >
1X
X (g) Ag (V (g) − X (g) γ ∗ ) =
dg
n
n

(35)

where the subscript (g) refers to the fact that we extract the rows corresponding to goup g.
Claim 1: Each component of dg ∈ Rp has an exponential-type tail with parameters (K1 , 12 ).
To show this, let us decompose the j-th component as follows:
X
X
∗
djg =
xij aik (vk − X >
(36)
k γ ).
i∈groupg k∈groupg

As the group g contains a bounded number of elements ng , it suffices to show that each xij aik (vk −
1
∗
X>
k γ ) has an exponential-type tail with parameters (K0 , 2 ) to prove the claim.
We note that
∗
−1/2
vk − X >
+ k + wk bg
k γ = zk hn
28

with var(bg ) = ψ ∈ Rq×q . The sub-Gaussianity of k and of each component of bg and the fact that
each component of wk is bounded imply by the properties of sub-Gaussian random variables that
∗
vk − X >
k γ is sub-Gaussian. Now, using that xij is sub-Gaussian and that aij is O(1) we obtain by
∗
Lemma 5.14 of Vershynin (2010) that xij aik (vk −X >
k γ ) has an exponential-type tail with parameters
1
(K0 , 2 ). This concludes the proof of claim 1.
n oN
, it follows from
Observe that djg ⊥ djg0 . Thanks to the exponential-type tail property of djg
g=1

Proposition 5.16 of Vershynin (2010) that there exist constants K2 and K3 such that for any x > 0
we have


N
X
(djg − Edjg ) > x ≤ K2 exp(−K3 x2 ).
(37)
P N −1/2
g=1

We compute

Edjg

i
h
−1/2
>
hZ (g) + (g) + W (g) bg )
= E X (g)j Ag (n
h
i
(i)
−1/2
∗
= E X>
A
n
hX
θ
g
(g)
(g)j
(ii)

= hn−1/2 trace(Ag )Σj: θ ∗

where (i) comes from the independence of X (g)j with u,  and W b and (ii) comes from the independence between xkj and xlj for k 6= l, this is,
h
i
X
X
X
E X>
E [xkj Akl xlj ] =
E [xkj Akk xkj ] = trace(Ag )Σj: .
(g)j Ag X (g) ) =
k∈groupg l∈groupg

k∈groupg

P
j
−1/2 N −1 trace(A)Σ . Now, observe that since A satisfies
Thus we obtain that N1 N
j:
g=1 E[dg ] = hn
the P -condition it holds that n−1 trace(A) is bounded. It follows that there exists a constant K4 such
that
N
1 X
n p −1
E[djg ] ≤ K4 trace(n−1 A)
n log p
(38)
N
N
g=1

Combining (35), (37) and (38) and using the union bound property and the inequality (a − b)2 ≥

29

a2
2

− b2 , we derive the following result:


> ηγ
P n−1 X > A(V − Xγ ∗ )
∞


−1 >
∗
≤ p P n X j A(V − Xγ ) > ηγ


N
X
n
= p P  N −1/2
djg > ηγ √ 
N
g=1


N
X
p
n
n
= p P  N −1/2
(djg − E[djg ]) > ηγ √ − K4 trace(n−1 A) √
n−1 log p
N
N
g=1
"

2 #
p
n
n
n−1 log p
≤ K2 p exp −K3 ηγ √ − K4 trace(n−1 A) √
N
N



1 2n
2
−1
2 n −1
≤ K2 p exp −K3 n
η
− K4 trace(n A)
n log p
2 γN
N


n
K3 2 n
1+K3 K42 trace(n−1 A)2 N
exp −
nηγ
.
= K2 p
2
N

By definition ofpA, trace(n−1 A) ≤ log2 n. This implies that for a constant
C1 large enough and

for ηγ ≥ C1 (log n) n−1 log p we have P n−1 X > A(V − Xγ ∗ ) ∞ ≤ ηγ with probability tending
to 1.
This ends the first part of the proof.
Let us now tackle √
the third constraint of (33) and show that
√ there exists a constant C2 > 0 such
∗
that for any µγ > C2 log n, we have kA(V − Xγ )k∞ ≤ C2 log n with probability tending to 1.
We first observe that each component of A(V − Xγ ∗ ) is sub-Gaussian because of the structure of A
and because components of V − Xγ ∗ are sub-Gaussian. By the union bound and the definition of a
sub-Gaussian variable, there exists a constant K4 > 0 such that for any x > 0


p
P kA(V − Xγ ∗ )k∞ > x log n ≤ n exp(1 − K5 x2 log n).
(39)

The result follows by taking a constant x large enough. To conclude, let
 us tackle the second
−1
−1
2
constraint of (33). Let us show that for any η̄γ ∈ 0, trace(n σ AP ) − c0 , we have


P n−1 V > A(V − Xγ ∗ ) < η̄γ → 0.
(40)
Working at the group level induced by the block-diagonal structure of A, we write
n−1 V > A(V − Xγ ∗ ) = n−1

N
X
j=1

V j Aj (V j − X j γ ∗ ) := n−1

Let us consider the following decomposition:

N
X

bj

j=1

∗
bj = (V j − X j γ ∗ )> Aj (V j − X j γ ∗ ) + γ ∗> X >
j Aj (V j − X j γ ) .
{z
} |
|
{z
}
bj,1

bj,2

Using the equality

V j − X j γ ∗ = W j bj + j + hn−1/2 (X j θ ∗ + uj )
30

and the independence between X j and j , bj and uj , we obtain that
∗
E[bj,2 ] = hn−1/2 γ ∗> E[X >
j Aj X j ]θ

(41)

≤ hn−1/2 trace(Aj ) kγ ∗ k2 kΣ∗ k2 kθ ∗ k2 = o(1).

(42)

Furthermore, we can show that the variance of bj,2 is bounded in N . Note that bj,1 is a quadratic
form. Using that E[V j − X j γ ∗ ] = 0 and Cov[V j − X j γ ∗ ] = σ2 P −1
j + o(1), it follows that E[bj,1 ] =
−1
2
trace(σ P j Aj ) + o(1) and
2 −1
Var[bj,2 ] = trace(σ2 P −1
j Aj σ P j Aj )

which is bounded in N . By consequence the mean of the random variable bj is given by E[bj ] =
trace(σ2 P −1
j Aj ) + o(1) and its variance is bounded in N . It follows from Markov’s law of large
numbers that
N
 p
1 X
bj − trace(σ2 P −1
A
)
→ 0.
j
j
N
j=1

This implies that n−1 V > A(V − Xγ ∗ ) = n−1 trace(σ2 P −1 A) + oP (1) which proves (40). This
ends the proof of Lemma 5.
b
b and θ.
The following lemma is used to derive the properties of the Dantzig-like estimators γ

Lemma 6. Let Y ∈ Rn , X ∈ Rn×p and C ⊆ Rp . Let A ∈ Rn×n satisfying the P -condition. Suppose
there exists β ∗ ∈ C such that n−1 X > A(Y − Xβ ∗ ) ∞ ≤ η. Let βb = argminβ ∈C kβk1 subject
to n−1 X > A(Y − Xβ) ∞ ≤ η. If s = kβ ∗ k0 and if there exists κ > 0 such that the restricted
eigenvalue condition holds for (s, κ, A1/2 X), this is

min

J0 ⊆{1,...,p}
|J0 |≤s

min

δ6=0
≤kδ J0 k1

δ J0c

A1/2 Xδ
2
√
≥κ
n kδ J0 k2

1

then for δ = βb − β ∗ it holds that kδk1 ≤ 8ηsκ−2 and n−1 δ > XAXδ ≤ 16η 2 sκ−2

Proof of Lemma 6. First note that β ∗ satisfies the constraints of the optimization problem
whose βb is the optimal solution. This implies that
βb

1

≤ kβ ∗ k1 .

(43)

Secondly observe that for the true active set S0 of β ∗ it holds that
β ∗ + δ S0c

1

= kβ ∗ k1 + δ S0c

1

(44)

.

Thirdly, combining (43), (44) and the triangle inequality |a| − |b| ≤ |a + b| with a = β ∗ + δ S0c and
b = δ S0 we obtain
(45)
δ S0c 1 ≤ kδ S0 k1 .

Indeed, we have

kβ ∗ k1 + δ S0c

1

− kδ S0 k1 = β ∗ + δ S0c

1

− kδ S0 k1 ≤ β ∗ + δ S0c + δ S0

31

1

= βb

1

≤ kβ ∗ k1 .

Now, using the restricted eigenvalue condition together with (45), we have
A1/2 Xδ

≥

2

√
nκ kδ S0 k2 .

(46)

We can then derive the following result:
n−1 A1/2 Xδ

2
2

= n−1 δ > X > AXδ
≤ kδk1 n−1 X > AX(βb − β ∗ )

(i)

∞

≤ 2η kδk1

√
≤ 4η s kδ S0 k2

(ii)

√
≤ 4η sκ−1 n−1/2 A1/2 Xδ

(iii)

2

√
where (i) follows from the definitioin of η, (ii) follows from (45) and kδ S0 k1 ≤ s kδ S0 k2 , and (iii)
√
follows from (46). This implies that n−1/2 A1/2 Xδ ≤ 4η sκ−1 and kδk1 ≤ 8ηsκ−2 , which ends
2
the proof of Lemma 6.
b)
Lemma 7. Under conditions of Theorem 4 and for σ
b2 = n−1 Pe (V − X γ

it holds that σ
b2 = n−1 trace(K) + oP (1).

2
2

and K = σ2 Pe P −1 Pe ,

Proof of Lemma 7. Using V = Xγ ∗ + hn−1/2 Z + W b +  and writing ∗ = Pe (W b + ), we
first observe that
b)
σ
b2 = n−1 Pe (V − X γ

2
2

b ) + hn−1/2 Pe Z
= n−1 Pe (W b + ) + Pe X(γ ∗ − γ
=n

−1

Pe (W b + )

2
2

2
2

2

+ oP (1)

= n−1 ∗> ∗ + oP (1)
b)
since n−1 Pe X(γ ∗ − γ

2

(47)

1/2
b)
≤ λmax (Pe ) n−1 Pe X(γ ∗ − γ

2
2

= oP (1) by (26) and by the fact that

λmax (Pe ) is bounded in n and since n−1/2 hn−1/2 Pe Z ≤ hn−1 Pe kZk2 = oP (1) by the sub2
2
Gaussianity of zi .
Let us define K = σ2 Pe P −1 Pe and observe that Var(∗ ) = K. Working at the group level induced
n oN
and corresponding covariance matrices {K j }N
by the mixed model, we consider ∗j
j=1 . We
j=1
oN
n
∗
∗
observe that ∗>
is a sequence of independent random variables with mean µj = E[∗>
j j ] =
j j
j=1

∗
trace(K j ) and Var[∗>
j j ] = 2trace(K j K j ) < ∞. The Markov’s law of large numbers implies that
N
N
1 X ∗> ∗
1 X
p
j j −
trace(K j ) → 0.
N
N
j=1

j=1

This together with (47) concludes the proof of Lemma 7.
32

b
Lemma 8. Under conditions of Theorem 4 and for K = σ2 Pe P −1 Pe and w = n−1/2 σ
bu−1 (Z − X θ),
>
−1
it holds that w Kw = n trace(K) + oP (1).

Proof of Lemma 8. Let us define w∗ = n−1/2 σu−1 (Z − Xθ ∗ ) and observe that E[w∗ ] = 0 and
Var(w∗ ) = n−1 I n and
o
σu n ∗
w − n−1/2 σu−1 X(θb − θ ∗ ) .
w=
σ
bu
We consider the following decomposition:

 2 
σu  ∗>
w> Kw =
w
Kw}∗ − 2n−1/2 σu−1 w∗> KX(θb − θ ∗ )
| {z
|
{z
}
σ
bu 

T1
T2




n
o>
−1 −2
∗
∗
+ n σu X(θb − θ ) KX(θb − θ )

|
{z
}

T3

First, notice that T1 is a quadratic form whose mean and variance are given by:
E[w∗> Kw∗ ] = trace(n−1 I n K) = n−1 trace(K)

and

Var[w∗> Kw∗ ] = 2trace(n−1 I n Kn−1 I n K) = 2n−2 trace(KK) = O(n−1 )

where the block diagonal structure of K is used to guarantee that the elements of KK are bounded
in n. Thus the random variable w∗> Kw∗ − n−1 trace(K) has mean 0 and variance converging to 0.
It follows from Chebyshev’s inequality that
T1 = n−1 trace(K) + oP (1).

(48)

Secondly, observe that T2 can be treated as T32 in the proof of Theorem 2. Indeed using that
w∗ = n−1/2 σu−1 u, we have
|T2 | = 2σu−2 n−1 u> KX(θb − θ ∗ )

≤ 2σu−2 n−1 kKuk2 X(θb − θ ∗ )

2

= OP (ηθ

where we use (26) and the sub-Gaussianity of Ku.
Thirdly, observe that
|T3 | = σu−2 n−1 K 1/2 X(θb − θ ∗ )

2
2

≤ σu−2 K 1/2

2
2

q
kθ ∗ k0 ) = oP (1)

n−1 X(θb − θ ∗ )

where the last equality follows from (26) and the fact that K 1/2

2
2

w Kw =



σu
σ
bu

2

2

= oP (1).

(50)

≤ λmax (K) is bounded in n.

Since σ
bu = σu +oP (1), the continuous mapping theorem implies that
with (48), (49), (50) and Slutzky’s theorem we obtain that
>

2

(49)

σu
σ
bu

p

→ 1. Hence, this together

(T1 + T2 + T3 ) = n−1 trace(K) + oP (1).

The following Lemma was presented in a preliminary version of Zhu and Bradic (2016) and is
33

used in the proof of Theorem 4. It is presented here for the sake of completeness.
Lemma 9. Let {wj }nj=1 and {hj }nj=1 be two independent sequences of random variables. Suppose
that
(i) hj and hj 0 are independent for j 6= j 0 , E[hj ] = 0, σh2 = E[h2j ] = O(1) and κ4 = E[(h2j − σh2 )2 ] =
O(1)
Pn
2
2
(ii)
j=1 wj = 1, max |wj hj | = oP (1) and max wj = oP (1).
j=1,...,n

Then,

Pn

j=1 wj hj

j=1,...,n

d

→ N (0, σh2 ).

The following lemma is an extension of Lemma 9 to the multivariate case. It is useful to obtain
b and
the normality of σ
bu−1 T1 in the proof of Theorem 4 where it is used with w = n−1/2 σ
bu−1 (Z − X θ)
e
h = P (W b + ).
N
Lemma 10. Let {wj }N
j=1 and {hj }j=1 be two independent sequences of random vectors with w j , hj ∈
Rnj and nj bounded in N . Suppose that

(i) hj and hj 0 are independent for j 6= j 0 , E[hj ] = 0, E[hj h>
j ] = K j with
λmax (K j ) = O(1)
and

lm
E[(hij hkj − Kjik )(hlj hm
j − Kj )] = O(1)

for 1 ≤ i, k, l, m ≤ nj .
PN
>
>
>
(ii)
j=1 w j w j = 1, max |w j hj | = oP (1) and max w j w j = oP (1).
j=1,...,N

Then,

j=1,...,N

PN

w>
d
j hj
→ N (0, 1).
>
j=1 w j K j w j

qP j=1
N

Proof of Lemma 10. Let FN,0 be the σ-algebra generated by W = {wj }N
and FN,j be
Pj j=1 >
j
2 =
>
the σ-algebra generated by W and {hk }k=1 . Define XN j = wj hj and SN j = k=1 wk hk , UN
j
P
Pj
N
>K w .
> h )2 and u2 =
w
(w
j j
j=1 j
N
k k
k=1
d

By Theorem 3.4 of Hall and Heyde (1980) it holdsh that SN jNi/UN jN → N (0, 1) if the following
2
2
2
conditions are satisfied: (a) maxj |XN j | = oP (1), (b) E maxj XN
j is bounded in N , (c) UN jN −uN =
P
PN
2
oP (1), (d) N
j=1 E [XN j |FN,j−1 ] = oP (1) and (e)
j=1 |E [XN j |FN,j−1 ]| = oP (1).
Observe that (a) holds from the assumptions of the Lemma. Secondly, notice that

 X
N
h
i
>
2
2
E max(wj hj ) ≤
E (w>
h
)
j
j
j

(51)

j=1

=

N
X
j=1

N
h h
ii X
h
i
>
2
E E (wj hj ) |FN,j−1 =
E w>
j K j wj .
j=1

34

(52)

1/2

Since w>
j K j wj = K j wj

2
2

, we derive that


N
N
h
i X
X
1/2
E w>
K
w
≤
E
Kj
j j
j
j=1

2
2

j=1

kwj k22





N
X

≤ max λmax (K j ) E 
w>
j w j < ∞,
j=1,...,N

(53)

(54)

j=1

P
>
where we use N
j=1 w j w j = 1 and the fact that maxj λmax (K j ) is bounded in N . Thus, (b) is
P
>
>
proven. To show (c), we define QN = N
j=1 w j (hj hj − K j )w j and need to show that QN = oP (1).
First, we derive that
E[Q2N |W ]

=

N X
N
X
j=1 k=1

(i)

=

N
X
j=1

(ii)

=

>
>
>
E[w>
j (hj hj − K j )w j w k (hk hk − K k )w k |W ]

>
>
>
E[w>
j (hj hj − K j )w j w j (hj hj − K j )w j |W ]

N
X
j=1



>
trace wj w>
j Cov[(hj hj − K j )w j |W ]

where (i) follows from the independence assumption of hj and (ii) follows from the mean formula of
2
a quadratic form and from E[(hj h>
j − K j )w j |W ] = 0. Let us now define the nj × nj matrix
wdiag
j



wj
0n
 j
= .
 ..

0nj

and notice that


. . . 0nj
. . . 0nj 

..
.. 
.
. 
. . . wj

0nj
wj
..
.
0nj

>

diag
Cov[(hj h>
K4,j wdiag
j − K j )w j |W ] = w j
j

where K4,j is a n2j × n2j matrix containing all the fourth moment elements of hj . By assumption (i),
all the elements of K4,j are bounded. It follows that
E[Q2N |W ] =
=

N
X
j=1

N
X
j=1

≤
≤

N
X
j=1

N
X

>

diag
w>
K4,j wdiag
wj
j wj
j

K4,j wdiag
wj
j
1/2

1/2 2

K4,j

2

wdiag
j

2
C n2j (w>
j wj )

j=1

35

2
2
2
2

kwj k22

≤C

n2j

max

j=1,...,N

(w>
j wj )

N
X

w>
j wj

j=1

= oP (1),

for a constant C, where we use
cycling property
of the trace, the boundedness of the fourth
P the
P
2
moments of hj , the property j cj ≤ max |cj |
j |cj | and the assumption (ii). By consequent, we
have
E[Q2N ] = E[E[Q2N |W ]] = oP (1)
which implies by Markov’s inequality that Q2N = oP (1). This proves (c).
Regarding (d) and (e), we observe that
h
i
>
E [XN j |FN,j−1 ] = E w>
j hj |FN,j−1 = w j E [hj |FN,j−1 ] = 0

by the independence assumption between hj and h1 , . . . , hj−1 and the independence between hj and
wj .
d

This guarantees that (d) and (e) hold. We thus obtain that SN jN /UN jN → N (0, 1). This together
with (c) and Slutsky’s Theorem implies that
PN

>
j=1 w j hj

qP
N

>
j=1 w j K j w j

d

→ N (0, 1).

References
Anderson, T. W., editor (1984). An Introduction to Multivariate Statistical Analysis. Wiley.
Athey, S., Imbens, G. W., and Wager, S. (2016). Approximate residual balancing: De-biased inference
of average treatment effects in high dimensions. arXiv preprint arXiv:1604.07125.
Belloni, A., Chernozhukov, V., Chetverikov, D., and Wei, Y. (2015). Uniformly valid postregularization confidence regions for many functional parameters in z-estimation framework. arXiv
preprint arXiv:1512.07619.
Belloni, A., Chernozhukov, V., Fernández-Val, I., and Hansen, C. (2017). Program evaluation and
causal inference with high-dimensional data. Econometrica, 85(1):233–298.
Belloni, A., Chernozhukov, V., and Kato, K. (2014). Uniform post-selection inference for least
absolute deviation regression and other z-estimation problems. Biometrika, 102(1):77–94.
Bickel, P. J., Ritov, Y., and Tsybakov, A. B. (2009). Simultaneous analysis of lasso and dantzig
selector. The Annals of Statistics, pages 1705–1732.
Bonnet, A., Gassiat, E., Lévy-Leduc, C., et al. (2015). Heritability estimation in high dimensional
sparse linear mixed models. Electronic Journal of Statistics, 9(2):2099–2129.
Breslow, N. E. and Clayton, D. G. (1993). Approximate inference in generalized linear mixed models.
Journal of the American Statistical Association, 88(421):9–25.
Bühlmann, P. and Van De Geer, S. (2011). Statistics for high-dimensional data: methods, theory and
applications. Springer Science & Business Media.
Cai, T. T., Guo, Z., et al. (2017). Confidence intervals for high-dimensional linear regression: Minimax
rates and adaptivity. The Annals of Statistics, 45(2):615–646.
Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., et al. (2016). Double machine
learning for treatment and causal parameters. arXiv preprint arXiv:1608.00060.
36

Chernozhukov, V., Hansen, C., and Spindler, M. (2015). Valid post-selection and post-regularization
inference: An elementary, general approach. Annual Review of Economics, 7:649–688.
Crainiceanu, C. M. and Ruppert, D. (2004). Likelihood ratio tests in linear mixed models with one
variance component. Journal of the Royal Statistical Society: Series B (Statistical Methodology),
66(1):165–185.
Fan, J. and Li, R. (2001). Variable selection via nonconcave penalized likelihood and its oracle
properties. Journal of the American statistical Association, 96(456):1348–1360.
Fan, Y. and Li, R. (2012). Variable selection in linear mixed effects models. Ann. Statist., 40(4):2043–
2068.
Ghosh, A. and Thoresen, M. (2016). Non-concave penalization in linear mixed-effects models and
regularized selection of fixed effects. arXiv preprint arXiv:1607.02883.
Goeman, J. J., Van Houwelingen, H. C., and Finos, L. (2011). Testing against a high-dimensional
alternative in the generalized linear model: asymptotic type i error control. Biometrika, pages
381–390.
Groll, A. and Tutz, G. (2014). Variable selection for generalized linear mixed models by l1-penalized
estimation. Statistics and Computing, pages 1–18.
Hall, P. and Heyde, C. (1980). Martingale Limit Theory and its Application. Probability and Mathematical Statistics: A Series of Monographs and Textbooks. Academic Press.
Heagerty, P. J. and Kurland, B. F. (2001). Misspecified maximum likelihood estimates and generalised
linear mixed models. Biometrika, 88(4):973.
Hobert, J. P. and Casella, G. (1996). The effect of improper priors on gibbs sampling in hierarchical
linear mixed models. Journal of the American Statistical Association, 91(436):1461–1473.
Hui, F. K., Müller, S., and Welsh, A. (2016). Joint selection in mixed models using regularized pql.
Journal of the American Statistical Association, (just-accepted).
Jankova, J., Van De Geer, S., et al. (2015). Confidence intervals for high-dimensional inverse covariance estimation. Electronic Journal of Statistics, 9(1):1205–1229.
Javanmard, A. and Montanari, A. (2014a). Confidence intervals and hypothesis testing for highdimensional regression. Journal of Machine Learning Research, 15(1):2869–2909.
Javanmard, A. and Montanari, A. (2014b). Confidence intervals and hypothesis testing for highdimensional regression. Journal of Machine Learning Research, 15:2869–2909.
Kenward, M. G. and Roger, J. H. (1997). Small sample inference for fixed effects from restricted
maximum likelihood. Biometrics, 53(3):983–997.
Lachos, V. H., Ghosh, P., and Arellano-Valle, R. B. (2010). Likelihood based inference for skewnormal independent linear mixed models. Statistica Sinica, 20(1):303–322.
Lindstrom, M. J. and Bates, D. M. (1988). Newton?raphson and em algorithms for linear mixed-effects
models for repeated-measures data. Journal of the American Statistical Association, 83(404):1014–
1022.
Litière, S., Alonso, A., and Molenberghs, G. (2007). Type i and type ii error under random-effects
misspecification in generalized linear mixed models. Biometrics, 63(4):1038–1044.
Mack, M., van Loon, A. P., and Hohmann, H.-P. (1998). Regulation of riboflavin biosynthesis
inbacillus subtilis is affected by the activity of the flavokinase/flavin adenine dinucleotide synthetase
encoded byribc. Journal of bacteriology, 180(4):950–955.
McCulloch, C. E. (1997). Maximum likelihood algorithms for generalized linear mixed models. Journal
of the American statistical Association, 92(437):162–170.
McGilchrist, C. (1994). Estimation in generalized mixed models. Journal of the Royal Statistical
Society. Series B (Methodological), pages 61–69.
Mörtl, S., Fischer, M., Richter, G., Tack, J., Weinkauf, S., and Bacher, A. (1996). Biosynthesis of
riboflavin lumazine synthase of escherichia coli. Journal of Biological Chemistry, 271(52):33201–
33207.
37

Müller, S., Scealy, J., Welsh, A., et al. (2013). Model selection in linear mixed models. Statistical
Science, 28(2):135–167.
Neyman, J. (1959). Optimal asymptotic tests of composite statistical hypotheses. Probability and
statistics, 57:213.
Ning, Y., Liu, H., et al. (2017). A general theory of hypothesis tests and confidence regions for sparse
high dimensional models. The Annals of Statistics, 45(1):158–195.
Quintana, F. A., Johnson, W. O., Waetjen, L. E., and Gold, E. B. (2016). Bayesian nonparametric
longitudinal data analysis. Journal of the American Statistical Association, 111(515):1168–1181.
Ren, Z., Sun, T., Zhang, C.-H., Zhou, H. H., et al. (2015). Asymptotic normality and optimalities in
estimation of large Gaussian graphical models. The Annals of Statistics, 43(3):991–1026.
Rohart, F., San-Cristobal, M., and Laurent, B. (2014). Fixed effects selection in high dimensional
linear mixed models. Computational Statistics and Data Analysis.
Rudelson, M. and Zhou, S. (2013). Reconstruction from anisotropic random measurements. IEEE
Trans. Information Theory, 59(6):3434–3447.
Ryzhov, I. O., Han, B., and BradicĞ, J. (2016). Cultivating disaster donors using data analytics.
Management Science, 62(3):849–866.
Schelldorfer, J., Buhlmann, P., and van de Geer, S. (2011). Estimation for high-dimensional linear
mixed-effects models using `1 -penalization. Scandinavian Journal of Statistics, 38(2):197–214.
Song, P. X.-K., Zhang, P., and Qu, A. (2007). Maximum likelihood inference in robust linear mixedeffects models using multivariate t distributions. Statistica Sinica, pages 929–943.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical
Society. Series B (Methodological), pages 267–288.
van de Geer, S., Bühlmann, P., Ritov, Y., and Dezeure, R. (2014). On asymptotically optimal
confidence regions and tests for high-dimensional models. The Annals of Statistics, 42(3):1166–
1202.
Van de Geer, S., Bühlmann, P., Ritov, Y., Dezeure, R., et al. (2014). On asymptotically optimal
confidence regions and tests for high-dimensional models. The Annals of Statistics, 42(3):1166–
1202.
Verbeke, G. and Lesaffre, E. (1996). A linear mixed-effects model with heterogeneity in the randomeffects population. Journal of the American Statistical Association, 91(433):217–221.
Verbeke, G. and Molenberghs, G. (2009). Linear Mixed Models for Longitudinal Data. Springer,
corrected edition.
Vershynin, R. (2010). Introduction to the non-asymptotic analysis of random matrices.
Vitreschak, A. G., Rodionov, D. A., Mironov, A. A., and Gelfand, M. S. (2002). Regulation of
riboflavin biosynthesis and transport genes in bacteria by transcriptional and translational attenuation. Nucleic acids research, 30(14):3141–3151.
Vogl, C., Grill, S., Schilling, O., Stülke, J., Mack, M., and Stolz, J. (2007). Characterization of
riboflavin (vitamin b2) transport proteins from bacillus subtilis and corynebacterium glutamicum.
Journal of bacteriology, 189(20):7367–7375.
Wang, L., Zhou, J., and Qu, A. (2012). Penalized generalized estimating equations for highdimensional longitudinal data analysis. Biometrics, 68(2):353–360.
Zhang, C.-H. and Zhang, S. S. (2014). Confidence intervals for low dimensional parameters in high
dimensional linear models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(1):217–242.
Zhang, D. and Davidian, M. (2001). Linear mixed models with flexible distributions of random effects
for longitudinal data. Biometrics, 57(3):795–802.
Zhu, Y. and Bradic, J. (2016). Linear hypothesis testing in dense high-dimensional linear models. to
appear in the Journal of the American Statistical Association: Theory and Methods.

38

