VERSITA TaTra AOuNTaiNS Mathematical Publications DOI: 10.2478/vl0127-012-0009-9 Tatra Mt. Math. Publ. 51 (2012), 83-90 DETECTION OF NON-AFFINE SHAPE OUTLIERS FOR MATCHED-PAIR SHAPE DATA Stanislav Katina ABSTRACT. Cleft lip/palate (CLP) is a relatively common birth defect so disfiguring that nowadays it is almost always corrected surgically as early as possible. The postnatal surgical correction does not, however, result in a normally growing upper jaw, but instead, owing to scar tissue, one that grows abnormally. It is important to decide if a clinical treatment group is homogeneous. The example involves data from digitally processed lateral X-ray films of 48 boys who have complete unilateral CLP but no other malformation. 22 landmarks were represented by their Procrustes shape coordinates, principal components of matched-pair differences were examined, and the distribution of the 48 shape changes was studied for outliers in the affine and non-affine subspaces of the full Procrustes shape and form space. To separate outliers from inliers we use bagplots. There are no outliers apparent in the affine subspace. In the non-affine subspaces, we found no outliers in the subspace of bending patterns at large scale but four outliers in the subspace of local changes at small scale. Almost the same outliers were found by form-space PC A. These latter are associated with possible creases of the corresponding thin-plate splines. In those cases we can use the same spline formalism to relax the outlying form to an inlier by optimal relaxation along the curve décolletage that weighs bending energy against Procrustes distance and stop relaxation on the fence. These maneuvers suggest a possibly novel and interesting fusion of the Procrustes-spline toolkit with outlier detection. They also have practical implications for craniofacial management of CLP follow-up as well as suggestive implications for outlier detection in applied craniometries and anthropometrics more generally. 1. Introduction Cleft lip/palate (CLP) is a relatively common birth defect (about 1 per 1000) so disfiguring that nowadays it is almost always corrected surgically as early © 2012 Mathematical Institute, Slovak Academy of Sciences. 2010 Mathematics Subject Classification: Primary 62H25, 62H35, 62P10; Secondary 6209, 62H11, 62G35. Keywords: shape outliers, relative warp analysis, matched-pair shape differences, bagplots, thin-plate splines, cleft lip/palate. This research was supported by VEGA [2/0038/12] and Welcome Trust [WT086901MA]. 83 STANISLAV KATI NA as possible. The postnatal surgical correction does not, however, result in a normally growing maxilla (upper jaw), but instead, owing to scar tissue, one that grows abnormally. For the clinical biostatistician it is important to decide if a clinical treatment group is homogeneous. In this paper we comment on this question both as an example of geometric morphometries (GMM) and as an example of the kind of question that craniofacial surgeons might propose in their role as physical anthropologists, students of the variation of human form. 2. Data The example involves data from digitally processed lateral X-ray films (under standard conditions—focus-object distance is 370 cm, object-film distance is 30 cm, magnification is 8.1%) of 48 boys who have complete unilateral CLP but no other malformation, born between the years 1985-1988. All were surgically corrected in infancy by the same team of surgeons at the Clinic of Plastic Surgery in Prague, Czech Republic. Primary cheiloplasty performed according to Tennison's method at an average age of 9 month was associated with periosteo-plasty (without nasal septum repositioning) using a 5-7mm wide and 15-20mm long periostal flap obtained from the lateral maxillary segment [8]. All patients underwent palatoplasty at an average age of 5 years, always by method of push-back with pharyngeal flap surgery. The corrective procedures were implemented in the patients with persisting soft tissue deformations like lip and nose, lip or nose [9]. Data are from follow-up at two ages, 10 and 15 years (measured at both times), within the adolescent growth spurt. 3. Detection and relaxation of shape outliers 22 cephalometric points (landmarks) were represented by their Procrustes shape coordinates [1], principal components of matched-pair differences were examined, and the distribution of the 48 shape changes was studied for outliers in the affine and the non-affine subspace of the full Procrustes shape and form space [3]. To separate outliers from inliers we use bagplots (bivariate boxplots) as supplied by R software [5]. Along curve decolletage, we relax the outliers by flattening the surface in the direction to Procrustes mean shape. Relative warp analysis (RWA) is modified PC A [1], where relative warps (RWs) are PCs with respect to the bending energy or inverse bending energy matrix. Let xp,i = Vec (Xp;i), i = 1,2, ...,n; n = 48 be a 2fc-vector of Procrustes shape coordinates, where Xp - - fx(1) : x(2)>l x(d) - (x{d) x{d) x(d)) d-12 84 DETECTION OF NON-AFFINE SHAPE OUTLIERS FOR MATCHED-PAIR SHAPE DATA Let ~x.d,i be 2fc-vectors (k = 22) of matched-pair differences of vectorized Procrustes shape coordinates, xP,i5,i = Vec(Xp)i5)i) and xp;i0,i = Vec(Xp)io,i), Sp, be the covariance matrix of the data xp,;i, Xp,i0 = (4!!o : *pio) = (xi, • • •, xfc)T be k x 2 matrix of mean Procrustes shape coordinates x7- of 10-year group, j = 1, 2,..., k. Let L / S lfc Xp.io \ , l11 l12 If 0 0 V xL oo/ T —1 _ / fcxfc fcx3 _ ' t 21 t 22 -^Xfe -^3x3 "P,10 where L is symmetric positive definite, the inverse of S exists as long as the landmarks are at least four in number, not all on one straight line, and also not in the same place (coincident); then inverse of L exists and is equal to L_1, SJS = 0(x"j -xfl); j, s = 1, 2,..., k, (x) = ||x||2 log(||x||2), V||x||2 > O^if ||x||2 = 0, (x) = 0; kxk matrix Be = L11 is called bending energy matrix ofXp^o, 2kx2k matrix B = 12x2 Be. There are these constrains of this matrix l^Bg = 0, XTBe = 0, so the rank of the bending energy matrix is k — 3. Then the non-zero eigenvalues of (B~)a^2 Sp, (B~)a^2 are lj with corresponding eigenvectors gj (PC loadings, RWs), (B~)a^2 = A- ^^jJjj is Moore-Penrose generalized Q; /2 inverse of Ba/2 RW scores are defined as r^- = gj (B ) x^^. The effect of the jth RW can be viewed by plotting Vec(Xp(Cj,j,a)) =Vec(Xp,1o)±cJBa/2g^1/2, rj=c$/2, c3 e R+ (1) for various values of rj G (0, max(|rij|)) (or some magnification of max(|r£j|); ^i ii alternatively, fixing Cj = 1, magnification of / ■ as standard deviation of PCj scores). The effect of the linear combination of RWi and RW2 can be viewed by plotting Vec(XP (ci, c2, a)) = Vec(XP,10) ± c^"/2^^2 ± ^V*/2^2. A PC summary of the shape data 2 2 Vec(Xp,15,, (a)) = Vec(XP,io) ± B^2 ^rljgj = Vec(XP,io) + J]gJ-gJxD)i. 3=1 3=1 If a = 1, bending patterns are at large scale, if a = —1, bending patterns are at small scale, and if a = 0, then we take B° = I as the 2k x 2k identity matrix and the procedure is exactly the same as classical PCA of matched-pair differences. 85 STANISLAV KATI NA Any value of a other than zero converts the analysis to RWA of shape covariances with respect to bending energy instead of Procrustes distance. To find the affine component we use linear regression model (d) (d) (d) o{d) id) d-12- i-12 n then is the affine component of Xp^. In affine subspace, Sp,^ stands for sample co-variance matrix of ~x.da,i = Vec(XJ4;i) — Vec(Xp;i0). Then PCA of Sda is called affine-subspace PCA. Let Xqj? = (Xp,:xsize), be an n x (2fc + 1) matrix with the rows equal to xdf,i=(^b,iM{CSi)) , i = l,2,...,n, and xLe=(ln(C5i),...,ln(C5n)). Let Spip be the covariance matrix of the data ~x.df,i- Then PCA of Sp,p is called form-space PCA. The first PC represents allometry—shape change during growth. Visualization of interpolated shape changes can be done via thin-plate spline (TPS) deformation grids [1], field of vectors (within the convex hull of reference shape Xpio, where longer vectors show stronger deformation in the specific direction of the shape change) superimposed with the grid of gray-scale rectangles with colors corresponding to the Procrustes distances (regions showing milder deformation are lighter, regions with stronger deformation are darker; the surface does not show the direction—but only the size—of some shape change). Let x = Vec (X) be a 2fc-vector of vectorized Procrustes shape coordinates. To relax a configuration x onto the mean configuration xm, seek the configuration xr that minimizes the weighted sums of two energies. One term is the bending energy (xr — xm) B (xr — xm) of xr considered as a deformation of the mean form. The other term (xr — x) (xr — x) is the sum of squares that corresponds to squared Procrustes distance between the relaxed landmark configuration xr and the vector x actually encountered. For any regularization parameter (relative weight) A, we seek the relaxed configuration xr that minimizes the weighted sum (xr — x) (xr — x) + A (xr — xm) B (xr — xm) of the energies. Setting the gradient of this expression to zero straightforwardly generates the solution xr = xr (A) = (I2fc + AB)"1 (x + ABxm). (2) As A varies, xr (A) traverses a smooth curve in the space of landmark configurations, named curve decolletage [2]. Note that when A = 0 then xr (0) = x, 86 DETECTION OF NON-AFFINE SHAPE OUTLIERS FOR MATCHED-PAIR SHAPE DATA in the absence of any penalty for deviation from the mean configuration xm, the best fit to any data is the selfsame data. In the other limit A —> oo, xr tends to the configuration of zero bending having the least Euclidean cost (ordinary least square fit to the data, in the affine subspace). Because B is a function of xm, this formulation is not symmetric in x and xm. By minimizing quadratic variation, the relaxation is trying to flatten the surface (to make the surface fiat, linear) to the maximum extend consistent with the data. To generalize the median to higher-dimensional settings, a variety of different maximum depth estimators have been introduced. They extend, for example, the halfplane location depth of T u k e y [7]. Let X\,..., Xn be independent and identically distributed bivariate random variables, Xi £ M2, i = 1, 2,..., n. For given observations x = (xi,..., xn)^ x^ £ M2, we write x^ = (xn, Xi2)T, where Xij are RW^ scores (here j = 1,2). The halfplane location depth of an arbitrary point 6 £ M2 relative to Xj is denned as di (6, x) = min#- # {i : Xj £ H} , i = 1,..., n, where H ranges over all closed halfplanes of which the boundary line passes through 9. The deepest location is denned as 6di = argmaxd/ (6*,x) and often called the Tukey median [6]. The depth region of depth s is denned as the set Ds of points 6 with di (6, x) > s. Equivalently, Ds is the intersection of all closed halfplanes that contains at least n — s + 1 observations, hence Ds is a bounded convex polygon, and Ds+i C Ds. The boundary of Ds is a convex polygon, which is called the contour of depth s (sth depth contour, sth hull). Therefore, each vertex of a depth contour is the intersection point of two lines, each passing through two observations. The univariate boxplot is based on the ranks since the box goes from the observation with the rank |_^J to that with the rank |~^p~|, and the central bar of the box is drawn at the median. A natural generalization of ranks to multivariate data is the notion of halfspace depth [7]. Using this concept in shape analysis, we propose a bivariate boxplot (bagplot) with the bag, that contains 50 % of the data points x^ a fence that separates inliers from outliers and a loop indicating the points outside the bag but inside the fence [6]. The loop plays the same role as the two whiskers in one dimension, so we could call this plot also ubag-and-bolster plot" to stress the analogy with "box-and-whisker plot". Like the univariate boxplot, bagplot visualizes several characteristics of the data—location, spread, correlation, skewness, and tails. The construction of bag B is as follows. Let #DS denote the number of points contained in Ds. One first determines the value s for which #DS < |_^J < #Ds-\ and then interpolates linearly between Ds and .Ds_i, relative to 9di, to obtain the set B (the bag B is a convex hull). The fence is obtained by inflating B, relative to 6di, by a factor 3 [8]. The loop contains all points between the bag and the fence and its outer boundary is the convex hull of the bag and the non-outliers. The points outside the fence are nagged as the outliers. 87 STANISLAV KATI NA Above mentioned Tukey depth and depth contoures can be directly applied to PCi and PC2 scores to detect possible outliers [4]. Furthermore, the bagplot fence can be used as stop-rule for relaxation of the configuration x, where this rule can be seen as multivariate winsorization. 4. Results and conclusions In the non-affine (non-uniform) subspaces using RWA, we find no outliers in the subspace of bending patterns at large scale but four outliers in the sub-space of local changes at small scale (Fig. 1; see also equation (1)). FIGURE 1. Results of RWA of matched-pair differences in non-affine sub-space of local changes with small scale (a = —1; RWi and rw2 scores with the bag and the fence of the bagplot (middle); TPS deformation grids and field of vectors superimposed with the grid of gray-scale rectangles with colors corresponding to the Procrustes distances of the mean shape Xpio to the Vec(Xp (cj , j, a)) (see equation (1)) indicated by the tips of the arrow heads in the middle figure). Almost the same outliers are found by form-space PC A (Fig. 2). These latter are associated with possible creases (extrema of directional derivative; [2]) 88 DETECTION OF NON-AFFINE SHAPE OUTLIERS FOR MATCHED-PAIR SHAPE DATA of the corresponding TPS. In those cases we can use the same spline formalism to relax the outlying form to an inlier by optimal relaxation along the curve decolletage that weighs bending energy against Procrustes distance (Fig. 3; see also equation (2)). FIGURE 2. Results of form-space PC A of matched-pair differences; PCi and pc2 scores with the bag and the fence of the bagplot (middle); TPS deformation grids and field of vectors superimposed with the grid of gray-scale rectangles with colors corresponding to the Procrustes distances of the mean shape Xpio to the Vec(Xp (ci, 1)) in PCi. FIGURE 3. Relaxation of Procrustes shape coordinates; TPS deformation grids and field of vectors superimposed with the grid of gray-scale rectangles with colors corresponding to the Procrustes distances of the shape Xp 10 29 to the shape Xp 15 29 (left) and to the final relaxed shape Xp 15 29 (right); curve decolletage of the shape Xp 15 29 (x—shape Xp,io,29, big •—shape Xpi5 29, small red •—relaxed shapes Xp 15,29; middle); the difference of Xp 15,29 and Xp 10,29 m the subspace of the first two RWs is visualized in the scatter-plot of Fig. 1 as the most right •. 89 STANISLAV KATI NA For biostatisticians, these maneuvers suggest a possibly novel and interesting fusion of the Procrustes-spline toolkit as the part of GMM with another theme, outlier detection, of the standard modern data-analytic toolkit. They also have practical implications for craniofacial management of CLP follow-up as well as suggestive implications for outlier detection in applied craniometries and anthropometrics more generally. Acknowledgements. The author thanks A. Bowman, F. L. Bookstein and P. O' H i g g i n s for comments and J. Velemínská for data acquisition, pre-processing, and valuable comments. REFERENCES [1] BOOKSTEIN, F. L.: Morphometric Tools for Landmark Data—Geometry and Biology. Cambridge University Press, Cambridge, 1991. [2] BOOKSTEIN, F. L.: Creases as local features of deformation grids, Med. Image Anal. 4 (2000), 93-110. [3] DRYDEN, I. L.—MARDIA, K. V.: Statistical Shape Analysis. John Willey, New York, 1999. [4] KATINA, S.: Geometrical Aspects of Statistical Size and Shape Analysis: Multivariate Statistical Models and Statistical Inference. Habilitation Thesis, Comenius University, Bratislava, 2008. [5] R DEVELOPMENT CORE TEAM: R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, 2011. [6] ROUSSEEUW, P. J.—RUTS, I.—TUKEY, J. W.: The bagplot: a bivariate boxplot, Amer. Statist. 53 (1999), 382-387. [7] TUKEY, J. W.: Mathematics and the picturing of data, Proc. Int. Congress Math. 2 (1975), 523-531. [8] VELEMÍNSKÁ, J.—KATINA, S.—ŠMAHEL, Z.—SEDLÁČKOVÁ, M.: Analysis of facial skeleton shape in patients with complete unilateral cleft lip and palate: Geometric morphometries, Acta Chirurgiae Plasticae 48 (2006), 26—32. [9] VELEMÍNSKÁ, J.—ŠMAHEL, Z.—KATINA, S.: Development prediction of sagittal intermaxillary relations in patients with complete unilateral cleft lip and palate during puberty, Acta Chirurgiae Plasticae 49 (2007), 41-46. Received October 31, 2011 School of Mathematics and Statistics The University of Glasgow University Gardens G12 8QW Glasgow SCOTLAND, UK Department of Applied Mathematics and Statistics Comenius University Mlynská dolina SK-842-48 Bratislava SLOVAKIA E-mail: Stanislav.katina@glasgow.ac.uk 90