Functional Additive Models on Manifolds of Planar Shapes and Forms

Abstract

The “shape” of a planar curve and/or landmark configuration is considered its equivalence class under translation, rotation, and scaling, its “form” its equivalence class under translation and rotation while scale is preserved. We extend generalized additive regression to models for such shapes/forms as responses respecting the resulting quotient geometry by employing the squared geodesic distance as loss function and a geodesic response function to map the additive predictor to the shape/form space. For fitting the model, we propose a Riemannian L₂-Boosting algorithm well suited for a potentially large number of possibly parameter-intensive model terms, which also yields automated model selection. We provide novel intuitively interpretable visualizations for (even nonlinear) covariate effects in the shape/form space via suitable tensor-product factorization. The usefulness of the proposed framework is illustrated in an analysis of (a) astragalus shapes of wild and domesticated sheep and (b) cell forms generated in a biophysical model, as well as (c) in a realistic simulation study with response shapes and forms motivated from a dataset on bottle outlines. Supplementary materials for this article are available online.

Keywords:

• Boosting
• Functional regression
• Shape analysis
• Tensor-product model
• Visualization

1 Introduction

In many imaging data problems, the coordinate system of recorded objects is arbitrary or explicitly not of interest. Statistical shape analysis (Dryden and Mardia 2016) addresses this point by identifying the ultimate object of analysis as the shape of an observation, reflecting its geometric properties invariant under translation, rotation and rescaling, or as its form (or size-and-shape) invariant under translation and rotation. This article establishes a flexible additive regression framework for modeling the shape or form of planar (potentially irregularly sampled) curves and/or landmark configurations in dependence on scalar covariates. A rich shape analysis literature has been developed for 2D or 3D landmark configurations—presenting for instance selected points of a bone or face—which are considered elements of Kendall’s shape space (see, e.g., Dryden and Mardia 2016). In many 2D scenarios, however, observed points describe a curve reflecting the outline of an object rather than dedicated landmarks (Adams, Rohlf, and Slice 2013). Considering outlines as images of (parameterized) curves shows a direct link to functional data analysis (FDA, Ramsay and Silverman 2005) and, in this context, we speak of functional shape/form data analysis. As in FDA, functional shape/form data can be observed on a common and often dense grid (regular/dense design) or on curve-specific often sparse grids (irregular/sparse design). While in the regular case, analysis often simplifies by treating curve evaluations as multivariate data, more general irregular designs gave rise to further developments in sparse FDA (e.g., Yao, Müller, and Wang 2005; Greven and Scheipl 2017), explicitly considering irregular measurements instead of pre-smoothing curves. To the best of our knowledge, we are the first to consider irregular/sparse designs in the context of functional shape/form analysis.

Shapes and forms are examples of manifold data. Petersen and Müller (2019) propose “Fréchet regression” for random elements in general metric spaces, which requires estimation of a (potentially negatively) weighted Fréchet mean for each covariate combination. Their implicit rather than explicit model formulation renders model interpretation difficult. More explicit model formulations have been developed for the special case of a Riemannian geometry. Besides tangent space models (Kent et al. 2001), extrinsic models (Lin et al. 2017) and models based on unwrapping (Jupp and Kent 1987; Mallasto and Feragen 2018), a variety of manifold regression models have been designed based on the intrinsic Riemannian geometry. Starting from geodesic regression (Fletcher 2013), which extends linear regression to curved spaces, these include MANOVA (Huckemann, Hotz, and Munk 2010), polynomial regression (Hinkle, Fletcher, and Joshi 2014), smoothing splines (Kume, Dryden, and Le 2007), regression along geodesic paths with nonconstant speed (Hong et al. 2014), or kernel regression (Davis et al. 2010) and Kriging (Pigoli, Menafoglio, and Secchi 2016). However, mostly only one metric covariate or categorical covariates are considered, possibly in hierarchical model extensions for longitudinal data (Muralidharan and Fletcher 2012; Schiratti et al. 2017). By contrast, Zhu et al. (2009), Shi et al. (2009), and Kim et al. (2014) generalize geodesic regression to regression with multiple covariates focusing on Symmetric Positive-Definite (SPD) matrix responses. Cornea et al. (2017) develop a general Generalized Linear Model (GLM) analogue regression framework for responses in a symmetric manifold and apply it to shape analysis. Recently, Lin, Müller, and Park (2020) proposed a Lie group additive regression model for Riemannian manifolds focusing on SPD matrices rather than shapes.

In FDA, there is a much wider range of developed regression methods (see overviews in Morris 2015; Greven and Scheipl 2017). Among the most flexible models are Functional Additive Models (FAMs) for (univariate) functional responses (in contrast to FAMs with functional covariates (Ferraty et al. 2011)) with either semi- or nonparametric approaches to model (a) response functions and (b) smooth covariate effects. For (a), nonparametric approaches formulate estimation problems in infinite-dimensional model spaces to motivate finite-dimensional representations or effectively evaluate curves on grids (e.g., Jeon and Park 2020). Semi-parametric approaches directly employ finite expansions in spline bases (Brockhaus, Scheipl, and Greven 2015), Functional Principal Component (FPC) bases (Morris and Carroll 2006) or both (Scheipl, Staicu, and Greven 2015), as well as wavelets (Meyer et al. 2015), sometimes directly expanding functions to model on coefficients and sometimes expanding only predictions while keeping the raw measurements. Nonparametric approaches are formulated in infinite-dimensional model spaces and effectively evaluate curves on grids or apply pre-smoothing techniques (e.g., Jeon and Park 2020). For (b), again semiparametric penalized spline basis approaches are employed (Scheipl, Staicu, and Greven 2015; Brockhaus, Scheipl, and Greven 2015), or local linear/polynomial (Müller and Yao 2008; Jeon et al. 2022) or other nonparametric kernel-based approaches (Jeon and Park 2020; Jeon, Park, and Van Keilegom 2021). Semi- and nonparametric approaches come with different theoretical and practical advantages, but similarities such as regarding asymptotic behavior are also known from scalar nonparametric regression (Li and Ruppert 2008). Advantages of the semi-parametric approach summarized in Greven and Scheipl (2017) include its appropriateness for sparse irregular functional data and its modular extensibility to functional mixed models (Scheipl, Staicu, and Greven 2015; Meyer et al. 2015) and nonstandard response distributions (Brockhaus, Scheipl, and Greven 2015; Stöcker et al. 2021). For bivariate or multivariate functional responses, which are closest to functional shapes/forms but without invariances, Rosen and Thompson (2009), Zhu, Li, and Kong (2012), Olsen, Markussen, and Raket (2018) consider linear fixed effects of scalar covariates, the latter also allowing for warping. Zhu et al. (2017), Backenroth et al. (2018) consider one or more random effects for one grouping variable, linear fixed effects and common dense grids for all functions. Volkmann et al. (2021) combine the FAM model class of Greven and Scheipl (2017) with multivariate FPC analysis (Happ and Greven 2018) to model multivariate (sparse) functional responses.

This article establishes an interpretable FAM framework for modeling the shape or form of planar (potentially irregularly sampled) curves and/or landmark configurations in dependence on scalar covariates, extending L₂-Boosting (Bühlmann and Yu 2003; Brockhaus, Scheipl, and Greven 2015) to Riemannian manifolds for model estimation. The three major contributions of our regression framework are: (i) We introduce additive regression with shapes/forms of planar curves and/or landmarks as response, extending FAMs to nonlinear response spaces or, vice versa, extending GLM-type regression on manifolds for landmark shapes both to functional shape manifolds and to include (nonlinear) additive model effects. (ii) We propose a novel Riemannian L₂-Boosting algorithm for estimating regression models for this type of manifold response, and (iii) a visualization technique based on tensor-product factorization yielding intuitive interpretations even of multi-dimensional smooth covariate effects for practitioners. Although related tensor-product model transformations based on higher-order SVD have been used, e.g., in control engineering (Baranyi, Yam, and Várlaki 2013), we are not aware of any comparable application for visualization in FAMs or other statistical models for object data. Despite our focus on shapes and forms, transfer of the model, Riemannian L₂-Boosting, and factorized visualization to other Riemannian manifold responses is intended in the generality of the formulation and the design of the provided R package manifoldboost (developer version on github.com/Almond-S/manifoldboost). The versatile applicability of the approach is illustrated in three different scenarios: an analysis of the shape of sheep astragali (ankle bones) represented by both regularly sampled curves and landmarks in dependence on categorical “demographic” variables; an analysis of the effects of different metric biophysical model parameters (including smooth interactions) on the form of (irregularly sampled) cell outlines generated from a cellular Potts model; and a simulation study with irregularly sampled functional shape and form responses generated from a dataset of different bottle outlines and including metric and categorical covariates.

In Section 2, we introduce the manifold geometry of irregular curves modulo translation, rotation and potentially rescaling, which underlies the intrinsic additive regression model formulated in Section 3. The Riemannian L²-Boosting algorithm is introduced in Section 4. Section 5 analyzes different data problems, modeling sheep bone shape responses (Section 5.1) and cell outlines (Section 5.2). Section 5.3 summarizes the results of simulation studies with functional shape and form responses. We conclude with a discussion in Section 6.

2 Geometry of Functional Shapes and Forms

Riemannian manifolds of planar shapes (and forms) are discussed in various textbooks at different levels of generality, in finite (Kendall et al. 1999; Dryden and Mardia 2016) or potentially infinite dimensions (Srivastava and Klassen 2016; Klingenberg 1995). Starting from the Hilbert space $\mathcal{Y}$ of curve representatives y of a single shape or form observation, we successively characterize its quotient space geometry under translation, rotation and rescaling including the respective tangent spaces. Building on that, we introduce Riemannian exponential and logarithmic maps and parallel transports needed for model formulation and fitting, and the sample space of (irregularly observed) functional shapes/forms.

