Estimating the Random Effect in Big Data Mixed Models
Michael Law Ya'acov Ritov
math.ST stat.ME stat.TH

	
	
	
	
	
		
	Estimating the Random Effect in Big Data Mixed Models
    Michael Law Ya'acov Ritov Supported in part by NSF Grant DMS-1712962.
    
		University of Michigan

		December 30, 2023
=========================================================================


		
		We consider three problems in high-dimensional Gaussian linear mixed models.  Without any assumptions on the design for the fixed effects, we construct an asymptotic F-statistic for testing whether a collection of random effects is zero, derive an asymptotic confidence interval for a single random effect at the parametric rate √(n), and propose an empirical Bayes estimator for a part of the mean vector in ANOVA type models that performs asymptotically as well as the oracle Bayes estimator.  We support our results with numerical simulations and provide comparisons with oracle estimators.  The procedures developed are applied to the Trends in International Mathematics and Sciences Study (TIMSS) data.
		
	
	
	
	

§ INTRODUCTION

	In the past two decades, there has been a lot of progress in the theory for high-dimensional linear models.  However, its close cousin, the high-dimensional linear mixed model, has received significantly less attention; it was not until the past decade until there were procedures for estimation.  Consider a linear mixed model given by
	
    Y = μ + Zν + Wγ + ϵ,

	with Z ∈^n× v, W∈^n× r, and Y,μ,ϵ∈^n.  In addition, we observe covariates X∈^n× p such that μ - Xβ is a small remainder term for some sparse vector β∈^p (see Section <ref> for a rigorous definition).  Therefore, the model may be written as
	
    Y = Xβ + Zν + Wγ +  + ϵ.

	Here, X is a design corresponding to the fixed effects and (Z,W) to the random effects.  We consider the setting where the random effects are low-dimensional, v+r < n, but the fixed effects are high-dimensional, p>n.  We have separated the random effects in two to emphasize later that we are interested in ν and view γ as nuisance parameters.  Various authors have considered different aspects of this problem.  
	
	The earliest work of schelldorfer2011 proposed an estimator for both β and the variance components using a lasso-type approach.  These types of approaches were later extended by several authors who considered estimation with both convex, such as groll2014, and non-convex penalties, such as wang2012.  There has also a growing literature on model selection in high-dimensional linear mixed models (cf. muller2013 and the references therein).  
	
	The problem of inference is slightly less well studied.  To the best of our knowledge, chen2015 was the first to consider hypotheses testing problems for random effects in ANOVA type settings and bradic2017 for fixed effects in high-dimensional models.  During the preparation of this manuscript, we became aware of the very recent work of li2019, who consider the problem of inference in high-dimensional linear mixed models.  In particular, they discuss inference for fixed effects and estimation of variance components.
	
	The goal of the present paper is to contribute to this growing literature on high-dimensional linear mixed models, both in terms of estimation and inference when the random effects and error are all Gaussian.  In particular, we consider three related problems:
	
	
		
  * Testing whether ν≡ 0.
		
  * Constructing confidence intervals for σ_ν^2 when ν∼_v( 0_v , σ_ν^2 I_v ).
		
  * Estimating using empirical Bayes in ANOVA Type Models.
	
		
	Our methodology is inspired by both low-dimensional linear mixed models as well as high-dimensional linear models.  Specifically, our approach to all three problems starts with considering a procedure in the corresponding low-dimensional problem and retrofitting it with tools and techniques from high-dimensional linear models to produce a procedure for the high-dimensional linear mixed model.
	
	

 §.§ Organization of the Paper

	We will end the current section with a description of the notation that we will adopt throughout the paper.  Sections <ref>, <ref>, and <ref> consider each of the three problems outlined in the Introduction respectively.  Each one starts with a description of the problem setup, a brief motivation from the low-dimensional problem, a description of the estimator that will be considered, and ends with some theoretical results.  In Section <ref>, we provide an overview of the computation of the estimators as well as simulations and a real data application.  We conclude with discussions in Section <ref>.  For the ease of presentation, we defer all proofs to Section <ref>.
	
	

 §.§ Notation

	Throughout, all of our variables will have a dependence on n.  However, when this will not cause confusion, we will suppress this dependence.  Let ‖·‖ denote the standard (unscaled) Euclidean norm with the dimension of the space being implicit from the vector and ‖·‖_0 the L_0-“norm” on ^p.  We will also use ‖·‖ to denote the ℓ_2 norm for matrices and ‖·‖_HS to denote the Hilbert-Schmidt norm.  Furthermore, λ_max(·) and λ_min(·) will denote the maximal and minimal eigenvalue respectively.  
	
	Consistent with other high-dimensional works, we will assume that β is a sparse vector.  There are various notions of sparsity but we will assume the slightly more general setting of weak sparsity.  For u ∈, we will let _u denote the collection of all models with the dimension of the fixed effects design equal to u.  For a model m∈_u, X_m will denote the n× u sub-matrix of X corresponding to the columns indexed by m.  Moreover, P_m will denote the projection onto X_m and m the projection onto the orthogonal complement.  We will consider the following definition of weak sparsity.
	
	
		The vector μ is said to satisfy the weak sparsity property relative to X with sparsity s at rate k if the set 
		
    𝒮_μ{m∈_s : ‖mμ‖^2 = o(k) }

		is non-empty.  A set S∈𝒮_μ is said to be a weakly sparse set for the vector μ.
			
	
	Then, we will let S∈𝒮_μ denote any weakly sparse set for μ.
	
	

