As a prominent dimension reduction method for multivariate linear regression, the envelope model has received increased attention over the past decade due to its modeling flexibility and success in enhancing estimation and prediction efficiencies. Several enveloping approaches have been proposed in the literature; among these, the partial response envelope model [57] that focuses on only enveloping the coefficients for predictors of interest, and the simultaneous envelope model [14] that combines the predictor and the response envelope models within a unified modeling framework, are noteworthy. In this article we incorporate these two approaches within a Bayesian framework, and propose a novel Bayesian simultaneous partial envelope model that generalizes and addresses some limitations of the two approaches. Our method offers the flexibility of incorporating prior information if available, and aids coherent quantification of all modeling uncertainty through the posterior distribution of model parameters. A block Metropolis-within-Gibbs algorithm for Markov chain Monte Carlo (MCMC) sampling from the posterior is developed. The utility of our model is corroborated by theoretical results, comprehensive simulations, and a real imaging genetics data application for the Alzheimer’s Disease Neuroimaging Initiative (ADNI) study.
Bayesian envelope modelMultivariate regressionReducing subspaceSimultaneous envelopePartial envelopeImaging geneticsU.S. National Science Foundation DMS-2013486DMS-1712418Wisconsin Alumni Research FoundationUniversity of Wisconsin-MadisonAlzheimer’s Disease Neuroimaging Initiative U01 AG024904W81XWH-12-2-0012National Institute on AgingNational Institute of Biomedical Imaging and BioengineeringAbbVieAlzheimer’s Drug Discovery FoundationBioClinica, Inc.BiogenBristol-Myers Squibb CompanyEli Lilly and CompanyF. Hoffmann-La Roche LtdGenentech, Inc.FujirebioGE HealthcareLundbeckMeso Scale Diagnostics, LLC.Novartis Pharmaceuticals CorporationPfizer Inc.ServierTakeda Pharmaceutical CompanyFoundation for the National Institutes of HealthC. Zhang’s work was partially supported by U.S. National Science Foundation grants DMS-2013486 and DMS-1712418, and provided by the University of Wisconsin-Madison Office of the Vice Chancellor for Research and Graduate Education with funding from the Wisconsin Alumni Research Foundation. Park’s research is partially supported by University of Wisconsin-Madison Office of the Vice Chancellor for Research and Graduate Education. Data collection and sharing for this project was funded by the Alzheimer’s Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904) and DOD ADNI (Department of Defense award number W81XWH-12-2-0012). ADNI is funded by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering, and through generous contributions from the following: AbbVie, Alzheimer’s Association; Alzheimer’s Drug Discovery Foundation; Araclon Biotech; BioClinica, Inc.; Biogen; Bristol-Myers Squibb Company; CereSpir, Inc.; Cogstate; Eisai Inc.; Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann-La Roche Ltd and its affiliated company Genentech, Inc.; Fujirebio; GE Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy Research & Development, LLC.; Johnson & Johnson Pharmaceutical Research & Development LLC.; Lumosity; Lundbeck; Merck & Co., Inc.; Meso Scale Diagnostics, LLC.; NeuroRx Research; Neurotrack Technologies; Novartis Pharmaceuticals Corporation; Pfizer Inc.; Piramal Imaging; Servier; Takeda Pharmaceutical Company; and Transition Therapeutics. The Canadian Institutes of Health Research is providing funds to support ADNI clinical sites in Canada. Private sector contributions are facilitated by the Foundation for the National Institutes of Health (www.fnih.org). The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer’s Therapeutic Research Institute at the University of Southern California. ADNI data are disseminated by the Laboratory for Neuro Imaging at the University of Southern California.Introduction
Consider the standard multivariate linear regression
Y=μY+βTX+ϵY∣X,
where the responses Y∈Rr and the centered predictors X∈Rp are both multivariate, the means of responses μY∈Rr, the regression coefficients β∈Rp×r, and the random error vector ϵY∣X∼Nr(0,ΣY∣X). It has applications in broad fields of scientific research, including economics [6], chemistry [35], agriculture [51], engineering [44], bioinformatics [34], etc. Specifically, we consider the application to the imaging genetics problem (IGP).
Thanks to the rapid technological development in medical imaging and genome sequencing, the diagnosis, prevention and treatment of mental illness have been greatly improved. Compared with traditional medical imaging analysis and genome-wide association studies (GWAS) that make use of only imaging or only genome data and were studied to a large extent [67, 63], another more promising yet more challenging study is to directly relate the imaging phenotypes (IPs) as responses to genotypes (GPs) as predictors of interest, while also considering the effects of some demographic covariates or other risk factors. This is called IGP and is less explored in the statistical community [42]. Because of the Central Dogma that is often stated as DNA makes RNA and RNA makes protein, genetic markers should be promising in revealing the abnormal protein expression in cortical or sub-cortical structures of brain, whose volumes should be more informative quantitative traits than simple disease status.
Some early practices in studying IGP focus on univariate [54] or voxel-wise type of methods [25], which explores the relationship between every (voxel, SNP) pair (SNP: Single-nucleotide polymorphism), or between each individual voxel and all genetic markers together, respectively. These types of methods enjoy simplicity and can handle high dimensional imaging or genomic datasets, but cannot leverage spatial dependence among voxels. To this end, several multivariate regression models have been proposed to jointly relate tens of brain summary measures, like the volumes of Regions of Interest (ROIs), to hundreds of SNPs [65, 68, 23, 64, 46]. Among them, [46] is the only work that employs the envelope model to IGP, to the best of our knowledge.
As a dimension reduction method developed in the past decade, the first envelope model [15] aims at providing an efficient estimator of the regression coefficients β in the multivariate linear regression (1.1), by assuming that there are linear combinations of the responses Y that are not related to the regression, and removing its variation as immaterial information (see Section 2.2 for more details). With a great success in achieving efficient estimation and modeling flexibility, the original envelope model has been adapted or extended to predictor envelope [12], groupwise envelope [46], partial response envelope [57], simultaneous envelope [14], sparsity learning [56], matrix-valued [18] or tensor-valued [39] responses or predictors, quantile regression [19], generalized linear models [11], etc. The envelope methods were recently reviewed by [37] and more details were introduced in [13]. The aforementioned envelope methods are all developed from the frequentist perspective. However, the literature addressing Bayesian envelope methods is rare, but still interesting for its capabilities of incorporating prior information into inference, and implementing posterior uncertainty quantification, compared with the unnatural prior information incorporation and the bootstrap or asymptotic variance in the frequentist setting.
[33] introduced the first framework for the Bayesian envelope model. However, it is hard to apply this framework to other envelope contexts, since it requires the specific formulation of the response envelope to construct the conjugacy. Furthermore, sampling from the generalized matrix Bingham distributions and the truncated inverse gamma distributions is required in this framework, which makes the MCMC algorithm computationally expensive. Later, [5] proposed a new flexible and handy Bayesian framework for the envelope model, which bypasses the manifold restriction on the basis matrix of the envelope space by an unconstrained matrix re-parameterization. Since all parameters under this new framework are either vectors, or unconstrained, or positive definite matrices, the computation is much more efficient without sampling from the generalized matrix Bingham distributions and the truncated inverse gamma distributions.
The rapid development in the envelope family requires a model to unify various perspectives, and a more clear choice to practitioners on which envelope model to apply, especially when we have a particular interest for a subset of predictors that could be either continuous or discrete or a mix of them. The existing simultaneous envelope model [14] integrates both the response and predictor envelopes within a unified modeling framework. It simultaneously enjoys the benefits from two envelope components, i.e., more efficient estimation of the regression coefficients offered by the response envelope, and improved prediction efficiency of the responses obtained by the predictor envelope, but it cannot focus on the predictors of main interest, and therefore cannot leverage the idea from the partial response envelope [57] to obtain further estimation efficiency. More importantly, without additional assumptions imposed, it is guaranteed to outperform the predictor or the response envelope only for the Normal predictors ([14], Propositions 2 and 5). However, the discrete predictors will not follow the Normal distribution and including them in the simultaneous envelope by implicitly regarding them as Normal will lead to biased estimates and invalid inference. The partial response envelope [57] focuses on enveloping only the coefficients of the predictors of main interest X1 (where in (1.1), X is partitioned into (X1T,X2T)T, and X1 denotes the predictor of main interest to us and X2 are nuisance covariates), for example the genetic markers in IGP. Compared with full response envelope [15] that indistinguishably envelopes the coefficients of all predictors X, enhanced estimation efficiency is possible by this model since the envelope space could be shrunk if only concentrating on the coefficients of X1, but it could not leverage the benefits of the predictor envelope [12], which is known to have the potential to increase prediction efficiency. The envelope model that is proposed in this paper could address these limitations while still providing a valid inference procedure.
To address the aforementioned limitations, and improve the efficiencies for the estimation of the regression coefficients of main interest (the coefficients of X1, which we denote as β1 in (3.2)) and the prediction of the responses (Y in (3.2)), this paper proposes a new envelope model, called the simultaneous partial envelope model (3.7)–(3.10), in Section 3. The proposed envelope model combines the partial response envelope model [57] and the simultaneous envelope model [14] within a more general modeling framework, by further partitioning the predictors of main interest X1 into X1C and X1D, the continuous part and the discrete part. Our proposed envelope model simultaneously imposes the partial response envelope structure on the coefficients of X1, and the partial predictor envelope only on the coefficients of X1C, instead of the coefficients of whole X1. By this construction, our proposed model could address the aforementioned limitations for the simultaneous and the partial response envelopes. Firstly, two envelope components that are simultaneously imposed on our model are partial envelopes, which both focus on the predictors of interest only. Secondly, it avoids the Normal requirement mentioned above for the discrete predictors in the simultaneous envelope model while still providing a valid inference, by only imposing the partial response envelope (i.e. no partial predictor envelope) on the coefficients of X1D. It is also worthy to note that X1C is Normal with means depending on X1D in our model. Lastly, to further improve the estimation and prediction efficiencies upon the partial response envelope, the simultaneous envelope structure is considered in our model for the coefficients of X1C, which turns out to enhance not only the estimation of the coefficients of X1C, but also that of the coefficients of X1D in our simulation (see the comparisons on the estimation of the coefficents of X1D between our method and the Bayesian partial response envelope in Table 3 and Figure 2), as a possible synergetic effect.
Our proposed model includes a number of important models in the envelope family as special cases, including the (partial) response envelope, the (partial) predictor envelope and the simultaneous envelope (see Table 1 for their relationship). In Section 4, we develop an efficient Bayesian inference procedure based on [5], which allows convenient prior information incorporation and variability quantification, compared with the frequentist envelope methods.
The contribution of this paper is fivefold.
We propose a new simultaneous partial envelope model (3.7)–(3.10) that unifies several existing envelope models as special cases. Our proposed model allows both predictors of interest X1 and nuisance covariates X2 to be either continuous or discrete or a mix of them. Compared with any envelopes that our model degenerates to, our method provides a more efficient while still valid estimator for the regression coefficients of X1 and improved prediction for Y in (3.2). These advantages address the limitations that are mentioned above for the simultaneous envelope model [14] and the partial response envelope model [57]. The corresponding block Metropolis-within-Gibbs algorithm is developed for the posterior inference.
We establish the theoretical properties including the posterior propriety for our model and the Harris ergodicity for our algorithm.
We are the first to investigate the performance of several popular dimension selection methods together, among the Bayesian envelope literature, to the best of our knowledge.
We apply our method to an imaging genetics analysis for the ADNI study. We show the improved prediction accuracy of our method compared with all other competitors, and obtain some scientifically meaningful findings from posterior selection. By incorporating the weak imaging genetics relationship as prior information, the prediction accuracy is further improved and the associated shrinkage effect is also investigated.
We build the R package SIMP that implements our algorithm, which is available at https://github.com/yanbowisc/SIMP.git.
Some notations frequently used in this paper are summarized here. For any a=1,2,…, let PS be a projection matrix onto a subspace S⊆Ra and PS⊥=Ia−PS be the projection matrix onto S⊥, where S⊥ denotes the orthogonal complement of S and Ia denotes the a-dimensional identity matrix. Let S+a×a denote the class of all a×a real valued positive definite matrices. Let 1n denote a vector of all ones with length n. For a matrix A∈Ra×a and a subspace E⊆Ra, AE is defined as AE={Aϵ:ϵ∈E}. Notations |A|, ‖A‖, vec(A) and span(A) denote the determinant, the spectral norm, the vectorization by columns and the column space of the matrix A respectively. For a1,a2,…,an∈R and a square matrix S∈Ra×a, we describe a diagonal matrix by either specifying its diagonal elements by the notation diag(a1,…,an), or using the notation diag{S} to indicate that the diagonal elements of the square matrix S will be taken out to constitute this diagonal matrix. Moreover, Sii denotes the i-th diagonal element of the square matrix S. For a random vector Y, the distribution of Y∣X is interpreted as the distribution of Y given the fixed value of X for a non-stochastic vector X, or the distribution of Y conditional on X for a random vector X. For a,b,c=1,2,…, and random vectors X∈Ra, Y∈Rb, we use Cov((X,Y)∣Z) to denote the a×b dimensional covariance matrix of X and Y conditional on a random (or given a non-stochastic) vector Z∈Rc. If multiple random (or non-stochastic) vectors are conditional on (or given), they are bracketed together after the vertical line, like Y∣(·) and Cov((X,Y)∣(·)). Also, := and X=DY indicate equating by definition and random vectors X and Y are equally distributed respectively. Lastly, for any x,y∈R, x≪y indicates that x is much less than y, while x≫y means that x is much greater than y.
The rest of this paper is organized as follows. Section 2 reviews two existing envelope models, which are prototypes of two envelope components of our proposed model. Section 3 introduces the definition and the formulation of our proposed simultaneous partial envelope model. Section 4 develops a Bayesian inference procedure for our proposed model, with detailed procedure of MCMC algorithm left to Appendix A. Section 5 establishes the theoretical properties of our model and algorithm. Section 6 investigates the issue of model selection for our Bayesian envelope model. Sections 7 and 8 conduct comprehensive simulation studies and investigate a real data application respectively. Section 9 concludes our paper by further discussions.
Review of the Response and the Predictor Envelope ModelsPreliminaries
In this section, we introduce two definitions and one proposition, which are necessary for the formal introduction of the response envelope model [15] and the predictor envelope model [12]. For convenience, the definition of reducing subspace [10] is firstly introduced.
([<xref ref-type="bibr" rid="j_nejsds23_ref_010">10</xref>], Chapter II, Definition 3.5).
A subspaceE⊆Rais called a reducing subspace ofΣ∈Ra×afor some integera=1,2,…, ifΣE⊆EandΣE⊥⊆E⊥.
If E is a reducing subspace of Σ, there is a crucial decomposition of Σ as illustrated in the following proposition.
Eis a reducing subspace ofΣif and only ifΣcould be decomposed asΣ=PEΣPE+PE⊥ΣPE⊥, wherePEandPE⊥are the projection matrices onto subspacesEandE⊥.
Next, for a positive definite matrix M, we introduce the general definition of theM-envelope of a subspaceS. It will be used in defining the response and the predictor envelope models in Section 2.2 and Section 2.3 and our proposed model in Section 3.
Let a matrixM∈S+a×aand a subspaceS⊂span(M). TheM-envelope ofSis the intersection of all reducing subspaces ofMthat containS, and is denoted byEM(S).
For convenience, we will frequently use the notation EM(S) to define envelope structures in the rest of this paper. Below we briefly review the response envelope model [15] in Section 2.2, and mention how this idea was adapted to define the predictor envelope model [12] in Section 2.3.
The Response Envelope Model
In this section, we review the definitions and the coordinate form of the response envelope model [15], which aims to improve the estimation efficiency of β in (1.1) by assuming some linear combinations of Y do not depend on X. Here, we assume X in (1.1) to be non-stochastic. The efficiency gains come from removing those redundant variation in Y by estimating the associated linear combination coefficients. Formally, the response envelope model is defined by the smallest subspace EY∣X that satisfies Condition 1.
([<xref ref-type="bibr" rid="j_nejsds23_ref_015">15</xref>], p. 928).
Assume in (1.1), the subspaceSY∣X⊂Rrsatisfies
PSY∣X⊥Y∣X=DPSY∣X⊥Y,
Cov((PSY∣XY,PSY∣X⊥Y)∣X)=0.
Condition 1(a) assumes the distribution of PSY∣X⊥Y does not depend on the value of X. Meanwhile, Condition 1(b) excludes the indirect effect of X on PSY∣X⊥Y via PSY∣XY, by assuming PSY∣XY and PSY∣X⊥Y are uncorrelated given X. Therefore, PSY∣X⊥Y only contains immaterial information, and PSY∣XY carries all the material information and possibly some extra immaterial information for the regression. Instead of using any satisfied SY∣X, the definition of the response envelope model ensures the uniqueness of EY∣X by taking it to be the smallest one, i.e. the intersection of all satisfied subspaces. At the same time, by using EY∣X rather than any SY∣X, immaterial information is guaranteed to be not contained in PEY∣XY anymore, and could be removed to the largest degree. So we call PEY∣XY the material part and PEY∣X⊥Y the immaterial part of Y respectively. Although intuitive, it is not easy to formulate the response envelope model under Condition 1. [15] derived the following Condition 2, which is equivalent to Condition 1. The response envelope model could be equivalently defined by Condition 2 with subspace EY∣X, for (1.1).
ΣY∣X=PSY∣XΣY∣XPSY∣X+PSY∣X⊥ΣY∣XPSY∣X⊥, i.e.SY∣Xis a reducing subspace ofΣY∣X.
Condition 2 assumes that SY∣X is a reducing subspace of ΣY∣X that contains span(βT). Recalling Definition 2, EY∣X is exactly EΣY∣X(span(βT)), i.e. the ΣY∣X-envelope of span(βT). Therefore, the response envelope model could be defined directly by EΣY∣X(span(βT)) for (1.1).
Suppose that the orthonormal basis matrices of EY∣X and EY∣X⊥ are Γ and Γ0 respectively. Then according to Condition 2 with subspace EY∣X, the coordinate form of the response envelope model is Y=μY+ΓηX+ϵY∣X,ΣY∣X=ΓΔΓT+Γ0Δ0Γ0T, where η, Δ=ΓTΣY∣XΓ and Δ0=Γ0TΣY∣XΓ0 carry the coordinates of βT relative to Γ, ΣY∣X relative to Γ and Γ0 respectively. The response envelope model could provide significant efficiency gains when ‖Δ‖≪‖Δ0‖, since there will be a large amount of immaterial information for removal under this scenario.
The Predictor Envelope Model
In this section, we review the predictor envelope model [12] by showing how the idea of the response envelope could be adapted to reduce the predictor space. The predictor envelope model has the potential to improve the prediction of Y in (1.1).
To define the predictor envelope, we assume X to be stochastic following N(0,ΣX). Then (1.1) is called the predictor envelope model if the envelope structure EΣX(span(β)) is imposed. Unlike the response envelope, EΣX(span(β)) is the intersection of all reducing subspaces of ΣX that contain span(β). Similarly, the material and immaterial parts of X could be defined. The immaterial part of X is assumed not to carry the information for the regression. Suppose the orthonormal basis matrices of EΣX(span(β)) and its orthogonal complement subspace are Υ and Υ0, then the coordinate form of the predictor envelope model is Y=μY+ψTΥTX+ϵY∣X,ΣX=ΥΞΥT+Υ0Ξ0Υ0T, where ψ, Ξ=ΥTΣXΥ and Ξ0=Υ0TΣXΥ0 carry the coordinates of β relative to Υ, ΣX relative to Υ and Υ0 respectively. The large efficiency gains could be obtained when ‖Ξ‖≫‖Ξ0‖.
Proposed Simultaneous Partial Envelope Model
As mentioned in Section 1, suppose that our predictors X∈Rp could be partitioned into three parts
X=X1X2=X1CX1DX2,
where X1=(X1CT,X1DT)T∈Rp1 is the predictors of main interest to us in the multivariate linear regression, with its continuous and discrete parts being X1C∈RpC and X1D∈RpD respectively (p1=pC+pD), while X2∈Rp2 denotes the set of predictors that is not of main interest (p=p1+p2). By separating X1D with X1C, we could avoid the Normality restrictions for X1D and hence whole X1, while still providing a valid simultaneous envelope estimator for the coefficients of X1C. Assuming predictors are not centered, the standard multivariate linear model can be written as
Y=μY+β1CT(X1C−μ1C)+β1DT(X1D−μ1D)+β2T(X2−μ2)+ϵY∣X,
where μY∈Rr, μ1C∈RpC, μ1D∈RpD and μ2∈Rp2 denote the unknown means of Y, X1C, X1D and X2 respectively, β1C∈RpC×r, β1D∈RpD×r and β2∈Rp2×r denote the unknown regression coefficients of X1C, X1D and X2, and ϵY∣X is the random error vector. Specifically, combining β1C and β1D, β1=(β1CT,β1DT)T denotes the coefficients of main interest to us. Assume X1C to be stochastic, and X1D and X2 to be non-stochastic, with sample means X‾1D and X‾2. Replacing μ1D and μ2 with X‾1D and X‾2 respectively, (3.1) is equal to (3.2)
Y=μY+β1CT(X1C−μ1C)+β1DT(X1D−X‾1D)+β2T(X2−X‾2)+ϵY∣X.
For the purpose of expositional simplicity, we implicitly assume the discrete predictors inX1DandX2to be quantitative variables in (3.1), (3.2) and later in (3.7). However, it is worthy to note that categorical predictors, either nominal or ordinal, are still applicable to our model by including their dummy variables intoX1DandX2instead. In either case, the centering for these two predictors serves to ensure thatμYcould be interpreted as the expectation (conditional on model parameters) of the sample means ofY.
Relationship between SIMP (3.7)–(3.10) and several other envelope models. The meanings of r, pC, pD, p1 and p2 are introduced in (3.1), and dX and dY are dimensions of EC∣D and EY∣1 respectively. The effective number of parameters for each model is also listed. Note for the response and the partial response envelopes, X1D is required to be non-stochastic only and not necessary to be discrete.
Conditions for (3.7)–(3.10)
Effective number of parameters
SIMP
N.A.
dY(dX+pD)+r(r+2p2+3)/2+pC(pC+3)/2
response envelope [15]
pC=p2=0
dYp+r(r+3)/2
predictor envelope [12]
pD=p2=0 and dY=r
r(r+2dX+3)/2+p(p+3)/2
simultaneous envelope [14]
pD=p2=0
dYdX+r(r+3)/2+p(p+3)/2
partial response envelope [57]
pC=0
dYp1+r(r+2p2+3)/2
partial predictor envelope [45]
p2=0 and dY=r
r(r+2dX+2pD+3)/2+pC(pC+3)/2
We firstly impose some distributional assumptions on (3.2). Assume ϵY∣X∼Nr(0,ΣY∣X). Given X1D, assume that X1C=μ1C+γT(X1D−X‾1D)+ϵC∣D for an unknown γ∈RpD×pC, where ϵC∣D∼NpC(0,ΣC∣D) and is independent of ϵY∣X. Rigorously speaking, μ1C in (3.2) should only be interpreted as the expectation of X‾1C (rather than the expectation of X1C) conditional on model parameters under this assumption, due to the non-stochasticity of X‾1D. However, we will keep using μ1C to avoid abusing notations.
To provide the efficiency gains on estimating β1 and predicting Y, we combine the advantages of the response envelope on the estimation and the predictor envelope on the prediction, by further imposing the following two partial envelope structures on (3.2) (by Definition 2) simultaneously:
Partial predictor envelope: EΣC∣D(L), i.e. the ΣC∣D-envelope of L and is shortened to EC∣D, with dimension dX(d≤dX≤pC),
Partial response envelope: EΣY∣X(R), i.e. the ΣY∣X-envelope of R and is shortened to EY∣1, with dimension dY(d1≤dY≤r),
where L and R denote span(β1C) and span(β1T), and d and d1 denote rank(β1C) and rank(β1) respectively. (3.2) is called the Simultaneous Partial Envelope Model (SIMP), if these distributional assumptions and envelope structures are imposed. Note that SIMP avoids the Normality requirement for the whole X1 (i.e., no other assumptions on X1D besides the non-stochasticity, and X1C is Normal with means depending on X1D, which relaxes the identical Normal distribution assumption of each observation of predictors in [14]) by imposing EY∣1 on the row space of whole β1 while EC∣D on the column space of only β1C.
It is worthy to note that SIMP is a broad class of envelopes that degenerates to several popular envelope models under certain conditions as listed in Table 1. Below, we introduce Conditions 3 and 4 in Section 3.1, which could equivalently define SIMP and will offer more intuitions on where the efficiency gains come from by assuming these two envelope components. For the convenience of developing the Bayesian inference procedure in Section 4, we give the coordinate form of SIMP in (3.7)–(3.10) in Section 3.3. For the preparation of Section 3.3, Section 3.2 will introduce a re-parameterization trick, to avoid the direct Bayesian inference on matrix parameters living in the Stiefel manifold.
Equivalent Conditions for <inline-formula id="j_nejsds23_ineq_296"><alternatives><mml:math>
<mml:mi mathvariant="normal">SIMP</mml:mi></mml:math><tex-math><![CDATA[$\mathrm{SIMP}$]]></tex-math></alternatives></inline-formula>
In this section, we introduce Conditions 3 and 4, which can equivalently define SIMP for (3.2). These two conditions will not be used in the formulation of SIMP in Section 3.3, but will offer some intuitions on the source of efficiency gains. The equivalence between the conditions in this section and the envelope structures that we have given for SIMP (by Definition 2) shares the same reason with the equivalence between Conditions 1 and 2 for the response envelope model in Section 2.2.
Assume in (3.2), the subspaceSC∣D⊂RpCsatisfies,
Cov((Y,PSC∣D⊥X1C)∣(PSC∣DX1C,X1D,X2))=0,
Cov((PSC∣DX1C,PSC∣D⊥X1C)∣X1D)=0.
Assume in (3.2), the subspaceSY∣1⊂Rrsatisfies,
PSY∣1⊥Y∣(X1,X2)=DPSY∣1⊥Y∣X2,
Cov((PSY∣1Y,PSY∣1⊥Y)∣(X1,X2))=0.
Condition 3 excludes the direct and indirect partial effects (via PSC∣DX1C) of PSC∣D⊥X1C on Y. The intersection of all subspaces that satisfy Condition 3 gives the same partial predictor envelope space EC∣D as defined previously by EΣC∣D(L). We call PEC∣DX1C and PEC∣D⊥X1C the material part and immaterial part of X1C. Given X1D and X2, the immaterial part of X1C is assumed not to affect Y by Condition 3. Similarly, the partial response envelope space EY∣1 previously defined by EΣY∣X(R) is also the intersection of all subspaces that satisfy Condition 4. We call PEY∣1Y and PEY∣1⊥Y the material part and immaterial part of Y. Given X2, the immaterial part of Y is assumed not to be affected by X1.
Re-Parameterization of the Basis Matrices
We only consider the case of 0<dX<pC and 0<dY<r in this section. When dX∈{0,pC} or dY∈{0,r}, we can simply take the orthonormal bases of envelope space EC∣D or EY∣1 and the orthogonal complement subspace to be either null or the identity matrix. When 0<dX<pC and 0<dY<r, let L∈RpC×dX, L0∈RpC×(pC−dX), R∈Rr×dY and R0∈Rr×(r−dY) denote the orthonormal bases of EC∣D, EC∣D⊥, EY∣1 and EY∣1⊥, respectively.
To avoid the difficulty of direct Bayesian inference on L, L0, R and R0, which all live in the Stiefel manifold, re-parameterization of L, L0, R and R0 is considered in the coordinate form of SIMP in Section 3.3. We illustrate the re-parameterization of L and L0 here, and this idea can be similarly exploited to re-parameterize R and R0. Let L1 be the matrix formed by the upper dX rows of L, and L2 be the matrix that contains the remaining rows. Without loss of generality, we assume that L1 is nonsingular. Otherwise, we could reorder the rows of L to make L1 nonsingular. Then,
L=L1L2=IdXAL1:=CdX(A)L1,
where A=L2L1−1 and CdX(A)=IdXA. [56] shows that A depends on L only through span(L) and there is a one-to-one correspondence between EC∣D and A. [8] further shows that if CdX(A) is a basis of EC∣D, then
DpC−dX(A)=−ATIpC−dX
is a basis of EC∣D⊥. Therefore, after normalization,
L(A)=CdX(A)(CdX(A)TCdX(A))−1/2
and
L0(A)=DpC−dX(A)(DpC−dX(A)TDpC−dX(A))−1/2
are a pair of orthonormal bases of EC∣D and EC∣D⊥ respectively. The re-parameterization
R(B)=CdY(B)(CdY(B)TCdY(B))−1/2
and
R0(B)=Dr−dY(B)(Dr−dY(B)TDr−dY(B))−1/2
by an unconstrained matrix B∈R(r−dY)×dY can be similarly constructed.
Coordinate Form of SIMP
Like (2.1)–(2.2) for the response envelope model and (2.3)–(2.4) for the predictor envelope model, we intend to give the coordinate form of SIMP explicitly as in (3.7)–(3.10), i.e., formulate two envelope structures EC∣D and EY∣1 by analytical expressions for β1C, β1D, ΣC|D and ΣY|X, and correspondingly rewrite (3.2), for the preparation of the Bayesian inference in Section 4.
Let β1C=UDVT be the singular value decomposition, where U∈RpC×d and V∈Rr×d are semi-orthogonal matrices, and D=diag(λ1,…,λd), with λi≥0, i=1,…,d, being d singular values of β1C where d=rank(β1C). Since L=span(β1C)=span(U) and R=span(β1T)=span(β1CT,β1DT)=span(V,β1DT), according to the definitions of EC∣D and EY∣1 by Definition 2 at the start of Section 3, U=LO and V=RP for some O∈RdX×d and P∈RdY×d. Hence the coordinate form of β1C is
β1C=(LO)D(RP)T:=LηCTRT,
where ηC=PDOT∈RdY×dX is the coordinate of β1CT relative to R and LT. Also,
β1DT=RηD,
for some ηD∈RdY×pD, which carries the coordinate of β1DT relative to R. Then the coordinate form of SIMP is given by Y=μY+R(B)ηCL(A)T(X1C−μ1C)+R(B)ηD(X1D−X‾1D)+β2T(X2−X‾2)+ϵY∣X,X1C=μ1C+γT(X1D−X‾1D)+ϵC∣D,ΣC∣D=L(A)ΩL(A)T+L0(A)Ω0L0(A)T,ΣY∣X=R(B)ΦR(B)T+R0(B)Φ0R0(B)T, where β1C=L(A)ηCTR(B)T and β1D=ηDTR(B)T, as illustrated in (3.5) and (3.6). This implies that the rows and columns of β1C depend on the partial predictor and partial response envelopes only, and the columns of β1D depend on the partial response envelope only. Note that Ω=L(A)TΣC∣DL(A), Ω0=L0(A)TΣC∣DL0(A), Φ=R(B)TΣY∣XR(B) and Φ0=R0(B)TΣY∣XR0(B) carry the coordinates of ΣC∣D relative to L(A) and L0(A), and ΣY∣X relative to R(B) and R0(B) respectively. Multiply R(B)T or R0(B)T on both sides of (3.7),
R(B)TY=R(B)TμY+ηC(L(A)T(X1C−μ1C))+ηD(X1D−X‾1D)+R(B)Tβ2T(X2−X‾2)+R(B)TϵY∣X,R0(B)TY=R0(B)TμY+R0(B)Tβ2T(X2−X‾2)+R0(B)TϵY∣X.
This pair of equations shows that X2 affects both material and immaterial parts of Y, X1D affects only the material part of Y, and only the material part of X1C affects the material part of Y only. The estimation of β1C is benefitted from removing the redundant variations of both R0(B)TY and L0(A)TX1C, whereas the efficiency gains for estimating β1D only come from removing the redundant variation of R0(B)TY.
Inherited from the (partial) predictor and the (partial) response envelope models, the partial predictor envelope component of SIMP offers large efficiency gains when ‖Ω‖≫‖Ω0‖, and significant advantages of the partial response envelope component require ‖Φ‖≪‖Φ0‖.
Bayesian Inference
In this section, we develop a Bayesian procedure for the statistical inference of SIMP, hence our method is called the Bayesian simultaneous partial envelope model. It is implemented through sampling from the posterior distribution of parameters. The determination of the posterior distribution requires both the likelihood function and prior distributions.
Suppose that we have n independent observations D={Yi,X1C,i,X1D,i,X2,i}i=1n from SIMP. Let Y=(Y1,…,Yn)T∈Rn×r, X1C=(X1C,1,…,X1C,n)T∈Rn×pC, X1D=(X1D,1,…, X1D,n)T∈Rn×pD and X2=(X2,1,…,X2,n)T∈Rn×p2 be data matrices. For notational simplicity, we consider standardized data matrices Y˜=Y−1nμYT, X˜1C=X1C−1nμ1CT, X˜1D=X1D−1nX‾1DT and X˜2=X2−1nX‾2T for Y, X1C, X1D and X2, respectively. Then, for fixed envelope dimensions dX and dY, the log-likelihood function for SIMP (3.7)–(3.10) is given by
l(Θ)=const−n2log(|Ω|)−n2log(|Ω0|)−12tr{(X˜1C−X˜1Dγ)(L(A)ΩL(A)T+L0(A)Ω0L0(A)T)−1(X˜1C−X˜1Dγ)T}−n2log(|Φ|)−n2log(|Φ0|)−12tr{(Y˜−X˜1CL(A)ηCTR(B)T−X˜1DηDTR(B)T−X˜2β2)(R(B)ΦR(B)T+R0(B)Φ0R0(B)T)−1(Y˜−X˜1CL(A)ηCTR(B)T−X˜1DηDTR(B)T−X˜2β2)T},
where const denotes a constant which does not depend on parameters of the model, Θ={μ1C, μY, β2, γ, ηC, ηD, A, B, Ω, Ω0, Φ, Φ0}.
To perform Bayesian inference on Θ related to SIMP, we put the following prior distributions on Θ:
π(μ1C,μY)∝1,β2∣Φ,Φ0,B∼MNp2,r(M−1Z,M−1,ΣY∣X),γ∣Ω,Ω0,A∼MNpD,pC(Λ−1F,Λ−1,ΣC∣D),Ω∼IWdX(ΨX,wX),Ω0∼IWpC−dX(ΨX0,wX0),Φ∼IWdY(ΨY,wY),Φ0∼IWr−dY(ΨY0,wY0),A∼MNpC−dX,dX(A0,KA,ΣA),B∼MNr−dY,dY(B0,KB,ΣB),ηC∣Φ∼MNdY,dX(WE−1,Φ,E−1),ηD∣Φ∼MNdY,pD(TQ−1,Φ,Q−1),
where MN and IW denote Matrix normal and Inverse-Wishart distributions. See Appendix E for the detailed introduction of these two distributions. Here we implicitly assume A0∈R(pC−dX)×dX, B0∈R(r−dY)×dY, Z∈Rp2×r, F∈RpD×pC, W∈RdY×dX, T∈RdY×pD, wX>dX−1, wX0>pC−dX−1, wY>dY−1, wY0>r−dY−1 and M∈S+p2×p2, Λ,Q∈S+pD×pD, ΨX,ΣA,E∈S+dX×dX, ΨY,ΣB∈S+dY×dY, ΨX0,KA∈S+(pC−dX)×(pC−dX), ΨY0,KB∈S+(r−dY)×(r−dY). The prior information for the parameters could be easily incorporated through these hyperparameters. For example, if we know a priori that the partial predictor envelope EC∣D is likely to be EˆC∣Dprior, then we could compute Aˆprior from the basis of EˆC∣Dprior by (3.3), and set A0 as Aˆprior. Our confidence about the prior of A is encoded in KA and ΣA. Specifically, the prior covariance matrices for the i-th row and the j-th column of A are KA,iiΣA and ΣA,jjKA respectively, where KA,ii and ΣA,jj are the i-th and j-th diagonal elements of KA and ΣA, and the small KA,ii’s or ΣA,jj’s can reflect our strong prior belief of the specific row(s) or column(s) of Aˆprior. The Bayesian inference is made via the Metropolis-within-Gibbs algorithm to generate the sample from the posterior distribution of the parameters Θ. A number of posterior samples at early iterations are discarded as burn-in. Suppose the retained posterior samples are D∗={Θ(s)}s=1S:={μ1C(s), μY(s), β2(s), γ(s), ηC(s), ηD(s), A(s), B(s), Ω(s), Ω0(s), Φ(s), Φ0(s)}s=1S with sample size S. The algorithm is provided in Appendix A, and the warm start initial estimator that we propose for the MCMC algorithm is detailed in Appendix B.
Theoretical Properties
In this section, we establish two theoretical properties of our method. We first establish the propriety of the joint posterior density of all parameters in Theorem 1, although the priors of μ1C and μY are improper. Proofs are left to Appendix D.
The posterior density of(μ1C,μY,β2,γ,ηC,ηD,A,B,Ω,Ω0,Φ,Φ0)with respect to the Lebesgue measure onRPC×Rr×Rp2×r×RpD×pC×RdY×dX×RdY×pD×R(pC−dX)×dX×R(r−dY)×dY×S+dX×dX×S+(pC−dX)×(pC−dX)×S+dY×dY×S+(r−dY)×(r−dY)is proper.
Theorem 2 establishes the Harris ergodicity of the Metropolis-within-Gibbs algorithm developed for the Bayesian inference of SIMP. This property ensures the asymptotic convergence of the Markov chain generated by our algorithm to the joint posterior distribution, irrespective of the starting points. This convergence may fail on a set of measure zero without the guarantee of this property. Despite measure zero, some choices for starting points in this pathological null set can arise naturally, as illustrated in [48] for example. Technical details including definitions for probabilistic terminologies and proofs are left to Appendix D.
Whenever0≤dX≤pCand0≤dY≤r, the Markov chain generated by the Metropolis-within-Gibbs algorithm for the posterior sampling ofΘinSIMP, as detailed in Appendix A, is Harris ergodic, i.e. (a) ϕ-irreducible with respect to some measure ϕ, (b) Aperiodic, and (c) Harris recurrent.
Selection of Envelope Dimensions
To effectively search the envelope dimensions, we first select the minimum possible dimension of two envelope components, which is the rank of β1C, and narrow down the range of dimensions to be searched to save the computational cost.
To determine d, we adapt the Bura-Cook estimator [4] to our model setting. This estimator includes a sequence of Chi-squared tests. For k=0,1,…,min(pC,r)−1, the test statistic is Tk=n∑j=k+1min(pC,r)φj, where φ1≥⋯≥φmin(pC,r) are eigenvalues of
βˆ1C,std={(n−pC−1)/n}1/2S1C1/2βˆRY∣1D,2∣1CSRY∣1D,2∣1C−1/2,
where S1C is the sample covariance matrix of X1C. The residuals of the regression from Y on X1D and X2 are regressed on X1C again, and the estimated regression coefficients and residual sample covariance matrix are βˆRY∣1D,2∣1C and SRY∣1D,2∣1C respectively. Hereinafter, if we don’t mention specifically for the method of a regression, it always refers to the frequentist ordinary least squares. [4] shows Tk follows the χ(pC−k)(r−k)2 distribution under the null hypothesis that rank(β1C) is equal to k. For each k=0,1,…,min(pC,r)−1, we compare Tk with the upper α quantile of the χ(pC−k)(r−k)2 distribution, where α is the pre-specified significance level, and default to be 0.05 in our numerical studies. We choose d as the first non-significant value of k or d=min(pC,r) otherwise. The Bura-Cook estimator exhibits excellent ability in selecting d, as illustrated in the third column of Table 2. Details of the simulation set-up and numerical performance are introduced in Section 7.1.
Selecting Envelope Dimensions
Dimensions dX and dY should be specified before fitting SIMP. In Section 7.1, we investigate the performance of Bayesian cross-validation (CV) and four information criteria, Akaike information criterion (AIC-MCMC), Bayesian information criterion (BIC-MCMC), Deviance information criterion (DIC) and Watanabe-Akaike information criterion (WAIC) in selecting these two envelope dimensions for SIMP. They are selected for their popularity in the Bayesian literature. Dimension selection is critical in envelope modeling. To the best of our knowledge, we are the first to investigate the performance of these interesting methods together for the dimension selection of the Bayesian envelope model, while previous Bayesian envelope literature simply chooses one of them directly [5, 38]. For each fixed dX and dY, AIC−MCMC=−2l(ΘˆdX,dYmax)+2K(dX,dY),BIC−MCMC=−2l(ΘˆdX,dYmax)+log(n)K(dX,dY), where l(ΘˆdX,dYmax) is the largest log-likelihood value attained by MCMC samples after burn-in, with envelope dimensions fixed at dX and dY for associated parameters, and K(dX,dY)=dY(dX+pD)+r(r+2p2+3)/2+pC(pC+3)/2 is the effective number of free parameters in SIMP. Note that the definitions of AIC-MCMC and BIC-MCMC here are slightly different from the ones that readers might be used to in the frequentist setting, where we replace the maximum likelihood estimates of the parameters by the empirical maximizers in the retained MCMC samples. Definitions of DIC, WAIC and Bayesian CV involve more notations, and hence are left to Appendix F. We choose dX and dY as the minimizers of either the minus of the average out-of-sample estimate of the log predictive density for the Bayesian CV (see details in Appendix F), or one of these four information criteria, among each pair of dX=dˆ,…,pC and dY=dˆ,…,r by the grid search, where dˆ is the estimate of d from the Bura-Cook estimator introduced in Section 6.1.
Percentages that the Bura-Cook estimator correctly identifies d (Column 3), and each of five methods correctly identifies dX or dY individually or the (dX,dY) pair (Columns 4-8) out of 500 repetitions. Among five methods, the best one for selecting (dX,dY) pair under each scenario is highlighted in bold face.
(dX,dY)
n
B-C
AIC-MCMC
BIC-MCMC
DIC
WAIC
5-fold CV
(2, 2)
150
0.99
0.37/0.65/0.29
0.85/0.98/0.84
0.00/0.20/0.00
0.00/0.24/0.00
0.58/0.56/0.26
300
1.00
0.35/0.63/0.25
0.92/0.99/0.91
0.00/0.19/0.00
0.00/0.20/0.00
0.54/0.61/0.27
500
1.00
0.44/0.68/0.35
0.96/0.99/0.95
0.00/0.15/0.00
0.00/0.24/0.00
0.53/0.63/0.30
(6, 2)
150
0.98
0.51/0.79/0.41
0.42/0.98/0.40
0.07/0.15/0.01
0.54/0.30/0.18
0.08/0.09/0.01
300
0.99
0.59/0.84/0.49
0.66/0.99/0.66
0.03/0.16/0.00
0.55/0.28/0.15
0.08/0.06/0.01
500
0.99
0.61/0.87/0.53
0.75/0.99/0.75
0.02/0.13/0.00
0.53/0.27/0.15
0.15/0.05/0.01
(2, 6)
150
0.99
0.50/0.79/0.39
0.99/0.99/0.98
0.00/0.36/0.00
0.04/0.71/0.03
0.41/0.01/0.01
300
1.00
0.55/0.81/0.45
0.99/1.00/0.99
0.00/0.31/0.00
0.02/0.72/0.01
0.40/0.08/0.04
500
1.00
0.63/0.81/0.51
1.00/1.00/0.99
0.00/0.31/0.00
0.01/0.69/0.01
0.44/0.16/0.08
(6, 6)
150
0.99
0.90/0.92/0.83
0.99/0.99/0.99
0.12/0.44/0.05
0.89/0.91/0.82
0.64/0.82/0.52
300
1.00
0.91/0.94/0.86
1.00/1.00/1.00
0.08/0.51/0.05
0.90/0.90/0.81
0.61/0.76/0.44
500
0.99
0.89/0.95/0.85
0.99/0.99/0.99
0.05/0.46/0.01
0.87/0.88/0.78
0.57/0.71/0.41
Simulation StudyPerformance of the Envelope Dimensions Selection
In this section, we investigate the numerical performance of the dimension selection methods for SIMP, which are introduced in Section 6.2. We fix r=pC=8, pD=p2=2, and we consider four cases where true dimensions (dX,dY) are (2,2), (6,2), (2,6) or (6,6), with their corresponding effective numbers of parameters being 112, 120, 128 and 152. Three sample sizes, 150, 300 and 500, are considered, given the magnitude of the effective numbers of parameters. We generate Ω, Ω0, Φ and Φ0 all as diagonal matrices, with their associated L1,1 norms being 100, 0.5, 1 and 10, and the assignment proportions to each diagonal element are generated from Dirichlet distributions with shape vectors of all 5, 5, 1 and 1 respectively. Entries of β2, γ, ηC, ηD are generated independently from Unif(−2,2), and those of A and B are from Unif(−1,1) independently. All elements of μ1C and μY are independently from Unif(0,10). Covariates of X1D are generated independently from discrete Unif{0,1,2}, and samples for X2 are independently from Np2(μ2,Σ2) with μ2=(2,5)T, and Σ2 as a realization from IWp2(Ip2,p2).
For our method, we specify the hyperparameters Z, F, W, T, A0 and B0 to be zero matrices, M, Λ, E and Q to be 10−6 times the identity matrices, KA, KB, ΣA and ΣB to be 106 times the identity matrices, wX, wX0, wY and wY0 to be the row numbers of Ω, Ω0, Φ and Φ0 respectively, and ΨX, ΨX0, ΨY and ΨY0 to be identity matrices multiplied by 10−6 and their respective degrees of freedom (correspondingly, wX, wX0, wY and wY0). The MCMC algorithm for our method is run for 20,000 iterations with the first 50% iterations as burn-in. For each case, their performance is averaged over 500 repetitions. Percentages that the Bura-Cook estimator correctly identifies d and each of five methods that correctly identifies true dX or dY individually or true (dX,dY) jointly are listed in Table 2.
Results in Table 2 reveal the advantage of BIC-MCMC in selecting dX and dY for SIMP, among various settings of true dX and dY and sample sizes, except when dX=6, dY=2 and the sample size is small. It suggests that the excellent performance of BIC-MCMC may need at least moderate sample sizes. AIC-MCMC and WAIC also perform relatively well in the large dX and large dY setting. These observations for BIC-MCMC that needs large sample size to respond, and the better performance of AIC-MCMC under large envelope dimensions are consistent with the previous findings for envelope models ([58]). Compared with underestimation, a little over-estimating of these two envelope dimensions is acceptable, since only some estimation efficiency, rather than the material information, is lost [5]. The right panel of Figure 1 illustrates this phenomenon, since it is observed that the estimation biases (i.e., the square root of the vertical gaps between the green and the red lines, due to the Bias-variance decomposition) are large for underestimated envelope dimensions, whereas with correct or overestimated envelope dimensions, the biases are nearly zero and the estimated variances only increase linearly. Although all of AIC-MCMC, DIC and WAIC show some tendency in over-estimating dX and dY in Table 2, the most significant one is DIC (see the left panel of Figure 1). Hence, AIC-MCMC and WAIC may also be considered in practice, especially for small sample sizes or the large dX and large dY setting. The Bayesian CV chooses envelope dimensions based on the prediction performance. Although its performance in selecting the true envelope dimensions is far from the best in our simulation, it is still recommended in practice if the prediction performance is our utmost concern and the computational burden is manageable.
Left: Selection frequencies of dX and dY (Rows 1 and 2) of five criteria, for different settings of true (dX,dY) (Columns 1–4). Colors represent the selection frequency; Right: MSE and estimated variances (Var) of βˆ1 from SIMP across all possible choices of dY (Row 1) or dX (Row 2) with dX or dY fixed at 2 (the true value) respectively. The horizontal blue dotted line in the figure that is at the bottom of the right panel shows the MSE of βˆ1 from the frequentist partial response envelope with the envelope dimension as 2. Here the results of two panels are based on the sample size as n=300 over 500 repetitions. Results with n=150 and n=500 are similar.
Estimation PerformanceSetting
In this section, we compare the estimation performance of SIMP with the Bayesian partial predictor envelope (PX-env), the Bayesian partial response envelope (PY-env) and several other classical envelope methods or multivariate regression approaches, including the frequentist predictor envelope (FX-env) [12], the frequentist response envelope (FY-env) [15], the frequentist simultaneous envelope (FS-env) [14], the frequentist partial response envelope (FPY-env) [57], the frequentist ordinary least squares (FOLS), the principal component regression (PCR) [32], the partial least squares regression (PLSR) [16], the canonical correlation analysis (CCA) [26] and the reduced rank regression (RRR) [29]. Note in Sections 7.2 and 8, PX-env and PY-env simply refer to SIMP with dY and dX fixed at r and pC respectively, rather than the Bayesian versions of the partial predictor and partial response envelopes as exactly indicated in Table 1, however their differences are small (see the simulation results between FPY-env and PY-env in Sections 7.2.2–7.2.5) but will facilitate the comparison. We choose the posterior mean as the point estimator here and later in Section 8 for all the Bayesian (envelope) approaches. It is a decision-theoretic estimator and could be obtained easily once the posterior samples are available. In comparison, the posterior median and the maximum a posteriori (MAP) estimators are not considered in our paper, since they are expected to perform similarly (see Table C.4 in Appendix C for a numerical evidence of the posterior median estimator) due to the bell-shaped posterior distribution that is observed in both simulation studies (see Appendix G) and the real data application (see Appendix H.2), and the MAP estimator requires additional computation for the optimization.
In this section, we keep almost all set-ups from Section 7.1 except the following. We test on nine cases, for each combination of r=3,8 or 15 with n=300,500 or 1000. The covariates of X1D are generated independently from discrete Unif{0,1}. Dimensions dX and dY are both fixed at 2. The specification of the hyperparameters (if still existed), the number of iterations and the burn-in proportion of the MCMC algorithm for PX-env, PY-env and SIMP in this section are the same as those for SIMP in Section 7.1.
To make comparison fair, for competitors that cannot account for the partial structures (i.e., all methods excluding SIMP, PX-env, PY-env and FPY-env), Y is regressed on X2 first, and the fitted residuals are the actual responses Yˇ for the regression on our interested predictors X1. This is a common way to adjust the effects of prognostic factors that are not of main interest for IGP [68]. Tuning parameters or envelope dimensions for all methods are assumed to be unknown. The envelope dimensions for all Bayesian envelope methods (SIMP, PX-env, PY-env) are chosen by BIC-MCMC, and those for the frequentist envelope methods (FPY-env, FX-env, FY-env and FS-env) are chosen by the traditional BIC. The numbers of components for PCR and PLSR are selected based on minimizing the Mean squared prediction error (MSPE) from 5-fold CV. For k=1,2,…,5, let Ck represent the k-th index set among five partition sets of n samples. For ℓ=1,2,…,p1, we determine the number of components as the ℓ that minimizes
MSPE(ℓ)=1n∑k=1K∑i∈Ck∑j=1r(yˇij−yˇˆij−k(ℓ))2.
For the sample i in Ck, yˇij represents the j-th response adjusted by X2, while yˇˆij−k(ℓ) is the prediction of yˇij by X1C,i and X1D,i, for either PCR or PLSR trained on all samples except those in Ck, and with number of components as ℓ. The rank for RRR is determined by the original Bura-Cook estimator [4]. As the estimation performance of μ1C, μY and β2 is very close among the competitors that we choose, and not of main interest to us, their results are not reported for compactness. For the true data generating mechanisms, we consider three models (M1-M3) from SIMP, and one model (M4) from RRR which is to investigate the model mis-specification.
{Ω,Ω0,Φ,Φ0}={10IdX,IpC−dX,IdY,5Ir−dY},
{Ω,Ω0,Φ,Φ0}={IdX,10IpC−dX,IdY,5Ir−dY},
{Ω,Ω0,Φ,Φ0}={10IdX,IpC−dX,5IdY,Ir−dY},
β1C=GCH, β1D=GDH, where GC∈RpC×d1, GD∈RpD×d1, H∈Rd1×r, with d1=1 and their elements are all generated independently from Unif(−1,1). Samples of X1C are generated from NpC(μ1C,ΣC), and ΣY∣X and ΣC are one realizations of IW(Ir,r) and IW(IpC,pC) respectively.
Note that we allow the true data generating models to be SIMP with ‖Ω‖≫‖Ω0‖ and ‖Φ‖≪‖Φ0‖ under (M1) or ‖Ω‖≪‖Ω0‖ and ‖Φ‖≪‖Φ0‖ under (M2) or ‖Ω‖≫‖Ω0‖ and ‖Φ‖≫‖Φ0‖ under (M3). (M4) assumes the true model to be RRR and is designed to test the robustness of SIMP under model mis-specification. Results under all scenarios are displayed from Section 7.2.2 to Section 7.2.5.
MSE comparison between SIMP and other 11 competitors for estimating β1C and β1D over 500 repetitions under the data generating mechanism (M1). The lowest MSE under each combination of r and n is in bold face.
r
n
FOLS
RRR
PCR
PLSR
CCA
FX-env
FY-env
FS-env
FPY-env
PX-env
PY-env
SIMP
3
300
β1C
0.16
0.07
0.16
0.17
3.95
0.10
0.05
0.05
0.05
0.02
0.05
0.02
β1D
0.55
0.27
0.63
0.69
13.52
0.98
0.18
0.33
0.16
0.23
0.17
0.08
500
β1C
0.09
0.04
0.09
0.09
3.81
0.05
0.03
0.02
0.03
0.01
0.03
0.01
β1D
0.34
0.16
0.34
0.36
12.87
0.38
0.10
0.14
0.09
0.14
0.10
0.05
1000
β1C
0.04
0.02
0.04
0.04
3.80
0.02
0.01
0.01
0.01
0.00
0.01
0.00
β1D
0.16
0.07
0.16
0.17
13.68
0.17
0.05
0.06
0.05
0.07
0.05
0.02
8
300
β1C
0.69
0.10
0.30
0.27
23.17
0.23
0.06
0.07
0.05
0.03
0.05
0.02
β1D
2.58
0.39
3.23
3.14
86.62
2.77
0.19
0.52
0.18
0.93
0.18
0.11
500
β1C
0.40
0.05
0.28
0.27
22.80
0.15
0.03
0.03
0.03
0.02
0.03
0.01
β1D
1.53
0.21
3.14
2.96
85.83
1.57
0.12
0.21
0.11
0.57
0.11
0.06
1000
β1C
0.20
0.02
0.21
0.22
22.82
0.08
0.02
0.01
0.01
0.01
0.02
0.01
β1D
0.74
0.09
1.5
1.57
84.69
0.89
0.05
0.07
0.05
0.27
0.05
0.03
15
300
β1C
1.45
0.16
0.16
0.11
30.38
0.10
0.09
0.13
0.08
0.06
0.08
0.05
β1D
4.74
0.38
1.39
1.37
103.70
1.31
0.19
0.80
0.15
1.87
0.15
0.17
500
β1C
0.86
0.09
0.16
0.10
30.18
0.08
0.05
0.06
0.04
0.04
0.04
0.03
β1D
2.86
0.21
1.34
1.32
103.07
1.31
0.10
0.42
0.09
1.13
0.09
0.09
1000
β1C
0.42
0.04
0.14
0.11
30.20
0.07
0.02
0.03
0.02
0.02
0.02
0.01
β1D
1.39
0.09
1.29
1.24
102.97
1.30
0.05
0.15
0.04
0.54
0.04
0.04
Average PSD across each coordinate of vec(β1C) (1-8r), followed by vec(β1D) ((8r+1)-10r), over 500 repetitions. Rows and columns indicate different sample sizes and dimensions of responses.
Under (M1), the mean squared errors (MSE) for β1C and β1D, two parameters that we are most interested in, are reported in Table 3. For almost every combination of r and n, SIMP always performs the best for the estimation of both β1C and β1D. PX-env or PY-env, two special cases of SIMP by setting dY=r or dX=pC, can achieve competitive performance for either only β1C, or whole β1 but is still less efficient than SIMP. Comparing these two models, we find the former obtains a slightly better performance for β1C while the latter achieves a much better performance for β1D. This observation is consistent with our assumptions for SIMP in (3.7), where we have assumed both the partial predictor and partial response envelope structures for β1C with a slightly stronger signal from the former component here (‖Ω‖/‖Ω0‖>‖Φ0‖/‖Φ‖), and only the partial response envelope structure is imposed on β1D. Although PY-env and SIMP shares the same partial response envelope structure for β1D, it is surprising that the latter outperforms the former for the estimation of β1D when r=3 and 8, which is possibly due to the reason that the more efficient estimator of β1C improves the estimation of R(B), and hence the estimation of β1D in SIMP. Meanwhile, it is noteworthy that PY-env and FPY-env show similar estimation performance for both β1C and β1D under all cases, but the former one allows more convenient posterior variability quantification and utilization of the prior information. FX-env, FY-env, FS-env and other five non-envelope methods cannot outperform their partial envelope counterparts and any partial envelope methods respectively.
To further illustrate the advantage of SIMP from the perspective of posterior variability, the average of the posterior standard deviations (PSD) of each coordinate of vec(β1C) and vec(β1D) is calculated, as the average of empirical standard deviations of posterior samples of vec(β1C) and vec(β1D) across 500 repetitions. The performance on average PSD is compared between SIMP, Bayesian linear regression (BLR) by implementing SIMP with dX=pC and dY=r, PX-env and PY-env under (M1). The envelope dimensions are assumed to be unknown and selected by BIC-MCMC for all three envelope methods. The results are reported in Figure 2. In terms of both β1C and β1D, three envelope methods can always significantly outperform BLR under all cases. SIMP and PX-env obtain the lowest average PSD for β1C, and are better than PY-env especially when r=3. This is as expected since we have imposed both partial predictor and partial response envelope structures on β1C but the signal of the former is stronger. As r becomes 8 or 15, the performance of PY-env is close to other two envelope methods. In terms of β1D, SIMP achieves the lowest PSD under all scenarios. When r=8 or 15, the performance of PY-env is close to SIMP, and is significantly better than PX-env. This is as expected as well, since the partial response envelope structure is imposed on β1D for both SIMP and PY-env, and the slight advantage of SIMP for the estimation of β1D is possibly due to the more efficient estimation of β1C, which is one synergetic effect of imposing two envelope components simultaneously. Although no envelope assumption is imposed on β1D for PX-env, the posterior variability for the estimator of β1D from PX-env is still smaller than that of BLR.
MSE comparison between SIMP and other 11 competitors for estimating β1C and β1D over 500 repetitions under the data generating mechanism (M2). The lowest MSE under each combination of r and n is in bold face.
r
n
FOLS
RRR
PCR
PLSR
CCA
FX-env
FY-env
FS-env
FPY-env
PX-env
PY-env
SIMP
3
300
β1C
0.06
0.05
0.06
0.06
1.31
0.05
0.04
0.02
0.04
0.06
0.04
0.03
β1D
0.32
0.46
0.34
0.35
7.11
0.33
0.47
0.12
0.45
0.31
0.48
0.12
500
β1C
0.04
0.03
0.04
0.04
1.29
0.03
0.01
0.01
0.01
0.04
0.01
0.02
β1D
0.20
0.22
0.20
0.20
7.04
0.21
0.07
0.07
0.07
0.19
0.08
0.06
1000
β1C
0.02
0.01
0.02
0.02
1.26
0.02
0.01
0.01
0.01
0.01
0.01
0.01
β1D
0.10
0.08
0.10
0.10
6.87
0.10
0.03
0.04
0.03
0.09
0.03
0.03
8
300
β1C
0.30
0.17
0.31
0.31
5.74
0.25
0.04
0.05
0.04
0.23
0.04
0.04
β1D
1.42
0.52
1.63
1.59
39.61
1.51
0.11
0.17
0.11
1.25
0.12
0.13
500
β1C
0.17
0.09
0.18
0.18
5.63
0.14
0.02
0.03
0.02
0.13
0.02
0.02
β1D
0.86
0.29
1.00
1.00
39.68
0.91
0.07
0.10
0.07
0.76
0.08
0.08
1000
β1C
0.08
0.04
0.08
0.08
5.42
0.07
0.01
0.01
0.01
0.06
0.01
0.01
β1D
0.41
0.13
0.41
0.43
38.19
0.41
0.03
0.04
0.03
0.37
0.04
0.04
15
300
β1C
0.60
0.42
0.59
0.55
5.19
0.43
0.10
0.16
0.10
0.46
0.12
0.13
β1D
2.32
0.29
1.83
1.19
49.09
1.36
0.09
0.40
0.09
2.04
0.11
0.17
500
β1C
0.36
0.25
0.37
0.34
5.00
0.27
0.06
0.09
0.06
0.28
0.06
0.07
β1D
1.42
0.16
1.56
1.13
49.09
0.97
0.05
0.20
0.06
1.24
0.06
0.09
1000
β1C
0.18
0.12
0.18
0.17
4.80
0.13
0.03
0.03
0.03
0.14
0.03
0.03
β1D
0.68
0.07
0.69
0.96
48.92
0.50
0.03
0.05
0.03
0.60
0.03
0.04
MSE comparison between SIMP and other 11 competitors for estimating β1C and β1D over 500 repetitions under the data generating mechanism (M3). The lowest MSE under each combination of r and n is in bold face.
Table 4 shows the performance under (M2). Since ‖Ω‖≪‖Ω0‖, the partial predictor envelope component cannot provide too many efficiency gains. Therefore, as expected, the best estimation performance comes from envelope models that contain the (partial) response envelope structure, including FY-env, FS-env, FPY-env, PY-env and SIMP. Although the best performance is mostly achieved by FY-env and FPY-env, the difference between these five methods is small. Since the partial response envelope structure is imposed on both β1C and β1D in (3.7), large efficiency gains could be observed in estimating both of these two parameters for these five methods, as expected. Among these five methods, the slightly worse performance of SIMP under some cases is largely due to the model selection. See the performance with the true envelope dimensions known in Table C.3 in Appendix C.
The estimation performance under (M3) is shown in Table 5. Without surprise, the best performance is obtained by envelope models containing the partial predictor envelope structure, i.e. PX-env and SIMP. Similarly, their estimation performances for β1C are close, and the tiny difference between them mainly comes from the model selection. It is worthy to notice that when r=8 or 15, some additional gains on estimating β1D could be offered by SIMP compared with PX-env, since the partial response envelope structure is assumed in SIMP for β1D, although the signal may not be strong (‖Φ‖≫‖Φ0‖), however, the synergetic effect of SIMP that is mentioned above could possibly enhance this advantage.
MSE comparison between SIMP and other 11 competitors for estimating β1C and β1D over 500 repetitions under the data generating mechanism (M4). The lowest MSE under each combination of r and n among all methods and within envelope methods are in bold face and blue respectively.
Results on (M4): Model Mis-Specification
To test the robustness of our method, an additional simulation under the model mis-specification scenario (M4) is studied. From Table 6, RRR estimator almost always produces the lowest MSE when r=3 or 8. Although under the mis-specified model, envelope methods still perform quite well. Within envelope methods, SIMP obtains a relatively better performance for both β1C and β1D when r is 8, and FX-env is better for β1D when r is 3 or 15. Note that when r=15, several envelope methods could even outperform the RRR estimator. When we are more cautious for the results displayed here, we will find that the performance of all envelope methods especially SIMP is largely contaminated by the “bad” envelope dimensions that are selected under this model mis-specification scenario, since the likelihoods of them might significantly depart from the truth and hence are misleading in model selection. See Table C.5 in Appendix C for the estimation results if we rule out this factor by using the “optimal” tuning parameter or envelope dimension for each method, if we pretend to know the true values of β1C and β1D. Under such scenario, the envelope methods almost always dominate other methods (even RRR), and SIMP is almost always the best one among envelope methods except when r=3. And when r=3, the performance of SIMP is very close to the best one. Therefore, the robustness of our method and other envelope methods to a general low rank model RRR is observed in our simulation.
Besides, in Appendix C, we also display additional numerical results for the estimation performance of FPY-env, PX-env, PY-env and SIMP with the true envelope dimensions known under (M1)–(M3) or RRR with the true rank known under (M4), and the results for p2=4 or pD=4 with other settings same as those in (M1).
Real Imaging Genetics Application on ADNI1Background
Data used in the preparation of this article were obtained from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). Led by Principal Investigator Michael W. Weiner, MD, ADNI was launched in 2003 as a public-private partnership, and lasts for 4 stages (ADNI1, ADNI-GO, ADNI2 and ADNI3) till now. The primary goal of ADNI has been to test whether serial magnetic resonance imaging (MRI), positron emission tomography (PET), genetic or other biological markers, clinical and neuropsychological assessments can be combined to measure the progression of mild cognitive impairment (MCI) and early Alzheimer’s disease (AD). This study has strict enrollment standards, follow-up and data checking protocols. All raw or preprocessed data is available through the Image and Data Archive (IDA) at https://ida.loni.usc.edu/, upon approval of the application. In this real data application, we want to study the genetic effects, including the effect of a well-known AD gene Apolipoprotein E (APOE) ϵ4, on the volumes of some cortical or sub-cortical structures in brain, while also accounting for the effects of some other prognostic factors by applying SIMP, and compare with several other methods in terms of the prediction performance. Following previous practices on analyzing ADNI data [68, 46], we focus on Caucasian participants in the ADNI1 phase, to reduce the population stratification.
Data Pre-processing Procedure
Imaging responses, i.e. Y in (3.2), we adopted are the volume measurements on 12 ROIs that are listed in [42] (i.e., the Amygdala, the Cerebral cortex, the Cerebral white matter, the Hippocampus, the Inferior lateral ventricle and the Lateral ventricle on two hemispheres), which are obtained from 1.5T MRI scan at the screening visit of ADNI1. All imaging data is downloaded from IDA, and has been pre-processed by FreeSurfer [21] already. Among 845 participants who took MRI scans, we only retain 593 participants who passed the quality control (QC). The log transformation is applied to reduce skewness.
MSPE from 5-fold CV for various envelope and non envelope methods. Envelope dimensions, number of components (ncomp) or rank selected by 5-fold CV are mentioned inside parenthesis.
Non envelope methods
Non-partial envelope methods
Partial envelope methods
FOLS
14.02
FX-env (dX=1)
8.34
PX-env (dX=1)
7.94
RRR (rank=1)
10.87
FY-env (dY=1)
9.77
FPY-env (dY=1)
9.82
PLSR (ncomp=1)
8.20
FS-env (dX=7, dY=1)
8.06
SIMP(dX=1, dY=3)
7.88
620,901 SNPs on all autosomes and sex chromosomes were genotyped for 757 participants in ADNI1, who are not contained in but are overlapped with the 845 participants who have MRI scans. Among 11,632 SNPs meta-analyzed by International Genomics of Alzheimer’s Project (IGAP) [36], 852 autosome SNPs are selected as our candidate SNPs for model fitting by checking two conditions: (1) Genotyped by ADNI1; (2) P-values from IGAP are smaller than 0.01. Developed by Shaun Purcell and Christopher Chang, PLINK [7] is an open source software for performing QC and some routine analyses in GWAS in an efficient manner, and is available at https://www.cog-genomics.org/plink/1.9/. Following [68], we perform the following two lines of QC on the selected SNP data by PLINK. The first line of QC includes (1) call rate check per subject and per SNP marker, (2) gender check, (3) sibling pair identification, (4) Hardy-Weinberg equilibrium test, and (5) population stratification check. The second line of QC removes SNPs with (1) more than 5% missing values, (2) minor allele frequency smaller than 10%, and (3) Hardy-Weinberg equilibrium p-value <10−6. This leaves us with 732 SNPs that passed QC. A free genotype imputation and haplotype phasing program IMPUTE2 ([27], [28]) is used to impute SNPs with missingness. We choose HapMap 3 as our reference panel for the imputation, since it shares the same genome build b36 of National Center for Biotechnology Information (NCBI) with ADNI1. We delete 39 SNPs that are not shown in the reference panel, which leads to 693 SNPs finally. To reduce dimensionality, the principal component analysis (PCA) is applied on these 693 SNPs. By thresholding the eigenvalues of its sample covariance matrix by 1, the first 186 principal components (PCs), which explains 87.19% of the total variation of these 693 SNP predictors, are selected as our X1C, the continuous part of our predictors of main interest in (3.2).
APOE ϵ4, taking values in {0,1,2}, is included in X1D as the only discrete (quantitative) genetic predictor in (3.2). We have incorporated six other important prognostic factors in X2 in (3.2), including gender, marital status, handness, age, years of education and intracranial volume (ICV). The intersection of participants from all types of abovementioned datasets leads to 498 samples in our analysis finally. All predictors and responses are standardized before the model fitting.
Prediction Performance
We first compare the prediction performance of SIMP with other two partial envelope methods (PX-env, FPY-env), three envelope models that cannot account for the partial structures (FX-env, FY-env, FS-env; we call them non-partial envelope methods), and three popular multivariate linear regression methods that are not based on the envelope models (FOLS, RRR and PLSR; we call them non envelope methods). Similar with Section 7.2, for all methods except the partial envelope methods (PX-env, FPY-env and SIMP), Y is adjusted as the residuals from the regression on X2 first, and the model fitting will be done by the genetic markers X1 and the adjusted responses. For the best prediction performance, all associated envelope dimensions for envelope methods, or the number of components for PLSR, or the rank for RRR are selected by minimizing MSPE from 5-fold CV. The MCMC algorithm for SIMP and PX-env is run for 20,000 iterations with first 50% iterations as burn-in, and their specifications for the hyperparameters are the same as those in Section 7.
Table 7 reveals that the best prediction performance is achieved by SIMP. All of the envelope methods in the last two columns significantly outperform FOLS and RRR. FS-env and SIMP achieve the best prediction performance within the non-partial envelope methods and partial envelope methods we considered respectively. This illustrates the advantages of combining the predictor and response envelopes in both non-partial and partial envelope contexts. Comparing partial envelope methods with their non-partial counterparts, PX-env and SIMP both obtain lower MSPEs than FX-env and FS-env respectively. The finding that FPY-env cannot outperform FY-env seems surprising, but is consistent with the previous theoretical result that if the envelope spaces of FPY-env and FY-env are equal, then the asymptotic variance of the estimator of β1 from FPY-env could not be smaller than that from FY-env (see the proposition 2 of [57]), and we have verified the estimated envelope spaces by FPY-env and FY-env are very close in this analysis. It is also interesting to note that FX-env shares similar MSPE with PLSR. Their connection has been discussed in [12].
Upper left: Indicator of significance of each (SNP, IP) pair or APOE ϵ4 with each IP. APOE ϵ4 corresponds to the last vertical line in red; Upper right: Indicator of significance of each (prognostic factor, IP) pair. Red signifies significance in the figures of the upper row. Bottom: Heatmap of the estimated regression coefficients matrix βˆSNPT between 693 SNPs and 12 IPs.
Posterior Estimation and Selection Performance
Both AIC-MCMC and BIC-MCMC select dX=87 and dY=1 for SIMP. Using dX=87 and dY=1, SIMP is fitted with the number of iterations, burn-in proportion and specification of the hyperparameters same as those in Section 7 and Section 8.3. It costs 14.82 hours for a run of 20, 000 iterations under these envelope dimensions, for a Intel Xeon E5-2680 v3 processor (2.50 GHz CPU).
Aiming at interpreting our results with some scientifically meaningful findings and due to the large variation that is explained by the PCs that we choose, we intend to approximate the posterior samples of βSNP∈R693×12, the regression coefficients between the standardized version of the original 693 SNP predictors (i.e., the 693 SNP predictors before PCA) and 12 IPs. For s=1,2,…,10,000, the s-th retained posterior sample β1C(s)∈R186×12 is mapped to βSNP(s), the s-th approximate posterior sample of βSNP by
βSNP(s)=(diag{SSNP})1/2·LD·(diag{S1C})−1/2·β1C(s),
where diag{SSNP}∈S+693×693 and diag{S1C}∈S+186×186 are diagonal matrices with diagonal elements being sample variances of the original 693 SNPs and 186 PCs respectively, and LD∈R693×186 is the loading matrix from PCA.
Out of 37 significant SNPs from SIMP (besides APOE ϵ4), 17 SNPs are reported here for the previously reported associations of them (or their RefSeq gene or their closest RefSeq gene with names indicated in Columns 2 and 5) with AD in the literature. In the literature, one previous study on the reported association is selected for each of these SNPs and displayed in Columns 3 and 6.
SNP
Gene
SNP
Gene
rs11685593
BIN1
[2]
rs10894473
NTM
[66]
rs7561528
BIN1
[24]
rs11064498
C1S
[61]
rs11706690
CHL1
[50]
rs757402
OAS2
[3]
rs10513391
P2RY14
[53]
rs2274736
PTPN21
[62]
rs2439538
TBC1D7
[17]
rs10498633
SLC24A4/RIN3
[60]
rs9381563
CD2AP
[30]
rs2554389
ADAMTSL3
[55]
rs2280231
NDUFS3
[49]
rs4265771
ADAMTSL3
[55]
rs7120548
MTCH2
[31]
rs17809911
CCDC102B
[41]
rs10501927
CNTN5
[24]
The corresponding (1−5.95×10−6)×100% two-sided credible intervals for βSNP, β1D and β2 are calculated, where the confidence level is determined by 0.95 with the Bonferroni correction for all 8,400 regression coefficients. Significance of any SNP, APOE ϵ4 or prognostic factor with any one of 12 IPs is determined by the exclusion of zero for the associated credible interval.
The upper-left panel of Figure 3 shows the significance of a few SNP and APOE ϵ4 with all 12 IPs. Besides the well-known AD gene APOE ϵ4, we identify another 37 significant SNPs that are possibly related to AD (see Appendix H.3 for the full list) under the control of the Family-wise error rate (FWER, the probability of reporting at least one false positives) at 0.05. Among them, 17 SNPs are reported in Table 8 for their appearance in the previous studies. Meanwhile, under the same control of FWER, the key role of the age for AD and the impact of ICV on the volumes of brain structures are confirmed at the same time, by the upper-right panel of Figure 3. The significance of SNPs seems to be consistent over all IPs but this is not observed for prognostic factors. This is expected since the rows of βˆ1CT (and also βˆSNPT, the posterior mean estimator of βSNPT) depend on the estimated partial response envelope only, which happens to contribute much weaker effects than the estimated partial predictor envelope in this dataset. Such pattern is not observed in βˆ2T, since we have not imposed any envelope structure on β2.
The heatmap of the estimated regression coefficients matrix βˆSNPT is shown in the bottom panel of Figure 3. The elements of βˆ1D are much larger in scale than those in βˆSNPT, hence are not displayed in this figure (The elements of βˆ1D are listed here instead. Amygdala.L: −0.19, Amygdala.R: −0.19, Cerebral cortex.L: −0.12, Cerebral cortex.R: −0.13, Cerebral white matter.L: −0.11, Cerebral white matter.R: −0.11, Hippocampus.L: −0.22, Hippocampus.R: −0.22, Inferior lateral ventricle.L: 0.18, Inferior lateral ventricle.R: 0.19, Lateral ventricle.L: 0.14, Lateral ventricle.R: 0.15). It is interesting that for each genetic predictor, the estimated coefficients related to each pair of brain measures on two hemispheres are close. This is reasonable since the bilateral correlations within all pairs of brain measures over two hemispheres are strong (see Appendix H.4 for the numerical evidence and [52] for another discussion on this bilateral correlation). The ROIs that have relatively large estimated effects with some SNPs are the Inferior lateral ventricle, the Hippocampus and the Amygdala (in both hemispheres). The Hippocampus and the Amygdala are related to memory and motor behavior respectively, and the importance of these three brain structures for AD have already been verified in the past studies [20, 47, 43].
Shrinkage Estimation by Incorporating the Prior Information of Weak Imaging Genetics Relationship
The past studies on AD have verified the weak effects of the SNP predictors on predicting AD outcomes, either the disease status [36] or the imaging phenotypes ([68, 46] and etc). This pattern of weak effects has also been verified by us from the bottom panel of Figure 3 in Section 8.4. Therefore, it is reasonable to incorporate the prior information of “weak signals in β1C” (and hence “weak signals in βSNP”) into our analysis, and it should also improve the estimation performance from the standpoint of the shrinkage effect offered for this high dimensional problem, and help to illustrate the advantage of our Bayesian envelope method on prior information incorporation, compared with the frequentist envelope methods. The envelope dimensions are fixed at dX=87 and dY=1. We adjust the (all equal) diagonal elements of the hyperparameter E (we set to be a diagonal matrix) in the prior of ηC from 10−6 to 106 while keeping the specification for other hyperparameters same as previous sections. This adjustment strengthens our prior belief of ηC (and hence β1C) to be a zero matrix gradually, since the prior mean of ηC is fixed at the zero matrix by assigning W to be the zero matrix, and the prior covariance matrices of rows of ηC are proportional to E−1. From the frequentist perspective, the Bayesian estimator (strictly speaking, the MAP estimator, but here we use the posterior mean estimator for simplicity) with increasingly stronger prior belief in the weak signals of ηC corresponds to an estimator with gradually stronger shrinkage to the zero matrix.
The left panel of Figure 4 shows improved prediction performance as the diagonal elements of E increase from 10−6 to 103, and the MSPE deteriorates slightly as the prior belief increases further beyond 103, possibly due to the bias that is caused by the over-strong prior information of weak signals, or the over-penalization from the frequentist perspective. The middle and right panels of Figure 4 display the regularization paths for the row sums and the column sums of βˆSNP respectively. The more significant shrinkage effects could be observed for the stronger prior belief of the weak signals in β1C (and hence in βSNP).
MSPE from the 5-fold CV (left), the row sums of βˆSNP (middle; each line corresponds to one SNP) and the column sums of βˆSNP (right; each line corresponds to one IP) with respect to the diagonal elements of E (in the log10 scale).
Discussion
In this paper, we propose a new unified Bayesian envelope model by integrating the partial predictor envelope and the partial response envelope under a convenient Bayesian framework. The proposed model degenerates to several well established envelope models in the literature under specific conditions. It addresses the limitations mentioned for the simultaneous and the partial response envelopes. Specifically, our method improves the efficiencies for estimating β1 and predicting Y in (3.2), and has no restrictions on X1 to be continuous or Normal, by a subtle construction of separating the discrete part from the whole predictors of interest. Compared with the frequentist envelope approaches, our method is more flexible in incorporating prior information and quantifying uncertainty through posterior distribution of parameters. Overall, our method is inclusive and could be regarded as a building block for the future theoretical research of the envelope model, and an ideal solution for practitioners who seek dimension reduction and want to apply envelope methods to their application problems, including almost any problems that could be formulated by the multivariate linear regression with predictors either continuous or discrete or a mix of them and even containing nuisance covariates, but are worried about which specific envelope model to use. Meanwhile, we are the first to investigate the performance of several popular dimension selection methods together, among the Bayesian envelope literature, to the best of our knowledge.
To be less conservative and make more meaningful discoveries, a well-designed multiple comparisons procedure to control the False discovery rate, rather than FWER, could be considered in our real data application. A natural extension of our work is to generalize SIMP to the generalized linear model setting. However, for the multivariate probit model as pointed in [9], the identification issue exists for the error covariance matrix, which is suggested to be restricted to the correlation matrix. Hence, the extension of SIMP to the generalized linear model setting might require some special strategies, for example the parameter expansion for data augmentation technique ([40], [59]) under a marginally uniform prior for the error correlation matrix [1].
Metropolis-within-Gibbs MCMC Algorithm for <inline-formula id="j_nejsds23_ineq_958"><alternatives><mml:math>
<mml:mi mathvariant="normal">SIMP</mml:mi></mml:math><tex-math><![CDATA[$\mathrm{SIMP}$]]></tex-math></alternatives></inline-formula>
In this section, we show how to generate a MCMC chain of length S on Θ from SIMP, for a pre-specified pair of dimensions (dX,dY)∈{1,…,pC−1}×{1,…,r−1}. We could arbitrarily choose the initial value Θ(0), but a warm start initial estimator provided in Appendix B is recommended for faster convergence. For each s=1,…,S, we iterate the following steps to update Θ.
Generate β2(s) from MNp2,r(M˜−1Z˜(s−1),M˜−1,ΣY∣X(s−1)), where
M˜=X˜2TX˜2+M,Z˜(s−1)=X˜2T(Y˜(s−1)−X˜1C(s−1)L(A(s−1))(ηC(s−1))TR(B(s−1))T−X˜1D(ηD(s−1))TR(B(s−1))T)+Z.
Generate γ(s) from MNpD,pC(Λ˜−1F˜(s−1),Λ˜−1,ΣC∣D(s−1)), where Λ˜=X˜1DTX˜1D+Λ, and F˜(s−1)=X˜1DTX˜1C(s−1)+F.
Generate (μ1C(s),μY(s))T from NpC+r((X‾1CY‾),Δ(s−1)/n) where Δ(s−1)=ΣC∣D(s−1)L(A(s−1))Ω(s−1)(ηC(s−1))TR(B(s−1))TR(B(s−1))ηC(s−1)Ω(s−1)L(A(s−1))TΣY∣X(s−1)+R(B(s−1))ηC(s−1)Ω(s−1)(ηC(s−1))TR(B(s−1))T.
Generate ηC(s) from MNdY,dX(W˜(s−1)(E˜(s−1))−1,Φ(s−1),(E˜(s−1))−1), where
E˜(s−1)=L(A(s−1))T(X˜1C(s−1))TX˜1C(s−1)L(A(s−1))+E,W˜(s−1)=((Y˜(s−1)−X˜2β2(s−1))R(B(s−1))−X˜1D(ηD(s−1))T)T(X˜1C(s−1))L(A(s−1))+W.
Generate ηD(s) from MNdY,pD(T˜(s−1)(Q˜)−1,Φ(s−1),(Q˜)−1), where Q˜=(X˜1D)TX˜1D+Q and T˜(s−1)=((Y˜(s−1)−X˜2β2(s−1))R(B(s−1))−X˜1C(s−1)L(A(s−1))(ηC(s−1))T)TX˜1D+T.
Generate Ω(s) from IWdX(Ψ˜X(s−1),w˜X), where
Ψ˜X(s−1)=L(A(s−1))T(X˜1C(s−1)−X˜1Dγ(s−1))T(X˜1C(s−1)−X˜1Dγ(s−1))L(A(s−1))+L(A(s−1))T(γ(s−1)−Λ−1F)TΛ(γ(s−1)−Λ−1F)L(A(s−1))+ΨX
and w˜X=n+pD+wX.
Generate Ω0(s) from IWpC−dX(Ψ˜X0(s−1),w˜X0), where
Ψ˜X0(s−1)=L0(A(s−1))T(X˜1C(s−1)−X˜1Dγ(s−1))T(X˜1C(s−1)−X˜1Dγ(s−1))L0(A(s−1))+L0(A(s−1))T(γ(s−1)−Λ−1F)TΛ(γ(s−1)−Λ−1F)L0(A(s−1))+ΨX0
and w˜X0=n+pD+wX0.
Generate Φ(s) from IWdY(Ψ˜Y(s−1),w˜Y), where
Ψ˜Y(s−1)=((Y˜(s−1)−X˜2β2(s−1))R(B(s−1))−X˜1C(s−1)L(A(s−1))(ηC(s−1))T−X˜1D(ηD(s−1))T)T((Y˜(s−1)−X˜2β2(s−1))R(B(s−1))−X˜1C(s−1)L(A(s−1))(ηC(s−1))T−X˜1D(ηD(s−1))T)+R(B(s−1))T(β2(s−1)−M−1Z)TM(β2(s−1)−M−1Z)R(B(s−1))+(ηC(s−1)−WE−1)E(ηC(s−1)−WE−1)T+(ηD(s−1)−TQ−1)Q(ηD(s−1)−TQ−1)T+ΨY
and w˜Y=n+dX+pD+p2+wY.
Generate Φ0(s) from IWr−dY(Ψ˜Y0(s−1),w˜Y0), where
Ψ˜Y0(s−1)=R0(B(s−1))T(Y˜(s−1)−X˜2β2(s−1))T(Y˜(s−1)−X˜2β2(s−1))R0(B(s−1))+R0(B(s−1))T(β2(s−1)−M−1Z)TM(β2(s−1)−M−1Z)R0(B(s−1))+ΨY0
and w˜Y0=n+p2+wY0.
Let h(s)(A) be the log full conditional density of A at the s-th iteration, i.e.,
h(s)(A)=const−12tr{(X˜1C(s−1)−X˜1Dγ(s−1))(L(A)Ω(s−1)L(A)T+L0(A)Ω0(s−1)L0(A)T)−1(X˜1C(s−1)−X˜1Dγ(s−1))T+(Y˜(s−1)−X˜1C(s−1)L(A)(ηC(s−1))TR(B(s−1))T−X˜1D(ηD(s−1))TR(B(s−1))T−X˜2β2(s−1))(ΣY∣X(s−1))−1(Y˜(s−1)−X˜1C(s−1)L(A)(ηC(s−1))TR(B(s−1))T−X˜1D(ηD(s−1))TR(B(s−1))T−X˜2β2(s−1))T+(L(A)Ω(s−1)L(A)T+L0(A)Ω0(s−1)L0(A)T)−1(γ(s−1)−Λ−1F)TΛ(γ(s−1)−Λ−1F)+ΣA−1(A−A0)TKA−1(A−A0)}.
We want to generate a Markov chain realization for A from the stationary density proportional to exp(h(s)(A)). We update A(s−1) to A(s) columnwisely via Metropolis steps as follows. Let Aj(s−1)∈RpC−dX denote the jth column of A(s−1), j=1,…,dX. Let τA>0 be a given tuning parameter and let {i1,…,idX} be a random permutation of {1,…,dX}. For each k=ij, j=1,…,dX,
Generate Ak∗ from NpC−dX(Ak(s−1),τA2IpC−dX).
Calculate the ratio
r=exp{h(A∗)}/exp{h(A(s−1))},
where A∗ is the resulting matrix after Ak(s−1) is replaced with Ak∗ in A(s−1).
Generate a binary indicator with success probability r. If a success is achieved, update Ak(s−1) to Ak∗ as the k-th column of A(s). Otherwise, retain Ak(s−1) as the k-th column of A(s).
Update A(s−1) to A(s) after all columns are investigated.
Let h(s)(B) be the log full conditional density of B at the s-th iteration, i.e.,
h(s)(B)=const−12tr{(Y˜(s−1)−X˜1C(s−1)L(A(s−1))(ηC(s−1))TR(B)T−X˜1D(ηD(s−1))TR(B)T−X˜2β2(s−1))(R(B)Φ(s−1)R(B)T+R0(B)Φ0(s−1)R0(B)T)−1(Y˜(s−1)−X˜1C(s−1)L(A(s−1))(ηC(s−1))TR(B)T−X˜1D(ηD(s−1))TR(B)T−X˜2β2(s−1))T+(R(B)Φ(s−1)R(B)T+R0(B)Φ0(s−1)R0(B)T)−1(β2(s−1)−M−1Z)TM(β2(s−1)−M−1Z)+ΣB−1(B−B0)TKB−1(B−B0)}.
The procedure to generate a Markov chain realization for B from the stationary density proportional to exp(h(s)(B)) is similar as that for A in Step 11.
Note that whendX=0(or,dX=pC), steps 5, 7 and 11 (or, steps 8 and 11) should be skipped,L(s)does not exist,L0(s)=IpC(or,L(s)=IpCandL0(s)does not exist) always and any terms involving eitherL(s)orΩ(s)orηC(s)(or, eitherL0(s)orΩ0(s)) orA(s)are omitted. Similarly, whendY=0(or,dY=r), steps 5, 6, 9 and 12 (or, steps 10 and 12) should be skipped,R(s)does not exist,R0(s)=Ir(or,R(s)=IrandR0(s)does not exist) always and any terms involving eitherR(s)orΦ(s)orηC(s)orηD(s)(or, eitherR0(s)orΦ0(s)) orB(s)are omitted.
In our actual implementation of this Metropolis-within-Gibbs algorithm, once any parameter is updated, it will be immediately used for the updating of other parameters in the same iteration. Even when we are updatingAandB, the previously updated columns will be used immediately to update the rest of columns.
The updating of each column ofAandBin steps 11 and 12 is implemented in a random order as illustrated. Instead, it could be in a deterministic order as well. Meanwhile, the order of updating different parameter blocks in an iteration could also be random without affecting the convergence.
MSE comparison for SIMP and other 11 estimators of β1C and β1D over 500 repetitions under the data generating mechanism (M5). The lowest MSE under each combination of r and n is in bold face.
r
n
OLS
RRR
PCR
PLSR
CCA
FX-env
FY-env
FS-env
FPY-env
PX-env
PY-env
SIMP
3
300
β1C
0.16
0.07
0.18
0.17
3.82
0.10
0.06
0.05
0.05
0.02
0.05
0.02
β1D
0.60
0.26
1.57
1.33
14.70
1.16
0.20
0.58
0.17
0.25
0.18
0.10
500
β1C
0.09
0.04
0.10
0.10
3.78
0.05
0.03
0.03
0.03
0.01
0.03
0.01
β1D
0.35
0.15
0.67
0.64
14.47
0.58
0.11
0.22
0.10
0.14
0.10
0.06
1000
β1C
0.05
0.02
0.05
0.05
3.85
0.02
0.01
0.01
0.01
0.01
0.01
0.01
β1D
0.17
0.07
0.17
0.21
14.98
0.19
0.05
0.07
0.05
0.06
0.05
0.03
8
300
β1C
0.69
0.11
0.59
0.52
22.98
0.27
0.07
0.08
0.06
0.04
0.06
0.03
β1D
2.46
0.34
4.29
4.08
84.50
3.31
0.24
0.54
0.18
0.90
0.18
0.12
500
β1C
0.41
0.06
0.42
0.41
22.82
0.19
0.04
0.04
0.03
0.02
0.03
0.02
β1D
1.46
0.19
1.88
1.90
83.58
1.75
0.13
0.21
0.10
0.55
0.10
0.07
1000
β1C
0.20
0.03
0.21
0.20
22.76
0.09
0.02
0.02
0.02
0.01
0.02
0.01
β1D
0.73
0.08
0.84
1.19
83.02
1.02
0.05
0.08
0.05
0.27
0.05
0.03
15
300
β1C
1.44
0.17
0.20
0.16
30.25
0.13
0.09
0.13
0.08
0.06
0.08
0.06
β1D
4.66
0.46
2.25
2.20
100.73
2.21
0.22
0.81
0.17
1.94
0.17
0.23
500
β1C
0.85
0.09
0.17
0.14
30.13
0.11
0.05
0.06
0.04
0.04
0.05
0.03
β1D
2.81
0.22
2.23
2.17
101.26
2.22
0.11
0.26
0.09
1.16
0.09
0.13
1000
β1C
0.42
0.04
0.17
0.13
30.18
0.10
0.02
0.03
0.02
0.02
0.02
0.01
β1D
1.38
0.10
2.18
2.10
101.38
2.11
0.05
0.10
0.05
0.57
0.05
0.04
MSE comparison for SIMP and other 11 estimators of β1C and β1D over 500 repetitions under the data generating mechanism (M6). The lowest MSE under each combination of r and n is in bold face.
r
n
OLS
RRR
PCR
PLSR
CCA
FX-env
FY-env
FS-env
FPY-env
PX-env
PY-env
SIMP
3
300
β1C
0.15
0.06
0.17
0.17
3.13
0.10
0.05
0.07
0.04
0.02
0.05
0.01
β1D
1.18
0.50
1.50
1.43
24.57
1.59
0.36
0.97
0.35
0.44
0.35
0.15
500
β1C
0.09
0.04
0.10
0.10
3.16
0.07
0.03
0.04
0.03
0.01
0.03
0.01
β1D
0.71
0.29
0.84
0.82
24.43
1.06
0.21
0.51
0.20
0.25
0.20
0.09
1000
β1C
0.04
0.02
0.05
0.05
3.14
0.04
0.01
0.02
0.01
0.00
0.01
0.00
β1D
0.34
0.14
0.36
0.40
23.84
0.43
0.10
0.16
0.10
0.13
0.10
0.04
8
300
β1C
0.69
0.09
0.64
0.62
18.75
0.31
0.06
0.08
0.05
0.03
0.05
0.02
β1D
4.61
0.76
3.53
3.37
122.38
2.96
0.36
0.80
0.33
1.87
0.33
0.20
500
β1C
0.41
0.05
0.40
0.36
18.62
0.20
0.03
0.04
0.03
0.02
0.03
0.01
β1D
2.75
0.40
2.47
1.70
122.52
1.85
0.20
0.43
0.20
1.12
0.20
0.11
1000
β1C
0.20
0.03
0.21
0.18
18.43
0.11
0.02
0.02
0.01
0.01
0.01
0.00
β1D
1.36
0.21
1.47
0.87
121.33
1.22
0.10
0.18
0.10
0.55
0.10
0.06
15
300
β1C
1.45
0.16
0.38
0.34
30.92
0.31
0.07
0.12
0.06
0.06
0.06
0.04
β1D
9.46
1.17
7.10
6.68
200.78
6.85
0.42
0.88
0.36
3.88
0.37
0.32
500
β1C
0.86
0.08
0.48
0.37
30.68
0.32
0.04
0.06
0.04
0.03
0.04
0.02
β1D
5.52
0.57
6.74
6.39
198.22
4.80
0.23
0.43
0.21
2.30
0.21
0.17
1000
β1C
0.42
0.03
0.45
0.38
30.49
0.21
0.02
0.03
0.02
0.02
0.02
0.01
β1D
2.79
0.26
3.26
2.20
198.11
2.35
0.11
0.20
0.10
1.14
0.10
0.07
Initial Estimator for MCMC Algorithm of <inline-formula id="j_nejsds23_ineq_1090"><alternatives><mml:math>
<mml:mi mathvariant="normal">SIMP</mml:mi></mml:math><tex-math><![CDATA[$\mathrm{SIMP}$]]></tex-math></alternatives></inline-formula>
In this section, we propose a warm start initial estimator for all parameters in SIMP, including μ1C, μY, β2, γ, ηC, ηD, A, B, Ω, Ω0, Φ and Φ0 if all existed (i.e. 0<dX<pC, 0<dY<r and r, pC, pD, p2>0). It will be used as the initial value Θ(0) for our implementation of the MCMC algorithm for SIMP, for faster convergence. Given the data matrices X1C, X1D, X2 and Y, and any valid envelope dimensions dX, dY as input, we implement the following steps.
Estimate μY(0)=YT1n/n and μ1C(0)=X1CT1n/n.
Center X1C, X1D, X2 and Y to be X˜1C=X1C−1n(μ1C(0))T, X˜1D=X1D−1nX‾1DT, X˜2=X2−1nX‾2T and Y˜=Y−1n(μY(0))T.
Perform the frequentist ordinary least squares for X˜1C on X˜1D, and get γ(0) and ΣC∣D(0) as the estimate for the regression coefficients matrix and the residual covariance matrix, respectively.
Perform the frequentist partial response envelope model with dimension dY (by the function penv() in the R library Renvlp) on
regressingY˜onX˜1C,X˜1DandX˜2,
where we treat X˜1C and X˜1D together as our predictors of interest in this partial response envelope model. From this fitting, β1C(0), β1D(0), β2(0), R(B(0)), R0(B(0)), Φ(0) and Φ0(0) could be directly obtained. Compute B(0) from R(B(0)) by the trick of re-parameterization, as illustrated in Section 3.2. Estimate ηD by ηD(0)=(β1D(0)R(B(0)))T.
Perform the frequentist predictor envelope model with dimension dX (by the function xenv() in the R library Renvlp) on the regression from Y˜∗ on X˜1C∗, where Y˜∗ and X˜1C∗ are the residual matrices from the following two frequentist ordinary least squares regressions respectively:
RegressY˜onX˜1DandX˜2,RegressX˜1ConX˜1DandX˜2.
Estimates L(A(0)), L0(A(0)) could be immediately obtained from this regression. Then, Ω(0) and Ω0(0) could be calculated as Ω(0)=L(A(0))TΣC∣D(0)L(A(0)) and Ω0(0)=L0(A(0))TΣC∣D(0)L0(A(0)), and A(0) could be computed from L(A(0)) by the idea in Section 3.2 similarly. Estimate ηC(0)=R(B(0))T(β1C(0))TL(A(0)).
Note when dX=0 (or, dX=pC), Ω, A and ηC (or, Ω0 and A) don’t exist. When dY=0 (or, dY=r), Φ, B, ηC and ηD (or, Φ0 and B) don’t exist. Therefore, it is not necessary to initialize them in the corresponding cases.
Additional Simulation Results for Section <xref rid="j_nejsds23_s_017">7.2</xref>
Tables C.1 and C.2 show the parameter estimation performance for β1C and β1D under exactly the same setting for (M1) in Section 7.2, except now p2=4 with μ2=(2,−1,−3,5) or pD=4 respectively. We call these two data generating mechanisms as (M5) and (M6). As Table C.1 reveals, MSE for both β1C and β1D from our method and our competitors remain roughly the same or become little inflated as p2 increases from 2 to 4, except few cases. Table C.2 suggests the MSE for β1D grow linearly with respect to pD, and the MSE for β1C remain almost the same for our method and most of our competitors, except few cases. Therefore, SIMP is observed to show more advantageous numerical performance with increased pD in our limited simulation studies.
Table C.3 summarizes the performance of FPY-env, PX-env, PY-env and SIMP with the true envelope dimensions known under (M1)–(M3), (M5) and (M6), or the performance of RRR with the true rank known under (M4). With the true envelope dimensions known, SIMP almost always obtains the best estimation performance under all scenarios except the model mis-specification case (M4).
MSE comparison between FPY-env, PX-env, PY-env and SIMP for estimating β1C and β1D over 500 repetitions, with the true envelope dimensions known under data generating mechanisms (M1)–(M3), (M5) and (M6), and the MSE of the RRR estimator with the true rank known under data generating mechanism (M4). “*” in the names of methods means that the true envelope dimension/rank is used. The lowest MSE for each combination of r and n under each true data generating mechanism is highlighted in bold face.
r
3
8
15
n
300
500
1000
300
500
1000
300
500
1000
β1C
β1D
β1C
β1D
β1C
β1D
β1C
β1D
β1C
β1D
β1C
β1D
β1C
β1D
β1C
β1D
β1C
β1D
M1
FPY-env*
0.05
0.16
0.03
0.09
0.01
0.05
0.05
0.18
0.03
0.11
0.01
0.05
0.08
0.15
0.04
0.09
0.02
0.04
PX-env*
0.01
0.22
0.01
0.14
0.00
0.07
0.03
0.92
0.02
0.57
0.01
0.27
0.06
1.87
0.04
1.12
0.02
0.55
PY-env*
0.05
0.16
0.03
0.09
0.01
0.05
0.05
0.18
0.03
0.11
0.01
0.05
0.08
0.15
0.04
0.09
0.02
0.04
SIMP*
0.01
0.07
0.01
0.04
0.00
0.02
0.02
0.09
0.01
0.05
0.00
0.03
0.04
0.09
0.03
0.05
0.01
0.03
M2
FPY-env*
0.02
0.10
0.01
0.06
0.01
0.03
0.04
0.11
0.02
0.07
0.01
0.03
0.10
0.09
0.06
0.05
0.03
0.03
PX-env*
0.05
0.30
0.03
0.19
0.01
0.09
0.23
1.25
0.13
0.76
0.06
0.37
0.46
2.04
0.28
1.24
0.14
0.59
PY-env*
0.02
0.10
0.01
0.06
0.01
0.03
0.04
0.11
0.02
0.07
0.01
0.03
0.11
0.09
0.06
0.05
0.03
0.03
SIMP*
0.02
0.09
0.01
0.06
0.01
0.03
0.04
0.11
0.02
0.06
0.01
0.03
0.11
0.09
0.06
0.05
0.03
0.03
M3
FPY-env*
0.23
0.81
0.13
0.46
0.06
0.23
0.22
0.80
0.13
0.49
0.06
0.23
0.23
0.70
0.13
0.41
0.07
0.21
PX-env*
0.02
0.33
0.01
0.20
0.00
0.10
0.02
0.47
0.01
0.28
0.01
0.14
0.03
0.66
0.02
0.40
0.01
0.20
PY-env*
0.23
0.81
0.13
0.46
0.06
0.23
0.22
0.80
0.13
0.49
0.06
0.23
0.23
0.69
0.13
0.41
0.07
0.21
SIMP*
0.01
0.30
0.01
0.19
0.00
0.09
0.02
0.30
0.01
0.18
0.01
0.09
0.03
0.30
0.02
0.19
0.01
0.09
M4
RRR*
0.39
0.05
0.23
0.03
0.11
0.01
0.32
0.04
0.17
0.02
0.08
0.01
9.04
1.00
3.75
0.40
1.56
0.16
M5
FPY-env*
0.05
0.17
0.03
0.10
0.01
0.05
0.06
0.18
0.03
0.10
0.02
0.05
0.07
0.16
0.04
0.09
0.02
0.05
PX-env*
0.01
0.22
0.01
0.13
0.00
0.06
0.03
0.90
0.02
0.55
0.01
0.27
0.06
1.94
0.03
1.16
0.02
0.57
PY-env*
0.05
0.17
0.03
0.10
0.01
0.05
0.06
0.17
0.03
0.10
0.02
0.05
0.07
0.16
0.04
0.09
0.02
0.05
SIMP*
0.01
0.08
0.01
0.05
0.00
0.02
0.02
0.09
0.01
0.05
0.01
0.03
0.04
0.09
0.02
0.06
0.01
0.03
M6
FPY-env*
0.04
0.35
0.03
0.20
0.01
0.10
0.05
0.33
0.03
0.20
0.01
0.10
0.06
0.36
0.04
0.21
0.02
0.10
PX-env*
0.01
0.40
0.00
0.24
0.00
0.12
0.03
1.86
0.02
1.13
0.01
0.55
0.06
3.87
0.03
2.31
0.02
1.14
PY-env*
0.04
0.35
0.03
0.20
0.01
0.10
0.05
0.33
0.03
0.19
0.01
0.10
0.06
0.37
0.04
0.21
0.02
0.10
SIMP*
0.00
0.13
0.00
0.08
0.00
0.04
0.01
0.18
0.01
0.11
0.00
0.05
0.03
0.23
0.02
0.13
0.01
0.06
Table C.4 shows the performance of all three Bayesian envelope estimators (PX-env, PY-env and SIMP) of β1C and β1D under the data generating mechanism (M1), except now the posterior median rather than the posterior mean is chosen as the point estimator. Comparing the results in this table and Table 3, we empirically find no observable difference in performance between the posterior mean and the posterior median as the point estimator under (M1).
MSE comparison for three Bayesian envelope estimators (PX-env, PY-env and SIMP) on estimating β1C and β1D over 500 repetitions under the data generating mechanism (M1), except the posterior median rather than the posterior mean is chosen as the point estimator. The lowest MSE under each combination of r and n is in bold face. The tiny differences between the results in this table and those in Table 3 are colored in red.
MSE comparison between SIMP and other 11 competitors for estimating β1C and β1D over 500 repetitions under the data generating mechanism (M4). Different from Table 6, this table displays the estimation performance of 12 methods, if we pretend that the true values of the parameters β1C and β1D are known, and the “optimal” tuning parameters or the envelope dimensions are selected by minimizing the ℓ2 distances between the estimates and the true values of β1C and β1D separately, instead of applying the model selection methods that are introduced in Section 7.2.1. The lowest MSE under each combination of r and n among all methods and within envelope methods are in bold face and blue respectively.
Theoretical Proofs
The integrability of the posterior density could be immediately obtained following the proof of Theorem 2. □
To prove the Harris ergodicity of the Markov chain generated by our Metropolis-within-Gibbs algorithm, we need to prove 3 properties, namely, (a) ϕ-irreducibility with respect to some measure ϕ, (b) Aperiodicity, and (c) Harris recurrence. Firstly, we prove ϕ-irreducibility which requires that the Markov chain could reach any set with positive measure from any state in the state space. The rigorous definition is given in Definition 3. (<italic>ϕ</italic>-irreducibility).
For a Markov chain(Xn)n=1∞, suppose that the state space isXwith σ-algebraF, and the n-step transition kernel isκn:X×F→[0,1]. Given a nonzero measure ϕ, the Markov chain is called ϕ-irreducible, if for every statex∈Xand every subsetA∈Fwithϕ(A)>0, there exists an integer n such thatκn(x,A)>0.
To prove the ϕ-irreducibility of the Markov chain generated by the Metropolis-within-Gibbs algorithm that we developed, we just need to check that the proposal densities and the acceptance probabilities for updating all parameter blocks are positive, then the ϕ-irreducibility is proved under 1-step transition (take n=1). For 1-dimensional MCMC, it is obvious that,
κ1(x,A)=∫Aq(x,y)a(x,y)ϕ(dy)>0for∀x∈X,A∈B(X)
as long as the proposal density q(x,y)>0 and the acceptance probability a(x,y)>0 always. To extend this conclusion to the multi-dimensional case, the technique of induction could be applied. See details in Section 4.1.8 of [22]. Among 12 parameter blocks {μ1C,μY,β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD}, if all existed, the updating of A and B requires the Metropolis steps, and that for the rest of parameter blocks only needs normal Gibbs steps by the conjugacy. Since the exponential function and the vector Normal density function are always positive, we can conclude that the acceptance probabilities and the proposal densities of the Metropolis steps, if present, are always positive. For the Gibbs steps, the proposal distributions are Matrix normal and Inverse-Wishart distributions, both with densities always positive. The acceptance probabilities are always 1, which is positive as well. This completes the proof for the ϕ-irreducibility.
Next, we intend to show the Aperiodicity, whose definition for a ϕ-irreducible Markov chain is given in Definition 4. (Aperiodicity).
For a ϕ-irreducible Markov chain(Xn)n=1∞, suppose that the state space isXwith σ-algebraF, the (1-step) transition probability isκ1:X×F→[0,1], and the stationary probability distribution isπ(·). The period of(Xn)n=1∞is defined as the largest positive integer T for which there exist disjoint subsetsA1,A2,…,AT∈Fwithπ(Ai)>0, such thatκ1(x,Ai+1)=1for allx∈Ai(1≤i≤T−1), andκ1(x,A1)=1for allx∈AT. The Markov chain is called Aperiodic ifT=1.
Aperiodicity could be immediately obtained, since every set with nonzero probability (under the stationary probability distribution π) in the σ-algebra of the state space could be reached from any point in the state space through one step transition in the Markov chain.
Lastly, we intend to prove the Harris recurrence. The definition of the Harris recurrence for a ϕ-irreducible Markov chain is given in Definition 5. (Harris recurrence).
For a ϕ-irreducible Markov chain(Xn)n=1∞, suppose that the state space isXwith σ-algebraF, and the stationary probability distribution isπ(·). The Markov chain is called Harris recurrent, if for allA⊆Xwithπ(A)>0and any starting pointx∈X, we have the conditional probabilityP(Xn∈Ainfinitelyofteninn∣X1=x)=1.
When dX∈{0,pC} and dY∈{0,r}, A and B don’t exist, then the Metropolis-within-Gibbs algorithm is actually a normal Gibbs sampler. Then Harris recurrence is directly implied by ϕ-irreducibility ([48], Corollary 13). When either 0<dX<pC or 0<dY<r, there exists Metropolis step for either A or B. The Harris recurrence could be obtained by proving the Lebesgue integrability of the posterior density over every combination of parameter blocks ([48], Corollary 19). For simplicity, we assume 0<dX<pC and 0<dY<r. Hence we only need to prove the posterior density f(μ1C,μY,β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D) is Lebesgue integrable over every k(1≤k≤12) elements from {μ1C,μY,β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD}. Theorem 2 under the other two cases (1) dX∈{0,pC} and 0<dY<r; (2) 0<dX<pC and dY∈{0,r} could be similarly proved with simple modifications of the following proof. Theorem 1 should be immediately proved by letting k to be the maximal value 12, i.e., the Lebesgue integrability over all parameter blocks is considered. Below we show the proof. The unnormalized posterior density is
f(μ1C,μY,β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D)∝f0(μ1C,μY,β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D)=exp{−12tr{(X˜1C−X˜1Dγ)(L(A)ΩL(A)T+L0(A)Ω0L0(A)T)−1(X˜1C−X˜1Dγ)T}}exp{−12tr{(Y˜−X˜1CL(A)ηCTR(B)T−X˜1DηDTR(B)T−X˜2β2)(R(B)ΦR(B)T+R0(B)Φ0R0(B)T)−1(Y˜−X˜1CL(A)ηCTR(B)T−X˜1DηDTR(B)T−X˜2β2)T}}|Ω|−(n+pD+wX+dX+1)/2|Ω0|−(n+pD+wX0+pC−dX+1)/2|Φ|−(n+p2+wY+dY+dX+pD+1)/2|Φ0|−(n+p2+wY0+r−dY+1)/2exp{−12tr(ΨXΩ−1)}exp{−12tr(ΨX0Ω0−1)}exp{−12tr(ΨYΦ−1)}exp{−12tr(ΨY0Φ0−1)}exp{−12tr{ΣC∣D−1(γ−Λ−1F)TΛ(γ−Λ−1F)}}exp{−12tr{ΣY∣X−1(β2−M−1Z)TM(β2−M−1Z)}}exp{−12tr{ΣA−1(A−A0)TKA−1(A−A0)}}exp{−12tr{ΣB−1(B−B0)TKB−1(B−B0)}}exp{−12tr{E(ηC−WE−1)TΦ−1(ηC−WE−1)}}exp{−12tr{Q(ηD−TQ−1)TΦ−1(ηD−TQ−1)}}.
Specifically, in the following term, what inside trace function is a product of two positive semi-definite matrices, so
exp{−12tr{(Y−X˜1CL(A)ηCTR(B)T−X˜1DηDTR(B)T−X˜2β2)(R(B)ΦR(B)T+R0(B)Φ0R0(B)T)−1(Y−X˜1CL(A)ηCTR(B)T−X˜1DηDTR(B)T−X˜2β2)T(In−1nJn)}}≤1,
where Jn=1n1nT. Therefore,
∫exp{−12tr{(Y−1nμYT−X˜1CL(A)ηCTR(B)T−X˜1DηDTR(B)T−X˜2β2)(R(B)ΦR(B)T+R0(B)Φ0R0(B)T)−1(Y−1nμYT−X˜1CL(A)ηCTR(B)T−X˜1DηDTR(B)T−X˜2β2)T}}dμY=(2πn)r/2|Φ|1/2|Φ0|1/2exp{−12tr{(Y−X˜1CL(A)ηCTR(B)T−X˜1DηDTR(B)T−X˜2β2)(R(B)ΦR(B)T+R0(B)Φ0R0(B)T)−1(Y−X˜1CL(A)ηCTR(B)T−X˜1DηDTR(B)T−X˜2β2)T(In−1nJn)}}≤(2πn)r/2|Φ|1/2|Φ0|1/2.
Meanwhile,
exp{−12tr{(Y˜−X˜1CL(A)ηCTR(B)T−X˜1DηDTR(B)T−X˜2β2)(R(B)ΦR(B)T+R0(B)Φ0R0(B)T)−1(Y˜−X˜1CL(A)ηCTR(B)T−X˜1DηDTR(B)T−X˜2β2)T}}≤1.
We can conclude
f0(μ1C,μY,β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D)≤c1f1(0)(μ1C,β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D),
and
∫f0(μ1C,μY,β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D)dμY≤c1f1(1)(μ1C,β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D),
where c1=1+(2πn)r/2 and for Ξ=0,1,
f1(Ξ)(μ1C,β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D)=exp{−12tr{(X˜1C−X˜1Dγ)(L(A)ΩL(A)T+L0(A)Ω0L0(A)T)−1(X˜1C−X˜1Dγ)T}}|Ω|−(n+pD+wX+dX+1)/2|Ω0|−(n+pD+wX0+pC−dX+1)/2|Φ|−(n+p2+wY+dY+dX+pD+1−Ξ)/2|Φ0|−(n+p2+wY0+r−dY+1−Ξ)/2exp{−12tr(ΨXΩ−1)}exp{−12tr(ΨX0Ω0−1)}exp{−12tr(ΨYΦ−1)}exp{−12tr(ΨY0Φ0−1)}exp{−12tr{ΣC∣D−1(γ−Λ−1F)TΛ(γ−Λ−1F)}}exp{−12tr{ΣY∣X−1(β2−M−1Z)TM(β2−M−1Z)}}exp{−12tr{ΣA−1(A−A0)TKA−1(A−A0)}}exp{−12tr{ΣB−1(B−B0)TKB−1(B−B0)}}exp{−12tr{E(ηC−WE−1)TΦ−1(ηC−WE−1)}}exp{−12tr{Q(ηD−TQ−1)TΦ−1(ηD−TQ−1)}}.
Next, we just need to show f1(Ξ)(μ1C,β2,γ,Ω,Ω0,Φ,Φ0, A, B, ηC, ηD∣D) is integrable with respect to any 1≤k≤11 elements in {μ1C,β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD}, for any value of Ξ. Similarly, note that
∫exp{−12tr{(X1C−1nμ1CT−X˜1Dγ)(L(A)ΩL(A)T+L0(A)Ω0L0(A)T)−1(X1C−1nμ1CT−X˜1Dγ)T}}dμ1C=(2πn)pC/2|Ω|1/2|Ω0|1/2exp{−12tr{(X1C−X˜1Dγ)(L(A)ΩL(A)T+L0(A)Ω0L0(A)T)−1(X1C−X˜1Dγ)T(In−1nJn)}}≤(2πn)pC/2|Ω|1/2|Ω0|1/2,
since
exp{−12tr{(X1C−X˜1Dγ)(L(A)ΩL(A)T+L0(A)Ω0L0(A)T)−1(X1C−X˜1Dγ)T(In−1nJn)}}≤1.
Together with the fact that
exp{−12tr{(X˜1C−X˜1Dγ)(L(A)ΩL(A)T+L0(A)Ω0L0(A)T)−1(X˜1C−X˜1Dγ)T}}≤1,
we get for any Ξ=0,1,
f1(Ξ)(μ1C,β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D)≤c2f2(Ξ,0)(β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D),
and
∫f1(Ξ)(μ1C,β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D)dμ1C≤c2f2(Ξ,1)(β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D),
where c2=1+(2πn)pC/2 and for Ξ=0,1, ζ=0,1,
f2(Ξ,ζ)(β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D)=|Ω|−(n+pD+wX+dX+1−ζ)/2|Ω0|−(n+pD+wX0+pC−dX+1−ζ)/2|Φ|−(n+p2+wY+dY+dX+pD+1−Ξ)/2|Φ0|−(n+p2+wY0+r−dY+1−Ξ)/2exp{−12tr(ΨXΩ−1)}exp{−12tr(ΨX0Ω0−1)}exp{−12tr(ΨYΦ−1)}exp{−12tr(ΨY0Φ0−1)}exp{−12tr{ΣC∣D−1(γ−Λ−1F)TΛ(γ−Λ−1F)}}exp{−12tr{ΣY∣X−1(β2−M−1Z)TM(β2−M−1Z)}}exp{−12tr{ΣA−1(A−A0)TKA−1(A−A0)}}exp{−12tr{ΣB−1(B−B0)TKB−1(B−B0)}}exp{−12tr{E(ηC−WE−1)TΦ−1(ηC−WE−1)}}exp{−12tr{Q(ηD−TQ−1)TΦ−1(ηD−TQ−1)}}.
Then it is enough to show f2(Ξ,ζ)(β2,γ,Ω,Ω0,Φ,Φ0,A,B, ηC, ηD∣D) is integrable over every 1≤k≤10 members in {β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD}, for any value of Ξ and ζ. Observe
∫exp{−12tr{E(ηC−WE−1)TΦ−1(ηC−WE−1)}}dηC=(2π)dXdY/2|E|−dY/2|Φ|dX/2
and
exp{−12tr{E(ηC−WE−1)TΦ−1(ηC−WE−1)}}≤1.
Hence for any Ξ=0,1, ζ=0,1,
f2(Ξ,ζ)(β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D)≤c3f3(Ξ,ζ,0)(β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηD∣D),
and
∫f2(Ξ,ζ)(β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηC,ηD∣D)dηC≤c3f3(Ξ,ζ,1)(β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηD∣D),
with c3=1+(2π)dXdY/2|E|−dY/2. For Ξ, ζ, λ=0,1,
f3(Ξ,ζ,λ)(β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηD∣D)=|Ω|−(n+pD+wX+dX+1−ζ)/2|Ω0|−(n+pD+wX0+pC−dX+1−ζ)/2|Φ|−(n+p2+wY+dY+pD+(1−λ)dX+1−Ξ)/2|Φ0|−(n+p2+wY0+r−dY+1−Ξ)/2exp{−12tr(ΨXΩ−1)}exp{−12tr(ΨX0Ω0−1)}exp{−12tr(ΨYΦ−1)}exp{−12tr(ΨY0Φ0−1)}exp{−12tr{ΣC∣D−1(γ−Λ−1F)TΛ(γ−Λ−1F)}}exp{−12tr{ΣY∣X−1(β2−M−1Z)TM(β2−M−1Z)}}exp{−12tr{ΣA−1(A−A0)TKA−1(A−A0)}}exp{−12tr{ΣB−1(B−B0)TKB−1(B−B0)}}exp{−12tr{Q(ηD−TQ−1)TΦ−1(ηD−TQ−1)}}.
Therefore, the task right now is to prove f3(Ξ,ζ,λ)(β2,γ, Ω, Ω0, Φ, Φ0, A, B, ηD∣D) is integrable over any 1≤k≤9 elements in {β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηD}, for any value of Ξ, ζ and λ. Since
∫exp{−12tr{Q(ηD−TQ−1)TΦ−1(ηD−TQ−1)}}dηD=(2π)dYpD/2|Q|−dY/2|Φ|pD/2
and
exp{−12tr{Q(ηD−TQ−1)TΦ−1(ηD−TQ−1)}}≤1,
we get for any Ξ, ζ, λ=0,1,
f3(Ξ,ζ,λ)(β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηD∣D)≤c4f4(Ξ,ζ,λ,0)(β2,γ,Ω,Ω0,Φ,Φ0,A,B∣D),
and
∫f3(Ξ,ζ,λ)(β2,γ,Ω,Ω0,Φ,Φ0,A,B,ηD∣D)dηD≤c4f4(Ξ,ζ,λ,1)(β2,γ,Ω,Ω0,Φ,Φ0,A,B∣D),
with c4=1+(2π)dYpD/2|Q|−dY/2 and for Ξ, ζ, λ, τ=0,1,
f4(Ξ,ζ,λ,τ)(β2,γ,Ω,Ω0,Φ,Φ0,A,B∣D)=|Ω|−(n+pD+wX+dX+1−ζ)/2|Ω0|−(n+pD+wX0+pC−dX+1−ζ)/2|Φ|−(n+p2+wY+dY+(1−τ)pD+(1−λ)dX+1−Ξ)/2|Φ0|−(n+p2+wY0+r−dY+1−Ξ)/2exp{−12tr(ΨXΩ−1)}exp{−12tr(ΨX0Ω0−1)}exp{−12tr(ΨYΦ−1)}exp{−12tr(ΨY0Φ0−1)}exp{−12tr{ΣC∣D−1(γ−Λ−1F)TΛ(γ−Λ−1F)}}exp{−12tr{ΣY∣X−1(β2−M−1Z)TM(β2−M−1Z)}}exp{−12tr{ΣA−1(A−A0)TKA−1(A−A0)}}exp{−12tr{ΣB−1(B−B0)TKB−1(B−B0)}}.
Then we intend to show f4(Ξ,ζ,λ,τ)(β2,γ,Ω,Ω0,Φ,Φ0, A, B∣D) is integrable over any 1≤k≤8 elements in {β2,γ,Ω,Ω0,Φ,Φ0,A,B}, for any value of Ξ, ζ, λ and τ. Since
∫exp{−12tr{ΣY∣X−1(β2−M−1Z)TM(β2−M−1Z)}}dβ2=(2π)p2r/2|M|−r/2|Φ|p2/2|Φ0|p2/2,
and
exp{−12tr{ΣY∣X−1(β2−M−1Z)TM(β2−M−1Z)}}≤1,
we get for any Ξ, ζ, λ, τ=0,1,
f4(Ξ,ζ,λ,τ)(β2,γ,Ω,Ω0,Φ,Φ0,A,B∣D)≤c5f5(Ξ,ζ,λ,τ,0)(γ,Ω,Ω0,Φ,Φ0,A,B∣D),
and
∫f4(Ξ,ζ,λ,τ)(β2,γ,Ω,Ω0,Φ,Φ0,A,B∣D)dβ2≤c5f5(Ξ,ζ,λ,τ,1)(γ,Ω,Ω0,Φ,Φ0,A,B∣D),
where c5=1+(2π)p2r/2|M|−r/2 and for Ξ, ζ, λ, τ, χ=0,1,
f5(Ξ,ζ,λ,τ,χ)(γ,Ω,Ω0,Φ,Φ0,A,B∣D)=|Ω|−(n+pD+wX+dX+1−ζ)/2|Ω0|−(n+pD+wX0+pC−dX+1−ζ)/2|Φ|−(n+wY+dY+(1−χ)p2+(1−τ)pD+(1−λ)dX+1−Ξ)/2|Φ0|−(n+wY0+r−dY+(1−χ)p2+1−Ξ)/2exp{−12tr(ΨXΩ−1)}exp{−12tr(ΨX0Ω0−1)}exp{−12tr(ΨYΦ−1)}exp{−12tr(ΨY0Φ0−1)}exp{−12tr{ΣC∣D−1(γ−Λ−1F)TΛ(γ−Λ−1F)}}exp{−12tr{ΣA−1(A−A0)TKA−1(A−A0)}}exp{−12tr{ΣB−1(B−B0)TKB−1(B−B0)}}.
Now it is suffice to show f5(Ξ,ζ,λ,τ,χ)(γ,Ω,Ω0,Φ,Φ0,A,B∣D) is integrable with respect to any 1≤k≤7 members in {γ,Ω,Ω0,Φ,Φ0,A,B}, for any value of Ξ, ζ, λ, τ and χ. Applying the similar trick as before, since
∫exp{−12tr{ΣC∣D−1(γ−Λ−1F)TΛ(γ−Λ−1F)}}dγ=(2π)pDpC/2|Λ|−pC/2|Ω|pD/2|Ω0|pD/2,
and
exp{−12tr{ΣC∣D−1(γ−Λ−1F)TΛ(γ−Λ−1F)}}≤1,
we conclude for any Ξ, ζ, λ, τ, χ=0,1,
f5(Ξ,ζ,λ,τ,χ)(γ,Ω,Ω0,Φ,Φ0,A,B∣D)≤c6f6(Ξ,ζ,λ,τ,χ,0)(Ω,Ω0,Φ,Φ0,A,B∣D),
and
∫f5(Ξ,ζ,λ,τ,χ)(γ,Ω,Ω0,Φ,Φ0,A,B∣D)dγ≤c6f6(Ξ,ζ,λ,τ,χ,1)(Ω,Ω0,Φ,Φ0,A,B∣D),
with c6=1+(2π)pDpC/2|Λ|−pC/2 and for Ξ, ζ, λ, τ, χ, ν=0,1,
f6(Ξ,ζ,λ,τ,χ,ν)(Ω,Ω0,Φ,Φ0,A,B∣D)=|Ω|−(n+wX+dX+pD(1−ν)+1−ζ)/2|Ω0|−(n+wX0+pC−dX+pD(1−ν)+1−ζ)/2|Φ|−(n+wY+dY+(1−χ)p2+(1−τ)pD+(1−λ)dX+1−Ξ)/2|Φ0|−(n+wY0+r−dY+(1−χ)p2+1−Ξ)/2exp{−12tr(ΨXΩ−1)}exp{−12tr(ΨX0Ω0−1)}exp{−12tr(ΨYΦ−1)}exp{−12tr(ΨY0Φ0−1)}exp{−12tr{ΣA−1(A−A0)TKA−1(A−A0)}}exp{−12tr{ΣB−1(B−B0)TKB−1(B−B0)}}.
Notice this is a constant multiple of the product of densities for IWdX(ΨX,n+wX+pD(1−ν)−ζ), IWpC−dX(ΨX0,n+wX0+pD(1−ν)−ζ), IWdY(ΨY,n+wY+p2(1−χ)+pD(1−τ)+dX(1−λ)−Ξ), IWr−dY(ΨY0,n+wY0+p2(1−χ)−Ξ), MNpC−dX.dX(A0,KA,ΣA) and MNr−dY.dY(B0,KB,ΣB). These Inverse-Wishart distributions are well-defined by checking
n+wX+pD(1−ν)−ζ>dX−1,n+wX0+pD(1−ν)−ζ>pC−dX−1,n+wY0+p2(1−χ)−Ξ>r−dY−1,
and
n+wY+p2(1−χ)+pD(1−τ)+dX(1−λ)−Ξ>dY−1,
for any value of Ξ, ζ, λ, τ, χ, ν=0,1, since wX>dX−1, wX0>pC−dX−1, wY>dY−1, wY0>r−dY−1, pD, p2, dX≥0 and n≥1. This completes the proof. □
Review of the Matrix Normal and the Inverse-Wishart Distributions
We introduce the Matrix normal and the Inverse-Wishart distributions here, for the description of the prior and the posterior distributions:
Let X∼MNa,b(M,U,V) denote the Matrix normal distribution with location parameter M∈Ra×b and scale parameters U∈Ra×a and V∈Rb×b such that U and V are positive definite, and it has the density function
exp{−tr(V−1(X−M)TU−1(X−M))/2}(2π)ab/2|V|a/2|U|b/2.
Let X∼IWa(X0,ν) denote the Inverse-Wishart distribution with scale matrix parameter X0∈S+a×a and degrees of freedom parameter ν>a−1. The density function is given by
|X0|ν/22νa/2Γa(ν/2)|X|−(ν+a+1)/2exp{−tr(X0X−1)/2}.
Here Γa(·) denotes the Multivariate gamma function, which is defined by
Γa(u)=πa(a−1)/4∏j=1aΓ{u+(1−j)/2},foru>0,
where a is a positive integer, and Γ{·} is the usual Gamma function.
Review of DIC, WAIC and Bayesian CV
We introduce the definitions of two information criteria, DIC and WAIC, and Bayesian CV here. Let ΘˆdX,dYmean denote the estimated posterior mean of all parameters Θ among retained posterior samples, when fixing the envelope dimensions at dX and dY. Recall that in (4.1), we denote l(Θ) as the log-likelihood function of Θ conditional on all observations D. For convenience, for i=1,2,…,n, we further denote f(Di∣Θ)=exp{li(Θ)} as the density function value of the observation i when using parameter Θ, and its log-transformation is denoted as li(Θ), where Di={Yi,X1C,i,X1D,i,X2,i}. Then DIC and WAIC are defined as
DIC=−2l(ΘˆdX,dYmean)+2{l(ΘˆdX,dYmean)−1S∑s=1Sl(ΘdX,dY(s))},WAIC=−2∑i=1nlog{1S∑s=1Sf(Di∣ΘdX,dY(s))}+2∑i=1n{1S−1∑s=1S(li(ΘdX,dY(s))−1S∑s=1Sli(ΘdX,dY(s)))2},
where {ΘdX,dY(s)}s=1S are the retained posterior samples, and we use the subscript for Θ to highlight that the envelope dimensions are fixed at dX and dY. Compared with AIC−MCMC as defined in (6.1), DIC replaces ΘˆdX,dYmax and K(dX,dY) with ΘˆdX,dYmean and a bias correction term based on data respectively, while WAIC replaces the maximized log-likelihood function by a term called the computed log pointwise posterior predictive density, and estimates K(dX,dY) by another bias correction term.
To introduce Bayesian CV, suppose n samples are partitioned into K subsets, with corresponding index sets as {Ck}k=1K. For fixed dX and dY, and each k=1,2,…,K, we run the MCMC algorithm on all observations D except those in Ck, and denote the retained posterior samples as {ΘdX,dY−k,(s)}s=1S. Instead of computing the MSPE that is usually used in the frequentist setting, Bayesian CV aims at finding dX and dY that minimize the minus of the average out-of-sample estimate of the log predictive density, with expression
−1K∑k=1Klog{1S∑s=1S∏i∈Ckf(Di∣ΘdX,dY−k,(s))}.
MCMC Diagnostic and Posterior Density Plots for Section <xref rid="j_nejsds23_s_017">7.2</xref>
The trace plots (left), the autocorrelation plots (middle; based only on Chain 1) and the evolution of the (median and 97.5% upper confidence limit of) Gelman-Rubin’s shrink factor as the number of iterations increases (right; PSRF represents the point estimate, i.e. the potential scale reduction factor) for a randomly selected and representative element of each of βˆ1C, βˆ1D, βˆ2, ΣˆC∣D and ΣˆY∣X from SIMP (correspondingly from row 1 to row 5), under the data generating mechanism (M1) of Section 7.2, with r=15 and n=300 (the most challenging case).
The (marginal) empirical posterior density plots for a randomly selected and representative element of each of βˆ1C and βˆ1D from SIMP, under the data generating mechanism (M1) of Section 7.2, with r=15 and n=300. These density plots under other simulated settings are similar.
Figure G.1 shows the convergence of the proposed MCMC algorithm by showing the traceplots, the autocorrelation plots and the evolution of the Gelman-Rubin’s shrink factor for a randomly selected and representative element of each of βˆ1C, βˆ1D, βˆ2, ΣˆC∣D and ΣˆY∣X from SIMP, under the data generating mechanism (M1) of Section 7.2, with r=15 and n=300 (the most challenging case). All of those plots indicate the convergence and weak autocorrelation between MCMC samples.
Figure G.2 displays the empirical posterior density plots of the randomly selected and representative element of each of βˆ1C and βˆ1D under the same setting. These (estimated) posterior distributions are bell-shaped and the shape is very common among the (estimated) posterior distributions of all elements in βˆ1C and βˆ1D.
Plots of Exploring the Population Stratification for Section <xref rid="j_nejsds23_s_023">8</xref>
Plots of all selected participants with their i-th SNP PC and (i+1)-th SNP PC as X and Y coordinates respectively, (i=1,3,5,7,9,11 for the upper left, upper middle, upper right, lower left, lower middle and lower right panels. The percentage of the total variation that is explained by that PC is included in the parenthesis) for the real data application. No population stratification is observed in all six plots.
MCMC Diagnostic and Posterior Density Plots for Section <xref rid="j_nejsds23_s_027">8.4</xref>
The trace plots (left), the autocorrelation plots (middle; based only on Chain 1) and the evolution of the (median and 97.5% upper confidence limit of) Gelman-Rubin’s shrink factor as the number of iterations increases (right; PSRF represents the point estimate, i.e. the potential scale reduction factor) for a randomly selected and representative element of each of βˆ1C, βˆ1D, βˆ2, ΣˆC∣D and ΣˆY∣X from SIMP (correspondingly from row 1 to row 5) for the real data application.
The (marginal) empirical posterior density plots for a randomly selected and representative element of each of βˆ1C and βˆ1D for the real data application.
Figures H.2 and H.3 are the MCMC diagnostic plots and the (marginal) empirical posterior density plots for the real data application. The observations are similar to those in Appendix G.
Besides APOE ϵ4, the following 37 SNPs are found to be significant with all 12 IPs from SIMP. The RefSeq gene is the associated (or the nearest if known) gene for the SNP in the NCBI Reference Sequence Database, where “-” means that no such associated gene is found by us.
SNP
RefSeq gene
SNP
RefSeq gene
rs10493854
-
rs7120548
MTCH2
rs12059556
-
rs10501927
CNTN5
rs1229119
-
rs489243
-
rs6428182
-
rs10894473
NTM
rs10921394
-
rs11064498
C1S
rs12138394
-
rs10877700
-
rs12756397
-
rs757402
OAS2
rs11685593
LOC105373605/BIN1
rs845757
SPATA7
rs7561528
LOC105373605/BIN1
rs2274736
PTPN21
rs11706690
CHL1
rs10498633
SLC24A4/RIN3
rs10513391
P2RY14
rs2554389
ADAMTSL3
rs1623177
-
rs437649
LOC105371501
rs10515893
-
rs470268
-
rs2439538
TBC1D7
rs17809911
CCDC102B
rs9381563
CD2AP
rs9959820
-
rs1516895
-
rs4809760
SLC9A8
rs11038268
LOC105376650
rs2281223
SLC9A8
rs2280231
NDUFS3
Boxplot of Bilateral Correlations for Section <xref rid="j_nejsds23_s_027">8.4</xref>
Boxplot of the Pearson correlation coefficients within all 6 pairs of brain measures over two hemispheres. Strong bilateral correlations are observed in this figure.
Acknowledgements
The authors thank the Editor, Associate Editor and two reviewers for insightful comments. Data used in preparation of this article were obtained from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). As such, the investigators within the ADNI contributed to the design and implementation of ADNI and/or provided data but did not participate in analysis or writing of this report. A complete listing of ADNI investigators can be found at: http://adni.loni.usc.edu/wp-content/uploads/how_to_apply/ADNI_Acknowledgement_List.pdf.
ReferencesBarnard, J., McCulloch, R. and Meng, X.-L. (2000). Modeling covariance matrices in terms of standard deviations and correlations, with application to shrinkage. Statistica Sinica1281–1311.Boada, M., Antunez, C., Ramírez-Lorca, R., DeStefano, A. L., Gonzalez-Perez, A., Gayán, J., López-Arrieta, J., Ikram, M. A., Hernández, I., Marin, J.et al. (2014). ATP5H/KCTD2 locus is associated with Alzheimer’s disease risk. Molecular Psychiatry19(6) 682–687.Broce, I. J., Tan, C. H., Fan, C. C., Jansen, I., Savage, J. E., Witoelar, A., Wen, N., Hess, C. P., Dillon, W. P., Glastonbury, C. M.et al. (2019). Dissecting the genetic relationship between cardiovascular risk factors and Alzheimer’s disease. Acta Neuropathologica137(2) 209–226.Bura, E. and Cook, R. D. (2003). Rank estimation in reduced-rank regression. Journal of Multivariate Analysis87(1) 159–176.Chakraborty, S. and Su, Z. A comprehensive Bayesian framework for envelope models. Technical Report.Chamberlain, G. (1982). Multivariate regression models for panel data. Journal of Econometrics18(1) 5–46.Chang, C. C., Chow, C. C., Tellier, L. C., Vattikuti, S., Purcell, S. M. and Lee, J. J. (2015). Second-generation PLINK: Rising to the challenge of larger and richer datasets. GigaScience4(1) 13742–015.Chen, T., Su, Z., Yang, Y., Ding, S.et al. (2020). Efficient estimation in expectile regression using envelope models. Electronic Journal of Statistics14(1) 143–173.Chib, S. and Greenberg, E. (1998). Analysis of multivariate probit models. Biometrika85(2) 347–361.Conway, J. (1990) A course in functional analysis. New York: Springer.Cook, R. D. and Zhang, X. (2015). Foundations for envelope models and methods. Journal of the American Statistical Association110(510) 599–611.Cook, R. D., Helland, I. S. and Su, Z. (2013). Envelopes and partial least squares regression. Journal of the Royal Statistical Society: Series B75(5) 851–877.Cook, R. D. (2018). An introduction to envelopes: Dimension reduction for efficient estimation in multivariate statistics.Cook, R. D. and Zhang, X. (2015). Simultaneous envelopes for multivariate linear regression. Technometrics57(1) 11–25.Cook, R. D., Li, B. and Chiaromonte, F. (2010). Envelope models for parsimonious and efficient multivariate linear regression (with discussion). Statistica Sinica20927–1010.De Jong, S. (1993). SIMPLS: An alternative approach to partial least squares regression. Chemometrics and Intelligent Laboratory Systems18(3) 251–263.Dibble, C. C., Elis, W., Menon, S., Qin, W., Klekota, J., Asara, J. M., Finan, P. M., Kwiatkowski, D. J., Murphy, L. O. and Manning, B. D. (2012). TBC1D7 is a third subunit of the TSC1-TSC2 complex upstream of mTORC1. Molecular Cell47(4) 535–546.Ding, S. and Dennis Cook, R. (2018). Matrix variate regressions and envelope models. Journal of the Royal Statistical Society: Series B (Statistical Methodology)80(2) 387–408.Ding, S., Su, Z., Zhu, G. and Wang, L. (2020). Envelope quantile regression. Statistica Sinica31(1) 79–105.Dubois, B., Hampel, H., Feldman, H. H., Scheltens, P., Aisen, P., Andrieu, S., Bakardjian, H., Benali, H., Bertram, L., Blennow, K.et al. (2016). Preclinical Alzheimer’s disease: Definition, natural history, and diagnostic criteria. Alzheimer’s & Dementia12(3) 292–323.Fischl, B. (2012). FreeSurfer. Neuroimage62(2) 774–781.Geyer, C. J. (1998). Markov chain Monte Carlo lecture notes. Course notes, Spring Quarter80.Greenlaw, K., Szefer, E., Graham, J., Lesperance, M., Nathoo, F. S. and Alzheimer’s Disease Neuroimaging Initiative (2017). A Bayesian group sparse multi-task regression model for imaging genetics. Bioinformatics33(16) 2513–2522.Harold, D., Abraham, R., Hollingworth, P., Sims, R., Gerrish, A., Hamshere, M. L., Pahwa, J. S., Moskvina, V., Dowzell, K., Williams, A.et al. (2009). Genome-wide association study identifies variants at CLU and PICALM associated with Alzheimer’s disease. Nature Genetics41(10) 1088–1093.Hibar, D. P., Stein, J. L., Kohannim, O., Jahanshad, N., Saykin, A. J., Shen, L., Kim, S., Pankratz, N., Foroud, T., Huentelman, M. J.et al. (2011). Voxelwise gene-wide association study (vGeneWAS): Multivariate gene-based association testing in 731 elderly subjects. NeuroImage56(4) 1875–1891.Hotelling, H. (1936). Relations between two sets of variates. Biometrika28(3-4) 321–377.Howie, B. N., Donnelly, P. and Marchini, J. (2009). A flexible and accurate genotype imputation method for the next generation of genome-wide association studies. PLOS Genetics5(6) 1000529.Howie, B., Fuchsberger, C., Stephens, M., Marchini, J. and Abecasis, G. R. (2012). Fast and accurate genotype imputation in genome-wide association studies through pre-phasing. Nature Genetics44(8) 955–959.Izenman, A. J. (1975). Reduced-rank regression for the multivariate linear model. Journal of Multivariate Analysis5(2) 248–264.Jansen, I. E., Savage, J. E., Watanabe, K., Bryois, J., Williams, D. M., Steinberg, S., Sealock, J., Karlsson, I. K., Hägg, S., Athanasiu, L.et al. (2019). Genome-wide meta-analysis identifies new loci and functional pathways influencing Alzheimer’s disease risk. Nature Genetics51(3) 404–413.Karch, C. M., Ezerskiy, L. A., Bertelsen, S., (ADGC), A. D. G. C. and Goate, A. M. (2016). Alzheimer’s disease risk polymorphisms regulate gene expression in the ZCWPW1 and the CELF1 loci. PLOS One11(2) 0148717.Kendall, M. G. (1957). A course in multivariate analysis. Technical Report.Khare, K., Pal, S. and Su, Z. (2017). A Bayesian approach for envelope models. The Annals of Statistics45(1) 196–222.Kim, S., Sohn, K.-A. and Xing, E. P. (2009). A multivariate regression approach to association analysis of a quantitative trait network. Bioinformatics25(12) 204–212.Kumar, N., Bansal, A., Sarma, G. and Rawal, R. K. (2014). Chemometrics tools used in analytical chemistry: An overview. Talanta123186–199.Lambert, J.-C., Ibrahim-Verbaas, C. A., Harold, D., Naj, A. C., Sims, R., Bellenguez, C., Jun, G., DeStefano, A. L., Bis, J. C., Beecham, G. W.et al. (2013). Meta-analysis of 74,046 individuals identifies 11 new susceptibility loci for Alzheimer’s disease. Nature Genetics45(12) 1452–1458.Lee, M. and Su, Z. (2020). A review of envelope models. International Statistical Review88(3) 658–676.Lee, M., Chakraborty, S. and Su, Z. (2022). A Bayesian approach to envelope quantile regression. Statistica Sinica321–19.Li, L. and Zhang, X. (2017). Parsimonious tensor response regression. Journal of the American Statistical Association112(519) 1131–1146.Liu, J. S. and Wu, Y. N. (1999). Parameter expansion for data augmentation. Journal of the American Statistical Association94(448) 1264–1274.Naj, A. C., Leonenko, G., Jian, X., Grenier-Boley, B., Dalmasso, M. C., Bellenguez, C., Sha, J., Zhao, Y., van der Lee, S. J., Sims, R. et al. (2021). Genome-wide meta-analysis of late-onset Alzheimer’s disease using rare variant imputation in 65,602 subjects identifies novel rare variant locus NCK2: The International Genomics of Alzheimer’s Project (IGAP). medRxiv.Nathoo, F. S., Kong, L., Zhu, H. and Alzheimer’s Disease Neuroimaging Initiative (2019). A review of statistical methods in imaging genetics. Canadian Journal of Statistics47(1) 108–131.Nestor, S. M., Rupsingh, R., Borrie, M., Smith, M., Accomazzi, V., Wells, J. L., Fogarty, J., Bartha, R. and Alzheimer’s Disease Neuroimaging Initiative (2008). Ventricular enlargement as a possible measure of Alzheimer’s disease progression validated using the Alzheimer’s disease neuroimaging initiative database. Brain131(9) 2443–2454.Noorossana, R., Eyvazian, M., Amiri, A. and Mahmoud, M. A. (2010). Statistical monitoring of multivariate multiple linear regression profiles in phase I with calibration application. Quality and Reliability Engineering International26(3) 291–303.Park, Y., Su, Z. and Chung, D. (2022). Envelope-based partial partial least squares with application to cytokine-based biomarker analysis for COVID-19. Statistics in Medicine 1–15.Park, Y., Su, Z. and Zhu, H. (2017). Groupwise envelope models for imaging genetic analysis. Biometrics73(4) 1243–1253.Poulin, S. P., Dautoff, R., Morris, J. C., Barrett, L. F., Dickerson, B. C. and Alzheimer’s Disease Neuroimaging Initiative (2011). Amygdala atrophy is prominent in early Alzheimer’s disease and relates to symptom severity. Psychiatry Research: Neuroimaging194(1) 7–13.Roberts, G. O. and Rosenthal, J. S. (2006). Harris recurrence of Metropolis-within-Gibbs and trans-dimensional Markov chains. The Annals of Applied Probability16(4) 2123–2139.Rosenthal, S. L., Barmada, M. M., Wang, X., Demirci, F. Y. and Kamboh, M. I. (2014). Connecting the dots: Potential of data integration to identify regulatory SNPs in late-onset Alzheimer’s disease GWAS findings. PLOS One9(4) 95152.Rowe, B. Bayesian variable selection methods for genome wide association studies with categorical phenotypes (2020). PhD thesis, University of Nevada, Las Vegas.Selige, T., Böhner, J. and Schmidhalter, U. (2006). High resolution topsoil mapping using hyperspectral image and field data in multivariate regression modeling procedures. Geoderma136(1-2) 235–244.Song, Y., Ge, S., Cao, J., Wang, L. and Nathoo, F. S. (2022). A Bayesian spatial model for imaging genetics. Biometrics78(2) 742–753.Squillario, M., Tomasi, F., Tozzo, V., Barla, A., Uberti, D. and Alzheimer’s Disease Neuroimaging Initiative (2018). A 3-fold kernel approach for characterizing late-onset Alzheimer’s disease. bioRxiv 397760.Stein, J. L., Hua, X., Lee, S., Ho, A. J., Leow, A. D., Toga, A. W., Saykin, A. J., Shen, L., Foroud, T., Pankratz, N.et al. (2010). Voxelwise genome-wide association study (vGWAS). NeuroImage53(3) 1160–1174.Strickland, S. L., Reddy, J. S., Allen, M., N’songo, A., Burgess, J. D., Corda, M. M., Ballard, T., Wang, X., Carrasquillo, M. M., Biernacka, J. M.et al. (2020). MAPT haplotype–stratified GWAS reveals differential association for AD risk variants. Alzheimer’s & Dementia16(7) 983–1002.Su, Z., Zhu, G., Chen, X. and Yang, Y. (2016). Sparse envelope model: Efficient estimation and response variable selection in multivariate linear regression. Biometrika103(3) 579–593.Su, Z. and Cook, R. D. (2011). Partial envelopes for efficient estimation in multivariate linear regression. Biometrika98133–146.Su, Z. and Cook, R. D. (2013). Estimation of multivariate means with heteroscedastic errors using envelope models. Statistica Sinica23(1) 213–230.Talhouk, A., Doucet, A. and Murphy, K. (2012). Efficient Bayesian inference for multivariate probit models with sparse inverse correlation matrices. Journal of Computational and Graphical Statistics21(3) 739–757.Tan, M.-S., Yang, Y.-X., Xu, W., Wang, H.-F., Tan, L., Zuo, C.-T., Dong, Q., Tan, L., Suckling, J. and Yu, J.-T. (2021). Associations of Alzheimer’s disease risk variants with gene expression, amyloidosis, tauopathy, and neurodegeneration. Alzheimer’s Research & Therapy13(1) 1–11.Torvell, M., Carpanini, S. M., Daskoulidou, N., Byrne, R. A., Sims, R. and Morgan, B. P. (2021). Genetic Insights into the Impact of Complement in Alzheimer’s Disease. Genes12(12) 1990.Veera Manikandan, R. and Anand, R. (2015). P2-012: A genome wide scan for genetic variations with inverse association between Alzheimer’s disease and breast cancer. Alzheimer’s & Dementia11(7S_Part_10) 485–485.Visscher, P. M., Wray, N. R., Zhang, Q., Sklar, P., McCarthy, M. I., Brown, M. A. and Yang, J. (2017). 10 years of GWAS discovery: Biology, function, and translation. The American Journal of Human Genetics101(1) 5–22.Vounou, M., Nichols, T. E., Montana, G. and Alzheimer’s Disease Neuroimaging Initiative (2010). Discovering genetic associations with high-dimensional neuroimaging phenotypes: A sparse reduced-rank regression approach. NeuroImage53(3) 1147–1159.Wang, H., Nie, F., Huang, H., Risacher, S. L., Saykin, A. J., Shen, L. and Alzheimer’s Disease Neuroimaging Initiative (2012). Identifying disease sensitive and quantitative trait-relevant biomarkers from multidimensional heterogeneous imaging genetics data via sparse multimodal multitask learning. Bioinformatics28(12) 127–136.White, C. C., Yang, H.-S., Schneider, J. A., Bennett, D. A., De Jager, P. L. and group, C. A. F. (2021). A genome-wide investigation of clinicopathologic endophenotypes uncovers a new susceptibility locus for tau pathology at Neurotrimin (NTM). Alzheimer’s & Dementia17051682.Zhang, C. and Yu, T. (2008). Semiparametric detection of significant activation for brain fMRI. The Annals of Statistics36(4) 1693–1725.Zhu, H., Khondker, Z., Lu, Z. and Ibrahim, J. G. (2014). Bayesian generalized low rank regression models for neuroimaging phenotypes and genetic markers. Journal of the American Statistical Association109(507) 977–990.