Fast and reliable confidence intervals for a variance component or proportion
Yiqiao Zhang Karl Oskar Ekvall Aaron J. Molstad

Outage Probability Analysis of  Wireless Paths with Faulty Reconfigurable Intelligent Surfaces
    Mounir Bensalem and Admela Jukan
Technische Universität Braunschweig, Germany;
{mounir.bensalem,  a.jukan}@tu-bs.de







    Version of April 23, 2024
==============================================================================================================================






    We show that confidence intervals for a variance component or proportion,
    with asymptotically correct uniform coverage probability, can be obtained by
    inverting certain test-statistics based on the score for the restricted
    likelihood. The results apply in settings where the variance or proportion
    is near or at the boundary of the parameter set. Simulations indicate the
    proposed test-statistics are approximately pivotal and lead to confidence
    intervals with near-nominal coverage even in small samples. We illustrate
    our methods' application in spatially-resolved transcriptomics where we
    compute approximately 15,000 confidence intervals, used for gene ranking, in
    less than 4 minutes. In the settings we consider, the proposed method is
    between two and 28,000 times faster than popular alternatives, depending on
    how many confidence intervals are computed.






§ INTRODUCTION

Variance component models are widely used in applied
statistics. Often, researchers are interested in quantifying how much of the
variability in the data is due to a specific source, possibly after controlling
for, or conditioning on, other sources of variability. In modern genetics and
genomics, for example, it is common to assume a vector Y ∈ℝ^n of
responses for n sample units satisfies

    Y ∼𝒩_n(Xβ, σ_g^2 K + σ_e^2 I_n),

with X ∈ℝ^n× p a design matrix, K ∈n× n a
known positive semi-definite matrix, and parameter θ = (β, σ_g^2,
σ_e^2) ∈p× [0, ∞) × (0, ∞)
<cit.>.  In statistical genetics, when
elements of Y are phenotypes for n individuals, K can be a genetic
relatedness matrix quantifying similarities between individuals' common SNP
genotypes <cit.>. Then, the quantity of inferential interest
is the heritability <cit.>

    h^2 = σ_g^2/σ_g^2 + σ_e^2∈ [0, 1).

In spatially resolved transcriptomics <cit.>, when elements of
Y are a gene's expression at n distinct two- or three-dimensional spatial
positions in a tissue sample and K is a matrix encoding distances between
positions, h^2 is known as the fraction of spatial variance
<cit.>.

Though (<ref>) may seem a simple model, inference on h^2 can be
complicated, especially when it is nearly, or equal to, zero or one, as those
are boundary points of the parameter set. It is well-known that common
test-statistics, such as Wald and likelihood ratio statistics, have non-standard
distributions on boundary points. The distributions can often be obtained by
simulation <cit.> or approximated, either using asymptotic theory
<cit.>, or by other means
<cit.>. Thus, it is often possible to construct tests with
(asymptotically) correct pointwise size, and the likelihood ratio test has
recently been suggested for use in related settings
<cit.>. However, creating a reliable confidence region by
inverting such a test-statistic is more complicated, for at least two reasons. 

First, one in general has to compute the statistic for every relevant h^2.
Thus, methods which are designed specifically for testing h^2 = 0 are
inapplicable, and simulation-based methods can be prohibitively time-consuming.
Secondly, even when inverting a test with asymptotically correct pointwise size
is computationally feasible, the resulting confidence region can have poor
coverage properties in finite samples. Specifically, the actual coverage
probability often depends substantially on how close the true parameter is to
the boundary. To illustrate the issues, Fig. <ref> shows
coverage probabilities for different values of h^2 for Wald, likelihood ratio,
and proposed confidence intervals defined in Section  <ref>. For
interior points of the parameter set, the confidence intervals are obtained by
inverting test-statistics using quantiles of the chi-square distribution with
one degree of freedom, as suggested by classical theory. Notably, the coverage
probabilities for Wald and likelihood ratio intervals are substantially
different from the nominal level when h^2 is near zero or one. The reason is
that, when the true parameter is close to the boundary, the distributions of the
test-statistics are close to what they are on the boundary. Thus, even when the
parameter is in the interior, classical asymptotic theory for interior points
provides a poor approximation. Removing the boundary points from the parameter
set or using a larger sample size is not sufficient to address the issue, as
Fig. <ref> illustrates; loosely speaking, for every sample
size, there are interior points close enough to the boundary to cause problems.
More formally, test-statistics with different asymptotic distributions on
boundary and interior points lead to asymptotically incorrect uniform
coverage probability <cit.>.



For reasons to become apparent, we focus on test-statistics based on the score
and information of the restricted likelihood. By doing so we are able to
construct a (1- α) confidence region ℛ_n(α) that
satisfies, for any compact C ⊆ [0, 1] and as n→∞,

    sup_θ∈ C|_θ{h^2 ∈ℛ_n(α)} - 
            (1 - α)|→ 0.

That is, the proposed confidence region, which is typically an interval, has
asymptotically correct uniform coverage probability on compact sets. This
property leads to near-nominal coverage of every practically relevant h^2,
including h^2 near or equal to zero, as is illustrated in
Fig. <ref>. We emphasize that when focusing on uniform
inference in the presence of boundary points, conclusions about score
test-statistics in general do not carry over to Wald and likelihood ratio
test-statistics. Indeed, our results show, in a sense to be made precise, that
certain score-statistics are asymptotically pivotal, but Wald and likelihood
ratio test-statistics are not. The score statistics can therefore be inverted
using the same quantiles for every point in the parameter set, leading to a
confidence interval that is easy to implement and orders of magnitude faster to
compute than competing ones. In a typical example we consider, our confidence
intervals take less than a second to compute, whereas the competitors often take
500 seconds, and sometimes much longer.

Previously, the score standardized by expected information has been shown to
give confidence regions with asymptotically correct uniform coverage probability
in settings with singular Fisher information
<cit.>. We build on ideas from that literature
here, but the setting, and hence the arguments required for the main results,
are quite different. In particular, here, the information matrix is positive
definite under weak conditions (Theorem <ref>). There are also two
nuisance parameters, β and σ^2 = σ_g^2 + σ_e^2. While
dealing with β is relatively straightforward since β and (h^2,
σ^2) are orthogonal in the Fisher information sense, our results are, to
the best of our knowledge, the first of their kind with non-orthogonal nuisance
parameters. In contrast to classical asymptotic theory for inference in mixed
models, we also let the number of predictors grow with the sample size. This is
an important development as the restricted likelihood is typically used
precisely when the number of predictors is smaller than the number of
observations, but large enough to induce substantial bias in the inference on
h^2.



§ MODEL
 

Let Σ = _θ(Y) = σ_g^2K + σ_e^2I_n. A log-likelihood
corresponding to (<ref>) is

    l(β, σ_g^2, σ_e^2) = l(β, σ_g^2, σ_e^2; Y, X, K) =
    -1/2{log|Σ|+(Y-Xβ)^Σ^-1(Y-β)},

the domain of which is the parameter set p× [0, ∞) × (0,
∞). It is important that (<ref>) allows σ^2_g =
0, which is often of practical interest. When K is positive definite, points
with σ^2_e = 0 and σ_g^2 > 0 can be added to the parameter set, but
with semi-definite K, (<ref>) requires
σ_e^2 > 0. Assuming X has full column rank, the restricted
log-likelihood, often used for inference on the covariance parameters, is

    l_R(σ_g^2, σ_e^2) = -1/2{log|Σ| + 
        log| X^Σ^-1X| + (Y - Xβ̃)^Σ^-1
        (Y - Xβ̃)},