§ A HIGH-DIMENSIONAL F-TEST: TESTING Ν≡ 0


	

 §.§ Model and Motivation

	Consider the high-dimensional linear mixed model given in equation (<ref>).  For the fixed effects X, we will make no assumption on the design besides independence with (ν,γ,ϵ).  For the random effects, we will similarly make no assumption on (W,Z), except for the fact that (W,Z) is independent of (ν,γ,ϵ).  We will assume that ϵ∼_n(0_n,σ_ϵ^2 I_n) for some positive constant σ_ϵ^2>0 and ν∼_v(0_v,Ψ) for some symmetric positive semi-definite matrix Ψ.  No assumption will be necessary on γ as the nuisance parameters will be projected out in the first stage.
	
	We are interested in the hypotheses testing problem
	
    H_0:  λ_max(Ψ) = 0           H_1:  λ_max (Ψ) > 0.

	
	Suppose temporarily that we are in the low-dimensional setting with s = p, p + v + r < n, and ≡ 0.  Then, in this low-dimensional problem, with ν and ϵ both distributed as Gaussian, the standard procedure for testing whether ν≡ 0 is through the F-test.  Letting Ã denote the projection onto Z ⊖ (X,W) with n_ã = (Ã) and B̃ denote the projection onto (X,Z,W)^⊥ with n_b̃=(B̃), we have that 
	
    F = ‖Ã Y ‖^2/n_ã/‖B̃ Y ‖^2/n_b̃.

	
	Then, under the null hypothesis, the above statistic has an F_n_ã,n_b̃ distribution.  It is known that the F-test has certain optimality properties for certain classes of hypotheses testing problems (cf jiang2007 and the references therein for more details).
	
	The main obstacle to directly using this F-test in the high-dimensional setting is removing the contribution of the fixed effects.  One possibility is to perform model selection and choose the relevant covariates from X and then use the low-dimensional F-test.  chen2015 consider a similar problem in the high-dimensional setting and they use a SCAD based approach for variable selection.  As a consequence, they require p=o(√(n)).  Instead, we leverage the fact that a projection onto a particular space is a regression onto a design whose columns span that same space.
	
	Expanding both the numerator and the denominator of the classical F-statistic, we have that
	
    ÃY = ÃZν + Ãϵ
       B̃Y = B̃ϵ.

	
	In both matrices above, they project onto the orthogonal complement of W, which may still be achieved in the high-dimensional problem since W is a low-dimensional matrix.  We may find two matrices, A and B, such that
	
    AY = AXβ + AZν + Aϵ,
       BY = BXβ + Bϵ,

	where A ∈^n_a× n and B ∈^n_b× n have orthonormal rows that are mutually orthogonal (ie. AB^ = 0_n_a× n_b).  For simplicity, we will write 
	
    Q [ A; B ].

	Rewriting, we have that
	
    [ AY; BY ] = QY = QXβ + QZν + Qϵ = [ AXβ + AZν + Aϵ;       BXβ + Bϵ ].

	If QX was low-dimensional, to obtain the projection of QY onto the orthogonal complement of QX, this is equivalent to finding the residuals of QY using the covariates QX; this yields QY - QXβ̂.  Then, we have that
	
    [ AY - AXβ̂; BY - BXβ̂ ]   = QY - QXβ̂ = QXβ - QXβ̂ + QZν + Qϵ
       = [ AXβ - AXβ̂; BXβ - BXβ̂ ] + [ AZν + Aϵ;       Bϵ ].

	
	Hence, this recasts the problem into one of high-dimensional prediction, for which there have been many procedures suggested to estimate QXβ, such as the lasso and exponential weighting (cf. tibshirani1996 and leung2006 respectively).  Therefore, we propose using a plug-in estimator for QXβ using exponential weighting of all models of a particular size and then consider the resultant residuals.  This idea, under mild assumptions, will provide an asymptotic F-test.
	
	

 §.§ Estimator

	
	Rather than using exponential weighting for AY and BY separately, we will estimate the mean vector jointly; the advantages to this approach are twofold.  On the one hand, it will significantly simplify our implementation as we will only need a single Markov chain (see Section <ref>) and, on the other, the estimation seems to perform better empirically since we pool more observations together.  Therefore, we will fix a sequence of constants u∈.  Let  denote the least-squares estimator of β using the model m∈_u with covariates QX_m.  We will define our set of exponential weights by
	
    w_m,Qexp( -1/α‖ QY - QX‖^2 )/∑_k∈_uexp( -1/α‖ QY - QXβ̂_k,Q‖^2 ),

	where α≥ 6λ_max( QZZ^ Q^ + σ_ϵ^2 I_n_a+n_b).  We will estimate β by 
	
    ∑_m∈_u w_m,Q.

	The bound on α is to ensure the consistency result of Lemma <ref> under both the null and the alternative hypotheses.  The subscript Q is to emphasize that the estimator is constructed only using the Q projected sub-data from (<ref>).  Then, the corresponding F-statistic will be
	
    ‖ AY - AX‖^2/n_a/‖ BY - BX‖^2/n_b.

	Similar to the classical F-statistic, we will reject the null hypothesis for large values of .  In particular, for a value δ∈(0,1), let F_n_a,n_b,δ denote the δ upper quantile of an F_n_a,n_b distribution.  Then, we will consider tests of the form
	
    1( > F_n_a,n_b,δ).

	
	

 §.§ Assumptions

	We will assume without the loss of generality that the columns of X have squared norm that is 𝒪(n).  Moreover, we will assume that
	
		
  *  The vector ν satisfies
		
    ν∼_v(0_v,Ψ).

		Furthermore, λ_max( QZ Ψ Z^ Q^ + σ_ϵ^2 I_n ) = 𝒪(1).
		
  *  The sequence of constants u satisfies
		
    lim inf_n→∞ (u-s) ≥ 0,
       lim_n→∞ulog(p)/n_a ∧ n_b→ 0.

	
	
		
		
  *  The sequence of constants u and the number of observations in the reduced models, n_a and n_b, satisfy
		
    lim inf_n→∞ (u-s) ≥ 0,
       lim_n→∞ulog(p)/√(n) = 0,
       n_a ≍ n_b ≍ n.

	
	
	
		There are two assumptions regarding the sparsity and the number of random effects.  The first assumption, <ref>, is quite weak in terms of the sparsity level and is used to establish the distribution of the estimator under the null hypothesis.  Considering the models in equations (<ref>) and (<ref>), which have n_a and n_b observations respectively, the sparsity only needs to be of smaller order than the number of observations in each of the two reduced models.  This is known to be the optimal rate at which the mean vector can be estimated in a high-dimensional linear model (for example, cf. raskutti2011 or rigollet2011).  The requirement that lim inf (u-s)≥ 0 is only to ensure that the models we consider for exponential weighting are sufficiently large to remove the fixed effects. 
		
		The stronger assumption, <ref>, is used when considering the power of the test under a specific contiguous alternative (see equation (<ref>)), which gives the parametric rate of √(n).  However, in cases where we have a different sequence of alternatives, the assumption may be weakened accordingly.  Note that this is similar to the setting of the de-biased lasso, which requires a ultra sparsity to achieve the parametric rate in testing a single covariate (cf. cai2017).  
	
	
	

 §.§ Main Results for 

	We start by stating the asymptotic distribution of  under the null hypothesis.
	
		Consider the linear mixed model given by equation (<ref>) and the hypotheses testing problem from equation (<ref>).  Assume <ref>.  Under the null hypothesis and α≥ 4σ_ϵ^2,
		
    = F + (1),

		where F∼ F_n_a,n_b.
	
	
	
	This naturally leads to the question regarding the power of the testing procedure.  In particular, between what contiguous alternatives can we distinguish?  We will consider the following contiguous testing problem.
	
    H_0:  λ_max(AZΨ Z^ A^) = 0           H_1:  λ_min(AZΨ Z^ A^) = h/√(n).

	
	
		Suppose that ν corresponds to a single random effect and the design is balanced, with ν∼_v(0_v,σ_ν^2 I_v).  Then, the above hypotheses becomes
		
    H_0:  σ_ν^2 = 0           H_1:  σ_ν^2 = h/√(n),

		which is a standard hypotheses testing problem, such as in the balanced one-way random effects model.  In this model, in the low-dimensional setting, the rate of √(n) is optimal.
	
	
	
		Consider the linear mixed model given by equation (<ref>) and the hypotheses testing problem from equation (<ref>).  Assume further <ref> and <ref>.  Fix a value of δ>0.  Under the alternative hypothesis with h>0 sufficiently large (not depending on n) and α≥ 6λ_max( QZΨ Z^ Q^ + σ_ϵ^2 I_n_a+n_b), the sum of type I and type II error for the test statistic  is between zero and one.
	
	
	To prove the theorem, we will need the following lemma that provides consistency of  under the setting of correlated errors, which may be of interest in its own right.
	
	
		Consider a high-dimensional linear model given by
		
    Y = μ + ξ,

		with ξ∼_n(0_n,Σ) for some covariance matrix Σ satisfying λ = λ_max(Σ) = 𝒪(1).  Assume that μ is weakly sparse relative to X with sparsity s at rate n and ‖ Xβ‖^2 = 𝒪(n).  Assume further log(|_u|) = o(n).  Letting  denote the least-squares estimator for β using the covariates X_m, define the exponential weights as
		
    w_m exp( -1/α‖ Y - X‖^2 )/∑_k∈_uexp( -1/α‖ Y - Xβ̂_k ‖^2 ),

		with α > 6λ.  Then,
		
    1/n‖∑_m∈_u w_m X - μ‖^2 → 0.

		If in addition μ is weakly sparse relative to X with sparsity s at rate √(n) and log(|_u|) = o(√(n)), then
		
    1/√(n)‖∑_m∈_u w_m X - μ‖^2 → 0.

	
	
	
		The above lemma slightly generalizes Theorem 5 of leung2006 to heteroskedastic and dependent, but still Gaussian, errors.  However, in the setting where the errors are independent and identically distributed, the result of leung2006 is stronger as their bound for α is smaller.
	
	
	
	