To make use of complex arithmetic, we identify the two-dimensional plane with the complex numbers, $R^2\cong \mathrm{ }C$, and consider a planar curve to be a function $y:\mathrm{ }R\supset \mathcal{T}\to \mathrm{ }C$, element of a separable complex Hilbert space $\mathcal{Y}$ with a complex inner product $〈·,·〉$ and corresponding norm $\Vert ·\Vert$. This allows simple scalar expressions for the group actions of translation $\text{Trl}=\{y\stackrel{{\text{Trl}}_{\gamma }}{\to }y+\gamma \mathrm{}\mathbb{1}:\gamma \in C\}$ with $\mathbb{1}\in \mathcal{Y}$ canonically given by $\mathbb{1}:t\mapsto \frac{1}{\Vert t\mapsto 1\Vert }$ the real constant function of unit norm; rescaling $\text{Scl}=\{y\stackrel{{\text{Scl}}_{\lambda }}{\to }\lambda ·(y-0_y)+0_y:\lambda \in \mathrm{ }{R}^{+}\}$ around the centroid $0_y=〈\mathbb{1} ,y〉\mathbb{1}$ (which we consider more natural than using $0$, the zero element of $\mathcal{Y}$, mostly chosen in the literature); and rotation $\text{Rot}=\{y\stackrel{{\text{Rot}}_u}{\to }u·(y-0_y)+0_y:u\in S^1\}$ around $0_y$ with $S^1=\{u\in C:|u|=1\}=\{\hspace{0.17em}\exp \hspace{0.17em}(\omega \sqrt{-1}):\omega \in R\}$ reflecting counterclockwise rotations by ω radian measure. Concatenation yields combined group actions G as direct products, such as the rigid motions $G=\text{Trl}\times \text{Rot}=\{{\text{Trl}}_{\gamma }°{\text{Rot}}_u:\gamma \in C,u\in \mathrm{ }{S}^1\}\cong C\times S^1$ (see Section S.1.1, supplementary materials for more details). The two real-valued component functions of y are identified with the real part $\text{Re}\mathrm{}(y):\mathcal{T}\to \mathrm{ }R$ and imaginary part $\text{Im}\mathrm{}(y):\mathcal{T}\to \mathrm{ }R$ of $y=\text{Re}\mathrm{}(y)+\text{Im}\mathrm{}(y)\sqrt{-1}$. While the complex setup is used for convenience, the real part of $〈·,·〉$ constitutes an inner product $\text{Re}\mathrm{}(〈y_1,y_2〉)=〈\text{Re}\mathrm{}(y_1),\text{Re}\mathrm{}(y_2)〉+〈\text{Im}\mathrm{}(y_1),\text{Im}\mathrm{}(y_2)〉$ for $y_1,y_2\in \mathcal{Y}$ on the underlying real vector space of planar curves. Typically $\text{Re}\mathrm{}(y),\text{\hspace{0.17em}Im}\mathrm{}(y)$ are assumed square-intregrable with respect to a measure ν and we consider the canonical inner product $〈y_1,y_2〉=\int y_1^{†}{y}_{2}d\nu$ where $y^{†}$ denotes the conjugate transpose of y, that is, $y^{†}(t)=\text{Re}\mathrm{}(y)(t)-\text{Im}\mathrm{}(y)(t)\sqrt{-1}$ is simply the complex conjugate, but for vectors $\mathbf{y}\in C^k$, the vector ${\mathbf{y}}^{†}$ is also transposed. For curves, we typically assume ν to be the Lebesgue measure on $\mathcal{T}=[0,1]$; for landmarks, a standard choice is the counting measure on $\mathcal{T}=\{1,\dots ,k\}$.

The ultimate response object is given by the orbit ${[y]}_G=\{g(y):g\in G\}$ (or short $[y]$) of $y\in \mathcal{Y}$, the equivalence class under the respective combined group actions G: with $G=\text{Trl}\times \text{Rot}\times \text{Scl},[y]={[y]}_{\text{Trl}\times \text{Rot}\times \text{Scl}}=\{\lambda u y+\gamma \mathbb{1}:\lambda \in \mathrm{ }{R}^{+},u\in S^1,\gamma \in \mathrm{ }C\}$ is referred to as the shape of y and, for $G=\text{Trl}\times \text{Rot}$, $[y]={[y]}_{\text{Trl}\times \text{Rot}}=\{uy+\gamma \mathrm{}\mathbb{1}:u\in S^1,\gamma \in C\}$ as its form or size-and-shape. ${\mathcal{Y}}_{/G}=\{{[y]}_G:y\in \mathcal{Y}\}$ denotes the quotient space of $\mathcal{Y}$ with respect to G. The description of the Riemannian geometry of ${\mathcal{Y}}_{/G}$ involves, in particular, a description of the tangent spaces $T_{[y]}{\mathcal{Y}}_{/G}$ at points $[y]\in {\mathcal{Y}}_{/G}$, which can be considered local vector space approximations to ${\mathcal{Y}}_{/G}$ in a neighborhood of $[y]$. For a point q in a manifold $\mathcal{M}$ the tangent vectors $\beta \in T_q\mathcal{M}$ can, i.a., be thought of as gradients $\dot{c}(0)$ of paths $c:R\supset (-\delta ,\delta )\to \mathcal{M}$ at 0 where they pass through $c(0)=q$. Besides their geometric meaning, they will also play an important role in the regression model, as additive model effects are formulated on tangent space level. Choosing suitable representatives ${\tilde{y}}^G\in {[y]}_G\subset \mathcal{Y}$ (or short $\tilde{y}$) of orbits ${[y]}_G$, we use an identification of tangent spaces with suitable linear subspaces $T_{{[y]}_G}{\mathcal{Y}}_{/G}\subset \mathcal{Y}$.

Form geometry: Starting with translation as the simplest invariance, an orbit ${[y]}_{\text{Trl}}$ can be one-to-one identified with its centered representative ${\tilde{y}}^{\text{Trl}}=y-〈y,\mathrm{}\mathbb{1}〉 \mathrm{}\mathbb{1}$ yielding an identification ${\mathcal{Y}}_{/\text{Trl}}\cong \{y\in \mathcal{Y}:〈y,\mathrm{}\mathbb{1} 〉=0\}$ with a linear subspace of $\mathcal{Y}$. Hence, also $T_{[y]}{\mathcal{Y}}_{/\text{Trl}}=\{y\in \mathcal{Y}:〈y,\mathrm{}\mathbb{1} 〉=0\}$. For rotation, by contrast, we can only find local identifications with Hilbert subspaces (i.e., charts) around reference points ${[p]}_{\text{Trl}\times \text{Rot}}$ we refer to as “poles”. Moreover, we restrict to $y,p\in {\mathcal{Y}}^{*}=\mathcal{Y}\setminus {[0 ]}_{\text{Trl}}$ eliminating constant functions as degenerate special cases in the translation orbit of zero. For each ${[y]}_{\text{Trl}\times \text{Rot}}$ in an open neighborhood around ${[p]}_{\text{Trl}\times \text{Rot}}$ which can be chosen with $〈{\tilde{y}}^{\text{Trl}},{\tilde{p}}^{\text{Trl}}〉\ne 0$, y can be uniquely rotation aligned to p, yielding a one-to-one identification of the form ${[y]}_{\text{Trl}\times \text{Rot}}$ with the aligned representative given by ${\tilde{y}}^{\text{Trl}\times \text{Rot}}=\frac{〈{\tilde{y}}^{\text{Trl}},{\tilde{p}}^{\text{Trl}}〉}{|〈{\tilde{y}}^{\text{Trl}},{\tilde{p}}^{\text{Trl}}〉|}{\tilde{y}}^{\text{Trl}}={\text{argmin}}_{y\prime \in {[y]}_{\text{Trl}\times \text{Rot}}}\Vert y\prime -p\Vert$ (compare Figure 1). While ${\tilde{y}}^{\text{Trl}\times \text{Rot}}$ depends on p, we omit this in the notation for simplicity. All ${\tilde{y}}^{\text{Trl}}$ rotation aligned to ${\tilde{p}}^{\text{Trl}}$ lie on the hyper-plane determined by $\text{Im}\mathrm{}(〈{\tilde{y}}^{\text{Trl}},{\tilde{p}}^{\text{Trl}}〉)=0$ (Figure 1), which yields $T_{[p]}{\mathcal{Y}}_{/\text{Trl}+\text{Rot}}^{*}=\{y\in \mathcal{Y}:〈y,\mathrm{}\mathbb{1}〉=0,\text{\hspace{0.17em}Im}\mathrm{}(〈y,p〉)=0\}$ with normal vectors ${\zeta }^{(1)}=\mathrm{}\mathbb{1},{\zeta }^{(2)}=\sqrt{-1} \mathbb{1},{\zeta }^{(3)}=\sqrt{-1} p$. Note that, despite the use of complex arithmetic, $T_{[p]}{\mathcal{Y}}_{/\text{Trl}\times \text{Rot}}^{*}$ is a real vector space not closed under complex scalar multiplication. The geodesic distance of ${[y]}_{\text{Trl}\times \text{Rot}}$ to the pole ${[p]}_{\text{Trl}\times \text{Rot}}$ is given by $d({[y]}_{\text{Trl}\times \text{Rot}},{[p]}_{\text{Trl}\times \text{Rot}})=\Vert {\tilde{y}}^{\text{Trl}\times \text{Rot}}-{\tilde{p}}^{\text{Trl}}\Vert =\underset{y\prime \in {[y]}_{\text{Trl}\times \text{Rot}},p\prime \in {[p]}_{\text{Trl}\times \text{Rot}}}{\text{argmin}}\Vert y\prime -p\prime \Vert$. It reflects the length of the shortest path (i.e., the geodesic) between the forms and the minimum distance between the orbits as sets.

Fig. 1 Left: Quotient space geometry: assuming p and y centered, translation invariance is not further considered in the plot; given pole representative p, we express $y=\frac{\text{Re}\mathrm{}(〈p,y〉)}{\Vert p{\Vert }^2}p+\frac{\text{Im}\mathrm{}(〈p,y〉)}{\Vert p{\Vert }^2}ip+(y-\frac{〈p,y〉}{\Vert p{\Vert }^2}p)\in \mathcal{Y}$ in its coordinates in p and ip direction, subsuming all orthogonal directions in the third dimension. In this coordinate system, the rotation orbit ${[y]}_{\text{Rot}}$ corresponds to the dotted horizontal circle, and is identified with the aligned $\tilde{y}:={\tilde{y}}^{\text{Rot}}$ in the half-plane of p; ${[y]}_{\text{Rot}\times \text{Scl}}$ is identified with the unit vector ${\tilde{y}}^{\text{Rot}\times \text{Scl}}=\frac{\tilde{y}}{\Vert \tilde{y}\Vert }$ projecting $\tilde{y}$ onto the hemisphere depicted by the vertical semicircle. Form and shape distances between $[p]$ and $[y]$ correspond to the length of the geodesics $c(\tau )$ (thick lines) on the plane and sphere, respectively. Right: Geodesic line $c(\tau )$ between $p=c(0)$ and $p\prime =c(1)$, Log-map projecting y to $\epsilon \in T_p\mathcal{M}$, parallel transport ${\text{Transp}}_{pp\prime }$ forwarding ε to $\epsilon \prime \in T_{p\prime }\mathcal{M}$, and Exp-map projecting $\epsilon \prime$ onto $\mathcal{M}$ visualized for a sphere. Tangent spaces, identified with subspaces of the ambient space, are depicted as gray planes above the respective poles. The parallel transport preserves all angles between tangent vectors and identifies $\dot{c}(0)\cong \dot{c}(1)$.

