Probabilistic size-and-shape functional mixed models
Authors: Fangyi Wang
Karthik Bharath
Oksana Chkrebtii
Sebastian Kurtek
Venue: NeurIPS 2024 poster (2024)
Keywords: statistical shape analysis
size-and-shape perturbation model
Bayesian random effects model
norm-preserving transformation
phase function
Primary area: probabilistic_methods
TL;DR: We develop a Bayesian mixed effects model for inference on the functional size-and-shape mean via a combination of size-and-shape preserving and altering random effects.

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.

Probabilistic size-and-shape functional mixed models
Fangyi Wang Department of Statistics The Ohio State University Columbus, OH, 43210 wang. 15022150 edu
Karthik Bharath
School of Mathematical Sciences University of Nottingham Nottingham, UK, NG7 2RD
Karthik. Berecckk.BBaaaahhhotnnnngge ac .k
Oksana Chkrebtii Department of Statistics
Sebastian Kurtek
Department Statistics
The Ohio State University Columbus, OH, 43210 oksana@stat.. osu. edu
The Ohio State University
Columbus, OH, 43210
kurtek. 1@stat.ttcccmmmm edu
Abstract
The reliable recovery and uncertainty quantification of fixed effect function in a functional mixed model, for modelling population- and object-level variability in noisily observed functional data, is notoriously challenging task: variations along the z and y axes are confounded with additive measurement error, and cannot in general be disentangled. The question then be reliably recovered becomes important. We demonstrate that is possible to recover the size-and-shape of square-integrable y under Bayesian functional mixed model. The size-and-shape of is geometric property invariant to a family of space-time unitary transformations, viewed rotations of the Hilbert space, that jointly transform the æ and axes. A random object-level unitary transformation then captures size-and-shape preserving deviations μ from an individual function, while a random linear term and measurement error capture size-and-shan and-shape altering deviations. The model is regularized by appropriate priors on the unitary transformations, posterior summaries of which may then be suitably interpreted optimal data-driven rotations fixed orthonormal basis for the Hilbert space. Our numerical experiments demonstrate utility of the proposed model, and superiority over the current state-of-the-ase
what properties of μ may
Introduction
Consider a sample of functions from the much-studied Berkeley growth data in Figure T [a), where growth rate curves from measurements heights in centimeters of 54 girls and 39 boys from age I 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
= -
1,...,
comprising three component real-valued functions on [0, 1]: (i) fixed population-level function that represents an average, or representative, change in growth rate; (ii) an individual- or object-level
38th Conference on Neural Information Processing Systems (NeurIPS 2024).

---