§ CONFIDENCE INTERVALS FOR Σ_Ν^2

	
	

 §.§ Model and Motivation

	
	
	In the previous section, we considered the problem of testing a collections of random effects.  However, it is often of interest to construct confidence intervals for the variance of a particular random effect.  Suppose that Ψ = σ_ν^2 I_v.  In the low-dimensional setting, there have been many procedures suggested to construct confidence intervals, from likelihood based approaches to F-test inversions (for example, see jiang2007 and the references therein for a non-exhaustive list).
	
	We take an approach similar to variance estimation in linear models.  Consider the two reduced models from Section <ref>, which are reproduced below for convenience:
	
    AY = AXβ + AZν + Aϵ,
       BY = BXβ + Bϵ,

	
	The second model is a standard high-dimensional linear model with Bϵ∼_n_b(0_n_b,σ_ϵ^2 I_n_b).  In the low-dimensional setting, σ_ϵ^2 would be estimated by mean squared error, which is computed through projecting BY onto the orthogonal complement of BX.  Hence, we may achieve a similar result by viewing this as a prediction problem and regressing out BX.  For the other model, note that we have a high-dimensional linear model with noise AZν + Aϵ∼_n_a( 0_n_a , σ_ν^2 AZZ^ A^ + σ_ϵ^2 I_n_a).  Since the matrix A and Z are both observed, we may whiten the data with respect to ν by taking the Spectral Decomposition of AZZ^ A^ = V D V^ and left multiplying by D^-1/2 V^.  Here, D is diagonal and V is unitary.  This yields
	
    D^-1/2V^ A Y = D^-1/2V^ AXβ + D^-1/2V^ AZν + D^-1/2V^ A ϵ .

	Then, the above can be viewed as a high-dimensional linear model with noise given by D^-1/2V^ AZν + D^-1/2V^ A ϵ∼_n_a( 0_n_a , σ_ν^2 I_n_a + σ_ϵ^2D^-1).  Since D^-1 is observed and σ_ϵ^2 can be estimated, this will provide an estimator of σ_ν^2, which we may then use to construct a pivot.
	
	

 §.§ Estimator

	As outlined previously, we will need to construct two estimators, one for σ_ϵ^2 and one for σ_ν^2.  Again, we will jointly estimate the mean vector Qμ.  Therefore, we will start by defining
	
    1/n_b‖ BY - BX‖^2.

	Then,  will be given by
	
    1/(D^-1)‖ D^-1/2V^ AY - D^-1/2 V^ AX‖^2 - .

	
	

 §.§ Assumptions

	In addition to the assumptions from Section <ref>, we will need two additional assumptions regarding the distribution of the values of the diagonal matrix D.
	
	
		
  *  Letting d=(D), the norms satisfy
		
    lim sup_n→∞‖σ_ν^2 1_n_a + σ_ϵ^2 d ‖_∞/‖σ_ν^2 1_n_a + σ_ϵ^2 d ‖_2→ 0.


		
  *  The sample sizes n_a and n_b satisfy
		
    lim_n→∞n‖σ_ν^2 1_n_a + σ_ϵ^2 d ‖_2^2/((D^-1))^2→σ_a^2>0,
       lim_n→∞2σ_ϵ^4n/n_b→σ_b^2>0,
       lim_n→∞n_a/(D^-1)→ > 0.

	
	
	
		Assumption <ref> on the matrix D is to ensure that the design is sufficiently well balanced and that Z is distinct enough from W.  This assumption is satisfied in general settings such as relatively balanced ANOVA models.

		Assumption <ref> is a technical regularity condition in view of the triangular array framework.  By possibly reducing to a convergent sub-sequence, we may always satisfy assumption <ref>.
	
	
	

 §.§ Main Results for 

	We will start by stating the asymptotic distribution of .
	
	
		Consider the linear mixed model given by equation (<ref>), with the reduced models given in equations (<ref>) and (<ref>).  Assume further <ref> with Ψ = σ_ν^2 I_v, <ref>, <ref>, and <ref>.  Then,
		
    √(n)(  - σ_ν^2 ) (0,σ_a^2 + σ_b^2).

	
	
	From the previous theorem, we immediately have that the following is a (1-α) confidence interval for σ_ν^2.
	
    ( 1/(  - √(σ_a^2+σ_b^2)/√(n) z_α/2) , 1/(  + √(σ_a^2+σ_b^2)/√(n) z_α/2) ).

	
	In general, the values of σ_a^2, σ_b^2, and  are unknown quantities and need to be estimated.  The following proposition yields consistent estimators for each of those quantities.
	
		Under the assumptions of Theorem <ref>,
		
    2n∑_k=1^n_a( /d̂ + d_k^-1)^2 /((D^-1))^2σ_a^2,
       2n()^2/n_bσ_b^2,
       d̂n_a/(D^-1).

	
	
	Therefore, we may plug in the values of , , and  into equation (<ref>) to obtain an asymptotic confidence interval.
	
	

§ EMPIRICAL BAYES IN ANOVA TYPE MODELS

	
	The motivating example of this problem framework is in terms of the Rasch model, originally proposed by rasch1960.  The model that we consider is slightly different than the classical Rasch model in that we have Gaussian responses as opposed to binary responses.  
	
	As an example of this model, the data that we will consider in Section <ref> is from the Trends in Mathematics and Sciences Study (TIMSS), an international study conducted every four years to measure fourth and eighth grade student achievement in mathematics and science.  We will only consider data from the year 2015.  Countries randomly sample a collection of nationally representative schools to take standardized examinations in both mathematics and science, with questions being either multiple choice or constructed response.  Then, each student within schools takes only a subset of the questions on the exams but all questions are answered by some students in each school.  In addition to recording student responses, the data also contains background covariates for schools.  martin2016 provides a more detailed description of the methods and procedures employed by TIMSS and more general information about TIMSS is available in mullis20years.
	
	For our analysis, we will only consider multiple choice questions and analyze on the level of school rather than students.  To construct a response variable for school, we compute the proportion of questions answered correctly by students in that school.  Note that, unlike the classical Rasch model, we assume a linear model and, for all schools, we have answers for all questions.  Thus, by the Central Limit Theorem, our response Y will be approximately Gaussian.  The fixed effects design X will be the background covariates for the school and the random effects design Z will be an indicator for the country, with ν corresponding to the unobserved variability of the countries.  In this example, we do since we have averaged over questions, we do not have any nuisance random effects.  The problem that we will consider in this section is ranking the countries based on mathematical ability and trying to estimate the average number of questions that any particular country will answer correctly.  That is, we would like to estimate μ+Zν for all countries.
	
	

 §.§ Model and Motivation

	The general problem framework that we consider is for K-factor ANOVA models.  However, we will derive the results in the setting when K=2.  That is, we will consider the model
	
    Y = μ + Zν + Wγ + ϵ.

	We do not assume that the design is fully crossed in the random effects.  The goal in the problem is to estimate a subset of the mean vector, ημ + Zν, since we view the random effects W as nuisance.  However, as the sample size increases, the number of observations per group stays bounded.  In the context of the motivating data example, each country still only answers a finite number of questions as we increase the sample size.  Therefore, it is not possible to consistently estimate the entire vector η.  A standard approach in the low-dimensional setting would be to use an empirical Bayes estimator by placing a Gaussian prior on both ν and γ (for example, see brown2018), which would then transform the problem into a standard high-dimensional linear mixed model.  Therefore, we will use a _v(0_v,σ_ν^2 I_v) and _r(0_r , σ_γ^2 I_r) prior on ν and γ respectively.
	
	Since we will need to estimate both σ_ν^2 and σ_γ^2 for the prior, our estimator for σ_γ^2 will be analogous to  from Section <ref>.  To this end, we will need an additional matrix C such that 
	
    CY = CXβ + CWγ + Cϵ.

	Further, C may be chosen such that C has orthonormal rows that are mutually orthogonal with (A^,B^)^.  We will let n_c denote the number of rows of C.  Similarly, we will need to consider the Spectral Decomposition of CWW^ C^ = ΓΛΓ^.  Similarly, we have that Λ is diagonal and Γ is unitary.  Then, we may take the transformation
	
    Λ^-1/2Γ^ CY = Λ^-1/2Γ^ Xβ + Λ^-1/2Γ^ CWγ + Λ^-1/2Γ^ Cϵ

	to whiten the data with respect to γ.
	
	

 §.§ Estimator

	To estimate μ, we will define  to be the exponentially weighted estimator for model (<ref>), with μ̂_EW X.
	
	From Theorem <ref>, we immediately have that σ_ν^2.  This suggests the empirical prior _v( 0_v , /d̂) for ν.  However, we will not use this estimator, but a slightly modified version by jointly estimating the vector μ using all of the observations.  Before we define a modified version of , we will first need to adapt our estimator of σ_ϵ^2, which will be given by
	
    1/n_b‖ BY - BX‖^2.

	Then, we may define
	
    1/(D^-1)‖ D^-1/2V^ AY - D^-1/2 V^ AX‖^2 - .

	We may similarly define  as
	
    1/(Λ^-1)‖Λ^-1/2Γ^ CY - Λ^-1/2Γ^ CX‖^2 - .

	Then, under some assumptions, it will follow that σ_γ^2 analogous to Theorem <ref> where  is defined in assumption <ref>.  Hence, we will estimate σ_ν^2 by /d̂, σ_γ^2 by /λ̂, and σ_ϵ^2 by .  This will give the following empirical Bayes estimator for η,
	
    μ̂_EW + /d̂ ZZ^( /d̂ ZZ^ + /λ̂WW^ +  I_n )^-1( Y - μ̂_EW).


	To compare our estimator, we will consider an oracle that has access to μ, σ_ν^2, σ_γ^2, and σ_ϵ^2.  Then, this oracle will use the Bayes estimator for η (see Lemma <ref>), given by
	
    μ + σ_ν^2 ZZ^( σ_ν^2 ZZ^ + σ_γ^2 WW^ + σ_ϵ^2 I_n )^-1( Y - μ).

	
	
	

 §.§ Assumptions

	Unlike the previous section, we will not need to establish the asymptotic distribution of , rather we only need the estimator to be consistent.  Accordingly, we may weaken our assumptions to the following
	
	
		
  *  The matrices D and Λ satisfy
		
    lim_n→∞( D^-2)/( D^-1)^2→ 0          lim_n→∞( Λ^-2)/( Λ^-1)^2→ 0
       lim_n→∞n_a/( D^-1)→         lim_n→∞n_c/( Λ^-1)→.

		
		
  *  The sequence of constants u satisfies
		
    lim inf_n→∞ (u-s) ≥ 0,
       lim_n→∞ulog(p)/n_a ∧ n_b ∧ n_c→ 0.

		
		
  *  The design matrices satisfy
		
    ‖σ_ν^2 ZZ^ + σ_γ^2 WW^ + σ_ϵ^2 I_n ‖^2 = 𝒪(1).

	
	
	
		Assumption <ref> is weaker than the combination of assumptions <ref> and <ref>.  This will ensure consistency in estimating σ_ν^2, σ_γ^2, and σ_ϵ^2.
		
		The next assumption, <ref>, is identical to assumption <ref> but with the additional constraint on n_c to ensure consistent prediction of Λ^-1/2Γ^ Cμ in equation (<ref>).
		
		The last assumption, <ref>, is equivalent to, in the context of the motivating example, each school only answering 𝒪(1) number of questions and each question being answered by 𝒪(1) number of schools.
	
	
	
	

 §.§ Main Results for 

	We will start this section by proving that , , and  are all consistent estimators.
	
		Consider the linear mixed model given by equation (<ref>), with the reduced models given in equations (<ref>), (<ref>), and (<ref>).  Then, assuming <ref> and <ref>,
		
    σ_ν^2,
    			         σ_γ^2,
    			      σ_ϵ^2.

	
	
	The following is a standard lemma regarding the empirical Bayes estimators in this problem setup, which we will prove for the sake of completeness.
	
		For a fixed vector μ∈^n and fixed values σ_ν^2>0, σ_γ^2>0, and σ_ϵ^2 >0, the Bayes estimator of η is given by
		
    ( η | Y ) = μ + σ_ν^2 ZZ^( σ_ν^2 ZZ^ + σ_γ^2 WW^ + σ_ϵ^2 I_n )^-1( Y - μ).

	
	
	We conclude this section with the main result regarding , that the empirical Bayes estimator performs nearly as well as the oracle Bayes estimator  asymptotically.
	
		Consider the linear mixed model given in (<ref>).  Assume <ref>, <ref>, and <ref>.  Then,
		
    1/n( ‖ - η‖^2 - ‖ - η‖^2 ) = (1).

	
	
	