where β̃ = β̃(σ_g^2, σ_e^2) = (X^Σ^-1X)^-1X^Σ^-1 Y is the partial maximizer of l(β,
σ^2_g, σ^2_e) in β, for fixed σ_g^2 and σ^2_e.
The domain of l_R is [0, ∞) × (0, ∞). The assumption that X
have full column rank is for notational convenience; one could define l_R
similarly as long as the column rank is less than n. Specifically, if the
column rank of X is p'< n and V ∈n× (n - p') is a
semi-orthogonal matrix satisfying V^ X = 0, then the restricted likelihood is
the likelihood of V^ Y <cit.>.

By spectral decomposition, K = OΛ O^ for some orthogonal O ∈n× n and Λ = (λ_1, …, λ_n); here and
elsewhere, eigenvalues are in decreasing order unless otherwise stated. It
follows that Σ = O (σ^2_gΛ + σ^2_e I_n)O^ and, hence,

    l(β, σ^2_g, σ^2_e) = -1/2∑_i = 1^n {log(λ_iσ^2_g + σ^2_e) + 
            (O^ Y - O^ X β)_i^2/λ_i σ^2_g  + σ^2_e},

where (·)_i means the ith element; and 

    l_R(σ^2_g, σ^2_e) = l(β̃, σ^2_g, σ^2_e) - log| (O^ X)^ (σ^2_g Λ + σ^2_e I_n)^-1 O^ X| / 2.
 
Thus, when developing theory, we may upon replacing Y and X by O^ Y and
O^ X, respectively, assume diagonal covariance matrices K = Λ
(equivalently, O = I_n) and Σ = σ^2_g Λ + σ^2_e I_n
without loss of generality. When doing so, the observations are ordered so that
the ith response Y_i corresponds to the ith eigenvalue λ_i.

At the expense of more notation, the assumption on Σ could be relaxed to
Σ = σ^2_G K + σ^2_e W for some known positive definite W ∈n× n. Indeed, in that case Y, X, and K can be replaced by,
respectively, W^-1/2Y, W^-1/2X, and W^-1/2K W^-1/2 before
repeating the above argument to get a diagonal Σ.

In several motivating applications, the parameter of interest is h^2 =
σ^2_g / (σ^2_g + σ^2_e). Therefore, we consider a common
reparameterization in terms of h^2 and σ^2 = σ^2_g + σ^2_e.
Then, Σ = σ^2 O{h^2Λ + (1 - h^2)I_n}O^. Overloading the
letter l and assuming O = I_n without loss of generality, the log-likelihood
becomes

    l(β, h^2, σ^2) = -1/2∑_i = 1^n {log(σ^2) + log(λ_i h^2 + 1 - h^2) + 
        (Y_i - X_i^β)^2/σ^2(λ_i h^2 + 1 - h^2)},

where X_i^ is the ith row of X. The domain of l is the parameter set
p× [0, 1) × (0, ∞). The corresponding restricted
log-likelihood, with domain [0, 1) × (0, ∞), is

    l_R(h^2, σ^2) = l(β̃, h^2, σ^2) +
        p/2log (σ^2) - 1/2log| X^{h^2Λ + (1 - h^2)I_n}^-1X|.




 §.§ Test-statistics and confidence regions


There are many potential test-statistics for inference based on l(β, h^2,
σ^2) and l_R(h^2, σ^2). However, most common ones, like Wald and
likelihood ratio statistics, do not give confidence regions satisfying
(<ref>) in boundary settings in general. Therefore, we focus on
score statistics standardized to have mean zero and diagonal covariance matrix.
Such test-statistics have previously been shown to lead to
(<ref>) in settings where, unlike here, the Fisher information
is singular <cit.>. Because our goal is to
construct confidence regions for the covariance parameters, effectively treating
β as a nuisance parameter, we primarily work with the restricted
likelihood.

Let V_1 = V_1(h^2, σ^2) = ∂Σ /
∂ h^2 = σ^2O(Λ-I_n)O^ and V_2 = V_2(h^2, σ^2) =
∂Σ / ∂σ^2 = O(h^2Λ + (1 - h^2)I_n)O^. Similarly,
let U_1(h^2, σ^2) = ∂ l_R(h^2, σ^2) / ∂ h^2 and
U_2(h^2, σ^2_e) = ∂ l_R(h^2, σ^2) / ∂σ^2, so
the restricted score is U(h^2, σ^2) = [U_1(h^2, σ^2), U_2(h^2,
σ^2)]^. Routine calculations show, for j∈{1, 2},

    U_j(h^2, σ^2)    =
        1/2(Y - Xβ̃)^Σ^-1V_j Σ^-1 (Y - Xβ̃) 
        -1/2(Σ^-1V_j)
          + 1/2(Σ^-1X(X^Σ^-1X)^-1X^Σ^-1V_j).

Observe U_j is defined on [0, 1)×(0, ∞) and has a continuous
extension to points where h^2 = 1, σ^2 > 0 if K is positive definite. The expression
simplifies substantially for j = 2 since V_2 = Σ / σ^2.
Specifically, assuming O = I_n,

    U_2(h^2, σ^2) = -n - p/2 σ^2 +
        1/2 σ^4∑_i = 1^n (Y_i - X_i^β̃)^2/h^2 λ_i + 1 - h^2.

Let ℐ(h^2, σ^2) denote the restricted information matrix, that
is, the covariance matrix of U(h^2, σ^2) when h^2 and σ^2 are
the true parameters. We next give an expression for ℐ(h^2, σ^2)
for easy reference. The proof, along with those of other formally stated
results, is in Appendix <ref>. 

Let P = Σ^-1/2X(X^Σ^-1X)^-1X^Σ^-1/2, Q = I_n - P,
D = {h^2 Λ + (1 - h^2)I_n }^-1(Λ - I_n), and H = ODO^.

 
    If (<ref>) holds for some X with full column rank, then the
    restricted information matrix has elements
    
    ℐ_11(h^2, σ^2) = 1/2(Q H QH);   ℐ_12(h^2, σ^2) = 1/2σ^2(Q H);  ℐ_22(h^2, σ^2)     = n - p/2σ^4.



For inference on (h^2, σ^2) jointly, we consider the (quadratic)
restricted score standardized by the restricted Fisher information:

    T^R_n(h^2, σ^2) = U(h^2, σ^2)^ℐ^-1(h^2, σ^2)U(h^2, σ^2),

For inference on h^2 only, we use

    T^R_n(h^2) = U_1(h^2, σ̃^2)^2
        ℐ^11(h^2, σ̃^2),

where ℐ^11(h^2, σ^2) is the leading element of ℐ^-1(h^2, σ^2) and σ̃^2 = σ̃^2
(h^2) is the partial maximizer of l_R(h^2, σ^2) in σ^2, for a
fixed h^2. By (<ref>), assuming O = I_n,

    σ̃^2 = 1/n - p∑_i = 1^n 
        (Y_i - X_i^β̃)^2/h^2 λ_i + 1 - h^2  = 1/n - p(Y - Xβ̃)^{h^2 Λ + (1 - h^2)I_n}^-1 (Y - Xβ̃).


When X has column rank less than n, σ̃^2 is greater than
zero with probability one, and hence a valid partial maximizer. Note
β̃ does not depend on σ^2, so the last display indeed
gives an explicit expression for σ̃^2. For inference on h^2,
we sometimes also consider S_n^R(h) = U_1(h^2, σ̃^2)
ℐ^11(h^2, σ̃^2)^1/2, that is, the signed square root
of (<ref>).


The proposed test-statistics require a positive definite ℐ(h^2,
σ^2), which the following result establishes. Notably,
ℐ^11(h^2, σ^2) is positive if ℐ(h^2, σ^2) is
positive definite. To understand the next result, notice that
Q has rank n-p, so QHQ has at most n-p non-zero eigenvalues, and recall
V ∈n × (n - p) is a semi-orthogonal matrix such that V^ X =
0.

 
   The information matrix ℐ(h^2, σ^2) is singular if and only
   if the following equivalent conditions hold: (i) the n-p possibly non-zero
   eigenvalues of QHQ are identical; (ii) V^ K V = c I_n-p for some c≥ 0.