Shape geometry: To account for scale invariance in shapes ${[y]}_{\text{Trl}\times \text{Rot}\times \text{Scl}}$, they are identified with normalized representatives ${\tilde{y}}^{\text{Trl}\times \text{Rot}\times \text{Scl}}=\frac{{\tilde{y}}^{\text{Trl}\times \text{Rot}}}{\Vert {\tilde{y}}^{\text{Trl}\times \text{Rot}}\Vert }$. Motivated by the normalization, we borrow the well-known geometry of the sphere $S=\{y\in \mathcal{Y}:\Vert y\Vert =1\}$, where $T_{p}S=\{y\in \mathcal{Y}:\text{Re}\mathrm{}(〈y,p〉)=0\}$ is the tangent space at a point $p\in S$ and geodesics are great circles. Together with translation and rotation invariance, the shape tangent space is then given by $T_{[p]}{\mathcal{Y}}_{/\text{Trl}\times \text{Rot}\times \text{Scl}}^{*}=T_{[p]}{\mathcal{Y}}_{/\text{Trl}\times \text{Rot}}^{*}\cap T_{p}S=\{y\in \mathcal{Y}:〈y,\mathbb{1}〉=0,〈y,p〉=0\}$ with normal vector ${\zeta }^{(4)}=p$ in addition to ${\zeta }^{(1)},{\zeta }^{(2)},{\zeta }^{(3)}$ above. The geodesic distance $d({[p]}_{\text{Trl}\times \text{Rot}\times \text{Scl}},{[y]}_{\text{Trl}\times \text{Rot}\times \text{Scl}})=\arccos |〈{\tilde{y}}^{\text{Trl}\times \text{Rot}\times \text{Scl}},{\tilde{p}}^{\text{Trl}\times \text{Rot}\times \text{Scl}}〉|$ corresponds to the arc-length between the representatives. This distance is often referred to as Procrustres distance in statistical shape analysis.

