Probabilistic size-and-shape functional mixed models

Authors: fangyi  wang
karthik  bharath
oksana  chkrebtii
sebastian  kurtek
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 16:21:12

Abstract: The reliable recovery and uncertainty quantification of a fixed effect function $\mu$ in a functional mixed model, for modeling population- and object-level variability in noisily observed functional data, is a notoriously challenging task: variations along the $x$ and $y$ axes are confounded with additive measurement error, and cannot in general be disentangled. The question then as to what properties of $\mu$ may be reliably recovered becomes important. We demonstrate that it is possible to recover the size-and-shape of a square-integrable $\mu$ under a Bayesian functional mixed model. The size-and-shape of $\mu$ is a geometric property invariant to a family of space-time unitary transformations, viewed as rotations of the Hilbert space, that jointly transform the $x$ and $y$ axes. A random object-level unitary transformation then captures size-and-shape preserving deviations of $\mu$ from an individual function, while a random linear term and measurement error capture size-and-shape altering deviations. The model is regularized by appropriate priors on the unitary transformations, posterior summaries of which may then be suitably interpreted as optimal data-driven rotations of a fixed orthonormal basis for the Hilbert space. Our numerical experiments demonstrate utility of the proposed model, and superiority over the current state-of-the-art.

Main paper: 
 introduction: 
 Consider a sample of functions from the much-studied Berkeley growth data in Figure 1(a) , where growth rate curves from measurements on heights in centimeters of 54 girls and 39 boys from age 1 to 18 are plotted on a rescaled time axis. It appears that individuals experience different numbers of small and large growth spurts that differ in magnitude and timing. Any reasonable generative model for inference on distributional aspects of the functions will need to account for the fact that the sample size (n = 93) is smaller than the dimension of each observation, and it is thus common to consider tractable parametric probabilistic models. A popular choice is a functional mixed model √ γ has the same size-and-shape as f as described in Section 2. (d) PQRST complexes. (e) PQRST pattern: P wave (first max), QRS complex (sharp min-max-min) and T wave (last max) [Pham et al., 2023] .EQUATIONrandom function v i that represents deviations from µ; and, (iii) a zero-mean random measurement error process ϵ i . The growth rate function f i (t) of individual i at time t differs pointwise from an average µ(t) via a smooth random translation v i (t) observed with additive error ϵ i (t).A function with any number of peaks representing spurts may be generated under the model, and as such exact recovery of µ is not possible without further constraints on v i and ϵ i . The issue is exacerbated when f i in (1) is changed to f i •γ i by incorporating phase variability, or x-axis variability, via a time dilation/contraction object-level increasing function γ i : [0, 1] → [0, 1], which models variability in timings of the spurts. Even for the simpler model, f i • γ i = σµ + ϵ i , for a fixed scalar σ > 0, recovery of µ is possible only when the measurement error process ϵ i has a rank one covariance operator [Kurtek et al., 2011, Chakraborty and Panaretos, 2021] .The crux of the problem in reliable recovery of the population-level function µ lies in the fact that its topological features such as number of critical points and their nature (non-degenerate or otherwise), while preserved by γ i , may be altered due to addition of v i and ϵ i . Consequently, a growth spurt (peak) in the average µ at time t (or γ i (t)) may either be preserved or destroyed. From a modeling perspective, uncertainty quantification of µ in (1) within a Bayesian framework is unreliable in the presence of phase γ i , regardless of how informative the prior distributions on v i and γ i are. Figure 1(d) &(e) show an additional motivating dataset of PQRST complexes segmented from an electrocardiogram (ECG) signal and a typical PQRST pattern, respectively. Summarily, under the functional mixed model (1), the topological shape of µ may be irrevocably altered, and thus ordinate values of µ may not be reliably inferred through probabilistic modeling. The goal of the paper is to investigate if the situation may be salvaged by considering a complementary geometric property of µ, its size-and-shape, characterized by a joint transformation of its range and domain by the phase function γ.Contributions. We focus on sampling from and summarizing the posterior distribution of a fixed effect function µ in a functional mixed model with random object-level phase and amplitude components, without a finite-rank covariance assumption on the error process. We do not make any claims on efficiency of the computational algorithms, but emphasize that the goal is to thoroughly investigate the novel size-and-shape perspective in inference on µ. Thus, our contributions are as follows.1. We propose a mixed model in a functional Hilbert space for the size-and-shape of a squareintegrable fixed effect µ by considering an isometric action of the infinite-dimensional group of phase functions Γ. The geometric property of µ thus preserved is referred to as its size-and-shape (Section 2), which we demonstrate may be reliably inferred under the model. To our knowledge, this is the first paper to consider a Bayesian functional mixed effects model with an unrestricted form of phase variation, and employ the novel perspective of inferring the size-and-shape of µ.2. Informative priors on the phase functions γ regularize the posterior of µ, and sampling is assisted by exploiting the group structure of the phase functions while exploring the parameter space (Section 3 and Appendix D).3. Isometric action of γ engenders a unitary transformation of the Hilbert space; the class of unitary transformations indexed by γ constitutes rotations of coordinates of the Hilbert space. Upon expressing µ and v i in a suitable orthonormal basis of the Hilbert space, inference for γ i translates to an automatic identification of a data-optimal rotation of the chosen basis that best captures population-and object-level variations, and performs better when compared to an empirical functional principal component analysis (FPCA) basis (Appendices B and G). 2 4. We carry out extensive numerical experiments to investigate utility of the proposed model, and demonstrate that the posterior mean of µ under our model better captures the properties of µ than the estimate given by current state-of-the-art [Claeskens et al., 2021] (Section 4 and Appendix I).Related work. In absence of phase functions as part of the random effect, approaches have broadly focused on two types of dimension reduction techniques for estimating the fixed effect: (i) empirical orthonormal basis from FPCA [e.g., Yao et al., 2005] , and (ii) pre-specified basis functions [e.g., Rice and Wu, 2001 , Guo, 2002 , Chen and Wang, 2011 , Zhu et al., 2011 , Morris and Carroll, 2006 , Huo et al., 2023] . Some recent works have extended methodology to multivariate functional mixed effect models [Volkmann et al., 2023] , scalar on function regression [Liu et al., 2017] , spatial-temporal variation [Zhu et al., 2019] , and generalized functional mixed effect models [St. Ville et al., 2022] .In the presence of phase functions, the notable works are by Claeskens et al. [2021] and Raket et al. [2014] , where the former provides sufficient conditions for exact recovery of the fixed effect when the error process is not assumed to have phase variation; in our numerical experiments, we compare our results to the estimator proposed by Claeskens et al. [2021] , which represents the state-of-the-art.We are unaware of work in literature on Bayesian approaches to functional mixed models with unrestricted, nonparametric phase functions to jointly model amplitude and phase variations, and not sequentially following pre-processing via registration. Of relevance, however, is the work by Schiratti et al. [2017] , where manifold-valued curves with reparameterization variation are considered.2 Phase functions and size-and-shape-preserving transformationsWithout loss of generality, we assume that all functions are observed on a fixed domain[0, 1]. Let L 2 ([0, 1], R) = {f : [0, 1] → R| 1 0 |f (t)| 2 dt < ∞} (henceforth simply referred to as L 2) denote the representation space of interest, i.e., the space of real-valued square-integrable functions on [0, 1]. The space L 2 when endowed with the norm ∥f ∥ := [1 0 |f (t)| 2 dt] 1/2 coming from the inner product ⟨f, g⟩ := 1 0 f (t)g(t)dt is a Hilbert space. Additionally, let Γ = {γ : [0, 1] → [0, 1] γ(0) = 0, γ(1) = 1, γ > 0}denote the group of orientation-preserving diffeomorphisms ( γ is the time derivative of γ) of [0, 1]. In the functional data analysis literature, Γ is used to model phase variability in functional observations. Thus, elements of Γ will henceforth be referred to as phase functions. Note that Γ is a Lie group with composition (γ 1 , γ 2 ) → γ 1 • γ 2 as the group operation, where • is function composition. This implies that (i) Γ is an infinite-dimensional smooth manifold, (ii)EQUATIONiii) it contains the identity element γ id (t) = t such that γ • γ id = γ for any γ ∈ Γ, and (iv) for any γ ∈ Γ there exists γ -1 ∈ Γ such that γ • γ -1 = γ id . The group structure of Γ plays a pivotal role when defining prior and proposal distributions for phase functions; this is further elucidated in Section 3.2 and Appendix D. Importantly, the group Γ can act on the function space L 2 from the right, engendering maps L 2 × Γ → L 2 , in different ways, thus resulting in different notions of phase variation in functional data.1. Value-preserving action. The value-preserving mapping is defined via composition:f • γ, f ∈ L 2 , γ ∈ Γ.This action is referred to as the time warping of a function since f and f • γ traverse the same exact y axis values, but at different times (x axis values). It is commonly used for alignment or registration of prominent features in functional data, e.g., local extrema, a process traditionally referred to as amplitude-phase separation. 2. Area-preserving action. The area-preserving mapping is defined as(f • γ) γ, f ∈ L 2 , γ ∈ Γ.This action is commonly used for statistical analysis of probability density functions. 3. Norm-preserving action. The norm-preserving mapping is defined as(f, γ) := (f • γ) √ γ, f ∈ L 2 , γ ∈ Γ.An important property of this action is that it preserves the L 2 norm of a function: ∥f ∥ = ∥(f, γ)∥. The action has been profitably used in the problem of function alignment/registration, wherein a desideratum is a cost function invariant to a simultaneous action of the group Γ, or time warping [Srivastava and Klassen, 2016, Chapter 4] .γ : L 2 → L 2 , f → D γ (f ) := (f • γ)√ γ arising from the norm-preserving action for a fixed γ ∈ Γ is a surjective isometry, and hence unitary. Thus, D γ , for a fixed γ ∈ Γ, is an infinite-dimensional rotation in the Hilbert space L 2 , and the group D := {D γ , γ ∈ Γ} of rotations of L 2 plays a prominent role in classical white noise calculus [Hida, 2015] .A general space-time diffeomorphism (t, x) → (σ 1 (t, x), σ 2 (t, x)) preserves the topology of the graph {(t, f (t))} of f viewed as a subset of [0, 1] × R. The operator D γ is associated with a special space-time diffeomorphism with σ 1 (t, x) = γ(t) and σ 2 (t, x) = x γ(t). In this sense, the size and shape of f , as it relates to its norm and its graph, is preserved by operators in D.To better understand how D γ transforms f , consider Figure 1 (b)&(c). As seen with the red curve in (c), the value-preserving action of γ preserves the image of t → f (t), where ordinate values of f are relabeled in an order-preserving way and no new values are created/destroyed. A more classical notion of the shape of f is thus preserved under the value-preserving action in that the topology of the level sets of f is unaltered. In contrast, the norm-preserving action shown with the yellow curve in (c) alters the image of t → f (t), and may create new critical points, but in a way such that its size-and-shape in the sense alluded to above is preserved. As such, the norm-preserving action is size-and-shape preserving, i.e., f and D γ (f ) for any γ ∈ Γ are equivalent in their size-and-shape.Use of the operator D γ is particularly useful for modeling phase variation in model ( 1):(i) Dimension reduction of the fixed and random effects is implemented via an orthonormal basis system, and operators in D allow us to rotate these basis systems to better align with observed data (see Appendix B). This provides modeling flexibility that is often needed due to the use of a finite, and potentially small, number of basis functions to represent these model components. (ii) The norm preserving action can be viewed as a combination of value-preserving warping and an associated local scaling of function values, i.e., when γ(t) > 1 ( γ(t) < 1) the function value f (γ(t)) is warped to the right (left) relative to f (t) and additionally rescaled by a factor of γ(t). (iii) The equivalence class of functions having the same size-and-shape under the norm-preserving action is 'larger' than the one under the value-preserving action in the following sense. The map π : L 2 → L 2 / ∼ takes a function to its equivalence class determined by equivalence relation ∼. It can be shown that the measure of an equivalence class, under the pushforward of a non-degenerate measure with support in L 2 under π, is larger when ∼ pertains to norm-preserving action as opposed to value-preserving one (see Example 6.5.2 in Bogachev [1998] for an example that illustrates this idea). In practice, the larger equivalence class under the norm-preserving action helps with reliable recovery of the size-and-shape of the fixed effect µ, in contrast to the more traditional notion of shape of µ under the value-preserving action, wherein the additive linear term v i + ϵ i more easily alters the shape of µ.3 Bayesian size-and-shape functional mixed modelThe operator D γ related to the norm-preserving action of Γ is an isometry of L 2 . The use of Euclidean isometries arising from translations and rotations as size-and-shape preserving has a long history in statistical shape analysis of landmark configurations [e.g., Dryden and Mardia, 2016] . To elaborate, let X i ∈ R K×2 denote a set of K landmarks, i.e., X i is a matrix that contains the coordinates of K points in R 2 , and SO(2) = {R ∈ R 2×2 |R T R = RR T = I, det(R) = 1} denote the rotation group with matrix multiplication as the group operation. Consider the following perturbation (generative) model for the size-and-shape of an object represented via a landmark configuration X i :EQUATIONwheret i ∈ R 2 is a random translation, 1 ∈ R K is a vector of 1s, R i ∈ SO(2) is a random rotation,M ∈ R K×2 is a fixed template, and E i ∈ R K×2 is the error. Using this model, it is of primary interest to estimate the fixed effect (size-and-shape) M based on observed landmark configurations X 1 , . . . , X n . Transformation of the model components M and E i using rotations R i (and translations) in order to align them to the observation X i may be viewed as an isometric, or norm-preserving, transformation of the coordinate system for (M + E i ), so that only the size-and-shape of M can be reliably recovered, but not M [Lele and McCulloch, 2002] .Inspired by size-and-shape analysis of landmark data, and in contrast to existing works using the valuepreserving action [Raket et al., 2014 , Claeskens et al., 2021] , our proposed functional mixed effects model utilizes the norm-preserving action of Γ on L 2 . Analogous to the model in (2), for a function f i ∈ L 2 , one can define a functional perturbation modelf i = D γi (µ + ϵ i ) = [(µ + ϵ i ) • γ i ] √γi for a template function µ with random phase functions γ i , wherein the L 2 coordinate system of the model components µ and ϵ i is aligned to the coordinate system of the function f i via a rotation using γ i .To the fixed effect, or signal plus noise model, one can additionally introduce an object-level random effect to increase modeling flexibility, resulting in:f i = D γi (µ + v i + ϵ i ) = [(µ + v i + ϵ i ) • γ i ] γi .(3) Thus, γ i and v i act as size-and-shape preserving and altering random effects, respectively. As in the landmark case, the primary goal of interest is to estimate µ and quantify its uncertainty using independent observations f 1 , . . . , f n , and it is possible to reliably infer only the size-and-shape of µ.While the model in (3) is written in terms of infinite-dimensional objects, functional data is often observed under a common discretization at time pointst 1 = 0 < t 2 < • • • < t T -1 < t T = 1.The assumption of a common discretization for all functional observations is adopted for clarity and succinctness in presenting the proposed Bayesian model; it can be easily relaxed to accommodate more general scenarios. The model is quite flexible and more realistic assumptions (e.g., arbitrary error dependence structure such as Matérn covariance, or non-Gaussian) may easily be incorporated at some computational cost. We instead focus on the phase component, and demonstrate the utility of the size-and-shape perspective in inferring geometric properties of µ. Let f i = (f i (t 1 ), . . . , f i (t T )) ⊤ ∈ R T , i = 1, . . . ,EQUATIONLetting {ϕ k } ∞ k=1 and { φk } ∞ k=1 denote two orthonormal basis systems for L 2 , we represent the model components µ and v i as linear combinations of basis functions. For dimension reduction, we define µ := B f k=1 a k ϕ k and v i := Br k=1 c i,k φk , i.e., the two basis sets are truncated to B f and B r basis functions for the fixed and (size-and-shape altering) random effects, respectively. We assume that the random effect coefficients are independent and identically normally distributed (iid) c i,k ∼N (0, σ 2 c ). The discretized error is assumed to be iid ϵ i (t j )∼N (0, σ 2 ). 