The proof of Theorem <ref>
is relatively straightforward using the expressions in Proposition <ref>. 
The key observation is that the determinant of ℐ(h^2, σ^2)
is

    1/2σ^4{ (n-p)(QHQ QHQ) - (QHQ)^2},

which is strictly positive by Jensen's inequality applied to the n-p possibly
non-zero eigenvalues of QHQ unless they are all identical, in which case the
determinant is zero.

Theorem <ref> is particularly intuitive if p = 0, so that Q = V=
I_n. Then the eigenvalues of QHQ are the diagonal elements of D, which are
identical if and only if the diagonal elements of Λ are; that is, if and
only if Λ, and hence K, is proportional to the identity. The
parameters are unidentifiable for such K, so we will generally assume K is
not proportional to the identity, leading to a positive definite information
matrix.

By inverting the proposed test-statistics using quantiles from chi-squared distributions, we get the confidence regions

    ℛ_n(α)    = {h^2 ∈ [0, 1] : T_n^R(h^2) ≤ q_1 - α, 1};
    ℛ^2_n(α)    = {(h^2, σ^2) ∈ [0, 1] × (0, ∞) : T_n^R(h^2, σ^2) ≤ q_1 - α, 2}.

 When h^2 = 1 is not in the parameter set, we take T^R_n(1) = T^R_n(1, σ^2) = ∞ so the confidence regions are well defined.



§ ASYMPTOTIC DISTRIBUTIONS
 

We assume K_n = Λ_n = (λ_n1, …, λ_nn) in this
section without loss of generality, and hence Σ_n = σ^2_n(I_n +
Λ_n) and H_n = D_n are diagonal. Here, a subscript n has been added
to σ^2 and other quantities since we will consider a sequence
{(β_n, h^2_n, σ^2_n)} of parameters indexed by the sample size,
and a corresponding sequence of models. Similarly, the number of predictors p =
p_n can change with n. By considering a sequence of (true) parameters, we can
obtain the asymptotic distribution of T^R_n(h^2_n, σ^2_n) and
T_n^R(h_n^2) under that sequence. This is useful because
(<ref>) holds for ℛ_n(α) in (<ref>) if and only if the asymptotic distribution of
T_n^R(h_n^2) is the same under any convergent sequence of parameters
<cit.>, and similarly for ℛ^2_n(α). This result essentially follows from the
equivalence of continuous convergence and uniform convergence on compact sets,
applied to the sequence of functions defined by θ↦_θ{T_n^R(h^2) ≤ q_1-α, 1}. More generally, examining the distribution under a sequence of
parameters tending to a point, can lead to a better understanding of the
behavior of the test-statistic near that point. Thus, for our purposes,
sequences tending to boundary points are of particular interest. The next result
is our first along these lines.

To state the result, let γ_i(·) denote the ith eigenvalue in
decreasing order.

 
    Assume that, for every n ∈{1, 2, …}, Y satisfies
    (<ref>) with an X with full column rank, K_n = Λ_n,
    and parameters (β_n, h^2_n, σ^2_n). If, as n→∞, (i) p/n → 0, (ii) ∑_i = 1^p γ_i(D_n^2) / (D_n^2) → 0, (iii) lim sup_n→∞{n^-1(D_n)^2 / (D_n^2)} < 1; then, in distribution,







    
    T_n^R(h_n^2, σ^2_n)→χ_2^2.



The proof of Lemma <ref> establishes (<ref>) by
showing asymptotic bivariate normality of ℐ_n(h_n^2,
σ_n^2)^-1/2U_n(h_n^2, σ^2_n). The key observation is that the
elements of U_n(h_n^2, σ^2_n) are quadratic forms in multivariate normal
vectors, centered to have mean zero. Therefore, their distributions are those of
linear combinations of independent chi-square random variables. The conditions
of the theorem ensure Lyapunov's central limit theorem applies to those
linear combinations.

Conditions (ii) and (iii) concern Λ_n and h_n^2, but neither of
β_n and σ^2_n appears in the conditions. Thus,
(<ref>) can hold for sequences one may not expect, such as
those where σ^2_n→ 0. This suggests the behavior of our
test-statistic is insensitive to the true σ^2, and may perform well
even if σ^2 is near the boundary, that is, near zero.

Condition (ii) says the sum of the p greatest eigenvalues of D_n^2 should
be negligible compared to the sum of all its eigenvalues. An eigenvalue, or
diagonal element, of D_n^2 will be relatively large if either (a) the
corresponding element of Λ_n is relatively close to zero and h^2
is close to one, or (b) the corresponding element of Λ_n is
relatively large. Thus, condition (ii) suggests the test-statistic may not
behave well if a few elements of O^ Y have much
larger variance than all other ones, or if some elements have very small
variance and most if it is attributable to K.

Condition (iii) essentially says that, when comparing the squared average of
the eigenvalues of D_n to the average of the squared eigenvalues, Jensen's
inequality is asymptotically strict; there has to be non-negligible
variation in the diagonal elements of D_n. This is intuitive as, if all
the diagonal elements of D_n are the same, then Λ_n is
proportional to the identity and h_n^2 unidentifiable, as noted in the
remark following Theorem <ref>. 

Conditions (ii) and (iii) could be replaced by slightly weaker ones by adding
assumptions on how X relates to K, but that would substantially increase
notational burden.


Our next result shows that the proposed test-statistic for h^2 has the desired
asymptotic distribution under the conditions of Lemma <ref>.
Substantial work is needed to establish the result since σ^2 is a
nuisance parameter for the purposes of inference on h^2 only. In particular,
the test-statistic for h^2 cannot depend on or be evaluated at the true
parameter, or null hypothesis, which was possible in Lemma <ref>.

 
Under the conditions of Lemma <ref>, T_n^R(h_n^2) →χ^2_1 in distribution as n→∞.


Lemma <ref> is proven by showing that T_n(h_n^2) has the same
asymptotic distribution as {U_n1 - ℐ_n12 U_n2 /
ℐ_n22}^2 ℐ_n^(11), with all quantities evaluated at
the true (h_n^2, σ^2_n) and ℐ_n^(11) denoting the leading
element of ℐ_n^-1. This result is consistent with what is often
obtained for nuisance parameters in classical settings, with fixed interior
parameters <cit.>. This is remarkable because those results in
general do not apply under sequences of parameters or with boundary points.

Our next result gives intuitive conditions on model parameters and K which
ensure the conditions in Lemmas <ref> and <ref>
hold. The conditions are not necessary in general, but they are more transparent, and
easier to assess in practice, than those in the lemmas.


    Assume that, for every n ∈{1, 2, …}, Y satisfies
    (<ref>) with an X with full column rank and parameters
    (β_n, h^2_n, σ^2_n). If (i) lim_n→∞ p_n/n = 0; (ii^*) lim sup_n→∞λ_n1 < ∞; (iii^*) lim inf_n→∞λ_nn > 0 or
        lim sup_n→∞h^2_n < 1; (iv^*) there exists an ϵ > 0 such that for any c ∈
        [0,∞), with k^c_n denoting the number of eigenvalues of
        Λ_n in [c-ϵ, c + ϵ], it holds that
        lim sup_n→∞k_n^c / n < 1,
 
	
	
	
	
	
	
	
	
	
	
 
	
	
	
	
 
 
 
	










    then, in distribution as n→∞,
    
    T_n^R(h_n^2, σ^2_n) →χ_2^2;   T_n^R(h_n^2)→χ_1^2.



