- Email: [email protected]
Rejoinder to the Discussion
Journal of Statistical Planning and Inference 145 (2014) 42–48
Contents lists available at ScienceDirect
Journal of Statistical Planning and Inference journal homepage: www.elsevier.com/locate/jspi
Rejoinder to the Discussion Rabi Bhattacharya a,n, Vic Patrangenaru b a b
Department of Mathematics, The University of Arizona, Tucson, AZ 85721, USA Department of Statistics, The Florida State University, Tallahassee, FL 32304, USA
The discussants have raised important issues on the subject matter of the present article. We are deeply appreciative of their valuable inputs and thank them for sharing their insights with us. As one of the early contributors to the theory of nonparametric statistics on manifolds, and as one whose research has made significant impacts on various aspects of statistics, Professor Beran is in a unique position to provide a broad perspective of the area under discussion. With regard to the need of developing appropriate analogs of multivariate analysis on manifolds that he rightly stresses, two contributions deserve particular attention. The first one by Huckemann et al. (2010b) has provided a theory of generalized principal component analysis (GPCA) on Riemannian manifolds and, more generally, on any quotient of a Riemannian manifold under the action of a Lie group of isometries. These cover, in particular, Kendall's shape spaces ∑km (for all m Z 2). Here the generalized principal components (GPC-s) are generalized geodesics. In particular, on a Riemannian manifold, the first GPC is the geodesic that has the smallest expected squared distance from a random variable X on the manifold with the given distribution Q. The second GPC is the geodesic with the smallest expected squared distance from X among all geodesics which intersect the first GPC orthogonally, and so on. Next, the article by Huckemann et al. (2010a) lays some groundwork for MANOVA on Riemannian manifolds. In the absence of the usual translation providing linearity or the additive structure of classical MANOVA on Euclidean spaces, on a Riemannian manifold M one employs parallel transport of vector fields along geodesics, made possible by the Levi-Civita connection. Specifically, in a two-factor model the r-th observation P ij;r at levels i and j of factors 1 and 2, respectively, is lifted to the tangent space T μij M by the log (inverse exponential) map, where μij is the intrinsic mean of the distribution of P ij;r . Because for different i, j the tangent spaces are not naturally related, log μij P ij;r is then transported to a common tangent space T ν M for some suitable ν (for example, a pooled intrinsic mean), by parallel translation along the geodesic joining μij and ν. The resulting vectors ϵij;r can play the role of the residuals, or i.i.d. errors, in the linear model. On functional inference mentioned by Professor Beran (2013), nonparametric density estimation, regression and classification have been carried out on manifolds in Bhattacharya and Dunson (2010) and Bhattacharya and Bhattacharya (2012), Chapters 13 and 14, via nonparametric Bayes theory. For density estimation, for example, this requires the imposition of a Dirichlet prior on the space of mixtures (probabilities) of a suitable parametric family of densities. The class of all densities obtained by such mixtures is required to be dense (say, in L1) in the class of all densities. In simulations and data examples, the Bayes posterior so obtained seems to perform at least as well as, and often better than, frequentist estimates by kernel methods such as that of Pelletier (2005). As to the two methods of bootstrapping mentioned by Professor Beran for computing or estimating p-values in twosample tests, while Bhattacharya and Patrangenaru (2005) use bootstrap confidence regions, both the confidence region based method and one based on simulations directly of the null distribution are considered in Bhattacharya and Bhattacharya (2012, Chapter 4). It may, however, be pointed out that in the context manifolds which are often high dimensional, and with samples of rather moderate sizes that are normally available, bootstrapping of the chisquare statistic runs into problems due to degeneracy of the bootstrapped covariance matrix. Also, even when one has non-degenerate bootstrap covariance, bootstrap methods do not seem particularly suitable to estimate the extremely small p-values such as obtained in Applications 3.3, 3.4 and 4.2. One may of course question the validity of the classical chisquare approximation of such small probabilities. But these are used here (and in Bhattacharya and Bhattacharya, 2012, Chapter 2) mostly to indicate
n
Corresponding author. E-mail address: [email protected] (R. Bhattacharya).
0378-3758/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jspi.2013.08.006
R. Bhattacharya, V. Patrangenaru / Journal of Statistical Planning and Inference 145 (2014) 42–48
43
extremely small orders of p-values (beyond the usual domain of bootstrap calculations) which may be justified by the order of approximation of tail probabilities provided by deep inequalities such as that of von Bahr (1967, Theorem 1). We share some of Professor Beran's reservations about statistical inference that is commonly practiced, including the impossibility of drawing a truly random number. By its very nature statistical inference is inductive and deals with highly illposed inverse problems, unlike most mathematics which is deductive. As the science underlying a physical phenomenon is understood better the level of uncertainty is reduced and inference becomes more accurate. This is happening in biology. However, as statistical mechanics teaches us, certain uncertainty is inherent in nature. As illustrated in the discussion, it is important to realize the limitations of statistical methods. However, even the commonly used methods, and many nonparametric procedures among them, will always remain indispensable in analyzing a host of phenomena. We next turn to the discussion by Dryden et al. (2013). Professors Dryden and Le are among the distinguished group of scholars including J. Kent, K. Mardia, C. Goodall and M. Prentice, who developed the theory of parametric inference on manifolds, following pioneering work by D.G. Kendall, F. Bookstein and H. Ziezold. Their recent work, some of which is surveyed in their discussion, has contributed significantly to nonparametric inference based on Fréchet type means. The present discussion by them provides a number of valuable insights concerning the notion of the mean on a manifold and the asymptotics of its empirical version. As to the question “In praise of which mean?” in Dryden et al. (2013), let us first recall that in Bhattacharya and Patrangenaru (2003, 2005) and Bhattacharya and Bhattacharya (2012), as surveyed here, a general theory is developed on the asymptotics of intrinsic and extrinsic Fréchet means and their applications to nonparametric inference on general abstract manifolds, not just for Kendall's shape spaces. For example, it applies to Stiefel manifolds and spaces of affine and projective shapes (Mardia and Patrangenaru, 2005; Patrangenaru et al., 2010; Bhattacharya and Bhattacharya, 2012, Chapters 9–11). It encompasses Fréchet minimizers of expected squared distances with respect to arbitrary distances ρ, subject of course to appropriate constraints. The topological constraint (2.1), namely, “closed bounded subsets are compact”, is identified as the appropriate one for proving consistency under uniqueness (Theorem 2.1). In addition, some smoothness of the Fréchet function and nondegeneracy of its averaged Hessian are required for deriving the CLT and tests and confidence regions based on it (see Theorems 2.4, 2.12, 2.15–2.17). As pointed out in Bhattacharya and Patrangenaru (2005) and in Remark 2.5 here, Theorem 2.4 applies to both intrinsic and extrinsic sample Fréchet means. As we will specify later, an extension of this theorem due to Bhattacharya and Lin (2013) removes some unnecessary constraints and yields a general CLT for both intrinsic and extrinsic Fréchet means extending Corollary 2.6 and Proposition 2.10 on manifolds, and it also applies to several important examples of stratified spaces. It has been pointed out in Bhattacharya and Bhattacharya (2012), Sections 3.2, 8.5, that, under uniqueness, consistency holds with the squared distance ρ2 replaced by ρβ for β Z1 in the definition of the Fréchet function, and the CLT for the sample Fréchet mean may also hold for β 4 2. Why is then the case β ¼ 2 special? First, for extrinsic means as defined by us, precise (necessary and sufficient) conditions for uniqueness of the Fréchet minimizer may often be obtained, and an explicit analytical formula for it be found. This is not the case for β a 2. In particular, this clarifies the issue of consistency, or lack of it, of Procrustes type estimators (“partial”, “full”, “generalized”) for what was defined as the population mean shape in the earlier literature (see, e.g., Dryden and Mardia, 1998, p. 88), referred to as “the shape of the mean” by the present discussants. While a Procrustes estimator may be an extrinsic empirical (or, sample-) Fréchet mean, the quantity it was supposed to estimate, namely the mean shape as defined in the past literature, is not the corresponding Fréchet mean, leading to inconsistency unless some extraneous restriction such as isotropy is imposed. Also, for the intrinsic mean, for β o 2, especially for β ¼ 1, the uniqueness condition has been proven only under more restrictive conditions than for β ¼ 2 (See Afsari, 2011). Secondly, it is not clear that the smoothness condition required for the CLT holds in general for β o2, even when uniqueness of the minimizer is assured. Also, it is worth emphasizing again what is stated in our Introduction, namely, the embeddings that we use are equivariant under the action of a large group preserving a great deal of the geometry and this is undoubtedly a reason for the essential equivalence between the extrinsic and intrinsic inferences displayed in the Applications in the article (also see Chapter 2 of Bhattacharya and Bhattacharya, 2012). The difference between the two sample means and of estimates and test statistics (and p-values) based on them is of smaller magnitude by several orders than the average difference between the extrinsic and intrinsic distances computed for pairs of sample observations. This shows, first of all, that one may think of extrinsic inference as a computationally convenient very good approximation to intrinsic inference, provided the embedding for the extrinsic analysis is equivariant under a large group H, perhaps especially when H is a group of isometries on the manifold if endowed with the metric tensor induced by the embedding. Secondly, it provides some confidence in the validity of using the intrinsic mean in most real data examples, without being bothered by the vexing theoretical problem of uniqueness. On reflection-similarity shape spaces, under the MDS or Schoenberg embedding (see (4.2) in the present article), Dryden et al. (2013) provide significant results on the uniqueness and, consequently, consistency of empirical Fréchet means under a family of distances dα ðα Z 1=2Þ, with d1 being the induced Euclidean distance used in our definition of the extrinsic mean under the embedding. They have also provided explicit formulas for the corresponding Fréchet means. The case α ¼ 1 was solved by Bhattacharya (2008), who also provided the CLT for the empirical Fréchet mean for this case. It should be borne in mind that such a family of distances, based on a spectral decomposition, is not generally available for arbitrary embeddings of other manifolds. We now turn to the last section of the discussion by Dryden et al. (2013) on the limiting distribution of the intrinsic sample mean. It may be noted that while a general CLT for intrinsic sample means was derived in Bhattacharya and Patrangenaru (2005), a study linking the curvature of a Riemannian manifold and the dispersion of the limiting Normal
44
R. Bhattacharya, V. Patrangenaru / Journal of Statistical Planning and Inference 145 (2014) 42–48
distribution in the CLT appeared in Bhattacharya and Bhattacharya (2008). In the latter article an explicit expression of the dispersion matrix in terms of curvature is given for manifolds with constant sectional curvature, and a lower bound in terms of sectional curvature for arbitrary Riemannian manifolds. These results are extended in a fine article by Kendall and Le (2011), and they are stated in the discussion by Dryden et al. (2013). In particular, dispersion is expressed in terms of curvature for Kahler (or, complex-) manifolds with constant holomorphic curvature as well, with an explicit expression in the case of Kendall's planar shape space. A different form of this dispersion also appears in Bhattacharya and Bhattacharya (2008, Theorem 3.1) and remarks following it, but without linking this to curvature. As emphasized in the Introduction and elsewhere in our work, the Fréchet mean is used as a nonparametric tool, not so much for inference on the mean itself as it is for discrimination among distributions on manifolds, especially in two-sample problems. In exceptional cases, as in the case of directional statistics, the mean may be the main point of interest. Estimating the true direction of the earths magnetic poles in past geological time scales is an example, where the remanent fossil magnetism provides data with the true value of the direction corrupted by random errors or, what one may call, measuremental errors (Application 2.14). For testing solely via empirical Fréchet means if two random samples on a manifold are from different distributions, there is of course the usual theoretical problem of identification (namely, distributions with the same mean may be different). But in most cases of interest, such as shape spaces, the manifold is of high dimension and the Fréchet mean is shown to discriminate different distributions in two-sample problems quite well (Applications 3.3, 3.4, 4.2), often far better than parametric tests designed to test directly if the distributions are different, indicating perhaps that the parametric tests so carried out are highly misspecified. Because of the well known curse of dimensionality, in these cases the effectiveness of more elaborate nonparametric functional inference for discrimination requires much larger sample sizes than are generally available. Hence our phrase “in praise of the mean”. We do recognize the usefulness of any so-called mean which carries out this task of discrimination well. We now turn to the discussion by Huckemann (2013). Professor Huckeman's important contributions to multivariate analysis are mentioned at the beginning of the Rejoinder. In the discussion he proposes to use the Ziezold mean, following Ziezold (1994), which is defined on a quotient M ¼ N\G, where N is a Riemannian manifold embedded in a Euclidean space, and G is a Lie group of isometries on N which is also assumed to be an isometry on N with respect to the Euclidean distance on it. The Ziezold mean is the Fréchet mean with respect to the distance ρð½x; ½yÞ ¼ inf f J gx′y′ J : g A Gg, where x′; y′ are arbitrary elements in N of the orbits ½x, ½y A M, respectively, and J J is the Euclidean norm. If the action is free and proper then M is a manifold and one can define a Riemannian metric on it by Riemannian submersion of N into M, and the intrinsic distance would be different from the Ziezold distance. In general, M need not be a manifold. For example, M can be ∑km for any m Z 2. On ∑km the Ziezold distance is also known as the partial Procrustes distance (See Dryden and Mardia, 1998, p. 65). The comparisons between the Ziezold mean and the extrinsic mean for the 3D reflection shape by data analysis and simulation in Examples 3.1 and 3.2 are valid and indicate that if the underlying distributions assign positive probability near, or on, the singular part of the space then the extrinsic mean based on Schoenberg embedding is not going to be very effective as a tool for discrimination. We expect that the intrinsic mean on the whole space ∑k3 , however, will be at least as effective as the Ziezold mean in such cases. There remains the question of uniqueness of the Ziezold mean. The expectation that uniqueness of the Fréchet minimizer holds more broadly for this distance than for the intrinsic distance may be valid. But it needs to be explored. Professor Huckemann is quite right in emphasizing that the uniqueness problem is one of the main issues confronting the nonparametric theory using Fréchet means. It is here that the extrinsic mean holds an edge over other means. Still one expects uniqueness for both the intrinsic and Ziezold means to hold outside a relatively small set of distributions. On the discussion by Professor Kim et al. The recently National Institutes of Health sponsored Human Microbiome Project (HMP), aims to “characterize microbial communities found at several different sites on the human body, including nasal passages, oral cavities, skin, gastrointestinal tract, and urogenital tract, and to analyze the role of these microbes in human health and disease”. From a general statistician perspective, the HMP is directed toward understanding correlations between changes in the microbiome and human health. From a biologistʼs perspective, the HMP may be regarded as one small project of the larger prokaryote–eukaryotes relationship, and from a higher general scientific view, this relationship is more meaningful in presence of an understanding of and closeness to Mother Nature. Recall that viruses are the tiniest life forms known to date and they replicate inside living cells; archaea are unicellular microorganisms, microbes that have no cell nucleus or any other membrane-bound organelles within their cells, and bacteria are all other nucleus free cellular organisms (prokaryote). The organisms whose cells do have a nucleus are called eukaryotes. However as noted in the discussion by Kim et al. (2013), sequencing and classification are based on DNA structure and not on cellular organization of an organism. Classifying life forms is the focus of taxonomy, and the DNA tags attached to an OTU are commonly used in such classifications in recent times (see e.g. Ciccarelli et al., 2006, for the left hand side of Fig. 1). The suggestion made by Kim et al. (2013), of sketching possible directions that have not appeared in the literature for analyzing microbiome data, parallels similar projects of using phylogenetic analysis of DNA sequencing mostly associated with computational methods in Biology to study the microbiome; these include Holmes (2012), where in particular part of the human skin microbiome was discussed. In Fig. 1 right hand side, the human skin microbiome with types and proportions of prokaryote at various parts is schematically displayed. Since proportions and mean types of bacteria play a key role in general in the health state of eukaryotes, the suggestion of running a Fréchet mean analysis for microbiome data is welcome. As noted in the discussion by Dryden et al. (2013), the Fréchet mean depends on the distance used on the sample space. This is in particular true in the case of operational taxonomic units (OTUs) that are represented as leaves on phylogenetic trees, points on the tree spaces Tp.
R. Bhattacharya, V. Patrangenaru / Journal of Statistical Planning and Inference 145 (2014) 42–48
45
Fig. 1. Left: An oversimplified color-coded version of the “tree of life”: prokaryote (blue), archaea (green) and eukaryotes (red). Right: Human Skin Microbiome. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
Fig. 2. Left: T4, space of trees with four leaves. Right: Petersen graph.
We described in detail the tree spaces, and commented on their stratification. In the low dimensional cases p r4, the topological structure is well known as presented in our paper. For example T4 can be representation as a surface with singularities obtained from the polyhedral surface given on the right hand side of Fig. 2 by identifying the edges labeled with the same letter and then embedded in R4 : The intersection of the 15 quadrants with the unit L1 in R4 is given by the 15 edges of the Petersen graph (see left side of Fig. 2). This graph provides a convenient representation for data on T4. The Fréchet mean of a sample of size 10 of simulated phylogenetic trees obtained with respect to this embedding for the concatenated eukaryote Parkinsea RNA tree with four leafs data was computed in Ellingson et al. (2012). There they show the resulting Fréchet sample mean tree and its relative location in the tree space in terms of the Petersen graph. Computational algorithms for intrinsic sample means on Tp were given by Owen and Provan (2011). Recently Ellingson et al. (2012) also considered some other computational examples to illustrate the behavior of extrinsic means for simulated data on T4. For the asymptotic distribution of the intrinsic sample means on this space, we refer to Barden et al. (2013). These early results along the line of the Fréchet mean analysis on phylogenetic tree spaces suggested by Professor Kim et al. (2013) are encouraging for more computationally intensive examples arising in the prokaryote–eukaryote relationship. In conclusion we state a recent general CLT for empirical or sample Fréchet means on manifolds and stratified spaces, since asymptotic distribution theory is a major theme underlying the discussions. Let ðS; ρÞ be a metric space and Q a probability measure on its Borel sigma-field.
46
R. Bhattacharya, V. Patrangenaru / Journal of Statistical Planning and Inference 145 (2014) 42–48
Assume that the Fréchet function F under squared distance is finite on S, and that F has a unique minimizer μ ¼ argminp FðpÞ—the Fréchet mean of Q (with respect to the distance ρ). Under broad conditions, the Fréchet sample mean μn of the empirical distribution Q n ¼ 1=n∑nj¼ 1 δY j based on independent S-valued random variables Yj ðj ¼ 1; …; nÞ with common distribution Q is a consistent estimator of μ. That is, μn -μ almost surely, as n-1. Here μn may be taken to be any measurable selection from the (random) set of minimizers of the Fréchet function of Qn, namely F n ðpÞ ¼ ð1=nÞ∑nj¼ 1 ρ2 ðp; Y j Þ. We make the following assumptions:
(A1) The Fréchet mean μ of Q is unique. (A2) μA G, where G is a measurable subset of S, and there is a homeomorphism ϕ : G-U, where U is an open subset of Rs for some s Z1 and G is given its relative topology on S. Also, x↦hðx; qÞ≔ρ2 ðϕ1 ðxÞ; qÞ
ð1Þ
is twice continuously differentiable on U, for every q outside a Q-null set. (A3) μn A G almost surely, for all sufficiently large n. (A4) Let Dr hðx; qÞ ¼ ∂hðx; qÞ=∂xr ; r ¼ 1; …; s. Then EjDr hðϕðμÞ; Y 1 Þj2 o 1; EjDr;r′ hðϕðμÞ; Y 1 Þj o1
for r; r′ ¼ 1; …; s:
ð2Þ
(A5) Let ur;r′ ðϵ; qÞ ¼ supfjDr;r′ hðθ; qÞDr;r′ hðϕðμÞ; qÞj : jθϕðμÞj oϵg. Then Ejur;r′ ðϵ; Y 1 Þj-0
as ϵ-0
for all 1 r r; r′ r s:
ð3Þ
(A6) The matrix Λ ¼ ½EDr;r′ hðϕðμÞ; Y 1 Þr;r′ ¼ 1;…;s is nonsingular. Theorem 1. Under Assumptions (A1)–(A6), L
n1=2 ½ϕðμn ÞϕðμÞ-Nð0; Λ1 CΛ1 Þ;
ð4Þ
as n-1;
where C is the covariance matrix of fDr hðϕðμÞ; Y 1 Þ; r ¼ 1; …; sg. This theorem extends and improves the CLT-s given in Bhattacharya and Patrangenaru (2005) and Bhattacharya and Bhattacharya (2012), Chapters 3 and 4. It is derived in Bhattacharya and Lin (2013). It has the following immediate consequences. Corollary 2 (CLT for Intrinsic Means on a Manifold). Let (M,g) be a d-dimensional complete Riemannian manifold with metric tensor g and geodesic distance ρg . Suppose Q is a probability measure on M with intrinsic mean μI , which is the unique minimizer of the Fréchet function for the distance ρg . Let ϕ ¼ Exp μ1 be the inverse exponential, or log-, function at μI defined on a I neighborhood G of μ ¼ μI onto its image U in the tangent space T μI ðMÞ. Assume that the Assumptions (A4)–(A6) hold. Then, with s ¼d, the CLT (4) holds for the intrinsic sample mean μn ¼ μn;I . Remark 3. Corollary 2 improves Theorem 2.4 and Corollary 2.6. For the case of the extrinsic mean, let M be a d-dimensional differentiable manifold, and J : M-EN an embedding of M into an N-dimensional Euclidean space. Assume that J(M) is closed in EN, which is always the case, in particular, if M is compact. The extrinsic distance dE;J on M is defined as dE;J ðp; qÞ ¼ J JðpÞJðqÞ J for p; q A M, where J J denotes the Euclidean norm of EN. The image μ in J(M) of the extrinsic mean μE;J is then given by μ ¼ PðmÞ, where m is the usual mean of Q ○J 1 thought of as a probability on the Euclidean space EN, and P is the orthogonal projection defined on an N-dimensional neighborhood V of m into J(M) minimizing the Euclidean distance between m and J(M). If the projection P is unique on V then the projection μn ¼ Pðmn Þ of the Euclidean mean mn ¼ ∑nj¼ 1 JðY j Þ=n on J(M) is, with probability tending to one as n-1, unique and lies in an open neighborhood G of μ ¼ PðmÞ in J(M). Theorem 1 immediately implies the following result of Bhattacharya and Patrangenaru (2005) (Also see Bhattacharya and Bhattacharya, 2012, Proposition 4.3). Corollary 4 (CLT for Extrinsic Means on a Manifold). Assume that P is uniquely defined in a neighborhood of the N-dimensional Euclidean mean m of Q ○J 1 . Let ϕ be a diffeomorphism on a neighborhood G of μ ¼ PðmÞ in J(M) onto an open set U in Rd . Assume (A4)–(A6). Then, using the notation of (4), pffiffiffi pffiffiffi L n½ϕðμn ÞϕðμÞ ¼ n½ϕðPðmn ÞÞϕðPðmÞÞ-Nð0; Λ1 CΛ1 Þ;
as n-1:
R Remark 5. Note that Theorem 1 easily extends to Fréchet functions FðpÞ ¼ ρβ ðp; qÞQ ðdqÞ for arbitrary distances ρ and arbitrary β Z1, as long as the Assumptions (A1)–(A6) hold. See Bhattacharya and Bhattacharya (2012, Section 3), for questions concerning uniqueness of the Fréchet mean and consistency of its sample estimate.
R. Bhattacharya, V. Patrangenaru / Journal of Statistical Planning and Inference 145 (2014) 42–48
47
Example 6 (Kendall's shape space ∑km ; m 4 2Þ). We now turn to applications of Theorem 1 to the so-called stratified spaces S which are made up of several subspaces of different dimensions. One familiar example for the application is Kendall's shape space ∑km of k landmarks in Rm ðk 4mÞ. Each element of this space is the “shape” of a set of k points in Rm (not all the same), also called a k-ad, obtained modulo translation, scaling and rotation of the k-ad. After translation and scaling the k-ads lie in (and fill out) a preshape sphere Smkm1 . The shape space is then viewed as ∑km ¼ Smkm1 =SOðmÞ. For m 4 2 the shape space splits into strata of different dimensions (See Kendall et al., 1999, or Bhattacharya and Bhattacharya, 2012). For simplicity, consider the case m ¼3. One may split ∑k3 into two strata. The larger stratum S1 corresponds to shapes of non-collinear kads. The shape of a k-ad is then identified with the orbit of its preshape under rotation. Thus each orbit has dimension 3. If Q has a support contained in a geodesic ball in S1 of sufficiently small radius, then a unique Fréchet minimizer μ lying in S1 exists, and Theorem 1 applies with s ¼ d ¼ 3k7: Here we let μ be a “local” Fréchet mean, i.e., the minimization is only over the geodesic ball in which the support of Q lies (See Bhattacharya and Bhattacharya, 2008; Kendall, 1990). The stratum S0 may be given the structure of a differentiable manifold of dimension k2. A point in the preshape space of S0 may be labeled as ðλ; u1 ; …; uk1 Þ, where λ is a line in R3 through the origin (i.e., a point in RP 2 ), and the translated and scaled points of the kad lying on λ are expressed by their numerical positions along the line, so that ðu1 ; …; uk1 Þ belongs to Sk2 . The preshape space has therefore the dimension k. Since the orbit of a point of this space under SOð3Þ is of dimension 2, S0 has dimension k2. If Q assigns all its mass to S0, and the Fréchet minimizer is unique under the intrinsic distance of ∑3k , then the sample Fréchet mean belongs to S0 and Theorem 1 holds with s ¼ k2: If, however, Q ðS1 Þ 40, and the intrinsic mean μ exists as the unique minimizer, then, with probability tending to one as n-1, the sample Fréchet mean belongs to S1, and Theorem 1 holds with s ¼ d ¼ 3k7. The Theorem also applies to several examples where S is a space of non-positive curvature (NPC), which is not in general a differentiable manifold, but has a metric with properties of a geodesic distance (namely, minimum length of curves between points) and which is also somewhat analogous to differentiable manifolds of non-positive curvature. These spaces were originally studied by A.D. Alexandrov and developed further by Yu. G. Reshetnyak and M. Gromov (see Bridson and Häfliger, 1999; Sturm, 2003, for detailed treatments). Unlike differentiable manifolds of positive curvature where uniqueness of the intrinsic mean is known only under very restrictive conditions (see Karcher, 1977; Kendall, 1990; Afsari, 2011), on an NPC space the Fréchet mean is always unique, if the Fréchet function is finite (Sturm, 2003). An exposition of the CLT for stratified NPC spaces such as considered in Theorem 6.2 and, more generally, to the Open Book analyzed in Hotz et al. (2013), is given in Bhattacharya and Lin (2013). The Theorem is also applicable to other NPC spaces such as considered in Barden et al. (2013). References Afsari, B., 2011. Riemannian Lp center of mass: existence, uniqueness, and convexity. Proceedings of the American Mathematical Society 139, 655–673. Barden, D., Le, H., Owen, M., 2013. Central limit theorems for Fréchet means in the space of phylogenetic trees. Electronic Journal of Statistics, (to appear). Beran, R., 2013. Discussion: “Statistics on manifolds and landmarks based image analysis: a nonparametric theory with applications”, this issue. Bhattacharya, A., 2008. Statistical analysis on manifolds: a nonparametric approach for inference on shape spaces. Sankhya Series A 70 (2), 223–266. Bhattacharya, A., Bhattacharya, R., 2008. Statistics on Riemannian manifolds: asymptotic distribution and curvature. Proceedings of the American Mathematical Society 136, 2957–2967. Bhattacharya, A., Bhattacharya, R., 2012. Nonparametric Inference on Manifolds: With Applications to Shape Spaces. IMS Monograph Series, no. 2, Cambridge University Press. Bhattacharya, A., Dunson, D., 2010. Nonparametric Bayesian density estimation on manifolds with applications to planar shapes. Biometrika 97, 851–865. Bhattacharya, R., Lin, L., 2013. A Central Limit Theorem for Fréchet Means. 〈http://arxiv.org/abs/1306.5806〉. Bhattacharya, R.N., Patrangenaru, V., 2003. Large sample theory of intrinsic and extrinsic sample means on manifolds-I. Annals of Statistics 31, 1–29. Bhattacharya, R.N., Patrangenaru, V., 2005. Large sample theory of intrinsic and extrinsic sample means on manifolds-II. Annals of Statistics 33, 1225–1259. Bridson, M.R., Häfliger, A., 1999. Metric Spaces of Non-Positive Curvature. Grundlehren der mathematischen Wissenschaften, Springer. Ciccarelli, F., Doerks, T., von Mering, C., Creevey, C., Snel, B., Bork, P., 2006. Toward automatic reconstruction of a highly resolved tree of life. Science 311 (5765), 1283–1287. Dryden, I.L., Le, H., Preston, S., Wood, A., 2013. Discussion: “Statistics on manifolds and landmarks based image analysis: a nonparametric theory with applications”, this issue. Dryden, I.L., Mardia, K.V., 1998. Statistical Shape Analysis. Wiley, New York. Ellingson, L., Hendriks, H., Patrangenaru, V., Valentin, P., in press. CLT on low dimensional stratified spaces. In: Lahiri, S., Pispinis, D. (Eds.), Proceedings of the First INSPS Conference-Chalkidiki, Greece, 2012. Springer. Holmes, S., 2012. Bacterial trees in the human microbiome. 〈http://beta.mbi.ohio-state.edu/video?view=speakers&id=5111&item=Susan%20Holmes〉. Hotz, T., Huckemann, S., Le, H., Marron, J., Mattingly, J., Miller, E., Nolen, J., Owen, M., Patrangenaru, V., Skwerer, S., 2013. Sticky central limit theorems on open books. Advances in Applied Probability, in press. Huckemann, S., 2013. Discussion: “Statistics on manifolds and landmarks based image analysis: a nonparametric theory with applications”, this issue. Huckemann, S., Hotz, T., Munk, A., 2010a. Intrinsic MANOVA for Riemannian manifolds with an application to Kendall's spaces of planar shapes. IEEE Transactions on Pattern Analysis and Machine Intelligence 32 (4), 593–603. “Spotlight Paper” for this issue with its “Special Section on Shape Analysis and its Applications in Image Understanding”. Huckemann, S., Hotz, T., Munk, A., 2010b. Intrinsic shape analysis: geodesic PCA for Riemannian manifolds modulo isometric lie group actions (with discussions). Statistica Sinica 20, 1–100. Karcher, H., 1977. Riemannian center of mass and mollifier smoothing. Communications on Pure and Applied Mathematics 30, 509–554. Kendall, D., Barden, D., Carne, T., Le, H., 1999. Shape and Shape Theory. Wiley, New York. Kendall, W., 1990. Probability, convexity, and harmonic maps with small image I: uniqueness and fine existence. Proceedings of London Mathematical Society 61, 371–406. Kendall, W., Le, H., 2011. Limit theorems for empirical Fréchet means of independent and non-identically distributed manifold-valued random variables. Brazilian Journal of Probability and Statistics 25, 323–352.
48
R. Bhattacharya, V. Patrangenaru / Journal of Statistical Planning and Inference 145 (2014) 42–48
Kim, P., Pinder, S., Rush, S., 2013. Discussion: “Statistics on manifolds and landmarks based image analysis: a nonparametric theory with applications.” this issue. Mardia, K., Patrangenaru, V., 2005. Directions and projective shapes. Annals of Statistics 33, 1666–1699. Owen, M., Provan, S., 2011. A fast algorithm for computing geodesic distances in tree spaces. Computational Biology and Bioinformatics 8, 2–13. Patrangenaru, V., Liu, X., Sugathadasa, S., 2010. Nonparametric 3d projective shape estimation from pairs of 2d images. i. in memory of W.P. Dayawansa. Journal of Multivariate Analysis 101, 11–31. Pelletier, B., 2005. Kernel density estimation on Riemannian manifolds. Statistics and Probability Letters 73, 297–304. Sturm, K., 2003. Probability measures on metric spaces of nonpositive curvature. In: Auscher, P., Coulhon, T., Grigor'yan, A. (Eds.), Heat Kernels and Analysis on Manifolds, Graphs, and Metric Spaces, Contemporary Mathematics, vol. 338, AMS, 2003, pp. 357–390. von Bahr, B., 1967. On the central limit theorem on Rk. Arkiv för Matematik 7, 61–69. Ziezold, H., 1994. Mean figures and mean shapes applied to biological figure and shape distributions in the plane. Biometrical Journal 36, 491–510.
Contents lists available at ScienceDirect
Journal of Statistical Planning and Inference journal homepage: www.elsevier.com/locate/jspi
Rejoinder to the Discussion Rabi Bhattacharya a,n, Vic Patrangenaru b a b
Department of Mathematics, The University of Arizona, Tucson, AZ 85721, USA Department of Statistics, The Florida State University, Tallahassee, FL 32304, USA
The discussants have raised important issues on the subject matter of the present article. We are deeply appreciative of their valuable inputs and thank them for sharing their insights with us. As one of the early contributors to the theory of nonparametric statistics on manifolds, and as one whose research has made significant impacts on various aspects of statistics, Professor Beran is in a unique position to provide a broad perspective of the area under discussion. With regard to the need of developing appropriate analogs of multivariate analysis on manifolds that he rightly stresses, two contributions deserve particular attention. The first one by Huckemann et al. (2010b) has provided a theory of generalized principal component analysis (GPCA) on Riemannian manifolds and, more generally, on any quotient of a Riemannian manifold under the action of a Lie group of isometries. These cover, in particular, Kendall's shape spaces ∑km (for all m Z 2). Here the generalized principal components (GPC-s) are generalized geodesics. In particular, on a Riemannian manifold, the first GPC is the geodesic that has the smallest expected squared distance from a random variable X on the manifold with the given distribution Q. The second GPC is the geodesic with the smallest expected squared distance from X among all geodesics which intersect the first GPC orthogonally, and so on. Next, the article by Huckemann et al. (2010a) lays some groundwork for MANOVA on Riemannian manifolds. In the absence of the usual translation providing linearity or the additive structure of classical MANOVA on Euclidean spaces, on a Riemannian manifold M one employs parallel transport of vector fields along geodesics, made possible by the Levi-Civita connection. Specifically, in a two-factor model the r-th observation P ij;r at levels i and j of factors 1 and 2, respectively, is lifted to the tangent space T μij M by the log (inverse exponential) map, where μij is the intrinsic mean of the distribution of P ij;r . Because for different i, j the tangent spaces are not naturally related, log μij P ij;r is then transported to a common tangent space T ν M for some suitable ν (for example, a pooled intrinsic mean), by parallel translation along the geodesic joining μij and ν. The resulting vectors ϵij;r can play the role of the residuals, or i.i.d. errors, in the linear model. On functional inference mentioned by Professor Beran (2013), nonparametric density estimation, regression and classification have been carried out on manifolds in Bhattacharya and Dunson (2010) and Bhattacharya and Bhattacharya (2012), Chapters 13 and 14, via nonparametric Bayes theory. For density estimation, for example, this requires the imposition of a Dirichlet prior on the space of mixtures (probabilities) of a suitable parametric family of densities. The class of all densities obtained by such mixtures is required to be dense (say, in L1) in the class of all densities. In simulations and data examples, the Bayes posterior so obtained seems to perform at least as well as, and often better than, frequentist estimates by kernel methods such as that of Pelletier (2005). As to the two methods of bootstrapping mentioned by Professor Beran for computing or estimating p-values in twosample tests, while Bhattacharya and Patrangenaru (2005) use bootstrap confidence regions, both the confidence region based method and one based on simulations directly of the null distribution are considered in Bhattacharya and Bhattacharya (2012, Chapter 4). It may, however, be pointed out that in the context manifolds which are often high dimensional, and with samples of rather moderate sizes that are normally available, bootstrapping of the chisquare statistic runs into problems due to degeneracy of the bootstrapped covariance matrix. Also, even when one has non-degenerate bootstrap covariance, bootstrap methods do not seem particularly suitable to estimate the extremely small p-values such as obtained in Applications 3.3, 3.4 and 4.2. One may of course question the validity of the classical chisquare approximation of such small probabilities. But these are used here (and in Bhattacharya and Bhattacharya, 2012, Chapter 2) mostly to indicate
n
Corresponding author. E-mail address: [email protected] (R. Bhattacharya).
0378-3758/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jspi.2013.08.006
R. Bhattacharya, V. Patrangenaru / Journal of Statistical Planning and Inference 145 (2014) 42–48
43
extremely small orders of p-values (beyond the usual domain of bootstrap calculations) which may be justified by the order of approximation of tail probabilities provided by deep inequalities such as that of von Bahr (1967, Theorem 1). We share some of Professor Beran's reservations about statistical inference that is commonly practiced, including the impossibility of drawing a truly random number. By its very nature statistical inference is inductive and deals with highly illposed inverse problems, unlike most mathematics which is deductive. As the science underlying a physical phenomenon is understood better the level of uncertainty is reduced and inference becomes more accurate. This is happening in biology. However, as statistical mechanics teaches us, certain uncertainty is inherent in nature. As illustrated in the discussion, it is important to realize the limitations of statistical methods. However, even the commonly used methods, and many nonparametric procedures among them, will always remain indispensable in analyzing a host of phenomena. We next turn to the discussion by Dryden et al. (2013). Professors Dryden and Le are among the distinguished group of scholars including J. Kent, K. Mardia, C. Goodall and M. Prentice, who developed the theory of parametric inference on manifolds, following pioneering work by D.G. Kendall, F. Bookstein and H. Ziezold. Their recent work, some of which is surveyed in their discussion, has contributed significantly to nonparametric inference based on Fréchet type means. The present discussion by them provides a number of valuable insights concerning the notion of the mean on a manifold and the asymptotics of its empirical version. As to the question “In praise of which mean?” in Dryden et al. (2013), let us first recall that in Bhattacharya and Patrangenaru (2003, 2005) and Bhattacharya and Bhattacharya (2012), as surveyed here, a general theory is developed on the asymptotics of intrinsic and extrinsic Fréchet means and their applications to nonparametric inference on general abstract manifolds, not just for Kendall's shape spaces. For example, it applies to Stiefel manifolds and spaces of affine and projective shapes (Mardia and Patrangenaru, 2005; Patrangenaru et al., 2010; Bhattacharya and Bhattacharya, 2012, Chapters 9–11). It encompasses Fréchet minimizers of expected squared distances with respect to arbitrary distances ρ, subject of course to appropriate constraints. The topological constraint (2.1), namely, “closed bounded subsets are compact”, is identified as the appropriate one for proving consistency under uniqueness (Theorem 2.1). In addition, some smoothness of the Fréchet function and nondegeneracy of its averaged Hessian are required for deriving the CLT and tests and confidence regions based on it (see Theorems 2.4, 2.12, 2.15–2.17). As pointed out in Bhattacharya and Patrangenaru (2005) and in Remark 2.5 here, Theorem 2.4 applies to both intrinsic and extrinsic sample Fréchet means. As we will specify later, an extension of this theorem due to Bhattacharya and Lin (2013) removes some unnecessary constraints and yields a general CLT for both intrinsic and extrinsic Fréchet means extending Corollary 2.6 and Proposition 2.10 on manifolds, and it also applies to several important examples of stratified spaces. It has been pointed out in Bhattacharya and Bhattacharya (2012), Sections 3.2, 8.5, that, under uniqueness, consistency holds with the squared distance ρ2 replaced by ρβ for β Z1 in the definition of the Fréchet function, and the CLT for the sample Fréchet mean may also hold for β 4 2. Why is then the case β ¼ 2 special? First, for extrinsic means as defined by us, precise (necessary and sufficient) conditions for uniqueness of the Fréchet minimizer may often be obtained, and an explicit analytical formula for it be found. This is not the case for β a 2. In particular, this clarifies the issue of consistency, or lack of it, of Procrustes type estimators (“partial”, “full”, “generalized”) for what was defined as the population mean shape in the earlier literature (see, e.g., Dryden and Mardia, 1998, p. 88), referred to as “the shape of the mean” by the present discussants. While a Procrustes estimator may be an extrinsic empirical (or, sample-) Fréchet mean, the quantity it was supposed to estimate, namely the mean shape as defined in the past literature, is not the corresponding Fréchet mean, leading to inconsistency unless some extraneous restriction such as isotropy is imposed. Also, for the intrinsic mean, for β o 2, especially for β ¼ 1, the uniqueness condition has been proven only under more restrictive conditions than for β ¼ 2 (See Afsari, 2011). Secondly, it is not clear that the smoothness condition required for the CLT holds in general for β o2, even when uniqueness of the minimizer is assured. Also, it is worth emphasizing again what is stated in our Introduction, namely, the embeddings that we use are equivariant under the action of a large group preserving a great deal of the geometry and this is undoubtedly a reason for the essential equivalence between the extrinsic and intrinsic inferences displayed in the Applications in the article (also see Chapter 2 of Bhattacharya and Bhattacharya, 2012). The difference between the two sample means and of estimates and test statistics (and p-values) based on them is of smaller magnitude by several orders than the average difference between the extrinsic and intrinsic distances computed for pairs of sample observations. This shows, first of all, that one may think of extrinsic inference as a computationally convenient very good approximation to intrinsic inference, provided the embedding for the extrinsic analysis is equivariant under a large group H, perhaps especially when H is a group of isometries on the manifold if endowed with the metric tensor induced by the embedding. Secondly, it provides some confidence in the validity of using the intrinsic mean in most real data examples, without being bothered by the vexing theoretical problem of uniqueness. On reflection-similarity shape spaces, under the MDS or Schoenberg embedding (see (4.2) in the present article), Dryden et al. (2013) provide significant results on the uniqueness and, consequently, consistency of empirical Fréchet means under a family of distances dα ðα Z 1=2Þ, with d1 being the induced Euclidean distance used in our definition of the extrinsic mean under the embedding. They have also provided explicit formulas for the corresponding Fréchet means. The case α ¼ 1 was solved by Bhattacharya (2008), who also provided the CLT for the empirical Fréchet mean for this case. It should be borne in mind that such a family of distances, based on a spectral decomposition, is not generally available for arbitrary embeddings of other manifolds. We now turn to the last section of the discussion by Dryden et al. (2013) on the limiting distribution of the intrinsic sample mean. It may be noted that while a general CLT for intrinsic sample means was derived in Bhattacharya and Patrangenaru (2005), a study linking the curvature of a Riemannian manifold and the dispersion of the limiting Normal
44
R. Bhattacharya, V. Patrangenaru / Journal of Statistical Planning and Inference 145 (2014) 42–48
distribution in the CLT appeared in Bhattacharya and Bhattacharya (2008). In the latter article an explicit expression of the dispersion matrix in terms of curvature is given for manifolds with constant sectional curvature, and a lower bound in terms of sectional curvature for arbitrary Riemannian manifolds. These results are extended in a fine article by Kendall and Le (2011), and they are stated in the discussion by Dryden et al. (2013). In particular, dispersion is expressed in terms of curvature for Kahler (or, complex-) manifolds with constant holomorphic curvature as well, with an explicit expression in the case of Kendall's planar shape space. A different form of this dispersion also appears in Bhattacharya and Bhattacharya (2008, Theorem 3.1) and remarks following it, but without linking this to curvature. As emphasized in the Introduction and elsewhere in our work, the Fréchet mean is used as a nonparametric tool, not so much for inference on the mean itself as it is for discrimination among distributions on manifolds, especially in two-sample problems. In exceptional cases, as in the case of directional statistics, the mean may be the main point of interest. Estimating the true direction of the earths magnetic poles in past geological time scales is an example, where the remanent fossil magnetism provides data with the true value of the direction corrupted by random errors or, what one may call, measuremental errors (Application 2.14). For testing solely via empirical Fréchet means if two random samples on a manifold are from different distributions, there is of course the usual theoretical problem of identification (namely, distributions with the same mean may be different). But in most cases of interest, such as shape spaces, the manifold is of high dimension and the Fréchet mean is shown to discriminate different distributions in two-sample problems quite well (Applications 3.3, 3.4, 4.2), often far better than parametric tests designed to test directly if the distributions are different, indicating perhaps that the parametric tests so carried out are highly misspecified. Because of the well known curse of dimensionality, in these cases the effectiveness of more elaborate nonparametric functional inference for discrimination requires much larger sample sizes than are generally available. Hence our phrase “in praise of the mean”. We do recognize the usefulness of any so-called mean which carries out this task of discrimination well. We now turn to the discussion by Huckemann (2013). Professor Huckeman's important contributions to multivariate analysis are mentioned at the beginning of the Rejoinder. In the discussion he proposes to use the Ziezold mean, following Ziezold (1994), which is defined on a quotient M ¼ N\G, where N is a Riemannian manifold embedded in a Euclidean space, and G is a Lie group of isometries on N which is also assumed to be an isometry on N with respect to the Euclidean distance on it. The Ziezold mean is the Fréchet mean with respect to the distance ρð½x; ½yÞ ¼ inf f J gx′y′ J : g A Gg, where x′; y′ are arbitrary elements in N of the orbits ½x, ½y A M, respectively, and J J is the Euclidean norm. If the action is free and proper then M is a manifold and one can define a Riemannian metric on it by Riemannian submersion of N into M, and the intrinsic distance would be different from the Ziezold distance. In general, M need not be a manifold. For example, M can be ∑km for any m Z 2. On ∑km the Ziezold distance is also known as the partial Procrustes distance (See Dryden and Mardia, 1998, p. 65). The comparisons between the Ziezold mean and the extrinsic mean for the 3D reflection shape by data analysis and simulation in Examples 3.1 and 3.2 are valid and indicate that if the underlying distributions assign positive probability near, or on, the singular part of the space then the extrinsic mean based on Schoenberg embedding is not going to be very effective as a tool for discrimination. We expect that the intrinsic mean on the whole space ∑k3 , however, will be at least as effective as the Ziezold mean in such cases. There remains the question of uniqueness of the Ziezold mean. The expectation that uniqueness of the Fréchet minimizer holds more broadly for this distance than for the intrinsic distance may be valid. But it needs to be explored. Professor Huckemann is quite right in emphasizing that the uniqueness problem is one of the main issues confronting the nonparametric theory using Fréchet means. It is here that the extrinsic mean holds an edge over other means. Still one expects uniqueness for both the intrinsic and Ziezold means to hold outside a relatively small set of distributions. On the discussion by Professor Kim et al. The recently National Institutes of Health sponsored Human Microbiome Project (HMP), aims to “characterize microbial communities found at several different sites on the human body, including nasal passages, oral cavities, skin, gastrointestinal tract, and urogenital tract, and to analyze the role of these microbes in human health and disease”. From a general statistician perspective, the HMP is directed toward understanding correlations between changes in the microbiome and human health. From a biologistʼs perspective, the HMP may be regarded as one small project of the larger prokaryote–eukaryotes relationship, and from a higher general scientific view, this relationship is more meaningful in presence of an understanding of and closeness to Mother Nature. Recall that viruses are the tiniest life forms known to date and they replicate inside living cells; archaea are unicellular microorganisms, microbes that have no cell nucleus or any other membrane-bound organelles within their cells, and bacteria are all other nucleus free cellular organisms (prokaryote). The organisms whose cells do have a nucleus are called eukaryotes. However as noted in the discussion by Kim et al. (2013), sequencing and classification are based on DNA structure and not on cellular organization of an organism. Classifying life forms is the focus of taxonomy, and the DNA tags attached to an OTU are commonly used in such classifications in recent times (see e.g. Ciccarelli et al., 2006, for the left hand side of Fig. 1). The suggestion made by Kim et al. (2013), of sketching possible directions that have not appeared in the literature for analyzing microbiome data, parallels similar projects of using phylogenetic analysis of DNA sequencing mostly associated with computational methods in Biology to study the microbiome; these include Holmes (2012), where in particular part of the human skin microbiome was discussed. In Fig. 1 right hand side, the human skin microbiome with types and proportions of prokaryote at various parts is schematically displayed. Since proportions and mean types of bacteria play a key role in general in the health state of eukaryotes, the suggestion of running a Fréchet mean analysis for microbiome data is welcome. As noted in the discussion by Dryden et al. (2013), the Fréchet mean depends on the distance used on the sample space. This is in particular true in the case of operational taxonomic units (OTUs) that are represented as leaves on phylogenetic trees, points on the tree spaces Tp.
R. Bhattacharya, V. Patrangenaru / Journal of Statistical Planning and Inference 145 (2014) 42–48
45
Fig. 1. Left: An oversimplified color-coded version of the “tree of life”: prokaryote (blue), archaea (green) and eukaryotes (red). Right: Human Skin Microbiome. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
Fig. 2. Left: T4, space of trees with four leaves. Right: Petersen graph.
We described in detail the tree spaces, and commented on their stratification. In the low dimensional cases p r4, the topological structure is well known as presented in our paper. For example T4 can be representation as a surface with singularities obtained from the polyhedral surface given on the right hand side of Fig. 2 by identifying the edges labeled with the same letter and then embedded in R4 : The intersection of the 15 quadrants with the unit L1 in R4 is given by the 15 edges of the Petersen graph (see left side of Fig. 2). This graph provides a convenient representation for data on T4. The Fréchet mean of a sample of size 10 of simulated phylogenetic trees obtained with respect to this embedding for the concatenated eukaryote Parkinsea RNA tree with four leafs data was computed in Ellingson et al. (2012). There they show the resulting Fréchet sample mean tree and its relative location in the tree space in terms of the Petersen graph. Computational algorithms for intrinsic sample means on Tp were given by Owen and Provan (2011). Recently Ellingson et al. (2012) also considered some other computational examples to illustrate the behavior of extrinsic means for simulated data on T4. For the asymptotic distribution of the intrinsic sample means on this space, we refer to Barden et al. (2013). These early results along the line of the Fréchet mean analysis on phylogenetic tree spaces suggested by Professor Kim et al. (2013) are encouraging for more computationally intensive examples arising in the prokaryote–eukaryote relationship. In conclusion we state a recent general CLT for empirical or sample Fréchet means on manifolds and stratified spaces, since asymptotic distribution theory is a major theme underlying the discussions. Let ðS; ρÞ be a metric space and Q a probability measure on its Borel sigma-field.
46
R. Bhattacharya, V. Patrangenaru / Journal of Statistical Planning and Inference 145 (2014) 42–48
Assume that the Fréchet function F under squared distance is finite on S, and that F has a unique minimizer μ ¼ argminp FðpÞ—the Fréchet mean of Q (with respect to the distance ρ). Under broad conditions, the Fréchet sample mean μn of the empirical distribution Q n ¼ 1=n∑nj¼ 1 δY j based on independent S-valued random variables Yj ðj ¼ 1; …; nÞ with common distribution Q is a consistent estimator of μ. That is, μn -μ almost surely, as n-1. Here μn may be taken to be any measurable selection from the (random) set of minimizers of the Fréchet function of Qn, namely F n ðpÞ ¼ ð1=nÞ∑nj¼ 1 ρ2 ðp; Y j Þ. We make the following assumptions:
(A1) The Fréchet mean μ of Q is unique. (A2) μA G, where G is a measurable subset of S, and there is a homeomorphism ϕ : G-U, where U is an open subset of Rs for some s Z1 and G is given its relative topology on S. Also, x↦hðx; qÞ≔ρ2 ðϕ1 ðxÞ; qÞ
ð1Þ
is twice continuously differentiable on U, for every q outside a Q-null set. (A3) μn A G almost surely, for all sufficiently large n. (A4) Let Dr hðx; qÞ ¼ ∂hðx; qÞ=∂xr ; r ¼ 1; …; s. Then EjDr hðϕðμÞ; Y 1 Þj2 o 1; EjDr;r′ hðϕðμÞ; Y 1 Þj o1
for r; r′ ¼ 1; …; s:
ð2Þ
(A5) Let ur;r′ ðϵ; qÞ ¼ supfjDr;r′ hðθ; qÞDr;r′ hðϕðμÞ; qÞj : jθϕðμÞj oϵg. Then Ejur;r′ ðϵ; Y 1 Þj-0
as ϵ-0
for all 1 r r; r′ r s:
ð3Þ
(A6) The matrix Λ ¼ ½EDr;r′ hðϕðμÞ; Y 1 Þr;r′ ¼ 1;…;s is nonsingular. Theorem 1. Under Assumptions (A1)–(A6), L
n1=2 ½ϕðμn ÞϕðμÞ-Nð0; Λ1 CΛ1 Þ;
ð4Þ
as n-1;
where C is the covariance matrix of fDr hðϕðμÞ; Y 1 Þ; r ¼ 1; …; sg. This theorem extends and improves the CLT-s given in Bhattacharya and Patrangenaru (2005) and Bhattacharya and Bhattacharya (2012), Chapters 3 and 4. It is derived in Bhattacharya and Lin (2013). It has the following immediate consequences. Corollary 2 (CLT for Intrinsic Means on a Manifold). Let (M,g) be a d-dimensional complete Riemannian manifold with metric tensor g and geodesic distance ρg . Suppose Q is a probability measure on M with intrinsic mean μI , which is the unique minimizer of the Fréchet function for the distance ρg . Let ϕ ¼ Exp μ1 be the inverse exponential, or log-, function at μI defined on a I neighborhood G of μ ¼ μI onto its image U in the tangent space T μI ðMÞ. Assume that the Assumptions (A4)–(A6) hold. Then, with s ¼d, the CLT (4) holds for the intrinsic sample mean μn ¼ μn;I . Remark 3. Corollary 2 improves Theorem 2.4 and Corollary 2.6. For the case of the extrinsic mean, let M be a d-dimensional differentiable manifold, and J : M-EN an embedding of M into an N-dimensional Euclidean space. Assume that J(M) is closed in EN, which is always the case, in particular, if M is compact. The extrinsic distance dE;J on M is defined as dE;J ðp; qÞ ¼ J JðpÞJðqÞ J for p; q A M, where J J denotes the Euclidean norm of EN. The image μ in J(M) of the extrinsic mean μE;J is then given by μ ¼ PðmÞ, where m is the usual mean of Q ○J 1 thought of as a probability on the Euclidean space EN, and P is the orthogonal projection defined on an N-dimensional neighborhood V of m into J(M) minimizing the Euclidean distance between m and J(M). If the projection P is unique on V then the projection μn ¼ Pðmn Þ of the Euclidean mean mn ¼ ∑nj¼ 1 JðY j Þ=n on J(M) is, with probability tending to one as n-1, unique and lies in an open neighborhood G of μ ¼ PðmÞ in J(M). Theorem 1 immediately implies the following result of Bhattacharya and Patrangenaru (2005) (Also see Bhattacharya and Bhattacharya, 2012, Proposition 4.3). Corollary 4 (CLT for Extrinsic Means on a Manifold). Assume that P is uniquely defined in a neighborhood of the N-dimensional Euclidean mean m of Q ○J 1 . Let ϕ be a diffeomorphism on a neighborhood G of μ ¼ PðmÞ in J(M) onto an open set U in Rd . Assume (A4)–(A6). Then, using the notation of (4), pffiffiffi pffiffiffi L n½ϕðμn ÞϕðμÞ ¼ n½ϕðPðmn ÞÞϕðPðmÞÞ-Nð0; Λ1 CΛ1 Þ;
as n-1:
R Remark 5. Note that Theorem 1 easily extends to Fréchet functions FðpÞ ¼ ρβ ðp; qÞQ ðdqÞ for arbitrary distances ρ and arbitrary β Z1, as long as the Assumptions (A1)–(A6) hold. See Bhattacharya and Bhattacharya (2012, Section 3), for questions concerning uniqueness of the Fréchet mean and consistency of its sample estimate.
R. Bhattacharya, V. Patrangenaru / Journal of Statistical Planning and Inference 145 (2014) 42–48
47
Example 6 (Kendall's shape space ∑km ; m 4 2Þ). We now turn to applications of Theorem 1 to the so-called stratified spaces S which are made up of several subspaces of different dimensions. One familiar example for the application is Kendall's shape space ∑km of k landmarks in Rm ðk 4mÞ. Each element of this space is the “shape” of a set of k points in Rm (not all the same), also called a k-ad, obtained modulo translation, scaling and rotation of the k-ad. After translation and scaling the k-ads lie in (and fill out) a preshape sphere Smkm1 . The shape space is then viewed as ∑km ¼ Smkm1 =SOðmÞ. For m 4 2 the shape space splits into strata of different dimensions (See Kendall et al., 1999, or Bhattacharya and Bhattacharya, 2012). For simplicity, consider the case m ¼3. One may split ∑k3 into two strata. The larger stratum S1 corresponds to shapes of non-collinear kads. The shape of a k-ad is then identified with the orbit of its preshape under rotation. Thus each orbit has dimension 3. If Q has a support contained in a geodesic ball in S1 of sufficiently small radius, then a unique Fréchet minimizer μ lying in S1 exists, and Theorem 1 applies with s ¼ d ¼ 3k7: Here we let μ be a “local” Fréchet mean, i.e., the minimization is only over the geodesic ball in which the support of Q lies (See Bhattacharya and Bhattacharya, 2008; Kendall, 1990). The stratum S0 may be given the structure of a differentiable manifold of dimension k2. A point in the preshape space of S0 may be labeled as ðλ; u1 ; …; uk1 Þ, where λ is a line in R3 through the origin (i.e., a point in RP 2 ), and the translated and scaled points of the kad lying on λ are expressed by their numerical positions along the line, so that ðu1 ; …; uk1 Þ belongs to Sk2 . The preshape space has therefore the dimension k. Since the orbit of a point of this space under SOð3Þ is of dimension 2, S0 has dimension k2. If Q assigns all its mass to S0, and the Fréchet minimizer is unique under the intrinsic distance of ∑3k , then the sample Fréchet mean belongs to S0 and Theorem 1 holds with s ¼ k2: If, however, Q ðS1 Þ 40, and the intrinsic mean μ exists as the unique minimizer, then, with probability tending to one as n-1, the sample Fréchet mean belongs to S1, and Theorem 1 holds with s ¼ d ¼ 3k7. The Theorem also applies to several examples where S is a space of non-positive curvature (NPC), which is not in general a differentiable manifold, but has a metric with properties of a geodesic distance (namely, minimum length of curves between points) and which is also somewhat analogous to differentiable manifolds of non-positive curvature. These spaces were originally studied by A.D. Alexandrov and developed further by Yu. G. Reshetnyak and M. Gromov (see Bridson and Häfliger, 1999; Sturm, 2003, for detailed treatments). Unlike differentiable manifolds of positive curvature where uniqueness of the intrinsic mean is known only under very restrictive conditions (see Karcher, 1977; Kendall, 1990; Afsari, 2011), on an NPC space the Fréchet mean is always unique, if the Fréchet function is finite (Sturm, 2003). An exposition of the CLT for stratified NPC spaces such as considered in Theorem 6.2 and, more generally, to the Open Book analyzed in Hotz et al. (2013), is given in Bhattacharya and Lin (2013). The Theorem is also applicable to other NPC spaces such as considered in Barden et al. (2013). References Afsari, B., 2011. Riemannian Lp center of mass: existence, uniqueness, and convexity. Proceedings of the American Mathematical Society 139, 655–673. Barden, D., Le, H., Owen, M., 2013. Central limit theorems for Fréchet means in the space of phylogenetic trees. Electronic Journal of Statistics, (to appear). Beran, R., 2013. Discussion: “Statistics on manifolds and landmarks based image analysis: a nonparametric theory with applications”, this issue. Bhattacharya, A., 2008. Statistical analysis on manifolds: a nonparametric approach for inference on shape spaces. Sankhya Series A 70 (2), 223–266. Bhattacharya, A., Bhattacharya, R., 2008. Statistics on Riemannian manifolds: asymptotic distribution and curvature. Proceedings of the American Mathematical Society 136, 2957–2967. Bhattacharya, A., Bhattacharya, R., 2012. Nonparametric Inference on Manifolds: With Applications to Shape Spaces. IMS Monograph Series, no. 2, Cambridge University Press. Bhattacharya, A., Dunson, D., 2010. Nonparametric Bayesian density estimation on manifolds with applications to planar shapes. Biometrika 97, 851–865. Bhattacharya, R., Lin, L., 2013. A Central Limit Theorem for Fréchet Means. 〈http://arxiv.org/abs/1306.5806〉. Bhattacharya, R.N., Patrangenaru, V., 2003. Large sample theory of intrinsic and extrinsic sample means on manifolds-I. Annals of Statistics 31, 1–29. Bhattacharya, R.N., Patrangenaru, V., 2005. Large sample theory of intrinsic and extrinsic sample means on manifolds-II. Annals of Statistics 33, 1225–1259. Bridson, M.R., Häfliger, A., 1999. Metric Spaces of Non-Positive Curvature. Grundlehren der mathematischen Wissenschaften, Springer. Ciccarelli, F., Doerks, T., von Mering, C., Creevey, C., Snel, B., Bork, P., 2006. Toward automatic reconstruction of a highly resolved tree of life. Science 311 (5765), 1283–1287. Dryden, I.L., Le, H., Preston, S., Wood, A., 2013. Discussion: “Statistics on manifolds and landmarks based image analysis: a nonparametric theory with applications”, this issue. Dryden, I.L., Mardia, K.V., 1998. Statistical Shape Analysis. Wiley, New York. Ellingson, L., Hendriks, H., Patrangenaru, V., Valentin, P., in press. CLT on low dimensional stratified spaces. In: Lahiri, S., Pispinis, D. (Eds.), Proceedings of the First INSPS Conference-Chalkidiki, Greece, 2012. Springer. Holmes, S., 2012. Bacterial trees in the human microbiome. 〈http://beta.mbi.ohio-state.edu/video?view=speakers&id=5111&item=Susan%20Holmes〉. Hotz, T., Huckemann, S., Le, H., Marron, J., Mattingly, J., Miller, E., Nolen, J., Owen, M., Patrangenaru, V., Skwerer, S., 2013. Sticky central limit theorems on open books. Advances in Applied Probability, in press. Huckemann, S., 2013. Discussion: “Statistics on manifolds and landmarks based image analysis: a nonparametric theory with applications”, this issue. Huckemann, S., Hotz, T., Munk, A., 2010a. Intrinsic MANOVA for Riemannian manifolds with an application to Kendall's spaces of planar shapes. IEEE Transactions on Pattern Analysis and Machine Intelligence 32 (4), 593–603. “Spotlight Paper” for this issue with its “Special Section on Shape Analysis and its Applications in Image Understanding”. Huckemann, S., Hotz, T., Munk, A., 2010b. Intrinsic shape analysis: geodesic PCA for Riemannian manifolds modulo isometric lie group actions (with discussions). Statistica Sinica 20, 1–100. Karcher, H., 1977. Riemannian center of mass and mollifier smoothing. Communications on Pure and Applied Mathematics 30, 509–554. Kendall, D., Barden, D., Carne, T., Le, H., 1999. Shape and Shape Theory. Wiley, New York. Kendall, W., 1990. Probability, convexity, and harmonic maps with small image I: uniqueness and fine existence. Proceedings of London Mathematical Society 61, 371–406. Kendall, W., Le, H., 2011. Limit theorems for empirical Fréchet means of independent and non-identically distributed manifold-valued random variables. Brazilian Journal of Probability and Statistics 25, 323–352.
48
R. Bhattacharya, V. Patrangenaru / Journal of Statistical Planning and Inference 145 (2014) 42–48
Kim, P., Pinder, S., Rush, S., 2013. Discussion: “Statistics on manifolds and landmarks based image analysis: a nonparametric theory with applications.” this issue. Mardia, K., Patrangenaru, V., 2005. Directions and projective shapes. Annals of Statistics 33, 1666–1699. Owen, M., Provan, S., 2011. A fast algorithm for computing geodesic distances in tree spaces. Computational Biology and Bioinformatics 8, 2–13. Patrangenaru, V., Liu, X., Sugathadasa, S., 2010. Nonparametric 3d projective shape estimation from pairs of 2d images. i. in memory of W.P. Dayawansa. Journal of Multivariate Analysis 101, 11–31. Pelletier, B., 2005. Kernel density estimation on Riemannian manifolds. Statistics and Probability Letters 73, 297–304. Sturm, K., 2003. Probability measures on metric spaces of nonpositive curvature. In: Auscher, P., Coulhon, T., Grigor'yan, A. (Eds.), Heat Kernels and Analysis on Manifolds, Graphs, and Metric Spaces, Contemporary Mathematics, vol. 338, AMS, 2003, pp. 357–390. von Bahr, B., 1967. On the central limit theorem on Rk. Arkiv för Matematik 7, 61–69. Ziezold, H., 1994. Mean figures and mean shapes applied to biological figure and shape distributions in the plane. Biometrical Journal 36, 491–510.