choice of basis functions and likelihood: 
 Specifications of µ and v i require appropriate orthonormal basis functions, and we consider modified Fourier basis or the B-spline basis for µ, and the B-spline basis for v i . The modified Fourier basis contains the following elements { √ 3t, √ 3(1 -t), √ 2 cos(2πjt), √ 2 sin(2πjt); j = 1, 2, . . . ; t ∈ [0, 1]} and is subsequently orthonormalized via the Gram-Schmidt procedure under the L 2 metric; larger values of j yield basis functions with finer harmonics over the domain [0, 1] . This basis is used to specify the fixed effect function µ only as it is effective at representing global periodic trends and oscillations. On the other hand, the B-spline basis is defined locally via piecewise polynomials, and is therefore better at capturing local variation and finer function features; we also orthonormalize the B-spline basis via the Gram-Schmidt procedure under the L 2 metric. In some cases where the underlying µ is expected to exhibit more local features, we use the B-spline basis rather than the modified Fourier basis. Finally, one needs to choose the number of basis functions to model µ and v i . This choice depends on the application of interest. For example, one may use a larger number of basis functions for µ when the observations have a complex shared global structure. Effects of misspecifying the number of basis functions for µ and v i are studied in Appendix H.The main inferential tasks of interest are to estimate and assess uncertainty in (i) mean size-and-shape µ via the coefficient vector a = (a 1 , . . . , a B f ) ⊤ , (ii) variance of the size-and-shape altering random effect coefficients σ 2 c , and (iii) error variance σ 2 . Thus, to simplify the inference, we marginalize the likelihood with respect to the size-and-shape altering random effect. Let Φ i ∈ R T ×B f denote the matrix of evaluations of the fixed effect basis, whose (j, k)th entry is given by (ϕ k • γ i )(t j ) γi (t j ). One can view Φ i as a discretization of the rotated coordinate system, via the norm-preserving action that uses γ i , for µ. Similarly, let Φi ∈ R T ×Br be the matrix of evaluations of the random effect basis, whose (j, k)th entry is ( φk•γ i )(t j ) γi (t j ). Let a ∈ R B f and c i ∈ R Br , i = 1, .. . , n be the vectors of basis coefficients for µ and v i , i = 1, . . . , n, respectively. Finally, let ϵ γi i ∈ R T , i = 1, . . . , n be the vectors of γ i -transformed observation errors, with entries (ϵ i • γ i )(t j ) γi (t j ), j = 1, . . . , T that are independently distributed as N (0, σ 2 γi (t j )), conditional on γ i ; the additional time-dependent scaling of the error variance comes from the norm-preserving action. Using this simplified notation, we can rewrite the model in (4) in matrix form asEQUATIONwhich may be usefully compared to the size-and-shape perturbation model ( 2) for landmark configurations. With MVN used to denote the multivariate normal distribution, the distribution of f i , conditional on the vector of coefficients c i , isf i |c i ∼ MVN(Φ i a + Φi c i , σ 2 diag( γi (t j ))), where diag( γi (t j )) is a T × T diagonal matrix whose jth diagonal element is γi (t j ). Then, the following marginal distribution of f i is used to define the likelihood function (see Appendix C for full derivation):EQUATIONThe model parameters that need to be estimated are the fixed effect coefficients a, random effect variance σ 2 c , variance of observation error σ 2 , and individual phase functions γ i . For a, σ 2 c and σ 2 , we use weakly informative prior distributions:a ∼ MVN(0, 10000I B f ), σ 2 c ∼ IG(0.01, 0.01), σ 2 ∼ IG(0.01, 0.01) (IG is the inverse-gamma distribution).The prior distribution over the space of phase functions Γ, to model the shape-and-size preserving random effect, can be specified in different ways [e.g., Telesca and Inoue, 2008 , Bigot, 2013 , Lu et al., 2017] . In this work, we use two tractable prior distribution models on Γ to guard against confounding of inference between {γ i } and µ: (i) a one-parameter family; (ii) a nonparametric finite-dimensional family compatible with the time discretization. Importantly, the prior models rely on the group structure of Γ.• Prior Model 1 (PM1) on Γ. Each phase function γ i is defined via a single parameter α i (Section 5.2 in Srivastava and Klassen [2016] ) asEQUATIONThe phase functions defined in ( 7) form a one-dimensional subset of Γ, making posterior inference more tractable with significant reductions in computational complexity. This subset contains phase functions with γ(t) ∈ (0, 2) for all t ∈ (0, 1) and no inflection points. At the same time, we sacrifice flexibility in terms of the allowed rotations of L 2 . The α i , i = 1, . . . , n are assumed to be independent a priori following the Uniform(-1, 1) distribution. • Prior Model 2 (PM2) on Γ. A more flexible nonparametric prior model utilizes a point processbased prior distribution related to the Dirichlet process on Γ [Bharath and Kurtek, 2020] . Under discretized time, each phase function may be represented by a finite number of successive increments p(γi ) = (γ i (t 2 ) -γ i (0), . . . , γ i (t j ) -γ i (t j-1 ), . . . , γ i (1) -γ i (t Tγ -1 )) ∈ R Tγ -1 ; as T γ → ∞,the resulting prior has dense support in Γ. Then, the finite-dimensional prior distribution is placed on the vector of phase increments,EQUATIONwhere t = (t 2 , t 3 -t 2 , ..., 1 -t Tγ -1 ) is defined via T γ consecutive time points on [0, 1] and θ γ is a precision parameter. First, the time points defining t can be different from the time points at which the functional observations were recorded. In fact, we usually choose T γ to be relatively small, e.g., five or seven, to simplify the prior model. Second, a large value for θ γ regularizes the phase functions toward γ id . We use θ γ = 30 in all of our numerical experiments. The resulting γ i is a piecewise linear function with changes in the slope γi at t 2 , . . . , t Tγ -1 . This prior distribution is more flexible than the one in PM1. However, depending on the chosen number of discretization points T γ , the dimension of the phase parameter space can be much larger in this case, making posterior inference more challenging. Interested readers can refer to Matuk et al. [2022] for more details behind this choice of prior distribution on phase functions in the context of Bayesian modeling for functional data.Posterior inference on all parameters is conducted via Markov chain Monte Carlo (MCMC) sampling using the Metropolis-Hastings algorithm. To efficiently explore the parameter space, we use adaptive proposal distributions. We monitor MCMC convergence using standard diagnostic plots, e.g., trace plots and autocorrelation plots. Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [NA] Justification: Our paper does not include theoretical results, but we make modeling assumptions with appropriate justifications. Appendices C and D include detailed derivations for the marginal likelihood used in our model and the Metropolis-Hastings acceptance ratio for MCMC sampling.Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
