Confidence Tubes for Curves on SO(3) and Identification of Subject-Specific Gait Change after Kneeling

In order to identify changes of gait patterns, e.g. due to prolonged occupational kneeling, which is believed to be major risk factor, among others, for the development of knee osteoarthritis, we develop confidence tubes for curves following a Gaussian perturbation model on SO(3). These are based on an application of the Gaussian kinematic formula to a process of Hotelling statistics and we approximate them by a computible version, for which we show convergence. Simulations endorse our method, which in application to gait curves from eight volunteers undergoing kneeling tasks, identifies phases of the gait cycle that have changed due to kneeling tasks. We find that after kneeling, deviation from normal gait is stronger, in particular for older aged male volunteers. Notably our method adjusts for different walking speeds and marker replacement at different visits.


Introduction
There is overwhelming evidence that prolonged occupational kneeling (POK), e.g.floor tile laying, constitutes a major risk factor for the development of knee osteoarthritis, e.g.Cooper et al. (1994); Coggon et al. (2000); Rytter et al. (2009).Also, POK is a risk factor for the development of degenerative tears in medial menisci, e.g.Rytter et al. (2009).In order to identify hypothesized underlying changes of gait patterns, kneeling workers' and controls' gait has been compared by Gaudreault et al. (2013) and prolonged kneeling has been simulated and gait changes compared by Kajaks and Costigan (2015); Tennant et al. (2018).Also, dependence of kneeling effects due to footwear has been investigated by Tennant et al. (2015) and kneeling effects have been studied on cadavers with total knee arthroplasty Wilkens et al. (2007).
In order to assess the specifics of changes of gait patterns, the three dimensional rotational path in SO(3) of the relative motion of the tibia (larger lower leg bone) w.r.t. the femur (upper leg bone) is usually represented by the three Euler angles flexion/extension, adduction/abduction and internal/external rotation.Doing so, Gaudreault et al. (2013); Kajaks and Costigan (2015); Tennant et al. (2018) have found, among others, for each angle, loci of significant gait changes, without, however, addressing the issue of multiple testing, correlation of the sequential data and the effect of marker replacement.
In our approach, we address all of these issues, and in consequence, are able to test for subjectspecific changes of gait pattern.In application, we do this for pre-and post-kneeling, the method, however, is applicable for any change of condition (e.g.onset of otheoathritis) over a period of time, due to correcting for marker replacement.To this end, we recall a Gaussian perturbation model from Telschow et al. (2016) for curves on Lie groups and show that a Hotelling statistic for the corresponding process follows asymptotically (for vanishing variance) a Hotelling statistic that can be described by a Gaussian kinematic formula (GFK) from Taylor et al. (2005); Taylor (2006).In application to gait analysis, our method, relying on curves on the rotational group, takes advantage of simultaneously involving all three Euler angles in a canonical way.Moreover, as our test statics use maxima of stochastic processes, we resolve the multiple testing issue by providing for simultaneous confidence tubes over entire gait cycles.Further, sequential correlation is naturally modeled within the GFK approach by simulating quantiles from the empirical process.Indeed, simulations mimicking and going beyond the use case of low variance and high smoothness typical in gait analysis show that our method is well applicable.Then, for an experiment conducted in the School of Rehabilitation Science at McMaster University (Canada), for six out of eight healthy volunteers we identify individual changes of gait patterns after kneeling tasks.We find that after kneeling, deviation from normal gait is stronger, in particular for older aged male volunteers.
Remarkably, our method is also robust under specialist marker replacement, at a subsequent patient's visit, say.It is well known that Euler angle curves may considerable change after marker replacement and simple approaches subtracting average angles over gait cycles (cf.Kadaba et al. (1989)) have remained questionable, e.g.Delval et al. (2008);McGinley et al. (2009); Noehren et al. (2010); Røislien et al. (2012), also for other approaches, leading to the longstanding open problem of gait reproducibility, see Duhamel et al. (2004).
In a recent publication (Telschow et al. (2016)), we have developed a method to successfully correct for marker replacement by estimating a Lie group isometry, bringing two samples (each sample is a repeated measurement of the same person's gait with fixed marker placement) into optimal position to one another.Since volunteers will have different comfortable walking speeds at different visits, we have also corrected for a sample-specific time warping effect.This method is part of the tool chain developed in the present contribution which is available under www.stochastik.math.uni-goettingen.de/KneeMotionAnalyticsas an R-package.In particular, it contains all data and code used in this paper.

Testing Gaussian Perturbation Models on Lie Groups
Modulo Sample-Specific Spatio-Temporal Action The following is taken, from Telschow et al. (2016).It has been formulated for p = 1 and generalizes at once to arbitrary p ∈ N. Let G be a connected Lie group with Lie algebra g embedded in a suitable Euclidean space R m and Lie exponential Exp : g → G.With the unit interval I = [0, 1] we have the family Γ = C p ([0, 1], G) of p ∈ N times continuously differentiable curves on G.We assume in particular that G admits a bi-invariant Riemannian metric, a sufficient condition for which is that G is compact.
Definition 2.1.We say that a random curve γ ∈ Γ follows a Gaussian perturbation (GP) around a center curve γ 0 ∈ Γ if there is a g-valued zero-mean Gaussian process A t with a.s.C p paths, such that γ(t) = γ 0 (t)Exp A t for all t ∈ I . (1) The Gaussian process A t will be called the generating process.
This model, which is based on right multiplication with the exponential of the generating process is equivalent to one based on left multiplication and asymptotically (as the variance goes to zero) equivalent to one based on two-sided multiplication, cf.Telschow et al. (2016).Moreover, this model is invariant under the spatial action of the isometry group I(G) on G, where its connectivity component I 0 G of the identity element can be viewed as the analog of the orientation preserving Euclidean motions of a Euclidean space.Indeed, for G compact, semisimple with trivial center, we have G × G = I 0 G , cf.Telschow et al. (2016), which is the case for G = SO(3).

Confidence Tubes on G
Since G is connected by hypothesis, the inverse exponential is well defined on the complement in G of the cut locus of the unit element.Let Log : G → g denote a measurable extension.Further, since g is a linear space, let ι : R m → g be a suitable isomorphism and set Definition 3.1.Let γ 1 , . . ., γ N be a sample of a random curve γ ∈ Γ following a GP model around a center curve γ 0 and let γN be an estimator for γ 0 .Then x are called intrinsic population and sample residuals, respectively.
This gives rise to the following one-dimensional processes, where we assume that Ŝx,N From this we obtain at once simultaneous (1 − α)-confidence tubes for γ 0 , setting Theorem 3.2.Let γ 1 , . . ., γ N be a sample of a random curve γ ∈ Γ following a GP model around a center curve γ 0 .Let γN be an estimator for γ 0 and assume Ŝx,N t is non-singular for all t ∈ [0, 1].
The process Ĥx,N t from (5) serves as an approximation of the genuine Hotelling process determined by ( 1): where Among others, the following section makes this approximation explicit for the special case of G = SO(3).

GP Models and Approximating Confidence Tubes on SO(3)
For G, the compact and connected Lie group of three-dimensional rotations G = SO(3) we detail the above approximation.To this end, we first recall the structure of SO(3), extrinsic pointwise means as estimators γN and fundamental properties of corresponding GP models.

GP Models on SO(3)
The Lie group G = SO(3) comes with the Lie algebra g = so(3) = {A ∈ R 3×3 : A T = −A} of 3 × 3 skew symmetric matrices.This Lie algebra is a three-dimensional linear subspace of all 3 × 3 matrices and thus carries the natural structure of R 3 conveyed by the isomorphism ι : R 3 → so(3) given by This isomorphism exhibits at once the following relation We use the scalar product A, B := trace Bhattacharya and Patrangenaru (2003).Moreover, we denote with I 3×3 the unit matrix.As usual, A → Exp(A) denotes the matrix exponential which is identical to the Lie exponential and gives a surjection g → G. Due to skew symmetry, the following Rodriguez formula holds This yields that the Lie exponential is bijective on B π (0) = {A ∈ g : A F < π}.For a detailed discussion, see (Chirikjian and Kyatkin 2000, p. 121 then the following hold. (i) There is Ω ⊂ Ω measurable with P(Ω ) = 1 such that for every ω ∈ Ω there is N ω ∈ N such that for all N ≥ N ω , every ÊN (t) has a unique element γN (t), for all t ∈ [0, 1], and γN ∈ Γ; Corollary 4.2.With the notations and assumptions of Theorem 4.1 we have lim

Approximating Confidence Tubes on SO(3)
As the main result of this section we first show that in case of concentrated errors (as is typical in biomechanics, e.g.Rancourt et al. (2000)), the residual processes x N t and X N,n t from Definition 3.1 are approximatively the residuals of the generating Gaussian process (7) of the GP model.Then, we use this approximation to define an estimator for ĥγ,N,α based on the Gaussian kinematic formula for Hotelling processes (see Taylor and Worsley (2008)), which will be shown in Section 5, using simulations, to perform very reliably even if the sample sizes are small as it is usually the case in biomechanical gait analysis.
Theorem 4.3 (Approximations for Concentrated Errors).Let N ∈ N be fixed and γ 1 , . . ., γ N be a sample of a random curve γ ∈ Γ following a GP model around a center curve γ 0 .Additionally, assume that the generating Gaussian process Let γN (t) be a measurable selection of sample PEM curves.Then, for x N t and x N,n t from Definition 3.1, where Corollary 4.4 (Asymptotically genuine Hotelling process).With the assumptions and notations of Theorem 4.3 we obtain with O p (σ) uniformly over t ∈ [0, 1], if additionally cov a t = σ 2 Σ t with fixed and non-singular Σ t for all t ∈ [0, 1].
The following theorem gives the equivariance property of the simultaneous confidence tubes with respect to the group action on Γ of the group Theorem 4.5.Let γ 1 , . . ., γ N be a sample of a random curve γ ∈ Γ following a GP model around a center curve γ 0 with PEM curve γN .Moreover, let (ψ, φ) ∈ (G × G) × Diff + ( I ) be arbitrary and define the sample i.e., they can be derived from the simultaneous confidence tubes for γ 0 using γ 1 , ..., γ N and (ψ, φ) ∈ (G × G) × Diff + ( I ) only.
The Gaussian kinematic formula (GKF).Corollary 4.4 states that for concentrated errors the statistic H a t , which is the Hotelling T 2 statistic of a generating Gaussian process, approximates the statistic Ĥx,N t .Thus, in order to estimate the quantiles ĥγ,N,α for the process Ĥx,N t , derived from a GP model γ, we use the expected Euler characteristic heuristic (see Taylor et al. (2005)) and assume that where x(U) denotes the Euler characteristic (EC) of U ⊂ [0, 1].Although we cannot rigorously justify this approximation, our simulations in Section 5 show that this procedure works very well.
Under some additional technical assumptions on the generating Gaussian process A t = ι • a t given in Taylor (2006), it is shown in Taylor and Worsley (2008) that the expected EC of the with the so called Lipschitz-Killing curvatures The so called Euler characteristic densities ρ H j for j ∈ {1, 2} appearing in the GKF (15) can be computed from the EC densities of a T -process with N − 1 degrees of freedom via Roy's union intersection principle (cf.Taylor and Worsley (2008, Sec.3.1.))using the formula Here µ d (S 2 ) denotes the d-dimensional intrinsic volume of the two-sphere S 2 given by in Taylor and Worsley (2008, p. 23).In relation to the Stochastic Geometry literature, µ 0 gives twice the number of connected components and µ 2 gives the surface area of S 2 (e.g., Mecke and Stoyan (2000, p. 100)).Moreover, the EC densities of a T -process with (N − 1) degrees of freedom have the explicit representations , given in Taylor and Worsley (2007, p. 915).
Estimation of the quantile ĥγ,N,α .Using the GKF for Hotelling T 2 -processes together with the EC heuristic ( 14) yields which can be used if L 1 [0, 1] is known, to estimate the value ĥα,N,α for low probabilities α by solving 2ρ T 0 ( Thus, it remains to estimate the Lipschitz-Killing curvature L 1 [0, 1] .This has been achieved for Gaussian processes in R D , D ∈ N, in Taylor and Worsley (2007, Sect. 4) and Taylor and Worsley (2008), where they also proved that their estimator is consistent.By Theorem 4.3 the intrinsic residuals of a sample from a GP model γ are, in case of concentrated errors, close to the residuals of the generating Gaussian process A t = ι • a t .Since the estimator of Taylor and Worsley (2008, Equation (18)) is based only on the Gaussian residuals, we adapt their estimator by replacing their residuals by the intrinsic residuals given in Theorem 4.3 to obtain an estimator of the Lipschitz-Killing curvature L 1 [0, 1] .
For convenience we restate the resulting estimator.Let γ 1 , . . ., γ N be a sample of a GP model γ and assume the curves are observed at times 0 = t 1 < t 2 < ... < t K = 1.Then we define the matrix Further, denote by R d t k the d-th column of R t k and define the normalized residuals as for d ∈ {1, 2, 3} and k ∈ {1, . . ., K}.The estimator of the Lipschitz-Killing curvature L 1 [0, 1] is then given by

Simulations of Covering Rates
Since the estimation of the quantile ĥγ,N,α relies on an approximation for concentrated error processes given in Theorem 4.3, we study the actual covering rate of this method using simulations.
GP models used for simulation.Without loss of generality we may assume that our center curves satisfy γ 0 (t) = I 3×3 for all t ∈ [0, 1].Otherwise, multiply the sample with γ 0 (t) −1 .
From these error processes the generating Gaussian process A t of the GP model is constructed by the following formula A i,l,j,σ t = M j σε i,l 1,t , σε i,l 2,t , σε i,l for i ∈ {1, 2, 3}, j ∈ {1, 2}, l ∈ {1, 2, 3} and σ ∈ R >0 .Here we denote with ε i,l s,t for s = 1, 2, 3 independent realizations of ε i,l t t∈I .The matrices are introduced to include correlations among the coordinates.Moreover, ( 19) introduces different variances in the coordinates, since for j = 2 the second component has half the variance of the other two components.
Then (1 − α)-SCT are constructed using Theorem 3.2.Here the quantile ĥγ,N,α is estimated by equation ( 16) using the estimator (17) for the Lipschitz killing curvature.Afterwards it is checked whether γ 0 ≡ I 3×3 is contained in the SCT for all t ∈ T .This procedure is repeated M = 5000 times.The true covering rate is approximated by the relative frequency of the numbers of simulations, in which the constructed SCT contained the true center curve γ 0 for all t ∈ T .
Results of simulation of SCT for center curves of GP models.The results are reported in Table 1 and they convey a positive message: For a variance σ = 0.05, which is that of the data of the application in Section 6, the simulated covering rate is very close to (1 − α).Only in the case of the Ornstein-Uhlenbeck error process we have slightly too high covering rates.For higher variance (σ = 0.6) we underestimate the covering rate.This is expected, since the proposed estimator is designed for concentrated data and the map v → Log Exp(v) is only the identity on v < π and we have the inequality This implies that our estimated covariance matrix has smaller eigenvalues then the covariance matrix of the sample and hence our confidence sets will become smaller.This effect is more visible if the sample size is large, since more curves cross the cut locus.
6 Application: Assessing Kneeling Effects on Gait Study design.In a study conducted at the School of Rehabilitation Science (McMaster University, Canada), 8 volunteers (4 female , 4 male, for each gender, two aged 20-30 and two aged 50-60) with no previous knee injuries (external observation and subjective questioning revealed no obvious knee problems) with unremarkable knee kinematics motion have been selected.In the experiment retro-reflective markers were placed onto identifiable skin locations on upper and lower volunteers' legs by an experienced technician following a standard protocol.Eight cameras recorded the position of the markers and from their motions, a moving orthogonal frame E u (t) ∈ SO(3) describing the rotation of the upper leg w.r.t. the laboratory's fixed coordinate system was determined, and one for the lower leg, E l (t) ∈ SO(3), each of which was aligned near I 3×3 when the subject stood straight.As is common practice in clinical settings, subjects walked along a pre-defined 10 meter straight path at comfortable speed.For each of the following four sessions (A,B,C,D), for each subject a sample of N ≈ 12 (for details on N , see Table 2) repeated walks have been conducted and for every walk a single gait cycle γ(t) = E u (t)E l (t) T about half way through has been recorded, representing the motion of the upper leg w.r.t. the lower leg.
After each walk the volunteers stopped shortly and started again for the next 10 meter walk.Thus, by design the assumption of independence of recorded gait cycles is satisfied.0% 25% 50% 75% 100% A 11.00 12.00 12.00 13.00 14.00 B 12.00 12.00 13.00 13.25 14.00 C 9.00 11.75 12.00 12.25 14.00 D 9.00 11.00 12.00 12.25 13.00 Table 2: Reporting the quartiles of numbers of processed walks (gait cycles) of volunteers for each of the four sessions from Table 3.
The study consists of four sessions, each giving, as described above, a sample of walks for the left leg of each volunteer.Between samples A and B the markers were detached and placed again by the same technician following the same standard protocol.Hence the difference between these samples reflects the challenge of repeated reproducibility of gait patterns under clinical conditions.Before conducting the two sessions C and D markers where again replaced and the volunteers fulfilled a task of 15 minutes kneeling prior to data collection of session C and yet another 15 minutes kneeling prior to session D. This allows to study the effect of kneeling and prolonged kneeling on gait patterns.Table 3 gives an overview of the four sessions conducted.Sessions A and B have already been reported in Telschow et al. (2016) Estimation of P and Q and temporal alignment of the sample mean curves have been done as described in Section 2 and detailed in Telschow et al. (2016) in order to make the samples comparable.Indeed, by Theorem 4.5 the shape of the confidence sets does not depend on alignment correction.
In the following we report our findings, first in Table 4 using the permutation test from Telschow et al. (2016, Test 2.11) correcting for sample-specific group action.If we were not to correct for sample-specific group action, we would detect significant changes of gait for 6 out of the 8 volunteers, even for "A vs. B", where nothing changed but marker placement, cf.Telschow et al. (2016, Table 2) .The challenge dealt with in Telschow et al. (2016)  Results.In Table 4, we see significant (often highly significant) changes of gait of volunteers 2, 4, 6 and 7 after each of the kneeling tasks.Volunteers 3, 5 and 8 show no changes.Remarkably, these findings are consistent over marker replacement ("A vs. *" and "B vs. *") and only for Volunteer 1 the picture is unclear.
In order to locate changes of gait patterns, we apply our new test of simultaneous confidence tubes.In Table 5 we report the specific loci where 1 − α = 0.95 confidence tubes no longer overlap, using standard naming convention (e.g.Rodgers (1995)) as illustrated in Figure 1.Employing Euler angles, which are popular in the field, as a local chart of SO(3), the corresponding curves and specific loci of non-overlapping simultaneous confidence tubes are shown exemplary in Figures 2 for Volunteer 2 and in Figure 3 for Volunteer 6. Notably, non-overlapping confidence tubes have been determined in SO(3) and not in chart coordinates so that the chart representations only serve as an approximate visualization of the real situation which we cannot visualize.The other volunteers' (1, 3, 4 and 7) curves with loci of non-overlapping confidence tubes are shown in the appendix in Figure 4. Again, we see that Volunteers 5 and 8 feature no changes in gait pattern.Volunteer 7 reported physical pain after post-kneeling walking.Indeed, high variation in gait patterns corresponding to session D (red, in the left two displays of the bottom row in Figure 4) widened the corresponding confidence tubes such that changes of gait in Session D were not detected.
Combining Tables 4 and 5 and taking into account age and gender, we see that older age (volunteers with even numbers belong to age group 50 -60) favors a kneeling effect over young age (volunteers with odd numbers belong to age group 20 -30).As a surprise, the effect seems to be overall stronger for males.Having established a tool chain to study such effects, this experiment warrants larger studies.

Discussion
In conjunction with the permutation test and estimation of marker replacement effects from Telschow et al. (2016), with the test for simultaneous non-overlapping confidence tubes presented in this paper, we have developed a tool chain that can be used in clinical practice to assess changes of gait patterns and localize these.These are no longer based on (single) Euler angle representations, as are often used in the field, but take advantage of a Gaussian perturbation model defined in the Lie group of three dimensional rotations.Due to the conservation of moment, gait curves are naturally smooth, their variation over repeated walks is moderate and hence approximations via the Gaussian kinematic formula are rather accurate, as well as in theory as in practice.
In this study, with a small number of participants and a small number of repeated walks, we see that short kneeling tasks tend to affect gait patterns and it seems that older age and, possibly, male gender, favor this effect.We have made sure that this effect has not been caused by different marker placements.While specific loci of gait change depend on individuals, changes seem to occur least at local maxima of dominating flexion-extension, namely at MF and MS.
We believe that our results derived for G = SO(3) generalize to general connected Lie groups

Figure 1 :
Figure 1: Depicting standard naming convention for gait events with respect to the flexionextension angle.

Table 3 :
. Experiments conducted Dealing with the marker replacement effect.Replacing markers between sessions results in fixed and different rotations of the upper and lower leg, conveyed by suitable P, Q ∈ SO(3) such that E after u (t) = P E before u (t) and E after l

Table 4 :
was to design a test keeping the level, also under marker replacement.Reporting p-values (significant in bold face) obtained from the permutation test inTelschow et al. (2016, Test 2.11in the version of Remark 2.12) correcting for sample-specific group action.