§ IMPLEMENTATION AND SIMULATIONS

	

 §.§ Implementation

	
	In this section, we describe how to compute all three of the proposed estimators.  The computation can be split into two halves:
	
		
  * Obtaining the matrices A, B, C, D, V, Λ, and Γ with which to project the data.
		
  * Computing the estimators via exponential weighting.
	
	
	For the first part, the matrices A, B, and C can be obtained by adding the n standard basis vectors in ^n to (W,Z) to form a basis for ^n and applying the modified Gram-Schmidt procedure.  To compute D and V, we form the matrices AZZ^ A^ and directly apply the eigendecomposition to obtain both the eigenvalues and the eigenvectors.  An analogous procedure is used to compute Λ and Γ.
	
	For the second part, it may be further sub-divided into three parts:
	
		
  * Determining a bound α.
		
  * Determining a value of u.
		
  * Computing |_u| = p  u = 𝒪(p^u) models.
	
	To determine α, we only need an upper bound on the maximal eigenvalue of the noise vector.  This may simply be bounded by the variance of Y by choosing α = 4 ‖ Y ‖^2/n.  Though this will choose a conservative value of α, asymptotically this will have no effect.  For u, we may following the suggestion of rigollet2011 and apply the exponential screening prior.  They say that a covariate X_j is selected if β̂_ES,j > 1/n, where β̂_ES,j is the j'th entry of the exponential screening estimator.  Note that, after computing the matrix B, the model BY = BXβ + Bϵ is a standard high-dimensional linear model with the same level of sparsity as the original model, which is in the framework of rigollet2011.  Alternatively, we may apply cross-validation to tune the sparsity level, though this is very computationally intensive.
	
	Once both are selected, we still need to compute p  u least-squares estimators.  However, rigollet2011 and rigollet2012 proposed a Metropolis Hastings approach to approximating the exponential weighting estimator.  We will refer the interested reader to rigollet2012 for the details.

	

 §.§ Simulations

	

  §.§.§ Methods and Models

	
	We will consider the linear mixed model given by
	
    Y = Xβ + Zν + Wγ + ϵ,

	with n = 200, s=3, and p=500.  For simplicity, we will set μ = Xβ.  Since no assumption is made on the fixed effects design X, we will generate X from a multivariate Gaussian _n(0_n, Σ), where Σ is an equi-correlation matrix.  That is,
	
    Σ_i,j = 
    		
    			1  if  i=j
    ρ if  i≠ j

	for ρ∈{0,0.8}.  As our fixed effects design is perfectly symmetric, we will take the strongly sparse set S to be the first s components with β_j = 1 for j∈ S.  We will fix the value of σ_ϵ^2=1.  For simplicity, we will let ν and γ each correspond to a single random effect that are independent of each other.  Therefore, we will have that
	
    [ ν; γ; ϵ ]∼_v+r+n( 
    		[ 0_v; 0_r; 0_n ] , 
    		[ σ_ν^2 I_v    0_v× r    0_v× n;    0_r× v σ_γ^2 I_r    0_r× n;    0_n× v    0_n× r σ_ϵ^2 I_n ]).

	Finally, with regards to the random effects design, we will consider the setting where (v,r)∈{25,50}^2.  Since the number of observations will be less than the number of possible combinations, we will start by generating a fully crossed design and then down sample n rows without replacement of the fully crossed design.  To experiment with both the level and the power of our procedure, we will consider (σ_ν^2,σ_γ^2) ∈{0,1}^2.  Hence, we will consider each of the following settings
	
    (ρ,v,r,σ_ν^2,σ_γ^2) ∈{0,0.8}×{25,50}×{25,50}×{0,1}×{0,1}.

	
	For each setting, we will generate a single draw of the design for (X,W,Z).  Then, we will draw N=100 trials of (ν,γ,ϵ) and compute each of our estimators on these trials.  Whenever a variance estimate is negative, we will set the estimated value to zero.  When using exponential weighting, we will estimate the sparsity using the exponential screening prior of rigollet2011 on the model BY = BXβ + Bϵ.  Moreover, we will use the value of α = 4‖ Y ‖^2/n.
	
	All of our simulations are conducted in R.  For each of our three problems, we will compare our estimator with an oracle low-dimensional estimator as well as a low-dimensional version of our proposed high-dimensional statistic.  For the oracle estimators, in the setting of the F-test, we directly apply the classical low-dimensional F-test that has access to the true sparse set S, as given in equation (2.3) of jiang2007, which we will denote by .  For the confidence intervals, we fit the linear mixed models with the true sparse set S using lmer and applying the confint function.  We will denote this estimator as .  Finally, in the setting of estimation, we directly compute the oracle Bayes estimator  described in Section <ref>.  Collectively, things subscripted with “LD” will denote low-dimensional estimators.
	
	In addition to comparing with the low-dimensional estimators, we also construct low-dimensional versions of our proposed high-dimensional statistics.  To do so, we use the exact same statistic as in the high-dimensional setting but replace exponential weighting with least-squares using the sparse set S.  Hence, all of these statistics will be subscripted by “LS” to denote the usage of least-squares.  We make this comparison since all of our proposed statistics rely on two layers of asymptotics:  
	
		
  * In the prediction of the mean vector via exponential weighting.
		
  * In the convergence once the residuals are obtained.
	
	To differentiate between these two, we introduce an intermediate statistic that relies on least-squares, which we think of as low-dimensional versions of our statistics.  For example, recalling that β̂_S is the least-squares estimator of β using the covariates X_S, we will also consider the statistic
	
    ‖ AY - AXβ̂_S‖^2/n_a/‖ BY - BXβ̂_S‖^2/n_b∼ F_n_a,n_b + (1).

	Similarly, we will consider  and  given by
	
    1/n_b‖ Bϵ - BXβ̂_S‖^2,
       1/(D^-1)‖ D^-1/2V^ AY - D^-1/2 V^ AXβ̂_S ‖^2 - ,
       1/(Λ^-1)‖Λ^-1/2Γ^ CY - Λ^-1/2Γ^ CXβ̂_S ‖^2 - .

	Then, we will have that
	
    √(n)(  - σ_ν^2 ) (0,σ_a^2 + σ_b^2)

	and
	
    Xβ̂_S + /d̂ ZZ^( /d̂ ZZ^ + /λ̂WW^ +  I_n )^-1( Y - Xβ̂_S ).

	To compare the three procedures, we will consider the following metrics
	
		
  * Average Coverage:  The percentage of time the correct hypothesis is selected for F-tests or the percentage of time the true value of σ_ν^2 is in the confidence interval.
		
  * Average Length:  The average length of the confidence interval, taken as the upper endpoint minus the lower endpoint.
		
  * Average Loss:  The average squared Euclidean distance between the estimated vector η̂ and the true vector η divided by n.
	
	
	

  §.§.§ Results

	Tables (<ref>) - (<ref>) provide the results of our simulations.  We notice that, on average,  and  are slightly less conservative than  and  respectively.  However, as v increases from 25 to 50, we generally notice that type I error decreases and coverage increases, suggesting that the main problem that is asymptotic.  In particular, even with the correct low-dimensional covariates,  and  still have higher type I error and lower coverage but does outperform both  and , as expected.  However, the confidence intervals for  seem to be overly optimistic when the design is more correlated and the true value of σ_ν^2=1.  In regards to estimation of η, the loss for  is generally close to the oracle loss, as predicted by Theorem <ref>.
	
	Moreover, Figures (<ref>) and (<ref>) show some simulated confidence intervals for 10 trials.  The fact that  and  occasionally have shorter confidence intervals than profile likelihood might be attributable to the more optimistic nature of those two procedures.
	
	

	
	
	
	
	
	
		
	
	
	
	

 §.§ Real Data Application

	
	Following in the motivating example of Section <ref>, we consider the Trends in International Mathematics and Sciences Study dataset, which is freely available online.  To simplify our analysis, we only consider the mathematics questions.  After filtering out for complete cases on background covariates, we are left with 146 questions, r=43 unique countries, p=106 covariates, and 6808 schools.  Therefore, we had a total of n=6808 responses after averaging over the students and questions within the schools.  
	
	To demonstrate our methodology, we apply all three procedures on the data.  When applying our methodology, we tune the sparsity using the exponential screening prior of rigollet2011 and set α = 4‖ Y ‖^2/n.  The high-dimensional F-test rejected the null hypothesis that σ_ν^2=0 and a 95% confidence interval for σ_ν^2 is (0.0020,0.0051), which suggests that, even controlling for school background characteristics, the country of the school impacts mathematical ability.  For the last part, we define a country's background characteristics X to be the arithmetic average of all the schools' background characteristics within that country.  Then, applying the empirical Bayes procedure, we rank the countries based on the predicted number of questions they would answer correctly.  The top five countries in order from our analysis are South Korea, Hong Kong, Singapore, Chinese Taipei, and Japan.  Up to some reordering, our results are mostly consistent with the report of mullis2016 based on individual student data, who had the same top five countries.
	
	
	