We may now define the maps needed for the regression model formulation. Let $\tilde{y}$ and $\tilde{p}$ be shape/form representatives of $[y]$ and $[p]$ rotation aligned to the shape/form pole representative p. Generalizing straight lines to a Riemannian manifold $\mathcal{M}$, geodesics $c:(-\delta ,\delta )\to \mathcal{M}$ can be characterized by their “intercept” $c(0)\in \mathcal{M}$ and “slope” $\dot{c}(0)\in T_{c(0)}\mathcal{M}$. The exponential map ${\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_q:T_q\mathcal{M}\to \mathcal{M}$ at a point $q\in \mathcal{M}$ is defined to map $\beta \mapsto c(1)$ for c the geodesic with $q=c(0)$ and $\beta =\dot{c}(0)$. It maps $\beta \in T_q\mathcal{M}$ to a point ${\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_q(\beta )\in \mathcal{M}$ located $d(q,{\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_q(\beta ))=\Vert \beta \Vert$ apart of the pole q in the direction of β. On the form space ${\mathcal{Y}}_{/\text{Trl}\times \text{Rot}}$, the exponential map is simply given by ${\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{{[p]}_{\text{Trl}\times \text{Rot}}}(\beta )=\mathrm{ }{\left[{\tilde{p}}^{\text{Trl}\times \text{Rot}}+\beta \right]}_{\text{Trl}\times \text{Rot}}$. On the shape space ${\mathcal{Y}}_{/\text{Trl}\times \text{Rot}\times \text{Scl}}$, identification with exponential maps on the sphere yields ${\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{{[p]}_G}(\beta )={\left[\cos \hspace{0.17em}(\Vert \beta \Vert ){\tilde{p}}^G+\hspace{0.17em}\sin \hspace{0.17em}(\Vert \beta \Vert )\frac{\beta }{\Vert \beta \Vert }\right]}_G$ with $G=\text{Trl}\times \text{Rot}\times \text{Scl}$. In an open neighborhood $\mathcal{U},q\in \mathcal{U}\subset \mathcal{M},{\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_q$ is invertible yielding the ${\hspace{0.17em}\mathrm{Log}\hspace{0.17em}}_q:\mathcal{U}\to T_q\mathcal{M}$ map from the manifold to the tangent space at q. For forms, it is given by ${\hspace{0.17em}\mathrm{Log}\hspace{0.17em}}_{{[p]}_{\text{Trl}\times \text{Rot}}}({[y]}_{\text{Trl}\times \text{Rot}})={\tilde{y}}^{\text{Trl}\times \text{Rot}}-{\tilde{p}}^{\text{Trl}\times \text{Rot}}$ and, for shapes, by ${\hspace{0.17em}\mathrm{Log}\hspace{0.17em}}_{{[p]}_G}({[y]}_G)=d({[p]}_G,{[y]}_G)\frac{{\tilde{y}}^G-〈{\tilde{p}}^G,{\tilde{y}}^G〉{\tilde{p}}^G}{\Vert {\tilde{y}}^G-〈{\tilde{p}}^G,{\tilde{y}}^G〉{\tilde{p}}^G\Vert }$ with $G=\text{Trl}\times \text{Rot}\times \text{Scl}$. Finally, ${\text{Transp}}_{q,q\prime }:T_q\mathcal{M}\to T_{q\prime }\mathcal{M}$ parallel transports tangent vectors $\epsilon \mapsto \epsilon \prime$ isometrically along a geodesic $c(\tau )$ connecting q and $q\prime \in \mathcal{M}$ such that the slopes ${\text{Transp}}_{q,q\prime }(\dot{c}(q))=\dot{c}(q\prime )$ are identified and all angles are preserved. For shapes, ${\text{Transp}}_{{[y]}_G,{[p]}_G}(\epsilon )=\epsilon -〈\epsilon ,{\tilde{p}}^G〉\frac{{\tilde{y}}^G+{\tilde{p}}^G}{1+〈{\tilde{y}}^G,{\tilde{p}}^G〉}$, with $G=\text{Trl}\times \text{Rot}\times \text{Scl}$, takes the form of the parallel transport on a sphere replacing the real inner product with its complex analogue. For forms, it changes only the $\text{Im}\mathrm{}(〈\epsilon ,\tilde{p}〉)$ coordinate orthogonal to the real $\tilde{y}$-$\tilde{p}$-plane as in the shape case, while the remainder of ε is left unchanged as in a linear space. This yields ${\text{Transp}}_{{[y]}_G,{[p]}_G}(\epsilon )=\epsilon -\text{Im}\mathrm{}(〈{\tilde{p}}^G/\Vert {\tilde{p}}^G\Vert ,\epsilon 〉)\frac{{\tilde{y}}^G/\Vert {\tilde{y}}^G\Vert +{\tilde{p}}^G/\Vert {\tilde{p}}^G\Vert }{1+〈{\tilde{y}}^G/\Vert {\tilde{y}}^G\Vert ,{\tilde{p}}^G/\Vert {\tilde{p}}^G\Vert 〉}\sqrt{-1}$, with $G=\text{Trl}\times \text{Rot}$, for form tangent vectors. While equivalent expressions for the parallel transport in the shape case can be found, for example, in Dryden and Mardia (2016), Huckemann, Hotz, and Munk (2010), a corresponding derivation for the form case is given in Section S.1.2, supplementary materials including a discussion of the quotient space geometry in differential geometric terms.

Based on this understanding of the response space, we may now proceed to consider a sample of curves $y_1,\dots ,y_n\in \mathcal{Y}$ representing orbits $[y_1],\dots ,[y_n]$ with respect to group actions G. In the functional case, with the domain $\mathcal{T}=[0,1]$, these curves are usually observed as evaluations ${\mathbf{y}}_i={(y_i(t_{i1}),\dots ,y_i(t_{i{k}_i}))}^{\top }$ on a finite grid $t_{i1}<\dots <t_{i{k}_i}\in \mathcal{T}$ which may differ between observations. In contrast to the regular case with common grids, this more general data structure is referred to as irregular functional shape/form data. To handle this setting, we replace the original inner product $〈·,·〉$ on $\mathcal{Y}$ by individual ${〈y_i,y_i^{\prime }〉}_i={\mathbf{y}}_i^{†}{\mathbf{W}}_i\mathbf{y}{\prime }_i$ providing inner products on the k_i -dimensional space ${\mathcal{Y}}_i=C^{k_i}$ of evaluations ${\mathbf{y}}_i,{\mathbf{y}}_i^{\prime }$ on the same grid. The symmetric positive-definite weight matrix ${\mathbf{W}}_i$ can be chosen to implement an approximation to integration w.r.t. the original measure ν with a numerical integration measure ν_i such as given by the trapezoidal rule. Alternatively, ${\mathbf{W}}_i=\frac{1}{k_i}{\mathbf{I}}_{k_i}$ with $k_i\times k_i$ identity matrix ${\mathbf{I}}_{k_i}$ presents a canonical choice that is analog to the landmark case for $k_i\equiv k$. Moreover, data-driven ${\mathbf{W}}_i$ could also be motivated from the covariance structure estimated for (potentially sparse) $y_1,\dots ,y_n$ along the lines of Yao, Müller, and Wang (2005), Stöcker et al. (2022). While this is beyond the scope of this article, potential procedures are sketched in Section S.7, supplementary materials. With the inner products given for $i=1,\dots ,n$, the sample space naturally arises as the Riemannian product ${\mathcal{Y}}_{1/G}^{*}\times \cdots \times {\mathcal{Y}}_{n/G}^{*}$ of the orbit spaces, with the individual geometries constructed as described above.

3 Additive Regression on Riemannian Manifolds

Consider a data scenario with n observations of a random response covariate tuple $(Y,\mathbf{X})$, where the realizations of Y are planar curves $y_i:\mathcal{T}\to \mathrm{ }C,i=1,\dots ,n$, belonging to a Hilbert space $\mathcal{Y}$ defined as above and potentially irregularly measured on individual grids $t_{i1}<\dots <t_{i{k}_i}\in \mathcal{T}$. The response object $[Y]$ is the equivalence class of Y with respect to translation, rotation and possibly scale and the sample $[y_1],\dots ,[y_n]$ is equipped with the respective Riemannian manifold geometry introduced in the previous section. For $i=1,\dots ,n$, realizations ${\mathbf{x}}_i\in \mathcal{X}$ of a covariate vector X in a covariate space $\mathcal{X}$ are observed. X can contain several categorical and/or metric covariates.

For regressing the mean of $[Y]$ on $\mathbf{X}=\mathbf{x}$, we model the shape/form $[\mu ]$ of $\mu \in \mathcal{Y}$ as(1) $$\left[\mu \right]={\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{\left[p\right]}\left(h(\mathbf{x})\right)={\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{\left[p\right]}\left(\sum _{j=1}^J{h}_j(\mathbf{x})\right),$$(1) with an additive predictor $h:\mathcal{X}\to T_{[p]}{\mathcal{Y}}_{/G}^{*}$ acting in the tangent space at an “intercept” $[p]\in {\mathcal{Y}}_{/G}^{*}$. Generalizing an additive model “$Y=\mu +ϵ=p+h(\mathbf{x})+ϵ$” in a linear space, we implicitly define $[\mu ]$ as the conditional mean of $[Y]$ given $\mathbf{X}=\mathbf{x}$ by assuming zero-mean “residuals” ϵ. In their definition, we follow Cornea et al. (2017) but extend to the functional shape/form and additive case. We assume local linearized residuals ${\epsilon }_{[\mu ]}={\hspace{0.17em}\mathrm{Log}\hspace{0.17em}}_{[\mu ]}([Y])$ in $T_{[\mu ]}{\mathcal{Y}}_{/G}^{*}$ to have mean $E({\epsilon }_{[\mu ]})=0$, which corresponds to $E({\epsilon }_{[\mu ]}(t))=0$ for (ν-almost) all $t\in \mathcal{T}$. Here, we assume $[Y]$ is sufficiently close to $[\mu ]$ with probability 1 such that ${\hspace{0.17em}\mathrm{Log}\hspace{0.17em}}_{[\mu ]}$ is well-defined, which is the case whenever $〈\tilde{Y},\tilde{\mu }〉\ne 0$ for centered shape/form representatives $\tilde{Y}$ and $\tilde{\mu }$, an unrestrictive and common assumption (compare also Cornea et al. 2017). However, residuals ${\epsilon }_{[\mu ]}$ for different $[\mu ]$ belong to separate tangent spaces. To obtain a formulation in a common linear space instead, local residuals are mapped to residuals $ϵ={\text{Transp}}_{[\mu ],[p]}({\epsilon }_{[\mu ]})$ by parallel transporting them from $[\mu ]$ to the common covariate independent pole $[p]$. After this isometric mapping into $T_{[p]}{\mathcal{Y}}_{/G}^{*}$, we can equivalently define the conditional mean $[\mu ]$ via $E(ϵ)=0$ for the transported residuals ϵ.

${\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{[p]}$ maps the additive predictor $h(\mathbf{x})={\sum }_{j=1}^J{h}_j(\mathbf{x})\in T_{[p]}{\mathcal{Y}}_{/G}^{*}$ to the response space. It is analogous to a response function in GLMs but depends on $[p]$. Although covariate effects $h_j(\mathbf{x})$ often only depend on an individual covariate in x for each j, they might also depend on covariate combinations in general to allow (smooth) interactions. While other response functions could be used, we restrict to the exponential map here, such that the model contains a geodesic model (Fletcher 2013)—the direct generalization of simple linear regression—as a special case for $h(\mathbf{x})=\beta x_1$ with a single covariate x₁ and tangent vector β. Typically, it is assumed that h is centered such that $E(h(\mathbf{X}))=0$, and the pole $[p]$ is the overall mean of $[Y]$ defined, like the conditional mean, via residuals of mean zero.

3.1 Tensor-Product Effect Functions h_j

Scheipl, Staicu, and Greven (2015) and other authors employ Tensor-Product (TP) bases for functional additive model terms. This naturally extends to tangent space effects, which we model as$$h_j(\mathbf{x})=\sum _{r,l}{\theta }_j^{(r,l)} b_j^{(l)}(\mathbf{x}) {\partial }_r$$with the TP basis given by the pair-wise products of $m_{}$ linearly independent tangent vectors ${\partial }_r\in T_{[p]}{\mathcal{Y}}_{/G}^{*},r=1,\dots ,m_{},$ and m_j basis functions $b_j^{(l)}:\mathcal{X}\to \mathrm{ }R,l=1,\dots ,m_j,$ for the jth covariate effect depending on one or more covariates. The real coefficients can be arranged as a matrix ${\left\{{\theta }_j^{(r,l)}\right\}}_{r,l}={\mathbf{\Theta }}_j\in \mathrm{ }{R}^{m_{}\times m_j}$. Also for infinite-dimensional $T_{[p]}{\mathcal{Y}}_{/G}^{*}$ and a general nonlinear dependence on x, a basis representation approach requires truncation to finite dimensions m and m_j in practice. Choosing the bases to capture the essential variability in the data, their size can be extended with increasing data size and computational resources.

While, in principle, the basis ${\left\{{\partial }_r\right\}}_r$ could also vary across effects $j=1,\dots ,J$, we assume a common basis for notational simplicity, which presents the typical choice. Due to the identification of $T_{[p]}{\mathcal{Y}}_{/G}^{*}$ with a subspace of the function space $\mathcal{Y}$, the ${\left\{{\partial }_r\right\}}_r$ may be specified using a function basis commonly used in additive models: Let $b_0^{(l)}:\mathcal{T}\to \mathrm{ }R,l=1,\dots ,m_0$ be a basis of real functions, say a B-spline basis (other typical bases used in the literature include wavelet (Meyer et al. 2015) or FPC bases (Müller and Yao 2008)). Then we construct the tangent space basis as ${\partial }_r={\sum }_{l=1}^{m_0}(z_p^{(l,r)}+z_p^{(m_0+l,r)}\sqrt{-1})b_0^{(l)}$, employing the same basis for the 1- and $\sqrt{-1}$-dimension before transforming it with a basis transformation matrix ${\mathbf{Z}}_p={\left\{z_p^{(l,r)}\right\}}_{l,r}\in \mathrm{ }{R}^{2m_0\times m_{}}$ implementing the linear tangent space constraints $\text{Re}\mathrm{}(〈{\partial }_l,{\zeta }^{(r)}〉)=0$ (or the empirical version) for all ${\partial }_l$ and normal vectors ${\zeta }^{(1)},{\zeta }^{(2)},{\zeta }^{(3)}$ for forms and additionally ${\zeta }^{(4)}$ for shapes defining $T_{[p]}{\mathcal{Y}}_{/G}^{*}$ as described in Section 2. Thus, the tangent space basis dimension is $m_{}=2m_0-3$ for forms or m = 2m₀ − 4 for shapes (or could, in principle, be larger if the original basis already meets the constraints). For details on the construction of ${\mathbf{Z}}_p$ see Section S.1.3, supplementary materials. For closed curves, we additionally choose ${\mathbf{Z}}_p$ to enforce periodicity, that is, ${\partial }_r(t)={\partial }_r(t+t_0)$ for some $t_0\in \mathrm{ }R$ (compare Hofner, Kneib, and Hothorn 2016).

Given the tangent space basis, we may now modularly specify the usual additive model basis functions $b_j^{(l)}:\mathcal{X}\to \mathrm{ }R,l=1,\dots ,m_j$, for the jth covariate effect to obtain the full functional additive model “tool box” offered by, for example, Brockhaus, Scheipl, and Greven (2015). Typically, $b_j^{(l)}(\mathbf{x})=b_j^{(l)}(z)$ depending on an individual covariate, say on z, in $\mathbf{x}={(\dots ,z,\dots )}^{\top }$. But for a single covariate also multiple different effects can be specified and a single interaction effect depends on multiple covariates. A linear effect—linear in the tangent space—of the form $h_j(\mathbf{x})=\beta z$ of a scalar (typically centered) covariate z and $\beta \in T_{[p]}{\mathcal{Y}}_{/G}^{*}$ is simply implemented by a single function $b_j^{(1)}(\mathbf{x})=z$. A smooth effect of the generic form $h_j(\mathbf{x})(t)=f(z,t)$ can be implemented by choosing, for example, a B-spline basis $b_j^{(1)}(z),\dots ,b_j^{(m_j)}(z)$ (asymptotic properties of penalized B-splines and connections to kernel estimators are discussed, for example, by Wood, Pya, and Säfken (2016), Li and Ruppert (2008)). For a categorical covariate κ in x, with effect $h_j(\mathbf{x}):\{1,\dots ,K\}\to T_{[p]}{\mathcal{Y}}_{/G}^{*},\kappa \mapsto {\beta }_{\kappa }$, the basis ${\mathbf{b}}_j(\mathbf{x})={(b_j^{(1)}(\kappa ),\dots ,b_j^{(m_j)}(\kappa ))}^{\top }$ maps $\kappa \mapsto {\mathbf{e}}_{\kappa }$ to a usual contrast vector ${\mathbf{e}}_{\kappa }$ with the basis being of dimension $m_j=K-1$ just as in standard linear models. Here, we typically use effect-encoding to obtain centered effects. Moreover, TP interactions of the model terms described above, as well as group-specific effects and smooth effects with additional constraints (Hofner, Kneib, and Hothorn 2016) can be specified in the model formula, relying on the mboost framework introduced by Hothorn et al. (2010), which also allows to define custom effect designs. For identification of an overall mean intercept $[p]$, sum-to-zero constraints yielding ${\sum }_{i=1}^n{h}_j({\mathbf{x}}_i)=0$ for observed covariates ${\mathbf{x}}_i$ can be specified, and similar constraints can be used to distinguish linear from nonlinear effects and interactions from their marginal effects (Kneib, Hothorn, and Tutz 2009). Different quadratic penalties can be specified for the coefficients ${\mathbf{\Theta }}_j$, allowing to regularize high-dimensional effect bases and to balance effects of different complexity in the model fit (see, Section 4).

3.2 Tensor-Product Factorization

The multidimensional structure of the response objects makes it challenging to graphically illustrate and interpret additive model terms, in particular when it comes to nonlinear (interaction) effects, or when effect sizes are visually small. To solve this problem, we suggest to rewrite estimated TP effects ${\widehat{h}}_j$ with estimated coefficient matrix ${\widehat{\mathbf{\Theta }}}_j$ as$${\widehat{h}}_j(\mathbf{x})=\sum _{r=1}^{m_j^{\prime }}{\xi }_j^{(r)}{\widehat{h}}_j^{(r)}(\mathbf{x})$$factorized into $m_j^{\prime }=\min (m_j,m_0)$ components consisting of covariate effects ${\widehat{h}}_j^{(r)}:\mathcal{X}\to \mathrm{ }R,r=1,\dots ,m_j^{\prime },$ in corresponding orthonormal directions ${\xi }_j^{(r)}\in T_{[p]}{\mathcal{Y}}_{/G}^{*}$ with $〈{\xi }_j^{(r)},{\xi }_j^{(l)}〉=\mathrm{ }1(r=l)$, that is, 1 if r = l and 0 otherwise. Assuming $E(b_j^{(l)}{(\mathbf{X})}^2)<\infty ,l=1,\dots ,m_j$, for the underlying effect basis, the ${\widehat{h}}_j^{(r)}$ are specified to achieve decreasing component variances $v_j^{(1)}\ge \dots \ge v_j^{(m_j^{\prime })}\ge 0$ given by $v_j^{(r)}=\mathrm{ }E({\widehat{h}}_j^{(r)}{(\mathbf{X})}^2)$. In practice, the expectation over the covariates X and the inner product $〈.,.〉$ are replaced by empirical analogs (compare Corollary 3, supplementary materials). Due to orthonormality of the ${\xi }_j^{(r)}$, the component variances add up to the total predictor variance $\sum _{r=1}^{m_j^{\prime }}{v}_j^{(r)}=v_j=E(〈{\widehat{h}}_j(\mathbf{X}),{\widehat{h}}_j(\mathbf{X})〉)$. Moreover, the TP factorization is optimally concentrated in the first components in the sense that for any $l\le m_j^{\prime }$ there is no sequence of ${\xi }_{*}^{(r)}\in \mathcal{Y}$ and ${\widehat{h}}_{*}^{(r)}:\mathcal{X}\to R$, such that $E(\Vert {\widehat{h}}_j(\mathbf{X})-{\sum }_{r=1}^l{\xi }_{*}^{(r)}{\widehat{h}}_{*}^{(r)}(\mathbf{X}){\Vert }^2)<\mathrm{ }E(\Vert h_j(\mathbf{X})-{\sum }_{r=1}^l{\xi }_j^{(r)}{\widehat{h}}_j^{(r)}(\mathbf{X}){\Vert }^2)$, that is, the series of the first l components yields the best rank l approximation of ${\widehat{h}}_j$. The factorization relies on SVD of (a transformed version of) the coefficient matrix ${\widehat{\mathbf{\Theta }}}_j$ and the fact that it is well-defined is a variant of the Eckart-Young-Mirsky theorem (proof in Section S.2, supplementary materials).

Particularly when large shares of the predictor variance are explained by the first component(s), the decomposition facilitates graphical illustration and interpretation: choosing a suitable constant τ ≠ 0, an effect direction ${\xi }_j^{(r)}$ can be visualized by plotting the pole representative p together with ${\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_p(\tau {\xi }_j^{(r)})$ on the level of curves, while accordingly rescaled $\frac{1}{\tau }{\widehat{h}}_j^{(r)}(\mathbf{x})$ is displayed separately in a standard scalar effect plot. Adjusting τ offers an important degree of freedom for visualizing ${\xi }_j^{(r)}$ on an intuitively accessible scale while faithfully depicting ${\xi }_j^{(r)}{\widehat{h}}_j^{(r)}(\mathbf{x})$. When based on the same τ, different covariate effects can be compared across the plots sharing the same scale. We suggest $\tau ={\max }_j\sqrt{v_j}$, the maximum total predictor standard deviation of an effect, as a good first choice.

Besides factorizing effects separately, it can also be helpful to apply TP factorization to the joint additive predictor, yielding$$h(\mathbf{x})=\sum _{r=1}^{m\prime }{\xi }^{(r)}{\widehat{h}}^{(r)}(\mathbf{x})=\sum _{r=1}^{m\prime }{\xi }^{(r)}\left({\widehat{h}}_1^{(r)}(\mathbf{x})+\dots +{\widehat{h}}_J^{(r)}(\mathbf{x})\right),$$with $m\prime =\min (\sum _j{m}_j,m_{})$ and again ${\xi }^{(r)}\in T_{[p]}{\mathcal{Y}}_{/G}^{*}$ orthonormal and the corresponding variance concentration in the first components, but now determined w.r.t. entire additive predictors ${\widehat{h}}^{(r)}={\sum }_{j=1}^J{\widehat{h}}_j^{(r)}$ spanned by all covariate basis functions in the predictor. In this representation, the first component yields a geodesic additive model approximation where the predictor moves along a geodesic line $c(\tau )={\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{[p]}({\xi }^{(1)}\tau )$ with the signed distance $\tau \in R$ from $[p]$, modeled by a scalar additive predictor ${\widehat{h}}^{(1)}(\mathbf{x})$ composed of covariate effects analogous to the original model predictor. In Section 5, we illustrate its potential in three different scenarios.

4 Component-Wise Riemannian L₂-Boosting

Component-wise gradient boosting (e.g., Hothorn et al. 2010) is a step-wise model fitting procedure accumulating predictors from smaller models, so called base-learners, to built an ensemble predictor aiming at minimizing a mean loss function. To this end, the base-learners are fit (via least squares) to the negative gradient of the loss function in each step and the best fitting base-learner is added to the current ensemble predictor. Due to its versatile applicability, inherent model selection, and slow over-fitting behavior, boosting has proven useful in various contexts (Mayr et al. 2014). Boosting with respect to the least squares loss function $\mathcal{l}(y,\mu )=\frac{1}{2}{(y-\mu )}^2,y,\mu \in R$, is typically referred to as L₂-Boosting and simplifies to repeated refitting of residuals $\epsilon =y-\mu =-{\nabla }_{\mu }\mathcal{l}(y,\mu )$ corresponding to the negative gradient of the loss function. For L₂-Boosting with a single learner, Bühlmann and Yu (2003) show how fast bias decay and slow variance increase over the boosting iterations suggest stopping the algorithm early before approaching the ordinary (penalized) least squares estimator. Lutz and Bühlmann (2006) prove consistency of component-wise L²-Boosting in a high-dimensional multivariate response linear regression setting and Stöcker et al. (2021) illustrate in extensive simulation studies how stopping the boosting algorithm early based on curve-wise cross-validation applies desired regularization when fitting (even highly autocorrelated) functional responses with parameter-intense additive model base-learners and, thus, leads to good estimates even in challenging scenarios.

When generalizing to least squares on Riemannian manifolds with the loss $\frac{1}{2}{d}^2([y],[\mu ])$ given by the squared geodesic distance, the negative gradient $-{\nabla }_{[\mu ]}\frac{1}{2}{d}^2([y],[\mu ])={\hspace{0.17em}\mathrm{Log}\hspace{0.17em}}_{[\mu ]}([y])={\epsilon }_{[\mu ]}$ (compare e.g., Pennec 2006) corresponds to the local residuals ${\epsilon }_{[\mu ]}$ defined in Section 3. This analogy to L₂-Boosting motivates the presented generalization where local residuals are further transported to residuals ϵ in a common linear space.

Consider the pole $[p]$ known and fixed for now. Assuming its existence, we aim to minimize the population mean loss$${\sigma }^2(h)=E\left(d^2\left(\left[Y\right],{\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{\left[p\right]}\left(h(\mathbf{X})\right)\right)\right)$$with the point-wise minimizer $h^{\star }(\mathbf{x})=\underset{{}_{h:\mathcal{X}\to T_{[p]}{\mathcal{Y}}_{/G}^{*}}}{\text{argmin}}E(d^2([Y],{\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{[p]}(h(\mathbf{X})))|\mathbf{X}=\mathbf{x})$ minimizing the conditional expected squared distance. Fixing a covariate constellation $\mathbf{x}\in \mathcal{X}$, the prediction $[\mu ] = {\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{[p]}(h^{\star }(\mathbf{x}))$ corresponds to the Fréchet mean (Karcher 1977) of $[Y]$ conditional on $\mathbf{X} = \mathbf{x}$. In a finite-dimensional context, Pennec (2006) show that $E({\epsilon }_{[\mu ]}) =0$ for a Fréchet mean $[\mu ]$ if residuals ${\epsilon }_{[\mu ]}$ are uniquely defined with probability one. This indicates the connection to our residual based model formulation in Section 3. We fit the model by reducing the empirical mean loss ${\widehat{\sigma }}^2(h)=\frac{1}{n}{\sum }_{i=1}^n{d}_i^2([y_i],{\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{[p]}(h({\mathbf{x}}_i))),$ where we replace the population mean by the sample mean and compute the geodesic distances d_i with respect to the inner products ${〈·,·〉}_i$ defined for the respective evaluations of y_i .

A base-learner corresponds to a covariate effect $h_j(\mathbf{x})=\sum _{r,l}{\theta }_j^{(r,l)} b_j^{(l)}(\mathbf{x}) {\partial }_r,{\mathbf{\Theta }}_j={\left\{{\theta }_j^{(r,l)}\right\}}_{r,l}$, which is repeatedly fit to the transported residuals ${ϵ}_1,\dots ,{ϵ}_n$ by penalized least-squares (PLS) minimizing ${\sum }_{i=1}^n\Vert {ϵ}_i-h_j({\mathbf{x}}_i){\Vert }_i^2+{\lambda }_j\text{tr}({\mathbf{\Theta }}_j{\mathbf{P}}_j{\mathbf{\Theta }}_j^{\top })+{\lambda }_{}\text{tr}({\mathbf{\Theta }}^{\top }{\mathbf{P}}_{}\mathbf{\Theta })$. Via the penalty parameters ${\lambda }_j,{\lambda }_{}\ge 0$ the effective degrees of freedom of the base-learners are controlled (Hofner et al. 2011) to achieve a balanced “fair” base-learner selection despite the typically large and varying number of coefficients involved in the TP effects. The symmetric penalty matrices ${\mathbf{P}}_j\in R^{m_j\times m_j}$ and ${\mathbf{P}}_{}\in R^{m_{}\times m_{}}$ (imposing, e.g., a second-order difference penalty for B-splines in either direction) can equivalently be arranged as a $m_j{m}_{}\times m_j{m}_{}$ penalty matrix ${\mathbf{R}}_j={\lambda }_j({\mathbf{P}}_j\otimes {\mathbf{I}}_{m_{}})+{\lambda }_{}({\mathbf{I}}_{m_j}\otimes {\mathbf{P}}_{})$ for the vectorized coefficients $\text{vec}({\mathbf{\Theta }}_j)={({\theta }_j^{(1,1)},\dots ,{\theta }_j^{(m_{},1)},\dots ,{\theta }^{(m_{},m_j)})}^{\top }$, where ⊗ denotes the Kronecker product. The standard PLS estimator is then given by $\text{vec}({\widehat{\mathbf{\Theta }}}_j)={({\mathbf{\Psi }}_j+{\mathbf{R}}_j)}^{-1}{\mathit{\psi }}_j$ with ${\mathbf{\Psi }}_j={\sum }_{i=1}^n\{\text{Re}\mathrm{}({〈b_j^{(l)}({\mathbf{x}}_i){\partial }_r,b_j^{(l\prime )}({\mathbf{x}}_i){\partial }_{r\prime }〉}_i){\}}_{\begin{array}{c}(r,l)=(1,1),\dots ,(m_{},1),\dots ,(m_{},m_j)\\ (r\prime ,l\prime )=(1,1),\dots ,(m_{},1),\dots ,(m_{},m_j)\end{array}}\in \mathrm{ }{R}^{m_{}\mathrm{ }{m}_j\times m_{}\mathrm{ }{m}_j}$ and ${\mathit{\psi }}_j={\sum }_{i=1}^n\{\text{Re}({〈b_j^{(l)}({\mathbf{x}}_i){\partial }_r,{ϵ}_i〉}_i){\}}_{(r,l)=(1,1),\dots ,(m_{},1),\dots ,(m_{},m_j)}\in R^{m_{}\mathrm{ }{m}_j}$. In a regular design, using the functional linear array model (Brockhaus, Scheipl, and Greven 2015) can save memory and computation time by avoiding construction of the complete matrices. The basis construction of ${\left\{{\partial }_r\right\}}_r$ via a transformation matrix ${\mathbf{Z}}_p$ (Section 3.1) is reflected in the penalty by setting ${\mathbf{P}}_{}={\mathbf{Z}}_p^{\top }({\mathbf{I}}_2\otimes {\mathbf{P}}_0){\mathbf{Z}}_p$ with ${\mathbf{P}}_0$ the penalty matrix for the un-transformed basis ${\left\{b_0^{(r)}\right\}}_r$.

In each iteration of the proposed Algorithm 1, the best-performing base-learner is added to the current ensemble additive predictor $h(\mathbf{x})$ after multiplying it with a step-length parameter $\eta \in (0,1]$. Due to the additive model structure this corresponds to a coefficient update of the selected covariate effect. Accordingly, after repeated selection, the effective degrees of freedom of a covariate effect, in general, exceed the degrees specified for the base-learner. They are successively adjusted to the data. To avoid over-fitting, the algorithm is typically stopped early before reaching a minimum of the empirical mean loss. The stopping iteration is determined, for example, by resampling strategies such as bootstrapping or cross-validation on the level of shapes/forms.

Algorithm 1:

Component-wise Riemannian L²-Boosting

# Initialization:

Geometry: specify geometry (shape/form) and pole representative p

Hyper-parameters: Step-length $\eta \in (0,1]$, number of boosting iterations

Base-learners: $h_j(\mathbf{x})$ with penalty matrix ${\mathbf{R}}_j$ and initial coefficient matrix ${\mathbf{\Theta }}_j=\mathbf{0}$

for j = 1 to J do # Prepare penalized least-squares (PLS)

# set up $m_{} m_j\times m_{} m_j$ matrix: ${\mathbf{\Psi }}_j\leftarrow {\sum }_{i=1}^n{\left\{\text{Re}\mathrm{}({〈b_j^{(l)}({\mathbf{x}}_i){\partial }_r,b_j^{(l\prime )}({\mathbf{x}}_i){\partial }_{r\prime }〉}_i)\right\}}_{\begin{array}{c}(r,l)=(1,1),\dots ,(m_{},1),\dots ,(m_{},m_j)\\ (r\prime ,l\prime )=(1,1),\dots ,(m_{},1),\dots ,(m_{},m_j)\end{array}}$

end

repeat # boosting steps

for $i=1,\dots ,n$ do # Compute current transported residuals

$[{\mu }_i]\leftarrow {\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{[p]}(h({\mathbf{x}}_i))$

${\epsilon }_{[{\mu }_i]}\leftarrow {\hspace{0.17em}\mathrm{Log}\hspace{0.17em}}_{[{\mu }_i]}([y_i])$

${ϵ}_i\leftarrow {\text{Transp}}_{[{\mu }_i],[p]}({\epsilon }_{[{\mu }_i]})$

end

for $j=1,\dots ,J$ do # PLS fit to residuals

# $m_{} m_j$ vector: ${\mathit{\psi }}_j\leftarrow {\sum }_{i=1}^n{\left\{\text{Re}\mathrm{}({〈b_j^{(l)}({\mathbf{x}}_i){\partial }_r,{ϵ}_i〉}_i)\right\}}_{(r,l)=(1,1),\dots ,(m_{},1),\dots ,(m_{},m_j)}$

${\widehat{\mathbf{\Theta }}}_j={\left\{{\widehat{\theta }}_j^{(r,l)}\right\}}_{r,l}\leftarrow$ Solve(

$({\mathbf{\Psi }}_j+{\mathbf{R}}_j)\text{vec}(\mathbf{\Theta })={\mathit{\psi }}_j$)

end

$\widehat{\mathit{ȷ}}\leftarrow \underset{{}_{j\in \{1,\dots ,J\}}}{\text{argmin}}{\sum }_{i=1}^n\Vert {ϵ}_i-\sum _{r,l}{\widehat{\theta }}_j^{(r,l)}{b}_j^{(l)}(\mathbf{x}){\partial }_r{\Vert }_i^2$;

# Select base-learner

${\mathbf{\Theta }}_{\widehat{\mathit{ȷ}}}\leftarrow {\mathbf{\Theta }}_{\widehat{\mathit{ȷ}}}+\eta {\widehat{\mathbf{\Theta }}}_{\widehat{\mathit{ȷ}}}$; # Update selected model coefficients

until Stopping criterion (e.g., minimal cross-validation error)

The pole $[p]$ is, in fact, usually not a priori available. Instead we typically assume $[p]=\underset{{}_{q\in {\mathcal{Y}}^{*}}}{\text{argmin}}E(d^2([Y],[q]))$ is the overall Fréchet mean, also often referred to as Riemannian center of mass for Riemannian manifolds or as Procrustes mean in shape analysis (Dryden and Mardia 2016). Here, we estimate it as $[p]={\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{[p_0]}(h_0)$ in a preceding Riemannian L²-Boosting routine. The constant effect $h_0\in T_{[p_0]}{\mathcal{Y}}_{/G}^{*}$ in the intercept-only special case of our model is estimated with Algorithm 1 based on a preliminary pole $[p_0]\in {\mathcal{Y}}_{/G}^{*}$. For shapes and forms, a good candidate for p₀ can be obtained as the standard functional mean of a reasonably well aligned sample $y_1,\dots ,y_n\in \mathcal{Y}$ of representatives.

The proposed Riemannian L₂-Boosting algorithm is available in the R (R Core Team 2018) package manifoldboost (github.com/Almond-S/manifoldboost). The implementation is based on the package FDboost (Brockhaus, Rügamer, and Greven 2020), which is in turn based on the model-based boosting package mboost (Hothorn et al. 2010).

5 Applications and Simulation

5.1 Shape Differences in Astragali of Wild and Domesticated Sheep

In a geometric morphometric study, Pöllath, Schafberg, and Peters (2019) investigate shapes of sheep astragali (ankle bones) to understand the influence of different living conditions on the micromorphology of the skeleton. Based on a total of n = 163 shapes recorded by Pöllath, Schafberg, and Peters (2019), we model the astragalus shape in dependence on different variables, including domestication status (wild/feral/domesticated), sex (female/male/NA), age (juvenile/subadult/adult/NA), and mobility (confined/pastured/free) of the animals as categorical covariates. The sample comprises sheep of four different populations: Asiatic wild sheep (Field Museum, Chicago; Lay 1967; Zeder 2006), feral Soay sheep (British Natural History Museum, London; Clutton-Brock et al. 1990), and domestic sheep of the Karakul and Marsch breed (Museum of Livestock Sciences, Halle (Saale); Schafberg and Wussow 2010). Table S1 in Section S.3, supplementary materials shows the distribution of available covariates within the populations. Each sheep astragalus shape, $i=1,\dots ,n$, is represented by a configuration composed of 11 selected landmarks in a vector ${\mathbf{y}}_i^{\text{lm}}\in C^{11}$ and two vectors of sliding semi-landmarks ${\mathbf{y}}_i^{\mathrm{c}1}\in C^{14}$ and ${\mathbf{y}}_i^{\mathrm{c}2}\in C^{18}$ evaluated along two outline curve segments, marked on a 2D image of the bone (dorsal view). Several example configurations are displayed in Figure S1, supplementary materials. In general, we could separately specify smooth function bases for the outline segments $y_i^{\mathrm{c}1}$ and $y_i^{\mathrm{c}2}$, respectively. Due to their systematic recording, we assume, however, that not only landmarks but also semi-landmarks are regularly observed on a fixed grid, and refrain from using smooth function bases for simplicity. Accordingly, shape configurations can directly be identified with their evaluation vectors ${\mathbf{y}}_i={({\mathbf{y}}_i^{\text{lm}\top },{\mathbf{y}}_i^{\mathrm{c}1\top },{\mathbf{y}}_i^{\mathrm{c}2\top })}^{\top }\in C^{43}=\mathcal{Y}$, and the geometry of the response space ${\mathcal{Y}}_{/\text{Trl}\times \text{Rot}\times \text{Scl}}^{*}$ widely corresponds to the classic Kendall’s shape space geometry, with the difference that, considering landmarks more descriptive than single semi-landmarks, we choose a weighted inner product $〈{\mathbf{y}}_i,{\mathbf{y}}_i^{\prime }〉={\mathbf{y}}_i^{†}\mathbf{W}{\mathbf{y}}_i^{\prime }$ with diagonal weight matrix W with diagonal ${({\mathbf{1}}_{11}^{\top },\frac{3}{14}{\mathbf{1}}_{14}^{\top },\frac{3}{18}{\mathbf{1}}_{18}^{\top })}^{\top }$ assigning the weight of three landmarks to each outline segment. We model the astragalus shapes $[{\mathbf{y}}_i]\in {\mathcal{Y}}_{/\text{Trl}\times \text{Rot}\times \text{Scl}}^{*}$ as$$\left[{\mathit{\mu }}_i\right]={\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{\left[\mathbf{p}\right]}\left({\mathit{\beta }}_{{\text{status}}_i}+{\mathit{\beta }}_{{\text{pop}}_i}+{\mathit{\beta }}_{{\text{age}}_i}+{\mathit{\beta }}_{{\text{sex}}_i}+{\mathit{\beta }}_{{\text{mobility}}_i}\right)$$with the pole $[\mathbf{p}]\in {\mathcal{Y}}_{/G}^{*}$ specified as overall mean and the conditional mean $[{\mathit{\mu }}_i]\in {\mathcal{Y}}_{/\text{Trl}\times \text{Rot}\times \text{Scl}}^{*}$ depending on the effect coded covariate effects $x_{ij}\mapsto {\mathit{\beta }}_{x_{ij}}\in T_{[\mathbf{p}]}{\mathcal{Y}}_{/\text{Trl}\times \text{Rot}\times \text{Scl}}^{*}$. For identifiability, the population and mobility effects are centered around the status effect, as we only have data on different populations/mobility levels for domesticated sheep. All base-learners are regularized to one degree of freedom by employing ridge penalties for the coefficients of the covariate bases ${\left\{b_j^{(l)}\right\}}_l$ while the coefficients of the response basis (the standard basis for $C^{43})$ are left un-penalized. With a step-length of $\eta =0.1$, 10-fold shape-wise cross-validation suggests early stopping after 89 boosting iterations. Due to the regular design, we can make use of the functional linear array model (Brockhaus, Scheipl, and Greven 2015) for saving computation time and memory, which lead to 8 sec of initial model fit followed by 47 sec of cross-validation. To interpret the categorical covariate effects, we rely on TP factorization (Figure 2). The first component of the status effect explains about 2/3 of the variance of the status effect and over 50% of the cumulative effect variance in the model. In that main direction, the effect of feral is not located between wild and domestic, as might be naively expected. By contrast, the second component of the effect seems to reflect the expected order and still explains a considerable amount of variance. Similar to Pöllath, Schafberg, and Peters (2019), we find little influence of age, sex, and mobility on the astragalus shape. Yet, all covariates were selected by the boosting algorithm.

Fig. 2 Left: Shares of different factorized covariate effects in the total predictor variance. Right: Factorized effect plots showing the two components of the status effect (rows): in the right column, the two first directions ${\mathit{\xi }}_1^{(1)},{\mathit{\xi }}_1^{(2)}\in T_{[\mathbf{p}]}{\mathcal{Y}}_{/\text{Trl}+\text{Rot}+\text{Scl}}^{*}$ are visualized via line-segments originating at the overall mean shape (empty circles) and ending in the shape resulting from moving 1 unit into the target direction (solid circles; large: landmarks; small: semi-landmarks along the outline); in the left column, the status effect in the respective direction is depicted. As illustrated in the middle plot, an effect of 1 would correspond to the full extend of the direction shown to the right.

Visually, differences in estimated mean shapes are rather small, which is, in our experience, quite usual for shape data. With differences in size, rotation and translation excluded by definition, only comparably small variance remains in the observed shapes. Nonetheless, TP factorization provides accessible visualization of the effect directions and allows to partially order the effect levels in each direction.

5.2 Cellular Potts Model Parameter Effects on Cell Form

The stochastic biophysical model proposed by Thüroff et al. (2019), a Cellular Potts Model (CPM), simulates migration dynamics of cells (e.g., wound healing or metastasis) in two dimensions. The progression of simulated cells is the result of many consecutive local elementary events sampled with a Metropolis-algorithm according to a Hamiltonian. Different parameters controlling the Hamiltonian have to be calibrated to match real live cell properties (Schaffer 2021). Considering whole cells, parameter implications on the cell form are not obvious. To provide additional insights, we model the cell form in dependence on four CPM parameters considered particularly relevant: the bulk stiffness $x_{i1}$, membrane stiffness $x_{i2}$, substrate adhesion $x_{i3}$, and signaling radius $x_{i4}$ are subsumed in a vector ${\mathbf{x}}_i$ of metric covariates for $i=1,\dots ,n$. Corresponding sampled cell outlines y_i were provided by Sophia Schaffer in the context of Schaffer (2021), who ran underlying CPM simulations and extracted outlines. Deriving the intrinsic orientation of the cells from their movement trajectories, we parameterize $y_i:[0,1]\to \mathrm{ }C$, clockwisely relative to arc-length such that $y_i(0)=y_i(1)$ points into the movement direction of the barycenter of the cell. With an average of $k=\frac{1}{n}{\sum }_{i=1}^n{k}_i\approx 43$ samples per curve (after sub-sampling preserving 95% of their inherent variation, as described in Volkmann et al. 2021, supplement), the evaluation vectors ${\mathbf{y}}_i\in \mathrm{ }{C}^{k_i}$ are equipped with an inner-product implementing trapezoidal rule integration weights. Example cell outlines are depicted in Figure S4, supplementary materials. The results shown below are based on cell samples obtained from 30 different CPM parameter configurations. For each configuration, 33 out of 10.000 Monte Carlo samples were extracted as approximately independent. This yields a dataset of $n=990=30\times 33$ cell outlines.

As positioning of the irregularly sampled cell outlines y_i , $i=1,\dots ,n$, in the coordinate system is arbitrary, we model the cell forms $[y_i]\in {\mathcal{Y}}_{/\text{Trl}+\text{Rot}}^{*}$. Their estimated overall form mean $[p]$ serves as pole in the additive model$$\begin{array}{c}\left[{\mu }_i\right]={\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{\left[p\right]}(h({\mathbf{x}}_i))\\ ={\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{\left[p\right]}(\sum _j{\beta }_j{x}_{ij}+\sum _j{f}_j(x_{ij})+\sum _{j\ne \ddot{\mathit{ȷ}}}{f}_{j\ddot{\mathit{ȷ}}}(x_{ij},x_{i\ddot{\mathit{ȷ}}}))\end{array}$$where the conditional form mean $[{\mu }_i]$ is modeled in dependence on tangent-space linear effects with coefficients ${\beta }_j\in T_{[p]}{\mathcal{Y}}_{/\text{Trl}\times \text{Rot}}$ and nonlinear smooth effects f_j for covariate $j=1,\dots ,4$, as well as smooth interaction effects $f_{j\ddot{\mathit{ȷ}}}$ for each pair of covariates $j\ne \ddot{\mathit{ȷ}}$. All involved (effect) functions are modeled via a cyclic cubic P-spline basis ${\left\{b_0^{(r)}\right\}}_r$ with 7 (inner) knots and a ridge penalty, and quadratic P-splines with 4 knots for the covariates x_ij equipped with a second-order difference penalty for the f_j and ridge penalties for interactions. Covariate effects are mean centered and interaction effects $f_{j\ddot{\mathit{ȷ}}}(x_j,x_{\ddot{\mathit{ȷ}}})$ are centered around their marginal effects $f_j(x_j),f_{\ddot{\mathit{ȷ}}}(x_{\ddot{\mathit{ȷ}}})$, which are in turn centered around the linear effects ${\beta }_j{x}_j$ and ${\beta }_{\ddot{\mathit{ȷ}}}{x}_{\ddot{\mathit{ȷ}}}$, respectively. Resulting predictor terms involve 69 (linear effect) to 1173 (interaction) basis coefficients but are penalized to a common degree of freedom of 2 to ensure a fair base-learner selection. We fit the model with a step-size of $\eta =0.25$ and stop after 2000 boosting iterations observing no further meaningful risk reduction, since no need for early-stopping is indicated by 10-fold form-wise cross-validation. Due to the increased number of data points and coefficients, the irregular design, and the increased number of iterations, the model fit takes considerably longer than in Section 5.1, with about 50 initial minutes followed by 8 hr of cross-validation. However, as usual in boosting, model updates are large in the beginning and only marginal in later iterations, such that fits after 1000 or 500 iterations would already yield very similar results.

Observing that the most relevant components point into similar directions, we jointly factorize the predictor as $\widehat{h}({\mathbf{x}}_i)=\sum _r{\xi }^{(r)}{\widehat{h}}^{(r)}({\mathbf{x}}_i)$ with TP factorization. The first component explains about 93% of the total predictor variance (Figure S3, supplementary materials), indicating that, post-hoc, a good share of the model can be reduced to the geodesic model $[{\widehat{\mu }}_i]={\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{[p]}({\xi }^{(1)}{\widehat{h}}^{(1)}({\mathbf{x}}_i))$ illustrated in Figure 3. A positive effect in the direction ${\xi }^{(1)}$ makes cells larger and more keratocyte/croissant shaped, a negative effect—pointing into the opposite direction—makes them smaller and more mesenchymal shaped/elongated. The bulk stiffness $x_{i1}$ turns out to present the most important driving factor behind the cell form, explaining over 75% of the cumulative variance of the effects (Figure S2, supplementary materials). Around 80% of its effect are explained by the linear term reflecting gradual shrinkage at the side of the cells with increasing bulk stiffness.

Fig. 3 Center: the main direction ${\xi }^{(1)}$ of the model illustrated as vectors pointing from the overall mean cell form $[p]$ (gray curve) to the form ${\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{[p]}({\xi }^{(1)})$ (filled dots), which are both oriented as cells migrating rightwards. Left: Effects of the bulk stiffness $x_{i1}$ into the direction ${\xi }^{(1)}$. A vertical line from 0, corresponding to $[p]$, to 1, corresponding to the full extent of ${\xi }^{(1)}$, underlines the connection between the plots and helps to visually asses the amount of change for a given value of $x_{i1}$. Right: The overall effect of $x_{i1}$ and membrane stiffness $x_{i2}$, comprising linear, smooth and interaction effects, as a 3D surface plot. The heat map plotted on the surface shows only the interaction effect $f_{12}^{(1)}(x_{i1},x_{i2})$ illustrating deviations from the marginal effects, which are of particular interest for CPM calibration.

5.3 Realistic Shape and form Simulation Studies

To evaluate the proposed approach, we conduct simulation studies for both shape and form regression for irregular curves. We compare sample sizes $n\in \{54,162\}$ and average grid sizes $k=\frac{1}{n}{\sum }_{i=1}^n{k}_i\in \{40,100\}$ as well as an extreme case with k_i = 3 for each curve but n = 720, that is, where only random triangles are observed (yet, with known parameterization over $[0,1]$). We additionally investigate the influence of nuisance effects and compare different inner product weights. While important results are summarized in the following, comprehensive visualizations can be found in Section S.5, supplementary materials.

Simulation design: We simulate models of the form $[\mu ]={\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{[p]}({\beta }_{\kappa }+f_1(z_1))$ with overall mean $[p]$, a binary effect with levels $\kappa \in \{0,1\}$ and a smooth effect of $z_1\in [-60,60]$. We choose a cyclic cubic B-spline basis with 27 knots for $T_{[p]}{\mathcal{Y}}_{/G}^{*}$, placing them irregularly at 1/27-quantiles of unit-speed parameterization time-points of the curves. Cubic B-splines with four regularly placed knots are used for covariates in smooth effects. True models are based on the bot dataset from R package Momocs (Bonhomme et al. 2014) comprising outlines of 20 beer (κ = 0) and 20 whiskey (κ = 1) bottles of different brands. A smooth effect is induced by the 2D viewing transformations resulting from tilting the planar outlines in a 3D coordinate system along their longitudinal axis by an angle of up to 60 degree toward the viewer ($z_1=60$) and away ($z_1=-60$) (i.e., in a way not captured by 2D rotation invariance). Establishing ground truth models based on a fit to the bottle data, we simulate new responses $[y_1],\dots ,[y_n]$ via residual resampling (Section S.5, supplementary materials) to preserve realistic autocorrelation. Subsequently, we randomly translate, rotate and scale $y_1,\dots ,y_n\in \mathcal{Y}$ somewhat around the aligned shape/form representatives to obtain realistic samples.

The implied residual variance $\frac{1}{n}{\sum }_{i=1}^n\Vert {ϵ}_i{\Vert }_i^2=\frac{1}{n}{\sum }_{i=1}^n{d}_i^2([y_i],[{\mu }_i])$ on simulated datasets ranges around 105% of the predictor variance $\frac{1}{n}{\sum }_{i=1}^n\Vert h({\mathbf{x}}_i){\Vert }_i^2=\frac{1}{n}{\sum }_{i=1}^n{d}_i^2([{\mu }_i],[p])$ in the form scenario and around 65% in the shape scenario. All simulations were repeated 100 times, fitting models with the model terms specified above and three additional nuisance effects: a linear effect $\beta z_1$ (orthogonal to $f_1(z_1)$), an effect f₂ of the same structure as f₁ but depending on an independently uniformly drawn variable z₂, and a constant effect $h_0\in T_{[p]}{\mathcal{Y}}_{/G}^{*}$ to test centering around $[p]$. Base-learners are regularized to 4 degrees of freedom (step-length $\eta =0.1$). Early-stopping is based on 10-fold cross-validation.

Form scenario: In the form scenario, the smooth covariate effect f₁ offers a particularly clear interpretation. TP factorization decomposes the true effect into its two relevant components, where the first (major) component corresponds to the bare projection of the tilted outline in 3D into the 2D image plane and the second to additional perspective transformations (Figure 4). For this effect, we observe a median relative mean squared error $\text{rMSE}({\widehat{h}}_j)={\sum }_{i=1}^n\Vert {\widehat{h}}_j({\mathbf{x}}_i)-h_j({\mathbf{x}}_i){\Vert }_i^2/{\sum }_{i=1}^n\Vert h({\mathbf{x}}_i){\Vert }_i^2$ of about 3.7% of the total predictor variance for small data settings with n = 54 and $k=100$ (5.9% with k = 40), which reduces to 1.5% for n = 162 (for both k = 40 and k = 100). It is typical for functional data that, from a certain point, adding more (highly correlated) evaluations per curve leads to distinctly less improvement in the model fit than adding further observations (compare, e.g., also Stöcker et al. 2021). In the extreme k_i = 3 scenario, we obtain an rMSE of around 15%, which is not surprisingly considerably higher than for the moderate settings above. Even in this extreme setting (Figure 4), the effect directions are captured well, while the size of the effect is underestimated. Rotation alignment based on only three points (which are randomly distributed along the curves) might considerably differ from the full curve alignment, and averaging over these sub-optimal alignments masks the full extend of the effect. Still, results are very good given the sparsity of information in this case. Having a simpler form, the binary effect ${\beta }_{\kappa }$ is also estimated more accurately with an rMSE of around 1.5% for n = 54, k = 100 (1.9% for k = 40) and less than 0.8% for n = 162 (for both k = 40 and k = 100). The pole estimation accuracy varies on a similar scale.

Fig. 4 Left: First (row 1) and second (row 2) main components of the smooth effect $f_1(z_1)$ in the form scenario obtained from TP factorization. Normalized component directions are visualized as bottle outlines after transporting them to the true pole (gray solid outline). Underlying truth (solid lines/shaded areas) are plotted together with five example estimates for n = 162 and k = 100 (black solid lines) and the extremely sparse k_i = 3 setting (gray dashed lines). Center: Conditional means for both bottle types with fixed metric covariate $z_1=0$ in the shape scenario with n = 54 and k = 40. Five example estimates (black solid outlines) are plotted in front of the underlying truth (shaded areas). Right: rMSE of shown example estimates (jittered diamonds) contextualized with boxplots of rMSE distributions observed in respective simulation scenarios.

Shape scenario: Qualitatively, the shape scenario shows a similar picture. For k = 40, we observe median rMSEs of 2.8% (n = 54) and 2.2% (n = 162) for $f_1(z_1)$, and 1.5% and 0.6% for the binary effect ${\beta }_{\kappa }$. For k = 100, accuracy is again slightly higher.

Nuisance effects and integration weights: Nuisance effects in the model where generally rarely selected and, if selected at all, only lead to a marginal loss in accuracy. The constant effect is only selected sometimes in the extreme triangle scenarios, when pole estimation is difficult. We refer to Brockhaus et al. (2017), who perform gradient boosting with functional responses and a large number of covariate effects with stability selection, for simulations with larger numbers of nuisance effects and further discussion in a related context, as variable selection is not our main focus here. Finally, simulations indicate that inner product weights implementing a trapezoidal rule for numerical integration are slightly preferable for typical grid sizes ($k=40,100$), whereas weights of $1/k_i$ equal over all grid points within a curve gave slightly better results in the extreme k_i = 3 settings.

All in all, the simulations show that Riemannian L₂-Boosting can adequately fit both shape and form models in a realistic scenario and captures effects reasonably well even for a comparably small number of sampled outlines or evaluations per outline.

6 Discussion and Outlook

Compared to existing (landmark) shape regression models, the presented approach extends linear predictors to more general additive predictors including also, for example, smooth nonlinear model terms and interactions, and yields the first regression approach for functional shape as well as form responses. Moreover, we propose novel visualizations based on TP factorization that, similar to FPC analysis, enable a systematic decomposition of the variability explained by an additive effect on tangent space level. Yielding meaningful coordinates for model effects, its potential for visualization will be useful also for FAMs in linear spaces and also beyond our model framework, such as we exemplarily illustrate for the nonparametric approach of Jeon and Park (2020) in Section S.8, supplementary materials.

Instead of operating on the original evaluations ${\mathbf{y}}_i\in C^{k_i}$ of response curves y_i as in all applications above, another frequently used approach expands y_i , $i=1,\dots ,n$, in a common basis first, before carrying out statistical analysis on coefficient vectors (compare Ramsay and Silverman (2005), Morris (2015), and Müller and Yao (2008) for smoothing spline, wavelet or FPC representations in FDA or Bonhomme et al. (2014) in shape analysis). Shape/form regression on the coefficients is, in fact, a special case of our approach, where the inner product is evaluated on the coefficients instead of evaluations (Section S.6, supplementary materials).

The proposed model is motivated by geodesic regression. However, in the multiple linear predictor, a linear effect of a single covariate does, in general, not describe a geodesic for fixed nonzero values of other covariate effects. Or put differently, ${\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{[p]}(h_1+h_2)\ne {\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{{\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{[p]}(h_1)}(h_2)\ne {\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{{\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{[p]}(h_2)}(h_1)$ in general. Thus, hierarchical geodesic effects of the form ${\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{{\hspace{0.17em}\mathrm{Exp}\hspace{0.17em}}_{[p]}(h_1)}(h_2)$, relevant, i.a., in mixed models for hierarchical/longitudinal study designs (Kim et al. 2017), present an interesting future extension of our model. Moreover, an “elastic” extension based on the square-root-velocity framework (Srivastava and Klassen 2016) presents a promising direction for future research, as do other manifold responses.

Acknowledgments

We sincerely thank Nadja Pöllath for providing carefully recorded sheep astragalus data and important insights and comments, and Sophia Schaffer for running and discussing cell simulations and providing fully processed cell outlines.

Supplementary Materials

Supplementary material with further details is provided in an online supplement.

Disclosure Statement

The authors report there are no competing interests to declare.

Additional information

Funding

We gratefully acknowledge funding by grant GR 3793/3-1 from the German research foundation (DFG) and support by the Open Access Publication Fund of Humboldt-Universität zu Berlin.