(b)
(c)
(d)
30
20
-Original -Value
0.8
15
Norm prov.
0.8
10
10
0.4
0.2
0
0.5
0
0.5
0.5
0.5
Figure 1: (a) Berkeley growth rate curves. (b) Convex phase function 7. (c) One example function f from (a) (blue) transformed by salue-preserving action o y (red) and norm preserving action (f o 7) (yellow). Here, f o has the same classical notion of shape as f. whereas (f o /f has the same size-and-shape as as described in Section 2 (d) PQRST complexes. (c) PORST pattern: P wave (first max), QRS complex (sharp
min-max-min) and T wave (last max) [Pham et al. 2023
random function v; that represents deviations from μı; and, (iil) a zero-mean random measurement enror process 6g. The growth rate function ft(t) of individual i at time differs pointwise from an average p()) via a smooth random translation (1) observed with additive error (1)
A function with any number of peaks representing spurts may be generated under the model, and as such exact recovery of j is not possible without further constraints on v/ and EG. The issue is exacerbated when f; in ( is changed to f10% by incorporating phase variability, or z-axis variability, via a time dilation/contractioo object-level increasing function 7/4 : [0, 1] + [0, 1], which models variability in timings of the spurts. Even for the simpler model, fi % oy + G,, for fixed scalar o > O, recovery of μ is possible only when the measurement error process 6, has a rank one covariance operator [Kurtek al., 201T Chakraborty and Panaretos 2021].
The crux of the problem in reliable recovery of the population-level function ju lies the fact that its topological features such number of critical points and their nature (non-degenerate otherwise), while preserved by "/s may be altered due to addition of e; and 6,. Consequently, a growth spurt (peak) in the average perspective, uncertainty quantification of μ in (I) within a Bayesian framework is unreliable in the presence of phase Y. regardless of how informative the prior distributions on v; and % are. Figure (ngure[c))ce (e) show an additional motivating dataset of PQRST complexes segmented from an electrocardiogram (ECG) signal and typical PQRST pattern, respectively.
time (or (t)) may either be preserved destroyed. From modeling
Summarily, under the functional mixed model ( the topological shape of μ may be irrevocably altered, and thus ordinate values of f 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 μr, its size-and and-sppe,, characterized by joint transformation of its range and domain by the phase function 7,
Contributions. We focus on sampling from and summarizing the posterior distribution of a fixed effect function μ in functional mixed model with random object-level phase and amplitude compo- nents, without a finite-rank covariance assumption on the eтror process. We do not make any claims on efficiency of the computational algorithms, but emphasize that the goal is 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 square- integrable fixed effect by considering an isometric action of the ciinite-dimensionall group of phase functions r' The geometric property of μ thus preserved is referred to its size-and-shapp (Section 2), which we demonstrate may be reliably inferred under the model. To our knowledge, this is the first paper to consider Bayesian functional mixed effects model with 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)
Isometric action of engenders unitary transformation of the Hilbert space; the class of unitary transformations indexed by constitutes rotations of coordinates of the Hilbert space. Upon expressing y and v; a suitable orthonormal basis of the Hilbert space, inference for 7, translates to automatic identification of 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

---

4. We carry out extensive numerical experiments to investigate utility of the proposed model, and demonstrate that the posterior mean of j under our model better captures the properties of μ than the estimate given by current state-of-the-art [Claeskens et al., 202 l (Section 4 and Appendix D
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., 20055, and (ii) pre-specified basis functions [e.g.. Rice
and Wu. 2001 Guo 2002 Chen and Wang 2011 Zhu et al. 2011 |Moris and Carroll 2006 Huo
et al. 2023]. Some recent works have extended methodology to multivariate functional mixed effect 2017], apatial-tempont
models [Volkmann et al. 2023]. scalar on function regression [Liu et al.
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. [021]. which represents the slate-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- pre-prceessing via registration. Of relevance, however, is the work by Schiratti] it al. [2017). where munifold-valuee curves with reparameterization variation are considered.
2 Phase functions and size-and-shape-preseree transformations
Without loss of generality, we assume that all functions are observed on a fixed domain [0, 1]. Let (10, 1]. R) = {f [o, 1] + R fo dt) < <o)} (henceforth simply referred to as IL3) denote the representation space of interest, i.e., the space of real-valued gquare-integrable functions on [0, 1]. The space I.2 when endowed with the norm || || == 16 [ coming from the inner product (f, = fo (t) is Hilbert space. Additionally, let I {r [0, 1] - [0, (/))(0) 0, = 1, > 0} denote the group of orientationn-peesrrinng diffeoomopphsms (5: is the time derivative of 7) of [0, 1]. In the functional data analysis literature, is used to model phase variability in functional observations. Thus, elements of I will henceforth be referred to as phase functions. Note that is Lie group with composition (71, + 7/1 o 722 as the group operation, where o is function composition. This implies that (i) r is an finitee-iimnsiooall smooth manifold, (ii) (71 o 72) o 733 = 7/1 o (7/ o 7/)) for any 71, 722 73 E r, (iii) it contains the identity element Tha(t) such that y o Yid y for any 7 E r, and (iv) for any 7 I there exists y r such that 7 y = "/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 I can act on the function space IL,2 from the right, engendering maps IL2 x r + IL2, in different ways, thus resulting in different notions of phase variation in functional data.
1. Value-preservingg action. The value-preserviv mapping is defined via composition: foy, f E L2, 7 E T. This action is referred to the time warping of function since and I y traverse the same exact y axis values, but at different times (z 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.
Area-preserving action. The pree-ries mapping is defined as (Jo I L2, 7 P. This action is commonly used for statistical analysis of probability density functions.
3. Norm-preserwwe
mapping is defined as (f.>) := (f o 7))) f L2 E I. An important property of this action is that it preserves the L norm of a func- tion: |l = 7)ll The action has been profitably used in the problem of function align- ment/registration, wherein desideratum is cost function invariant to simultaneous action of the group r, or time warping [Srivastava and Klassen 2016 Chapter 4].
action. The norm-preservin
The operator D. : 2 + L3, f D+( = (f o v arising from the norm norm-preserving action for a fixed 7 I' surjective isometry, and hence unitary. Thus, D,, for fixed 7 E I, is an dimensioai rotation in the Hilbert space L?, and the group D == {Dz, 7 I') of rotations
of L plays prominent role in classical white noise calculus Hida 2015
3

---

A general space-time diffeomorphism (t,x) + graph {(t. f())] of f viewed as subset of [0, 1] x R. The operator D.y is associated with a special space-time diffeomorphism with 01 (f, r) = (t) and 0z(t, x) = x /5((). In this sense, the size and shape of f. as it relates to its norm and its graph, is preserved by operators in D.
x)) preserves the topology of the
To better understand how D., transforms f. consider Figure I [b)((c). As seen with the red curve in (c), the value- preserving action of y preserves the image of / + [(t), where ordinate values of 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 the level sets of of 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 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 (f) for any y E Г are equivalent in their size-an aare-and-shee and-shape.
action in that the topology of
Use of the operator D, is particularly useful for modeling phase variation in model
(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 aalueppreerrvig warping and an associated local scaling of function values, i.e., when +() > 1 (5(1) < 1) the function value f (t)) is warped to the right (left) relative to ft) and additionally rescaled by factor of v+(t).
(iii) The equivalence class of functions having the same side-and-shape under the norn-preserving action is 'larger' than the one under the malue-preserving action in the following sense. The map क: T.2 - / ~ takes a function to its equivalence class determined by equivalence relation ~ can be shown that the measure of equivalence class, under the pushforward of non-degenerate measure with support in L under n, is larger when ~ pertains to norm-present 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- porm-prregg action helps with reliable recovery of the size-and-shape of the fixed effect μ4, in contrast to the more traditional notion of shape of μ under the value- pelleepreserving action, wherein the additive linear term U; + c, more easily alters the shape of j..
3 Bayesian size-and-shape functional mixed model
The operator D., related to the norm- norm-preserving action of I' is an isometry of L22 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 config configurations [e.g.. Dryden and Mardia 2016). To elaborate, let X; E R.K xx2 denote set of K landmarks, i.e., X; is matrix that contains the coordinates of K points in 2 and S0(2) = {R E R2221 with matrix multiplication the group operation. Consider the following perturbation (generative) model for the size-and-shape of an object represented via a landmark configuration X;:
" R = RR" = I, det( = }} denote the rotation group
X; (M + E,)R, + 1
where t; E ? is a random translation, 1 E K is a vector of 1s, R, E SO(() is a random rotation, M E Kx2 is fixed template, and E, Rx2 is the error. Using this model, it is of primary M based on observed landmark configurations
interest to estimate the fixed effect ssize-and-shape Xı,..., Xn. Transformation of the model components M and E, using rotations RG (and translations) in order to align them to the observation X; may be viewed as an isometric, or norm-preserving. transformation of the coordinate system for (M + F.). so that only the size-and-shape of M can be reliably recovered, but not M Lele and McCulloch 2002).
Inspired by size-and-shapp shape analysis of landmark data, and in contrast to existing works using the value- preserving action |Raket et al. 2014 Claeskens et al. 202T], our proposed functional mixed effects model utilizes the norm- Norm-preserving action of I on L Analogous to the model in
( for a function fi E IL, one can define a functional perturbation model fa = D (μ + 4) = [(4 + G) %lVi% for a template function a with random phase functions 7s wherein the II2 coordinate system of the model components μ and e; is aligned to the coordinate system of the function f. via rotation using 7/

---

To the fixed effect, or signal plus noise model, one can additionally introduce an object-level random effect to increase modeling flexibility, resulting in:
Se Dry (u + v + e.) Uu+ i4 + 70l+%s
Thus, "/ and v 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 fi. and it is possible to reliably infer only the size-and-shape of μ.
While the model in (3) is written in terms of infinite-din dimeesioaia objects, functional data is often observed under a common discretization at time points t = 0 < tz < < tr- < Br 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 the phase component, and demonstrate the utility of the size-and-shape perspective in inferring geometric properties of j..
Let f: = (fi(t))..... ..... ftt.......tt..)......... )(( E R i= = 1...., represent the discretized function values for each observation. The discretized observation model is then given by
fc(t,) 1/0 + i + 4:) rel(tt)
-1 ....., n, j ,,..., T.
Letting (ox) and { k== denote two orthonormal basis systems for L22 we represent the model components μ and 2 as linear combinations of basis functions. For dimension reduction, we define J : >== 0x and v;- := Ex= ŌR i.e., the two basis sets are truncated to By and B, 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) N2,6~~ o?). The discretized error is assumed to be iid ((t))~~((0,'))
3.1 Choice of basis functions and likelihood
Specifications of and U; require appropriate orthonormal basis functions, and we consider modified Fourier basis or the B-spline basis for μu, and the B-spline basis for us. The modified Fourier basis contains the following elements { (V3t, v3(1- (5(1-)). t), v2 cos((2jjt) v2 sin(nnt): j = 1, 2, [0, 1] and is subsequently orthonormalized via the Gram-Schmidt procedure under the T.2 metrie: larger values of yield basis functions with finer harmonies over the domain [0. 1]. This basis is used to specify the fixed effect function μa 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 L2 metric. In some cases where the underlying f1 is expected to exhibit more local features, modified Fourier basis. Finally, one needs to choose the number of basis functions to model μ and the application of interest. For example, one may use larger number
: :.;
use the B B-spline basis rather than the
2g- This choice depends of basis functions for μ when the observations have a complex shared global structure. Effects of misspecifying the number of basis functions for μ and i 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 effect coefficients o². and (iii) error variance o². Thus, to simplify the inference, we marginalize the likelihood with respect to the size-and- -nd--aaddpppe altering random effect. Let 4 E R denote the matrix of evaluations of the fixed effect basis, whose (j. k)th entry is given by (dx o (tts)) vi.(t,)) One can view D; as a discretization of the rotated coordinate system, via the norm-preserving action that uses y,, for μ.. Similarly, let Φ; E TR'',, Br be the matrix of evaluations of the random effect basis, whose (j, k) entry is (ōx (t)) (t), Let a EB, and c, E RBB, i i= -I..... be the vectors of basis coefficients for μ and U2 i I,..., n. respectively. Finally, let c,"h E ERT, i= 1....... be the vectors of ry-transformed observation errors, with entries (E, (t,) V j= =, L,... ... T that are independently distributed as N NN, (t)), conditional on 7: the additional time-dependent scaling of the eтror variance comes from the norm- pereeereeniing action. Using this simplified notation, we can rewrite the model in in matrix form as
(......, aB,) (ii) variance of the size-and-shape altering random
f. =.a + $.C, + c7 i i= -1..., .,
(5)
5

---

which may be usefully compared to the size-and- shape perturbation model (2 for landmark con- figurations. With MVN used denote the multivariate normal distribution, the distribution of f.. conditional on the vector of coefficients C, is fi C; MVN a + $,,,, o² dfdiag(*sstt))) where diag(fr th (fa(t))) is T x T diagonal matrix whose jth diagonal element is (t(z)) Then, the follow- ing marginal distribution of f: is used to define the likelihood function (see Appendix C for full derívation):
fi ~ MVN(D, o² diag(")(t)))) +
+ +00
(6)
3.2 Prior distributions
The model parameters that need to be estimated are the fixed effect coefficients a, random effect variance o² variance of observation eгтоr o², and individual phase functions ". For a, o² and a² we use weakly informative prior distributions: a ~ MV, 10000I of [G 0.01), 2² ~ IG(0.01, 0.01) (IG is the inverse-gande gamma distribution).
The prior distribution over the space of phase functions I, to model the shape-and-size preserving random effect, can be specified in different ways [e.g., Telesca and Inoue 2008 [Bioot] 2013 Lu et al. 2017]. In this work, we use two tractable prior distribution models I' to guard against confounding of inference between {7x} and μ:: (i) one-parameter family; (ii) nonparametric finiteddimensional
family compatible with the time discretization. Importantly, the prior models rely on the group structure of r.
Prior Model (PM1) r. Each phase function 7/j is defined via single parameter o; (Section 5.2 in Srivastava and Klassen (0166) as
(t) t att(t - 1), a; E (-1,1), E [0, 1].
The phase functions defined in 7 form a one-dimensional subset of r, making posterior inference more tractable with significant reductions in computational complexity. This subset contains phase functions with (tt E (0, 2) for all & E (0, 1) and no inflection points. At the same time, we sacrifice flexibility terms of the allowed rotations of 2 The ai, , n are assumed be independent a priori following the Uniform(- -1, 1) distribution.
Prior Model 2 (PM2) on r. A more flexible nonparametric prior model utilizes a point process- based prior distribution related to the Dirichlet process on I (Bharath and Kurtek 2020). Under discretized time, each phase function may be represented by a finite number of successive in- crements p(7)) (tx( - (0)..... .tt,., - (tt--)),,... (1) ((tr,-1)) E RT,-1; as
Tyy + oo, the resulting prior has dense support r. Then, the finite-dimensional prior distribution is placed on the vector of phase increments,
p(7)) ind (fr,
(8)
where t = (12, t/- t,,..., ftr,,, -1) is defined via T, consecutive time points on [0, 1] and 0, 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, usually choose T, to be relatively small, e.g., five or seven, to simplify the prior model. Second, large value for 0., regularizes the phase functions toward "id- We use 02- = 30 in all of our numerical experiments. The resulting % is a piecewise linear function with changes in the slope 7/s at f2..... .2....,,..11. This prior distribution is more flexible than the one PMI. 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- Netropolis-Hastiins 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. Trace plots provided in Appendix [ indicate good convergence for all of the examples presented in this manuscript. Full details of the proposal distributions and the MCMC algorithm are in Appendices D and

---

(b)
(d)
3000
0.8
0
2000
2000
0.8
-10
Centered Post. Samp. Centered Post. Mean warpMix Estimate Truth
1000
0.4
10
0.2
-Post. Samp. -Post. Mean
0.08
0.1
0.12
Truth
0.5
0.5
Post. mean: 0.10
Post. mean: 0.28
0.5
3000
8000
5
0.8
2000
0
-6
Centered Post. Samp. Centered Post. Mean
1000
1000
0.4
warphiix Estimate Truth
0.2
Post. Samp Moan -Truth
-15
15
0.1
0.14
0.2
0.3
0.4
0
mean:
Figure 2: Row 1: Phase functions from PM1. Row 2: Phase functions from PM2. (a) Simulated data (n = 30). (b) Estimation of μr: ground truth (black), posterior samples (blue), posterior mean (red), warpMix esimate (yellow). (c) Histograms of posterior samples for o² and o². respectively (posterior mean in red; ground truth in black). (e) Estimation of phase function a randomly chosen observation: ground truth (black). posterior samples (blue), posterior mean (red).
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-ard frequentist functional mixed model [Claeskens et al. 202T]; we 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 al., 2016), which models functions under the popular square-root velocity transformation [Srivastava et al. 201 as realizations of Gaussian process centered at mean function μ.. Phase variation is incorporated via the gelleepressriing action of I with prior model that is the same as PM2 in our model.
The main object of interest for posterior inference is μ. Note however that, discussed earlier, i and (+( (μ 7 for any 7 E T are equivalent in terms of their size-and-shape. Thus, we first center all posterior samples for μl, via 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 the orbit centering in Srivastava al. (0111]. Assuming that we have N posterior samples of each 7/ i I,..., n,, we first compute = 1/=n)(n) E Σi= 1/. where /f is the jth posterior sample of the phase function 7s: due to convexity of I, note that 7 r. We then compute the centered posterior samples of μ,, (ft 7) V7. which may be visualized directly or used estimate posterior summaries, c.g., pointwise posterior mean and credible interval. In some cases, we also visualize the posterior samples and pointwise posterior mean for phase function Y: 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 y and %: we use uniform subsample of size 000 to ensure that all plots are easily readable.
4.1 Simulations
Example 1: data generated from our model. We first consider an example based on data simulated from model in (6). We use B, = 6 modified Fourier basis functions for μ and B, = 6 B-spline basis functions for each v;. The ground truth variances are o² 0.25 and o² 0.1. Then, to generate the data, we sample a ~ MVN(0, JB,), ), and consider two cases for phase functions based on PMI and zid PM2: (i) a; ~ (niform(-1. 1), and (2) p(7)) itd (30t t = (0, 0,,5. 0.5, 0.75, 1). The sample size in this simulated example is = 30.
F displays results based on data generated via PM1 for phase functions (row 1) and PM2 (row 2). The simulated data is shown in panel (a). Upon visual inspection, it is difficult to discern the underlying μ.. Panel (b) displays estimation results for pez ground truth in black, centered posterior samples in blue, centered posterior mean in red, and warpMix estimate in yellow. The proposed
7

---

(b)
1
0.5
0
Centered Centered Post. Mean warpMiix Eatimate -Truth
-Cantered Post. Semp.
Post. Samp.
0.5
-Centered Post. Mean
2
Centered Post. Mean
-warpMix Estimate
warpMixx Estimate
20 Figure 3: Comparison of estimation results based on Model 2-B and warpMix for (a) μr. (b) μlz and (c) pa- In each panel, we show the ground truth (black), centered posterior samples (bluc), centered posterior mean (red), and warpMix estimate (yellow).
-Truth
3
-Truth
0.5
0
0.5
0.5
Table I: Comparison of fixed effect estimation accuracy based on posterior mean from proposed Bayesian models and warpMix estimate. Smallest estimation errors are highlighted in bold.
warpMix
Model I-F
Model 1-B
Model 2-F
Model 2-B
p1
0.0179
0.0194
0.0151
0.0193
0.0182
/42
0.0394
0.0134
0.0235
0.0070
0.0152
0.0077
0.0195
0.0111
0.0044
0.0033
model is effective at recovering the underlying size-and-shape of jr, which contains two local minima and maxima. On the other hand, the warpMix estimate only contains 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 g² and o 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 phase function 7s corresponding to a randomly chosen observation. As before, posterior samples are shown in blue, posterior mean in red and the ground truth in black. both cases (phase function generated from PMI or PM2), we reliably recover the underlying ground truth.
Example 2: data generated from warpMix model. Next, we consider more challenging scenario for the proposed model and specifically focus on a comparison to warpMix. In this case, phase functions and random effect functions generated using the warpMix model; set o² 0.25 and o² = 0.0001. Comparison of estimation accuracy is based on three different fixed effect functions: Nn (t) = (3rr) + 3rr//44 /4, p2(() = ex--(tt-.233) /0.04 + exp (1-0.755) /0.02 and (ts(t) = cos (xx
c2st + T7 1/2),, (0, 1]. To generate the data, we use the preueerpeeeviigggg action as in the warpMix specification, and not the norm-presse com--preeervng action that is utilized in our model.
To specify our models, we use By 6 modified Fourier basis functions or B-spline basis functions for μ with PM1 or PM2 for phase functions with T. 7, and (ii) B, 6 B-spline basis functions for each v/- In total, 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 present estimation results for (a) μ1. (b) #2 and (c) μı based on Model 2-B. The ground truth is 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 - -ty). The [ 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 μ사 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 warpMix For /43, 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
j=1 evaluation criterion used in Claeskens et al. [021]]: A := ET= [A(t,) =(t,
f for our models is defined as the centered posterior mean. Table
8

---

(b)
(d)
20
Centered Post. Somp.
4000
4000
20
Rotaled fixed effect -Functionl evservation
-Centered Meen
-war wwrrphix Entimate
3000
3000
0.8
15
10
2000
2000
10
5
0.4
1000
1000
0
-Post. Samp.
-5
1.2
1.3
1.4
6.2
0.3
0.4
0.5
Moan
Post. Mean: 1.27
Post. Mean: 0.34
Centered Post. Samp. -Centered Poet. Mean
3000
3000
10
-Rttttddd fixed effect Functioral ovservation
0.8
2000
2000
1000
1000
0.4
Samp. -Post Mean
-2
0.18
0.2
0.22
0.24
8.01
0.015
0.02 0.025
0
22
Post. Mear: 0.21
Post. Msan: 0.02
0.5
Figure 4: Estimation results for Berkeley data (row 1) and PQRST complexes (row 2). (a) Posterior samples
(blue) and posterior mean (red) of μ. and warpMix estimate (yellow). The warpMix model was unable to yield an estimate of j for PQRST data. (b)&(c) Histograms of posterior samples for a² and o². respectively (posterior mean red). (d) Posterior samples (blue) and mean (red) of phase function for randomly chosen observation. (e) Observation corresponding (d) (black) with rotated posterior samples of ju (bluc).
4.2 Real data examples
We now consider application of the proposed modeling framework to (i) Berkeley growth rate functions (n = 93) (Srivastava et al. 20116) (Figure [a)), and (ii) PQRST complexes (n = 40)
(Kurtek et al. 2013) (Figure dd)). 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 patter,, 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 she-and-shape of p.
For Berkeley data, we use modified Fourier basis for μ with B 6, B-spline basis for each v; with B, 6, and PM1 on I; for PQRST data use B-spline basis for ju with B, 12, B-spline basis for each i; with B, = 6 to better model the sharp local features of the QRS complex (Figure ( T))), and the PM2 model on I 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. [0111]. The marginal posterior uncertainty for a 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 a² and o². 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 xg (f²), the posterior samples of ju rotated using corresponding posterior samples of the phase function, where i, index the observation and posterior sample, respectively. Note that the blue functions fit the observed black function well.
Row 2 in Figure 4 shows shows estimation results for the PQRST data. Panels (a)-(c) are the same as in the previous description. The estimated size-and-shape of μ resembles 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 r. Finally, panel (e) shows that the rotated posterior samples of i fit randomly chosen observation well. Appendix [ shows posterior mean and 95% credible interval estimates, as well as warpMix estimates, of μ for five additional datasets.

---

Discussion
Numerical experiments in Section 4 and Appendices A and demonstraat benefits of the proposed
mixed model in recovering the size-and-shape of 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 approximately 93 and 111 minutes, respectively, under PMs and 2 on I, based on = 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 U; a priori. Alternatively,
could treat the number of basis functions By and B, as random and estimate them. This, however, would require more advanced MCMC algorithms for posterior inference. Further,
the modified Fourier basis and B 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 ju 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 p. which we plan 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 sparselyfiregularly sampled and fragmented, or partially observed, functional data [Matuk et al.
2022] The norm-prees norm-preserving action of I may be used perform inference for the size- shapp-addeeeeeeeeeee of
fixed effect parameterized open/closed curve [0. 1] Rd
Srivastava et al. 2011a Kurtek
al. 2012 Finally, for two-dimensional parameterized surfaces J D: R² IR similar
norn--reservvin action of the reparameterization diffeomorphism group with elements 7 D
Jermyn et al. 2017] may be used to infer the mean size-and shape of fixed effect surface.
Functional data is arising as 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 [)]. The ideas we have explored the proposed framework can be applied more broadly 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.
Acknowledgments and Disclosure of Funding
This research was partially funded by NIH R37-CA211955 (to SK and KB), NSF DMS-2015377 EPSRC EP/V048104/ (to KB), and NSF CCF-1740761 and NSF DMS-2015226 (to SK). The authors have no competing interests to disclose.
10

---

References
K. Bharath and S. Kurtek. Distribution on warp maps for alignment of open and closed curves. Journal of the American Statistical Association, 115(531):1378 1392, 2020.
J. Bigot. Fréchet means of curves for signal averaging and application to ECG data analysis. Annals of Applied Statistics, 7(4):2384-2401, 2013.
V.I. Bogachev. Gaussian Measures. Number 62. American Mathematical Society, 1998.
M.M. Bronstein, J. Bruna, Cohen, and Veličkovié. Geometric Deep Learning: Grids, Groups, Graphs, Geodesics, and Gauges. arXiv, 2104. 13478, 2021.
A. Chakraborty and V.M. Panaretos. Functional registration and local variations: Identifiability, rank, and tuning. Bernoulli, 27:1 103 1130, 2021.
H. Chen and Y. Wang. A penalized spline approach to functional mixed effects model analysis. Biometrics, 67(3):611-870, 2011.
W. Cheng, I.L. Dryden, and X. Huang. Bayesian registration of functions and curves. Bayesian Analysis, 1122)447-475, 2016.
G. Claeskens, E. Devijver, and I. Gijbels. Nonlinear mixed effects modeling and warping for functional data using B-splines. Electronic Journal of Statistics, 15(2):5245 5282, 2021.
I.L. Dryden and K. Mardia. Statistical Shape Analysis: With Applications in R. Wiley, 2016.
W. Guo. Functional mixed effects models. Biometrics, (8(1))121-12)) 2002.
T. Hida. Stationary Stochastic Processes (MN-8), volume 8. Princeton University Press, 2015,
Huo, J.S. Morris, and H. Zhu. Ultra-fast approximate inference using variational functional mixed models. Journal of Computational and Graphical Statistics, 32(2):353-365 2023.
I.H. Jermyn, S. Kurtek, H. Laga, and A. Srivastava. Elastic Shape Analysis of Three-Dimensiooal Objects. Synthesis Lectures on Computer Vision. Springer Cham, 2017.
A. Kneip and J.O. Ramsay. Combining registration and fitting for functional models. Journal of the American Statistical Association, 103(483): 2008
Kurtek, A. Srivastava, and W. Wu. Signal estimation under random warp-ngs and nonlinear signal alignment. In Advances in Neural Information Processing Systems, volume 24, 2011.
S. Kurtek, A. Srivastava, E.P. Klassen, and Z. Ding. Statistical modeling of curves using shapes and related features. Journal of the American Statistical Association, 107(499): 152-1165, 2012.
S. Kurtek, W. Wu, G. Christensen, and A. Srivastava. Segmentation, alignment and statistical analysis of biosignals with application to disease classification. Journal of Applied Statistics, 40:1270-1222 2013.
S.R. Lele and C.E. McCulloch. Invariance, identifiability, and morphometrics. Journal of the American Statistical Association, 97(459):796-80)) 806, 2002.
B. Liu, L. Wang, and J. Cao. Estimating functional linear mixed-effects regression models. Computa- tional Statisties and Data Analysis, 106:153-166 164, 2017.
Y. Lu, R. Herbei, and S. Kurtek. Bayesian registration of functions with a Gaussian process prior. Journal of Computational and Graphical Statistics, 26(4))894-904, 2017.
J. Matuk, K. Bharath, O. Chkrebtii, and S. Kurtek. Bayesian framework for simultaneous registration and estimation of noisy, sparse, and fragmented functional data. Journal of the American Statistical 2022.
Association, 117(540):1964-19))
J.S. Morris and R.J. Carroll. Wavelet-based functional mixed models. Journal of the Royal Statistical Society: Series B, 68(2):179-199, 2006.

