Classification in Postural Style Based on Stochastic Process Modeling

2.1 Original dataset

The dataset was collected at the center for the study of sensorimotor functioning (CESEM) of Université Paris Descartes. It is composed of a cohort of n=70 subjects. Among the 70 subjects, 22 are hemiplegic (due to cerebrovascular accidents), 16 are vestibular, while the 32 remaining subjects are identified as normal based on an initial medical evaluation.

Each subject completed two protocols designed to evaluate potential deficits in maintaining posture. These protocols have been identified as the most informative among four similar protocols for classifying hemiplegic versus normal subjects in the earlier study [1]. Both protocols are divided into three phases: a first phase (lasting 15 s) with no postural perturbation, a second phase (lasting 35 s) with postural perturbations (either some muscular perturbations or a combination of muscular and visual perturbations), and a third phase (lasting 20 s) with no postural perturbation. We expect that the subject’s postural sway changes around the beginning and the end of the perturbation phase (around 15 and 50 s).

For each protocol, the center-of-maximal-pressure exerted by each foot on a force-platform is recorded at equispaced discrete times. Thus each protocol results in a trajectory (Li,Ri)i≤N, where (Li,Ri) is the observation at time ti=iδ for i=1,…,N=2,800 and δ=0.025 s. For each i=1,…,N, Li=(Li1,Li2)∈ℝ2 and Ri=(Ri1,Ri2)∈ℝ2 respectively correspond to the left and right foot (Table 1).

Following Chambaz and Denis [1], we derive from (Li,Ri)i≤N the one-dimensional trajectory C1:N (sometimes we will denote a n-tuple (x1,…,xn) by x1:n) characterized by

where γ is defined as the component-wise median value of (12(Li+Ri))i≤400 over the ten first seconds of the protocol. The process C1:N provides a relevant description of the sway of the body during the course of the protocol. In Figure 1, we display the trajectories C1:N corresponding to protocol 1 (left plot) and protocol 2 (right plot) as completed by a single hemiplegic subject. Figure 1 confirms the intuition that the subject’s postural sway is not necessarily instantaneously affected by the start and end of perturbations.

Figure 1

Two trajectories C1:N respectively corresponding to protocol 1 (left) and protocol 2 (right) undergone by a single hemiplegic subject. Vertical lines correspond to the estimated change points (τˆ1δ,τˆ2δ) as obtained based on the two (real) observed trajectories

2.2 Data modeling

We model the trajectory C1:N as the observation at discrete times of a stochastic process (C(t))t∈[δ,Nδ] characterized by a stochastic differential equation.

We model the effects of the perturbations by possible changes in the volatility and drift functions at two unknown change points. In order to account for the fact that subjects react differently, we also assume that the change points differ among subjects. We denote by (T1,T2) the change points of a trajectory, with T0=τ0δ=δ<T1=τ1δ<T2=τ2δ<T3=70 (τ0=1, τ1 and τ2 are integers). We assume that there exist two functions a, b, and two parameters (ϕ1,ϕ2,ϕ3) and (σ1,σ2,σ3) such that, for all t∈[Tk−1,Tk) and k=1,2,3,

where (W(t))t∈[δ,Nδ] is a standard Wiener process.

We now specify the parametric forms of the volatility and drift functions a and b. Because Figure 1 indicates that the variance of C1:N is not constant over time and because, in addition, the process C1:N takes positive values, we decide to rely on the classical Cox–Ingersoll–Ross (CIR) process [6], by setting a(x,ϕ=(λ,μ))=λ(μ−x) and b(x,σ)=σx (for all x≥0).

In summary, we model the trajectory C1:N as Ci=C(iδ)=C(ti) the observation at discrete times of the stochastic process (C(t))t∈[δ,Nδ] being characterized, for t∈[Tk−1,Tk)(k=1,2,3), by