Since Λ_n is known, condition (ii^*) can be made to hold, for example by
replacing K by K/λ_n1. Then, for the conclusions to hold, conditions (iii^*) and (iv^*) would need to hold for the new eigenvalues λ_ni/λ_n1,
i ∈{1, …, n}.  Condition (iii^*) illustrates that if K is full rank, then
reliable inference on h^2 near one is possible; this is formalized further in
Corollary <ref>. Condition (iv^*) prevents the eigenvalues of K from
concentrating in some neighborhood asymptotically. That is, it prevents K from
being essentially proportional to the identity asymptotically.

 
    Under condition (i), (ii^*), and (iv^*) of Theorem <ref>, ℛ_n(α) and ℛ^2_n(α) in (<ref>) have asymptotically correct uniform coverage probability on
    compact subsets of [0, 1) and [0, 1)× (0, ∞), respectively
    (c.f. (<ref>)). The conclusion can be strengthened to
    compact subsets of [0, 1]× (0, ∞) and [0, 1], respectively, if
    the first part of condition (iii^*) of Theorem <ref> also
    holds.






§ COMPUTATIONAL CONSIDERATIONS


In this section, we consider the computational aspects of our test-statistics.
Like popular competing methods, we require a one-time preprocessing step where
the eigendecomposition of K = OΛ O^ is used to construct Y_O = O^ Y
and X_O = O^ X. In practice it is common to use the same K for many
different choices of response vector Y, with separate inference on h^2 for
each response vector. Then, the pre-processing step is only performed once.

We show that, after the pre-processing, the number of floating point operations
required to evaluate the proposed test-statistics grows linearly in n, making
our methods feasible for large-scale applications, such as those described in
Section <ref>. More specifically, we show that, if p is
potentially growing with n, the number of operations required is O(np^2 +
p^3).

Recall V_2 = h^2Λ + (1 - h^2)I_n and define P̃ =
V_2^-1/2X_O(X_O^ V_2^-1X_O)^-1X_O^ V_2^-1/2, Q̃ = I_n -
P̃, β̃ = (X_O^⊤ V_2^-1X)^-1X_O^⊤ V_2^-1Y_O, and
R̃ =Y_O - X_Oβ̃. Let R̃_i denote the ith
element of R̃ and Q_ii the ith diagonal element of Q. With
these definitions the score function is

    𝑈(h^2, σ^2) = 1/2[   ∑_i=1^n D_ii(σ^-2R̃_i^2 /h^2 λ_i + 1 - h^2  -  Q̃_ii); -n - p/σ^2 +
        1/σ^4∑_i= 1^n R̃_i^2/h^2 λ_i + 1 - h^2 ].

Both terms in (<ref>) involve the summation of n terms.
Construction of R̃ requires O(np) operations; and computation of all
Q̃_ii  requires only O(np^2) operations.  We do not
compute or store the entire matrix Q̃. 

To compute the information efficiently, note Q = OQ̃O^⊤, so by
Proposition <ref>, 

    ℐ(h^2, σ^2) = 1/2σ^4([ σ^4  tr(Q̃D Q̃D)     σ^2  tr(Q̃D);     σ^2  tr(Q̃D)            n - p ]).
 Because tr(Q̃D) = ∑_i=1^n Q̃_ii D_ii, ℐ_12 can be
evaluated in O(n) operations using the Q̃_ii computed before.
Notice 2ℐ_11 =  tr(Q̃D Q̃D) =
 tr(D^2) - 2 tr(P̃D^2) +  tr(P̃DP̃D).
Because D is a diagonal matrix, evaluating tr(D^2) requires O(n)
operations. Similarly, tr(P̃D^2) = ∑_i=1^n (1 -
Q̃_ii)D_ii^2, so evaluating this term also requires only O(n)
operations.  Finally, by the cyclic property of the trace operator, tr(P̃DP̃D) =  tr(ÃB̃ÃB̃)
where Ã = (X_O^ V_2^-1X_O)^-1∈ℝ^p × p and
B̃ = X_O^ V_2^-1/2D  V_2^-1/2X_O ∈ℝ^p × p. Thus,
evaluating tr(ÃB̃ÃB̃) requires O(np^2
+ p^3) operations, with the p^3 coming from a Cholesky (or other)
decomposition of Ã^-1 to deal with the inverse, and the np^2 is
from matrix multiplications.

We have created an R package  for computing our test-statistics
and confidence intervals, available for download at
. 