---

B.-T. Pham. P.T. Le, T.-.. Tai, Y.-. Hsu, Y.-H. Li, and J.-C. Wang. Electroeardiogram heartbeat classification for arrhythmias and myocardial infarction. Sensors, 23(6):2993, 2023.
L.L. Raket, S. Sommer, and B. Markussen. A nonlinear mixed-effects model for simultaneous smoothing and registration of functional data. Puttern Recognition Letters, 38:1-7, 2014.
J.O. Ramsay, X. Wang. and R. Flanagan. A functional data analysis of the pinch force of human fingers. Journal of the Royal Statistical Society: Series C, 44(1):17-30, 1995.
A.R. Rao and M. Reimherr. Modern non-linear function-on-ffunction regression. Statistics and Computing, 33(6))130, 2023a.
A.R. Rao and M. Reimherr. Nonlinear functional modeling using neural networks. Journal of Computational and Graphical Statistics, (3(4):1248-12)) 2023b.
J.A. Rice and C.O. Wu. Nonparametric mixed effects models for unequally sampled noisy curves. Biometrics, 57(1):253-259, 2001.
J.-B. Schiratti, S. Allassonnière, O. Colliot, and S. Durrleman. A Bayesian mixed-effects model to learn trajectories of changes from repeated manifold-valie observations. Journal of Machine Learning Research, 18(133):1-33 2017.
G. Simehoni and S. Rosset. Integrating random effects in deep neural networks. Journal of Machine Learning Research, 24(156):1-57, 2022.
A. Srivastava and E.P. Klassen. Functional and Shape Data Analysis, volume Springer, 2016.
A. Srivastava, E.P. Klassen, S.H. Joshi, and L.H. Jermyn. Shape analysis of elastic curves Euclidean spaces. IEEE Trans. Pattern Analysis and Machine Intelligence, (3(7):1415-144) 2011a.
A. Srivastava, W. Wu, S. Kurtek, E.P. Klassen, and J.S. Marron. Registration of functional data using Fisher-Rao metric. arXiv, 110..3817v2, 201
M. St. Ville, A.W. Bergen, J.W. Baurley, J.D. Bible, and C.S. McMahan. Assessing opioid use disorder treatments in trials subject to non-adherence via functional generalized linear mixed- effects model. International Journal of Environmental Research and Public Health, 19(9):5456, 2022.
D. Telesca and L.T.T. Inoue. Bayesian hierarchical curve registration. Journal of the American Statistical Association, 103(481):328-44 2008.
A. Volkmann, A. Stöcker, F. Scheipl, and S. Greven. Multivariate functional additive mixed models. Statistical Modelling. 23(4):303-32 2023.
Yao, H.-G. Müller, and J.-L. Wang. Functional data analysis for sparse longitudinal data. Journal of the American Statistical Association, (00(470):577559, 2005.
H. Zhu, P.J. Brown, and J.S. Morris. Robust, adaptive functional regression in functional mixed model framework. Journal of the American Statistical Association, 106(495): 2011.
H. Zhu, K. Chen, X. Luo, Y. Yuan, and J.-L. Wang. FMEM: Functional mixed effects models for longitudinal functional responses. Statistica Sinica, (4)):2077 2019.
12