where the initial condition C(Tk−1) is the observation Cτk−1 at time Tk−1, λk, and μk and σk are positive constants. It is known that if 2μkλk/σk2≥1 then the process remains positive and admits a stationary distribution, namely a Gamma distribution with shape parameter 2μkλk/σk2 and scale parameter σk2/2λk [6]. We summarize all parameters as a single vector (τ1,τ2,θ1,θ2,θ3) with θk=(λk,μk,σk) for k=1,2,3. We assume that the changes, induced by the perturbation phase, necessarily affect the drift parameter through a change inμ. In other words, we assume that μ1≠μ2 and μ2≠μ3. Changes can also affect the other parameters.

3.1 A two-step estimating procedure

We define a two-step estimating procedure: we first estimate the sequence of change points, and then estimate (θ1,θ2,θ3) on each interval characterized by the previously estimated change points.

A great variety of methods have been proposed for the estimation of change points (among many others [7–9–11]. We choose to adopt a popular approach originally proposed by Bai [7]. The approach relies on a least-squares criterion and aims at detecting change points which affect the mean of a linear process. The estimator (τˆ1,τˆ2) of the sequence of change points (τ1,τ2) is defined as:

where Cˉk is the arithmetic mean of Cτk−1:τk−1. This makes sense because if (C(t))t∈[Tk−1,Tk) reaches a stationary regime then the random variables Cτk,…,Cτk−1 are identically distributed from a stationary distribution whose mean parameter is μk.

Fix k∈{1,2,3}. On each interval [Tˆk−1,Tˆk)=[τˆ1δ,τˆ2δ), we estimate θk by minimizing a contrast function based on the log-likelihood of the approximated discrete-time Euler–Maruyama scheme (see Kessler [12], for instance). The latter scheme with step size δ guarantees that, for i=τˆk−1,…,τˆk−1,

where (ηi)τˆk−1<i≤τˆk is a sequence of independent random variables with standard Gaussian distribution. By eq. (2), it is convenient to consider the following equivalent parametrization: θk=(θ1,k,θ2,k,θ3,k), with θ1,k=(1−δλk), θ2,k=δλkμk and θ3,k=σkδ. The estimator θˆk of parameter θk is defined as:

which actually yields closed-form expressions:

(the sums in the above expressions range over [τˆk−1,τˆk−1]).

Under mild conditions, and if the process reaches the stationary regime, then (τˆ1,τˆ2) consistently estimates (τ1,τ2) (see Lavielle [13], Lavielle and Ludeña [14], for instance). Furthermore, if the true change points (τ1,τ2) are known then, under another set of mild conditions, the estimators (θˆ1,θˆ2,θˆ3) consistently estimate (θ1,θ2,θ3) (see Kessler [12], Theorem 1). To the best of our knowledge, there is no satisfactory result in the literature regarding the joint estimation of the change points (τ1,τ2) and the parameter (θ1,θ2,θ3).

3.2 Application to the real dataset

We undertake a simulation study of the properties of the two-step estimating procedure presented in the previous section and summarize its results in the next section. Because we characterize our simulation scheme based on the results obtained when applying the latter procedure to the real dataset, we first present them here.

For each subject and protocol, we estimate (τ1,τ2,θ1,θ2,θ3) from the corresponding (real) observed trajectory. The results pertaining to the estimation of (τ1,τ2) are summarized in Table 2 (the mean and standard deviation of the estimates (τˆ1δ,τˆ2δ) computed over each group of subjects are provided) and illustrated in Figure 1. Results pertaining to the estimation of (θ1,θ2,θ3) are summarized in Table 3 (the mean and standard deviation of the estimates (θˆ1,θˆ2,θˆ3) computed over each group of subjects are provided).

Three features of Table 2 are worth commenting on. First, one notes that there is no significant difference across groups of subjects (however, the means of τˆ1δ and τˆ2δ are slightly shifted to the left in the group of hemiplegic subjects relative to the two other groups). Second, the means of τˆ1δ are close to the time of start of perturbations for all groups and protocols. As for the means of τˆ2δ, they are close to the end time in vestibular and normal subjects and both protocols. In hemiplegic subjects, one notes that the standard deviations are quite large and that the mean of τˆ2δ is shifted to the left relative to the end time for perturbations. This is due to the fact that for each protocol, 30% of hemiplegic subjects feature an estimator τˆ2δ close to 30 s.

Third, judging by the standard deviations, all hemiplegic subjects tend to react similarly by adjusting quickly to the perturbations undergone in protocol 2. Likewise, all normal subjects tend to react similarly by adjusting quickly to the end of perturbations undergone in protocol 2. In protocol 1, the large standard deviation associated to the mean value of τˆ2δ computed over the group of normal subjects reflects the fact that 20% of these subjects feature an estimator τˆ2δ close to 30 s.

Regarding Table 3, we consider in turn parameters θ1,k, θ2,k and θ3,k (k=1,2,3). First, the estimates θˆ1,k behave similarly across groups of subjects, protocols and intervals [Tˆk−1,Tˆk) for each k=1,2,3. In contrast, for each k=1,2,3, the estimates of θˆ2,k behave quite differently across groups of subjects, protocols and intervals [Tˆk−1,Tˆk). This is a promising feature for the sake of classifying subjects by group, which is the problem at stake here. As for the estimates θˆ3,k (k=1,2,3), for any given group of subjects and protocol, they behave quite similarly. However, for protocol 2, it seems that the estimates θˆ3,k (k=1,2,3) in normal subjects behave differently from their counterpart in hemiplegic or vestibular subjects. This is another promising feature.

3.3 Simulation study

We carry out a simulation study to evaluate the performances of the two-step estimating procedure presented in Section 3.1. We directly simulate the trajectory C1:N. Specifically,

1. we set (T1,T2)=(τ1δ,τ2δ)=(15,50);

2. we rely on the Euler scheme (2) with step size δ/10 to approximate the sampling of (C(t))t∈[δ,Nδ] from the distribution characterized by eq. (1) where, for each k∈{1,2,3}, θk equals the mean of its 22 estimates based on the 22 real trajectories associated to the 22 hemiplegic subjects and protocol 1 (see Section 3.2 and Table 3);

3. we conclude by sub-sampling (C(t))t∈[δ,Nδ] to derive C1:N.

Three comments on Table 4 are in order. First, regarding the estimation of (τ1,τ2), we note that the means of the estimated change points are very close to their respective true values. Moreover, the standard deviations are small. Second, regarding the estimation of (θ1,k,θ3,k) (k=1,2,3), we note that the means of the estimates of θ1,k and θ3,k are very close to their respective true values and that the standard deviations are small. Third, regarding the estimation of θ2,k (k=1,2,3), we emphasize that the means are quite apart from their respective true values. Moreover, the standard deviations are not small. Overall this indicates a poorer estimation of θ2,k (k=1,2,3) than of (θ1,k,θ3,k) (k=1,2,3). This is probably due to the fact that the time intervals [Tk−1,Tk) (k=1,2,3) are relatively narrow given the value of δ.

4.1 Summary measures

Most classification procedures based on trajectories involve a preliminary step of dimension reduction where the high-dimensional trajectories are summarized into a low-dimensional summary measure [4, 5]. Here we build a tailored finite-dimensional summary measure of every trajectory based on the estimates (τˆ1,τˆ2,θˆ1,θˆ2,θˆ3) derived from it via our two-step estimating procedure.

We argue in Section 3.2 that, overall, τˆ1, τˆ2, θˆ2,k and θˆ3,k (k=1,2,3) may be relevant for the sake of classifying subjects as hemiplegic, vestibular or normal. Because the ratio θˆ2,k/(1−θˆ1,k)=μˆk is easier to interpret than θˆ2,k (and since θˆ1,k varies very little across subjects and protocols) and because θˆ3,k=σˆkδ, we choose to define our finite-dimensional summary measure as (τˆ1,τˆ2,μˆ1,μˆ2,μˆ3,σˆ1,σˆ2,σˆ3). Hereafter, we denote X1 the latter vector derived from the trajectory associated to protocol 1 and X2 that derived from the trajectory associated to protocol 2.

4.2 Classification procedure

We actually construct two classifiers, ϕˆ1 and ϕˆ2, to classify subjects as hemiplegic, vestibular or normal based on either X1 or X2 for ϕˆ1, and on both (X1,X2) for ϕˆ2.

The generic observed data structure associated to a given subject is (X1,X2,Y), where Y∈{1,2,3} indicates the subject’s group (with convention Y=1 for hemiplegic, Y=2 for vestibular, and Y=3 for normal). We denote P0 the true distribution of (X1,X2,Y).

Let S be the set of all classifiers based on X=(X1,X2). The misclassification risk associated to S∈S is R(P0)(S)=EP01{S(X)≠Y}. Denote by R∗=minS∈SR(P0)(S) its minimum. It is achieved at the Bayes classifier S∗∈S, characterized by S∗(X)=argmaxy∈{1,2,3}P0Y=y|X. For j=1,2, we also introduce the Bayes classifier S∗j∈S which relies only on Xj: S∗j(X)=argmaxy∈{1,2,3}P0(Y=y|Xj).

We consider Dn={(Xi,Yi),i=1,…,n}, where Xi=(Xi1,Xi2) is the summary measure associated to the subject i and Yi indicates his group. Our objective is to construct, based on Dn, ϕˆ1 and ϕˆ2 which respectively estimate the better classifier among S1∗ and S2∗(i.e., argminS1∗,S2∗R(P0)(Sj∗)) and S∗. We choose to rely on the popular methodology of random forests [15], which have proved powerful in a variety of applications. This is made easy thanks to the

R

randomForest

5.1 Performance of the classification procedure

We evaluate the performances of the classification procedure by the leave-one-out rule and the 0.632+bootstrap method as described by Efron and Tibshirani [16]. We apply the 0.632 + bootstrap method in order to take into account the fact that the leave-one-out rule may result in overly optimistic error rates. Resorting to the leave-one-out rule is notably motivated by the relatively small sample size of our dataset. The results are reported in Table 5 (second row).

With its leave-one-out performance equal to 74% and its 0.632 + bootstrap error rate equal to 0.35, the best classifier is ϕˆ1, which involves one protocol only. For curiosity, we also evaluate the performances of ϕˆ1 and ϕˆ2 when systematically replacing (τˆ1δ,τˆ2δ) with (15,50) (the start and end times of the perturbation phase). We report the results in Table 5 (first row). Quite satisfactorily, we note that both ϕˆ1 and ϕˆ2 are not as good as the original versions: estimating the change points proves very relevant.

5.2 Performance of an extended classification procedure

It is tempting to extend our classification procedure by simply extending the definition of the summary measure which is at its core. Following Chambaz and Denis [1], we merely substitute (X1,U1) and (X2,U2) for X1 and X2 with Uj (j=1,2) derived from C1:N as Uj=(Cˉ1+−Cˉ1−,Cˉ2−−Cˉ1+,Cˉ2+−Cˉ2−) where

We refer to Chambaz and Denis [1] for a justification. Then we apply the extended classification procedure and report its performances in Table 5 (third row).

The performances achieve by both classifiers ϕˆ1 and ϕˆ2 are slightly better than the best performance obtained in Section 5.1: enriching the summary measures seems to provide relevant additional information for the sake of classifying subjects as hemiplegic, vestibular or normal.

5.3 Perspectives

We address the statistical challenge of classifying subjects as hemiplegic, vestibular or normal based on complex trajectories obtained through two experimental protocols which were designed to evaluate potential deficits in postural control. The classification procedure involves a dimension reduction step where the complex trajectories are summarized by finite-dimensional summary measures based on a stochastic process model for a real-valued trajectory. This allows us to retrieve from the trajectories information relative to their temporal dynamic. A leave-one-out evaluation yields a 79% performance of correct classification for a total of n=70 subjects, with 22 hemiplegic (31%), 16 vestibular (23%) and 32 normal (46%) subjects.

In future work, we will extend the classification procedure by introducing finite-dimensional summary measures based on a stochastic process model for the original trajectories in ℝ4. We will also draw advantage from our good understanding of the classification problem to tackle the closely related statistical challenge of clustering subjects according to postural style.