numerical experiments: 
 We present posterior inference results from the model described in Section 3 for simulated and real data. Throughout, we compare our results to those generated by warpMix, a state-of-the-art frequentist functional mixed model [Claeskens et al., 2021] ; we use the default parameter settings in warpMix. In addition, in Appendix A, we compare our results to those generated by a state-of-the-art Bayesian approach [Cheng et al., 2016] , which models functions under the popular square-root velocity transformation [Srivastava et al., 2011b] as realizations of a Gaussian process centered at a mean function µ. Phase variation is incorporated via the value-preserving action of Γ with a prior model that is the same as PM2 in our model.The main object of interest for posterior inference is µ. Note however that, as discussed earlier, µ and D γ (µ) = (µ • γ) √ γ for any γ ∈ Γ are equivalent in terms of their size-and-shape. Thus, we first center all posterior samples for µ, via an average of the estimated posterior means of object-level phase functions, to obtain size-and-shape representatives in their equivalence classes under the norm-preserving action; this centering is similar in spirit to the orbit centering in Srivastava et al. [2011b] . Assuming that we have N posterior samples of each γ i , i = 1, . . . , n, we first compute γ := 1/(nN )n i=1 N j=1 γj i, where γj i is the jth posterior sample of the phase function γ i ; due to convexity of Γ, note that γ ∈ Γ. We then compute the centered posterior samples of µ, (μ j • γ) γ, which may be visualized directly or used to estimate posterior summaries, e.g., pointwise posterior mean and credible interval. In some cases, we also visualize the posterior samples and pointwise posterior mean for a phase function γ i . For all examples in this section, we use N = 100, 000 with a burn-in period of 200, 000 iterations; for visualization of posterior samples for µ and γ i , we use a uniform subsample of size 1, 000 to ensure that all plots are easily readable.Example 1: data generated from our model. We first consider an example based on data simulated from model in (6). We use B f = 6 modified Fourier basis functions for µ and B r = 6 B-spline basis functions for each v i . The ground truth variances are σ 2 c = 0.25 and σ 2 = 0.1. Then, to generate the data, we sample a ∼ MVN(0, I B f ), and consider two cases for phase functions based on PM1 and PM2: (i) α i iid ∼ Uniform(-1, 1), and (2) p(γ i ) iid ∼ Dirichlet(30t), t = (0, 0.25, 0.5, 0.75, 1). The sample size in this simulated example is n = 30. model is effective at recovering the underlying size-and-shape of µ, which contains two local minima and maxima. On the other hand, the warpMix estimate only contains one local minimum and maximum, and is clearly an inaccurate representation of the size-and-shape of µ. Panels (c) and (d) show histograms of posterior samples for the variance parameters σ 2 and σ 2 c , respectively; the ground truth values are in black and the estimated posterior means in red. In both cases, we slightly overestimate both variance parameters. Finally, in (e), we show estimation results for a phase function γ i corresponding to a randomly chosen observation. As before, posterior samples are shown in blue, posterior mean in red and the ground truth in black. In both cases (phase function generated from PM1 or PM2), we reliably recover the underlying ground truth.Example 2: data generated from warpMix model. Next, we consider a more challenging scenario for the proposed model and specifically focus on a comparison to warpMix. In this case, phase functions and random effect functions are generated using the warpMix model; we set σ 2 c = 0.25 and σ 2 = 0.0001. Comparison of estimation accuracy is based on three different fixed effect functions: µ 1 (t) = {sin(3πt) + 3πt}/4, µ 2 (t) = exp -(t-0.25) 2 /0.04 + exp -(t-0.75) 2 /0.02 and µ 3 (t) = cos(2πt + π/2), t ∈ [0, 1]. To generate the data, we use the value-preserving action as in the warpMix specification, and not the norm-preserving action that is utilized in our model.To specify our models, we use (i) B f = 6 modified Fourier basis functions or B-spline basis functions for µ with PM1 or PM2 for phase functions with T γ = 7, and (ii) B r = 6 B-spline basis functions for each v i . In total, we consider four Bayesian models, 1-F, 1-B, 2-F and 2-B, where the number indexes the prior model on phase and the letter indexes the basis used to model µ.First, in Figure 3 , we present estimation results for (a) µ 1 , (b) µ 2 and (c) µ 3 based on Model 2-B. The ground truth is in black, centered posterior samples in blue, centered posterior mean in red, and warpMix estimate in yellow. Visually, both warpMix and Model 2-B are effective at recovering the underlying fixed effect functions. To quantitatively assess estimation accuracy, we adopt the evaluation criterion used in Claeskens et al. [2021] :∆ µ := T -1 j=1 [μ(t j ) -µ(t j )] 2 (t j+1 -t j ).The μ for our models is defined as the centered posterior mean. Table 1 reports the results with best performance highlighted in bold. In all three cases, one of our Bayesian models yields the lowest estimation error. For µ 1 , only Model 1-B outperforms warpMix, which is not surprising since µ 1 has the simplest structure. Nonetheless, the other three Bayesian models are competitive as well. For µ 2 , all four of our models significantly outperform warpMix, with Model 2-F yielding an error reduction of 82% as compared to warpMix. For µ 3 , the two models with PM2 on Γ outperform warpMix suggesting the need for flexible phase functions in this case. These results show that the proposed models are effective in estimating µ even when the underlying data generating process uses value-preserving warping, and supports claim (iii) at the end of Section 2. We now consider application of the proposed modeling framework to (i) Berkeley growth rate functions (n = 93) [Srivastava et al., 2011b] (Figure 1 (a)), and (ii) PQRST complexes (n = 40) [Kurtek et al., 2013] (Figure 1(d) ). Our primary interest in the Berkeley data lies in estimating an average pattern of growth spurts. However, the number and magnitudes of growth spurts for each child may differ quite markedly from those in the average growth rate function. Majority of the PQRST complexes exhibit similar pattern of local extrema and we are interested in inferring the average pattern, while accounting for variability in magnitudes of the extrema across observations. A functional mixed effects model is thus appropriate for both these data settings to reliably estimate the size-and-shape of µ.For Berkeley data, we use modified Fourier basis for µ with B f = 6, B-spline basis for each v i with B r = 6, and PM1 on Γ; for PQRST data we use B-spline basis for µ with B f = 12, B-spline basis for each v i with B r = 6 to better model the sharp local features of the QRS complex (Figure 1(e) ), and the PM2 model on Γ to allow for inflection points in estimated phase.Row 1 in Figure 4 shows estimation results for the Berkeley data. In (a), we show centered posterior samples of µ in blue with the centered posterior mean in red. We uncover two growth spurts, a small initial one followed by a larger pubertal one. This agrees with previous literature that has considered this data [Srivastava et al., 2011b] . The marginal posterior uncertainty for µ is very small throughout the domain. We also show the warpMix estimate, which was only able to recover one small growth spurt. Panels (b) and (c) show histograms of posterior samples of σ 2 and σ 2 c , respectively. Panel (d) displays the posterior samples (blue) and posterior mean (red) of a phase function for a randomly chosen observation. Panel (e) shows the observation corresponding to panel (d) in black. In addition, the blue functions are D γj i (μ j ), the posterior samples of µ rotated using corresponding posterior samples of the phase function, where i, j index the observation and posterior sample, respectively. Note that the blue functions fit the observed black function well.Row 2 in Figure 4 shows estimation results for the PQRST data. Panels (a)-(e) are the same as in the previous description. The estimated size-and-shape of µ resembles a PQRST complex as desired and contains sharp features that are representative of the observed data. Further, posterior uncertainty appears greater along the P and T waves than the QRS complex. The warpMix model failed to yield an estimate in this case, potentially due to lack of modeling flexibility in capturing the QRS complex and phase variation. The estimated phase function in (d) contains an inflection point motivating our use of PM2 on Γ. Finally, panel (e) shows that the rotated posterior samples of µ fit a randomly chosen observation well. Appendix I shows posterior mean and 95% credible interval estimates, as well as warpMix estimates, of µ for five additional datasets.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: We have a clear description of all model parameters, prior distributions and proposal distributions for all numerical experiments presented in Section 4. We provide a detailed MCMC algorithm in appendix E. We also included code as supplementary material and provide step-by-step instructions on how to reproduce the table and figures presented in the main part of the paper.• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