§ DISCUSSIONS

	In this paper, we considered three problems related to high-dimensional linear mixed models.  Throughout, we used exponential weighting to perform high-dimensional prediction due to its nice theoretical properties that do not require any assumption on the design matrix.  However, exponential weighting may be replaced with another procedure as long as we also satisfy the conditions for an oracle inequality for prediction, such as assuming compatibility for the lasso.  The theory for all three problems will remain unchanged.  In particular, the lasso has the advantage of being a convex optimization problem which is easily implemented by standard software.  
	
	The main drawback to our approach is the assumption that both the random effects and the error distributions are Gaussian, which is used to maintain independence after constructing various projections.  However, there are many situations when this is not a reasonable assumption and other procedures will need to be developed.  We view our contribution as a step towards these problems, similar to the development of the theory for linear mixed models in the low-dimensional setting.
	
	

§ PROOFS

	

 §.§ Proofs for Section <ref>

	
	
		Under the null hypothesis, we have that
		
    QY = Qμ + Qϵ,

		which is a standard high-dimensional linear model.  By assumption <ref> and Theorem 5 of leung2006,
		
    ‖ Aμ - AX‖^2    ≤‖ Qμ - QX‖^2 
       ≤‖ Q ‖^2 ‖μ - X‖^2 = (n_a).

		By the triangle inequality, we have that
		
    1/√(n_a)   ( ‖ Aϵ‖ - ‖ Aμ - AX‖)
    ≤   1/√(n_a)‖ AY - AX‖
       ≤1/√(n_a)( ‖ Aϵ‖ + ‖ Aμ - AX‖).

		By equations (<ref>) and (<ref>),
		
    1/√(n_a)‖ AY - AX‖ = 1/√(n_a)‖ Aϵ‖ + (1).

		Squaring both sides yields
		
    1/n_a‖ AY - AX‖^2 = 1/n_a‖ Aϵ‖^2 + (1).

		Similarly, for the denominator,
		
    1/n_b‖ BY - BX‖^2 = 1/n_b‖ Bϵ‖^2 + (1).

		Noting that the matrices A and B have orthonormal rows and are mutually orthogonal, independence of Aϵ and Bϵ follows from fact that ϵ is Gaussian.  Dividing finishes the proof.
	
	
	
		Temporarily, we will start by considering the numerator for .  Let Qξ QZν + Qϵ∼_v( 0_v,QZΨ Z^ Q^ + σ_ϵ^2 I_v ).  By an identical line of reasoning as in the proof of Theorem <ref>, it follows that
		
    1/n_a‖ AY - AX‖^2 = 1/n_a‖ Aξ‖^2 + (1/√(n)).

		We would like to remark that the (1/√(n)) as opposed to (1) in Theorem <ref> comes from assumption <ref> and we have applied Lemma <ref> of the present paper as opposed to Theorem 5 from leung2006.  Similarly, for the denominator,
		
    1/n_b‖ BY - BX‖^2 = 1/n_b‖ Bξ‖^2 + (1/√(n)) = 1/n_b‖ Bϵ‖^2 + (1/√(n)).

		Thus, this yields
		
    = ‖ Aξ‖^2/n_a/‖ Bϵ‖^2/n_b + (1/√(n))
       ≥λ_min(AZΨ Z^ A^) + σ_ϵ^2/σ_ϵ^2F_v,q + (1/√(n))
       ≥h/√(n) + σ_ϵ^2/σ_ϵ^2F_v,q + (1/√(n)),

		We would like to show that
		
    _H_0 + _H_1 = 
    			_H_0(  > F_n_a,n_b,δ) + _H_1( ≤ F_n_a,n_b,δ) < 1.

		It suffices to show that
		
    _H_1(  > F_n_a,n_b,δ) > δ.

		We will start by providing an upper bound on F_n_a,n_b,δ.  From Lemma 1 of laurent2000, it follows that, for a χ^2_d random variable,
		
    (χ^2_d > d + 2√(dx) + 2x) ≤exp(-x),
       (χ^2_d ≤ d - 2√(dy)) ≤exp(-y).

		Therefore, it follows that, for any x,y>0,
		
    ( F_n_a,n_b > 1 + 2√(x/n_a) + 2x/n_a/1 - 2√(y/n_b)) ≤exp(-x) + exp(-y).

		By choosing x,y>0 such that
		
    exp(-x) + exp(-y) ≤δ,

		then
		
    F_n_a,n_b,δ≤1 + 2√(x/n_a) + 2x/n_a/1 - 2√(y/n_b).

		Let a,b>0 be constants that will be chosen later.  Define the event  as
		
    {χ^2_n_b≤ n_b + 2√(n_bb) + 2b , χ^2_n_a > n_a - 2√(n_aa)}.

		Again, by Lemma 1 of laurent2000,
		
    ( ^) ≤exp(-a) + exp(-b).

		Now,
		
    _H_1(  > F_n_a,n_b,δ)    ≥_H_1(  > F_n_a,n_b,δ , ) 
       ≥_H_1( h/√(n) + σ_ϵ^2/σ_ϵ^2( 1 - 2√(a/n_a)/1 + 2√(b/n_b) + 2b/n_b)  
       + (1/√(n)) > 1 + 2√(x/n_a) + 2x/n_a/1 - 2√(y/n_b) , ).

		From assumption <ref>, we have that n_a ≍ n_b ≍ n.  Let
		
    lim inf_n→∞n_a/n = c_a,          lim inf_n→∞n_b/n = c_b.

		Then, assuming that
		
    h > 2σ_ϵ^2 ( √(bc_b) + √(ac_a) + √(yc_b) + √(xc_a)),

		it will follow that
		
    lim inf_n→∞_H_1(  > F_n_a,n_b,δ) ≥lim inf_n→∞() ≥ 1 - exp(-a) - exp(-b),

		which finishes the proof by tuning a,b>0 and h sufficiently large appropriately.
	
	
	
		Let P_m denote the projection onto X_m and m the projection onto the orthogonal complement.  Let r_m as
		
    r_m ‖mμ‖^2.

		We will choose an arbitrary sequence of weakly sparse sets S∈𝒮_μ satisfying
		
    r_S = o(n).

		Fix t>0 arbitrarily and define the set _t as follows
		
    _t { m∈_u : r_m ≤ tn }.

		Looking at the squared norm, we will start with the convexity of the norm and the triangle inequality to obtain that
		
    1/n‖∑_m∈_u w_m( P_mY - μ) ‖^2
    			≤   1/n∑_m∈_u w_m ‖ P_m Y - μ‖^2
    ≤   2/n∑_m∈_u w_m ( r_m + ‖ P_mξ‖^2 )
    ≤   2/n(∑_m∈_t w_m r_m + ∑_m∈_t^ w_m r_m ) 
       + 2/n∑_m∈_u w_m ‖ P_mξ‖^2.

		We will consider each of the three terms separately.  For the first term, it follows from the definition of _t that
		
    lim sup_n→∞2/n( ∑_m∈_t w_m r_m ) ≤ 2t.

		For the second term, we have that
		
    ( ∑_m∈_t^ w_m r_m )
    			= ∑_m∈_t^ r_m  w_m.

		Fix m∈_t^ temporarily.  Then, we have that
		
    w_m    ≤exp( -1/α( ‖m Y ‖^2 - ‖S Y ‖^2 ) )
       ≤exp( -1/α( r_m - r_S + 2ξ^mμ - 2ξ^Sμ + ‖ P_Sξ‖^2 - ‖ P_mξ‖^2 ) )
       ≤exp( -1/α( r_m - r_S + 2ξ^mμ - 2ξ^Sμ - ‖ P_mξ‖^2 ) ).

		Therefore, by the generalized Hölder inequality, 
		
    w_m ≤   exp(-1/α( r_m - r_S ) ) (exp( -6/αξ^mμ))^1/3
       ×(exp( 6/αξ^Sμ))^1/3(exp( 3/α‖ P_mξ‖^2 ) )^1/3.

		Computing each of the three Laplace transforms separately, since
		
    ξ^mμ∼(0 , ‖Σ^1/2mμ‖^2 ),

		it follows that
		
    (exp( -6/αξ^mμ) )^1/3 = exp( 6/α^2‖Σ^1/2mμ‖^2 ) 
    			≤exp( 6λ/α^2 r_m ).

		Similarly,
		
    (exp( 6/αξ^Sμ) )^1/3≤exp( 6λ/α^2 r_S ).

		For the other term, observe that
		
    ‖ P_mξ‖^2 ≼λχ^2_u,

		where ≼ denotes stochastic ordering.  Hence,
		
    exp( 3/α‖ P_mξ‖^2 )
    			≤exp( 3λ/αχ^2_u ) = ( 1 - 6λ/α)^-u/2.

		Since α > 6λ and m∈_t^, it will follow that
		
    w_m    ≤exp(-1/α(1 - 6λ/α) r_m + 6λ/α^2r_S - u/2log( 1 - 6λ/α) )
       ≤exp(- (1 - 6λ/α) tn/α + o(n) - u/2log( 1 - 6λ/α) ).

		By assumption, 
		
    lim sup_n→∞1/n‖mμ‖^2 ≤lim sup_n→∞1/n‖μ‖^2 ≤ C,

		for some constant C > 0.  Hence,
		
    lim sup_n→∞   2/n( ∑_m∈_t^ w_m r_m ) ≤ 2C lim sup_n→∞∑_m∈_t^ w_m
       ≤ 2Clim sup_n→∞exp(- (1 - 6λ/α) tn/α + o(n) - u/2log( 1 - 2λ/α) + log(|_t^|) ) 
       → 0,

		since log(|_u|) = o(n).  This implies that
		
    2/n( ∑_m∈_t^ w_m r_m ) → 0.

		For the last term, fix a value of K > 0 large to be determined later.  Define the event  as
		
    ⋂_m∈_u{‖ P_mξ‖^2 ≤ Klog(|_u|) + λ s }.

		We will first provide an upper bound on (^).  Indeed,
		
    ( ^) ≤∑_m∈_u( ‖ P_mξ‖^2 > Klog(|_u|) + λ s ).

		Let ϵ∼_n(0_n,I_n) and define Q_m Σ^1/2 P_mΣ^1/2.  Fixing m∈_u, we have by the Hanson-Wright inequality (Theorem 1.1 of rudelson2013) that
		
    ( | ‖ P_mξ‖^2 - ‖ P_mξ‖^2 | > Klog(|_u|) ) 
       =( | ϵ^ Q_m ϵ - ϵ^ Q_m ϵ| > Klog(|_u|) ) 
       ≤ 2 exp( -c min( K^2log^2(|_u|)/‖ Q_m ‖^2_HS , Klog(|_u|)/‖ Q_m ‖_op) ).

		for some universal constant c>0.  Note that
		
    ‖ Q_m ‖_op≤‖Σ^1/2‖^2_op‖ P_m‖_op = λ.

		Since Q_m is a real symmetric matrix, it admits a Spectral Decomposition with maximal eigenvalue bounded by λ and rank at most s.  Hence, this gives
		
    ‖ Q_m ‖_HS^2 ≤λ^2 s.

		Directly computing the expectation,
		
    ϵ^ Q_m ϵ = (Q_m) ≤λ s.

		This gives
		
    ( | ‖ P_mξ‖^2 - ‖ P_mξ‖^2 | > Klog(|_u|) ) 
       ≤ 2 exp( -cmin( K^2log^2(|_u|)/λ^2 s , Klog(|_u|)/λ) ).

		This implies that
		
    ∑_m∈_u   (  ‖ P_mξ‖^2 > Klog(|_u|) + λ s ) 
       ≤ 2exp( -cmin( K^2log^2(|_u|)/λ^2 s , Klog(|_u|)/λ) + log(|_u|) ).

		Now, 
		
    ( ∑_m∈_u w_m ‖ P_mξ‖^2 )
    			= ( ∑_m∈_u w_m ‖ P_mξ‖^2 : ) + ( ∑_m∈_u w_m ‖ P_mξ‖^2 : ^)

		Computing each term separately, we have, by the definition of , that
		
    ( ∑_m∈_u w_m ‖ P_mξ‖^2 : ) ≤ Klog(|_u|) + λ s = o(n).

		For the other term, note that
		
    ( ∑_m∈_u w_m ‖ P_mξ‖^2 : ^) 
    			≤( ‖ξ‖^2 : ^)
       ≤( (‖ξ‖^4 ) (^) )^1/2≤( λ^4 (n^2 + 2n )(^))^1/2→ 0

		where the second inequality is due to Cauchy-Schwarz and the limit follows from equation (<ref>), assuming that K>0 was chosen sufficiently large.  Therefore, we have shown that
		
    lim sup_n→∞1/n‖∑_m∈_u w_m( P_mY - μ) ‖^2 ≤ 2t.

		Since t>0 was arbitrary, this proves the first claim.  For the second claim, note that the proof is identical except for defining the set _t as
		
    _t { m∈_u : r_m ≤ t√(n)}.

		This finishes the proof.		
	
	
	

 §.§ Proofs for Section <ref>

	
		Let ξ D^-1/2V^ AZν + D^-1/2V^ Aϵ.  From the proof of Theorem <ref>, it is easy to see that
		
    1/n_a‖ D^-1/2V^ AY - D^-1/2V^ AX‖^2 = 1/n_a‖ξ‖^2 + (1/√(n_a)),
       1/n_b‖ BY - BX‖^2 = 1/n_b‖ Bϵ‖^2 + (1/√(n_b)).

		Therefore, we immediately have that
		
    √(n)(  - σ_ν^2 )
    			   = √(n)( 1/n_a‖ξ‖^2 - 1/n_b‖ Bϵ‖^2 - σ_ν^2 ) + (1).

		We will consider each random variable separately.  Indeed, by the Central Limit Theorem,
		
    √(n_b)( 1/n_b‖ Bϵ‖^2 - σ_ϵ^2 ) (0,2σ_ϵ^4 ),

		or, equivalently,
		
    √(n_b)/√(2σ_ϵ^4)√(n)√(n)( 1/n_b‖ Bϵ‖^2 - σ_ϵ^2 ) (0,1).

		For the other term, we will verify Lindeberg's condition.  Define
		
    μ_k ξ_k^2 = σ_ν^2 + d_k^-1σ_ϵ^2,
       σ_k^2 (ξ^2_k) = 2( σ_ν^2 + d_k^-1σ_ϵ^2 )^2 = 2μ_k^2,
       s_n_a^2 ∑_k=1^n_aσ_k^2.

		Let (ζ_k)_k=1^n_a be a sequence of i.i.d. χ^2_1 random variables.  Then, fixing δ >0, 
		
    lim sup_n_a→∞   1/s_n_a^2∑_k=1^n_a( (ξ_k^2 - μ_k)^2 : | ξ_k - μ_k | > δ s_n_a)
       ≤lim sup_n_a→∞1/s_n_a^2∑_k=1^n_aμ_k^2 ( (ζ_k^2 - 1)^2  : | ζ_k - 1 | > δ s_n_a/μ_k)
       ≤lim sup_n_a→∞1/s_n_a^2∑_k=1^n_aμ_k^2 (( ζ_k^2 - 1 )^4)^1/2(| ζ_k - 1 | > δ s_n_a/μ_k)^1/2
       ≤lim sup_n_a→∞ C ⋁_k=1^n_a(| ζ_k - 1 | > δ s_n_a/μ_k)^1/2→ 0,

		where C is the square root of the fourth central moment of a χ^2_1 random variable.  The second inequality is an application of Cauchy-Schwarz and the limit follows from assumption <ref>.  Therefore, by the Lindeberg Central Limit Theorem,
		
    1/s_n_a( ‖ξ‖^2 - (n_aσ_ν^2 + (D^-1)σ_ϵ^2) ) (0,1).

		Equivalently,
		
    (D^-1)/s_n_a√(n)√(n)( 1/(D^-1)‖ξ‖^2 - (n_a/(D^-1)σ_ν^2 + σ_ϵ^2) ) (0,1).

		By the independence of ξ and Bϵ and assumption <ref>, the desired result now follows by differencing.
	
	
	
	
		Indeed, by the Law of Large Numbers,
		
    σ_ϵ^2.

		Hence, it is clear that σ_b^2.  The fact that d̂ is assumption <ref>.  Finally, for , it follows from Theorem <ref> that
		
    σ_ν^2.

		Now, define ŝ_n_a as
		
    ŝ^2_n_a 2∑_k=1^n_a( /d̂ + d_k^-1)^2.

		It suffices to show that
		
    ŝ^2_n_a/s^2_n_a 1,

		where s^2_n_a is defined in the proof of Theorem <ref>.  Then, we have that
		
    ŝ^2_n_a/s^2_n_a   = ∑_k=1^n_a( /d̂ + d_k^-1)^2/∑_k=1^n_a( σ_ν^2 + d_k^-1σ_ϵ^2 )^2
       = ∑_k=1^n_a( (σ_ν^2 + (1))/( + o(1)) + d_k^-1(σ_ϵ^2 + (1) ) )^2/∑_k=1^n_a( σ_ν^2 + d_k^-1σ_ϵ^2 )^2
       = ∑_k=1^n_a( ( σ_ν^2 + d_k^-1σ_ϵ^2 ) + d_k^-1(1) + (1) )^2/∑_k=1^n_a( σ_ν^2 + d_k^-1σ_ϵ^2 )^2
       = 1 + 2∑_k=1^n_a( σ_ν^2 + d_k^-1σ_ϵ^2 ) (d_k^-1(1) + (1) ) /∑_k=1^n_a( σ_ν^2 + d_k^-1σ_ϵ^2 )^2 + ∑_k=1^n_a(d_k^-1(1) + (1) )^2/∑_k=1^n_a( σ_ν^2 + d_k^-1σ_ϵ^2 )^2.

		We will consider each of the two terms separately. Applying Cauchy-Schwarz to the first term, it follows that
		
    |2∑_k=1^n_a( σ_ν^2 + d_k^-1σ_ϵ^2 ) (d_k^-1 + (1) ) (1)/∑_k=1^n_a( σ_ν^2 + d_k^-1σ_ϵ^2 )^2|
    			≤2√(∑_k=1^n_a(d_k^-1(1) + (1) )^2)/√(∑_k=1^n_a( σ_ν^2 + d_k^-1σ_ϵ^2 )^2) 0.

		For the second term, note that
		
    ∑_k=1^n_a( d_k^-1(1) + (1) )^2/∑_k=1^n_a( σ_ν^2 + d_k^-1σ_ϵ^2 )^2≤(1) 1/σ_ϵ^4 0,

		which finishes the proof.
	
	
	

 §.§ Proofs for Section <ref>

	
	
		Indeed, the consistency of  is by the Law of Large Numbers since
		
    = 1/n_b‖ Bϵ‖^2 + (1).

		Let (ζ_k)_k=1^n_a be a sequence of i.i.d. χ^2_1 random variables.  For , we have that
		
    1/(D^-1)   ‖ D^-1/2V^ AY - D^-1/2 V^ AX‖^2
       =1/(D^-1)‖ D^-1/2V^ AZν + D^-1/2V^ Aϵ‖^2 + (1)
       1/(D^-1)∑_k=1^n_a( σ_ν^2 + d_k^-1σ_ϵ^2 ) ζ_k + (1)
       = n_a/(D^-1)σ_ν^2/n_a∑_k=1^n_aζ_k + σ_ϵ^2/(D^-1)∑_k=1^n_a d_k^-1ζ_k + (1).

		Considering each term separately, we have, by the Law of Large Numbers and assumption <ref>,
		
    n_a/(D^-1)σ_ν^2/n_a∑_k=1^n_aζ_k σ_ν^2.

		For the second term, by a direct variance calculation,
		
    ( σ_ϵ^2/(D^-1)∑_k=1^n_a d_k^-1ζ_k ) = σ_ϵ^4/( D^-1)^2∑_k=1^n_a d_k^-2(ζ_k) = 2σ_ϵ^4 (D^-2)/(D^-1)^2→ 0,

		where the limit uses assumption <ref>.  Hence,
		
    σ_ϵ^2/(D^-1)∑_k=1^n_a d_k^-1ζ_k σ_ϵ^2.

		Combining everything together shows that
		
    σ_ν^2.

		The proof for  is analogous and is omitted, which finishes the proof.
	
	
	
		By independence, the joint distribution of (ν,γ,ϵ) is given by
		
    [ ν; γ; ϵ ]∼_v+r+n( 
    			[ 0_v; 0_r; 0_n ] , 
    			[ σ_ν^2 I_v    0_v× r    0_v× n;    0_r× v σ_γ^2 I_r    0_r× n;    0_n× v    0_n× r σ_ϵ^2 I_n ]).

		Therefore, the joint distribution of Y and η is given by
		
    [ Y; η ]∼_2n( 
    			[ μ; μ ],
    			[ σ_ν^2 ZZ^ + σ_γ^2 WW^ + σ_ϵ^2 I_n                         σ_ν^2 ZZ^;                         σ_ν^2 ZZ^                         σ_ν^2 ZZ^ ]).

		By standard results on the conditional mean of η|Y, it follows that
		
    ( η | Y ) = μ + σ_ν^2 ZZ^( σ_ν^2 ZZ^ + σ_γ^2 WW^ + σ_ϵ^2 I_n)^-1(Y - μ).

		which finishes the proof.
	
	
	
		From Proposition <ref>, it follows that
		
    /d̂ = σ_ν^2 + δ_ν^2          /λ̂ = σ_γ^2 + δ_γ^2        = σ_ϵ^2 + δ_ϵ^2,

		where δ_ν^2, δ_γ^2, and δ_ϵ^2 are all (1).  We will write δ_∗^2 max{δ_ν^2, δ_γ^2, δ_ϵ^2}.  Define
		
    Σσ_ν^2 ZZ^ + σ_γ^2 WW^ + σ_ϵ^2 I_n,
       Δδ_ν^2 ZZ^ + δ_γ^2 WW^ + δ_ϵ^2 I_n.

		Note that Σ is positive definite and Δ is invertible almost surely.  Then, by the Matrix Inversion Lemma,
		
    ( Σ + Δ)^-1 = Σ^-1 - Σ^-1( Δ^-1 + Σ^-1)^-1Σ^-1.

		Some algebra will yield
		
    =   μ̂_EW + /d̂ ZZ^( /d̂ ZZ^ + /λ̂WW^ +  I_n )^-1( Y - μ̂_EW)
    
    			=   μ̂_EW + ( σ_ν^2 + δ_ν^2 ) ZZ^( Σ + Δ)^-1( Y - μ̂_EW)
    
    			=   (μ̂_EW - μ) + ( σ_ν^2 + δ_ν^2 ) ZZ^( Σ + Δ)^-1( μ - μ̂_EW) 
    			
       + δ_ν^2 ZZ^( Σ + Δ)^-1( Zν + Wγ + ϵ)
    			
       - σ_ν^2 ZZ^Σ^-1( Δ^-1 + Σ^-1)^-1Σ^-1( Zν + Wγ + ϵ)
    			
       +μ + σ_ν^2 ZZ^Σ^-1( Zν + Wγ + ϵ),

		where we have added and subtracted μ and applied the Matrix Inversion Lemma.  From Lemma <ref>, it follows that
		
    = μ + σ_ν^2 ZZ^Σ^-1( Zν + Wγ + ϵ).

		Define ξ -.  Therefore,
		
    1/n( ‖ - η‖^2 - ‖ - η‖^2 )
    			= 1/n( ‖ξ‖^2 + 2ξ^(  - η) ).

		We will prove each of the two terms on the right hand side are (1).  Before doing so, we will prove a few useful facts that will facilitate the remainder of the proof.
		
			
  * ‖ ZZ^‖^2 = 𝒪(1).
			
  * ‖ Zν + Wγ + ϵ‖^2 = (n) and ‖ Zν‖^2 = (n).
			
  * ‖Δ‖^2 = (1) and ‖( Σ + Δ)^-1‖^2 = (1).
		
		(I) is immediate since Σ is the sum of three positive semi-definite matrices, one of which being σ_ν^2 ZZ^, and ‖Σ‖^2 = 𝒪(1) by assumption <ref>.
		The first part of (II) follows from the fact that 
		
    ‖ Zν + Wγ + ϵ‖^2 ≼λ_max( Σ) χ^2_n = (n),

		where ≼ denotes stochastic ordering and λ_max( Σ) = 𝒪(1) by assumption <ref>.  The second part of (II) is identical except we apply (I) to bound the maximal eigenvalue.  Finally, for (III), the first part follows from δ_∗^2  0 and ‖Σ‖ = 𝒪(1).  For the second part, note that
		
    ‖( Σ + Δ)^-1‖^2 ≤‖Σ^-1‖^2 ‖( I_n + Σ^-1Δ)^-1‖^2 = 𝒪(1) ‖( I_n + Σ^-1Δ)^-1‖^2.

		Then, ‖Σ^-1Δ‖ 0.  Therefore,
		
    ‖( I_n + Σ^-1Δ)^-1‖^2
    			= ‖∑_k=0^∞(- Σ^-1Δ)^k‖^2
    			≤∑_k=0^∞‖(- Σ^-1Δ)^k‖^2 < ∞.

		This proves all three claims.
		We can now show that ‖ξ‖^2 = (n), which we will show in parts.  To that end, note that
		
    ‖(I_n - ( σ_ν^2 + δ_ν^2 ) ZZ^( Σ + Δ)^-1) ( μ̂_EW - μ) ‖^2
       ≤‖ I_n - ( σ_ν^2 + δ_ν^2 ) ZZ^( Σ + Δ)^-1‖^2 ‖μ̂_EW - μ‖^2
       = (n).

		In the last line, we have used (I) and (III) from above and Lemma <ref>.  Similarly,
		
    ‖   δ_ν^2 ZZ^( Σ + Δ)^-1( Zν + Wγ + ϵ) ‖^2
       ≤δ_ν^4 ‖ ZZ^‖^2 ‖( Σ + Δ)^-1‖^2 ‖ Zν + Wγ + ϵ‖^2 
       = (1).

		The last line follows from the fact that δ_ν^2  0 and (I), (II), and (III) from above.  For the last term, we will apply the Matrix Inversion Lemma again to obtain
		
    ( Δ^-1 + Σ^-1)^-1 = Δ - Δ( Δ + Σ)^-1Δ.

		Hence,
		
    ‖σ_ν^2 ZZ^Σ^-1( Δ^-1 + Σ^-1)^-1Σ^-1( Zν + Wγ + ϵ) ‖^2
       ≤σ_ν^4 ‖ ZZ^‖^2 ‖Σ^-1‖^4 ‖Δ‖^2 ‖ I_n - ( Δ + Σ)^-1Δ‖^2 ‖ Zν + Wγ + ϵ‖^2
       = (n),

		where (I), (II), and (III) have all been used from above.  Combining these three results with the triangle inequality, this proves that
		
    ‖ξ‖^2 = (n).

		For the other quantity, we have that
		
    ‖ - η‖^2    = ‖σ_ν^2 ZZ^Σ^-1( Zν + Wγ + ϵ) - Zν‖^2
       ≤ 2 ‖σ_ν^2 ZZ^Σ^-1( Zν + Wγ + ϵ) ‖^2 + 2‖ Zν‖^2
       ≤ 2 ‖σ_ν^2 ZZ^Σ^-1( Zν + Wγ + ϵ) ‖^2 + 2‖ Zν‖^2
       ≤ 2 σ_ν^4 ‖ ZZ^‖^2 ‖Σ^-1‖^2 ‖  Zν + Wγ + ϵ‖^2 + 2‖ Zν‖^2.

		Applying all three facts from above demonstrates that
		
    ‖ - η‖^2 = (n).

		Using the Cauchy-Schwarz inequality yields
		
    2|ξ^(  - η)| ≤ 2‖ξ‖‖ - η‖ = (n),

		which finishes the proof.
	
	
	
	
newapa
acm	

	