---

We include the following appendices as support for the content in the main paper:
A. Additional comparison of estimation accuracy with oftt--o-thhearr methods.
B. Additional motivation behind the norm-preserving action of T on L22
C. Detailed derivation of marginal likelihood in
D. Specifications all proposal distributions and derivation of the Metropolis-Hastine ratio for all model parameters.
E. Detailed Markov chain Monte Carlo algorithm sample from the posterior distribution.
F. MCMC diagnostic plots for Example 1 in Section 4.1 and the two real data examples in Section 4.2
G. Results of using empirical FPCA basis to model the fixed effect function
H. Results of sensitivity analyses assess effects of hyperparameter misspecification.
Estimation results for five additional real data examples.
A Additional comparison of estimation accuracy
Table reports a comprehensive quantitative evaluation, in terms of estimation accuracy for the fixed effect function μr, on the five simulated datasets used in Section 4. Rows 1-2 consider data simulated from our model under Prior Models (PM 1) and (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 (4) to those produced using warpMix Claeskens al. 1 and the Bayesian model proposed in Cheng et al. [2016) (BRFC). As seen in the table, our model outperforms warpMix and BRFC all of these simulation scenarios.
Table 2: Comparison of fixed effect estimation accuracy based on posterior mean from proposed Bayesian models, and BRFC (posterior mean) and warpMix estimates. Model 1-F to Model 2-B correspond to models where the number indicates the prior model phase and the letter corresponds the type of basis used model the effect function (F=Fourier, B=B-splines). Smallest errors are highlighted bold.
warpMix
BRFC
Model 1-F
Model 1-B
Model 2-F
Model 2-B
PMI-F
0.7972
2.5738
0.0452
0.2542
0.0746
0.4039
PM2-F
0.6417
2.1379
0.1152
0.2457
0.0539
0.2878
warpMix-Ju
0.0179
0.1627
0.0194
0.0151
0.0193
0.0182
warpMix- Mi
0.0394
0.0103
0.0134
0.0235
0,0070
0.0152
warpMix-44s
0.0077
0.0582
0.0195
0.0111
0,0044
0.0033
B Additional motivation behind norm-preserving action
Figure provvdds an additional example of the difference between the value-preserving and norm- preserving actions. Panel (a) shows two phase functions generated from Prior Model 1. Panel (b) shows the first six modified Fourier basis functions. In (c)&(d), we show the same basis functions after applying the value- gelleepreserving and norm- pre--preeriing actions using the phase functions in (a). Finally, in (e). we provide an example function formed using 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 well as the transformed red function contain four extrema both rows. On the other hand, the yellow functions have five extrema (there is small local maximum in the yellow function in the bottom row near 0.1).
13

---

(b)
(c)
(d)
0.5
-5
Origgnal VValue prav.
-2
-2
Norm prev.
0.5
0.5
0.5
0.5
0.5
0.5
-10
-Oniginal
-2
VVllue prov.
Norm prsv.
0
20
0.5
0.5
0.5
0.5
0.5
Figure 5: (a) Phase function. (b) Six modified Fourier basis functions. (c) The same basis functions as in (b) after value preserving action using phase function in (a). (d) Same as (c), but using norm- preserving action. (e) Function formed using the same linear combination of basis functions in (b) (blue), (c) (red) and (d) (yellow).
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- comm-preserving action. We vary the number of basis functions from 30 and average the residuals over all n = 93 Berkeley growth rate functions. Given a function f E L2, its projection onto the basis is =ド= م where å, J/' (t))) used in the projection, and oj. j I...., B are the basis functions. The residual for (i) is defined as |l. fll. while the residual for (ii) is min,er Il. (f //l. Overall, we expect the average residual from (ii) to be lower than the average residual from (i). However, interesting aspect of this experiment is to understand how much more expressive the basis when an additional rotation via the form-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 -,....,37 30, based on (i) and (ii) are shown in blue and red. As expected, (ii) yields lower average residuals, with very large gap for 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 small number of basis elements is used in the specification of in our model, the norm- weee-preeriing action allows for more flexible and efficient modeling, by rotating the coordinate system for f the observed data.
1.... B. B is the number of basis functions
1,.... 10. The inclusion of the norm-preen corm-preserving action becomes less valuable
Proj. Proj. & Norm prsv.
Avg. Residual $2
10
20
30
Num. of Basis
Figure 6: Average residual of (i) projection onto modified Fourier basis functions (blue), and (ii) projection followed by optimization over I under norm-preserving action (red).
14

---

C Derivation of marginal likelihood
Conditional on the coefficients of the size-and-shape altering random effect, fo follows a multivariate normal distribution:
fice ~ MVN( + 0, 0² (t))))
The expected value and variance of fi are
BjA] = E| + = + 2.a + [. D.a.
and
Var[fi] [E[aar f[]]]]] + Var[B[[f]ca]] o diagg(t:())) Varl 0.a)
didig (1,)) + [a] o'diag(fi (f:(4)) + .....
Thus, the marginal distribution of fi is
f. MVN o'a o
(9)
D Proposals and derivation of Metropolis- acceptance ratio
D.1 Proposal distributions
We use the following proposal distributions in the Markov chain Monte Carlo algorithm:
Fixed effect coefficients: aa (a Ea).
Variance of error process: (o²)can ~ TN((o?) 72, o0).
3. Variance of size-and-shape altering random effect: (o²) ~ TN((o2) - 0, 00).
4. Size-and-shape preserving random effect (phase functions) under Prior Model 1: agan Uniform(( or."r + 8).
Size-and-shape preserving random effect (phase functions) under Prior Model 2: Yean Your 7, where p(7) Dirichlet(at).
TN stands for the truncated normal distribution. Notably, the proposal in 5. for phase functions utilizes the group structure of T. The proposal covariance matrix Σn 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, of and v2 respectively, well the parameter and precision parameter for the size-and-shee preserving random effect under Prior Models and 2, respectively, are adapted according to the acceptance rate during the burn-in period.
D.2 Derivation of Metropolis- Metropolis-Hastings Hastings acceptance ratio
Let L(-I-) denote the likelihood function, and π (-I-), p(-) and g(-l.) denote the posterior, prior and proposal distributions. Then, the general form of the MH acceptance ratio for parameter 8k, an element of parameter vector B, with observed data denoted by f.
π (B 3 (f) ( (
}
Ps. = min
x (S (-A) If) q ggam |
Lf k)) @ran ) (66urr 18
min
(10)
Lf geur 3eLh) )p())(((3 ( (Bean [3gur)
where B(-k) is the vector of parameters excluding 3r. Note that, for the model proposed in this manuscript B (a, 03, 7)).
15

---

Fixed effect coefficients a. Let B-a denote the vector of parameters excluding a. Since the proposal a random walk centered at the current value, it is symmetric. Thus,
Pa = min {1, L(fla ,350a can
{ L(flaeur 3Ba0a (ord)p((es ((a
(11)
Variance of error process, o², and variance of shape-and-size altering random effect, o² o² and o have the same prior and proposal distributions, we derive the acceptance ratio pa* only, and note that Pea is the same. The proposal for these two variance parameters is a truncated normal distribution. Thus,
Since
I o (o")aur Pycan (((o2)our) (o2))an)
+ (02)
Φ5 (o2)can -(a Ryruuu
a(( (o?)cur)
((0c(o2)ann<o))
+ (at)eur
where (-) and (-) denote the probability density and cumulative distribution functions for the standard normal distribution. As result.
1, L(fl(o2)aan s ((("))(((")))e)) (ffl(22)aw , )pp.)))(( (((02) (((())))a))) (o2)cur)
Pat min
I(fl(o2)aa 3 (10 + (0
-mnn
(0)
(12)
(ff((22)cur Bar )p ((²2)cur) Φ (a
Size-and-shape preserving random effect under Prior Model 1. Recall that, under this prior model, (t) + a1 - t), a E (-1,1), e [0, The Unif(-1,1) prior on a; assigns zero mass to açan $ (-1,1) resulting in automatic rejection. In other cases, the proposal Unif(a. pur 8, atu + 6) symmetric. Thus,
Lfflagan con
Bour
Pa. = min
ko 21)
(13)
Lffla our Bour
Size-and-shape preserving random effect 7/ under Prior Model 2. Since
[ L(Shian Brur ( (l>55rr120 L(f(68" ( ) (p()) cur )9()) hn cur)
Pr: = min {
(14)
focus the derivation of g We omit the subscript to simplify notation. Let a(1) byurr
p(>) (T((2) - (00))) - (l). ((,,,))))(()2)))))) (1),,1)) (A,)) =: A(7)
(15)
The proposal distribution is A(7) Dirichlet(at), where t = (t2, 13 tz 1 tr = and has density
T -1
r(a)
II (A, j=1
(16)
I ((ate))
16

---

Next, we derive the density of (TTau) given (()) and Your, which appears in the denominator of the acceptance ratio. Since 7 = Your Yeans we have
(s) (5(tz), ................................ ((ts) ((t2)) ................................ R(1) ((tx,-1))
Your Your (Tean (12)); Your (Tca((at))))) - Your (Team(t2)) Your (Yeon (1)) - (tr (tr,-1)))
(7our (((aa)))) Your
As (1ran))
Tour (T)aa)))))
5=1
As(Tra)))
A (rran))
Your
Your
T
Ty-22
Σ (Tean))
Σ j=l
A (Team) =
Tour
Your
=: Yeur).
(17)
Therefore, the density of ( given A()) and Your the density of ( i.e., the density of Dirichlet(at), multiplied by the determinant of the Jacobian of The Jacobian of is given by
0g
0g
8A1
8Ar.
0ga 0AT,
0ga
8g 0
8A1
Ogn, -1
Og
8A1
21 Year (Ai (7ean))
0
11 Your
(E)=) 4 (Tean)
(18)
E jj=1 A, (Yean))
-1
Yowr
and is lower triangular matrix. Then, the determinant of this Jacobian is
T,-1
0g
Σ AT((n
(19)
A(A (Yean)
Yeur
The numerator of the acceptance ratio contains the density of A(Teur): given A(5-1) and Year, which the density of ((7-1) multiplied by the determinant of the Jacobian of == ((A(Ter))) 7-1), ) the transformation from A (Tter) to A(5-1). This determinant derived the same way as the determinant of the Jacobian of g. Finally, we have
I(a)
I (Aş(7
T
-1 11 Tcon
Ek-n Du (Tcar)
((7 can) ca (eur) )our)
IT, ((atg))
(20)
I(a)
I1 (A
fy- -1 i1 Your
Ek=1 Ao Al(Co)))
1=1
(((tg))
E Detailed Markov chain Monte Carlo algorithm
The detailed MCMC algorithm is given in Algorithm
Algorithm Bayesian Functional Mixed Effects Model
Input: Data: fi = f((tz), i = 1,..., n, j = 1, , T; prior hyperparameters: 7" (prior variance for a), ar, ba (shape and scale for o²), do.: b. (shape and scale for of), 0y, t (concentration and discretization in Prior Model 2 for 7); proposal hyperparameters: Σn (proposal covariance for a), + (proposal variance for a²), T (proposal variance for o?), 8 (proposal for Prior Model for 7), (concentration and discretization in Prior Model 2 proposal for 7); number of
17

---

burn-in iteration: No; total number of MCMC iterations: N; number of iterations between tuning
of proposal parameters: Nç; initial values: ag. (o")o: (02)o: (ai)o: I...., (Prior Model or (7i)o i = ., n (Prior Model 2).
Output: Posterior samples of all parameters: Ok:: (a²);: (o (as)k, 1,..., n (Prior Model 1) or (7)),, i = 1,
(Prior Model 2), h = 1,..., N Ng.
for j in 1:N do
if mj N+) = 0 and j < N₂ then Update La based on the empirical correlation matrix of the last Ne samples of a. end if
2. Propose a amn and compute Pa using
if U < Pa. U ~ Un, 1) then
3. Set aj+1 a can
else
Set aj+1 O.j-
end if
if mj,(j, No) = O and j < N₂ then
4. Update = based the acceptance rate of the last N, samples of g². end if
Propose (o2)can and compute Pa* using if U < pra, U ~ Unif(0. 1) then Set (o )+1 (02)
else
Set (a²) (o²)
end if
mjd(j, N₂) == 0 and j < N₂ then
7. Update To_ based on the acceptance rate of the last N, samples of o².
end
8. Propose (oz)can and compute Pog using 12
if U < Pog. U ~ Unif(0, 1) then
9. Set (a2) (o.)oan
else
9. Set (o2)i+1 (o
end if
for i in 1:m do
if Prior Model then
if' mod(j, N₁) == and 10. Update d based on the acceptance rate of the last N₃ samples of oq. end if
then
11. Propose agen and compute Pa. using 13
if U Pas, U Unif(0, 1) then
12. Set (a+)j+1 =a?". agan.
else
12. Set (at)j+1 (a:))-
end if
end if
Prior Model then
if mod(jj N,) = 0 and No then
10. Update a based on the acceptance rate of the last N₂ samples of 7/- end if
11. Propose zna and compute Pr, using if U P.c: U Unif(0, 1) then
12. Set (7:)++1 = 2eon.
else
12. Set (7))411 = (74)s-
end if
end if
end for
end for
18

---

(b)
(c)
(d)
-0.2
1.5
0.13
0.5
-0.4
0.12
0.11
0.4
-0.8
-0.6
0.1
0.3
-0.9
0.08
0.2
5
5
10
10
10
5
10
100
x10
-10
x10
-0.2
1.8
0.14
0.5
0.8
-0.4
1.2
0.12
D.4
0.6
0.11
0.4
-0.8
0.1
0.2
-1
0.8
000000
0.2
10
10
1n
v10
10
>1n
4.6
4.8
1.4
0.5
1.35
4.4
4.6
1.3
0.4
4.2
4.4
1.25
0.5
1.2
4
0
10
4.2
10
1.15
5
10
2
10
5
10
x10
x100
x101
x10
0.24
-0.1
0.025
0.8
-0.15
0.05
0.22
0.02
0.6
0.4
0.2
0.015
0.2
-01
0.25
0.18
0.01
5
10
10
10
10
x10
x10
x10
x10
x10
Figure 7: Row 1: Simulated data- row in Figure in Section Row 2: Simulated data row
2 in Figure Q Section 4.1 Row 3: Berkeley row in Figure in Section 4.2 Row 4: PQRST row 2 in Figure 4 in Section 4.2 Trace plots for the (a) first two fixed effect coefficients in a. respectively, (c) error process variance o², (d) variance of size-and-shape altering random effect o²: and (e) size-and-shape preserving random effect (phase function) for randomly chosen observation. Ground truth and posterior mean are marked in black and red, respectively.
F MCMC diagnostic plots for examples in Section 4
Figure 7 shows trace plots of 100, 000 MCMC iterations after the burn-in period for all model parameters. The first two rows correspond the simulated data examples considered in Figure in Section [4. 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 o² and o². respectively. Finally, panel (c) shows trace plots for the parameter a, when Prior Model 1 was used, or the parameter values (t2) (blue), (ttts) (yellow) and (ta) (purple) when Prior Model 2 was used (recall that 1- = 0.25, to 0.5 and t4 = 0.75). The choices of phase functions 7/2 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 horizontal red line. All trace plots suggest convergence to the stationary posterior distribution.
19

---

(b)
(c)
10
10
10
-10
10
-10
-20
-20
-20
0.5
0.5
0.5
10
5
0
0
-10
-10
-10
-15
0
0.5
0.5
0
0.5
30
20
20
20
0
0.5
1
0.5
0
0.5
0.5
0.5
0.5
Figure 8: Top to bottom: data simulated using Prior Model for phase functions, data simulated using Prior Model 2 for phase functions, Berkeley growth rate functions and PQRST complexes. (a) Functions (ʃi o 7,") / and average j (black). (b) Average j (black) and jı + Un (blue). (c) Same as (b), but using U2-
G Empirical FPCA basis for fixed effect function
Here, we provide empirical evidence behind the claim that the use of an FPCA basis deteriorates estimation performance for (see Contribution 3 in Section IT. We demonstrate this using four datasets considered Section the two simulated datasets from Example
Section 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 WEL2 min.e Σ min || (fi vx (ii)
compute the sample covariance K of
(fr o 7 7. where 7 arg min || (fr o /ill (assuming that each w; is sampled at f 0, f2, ......, 1-, (T), and (iii) apply singular value decomposition to the covariance matrix K = USU The columns U provide a data-driven orthonormal FPCA basis for I.?.
The number of FPCA basis functions B, used to specify in the proposed Bayesian model is selected based on % variation explained: 90% for the two simulated datasets (By 7 and B, for data simulated under Prior Models I and 2 for phase functions, respectively) and for the PQRST complexes (B, 6); 80% for the Berkeley growth rate functions B, 8. In all cases, we use B, = 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
First, in Figure 8 we show results of applying the FPCA procedure to each dataset. Panel (a) shows the functions (Si o 7/ / and zi black. Panels (b)&(c) display the variation in the data captured by the two leading FPCA basis functions: f in black with f + Uy, 1, 2 in blue (j in (b)
20

---

(b)
5
2500
4000
2000
3000
1500
2000
1000
1000
-10
-Centared Samp.
500
-Cerilered Meen
-15
Truth
8.1
0.15
0.28
0
2
0.5
Post. mean: 0.20
Post. mean:
3000
0
2000
2000
-10
Centered Post. Samp, Cantaradd Post. Moan
-15
0.1
0:15
0.2
D.6
0.8
1.2
0.5
Post. mean: 0.15
Post. mean: 0.75
3000
Centered Samp. -entered Post. Mean
3000
2000
10
1000
1000
M/-
1.6
10
12
Post. Mean: 1.56
Post. Mean: 9.72
Centered Post. Samp.
3000
3000
Canttred Post. Mean
2000
4
1000
1000
1/
-2
0.16 0.17 0.18 0.18
0.02
0.04
0.06
Post.
Figure 9: Top to bottom: data simulated using Prior Model for phase functions, data simulated using Prior Model 2 for phase functions, Berkeley growth rate functions and PQRST complexes. Rows 1&3: model with Prior Model 1 on T. Rows 2&4: model with Prior Model on Г. (a) Estimation of μ:: ground truth when available (black), centered posterior samples (blue), centered posterior mean (red). (c)&(c) Histograms of posterior samples for o² and og, respectively (posterior mean in red; ground truth black).
and in (c)). Notably, the FPCA basis appears to capture various sources of variation including noise, which undesirable if they are to be used to model the fixed effect function μ..
prouidesssviess estimation results for (a) (ground truth in black when available, centered posterior samples in blue, and centered posterior mean red), (b) o² (ground truth black when available. posterior mean in red), and (c) of (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 reasonable estimate of μ. This in turn results in overestimation of 2 and o? (posterior samples of o? 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.
H Sensitivity analysis for hyperparameter misspecification
assess sensitivity posterior inference to under- or over-specification three hyperparameters in the proposed model: By (number of basis functions used to model the fixed effect function μ)), B, (number basis functions used to model the size-and-shape altering random effect v;) and 0, (concentration hyperparameter in Prior Model (PM 2) on the size-and-shape preserving random
21

---

Bj 4, PM1
By 4, PM2
B, PM
B, 4, PM
0
6
5
-6
-5
-10
Contered Post, Mean
110
Oenieced Post Mene
-10
-ontered Post. Mean
-Centendd Fisi. Mean
-10
Coatered Post, Mous
-Trash
-Trith
-Trthh
Truth
Truth
-95%
-15
15
-15
-15
0.5
0.5
0
0.5
0.5
0.5
By = 6, PM 1
B = 6, PM2
B,
6, PM 1
B. = 6, PM 2
0, 30
0
-5
-10
-Gentered Post. Moon
10
-Genteed Fust Mea
-10
-Gentened Post Meon
11
-Centened Post. Mean
10
-Ceatered Post Meazi -Truth
-Tuuhh
-Bruk
-Truth
-Truth
-06% CI
-95% Cl
-0035 CI
-05% C1
-26% CI
-15
15
-15
16
0.5
0.5
By 10, PM
By = 10, PM 2
B,
10, PM
B, 10, PM
0g 100
0
-10
-Centered
10
-Genteeed Post Mene
10
-Cuntered Poet, Mean
-10
-Centered
10
Fost Meas!
-Truth
-Trolki
-Truth
-ruthh
-Truth
-96% CI
-90% CI
-05% CI
-05% CI
-96% CI
-15
15
15
0.5
0.5
0.5
0.5
0 0.5
Figure 10: Centered posterior mean (red) and 95% credible interval (dashed blue) for μ (ground truth, black) for different choices of By, B, and f, in Prior Model 2 (PM 2) on phase functions; PM 1 refers to Prior Model I on phasc. The data was generated using B, B. 6 and 0, 30 (for PM 2). Row I: Estimation results for under-specified values of hyperparameters. Row 2: Estimation results for correctly specified values of hyperparameters. Row 3: Estimation results for over-specified values hyperparameters.
effect or phase function %)). 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, modified Fourier basis functions for and B, = B-spline basis functions for each u; Each panel in Figure shows the centered estimate of the posterior mean for μ (red), associated 95% eredible interval (dashed blue), and the ground truth μ (black). Rows 1-3 show results for under-specified, correctly specified and ever-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 By, (ii) under- or over-specification of B,, and (iii) under- or wwer-peeiffication of 0y. On the other hand, when B, is under-specified. we are unable to accurately recover i 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 larger number of basis functions to model the fixed and shappe altering random effects.
I Additional real data examples
We provide estimation results for μ for several additional datasets that have been analyzed in previous literature:
1. Pinch force data (Ramsay al. 1995). which was also analyzed by Claeskens et al. 2021
2. Respiration data (Kuttek et al. 2013
3. Gait data Kurtek et al. 2013).
4. Signature acceleration data (Keepp and Ramsay 2008
5. Gene expression data SSrivastava et al. 201
22

---

(d)
3
3
3
3
2
2
2
2
1
0
0
0
0
1'0
0.5
0.5
0.5
0.5
2
2
2
2
-2
2
2
-2
0.5
10
0.5
10
0.5
0.5
-2
2
-2
-2
0.5
10
0.5
0.5
10
0.5
2
2
2
2
1.5
1.5
1.5
1.5
1
0.5
0.5
0.5
0.5
0
0
0
0
0.5
0.5
1
0.5
0
0.5
2
2
2
2
1
0
-1
-1
-2
-2
-2
-2
0.5
1
0.5
1 0
0.5
1
0.5
Figure 11: Estimation results for the pinch force, respiration, gait, signature acceleration and gene expression datasets (top bottom). (a) Data. (b)&(c) Centered posterior mean (black) and 95% credible interval (dashed blue) for when Prior Models and 2 are used for phase functions, respectively. (d) warpMix estimate.
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 By = basis functions for the pinch force and respiration datasets, and B, 12 basis functions for the gait,
23

---

signature acceleration and gene expression datasets. We use B, 6 B-spline basis functions for each 2 in all cases.
Figure aa) displays each dataset (top to bottom: pinch force, respiration, gait, signature acceleration, gene expression). Panels (b) and (c) in Figure TT show estimation results for l when Prior Models and are used for phase functions, respectively. We show the centered posterior mean in black with 95% credible interval displayed using blue dashed lines. The credible interval computed pointwise, using the 2.5% and 97.5% empirical quantiles of the centered posterior samples, (jr v j = 1...., N. The warpMix estimate is shown in Figure di Overall, the proposed Bayesian model is effective recovering μ.. Compared to warpMix, our estimates contain finer geometric features and more representative of the size-and-shape patterns in the data.
24

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We state the contributions of our paper in Section In summary, the contribu- tions are as follows.
(a) We infer the size-and-shape of a fixed effect function μ using Bayesian functional mixed effect model.
(b) We use informative prior distributions to model phase variation, which allows for regularization of the posterior distribution of μ.
(c) We interpret phase variation data-driven rotation of fixed orthonormal basis systems for the Hilbert representation space.
The claims in the abstract and introduction match the three points above, and are further elaborated in Section Empirical support via simulated and real data examples is provided Section
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
make certain modeling choices that can be interpreted as limitations: common discretization observed functional data, independence in discretized eπror process, Gaussian likelihood. These assumptions can be relaxed at some computation cost. We discuss additional limitations in Section 5 High computational cost associated with sampling from the posterior distribution is discussed there, with proposals of possible solutions.
Justification: In Section
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested few datasets with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
25

---

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 faireess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, 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 impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
Our paper does not include theoretical results, but we make modeling assump- tions with appropriate justifications. Appendices Caveendicessaandd include detailed derivations for the marginal likelihood used in our model and the MietropoiisHHasting acceptance ratio for MCMC sampling.
Justification:
Guidelines:
The answer NA means that the paper does include theoretical results. the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: We have clear description of all model parameters, prior distributions and proposal distributions for all numerical experiments presented in Section We provide detailed MCMC algorithm in appendix E] We also included code supplementary material and provide step-by-step instructions on how to reproduce the table and figures presented in the main part of the paper.
Guidelines:
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 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 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 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 replicate the results, access to a hosted model (c.g., in the case
26

---

of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm.
(b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g.. with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors welcome describe the particular way they provide for reproducibility, In the case of closed-source models, it may be that access to the model is limited in some way (e-g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: The supplementary material include data and code for reproducibility purposes. The README.1x file provides detailed step-by-step instructions to reproduce the results presented in the main part of the paper.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.. ec/
coblic/guides htte:/uuddoo/oooouiiiisssee CodeSubmissionno. aesionPolicy ooiicy for more details.
While encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ calps..c/public//gue gutps//eeee//ee/ CwwdeSummsssionnooao.
for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, 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 respectively.
Guidelines:
27

---

The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to 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, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Figures 2 and show posterior samples and posterior means allow for assessment of posterior uncertainty. Figure II in Appendix pxpendix[shows> shows 95% pointwise credible intervals to help with posterior uncertainty assessment.
Guidelines:
The answer NA means that the paper does include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
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).
should be clear whether the error bar is the standard deviation the standard error of the mean.
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 96% CI, if the hypothesis of Normality of errors 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 out of range (e.g. negative error rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: The computer resources are specified in detail in the README. file in the supplementary material and briefly summarized in Section 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 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 the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code Of Ethics
28

---

Question: Does the research conducted in the paper conform, in every respect. with the NeurIPS Code of Ethics
Answer: [Yes]
The research conducted in our paper does not involve human subjects. All of the real data examples considered in the main part of the paper and appendices use data that is publicly available.
Justification:
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g.. there special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification:
The main focus of our paper is to develop and demonstrate the effectiveness of a new Bayesian functional mixed model that can reliably recover the size-and-shape of fixed effect function μ.. While the proposed modeling framework may have future societal impacts when applied to more diverse datasets, all of the examples presented in the paper are proof-of-concept.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no socictal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g.. disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learss from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (c.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The proposed model and data used our paper does not have a high risk for misuse.
Guidelines:
29

---

The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We cite original papers for all of the datasets used in Section 4.2 and Appendix We also cite the paper associated with the R package warpMix 0.1.0, license GPL(>
3) htttt:://// attts::/ccran. aaaaan---eeeetttt wtg/web/gaaneeeeeeeeeeeeeeeeeeee wwtppaassswaapiiiiiaapppap for comparison purposes.
html that is used
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URI..
The name of the license (e.g., CC-BY 4.0) should be included for each
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: Our paper does not introduce new datasets or computing packages. The code included as supplementary material is for reproducibility purposes.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/mooel as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
30

---

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 research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Question:
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 research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the 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.
31


Reviews:
[{"review_id": "f8padNkl8G", "reviewer": "Reviewer_yD41", "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.", "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. \n \nThe 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. \n \nCould both basis functions be learned jointly with weight coefficients? \n \nIn 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. \n \nFigure 1 (b) and (c) are not cited in the introduction. \n \nThe 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. \n\nWhat about failure cases?", "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.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "+ A Bayesian functional mixed effects model with an unrestricted form of phase variation is presented. \n \n+ A novel perspective of inferring a size-and-shape function, a type of geometric constraint. \n \n+ In general, the paper is concise and direct. Most explanations and discussions are provided. It is easy to follow, up to some points in formulation.", "weaknesses": "- 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). \n \n- 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. \n \n- The authors could include a practical motivation of the study, where they could apply their approach and why it is important.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Notation, Figure 1, Sample Sizes** We agree that the notation is dense in certain sections. If accepted, we will try to simplify notation. We denote vectors and matrices in bold and functions and scalars in regular font. We will include dimensionality when appropriate. Figure 1 (b-c) shows the effect of the norm-preserving action $D_{\\gamma}$ on $f$ (referred to on line 26 in Section 2). These images are provided in Figure 1 due to the page limit. In all simulations, we use $n=30$ (line 289). Thus, sample sizes in simulations and real data examples are comparable. The $N = 100,000$ is the number of posterior samples that are used for inference (line 280).\n\n**Motivation** Incorporating random effects in neural networks is a popular technique to account for correlations in input data (e.g., Simchoni et al., Integrating Random Effects in Deep Neural Networks, 2023). In this context, our main motivation stems from ascertaining what may be reliably recovered by a Bayesian model when inputs are functions, and the objective is to estimate a population-level fixed effect function $\\mu$ in the presence of two different types of random effects: size-and-shape altering one ($v_i$) and size-and-shape preserving one ($\\gamma_i$, phase variation). Simulated and real data examples show that our model better recovers geometric features of $\\mu$ than the SOTA model warpMix. We expect the understanding gained via the simple Bayesian model to have important consequences when overparameterized neural networks are used with input functional data. We discuss and analyze two motivating datasets in Sec. 1 and 4.2. The proposed framework is general with many applications; results for other real datasets are in App. H.\n\n**Type of Basis and Joint Estimation** Our model allows use of *any* basis for the fixed effect $\\mu$ and the size-and-shape altering random effect $v_i$. *We clarify that, despite choosing the type of basis apriori, the norm-preserving action $D_{\\gamma}$ rotates the basis system toward the data, allowing us to learn a data-driven basis for $\\mu$. Thus, theoretically, the apriori choice of basis is unimportant for recovery of size-and-shape of $\\mu$, a crucial implication of our geometric formulation of the problem.* Below, we provide practical guidance on this issue.\n\nOrthonormality: The basis should be orthonormal since our model is motivated by viewing the norm-preserving action as a rotation of the basis system.\n\nData complexity, fixed effect: Choice of basis for $\\mu$ should be motivated by the structure of given data, e.g., we use a Fourier basis for growth rate data, but B-spline basis for PQRST complexes; this is based on prior belief of whether $\\mu$ contains periodicity (growth rate data) or sharp geometric features (PQRST complexes). We found that a functional PCA (FPCA) basis, estimated using the data, is ineffective in modeling $\\mu$ as it captures noise and small-scale variation (App. F).\n\nRandom effect: We prefer a basis with compact support, e.g., B-splines, to capture local variation and finer features as compared to the fixed effect.   \n\nWe appreciate the comment on joint estimation, which can be quite challenging: one needs to disentangle MCMC errors related to estimating the functional subspace (what is the optimal basis and how many?) from that of estimating the size-and-shape of $\\mu$. We will pursue this in the future for the fixed effect, which is the primary object of interest for inference, in two ways: by (i) placing a prior on the Stiefel manifold of orthonormal frames, or (ii) using ideas from Matuk et al., Bayesian Modeling of Nearly Mutually Orthogonal Processes, 2023.\n\n**Number of Basis Functions** This is chosen apriori, and can be informed by exploratory analyses, e.g., FPCA. We found that posterior inference is robust to (i) under- and over-specification of the number of basis functions for the size-and-shape altering random effect, and (ii) over-specification for the fixed effect. It is somewhat robust to under-specification for the fixed effect; see figure in the attached file. Thus, one may choose a large basis for the fixed and random effects. Alternatively, we could treat $B_f$ and $B_r$ as random and infer them. We leave this for the future since posterior inference becomes more complex: we may have to use Reversible Jump MCMC to allow the parameter space to change dimension.\n\n**warpMix Comparisons** We carefully chose simulations to ensure fair comparisons to warpMix. For quantitative evaluation, we use the criterion used in their work (line 318). Simulated Example 2 is particularly illuminating wherein the datasets are generated using the warpMix model, which uses a value-preserving action on $\\mu$ and $v_i$ rather than a norm-preserving one. Our model qualitatively/quantitatively outperforms warpMix in all considered simulations. We provide comparisons to warpMix on real data in Sec. 4.2 and App. H; our model recovers a fixed effect with more pronounced geometric features. Also, since we use a Bayesian model, we are able to quantify uncertainty via 95% credible intervals; this is not straightforward in warpMix (result is only a point estimate).\n\nComputationally, warpMix is more efficient: we rely on MCMC while warpMix uses iterative likelihood optimization. We discuss computational cost of our method in Section 5, lines 366-370 and mention that there is room for improvement in MCMC efficiency (most efficient MCMC was not our primary focus). For comparison, 10 iterations in warpMix take ~5 minutes as compared to ~111 minutes needed by MCMC for 100,000 posterior samples for inference.\n\n**Failure Cases/Limitations** As mentioned earlier, we cannot reliably recover the true fixed effect when the number of basis functions is under-specified. This is not unexpected since the ground truth $\\mu$ does not lie in the span of the basis functions used in the model. We discuss limitations in Section 5, including lack of theoretical support and computational cost.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "h6lw9DUMVs", "reviewer": "Reviewer_qjNQ", "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.", "questions": "## Minor comments / questions:\n\n\n\n1. 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?\n2. The introduction was somewhat confusing, as illustrated below. \n    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.\n    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? \n\n    To sum up, I think there are multiple problems motivated and highlighted along the way which somehow buries the key issue of phase variability. \n\n3. 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.\n4. Line 119: gamma in Gamma?", "limitations": "See above", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. The presentation over all is clear, self-contained, and precise. While this might sound standard in the past, I find these properties rare in papers these days, so I appreciate them even more. \n2. The formulation is carefully backed up by justifications for the choice of components (e.g., basis functions, prior distributions) as well as numerical experiments.", "weaknesses": "## Significance:\n\n\n\n1. 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,\n    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?\n    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”.\n2. 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. \n\n## Clarity:\n\n\n\n1. 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”? \n    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).\n    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.\n\n## Novelty:\n\n\nThe 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. \n\n\n1. 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.\n2. 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].\n\nSo 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.\n\n[A] R. I. Hartley, “In defense of the eight-point algorithm,” IEEE Transactions on Pattern Analysis and Machine Intelligence, 1997.\n\n[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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Significance, Size-and-Shape Preserving Transformations** Broadly, for *any* application involving functional data, our model has a two-fold motivation: (i) norm-preserving action $D_\\gamma$ on fixed effect $\\mu$ and size-and-shape preserving random effect $v_i$, expressed as linear combinations of orthonormal bases, allows us to rotate the basis systems to better align with observed data, and (ii) isometric group actions play a key role in finite and infinite-dimensional shape analysis to account for nuisance variation (Srivastava et al., Functional and Shape Data Analysis, 2016). For (i), modeling flexibility enabled by $D_\\gamma$ results in more parsimonious models, especially since $\\mu$ is represented via a finite, and potentially low-dimensional, linear combination of basis functions. This has applied implications since the type of basis is chosen apriori. The viewpoint in (ii) allows us to formulate models and carry-out computations in a geometric fashion.\n\n**Comparisons** Quantitative comparisons to warpMix on synthetic data generated via the warpMix model in Sec. 4.1 show that our model performs better in terms of recovery of $\\mu$. In Sec. 4.2, we consider two real datasets: growth rate functions and PQRST complexes. In the first, warpMix only recovers a single growth spurt while our model recovers two growth spurts, a smaller initial one and later pubertal one. In the second, warpMix fails to produce an estimate. We provide results for more real datasets in App. H. In general, warpMix oversmooths the estimate of $\\mu$ resulting in fewer geometric features than expected. Finally, posterior MCMC samples from our model can be used to quantify uncertainty, which is an important applied consideration (warpMix provides a point estimate).\n\n**Equivalence Classes** The example you provide demonstrates surjectivity of the operator from a value-preserving action, but such an operator is not unitary since it does not preserve the $\\mathbb L^2$ norm. Operator $D_\\gamma$ on the other hand is unitary for every $\\gamma \\in \\Gamma$. Let $[f]_n:=$\\{$f(\\gamma)\\sqrt{\\dot \\gamma},\\gamma\\in\\Gamma$\\} be the equivalence class of a function $f$ under the norm-preserving action and $[f]_v:=$\\{$f(\\gamma),\\gamma\\in\\Gamma$\\} be the equivalence class under the value-preserving one. Class $[f]_v$ contains functions that have the *same* sequence of ordinate values as $f$ - image $t \\mapsto f(t)$ is preserved. On the other hand, $[f]_n$ contains functions with new ordinate values (due to scaling factor $\\sqrt{\\dot \\gamma}$) that are originally not in $f$, and new extrema may be created. Intuitively, thus, $[f]_n$ is \"larger\" than $[f]_v$. A rigorous proof that compares sizes of two open sets (e.g., open balls) in two norm topologies induced by two quotient metrics (w.r.t. to the equivalence relation) can be provided. However, the following simple example with Brownian motion $W(t)$ illustrates the above point. Let $X(t)=W(c t)$ with $c>0$ be a time-changed process, corresponding to a value-preserving action by scaling of time on some interval in $(0,\\infty)$. Then, the law of $X$ is singular w.r.t. the law of $W$ for every value of $c$ except $c=1$. However, under the norm-preserving map, $X(t)=W(\\gamma(t))\\sqrt{\\dot \\gamma(t)}$ is absolutely continuous w.r.t. $W(t)$ for any $\\gamma \\in \\Gamma$ (Example 6.5.2 in Bogachev, Gaussian Measures, 1998).\n\n**Equivalent Solutions** Indeed, line 276 introduces a post-hoc constraint to choose an appropriate element of the equivalence class of $\\mu$ in a data-driven manner, ensuring that the resulting $\\mu$ is associated with identity phase $\\gamma_{id}(t)=t$; this step is necessary since functional data comes under arbitrary phase variation. This is commonly referred to as orbit centering where the average of phase functions is used to choose an orbit representative (Chapter 4 of Srivastava et al., Functional and Shape Data Analysis, 2016). *This choice is optimal with respect to minimizing the extrinsic Frechet variance of the $\\gamma_i$s.* An alternative that can be used in the presence of phase outliers is the Frechet median. We thank the reviewer for pointing us to the two references. For [A], we agree that choosing an appropriate equivalence class representative is difficult in the presence of varying levels of noise. In our model, we use finite-dimensional priors on phase, which help us constrain the size of equivalence classes. For [B], it is an interesting idea to recover the entire equivalence class instead of a representative. However, the parameter spaces for our model are (high or) infinite-dimensional (even though some dimension reduction is enforced via prior formulations). Under the one-dimensional Prior Model 1 on phase, one can potentially explore the entire equivalence class of the posterior mean of $\\mu$.\n\n**Frequentist Inference** Maximum likelihood estimation of parameters would yield the maximum aposteriori estimate corresponding to our model under appropriate non-informative prior choices. Despite higher computational cost, we favor Bayesian inference since it allows us to easily quantify posterior uncertainty.\n\n**Introduction** We agree that, while important, the main issue motivating our framework is not the data dimension. Instead, the motivation stems from the fact that children can undergo different numbers of growth spurts of different magnitudes at different times. This implies that the data contains phase variation, which is not accounted for in model (1) (an implicit assumption in (1) is that there is no time variation in growth spurts). Phase variation in our model acts as a size-and-shape preserving random effect that allows better alignment of the magnitude and timing of growth spurts across observations. If accepted, we will clarify our motivation. Methodological motivation of our model, as opposed to one such as warpMix with value-preserving phase variation, is given in lines 31-45.\n\nWe will use $\\gamma\\in\\Gamma$ in line 119.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "zVvO0qAT5i", "reviewer": "Reviewer_JKM1", "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.\n\n----\nPost-rebuttal: The rebuttal addressed most of my concerns and I would like to raise my score.", "questions": "1)\tNeural networks can also recover functions by sampling. What are the advantages of the proposed method compared to neural networks?\n2)\tHow can we recover arbitrary by a fixed model in Eq.(3)?", "limitations": "NA", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1）\tThis paper tries to solve a challenging task, recovering the size-and-shape of an unknown function by sampling. \n2）\tThe theoretical analysis is interesting, e.g., Derivation of Metropolis–Hastings acceptance ratio.\n3）\tThe authors also give the convergence analysis by numerical experiments.", "weaknesses": "1)\tIt seems that we should give different prior models for different functions.\n2)\tThe parameter sensitivity of $\\theta$ and the running time of the proposed method should be given.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Prior Models** We want to address this comment from two different perspectives. First, depending on the real data scenario, the choice of the type/number of basis functions in the model for the fixed effect function $\\mu$ and the size-and-shape altering random effect $v_i$ can be different, e.g., for growth rate functions we used the modified Fourier basis to model the fixed effect $\\mu$, while for the PQRST complexes we used B-splines. This choice is data dependent. At the same time, different types of basis functions can be used for the fixed effect function and the size-and-shape altering random effect. In general, we recommend the use of B-splines for the size-and-shape altering random effect since they have compact support and are able to capture local variation around the fixed effect function. In the end, the prior distributions are specified for the basis function coefficients.\n\n**Computational Cost and Sensitivity Analyses** Lines 366-370 provide information about the computational cost of MCMC sampling for the proposed Bayesian model. In short, we require ~111 minutes to yield 100,000 posterior samples for inference. The figure in the attached file contains sensitivity analyses with respect to (i) the number of basis functions used to model the fixed effect $\\mu$, (ii) the number of basis functions used to model the size-and-shape altering random effect $v_i$, and (iii) the value of the concentration hyperparameter $\\theta_\\gamma$ in Prior Model 2 for the phase functions (size-and-shape preserving random effect). Through these experiments, we show that posterior inference is robust in almost all scenarios, except when the number of basis functions to model the fixed effect is under-specified (this is not unexpected since the ground truth $\\mu$ does not lie in the span of the basis functions used in the model).\n\n**Neural Networks** We agree that neural networks can be used to approximate posterior distributions and their functionals; in other cases, they can aid in designing more efficient MCMC samplers by optimizing the proposal distribution such as in Li et al., A Neural Network MCMC Sampler That Maximizes Proposal Entropy, 2021. We have not explored these directions and leave them for future work. The primary goal of the current paper is to introduce a Bayesian functional mixed model, motivated by size-and-shape models in the context of shape analysis, that accounts for size-and-shape altering as well as size-and-shape preserving random effects. The novelty is in treating phase variation in functional data as a size-and-shape preserving random effect. We use standard MCMC for posterior inference, but plan to consider other approaches for approximating the posterior distribution under this Bayesian model in the future.\n\n**Generative Model (3)** In (3), we specify a generative model for functional data that incorporates a fixed effect function $\\mu$, a noise process $\\epsilon_i$, a size-and-shape altering random effect $v_i$, and a size-and-shape preserving random effect $\\gamma_i$. The ability to simulate arbitrary functional data based on the model in (3) depends on the assumptions, e.g., smoothness, etc., on each of these model components. That said, this is a very general model, and to the best of our knowledge, the first one to consider both, size-and-shape altering and size-and-shape preserving random effects.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "WdTyPHxZfx", "reviewer": "Reviewer_jjxT", "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.", "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": "All limitations are adequately addressed.", "rating": 5, "confidence": 1, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- Due to the strong modeling assumptions, on synthethic datasets, the method can recover the true signal even in situations where there is lots of noise and classical methods (such as GPs) may fail. \n- The methods seems to be useful to analyze the Berkley growth spurt data set.", "weaknesses": "- 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. \n- 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer's concern regarding accessibility to the broader machine learning community. If the manuscript is accepted, we plan to simplify notation as much as possible and provide more intuitive descriptions of some of the mathematical concepts. \n\nThere were several factors that motivated us to choose NeurIPS as the outlet for this work. \n\n**1.** 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, see e.g., Rao et al., Modern Non-linear Function-on-function Regression, 2023 and Rao et al., Nonlinear Functional Modeling Using Neural Networks, 2023. Furthermore, employing mixed models with random effects to better model correlated input data in neural networks is fast gaining traction, see e.g., Simchon et al., Integrating Random Effects in Deep Neural Networks, 2023.\n\n**2.** The field of shape analysis and more broadly geometric data analysis on quotient spaces, which provides part of the motivation behind the proposed modeling framework, falls under the umbrella of machine learning with broad applications in computer vision, graphics and biomedical imaging. \n\n**3.** In both, functional data analysis and shape analysis, geometry and invariance to nuisance variation play an important role and aid in model formulation and estimation. Indeed, these ideas can be applied more broadly in machine learning and there is increasing interest in the use of geometry and invariance in neural networks; see Bronstein et al., Geometric Deep Learning: Grids, Groups, Graphs, Geodesics, and Gauges, 2021 and references within.\n\n**4.** Our motivation stems from understanding the type of signal that can be recovered by a neural network or any probabilistic model from complex input data in the presence of nontrivial symmetries and geometric information. We have considered one specific instance of this, versions of which have been extensively studied in the statistics community.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "O1OhwRNY1E", "reviewer": "Reviewer_UDtq", "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.\n\nThe 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.\n\nEvaluation is shown both on simulated data and some real functional data examples.", "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.\n\nAdditional questions are on evaluation, which are described in the experiments section of the review.", "limitations": "The paper has a reasonable description of the limitations of the model and its assumptions in section 3.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper does a good job of introducing various temporal warping models, and their intuition on many example data. I found the mathematical and visual explanations of the differences between value preserving and norm preserving warps illuminating. The paper is well positioned to advance the field of functional data analysis as more and more applications get impacted by these methods.\n\nThe paper tries to motivate the presented mixed effects models by highlighting historical examples based on 2D shape analysis, which provides good context for the proposed developments. \n\nThe paper also has good strategies for actually implementing the infinite dimensional computations using projections onto basis functions such as B-splines.", "weaknesses": "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? \n\nFurther, 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. \n\nBeyond 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Noise Variance** The warpMix default noise variance parameter, \"sigmaEpsilonTilde\", which is set to $10^{-3}$ in the R implementation, is the variance of $\\theta_i$, a set of parameters in the model for the phase functions (Equation (6) in Claeskens et al., Nonlinear Mixed Effects Modeling and Warping for Functional Data, 2021). It is not the noise variance or the variance of the random effect ($\\sigma^2$ and $\\sigma^2_c$ in our model). In Simulated Example 2 in Sec. 4.1, we compare the proposed Bayesian model to warpMix. There, we actually simulate data using the warpMix model (with default parameter values). Note that, as in the warpMix model, the data is generated using the value-preserving action and not the norm-preserving action used in our model. Despite this, our model recovers the underlying fixed effect function $\\mu$ with higher accuracy than the warpMix model (see Table 1 in the manuscript and attached file). In this regard, all of our comparisons to warpMix are fair.\n\n**Novelty and Comparison to Existing Methods** To the best of our knowledge, the Bayesian functional mixed effects model proposed in our manuscript is the first to use the norm-preserving action as a size-and-shape preserving random effect on the original function space. The cited papers use the norm-preserving action on a transformed space, the space of square-root velocity functions, which corresponds to the value-preserving action on the original space (see Chapter 4 in Srivastava et al., Functional and Shape Data Analysis, 2016). In addition, none of those models incorporate size-and-shape altering random effects. The table in the attached file includes a quantitative comparison of the proposed modeling framework to the model of Cheng et al., Bayesian Registration of Functions and Curves, 2016 which is a mean+noise model and utilizes the aforementioned square-root velocity representation. As can be seen in this table, our model outperforms their approach, labeled BRFC, in terms of recovery of the mean (fixed effect) function $\\mu$ on all simulated data examples.\n\n**Evaluation of Estimation Results and Novelty** The question of whether the fixed effect function $\\mu$ can be reliably recovered is a highly studied problem within the functional data analysis literature, and sufficient conditions have only recently been identified. In this context, our conclusion that the geometric size-and-shape of $\\mu$ is recoverable is novel. This finding has clear implications in many machine learning applications with functional data as input to neural networks, which additionally incorporate a random effect in its architecture (Simchoni et al., Integrating Random Effects in Deep Neural Networks, 2023). However, incorporating downstream tasks, e.g., a classification layer, into the proposed model is something of interest and we plan to pursue this in future work. Our quantitative evaluations are based on simulated datasets where the true underlying fixed effect function $\\mu$ is known. In these cases, we outperform state-of-the-art competitors including warpMix and BRFC as seen in the table in the attached file. As mentioned by the reviewer, we rely on qualitative evaluations for the real data examples. Nonetheless, in most real data scenarios it is fairly clear that the fixed effect function estimated using the proposed model is \"better\" or more representative of patterns in the observed data than that recovered using warpMix; see Sec. 4.2 and App. H. For the PQRST complexes, the warpMix model fails to yield an estimate, and for many of the other datasets it tends to oversmooth expected prominent features of $\\mu$.\n\n**Alternative Model Formulation** The fixed effect function and size-and-shape altering random effect components of our probability model are specified using linear combinations of basis functions, *but not the discretized data itself*. The main motivation behind this choice is dimension reduction (sample size $n$ is typically smaller than number of time points per function), and this is the most common approach in specifying models for functional data. That said, an alternative approach would be to specify the model pointwise using parameter function evaluations as suggested by the reviewer. This, however, would drastically increase the dimension of the parameter space for the fixed effect function, complicating inference. The basis functions are only used to specify the prior distributions for the fixed effect function and the size-and-shape altering random effect. In that sense, basis functions allow us to enforce desired smoothness in the estimate of $\\mu$, which would be much more difficult if the fixed effect function was modeled pointwise. As mentioned by the reviewer, we tried to use a data-driven FPCA basis in our model; see App. F. However, we found that this basis contained small scale variation, even in the leading components, making it ineffective at modeling the fixed effect.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank all reviewers for their careful consideration of our manuscript and constructive comments.

**Significance and Motivation** Employing mixed models with random effects to better model correlated input data in neural networks is fast gaining traction (e.g., Simchon et al., Integrating Random Effects in Deep Neural Networks, 2023). In this context, when the inputs are functional data observed with *arbitrary* phase variability, it is important to understand what feature/property of the fixed effect function $\mu$ may be reliably recovered. Our main contribution is the proposal, and investigation, of a geometric reformulation of the problem: we use the phase component in the probabilistic model as a space-time unitary transformation, and then determine a data-driven optimal rotation of the model's coordinates to recover the size-and-shape of $\mu$, with corresponding uncertainty estimates. The property of the norm-preserving action being an infinite-dimensional rotation of $\mathbb{L}^2$ is key in the formulation of the proposed model; this action is an isometry under the $\mathbb{L}^2$ metric. Since the fixed effect function (and the size-and-shape altering random effect) are eventually represented using a finite-dimensional basis set, such rotations allow us to align the $\mathbb{L}^2$ coordinate system of model components to the coordinate system of each observation, thus resulting in a more parsimonious model. Quantitative and qualitative evaluations on simulated and real datasets confirm these claims.

**Implementation** Our main contribution is a novel geometric perspective on the problem of reliably recovering the size-and-shape of a fixed effect function $\mu$ using a Bayesian functional mixed model with *unconstrained* (infinite-dimensional) size-and-shape altering and preserving random effects, in the presence of noise. As such, the computational algorithm used for posterior inference via MCMC serves more as a proof of concept rather than a definitive (optimal) computational tool, but nevertheless outperforms the current state-of-the-art (SOTA). There are many alternative algorithms that can be used to approximate the posterior distribution including variational approaches, and neural networks for intractable likelihoods and MCMC (e.g., Li et al., A Neural Network MCMC Sampler That Maximizes Proposal Entropy, 2021).

**Additional Results in Attached File** The attached PDF file contains a figure and a table. In the figure, we show results of more extensive sensitivity analyses to misspecification of three hyperparameters: $B_f$ (number of basis functions used to model the fixed effect function $\mu$), $B_r$ (number of basis functions used to model the size-and-shape altering random effect $v_i$) and $\theta_\gamma$ (concentration hyperparameter in Prior Model 2 (PM 2) on the size-and-shape preserving random effect or phase function $\gamma_i$). Each plot shows the centered estimate of the posterior mean (red), 95% credible interval (dashed blue) and ground truth $\mu$ (black). Rows 1-3 show results for under-specified, correctly specified and over-specified values of hyperparameters, respectively. Posterior inference, in terms of the posterior mean for $\mu$ 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 $\theta_\gamma$. When $B_f$ is under-specified, we are unable to reliably recover the ground truth $\mu$ since it 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.
The table reports a more comprehensive quantitative evaluation, in terms of estimation accuracy for the fixed effect function $\mu$, on five simulated datasets. Rows 1-2 consider data simulated from our model under Prior Model 1 (PM 1) and 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 refers to the prior model on phase and the letters F or B refer to the Fourier or B-spline bases) to those produced using warpMix (Claeskens et al., Nonlinear Mixed Effects Modeling and Warping for Functional Data Using B-splines, 2021) and the mean+noise model proposed in Cheng et al., Bayesian Registration of Functions and Curves, 2016 (labeled as BRFC). As seen in the table, our model outperforms warpMix and BRFC in all of these simulation scenarios.

Meta-review: 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.
Decision: Accept (poster)