discussion: 
 Numerical experiments in Section 4 and Appendices A and I demonstrate benefits of the proposed mixed model in recovering the size-and-shape of a fixed effect function µ, while outperforming current state-of-the-art. What is lacking is theoretical support for the same, and this is work in progress.To evolve the MCMC algorithm in MATLAB R2021a for 300, 000 iterations yielding 100, 000 posterior samples after burn-in, on a computing server with 6 parallel Intel(R) Xeon(R) CPUs with 20GB of memory, the computing time is approximately 93 and 111 minutes, respectively, under PMs 1 and 2 on Γ, based on n = 30 functions discretized at T = 50 points. There is room for improvement of efficiency in the MCMC computations, and alternatives may be explored.We specify the number of basis functions for the fixed effect µ and the size-and-shape altering random effect v i a priori. Alternatively, one could treat the number of basis functions B f and B r as random and estimate them. This, however, would require more advanced MCMC algorithms for posterior inference. Further, we use the modified Fourier basis and B-spline basis in our model. Other orthonormal basis functions could also be used, but this choice is not crucial for reliable recovery of the size-and-shape of µ since the norm preserving action D γ rotates the basis system toward the data, allowing us to learn a data-driven basis for µ. An alternative approach would be to directly learn an appropriate subspace for µ, which we plan to consider in future work.Exciting and novel extensions of the proposed mixed modeling framework are readily available to handle more complex functional data. The proposed model may be easily modified to handle sparsely/irregularly sampled and fragmented, or partially observed, functional data [Matuk et al., 2022] . The norm-preserving action of Γ may be used to perform inference for the size-and-shape of a fixed effect parameterized open/closed curve µ : [0, 1] → R d , d > 1 [Srivastava et al., 2011a , Kurtek et al., 2012] . Finally, for two-dimensional parameterized surfaces f : D ⊂ R 2 → R 3 , a similar norm-preserving action of the reparameterization diffeomorphism group with elements γ : D → D [Jermyn et al., 2017 ] may be used to infer the mean size-and shape of a fixed effect surface.Functional data is arising as a common object in various applications, including computer vision and biomedical imaging, and functional data models are becoming increasingly important in machine learning [Rao and Reimherr, 2023a,b] . The ideas we have explored in the proposed framework can be applied more broadly to other regimes. For instance, employing mixed models with random effects to better model correlated input data in neural networks is fast gaining traction [Simchoni and Rosset, 2023] . There is also increasing interest in the use of geometry and invariance to nuisance transformations in neural networks [Bronstein et al., 2021] . Our framework provides understanding for the type of signal that can be recovered from complex input data in the presence of nontrivial symmetries and geometric information, which offers a new perspective on incorporating geometric constraints in various types of data settings and models.We include the following appendices as support for the content in the main paper: Table 2 reports a comprehensive quantitative evaluation, in terms of estimation accuracy for the fixed effect function µ, on the five simulated datasets used in Section 4.1. Rows 1-2 consider data simulated from our model under Prior Models 1 (PM 1) and 2 (PM 2) for phase functions. Rows 3-5 consider data simulated using the warpMix model with default parameter values. We compare estimation results produced using our model (columns Model 1-F through Model 2-B, where the number indicates the prior model on phase and the letters F (Fourier) or B (B-spline) correspond to the type of basis used to model the fixed effect function µ) to those produced using warpMix [Claeskens et al., 2021] and the Bayesian model proposed in Cheng et al. [2016] (BRFC). As seen in the table, our model outperforms warpMix and BRFC in all of these simulation scenarios. shows the first six modified Fourier basis functions. In (c)&(d), we show the same basis functions after applying the value-preserving and norm-preserving actions using the phase functions in (a). Finally, in (e), we provide an example function formed using a linear combination of the basis functions in (b) (blue), and the same linear combinations but using basis functions in (c) (red) and (d) (yellow), respectively. Note that the norm-preserving action, when combined with an orthonormal basis system, provides more flexibility in modeling. The original blue function as well as the transformed red function contain four extrema in both rows. On the other hand, the yellow functions have five extrema (there is a small local maximum in the yellow function in the bottom row near t = 0.1). In another example, we study the magnitude of residuals when data is (i) projected onto a fixed number of modified Fourier basis functions, and (ii) projected onto the same number of modified Fourier basis functions followed by optimization over phase functions under the norm-preserving action. We vary the number of basis functions from 1 to 30 and average the residuals over all n = 93 Berkeley growth rate functions. Given a function f ∈ L 2 , its projection onto the basis is f = B j=1 âj ϕ j , where âj = 1 0 f (t)ϕ j (t)dt, j = 1, . . . , B, B is the number of basis functions used in the projection, and ϕ j , j = 1, . . . , B are the basis functions. The residual for (i) is defined as ∥f -f ∥, while the residual for (ii) is min γ∈Γ ∥f -( f • γ) √ γ∥. Overall, we expect the average residual from (ii) to be lower than the average residual from (i). However, an interesting aspect of this experiment is to understand how much more expressive the basis is when an additional rotation via the norm-preserving action is allowed, especially when very few basis elements are used in the projection. The results are presented in Figure 6 . The average residuals, for B = 1, . . . , 30, based on (i) and (ii) are shown in blue and red. As expected, (ii) yields lower average residuals, with a very large gap for B = 1, . . . , 10. The inclusion of the norm-preserving action becomes less valuable when a large number of basis elements is used. This is also expected since a data-driven rotation of the coordinate system becomes less crucial when more and more coordinates, defined with respect to the basis, are included. This shows that, when a small number of basis elements is used in the specification of µ in our model, the norm-preserving action allows for more flexible and efficient modeling, by rotating the coordinate system for µ to the observed data. Conditional on the coefficients of the size-and-shape altering random effect, f i follows a multivariate normal distribution:f i |c i ∼ MVN(Φ i a + Φi c i , σ 2 diag( γi (t j ))).The expected value and variance of f i areE[f i ] = E[E[Φ i a + Φi c i |c i ]] = E[Φ i a + Φi c i ] = Φ i a + Φi E[c i ] = Φ i a. and Var[f i ] = E[Var[f i |c i ]] + Var[E[f i |c i ]] = σ 2 diag( γi (t j )) + Var[Φ i a + Φi c i ] = = σ 2 diag( γi (t j )) + Var[ Φi c i ] = σ 2 diag( γi (t j )) + σ 2 c Φi ΦT i .Thus, the marginal distribution of f i isEQUATIONWe use the following proposal distributions in the Markov chain Monte Carlo algorithm:1. Fixed effect coefficients: a can ∼ N (a cur , Σ a ).2. Variance of error process: (σ 2 ) can ∼ TN((σ 2 ) cur , τ 2 σ , 0, ∞). 3. Variance of size-and-shape altering random effect: (σ 2 c ) can ∼ TN((σ 2 c ) cur , τ 2 σc , 0, ∞). 4. Size-and-shape preserving random effect (phase functions) under Prior Model 1:α can i ∼ Uniform(α cur i -δ, α cur i + δ).5. Size-and-shape preserving random effect (phase functions) under Prior Model 2:γ can = γ cur • γ,where p(γ) ∼ Dirichlet(αt).TN stands for the truncated normal distribution. Notably, the proposal in 5. for phase functions utilizes the group structure of Γ. The proposal covariance matrix Σ a is adapted during the burn-in period according to the empirical correlation matrix of the samples of a. The proposal variances for the error process and the size-and-shape altering random effect coefficients, τ 2 σ and τ 2 σc , respectively, as well as the parameter δ and precision parameter α for the size-and-shape preserving random effect under Prior Models 1 and 2, respectively, are adapted according to the acceptance rate during the burn-in period.Let L(•|•) denote the likelihood function, and π(•|•), p(•) and q(•|•) denote the posterior, prior and proposal distributions. Then, the general form of the MH acceptance ratio for a parameter β k , an element of a parameter vector β, with observed data denoted by f , isEQUATION)where β (-k) is the vector of parameters excluding β k . Note that, for the model proposed in this manuscript β = (a, σ 2 c , σ 2 , γ).Fixed effect coefficients a. Let β -a denote the vector of parameters excluding a. Since the proposal is a random walk centered at the current value, it is symmetric. Thus,EQUATION)Variance of error process, σ 2 , and variance of shape-and-size altering random effect, σ 2 c . Since σ 2 and σ 2 c have the same prior and proposal distributions, we derive the acceptance ratio ρ σ 2 only, and note that ρ σ 2 c is the same. The proposal for these two variance parameters is a truncated normal distribution. Thus,q((σ 2 ) cur |(σ 2 ) can ) = 1 τ ϕ (σ 2 ) cur -(σ 2 ) can τ Φ (σ 2 ) can τ 1 {0<(σ 2 ) cur <∞} , q((σ 2 ) can |(σ 2 ) cur ) = 1 τ ϕ (σ 2 ) can -(σ 2 ) cur τ Φ (σ 2 ) cur τ 1 {0<(σ 2 ) can <∞} ,where ϕ(•) and Φ(•) denote the probability density and cumulative distribution functions for the standard normal distribution. As a result,EQUATIONSize-and-shape preserving random effect γ i under Prior Model 1. Recall that, under this prior model, γ i (t) = t + α i t(1 -t), α ∈ (-1, 1), t ∈ [0, 1]. The Unif(-1, 1) prior on α i assigns zero mass to α can i / ∈ (-1, 1) resulting in automatic rejection. In other cases, the proposal Unif(α cur i -δ, α cur i + δ) is symmetric. Thus,EQUATION)Size-and-shape preserving random effect γ i under Prior Model 2. SinceEQUATION)we focus on the derivation ofq(γ cur i |γ can i ) q(γ can i |γ cur i) . We omit the subscript i to simplify notation. Letp(γ) = (γ(t 2 ) -γ(0), γ(t 3 ) -γ(t 2 ), . . . γ(t Tγ -1 ) -γ(t Tγ -2 ), γ(1) -γ(t Tγ -1 )) = =: (∆ 1 (γ), ..., ∆ Tγ -1 (γ)) =: ∆(γ) (15)The proposal distribution is ∆(γ) ∼ Dirichlet(αt), where t = (t 2 , t 3 -t 2 , . . . , 1 -t Tγ -1 ) = (t (1) , . . . , t (Tγ -1) ), and has densityEQUATIONNext, we derive the density of ∆(γ can ) given ∆(γ) and γ cur , which appears in the denominator of the acceptance ratio. Since γ = γ -1 cur • γ can , we haveEQUATIONTherefore, the density of ∆(γ can ) given ∆(γ) and γ cur is the density of ∆(γ), i.e., the density of Dirichlet(αt), multiplied by the determinant of the Jacobian of g. The Jacobian of g is given byEQUATIONand is a lower triangular matrix. Then, the determinant of this Jacobian isEQUATION)The numerator of the acceptance ratio contains the density of ∆(γ cur ), given ∆(γ -1 ) and γ can , which is the density of ∆(γ -1 ) multiplied by the determinant of the Jacobian of h := h(∆(γ cur ); γ-1 ), the transformation from ∆(γ cur ) to ∆(γ -1 ). This determinant is derived in the same way as the determinant of the Jacobian of g. Finally, we haveEQUATIONThe detailed MCMC algorithm is given in Algorithm 1.Input: Data: f i = f i (t j ), i = 1, . . . , n, j = 1, . . . , T ; prior hyperparameters: τ 2 a (prior variance for a), a σ , b σ (shape and scale for σ 2 ), a σc , b σc (shape and scale for σ 2 c ), θ γ , t (concentration and discretization in Prior Model 2 for γ); proposal hyperparameters: Σ a (proposal covariance for a), τ 2 σ (proposal variance for σ 2 ), τ 2 σc (proposal variance for σ 2 c ), δ (proposal for α in Prior Model 1 for γ), α, t (concentration and discretization in Prior Model 2 proposal for γ); number of Figure 7 shows trace plots of 100, 000 MCMC iterations after the burn-in period for all model parameters. The first two rows correspond to the simulated data examples considered in Figure 2 in Section 4.1. The third and fourth rows correspond to the two real data examples considered in Figure 4 Section in 4.2 (third row: Berkeley; fourth row: PQRST). Panels (a) and (b) show trace plots for the first two coefficients (first two entries in the vector a of the fixed effect function µ. Panels (c) and (d) show trace plots for σ 2 and σ 2 c , respectively. Finally, panel (e) shows trace plots for the parameter α i when Prior Model 1 was used, or the parameter values γ i (t 2 ) (blue), γ i (t 3 ) (yellow) and γ i (t 4 ) (purple) when Prior Model 2 was used (recall that t 2 = 0.25, t 3 = 0.5 and t 4 = 0.75). The choices of phase functions γ i in panel (e) match those in Sections 4.1 and 4.2. For the simulated examples, we show the ground truth value of each parameter as a horizontal black line. For all examples, we show the posterior mean of each parameter, estimated using the 100, 000 samples shown in the trace plots, using a horizontal red line. All trace plots suggest convergence to the stationary posterior distribution. Here, we provide empirical evidence behind the claim that the use of an FPCA basis deteriorates estimation performance for µ (see Contribution 3 in Section 1). We demonstrate this using four datasets considered in Section 4: the two simulated datasets from Example 1 in Section 4.1 and the two real datasets from Section 4.2. To estimate the FPCA basis, we follow the following steps: (i) estimate the (centered) sample average defined as μ = arg min µ∈L 2 n i=1 min γ∈Γ ∥µ -(f i • γ i ) √ γi ∥ 2 , (ii) compute the sample covariance K of w i = (f i • γ * i ) γ * i -μ where γ * i = arg min γ∈Γ ∥µ -(f i • γ i ) √ γi ∥ 2 (assuming that each w i is sampled at t 1 = 0, t 2 , . . . , t T -1 , t T ), and (iii) apply singular value decomposition to the covariance matrix K = U SU ⊤ . The columns U provide a data-driven orthonormal FPCA basis for L 2 .The number of FPCA basis functions B f used to specify µ in the proposed Bayesian model is selected based on % variation explained: 90% for the two simulated datasets (B f = 7 and B f = 5 for data simulated under Prior Models 1 and 2 for phase functions, respectively) and for the PQRST complexes (B f = 6); 80% for the Berkeley growth rate functions B f = 8. In all cases, we use B r = 6 B-spline basis functions for the shape-and-size altering random effect. We use the same prior models for phase functions for each dataset as those that were used in Section 4. and j = 2 in (c)). Notably, the FPCA basis appears to capture various sources of variation including noise, which is undesirable if they are to be to model the fixed effect function µ.Figure 9 provides estimation results for (a) µ (ground truth in black when available, centered posterior samples in blue, and centered posterior mean in red), (b) σ 2 (ground truth in black when available, posterior mean in red), and (c) σ 2 c (ground truth in black when available, posterior mean in red). It is clear that the data-driven FPCA basis does not provide a good model for the fixed effect function µ. Only for the PQRST data example, the model yields a reasonable estimate of µ. This in turn results in overestimation of σ 2 and σ 2 c (posterior samples of σ 2 c are very large for the Berkeley data). This suggests that most variation is absorbed into the size-and-shape altering random effect and observation error.We assess sensitivity of posterior inference to under-or over-specification of three hyperparameters in the proposed model: B f (number of basis functions used to model the fixed effect function µ), B r (number of basis functions used to model the size-and-shape altering random effect v i ) and θ γ (concentration hyperparameter in Prior Model 2 (PM 2) on the size-and-shape preserving random effect or phase function γ i ). The data for this experiment is exactly the same as in Example 1 in Section 4.1, i.e., the data was generated from the proposed model with B f = 6 modified Fourier basis functions for µ and B r = 6 B-spline basis functions for each v i . Each panel in Figure 10 shows the centered estimate of the posterior mean for µ (red), associated 95% credible interval (dashed blue), and the ground truth µ (black). Rows 1-3 show results for under-specified, correctly specified and over-specified values of the three hyperparameters, respectively. Overall, posterior inference, in terms of the posterior mean for µ and its uncertainty as ascertained via the 95% credible interval, is very robust to (i) over-specification of B f , (ii) under-or over-specification of B r , and (iii) under-or over-specification of θ γ . On the other hand, when B f is under-specified, we are unable to accurately recover µ as seen in row 1. This is not unexpected since the ground truth µ does not lie in the subspace spanned by the specified basis functions. Thus, in general, we recommend specifying a larger number of basis functions to model the fixed and size-and-shape altering random effects.B f = 4, PM 1 B f = 4, PM 2 B r = 4, PM 1 B r = 4, PM 2 θ γ = 1 B f = 6, PM 1 B f = 6, PM 2 B r = 6, PM 1 B r = 6,We provide estimation results for µ for several additional datasets that have been analyzed in previous literature:1. Pinch force data [Ramsay et al., 1995] , which was also analyzed by Claeskens et al. [2021] .2. Respiration data [Kurtek et al., 2013] .3. Gait data [Kurtek et al., 2013] .4. Signature acceleration data [Kneip and Ramsay, 2008] 5. Gene expression data [Srivastava et al., 2011b] . We omit detailed descriptions of these datasets here and refer the interested reader to the associated references.To specify the model for the fixed effect function µ, we use the modified Fourier basis with B f = 6 basis functions for the pinch force and respiration datasets, and B f = 12 basis functions for the gait, signature acceleration and gene expression datasets. We use B r = 6 B-spline basis functions for each v i in all cases.Figure 11 (a) displays each dataset (top to bottom: pinch force, respiration, gait, signature acceleration, gene expression). Panels (b) and (c) in Figure 11 show estimation results for µ when Prior Models 1 and 2 are used for phase functions, respectively. We show the centered posterior mean in black with a 95% credible interval displayed using blue dashed lines. The credible interval is computed pointwise, using the 2.5% and 97.5% empirical quantiles of the centered posterior samples, (μ j • γ) γ, j = 1, . . . , N . The warpMix estimate is shown in Figure 11(d) . Overall, the proposed Bayesian model is effective in recovering µ. Compared to warpMix, our estimates contain finer geometric features and are more representative of the size-and-shape patterns in the data.• The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results? Answer: [Yes] Justification: All model hyperparameters (for prior distributions) and proposal parameters are specified in Section 3.2 and Appendix D, respectively. Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: Figures 2, 3 and 4 show posterior samples and posterior means to allow for assessment of posterior uncertainty. Figure 11 in Appendix I shows 95% pointwise credible intervals to help with posterior uncertainty assessment.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: The computer resources are specified in detail in the README.txt file in the supplementary material and briefly summarized in Section 5. The execution time is also discussed in Section 5.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)? Answer: [NA] Justification: Our paper does not involve crowdsourcing nor research with human subjects. Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix:  