Our confidence interval for h^2 can be obtained by evaluating T_n^R(h^2)
from at equally spaced points in [0,1) and comparing the values to quantiles
from the chi-square distribution with one degree of freedom. In practice,
however, it is typically possible, and faster, to do a search for the points
where T_n^R(h^2) crosses the desired quantile. More specifically, when
T_n^R(h^2) is strictly quasiconvex in h^2, then ℛ_n(α) is
an interval, so we need only search for points a such that T_n^R(a) = q_1,
1-α. While we have not been able to prove T_n^R(h^2) is quasiconvex in
general due to technical difficulties, we have not encountered an example where
it is not. Typically, T_n^R(h^2) looks roughly like one of the three cases in
Fig. <ref>. 
Thus, in our software we first find a point
a' such that T_n^R(a') < q_1, 1-α. More specifically, we start a
ternary search for the minimum of T_n^R and terminate the first time a
negative value is encountered. Then, under quasiconvexity, the intervals
[0,T_n^R(a')] and [T_n^R(a'),1) contain lower and upper bounds of the
interval ℛ_n(α), respectively. We then use a bisection search
to find b ∈ [0,T_n^R(a')] such that T_n^R(b) = q_1, 1-α, and
similarly for [T_n^R(a'),1). This procedure often requires less computing time
than an exhaustive grid search.

We have also found that, in practice, S_n^R(h^2) = {ℐ_n^11(h^2,
σ̃^2)}^1/2U_1(h^2, σ̃^2) is typically monotonely
decreasing in h^2, so a search can also be based on that statistic and
comparisons to standard normal quantiles; see Section <ref> for
an example with one-sided confidence intervals.





























































































































§ NUMERICAL STUDIES



 §.§ Finite-sample interval width and coverage probabilities

To supplement our asymptotic theory, we study how n, h^2, and the
eigenvalues of K affect the coverage probability and width of our confidence
intervals. We consider K_ij = ρ^|i-j| for ρ∈ (0, 1) and i, j
∈{1, …, n}. When ρ is close to one, K is approximately
1_n1_n^, where 1_n is an n-vector of ones. Thus, approximately,
λ_n1 = n and λ_ni = 0, i ≥ 2, which is inconsistent
with conditions (ii^*) and (iv^*) of Theorem <ref>. Conversely,
when ρ is close to zero, K is approximately an identity matrix, which is
inconsistent with condition (iv^*) of Theorem <ref>. Thus, we
expect our test-statistic to perform best when ρ is not too close to either
zero or one.

In each setting considered, we generate X ∈ℝ^n × 5 with
independent and identically distributed standard normal entries, and then
generate Y ∼𝒩_n(Xβ, h^2 K + ( 1- h^2) I_n) (i.e.,  σ^2
= 1) with β = 0. For every combination of the triplet (n, h^2, ρ),
we generate 10^4 independent realizations of Y and compute the corresponding
95% confidence intervals for h^2. 

The bottom panel of Fig. <ref> displays coverage probabilities
with ρ and n varying for four different values of h^2. For ρ∈{0.1. 0.5. 0.95}, our interval appears to have the coverage probability 0.95
or higher. Only when ρ = 0.999 is the evidence of under-coverage for some
combinations of n and h^2. 

The top panel of Fig. <ref> displays interval widths for different
values of ρ and h^2. When ρ = 0.1, K ≈ I_n, which means
h^2 is approximately unidentifable. Consequently, interval widths are too
large to be practically useful, despite having the nominal coverage level. With
ρ = 0.5 or 0.95, the widths are reasonably small and decrease as n
increases. When ρ = 0.999, the eigenvalues of K are approximately (n, 0,
…, 0), which appear to lead to longer confidence intervals. Nevertheless,
as n increases, interval widths decrease as expected. As n approaches 2000,
the widths are small enough to be useful in practice. 






 §.§ Comparison to simulation-based test-statistics

Existing methods for constructing confidence intervals for h^2 often rely on
simulation. For example, <cit.> uses the parametric bootstrap
to approximate the distribution of the REML estimator of h^2. In comparison to
a naive implementation of the parametric bootstrap, their implementation is very
efficient. Specifically, to estimate the distribution of the REML estimator when
the true heritability is h^2, their method requires (i) drawing N
independent samples from (<ref>) with σ^2 = 1 and null
hypothesis value h^2 given;  and (ii) for each sample, finding an interval
containing a local maximizer of the restricted log-likelihood. Since each sample
is needed for step (ii), <cit.> show that it is sufficient to
draw samples from 𝒩_n(0, I_n) in step (i) and modify step (ii)
accordingly.  The modified step (ii) is performed by evaluating the derivative
of the restricted log-likelihood at fixed grid points. However, for a given
h_2 the derivative must, in general, be evaluated at many grid points; and
this procedure must be repeated for all N samples. Constructing confidence
intervals thus requires repeating this procedure many times, which can be
computationally burdensome. 

In follow-up work, <cit.> proposed a stochastic
approximation to the confidence interval from  <cit.>. This
approximation estimates the upper and lower bounds of the confidence interval
using the Robbins–Monro algorithm. Like our test-statistic, after the same
preprocessing step, the operations required for constructing the intervals from
<cit.> and <cit.> scale linearly with
n. 

In Fig. <ref>, we compare the time needed to compute our interval
to that of <cit.> and
<cit.>. The figure indicates our confidence
intervals can be computed most efficiently of the three methods. For example,
then n = 100 and h^2 = 0.01 our method took on average 3.7 × 10^-4
seconds compared to 3.11 seconds and 45.38 seconds for FIESTA and ALBI,
respectively. When n = 3000 and h^2 = 0.01, our method took 0.07 seconds on
average, whereas FIESTA and ALBI took 511.18 seconds and 1908.68 seconds,
respectively. ALBI and FIESTA are written in Python, while our software is
written in a combination of R and C++. All timings were obtained simultaneously
on the University of Florida's high-performance computing cluster, HiPerGator.

We also compared the width of our confidence intervals to those of
<cit.>.  We found that under a broad range of settings, the
distribution of our interval's width was similar to that of
<cit.>. In our simulations with n = 1000, we recorded the
standardized difference in widths (defined as <cit.>'s width
minus our interval's width, divided by our interval's width). On average the
standardized differences were 0.24, -0.03, and 0.005  when h^2 was 0.01,
0.1, and 0.5, respectively. Median standardized differences were -0.07,
-0.03, and 0.005. 






































 §.§ Gene ranking in spatially-resolved transcriptomics
 
In Section <ref>, we described one application of our method in
the analysis of spatially-resolved transcriptomic data. The rightmost panels of
Fig. <ref> display data generated in a spatially-resolved RNA
sequencing experiment performed on a mouse brain tissue sample slice. In such an
experiment, one measures the expression of thousands of genes (the number of RNA 
reads aligning with each gene) from cells in 50 micron diameter spatial locations (known as
spots) across the tissue sample slice. These locations are represented by the
colored circles in the rightmost panels of Fig. <ref>.

In this context, for one gene, h^2 is the fraction of variability in
expression attributed to the spatial orientation of the tissue sample. In a
typical analysis, one tests h^2 for more than 10,000 genes. Specially,
after normalization, it is common to assume that for the ℓth gene, 

    Y_ℓ∼𝒩_n (X β, σ^2{ h_ℓ^2 K +(1-h_ℓ^2) I_n}),
      ℓ∈{1, …, d},

where h_ℓ^2 is the fraction of spatial variability in the ℓth gene,
n is the number of spots in the tissue, d is the number of genes in the
study, and Y_ℓ and Y_ℓ' are independent for all ℓ≠ℓ'.
Typically K_ij= exp(-s_i-s_j/a) for some a > 0, where s_i is the
two dimensional spatial coordinates of the ith spot and · is the
Euclidean norm. 

To demonstrate the usefulness of our method in this application, we analyzed two
slices of a mouse's brain tissue (anterior and posterior regions), which are
available as the stxbrain data in the  R package.
These data were generated on the 10X Visium platform. Assuming
(<ref>) with a= 0.02 (chosen based on results from
<cit.>), after quality control of filtering out low-expressed
genes and boundary spots in the tissues, n = 2380 for d = 15117 genes in the
anterior slice and n=3011 for d = 14738 genes in the posterior slice.

Using our test statistics, we can provide a confidence interval for each
h_ℓ^2. Moreover, in order to rank genes, we can also provide one-sided
confidence intervals for h_ℓ^2 by comparing U_1(h^2, σ̃^2)
{ℐ^11(h^2, σ̃^2)}^1/2 to quantiles of the standard
normal distribution. 

Because each of the genes has the same K matrix, we can perform a
transcriptome-wide analysis very efficiently. After a one-time
eigendecomposition of K, the remaining processes took around 3 minutes and 4
minutes for the two tissues, respectively. For comparison, ALBI took
approximately 9 minutes and 13 minutes for the two tissues, including the time
needed to find the necessary maximum likelihood estimates of each h_ℓ^2,
which was around 3 minutes. FIESTA, which is generally slower than ALBI when
computing many confidence intervals, timed out after 96 hours.













 §.§ Acknowledgements

The authors thank Yu-Ru Su for a helpful conversation, and thank Regev Schweiger for assistance with , , and for providing helpful references. 





apalike






§ APPENDIX
 



    We first derive a convenient expression for the stochastic part of the
    restricted score function. Recall P = Σ^-1/2 X (X^T
    Σ^-1X)^-1X^T Σ^-1/2 and write Y - Xβ̃ = Y -
    X(X^TΣ^-1X)^-1X^T Σ^-1Y = (I_n - X(X^TΣ^-1X)^-1X^T
    Σ^-1)Y = (I_n - Σ^1/2 P Σ^-1/2)Y. Thus,
    
    (Y - X β̃)^T Σ^-1 V_j Σ^-1 (Y - X β̃)    = Y^T(I_n - Σ^-1/2 P Σ^1/2)Σ^-1V_j Σ^-1(I_n - Σ^1/2 P Σ^-1/2)Y
       = Y^Σ^-1/2(I_n - P)Σ^-1/2V_j Σ^-1/2(I_n - P)Σ^-1/2Y 
       = Y^Σ^-1/2Q Σ^-1/2 V_j Σ^-1/2Q Σ^-1/2Y.

    Since Q Σ^-1/2X = Σ^-1/2X - Σ^-1/2 X (X^T Σ^-1X)^-1X^TΣ^-1X = Σ^-1/2X - Σ^-1/2X = 0, the preceeding display has the same distribution as
    
    ξ_n^ QΣ^-1/2 V_j Σ^-1/2Q ξ_n,

    where ξ_n ∼(0, I_n) and the distribution holds jointly for j ∈{1, 2}. For j = 1 we get
    
    QΣ^-1/2 V_1 Σ^-1/2Q   = Q Σ^-1/2 O σ^2_n(Λ - I_n)O^TΣ^-1/2Q 
        = QO{h^2 Λ + (1 - h^2)I_n }^-1/2(Λ - I_n) {h^2 Λ + (1 - h^2)I_n }^-1/2O Q 
       = Q H Q,

    and, similarly, QΣ^-1/2 V_2 Σ^-1/2Q = Q Σ^-1/2(Σ / σ^2)Σ^-1/2Q = Q / σ^2. It follows that
    
    ℐ_11(h^2, σ^2)    = 1/2(QHQHQ) = 1/2(QHQH) 
    ℐ_12(h^2, σ^2)    = 1/2(Q HQ Q / σ^2)
            = 1/2σ^2(QH) 
    ℐ_22(h^2, σ^2)    = 1/2σ^4(Q) = n - p/2σ^4.




    We start with some useful observations. First, ℐ_11(h^2,
   σ^2) = (QHQ QHQ)/2 and ℐ_12(h^2, σ^2) =
   (QHQ)/(2σ^2) since Q is idempotent. Secondly, the eigenvalues of
   QHQ are real since QHQ is symmetric, and hence the eigenvalues of QHQ
   QHQ are positive as the squares of real numbers.
   
   Now, suppose QHQ has n-p eigenvalues equal to some c∈, with the
   remaining p eigenvalues equal to zero. Then the determinant
   
    |ℐ(h^2, σ^2)| = 1/4σ^4{(n-p)
                (QHQ QHQ) - (QHQ)^2}

   is zero since (QHQQHQ) = (n-p)c^2 and (QHQ) = (n-p)c. Thus,
   ℐ(h^2, σ^2) is singular.

   To establish the other direction, suppose instead that not all of the n-p
   possibly non-zero eigenvalues of QHQ are identical. Then Jensen's
   inequality gives
   
    1/n - p∑_i = 1^n γ_i(QHQ)^2 > {1/n - p∑_i = 1^n γ_i(QHQ) }^2,

   where γ_i(·) is the ith largest eigenvalue. Thus, |ℐ(h^2, σ^2)| > 0. Moreover, at least one eigenvalue of QHQQHQ is strictly positive, so ℐ_11(h^2, σ^2) > 0. Thus, ℐ(h^2, σ^2) is positive definite by Sylvester's criterion.































































 
    If (<ref>) holds for some X with full column rank, then
    
    (n - p)  σ̃^2(h^2) / σ^2_n∼χ^2_n - p.




    Assuming O = I_n, by (<ref>) and proof of Proposition <ref>,
    
    (n - p)  σ̃^2(h^2) / σ^2_n   = (Y - Xβ̃)^Σ^-1 (Y - Xβ̃) 
       = Y^(I_n - Σ^-1/2 P Σ^1/2) Σ^-1 (I_n - Σ^1/2 P Σ^-1/2)Y
       = Y^Σ^-1/2 (I_n-P) Σ^1/2Σ^-1Σ^1/2 (I_n-P) Σ^-1/2 Y
       = Y^Σ^-1/2 (I_n-P) (I_n-P) Σ^-1/2 Y
       = Y^Σ^-1/2 Q Σ^-1/2 Y.

    Since Q Σ^-1/2X = 0 given by proof of Proposition <ref>, Q Σ^-1/2 Y ∼ 𝒩(0, Q), then the above quadratic form has the same distribution as ξ_n^ Q ξ_n, where ξ_n ∼ 𝒩(0, I_n).

    By spectral decomposition, Q = C Ω C^T for some orthogonal C ∈n× n and Ω = (1, …, 1,0, …, 0) with rank n-p. Since C^ξ_n ∼ 𝒩(0, I_n), we can write
    
    ξ_n^ Q ξ_n =ξ_n^ C Ω C^ξ_n 
            = (C^ξ_n)^Ω (C^ξ_n) 
            = ∑_i=1^n-p (C^ξ_n)_i^2 ∼χ^2_n-p,

    where (C^ξ_n)_i is the ith element of C^ξ_n.




    We write ℐ_n = ℐ_n(h_n^2, σ^2_n) and U_n =
    U(h_n^2, σ^2_n) for simplicity. By continuous mapping theorem, it
    suffices to show
    
    𝒥_n^-1/2𝒟_n^-1U_n ⇝(0, I_2),

    where ℐ_n = 𝒟_n 𝒥_n𝒟_n,
    𝒥_n is the correlation matrix corresponding to ℐ_n,
    and 𝒟_n is diagonal with positive entries. By definition,
    𝒥_n is the covariance matrix of 𝒟^-1_n U_n. Also,
    from Proposition <ref>, 𝒟_n11 = ‖
    QHQ‖_F/√(2) and 𝒟_n22 = √((n -
    p))/(√(2)σ^2_n).

    To prove (<ref>), it suffices to show that (i) the
    eigenvalues of 𝒥_n are bounded away from zero and infinity
    asymptotically, and (ii) for any unit-length w∈2,
    
    w^𝒟_n^-1U_n/(w^𝒟_n^-1U_n)^1/2 = w^𝒟_n^-1U_n/(w^𝒥_n w)^1/2⇝(0, 1).

    We start with (ii). By arguments in the proof of Proposition <ref>,
    
    √(2)w^T 𝒟^-1_n U_n    ∼ w_1 ξ_n^ (Q_nH_nQ_n / ‖ Q_n H_n Q_n‖_F) ξ_n - w_1 (Q_n H_n Q_n / ‖ Q_n H_n Q_n‖_F) 
          + w_2 ξ_n^ Q_n ξ_n / √(n - p) - w_2 (Q_n / √(n - p)) 
       = ξ_n^(w_1 Q_nH_n Q_n/‖ Q_n H_n Q_n ‖_F + w_2 Q_n/√(n - p))ξ_n - w_1(Q_nH_n Q_n/‖ Q_n H_n Q_n ‖_F)
           - w_2 √(n - p)
       = ξ_n^ K_n ξ_n - (K_n),

    where K_n is defined by the last equality. Using spectral decomposition of
    K_n and rotational invariance of the normal distribution,
    
    ξ_n^ K_n ξ_n - (K_n)
                ∼∑_i = 1^n-pγ_i(K_n)(ξ_ni^2 - 1).

    We establish (<ref>) by verifying Lyapunov's condition
    <cit.> holds for the right-hand side in the
    last display. Let
    
    a_n = ∑_i = 1^n [{γ_i(K_n)(ξ_ni^2 - 1)}^2] = 
             {(ξ_11^2 - 1)^2}∑_i = 1^n γ_i(K_n)^2 =  
             {(ξ_11^2 - 1)^2}‖ K_n‖_F^2

    and
    
    b_n    = ∑_i = 1^n [{γ_i(K_n)(ξ_ni^2 - 1)}^4] = 
           {(ξ_11^2 - 1)^4}∑_i = 1^n γ_i(K_n)^4 ≤{(ξ_11^2 - 1)^4}‖ K_n‖^2‖ K_n‖_F^2.

    Lyapunov's condition says (<ref>) holds if b_n / a_n^2 →
    0, which holds if and only if ‖ K_n‖ / ‖ K_n‖_F → 0.
   We next show ‖ K_n‖→ 0 and lim inf_n→∞‖ K_n‖_F > 0.

    For the former we have, by the triangle inequality and submultiplicativity
    of the spectral norm,
    
    ‖ K_n‖   ≤| w_1|‖ Q_nH_nQ_n‖/‖ Q_nH_nQ_n‖_F + | w_2|‖ Q_n‖/√(n - p) = | w_1|‖ H_n‖/‖ Q_nH_nQ_n‖_F + | w_2|1/√(n - p)
          =| w_1|γ_1(D_n^2)^1/2/‖ Q_nH_nQ_n‖_F + | w_2|1/√(n - p),

    which tends to zero by assumption 2 since, by
    H.1.h of <cit.>,
    
    ‖ Q_n H_n Q_n‖_F^2    = (Q_nH_n Q_n H_n) 
       = (H_nQ_n H_n - P_n H_nQ_n H_n) 
       = (H_n^2 - H_nP_nH_n - P_nH_n^2 + P_nH_nP_n H_n) 
       ≥(H_n^2) - 2 (P_n H_n^2) 
       ≥(D_n^2) - 2 ∑_i = 1^pγ_i(D_n^2).

    
    To show lim inf_n→∞‖ K_n‖_F > 0 we consider two cases.
    First, if | w_1|≠| w_2|, then by reverse triangle
    inequality,
    
    ‖ w_1Q_n H_n Q_n /‖ Q_n H_n Q_n ‖_F +  w_2Q_n/√(n - p) + ‖_F 
                  ≥|‖ w_1Q_n H_n Q_n /‖ Q_n H_n Q_n ‖_F‖_F -  ‖ w_2Q_n/√(n - p)‖_F|
       = || w_1| - | w_2||,

    which is greater than zero and does not depend on n. If instead |w_1| =
    |w_2|,
    
    ‖ w_1Q_nH_nQ_n/‖ Q_nH_nQ_n‖_F +  w_2Q_n/√(n - p)‖_F^2 
                 = w_1^2 ‖Q_n/√(n - p)±Q_nH_nQ_n/‖ Q_n H_n Q_n‖_F‖_F^2 
       = w_1^2 ( Q_n^2/n - p + 
                (Q_nH_nQ_n)^2/‖ Q_n H_n Q_n‖_F^2± 2 
                Q_n/√(n - p) Q_n H_n Q_n/‖ Q_n H_n Q_n ‖_F)
        = 2 w_1^2  {1 ±(Q_n/√(n - p) Q_n H_n Q_n/‖ Q_n H_n Q_n ‖_F)},

    so, since Q_n is idempotent, it
    suffices to show
    
    lim sup_n→∞|(Q_n H_n)|/√(n - p)‖ Q_n H_n Q_n‖_F < 1.

    We bound the numerator first. By the triangle inequality,
    
    |(Q_n H_n)| = |(H_n) - (P_n H_n)|≤|(H_n)| + |(P_n H_n)|.

    Applying H.1.g and H.1.h of <cit.> (note their comment regarding positive
    semi-definiteness) to (P_n H_n), we have
    
    ∑_i = 1^pγ_n - p + i( D_n) ≤(P_nH_n) ≤∑_i = 1^pγ_i( D_n),

    which gives an upper bound of | (P_n
        H_n)|. Using this bound followed by Jensen's inequality gives
    
    | (P_n
            H_n)|   ≤max{|∑_i = 1^pγ_i( D_n)|, 
            |∑_i = 1^pγ_n - p + i( D_n) |}
       ≤∑_i = 1^p√(γ_i(D_n^2))
       ≤(p∑_i = 1^pγ_i(D_n^2))^1/2.

    We thus have that the ratio to be bounded is no greater than
    
    |(D_n)| + 
            {p∑_i = 1^pγ_i(D_n^2))}^1/2/√(n - p){(D_n^2) - 2 ∑_i = 1^pγ_i(D_n^2)}^1/2 =
            |(D_n)|/√(n - p){(D_n^2) - 
                2 ∑_i = 1^pγ_i(D_n^2)}^1/2 + o(1).

    By assumption 2, and since n / (n - p) → 1 by assumption 1,
    the last right-hand side is
    
    {1 + o(1)}n^-1/2|(D_n)|/(D_n^2)^1/2 + o(1),

    the upper limit of which is less than one by assumption 3, which establishes
    (<ref>).
    
    It remains only to establish (i). First, since 𝒥_n is a
    correlation matrix, it's entries are no greater than one in absolute value,
    and hence γ_1(𝒥_n) = ‖𝒥_n‖≤‖𝒥_n‖_F ≤ 2. To bound γ_2(𝒥_n) away from
    zero, suppose for contradiction lim inf_n→∞γ_2(𝒥_n) = 0. Then we can pick a subsequence along
    which γ_2(𝒥_n) → 0. Since ‖𝒥_n‖_F is
    bounded, we may also, by the Bolzano–Weierstrass property, pick the
    subsequence in such a way that 𝒥_n converges to some positive
    semi-definite 𝒥. Since γ_2(𝒥_n) → 0,
    𝒥 has at least one vanishing eigenvalue. Let w be a
    corresponding eigenvector. Then, along the subsequence, w^𝒥_n w → 0. But, by the arguments following
    (<ref>), w^𝒥_n w =  ‖ K_n‖_F^2,
    and we already proved lim inf_n→∞‖ K_n‖_F > 0, which
    gives the desired contradiction. Thus, (i) holds and the proof is complete.
 




    Define a remainder R_n1 to be dealt with later by
    
    U_n1(h_n^2, σ̃^2_n) = U_n1(h_n^2, σ^2_n) + ∇^2_12l_R(h_n^2, σ_n^2) (σ̃^2_n - σ^2_n) + R_n1.

    Since
    
    U_n2(h^2_n, σ^2_n) = -n - p/2σ^2_n + (n - p) σ̃^2_n/2σ^4_np,

    we have σ̃^2_n -  σ^2_n =  2 σ^4_np
    U_n2(h_n^2, σ^2_n) / (n - p) = U_n2 / ℐ_n22(h_n^2, σ^2_n). For the remainder of the proof, we omit
    the arguments when they are the true parameters (h_n^2, σ^2_n). Then
    
    U_n1(h_n^2, σ̃^2_n)    = U_n1 + (∇^2_12l_R)  U_n2 / ℐ_n22 + R_n1
            =      U_n1 -ℐ_n12ℐ_n22^-1 U_n2 + R_n1 + R_n2,

    where R_n2 = (∇^2_12l_R + ℐ_n12)
    U_n2/ℐ_n22. We complete the proof by showing (i) R_n1^2ℐ_n^(11)→ 0 in probability,
    (ii) R_n2^2 ℐ_n^(11)→ 0 in probability, and (iii)
    
    {U_n1 - ℐ_n12 U_n2 /  ℐ_n22}
            (ℐ_n^(11))^1/2⇝(0, 1).

    It will be useful to note
    
    ℐ_n^(11) = (ℐ_n11 - ℐ_n21^2/ℐ_n22)^-1 = ℐ_n22/ℐ_n11ℐ_n22 - ℐ_n21^2,

    which in fact does not depend on σ_n^2. We begin with (ii) and have
    
    ∇_12^2 l_R = ∇_1 U_n2   = - 1/2 σ^4_np∑_i = 1^n (λ_i - 1)(Y_i - X_i^Tβ̃)^2/(h_n^2 λ_i + 1 - h_n^2)^2
        = -1/2σ^2_n(Y - Xβ̃)^Σ_n^-1/2 D_n Σ_n^-1/2 (Y - Xβ̃)

    From the proof of Proposition <ref>, Y - Xβ̃ = (I_n - Σ_n^1/2PΣ^-1/2)Y = Σ_n^1/2(I_n - P_n)Σ_n^-1/2Y = Σ_n^1/2Q_n Σ_n^-1/2 Y, and so
    
    ∇_12^2 l_R    = -1/2σ^2_n Y^Σ_n^-1/2Q_n D_n Q_n Σ_n^-1/2Y.

    Thus, -R_n2(ℐ_n^(11))^1/2 is
    
    1/2 n^-1/2σ^2_n(ℐ_n11ℐ_n22 - ℐ_n12^2)^1/21/n^1/2(Y^Σ_n^-1/2Q_n D_n Q_n Σ_n^-1/2Y - (Q_n D_n) ) U_n2 / √(ℐ_n22),

    which we can write as I^-1/2×II×III by defining
    
    I   = 4 n^-1σ^4_np(ℐ_n11ℐ_n22 - ℐ_n12^2) = (Q_nD_nQ_nD_n) {1 - p / n - n^-1(Q_nD_n)^2/(Q_nD_nQ_nD_n)}; 
    II   = 1/n^1/2(Y^Σ_n^-1/2Q_n D_n Q_n Σ_n^-1/2Y - (Q_n D_n) ); 
    III   = U_n2/√(ℐ_n22).

    By calculations in the proof of Lemma <ref>, the
    term of I in curly brackets is bounded away from zero
    asymptotically (c.f. condition 3 of that lemma), so I≥ c
    (Q_nD_nQ_nD_n) for some c > 0 and all large enough n. Next,
    II has mean zero and variance (Q_nD_nQ_nD_n) / n. Thus,
    I^-1/2×II has mean zero and variance no greater
    than 1 / (cn), and hence it tends to zero in mean square and probability.
    By Chebyshev's inequality, III is O_(1), and hence
    I^-1/2×II×III→ 0 in
    probability, which completes the proof of (ii).

    To show (i), use that, by Taylor's theorem with Lagrange-form remainder, for some
    σ̅^2_np between σ^2_n and σ̃_np^2,
    
    R_n1 = ∇^3_122l_R(h_n^2, σ̅^2_np)
                (σ̃_np^2 - σ_n^2)^2/2.

    Using calculations similar to when deriving ∇^2_12l_R, we get
    
    ∇^3_122l_R(h_n^2, σ̅^6_np) =
            1/2σ̅^6_np∑_i = 1^n
            (λ_i - 1)(Y_i - X_i^Tβ̃)^2/(h_n^2 λ_i + 1 - h_n^2)^2
            = σ^2_n/2σ̅^6_np Y^Σ_n^-1/2Q_n D_n Q_n Σ_n^-1/2Y.

    Thus, 2 R_n1 (ℐ_n^(11))^1/2 is
    
    σ_n^6/σ̅^6_np1/2 σ^4_npℐ_n22^1/2/(ℐ_n11ℐ_n22 - ℐ_n12^2)^1/2 Y^Σ_n^-1/2Q_n D_n Q_n Σ_n^-1/2Y (U_n2 / ℐ_n22)^2 
       = σ_n^6/σ̅^6_npn^-1/2/2σ^2_nn^-1/2(ℐ_n11ℐ_n22 - ℐ_n12^2)^1/2 Y^Σ_n^-1/2Q_n D_n Q_n Σ_n^-1/2Y U_n2^2/ℐ_n22ℐ_n22^1/2/ℐ_n22σ^2_n
       =  σ_n^6/σ̅^6_npn^-1/2/I^1/2Y^Σ_n^-1/2Q_n D_n Q_n Σ_n^-1/2Y ×III^2 ×√(2)/√(n - p)
       =  σ_n^6/σ̅^6_np×I^-1/2× (II + (Q_n D_n)/√(n)) ×III^2 ×√(2)/√(n - p).

    Note |σ̅^2_np / σ^2_n - 1|≤|σ̃^2_n / σ^2_n - 1| and
    σ̃^2_n/σ^2_n∼χ^2_n-p/(n-p), which tends to 1 in
    probability. Thus, σ̅^2_np / σ^2_n tends to one in
    probability and, by the arguments used to establish (ii), it suffices to
    show |(Q_nD_n)/n |×I^-1/2→ 0. But
    Jensen's inequality applied to the sum of  the n real eigenvalues of QDQ, which are the same as those of QD, gives |(Q_nD_n)/n |≤(Q_n
    D_n Q_n D_n)^1/2 n^-1/2 and in the proof of (ii) we argued I≥  c (Q_nD_nQ_nD_n), so
    
    |(Q_nD_n)/n |×I^-1/2≤ (cn)^-1/2→ 0.


    Finally, to show (iii), let a_n = (ℐ_n^(11))^1/2 [1,
    -ℐ_n12 / ℐ_n22]^, so that what we want to show
    is a_n^ U_n ⇝(0, 1). A direct calculation shows the
    variance of a_n^ U_n is a_n^ℐ_n a_n = 1. Let b_n =
    ℐ_n^1/2 a_n. Then ‖ b_n‖ = 1 and a_n^ U_n =
    b_n^ℐ_n^-1/2U_n. The proof of Lemma <ref>
    shows ℐ_n^-1/2U_n ⇝(0, I_2). Moroever, by
    the Bolzano–Weierstrass Theorem and the subsequence principle, we may
    assume b_n → b for some b with ‖ b‖ = 1. The Cramér–Wold
    Theorem says b^ℐ_n^-1/2U_n ⇝(0, 1),
    and therefore b_n^ℐ_n^-1/2U_n = b^ℐ_n^-1/2U_n + (b_n - b)^ℐ_n^-1/2U_n = b^ℐ_n^-1/2U_n + o_(1) ⇝(0, 1) by Slutsky's
    Theorem, which completes the proof.




    
        Condition 1 is the same as in Lemma <ref>, so it suffices to
        verify the remaining conditions of Lemma <ref>. First, condition
        3 ensures
        
    lim inf_n→∞min_i ∈{1, …, n}{h_n^2(λ_ni-1)+1}>0,

        so together with condition 2 we have, for some c_max∈ (0, ∞),
        
    lim sup_n→∞√(γ_1(D^2_n)) = lim sup_n→∞max_i ∈{1, …, n}|λ_ni-1/h_n^2(λ_ni-1)+1| ≤ c_max.

        
       Condition 2 also implies we may upon increasing c_max if needed have, for any λ_ni and λ_nj,
        
    0 < (h_n^2(λ_ni-1)+1)(h_n^2(λ_nj-1)+1) 
                     ≤max{1,λ_ni}max{1,λ_nj}≤ c_max^2.

       Thus, if |λ_ni - λ_nj| > ϵ,
        
    |λ_ni-1/h_n^2(λ_ni-1)+1-λ_nj-1/h_n^2(λ_nj-1)+1|
                       =|λ_ni- λ_nj|/(h_n^2(λ_ni-1)+1)(h_n^2(λ_nj-1)+1)
       > ϵ/c^2_max.

        Thus, if two diagonal elements of Λ_n are more than ϵ
        apart, then the corresponding diagonal elements of D_n are more than
        ϵ / c_max^2 apart.
        
        Applying (<ref>) with λ_nj = 1 shows the number of
        eigenvalues of D_n in [-ϵ/(c_max^2), ϵ /
        (c_max^2)] is at most k_n^1, the number of eigenvalues of
        Λ_n in [1-ϵ,1+ϵ]. Consequently, together with
        condition 1 we have
        
    ∑_i=1^pγ_i(D_n^2)/∑_i=1^n γ_i(D_n^2)≤p c_max^2/(n - k^1_n) ϵ^2 / c_max^4→ 0,

        which verifies condition 2 of Lemma <ref>.

        Let γ̅_n = 1/n∑_i=1^nγ_i(D_n) and let
        r_n be the number of eigenvalues of D_n in [γ̅_n -
        ϵ / (4 c^2_max), γ̅_n + ϵ / (4 c^2_max)]. Then
        
    ∑_i=1^nγ_i(D^2_n)-nγ̅_n^2
                     = ∑_i=1^n{γ_i(D_n)-γ̅_n}^2 
                     ≥ϵ^2/16c^4_max(n - r_n).

        We are done if lim sup_n→∞ r_n / n < 1 since, then, condition
        3 of Lemma <ref> is satisfied:
        
    lim sup_n→∞(D)^2/n (D^2) = 1 - 
                lim inf_n→∞∑_i=1^nγ_i(D^2_n)
                -nγ̅_n^2/∑_i=1^nγ_i(D^2_n)≤ 1 
                - lim inf_n→∞ϵ^2 n /16 c_max^6(n - r_n) < 1.


        To see r_n / n is indeed bounded away from one, observe γ̅_n is
        bounded by (<ref>), so we may, by the
        Bolzano–Weierstrass property and passing to a subsequence if necessary,
        assume γ̅_n converges to some γ_0 ∈. Thus, for
        all large enough n, all the eigenvalues of D_n in [γ̅_n -
        ϵ / (4 c^2_max), γ̅_n + ϵ / (4 c^2_max)]
        are also in [γ_0 - ϵ / (2 c^2_max), γ_0 + ϵ
        / (2 c^2_max)]. But then, by (<ref>), the corresponding
        λ_ni are in [λ_0 - ϵ, λ_0 + ϵ] for
        some fixed λ_0. For example, we can fix any i and n such
        that  γ_i(D_n) ∈ [γ_0 - ϵ / (2 c^2_max),
        γ_0 + ϵ / (2 c^2_max)] and let λ_0 be the
        corresponding λ_ni. Thus, r_n ≤ k_n^λ_0 and the
        proof is completed since lim sup_n→∞k_n^λ_0/n < 1
        by condition 4.
    
    








































































        






























































































        

















        








































































































    

























