Decision: The reviewers found that the paper introduces an interesting infinite-dimensional Bayesian mixed model, offering a new perspective on incorporating geometric constraints. They acknowledged the efforts made during the rebuttal, which provided more convincing results. However, they raised concerns about the clarity of the writing and suggested several additions to better motivate the exposition of the method. These revisions should be incorporated into the final version to enhance the quality of the manuscript. I recommend acceptance of this paper.

Reviews: review: paper_summary: paper_summary: In this paper is proposed a mixed model in a functional Hilbert space for a size-and-shape (a type of geometric property) of a square-integrable fixed effect. To this end, the authors consider an isometric action of the infinite-dimensional group of phase functions. Synthetic and real experiments and comparisons with respect to another method in literature are reported, demonstrating how a posterior mean, and under the proposed model, captures the main properties of the function.
strength_weakness: strength_weakness: - I think the notation could be improved, representing, for example, the use of vectors, matrices and scalars differently. In some cases, including dimensionality may help readers. Moreover, some symbols are different along text (see the transpose operator). 
 
- A couple of linear combinations of basis functions are assumed. The rank of every subspace is Bf and Br. This type of formulation has been proposed in previous approaches, following an LDA style. The rank of every subspace is known a prior. 
 
- The authors could include a practical motivation of the study, where they could apply their approach and why it is important.
questions: questions: How is the rank of both subspaces fixed? By hand? How sensitive is the solution with respect to these parameters? Could Bf and Br be inferred automatically? This seems to be a key factor in the model, but the authors do not provide any discussion in this line. 
 
The authors propose to use a B-spline basis, so a pre-defined one. Why this type of function? Note that we can find many piecewise polynomials in literature. 
 
Could both basis functions be learned jointly with weight coefficients? 
 
In real data, n=93, but when the authors use synthetic one, n=100,000. I cannot understand why the difference needs to be so big. In my opinion, if the goal is to sort out the problem in scenarios where the amount of data is limited, this should be also considered in the synthetic scenario. 
 
Figure 1 (b) and (c) are not cited in the introduction. 
 
The method consistently outperforms warpMix on synthetic data. While I like this comparison, I think the authors should consider more datasets for this evaluation. In addition to that, what about computational cost? I would like to see the trade-off accuracy vs. cost for both methods. 

What about failure cases?
limitations: limitations: The authors do not analyze properly the limitations of the paper. In any case, I do not think the work in this paper can produce a negative social impact.
paper_summary: paper_summary: The paper studies the problem of recovering a fixed effect function µ in functional mixed models, where measurement errors and object-level phase variations make the task difficult. It focuses on disentangling the size-and-shape characteristics of µ, which remain invariant under certain transformations. The authors hypothesize that it is feasible to reliably recover the size-and-shape of µ using a Bayesian functional mixed model framework. The formulation is justified by numerical experiments on synthetic and real data.
strength_weakness: strength_weakness: ## Significance:



1. The main proposal of this paper is to integrate phase variations into the model in a phase-preserving manner (equations 3 and 4). Nonetheless, I am struggling to understand what is the significance of such a model. More precisely,
    1. Is this a reasonable model for real applications? Namely, are there some practical applications in which we have reasons to believe the phase functions act on the fixed effect mu in a phase-preserving way?
    2. What are the main conceptual benefits and drawbacks of phase-preserving compared to value-preserving actions? The only answer I seem to find is in lines 142-147, but I do not understand the statement there as I detail in “Clarity”.
2. The comparison of the proposed method with warpMix on real data (Figure 4) seems limited, and this comment relate to the one above on real applications. It is hard to say which method is better when there is no ground-truth - the two methods give different estimates of mu, but it is hard to say which one makes more sense. 

## Clarity:



1. Lines 142-147: Why is “The equivalence class of functions having the same size-and-shape under the norm-preserving action is in an appropriate sense ‘larger’ than the one under the value-preserving action”? 
    1. The explanation following the colon is not sufficient, since the same could be said for value-preserving action: for any g and gamma(t), one can find f(t) := g(gamma^{-1}(t)) such that f(gamma(t)) = g(gamma^{-1}(gamma(t))) = g(t).
    2. The explanation in footnote 1 is a bit hand-wavy and it would be great if rigorous proof could be provided or referred to.

## Novelty:


The idea that the parameterization of the model is not minimal and therefore there is the ambiguity of having an equivalence class of solutions has been there for a long time, as the authors also point out in equation (2). Since (2) is an example in geometric vision which i am a bit more familiar with, I would like to point out more high-level connections here. 


1. A basic idea to resolve the ambiguity is to enforce constraints that favors certain elements in each equivalence class than others. This is more or less what the authors are doing near line 76. Note however that: while two sets of parameters can be equivalent in the sense that if there were no noise then the two sets of parameters explain the observation equally well, in the noisy case the two sets of parameters could handle noise very differently. Therefore, one needs to be smart on which elements in each equivalence class to pick. For example, the work of [A] argues that picking a particular transformation (normalization) in the parameters improves numerical stability of the estimation in the noisy case.
2. Following the same line, two sets of “equivalent” parameters could also be favored differently when the data are corrupted by outliers. An alternative idea to resolve the ambiguity is to estimate the entire equivalent class (rather than finding an element of it), which is for example the path explored in the work of [B].

So far i do not see a strong need for the current paper to exhaust all options, but it would be great if the authors could discuss on how optimal the choice of line 76 is, what might be potential alternatives, etc.

[A] R. I. Hartley, “In defense of the eight-point algorithm,” IEEE Transactions on Pattern Analysis and Machine Intelligence, 1997.

[B] Ding, T., Yang, Y., Zhu, Z., Robinson, D. P., Vidal, R., Kneip, L., & Tsakiris, M. C, “Robust homography estimation via dual principal component pursuit”, In the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2020.
questions: questions: ## Minor comments / questions:



1. Could you comment on how does using MCMC to get parameters from the posterior compare to using optimization, e.g., by writing down the log likelihood and then maximize it over the parameters?
2. The introduction was somewhat confusing, as illustrated below. 
    3. The first paragraph motivates the reader with the Berkeley growth data example, saying that a problem for modeling is that the sample size is smaller than the dimension observation, so one must assume prior / constraints in the hypothesis class. Yet the model in equation (1) is rather general and does not assume much prior for recovery.
    4. Now the second paragraph makes the model even more complicated and harder to estimate by introducing phase variability. In particular, why isn’t model (1) sufficient and why does phase variability make sense? 

    To sum up, I think there are multiple problems motivated and highlighted along the way which somehow buries the key issue of phase variability. 

3. The reviewer is grateful that some basics on phase functions and size-and-shape-preserving transformations are provided in section 2. Nonetheless, when I read Section 2 (before going through the rest of the paper) it is unclear why are these properties important - are they going to contribute to the solution? Which key properties allow you to do something that previously was not done? This for example can be said in the preamble of section 2 to make the flow smoother.
4. Line 119: gamma in Gamma?
limitations: limitations: See above
paper_summary: paper_summary: For the problem of the reliable recovery of a fixed effect function $\mu$, this paper focuses on sampling from and summarizing the posterior distribution of a fixed effect function $\mu$ in a functional mixed model with random object-level phase and amplitude components, without a finite-rank covariance assumption on the error process.

----
Post-rebuttal: The rebuttal addressed most of my concerns and I would like to raise my score.
strength_weakness: strength_weakness: 1)	It seems that we should give different prior models for different functions.
2)	The parameter sensitivity of $\theta$ and the running time of the proposed method should be given.
questions: questions: 1)	Neural networks can also recover functions by sampling. What are the advantages of the proposed method compared to neural networks?
2)	How can we recover arbitrary by a fixed model in Eq.(3)?
limitations: limitations: NA
paper_summary: paper_summary: The paper considers uncertainty quantification for one-dimensional regression tasks, where the observations are noisy. A rather advanced additive model considering invariances under space-time unitary transformations is proposted. Numerical experiments demonstrate the superiority of the approach over other state-of-the-art methods.
strength_weakness: strength_weakness: - It is unclear whether the work will be of interest to the wider machine learning audience and the considered topic and datasets feel more of a niche. 
- Despite the well-written introduction, the paper was quite difficult to understand due to the advanced math and notations -- after reading it I am still quite unsure why I should really care about functional mixed models and their Bayesian inference.
questions: questions: I am slightly puzzled as of why NeurIPS was chosen as a venue for this paper - the work may not get the appreciation here it deserves. I spent quite some time on trying to understand the paper, but its motivations and applications to wider machine learning remain unclear. Perhaps the paper may be better suited towards a statistics journal, such as the ones referenced often in the paper?   If more general motivations are clarified, I am considering to change my score, but accessibility issues are remaining.
limitations: limitations: All limitations are adequately addressed.
paper_summary: paper_summary: The paper deals with estimation functional data, i.e. time-series represented as functions – under a specific observation model, that aims to separately account for a fixed effect (modeled as a mean function) and other variations (modeled as noise functions), with the added confounding variable being a norm-preserving warping function, that leads to observed data.

The paper proceeds to develop an Bayesian approach to estimating the model parameters, under certain assumed distributions for the noise terms, and certain basis functions to represent the functions in a finite dimensional space.

Evaluation is shown both on simulated data and some real functional data examples.
strength_weakness: strength_weakness: Comparisons are limited to a publicly available implementation of the technique known as warpMiX. Default parameters are used as given in the R implementation. At the very least some attempt at matching noise variance assumed in the warpmix model and the proposed model should be attempted. Default noise variance parameter in warpmix is at the level of 10^-3 whereas the noise variances for the proposed model are much higher at sigma_c^2 = 0.25, sigma^2 = 1 etc… Are these comparable in any sense, if not how should this comparison be made more fairly? 

Further, the paper cites a wealth of literature that uses a model similar to the one used, at least as far as the assumed norm preserving action model. There are papers by the cited authors that use these models for estimating mean functions from a stack of observed functions, mainly by Kurtek, Srivastava and colleagues. Many of these references are also cited in the paper, but some attempt should be made to compare to these approaches. 

Beyond simulation, for the task of recovery of functions (and functional parameters) under an assumed observation model, evaluation should include some real downstream task. Otherwise, we are only faced with visual assessments of quality. For instance, the results as shown in figure 4 are claimed to be better in the sense that the estimate is ‘as desired and contains sharp features that are representative of the observed data’ (line 351, pg 9). These are at best subjective assessments, and it is unclear if such an estimate is useful in some quantifiable way. There are many ways to take this to the next step of evaluation, to see if some downstream task such as classification shows improvement after the recovery of the needed parameters.
questions: questions: I would like to know how well the method could be implemented without using basis function approximations – i.e. representing functions simply by their samples. What gets simplified and what gets complicated. What are the practical implications of assuming a specific family of basis functions such as B-Splines. I did notice the data-driven functional PCA basis was tried as well, in the appendix, which seemed to perform worse than B-splines.

Additional questions are on evaluation, which are described in the experiments section of the review.
limitations: limitations: The paper has a reasonable description of the limitations of the model and its assumptions in section 3.

Comments: title: Thank you for your rebuttal
comment: Thank you for considering our rebuttal. If accepted, we will make sure to list future work related to rank/joint estimation as future work.
Thank you for the clarifications!

I am happy with the answers that have been provided. In any case, although recovering the rank is a complex problem, this point should not be hidden and must be included as part of the future work.
Dear Reviewers,

With the discussion period coming to an end tomorrow, we wanted to thank you again for providing constructive comments and considering our rebuttal. We are happy to answer any remaining questions.
Thank you for your additional questions.

**Clarity 1.1, Significance 1.2 / Equivalence Classes**: You are indeed right that $[f]_n$ need not be larger than $[f]_v$ for every function $f$ - that necessarily depends on the function $f$. Our claim on sizes of equivalence classes, however, is with respect to the measure induced on the quotient $\mathbb{L}^2/\sim_n$ compared to the quotient $\mathbb{L}^2/\sim_v$ starting from a non-degenerate measure on $\mathbb{L}^2$; here $\sim_n$ is the equivalence relation on $\mathbb{L}^2$ induced by the norm-preserving action, and similarly $\sim_v$ for the value-preserving action. The example in our response illustrates this when the measure on $\mathbb{L}^2$ is the Wiener measure (with support on continuous functions): its push forward under $\sim_v$ is effectively singular with respect to the Wiener measure, and thus prescribes zero mass to $[f]_v$; in contrast, the push forward of the Wiener measure under the quotient map corresponding to $\sim_n$ is absolutely continuous with respect to the Wiener measure, and prescribes positive mass to $[f]_n$ for every $f$, since the norm-preserving action is an isometry of $\mathbb{L}^2$. Indeed, this does not preclude the possibility of constructing a measure on $\mathbb{L}^2$ that violates the preceding claim; but, in our opinion, such measures are not “natural” since they are not compatible with isometries of $\mathbb{L}^2$, which in this case is the norm-preserving action.

Further, ignoring the random effect, if the fixed effect function is mis-specified with fewer basis functions than required to capture the correct number of extrema, the value-preserving action cannot rectify this, whereas the norm-preserving one can. This has modeling significance since, in a way, the norm-preserving action makes the model more robust. Of course, the size-and-shape altering random effect can introduce new extrema, but only on the individual level and not the population level.

The discussion has brought to light the need to be clearer with the claim on sizes of the equivalence classes in the paper, and we thank you for this.

**Significance 1.1**:

1. Using the Hilbert space $\mathbb L^2([0,1],\mathbb{R})$ as the representation space for the observed functional data is popular, mainly due to availability of an orthonormal basis that enables working with the basis coefficients and facilitates dimension reduction. Given this, modeling the presence of phase variability via the $\mathbb L^2$-norm preserving action of $\gamma$ allows us to use the geometry of the quotient space to develop computational tools; the value-preserving action, in contrast, preserves the uniform $\mathbb L^\infty$-norm $\sup_t|f(t)|$. Moreover, an action preserving any $\mathbb L^p$ norm for $1 \leq p <\infty$ can be converted to an action that preserves the $\mathbb L^2$-norm via the pointwise transformation $f(t)\mapsto \text{sgn}(f(t))|f(t)|^{p/2}$ for every $t \in [0,1]$. The norm-preserving action allows more flexibility in modeling the fixed effect function $\mu$ as compared to the value-preserving action since time warping is accompanied by local rescaling. While a definitive answer to your question is difficult, in many applications including biology, medicine and biometrics, the size-and-shape of $\mu$ is of primary interest for inference. One application that motivates our model, as opposed to one that utilizes a norm-preserving action such as warpMix, is in recovery of the underlying number and magnitude of growth spurts based on a sample of growth rate curves (from the Berkeley growth data). It is well known in this context that children undergo different numbers of growth spurts with different magnitudes occurring at different times. This implies that the data contains phase variation, and the primary task of interest is recovery of the size-and-shape (magnitude and number) of the population level fixed effect function $\mu$. Our estimate, which contains two growth spurts, a smaller initial one and a larger pubertal one later on, agrees with expectations, since it has been shown that children undergo more than one growth spurt during the growth process. The warpMix estimate on the other hand only contains only a single growth spurt.

2. It is commonplace in geometric statistics and statistical shape theory to incorporate the action of a nuisance transformation as an isometric group action (see, e.g., the books Shape and Shape Theory by Kendall et al. and Functional and Shape Data Analysis by Srivastava et al.). For example, when size-and-shape of a landmark configuration is of interest, since rotations and translations do not change the configuration's size-and-shape, they are nuisance, and due to use of the Euclidean norm on the landmarks, act isometrically on the configuration. Model (2) in our paper uses this in its formulation, and inspires our model for functional data in the presence of phase variation.
Thanks to the reviewers for the response! I am clear on all except the following concerns:

**Clarity 1.1**, **Significance 1.2** / **Equivalence Classes**: I am unsure how this intuition helps conclude that $[f]_n$ is larger than $[f]_v$. From your argument, I can agree that functions in $[f]_v$ satisfy the property you mentioned (i.e. same sequence of ordinate values), and that functions in $[f]_n$ do not satisfy that property. However, functions in $[f]_n$ could apriori satisfy some other properties that functions in $[f]_v$ do not satisfy, making $[f]_n$ small. 
- I do not have any clue if that is indeed the case, but your argument does not prevent such cases from happening (if I'm understanding it correctly), which is why I could not understand the intuition. Do I miss something here?

**Significance 1.1**: I guess I can agree on what (i) and (ii) do, i.e., what norm-preserving action and isometric group actions do, but that's not what I was trying to ask. Maybe I am not understanding your response. My original questions are regarding concrete applications, in which we know the better way to model the signals is by norm-preserving actions rather than the alternatives (e.g., value-preserving actions or area-preserving actions)
- Are there such applications? Why are norm-preserving actions a better model (other than empirical numbers showing its better)?
- Maybe you already had them in the paper/response: are you saying that equation (2) is an example? Or time warping?
Thank you again for providing constructive comments during the review period and for considering our rebuttal. If accepted, we will make sure to provide the necessary clarifications in the camera-ready version of the manuscript.
I am satisfied by the rebuttal and have upgraded my decision to Accept. I suggest summarizing the key clarifications in the camera-ready draft if the final decision is indeed an Accept.