Clustering of Directions Improves Goodness of Fit in Kinematic Data Collected in the Transverse Plane During Robot-Assisted Rehabilitation of Stroke Patients

The kinematic character of hand trajectory in reaching tasks varies by movement direction. Often, direction is not included as a factor in the analysis of data collected during multi-directional reach tasks; consequently, this directionally insensitive model (DI) may be prone to type-II error due to unexplained variance. On the other hand, directionally specific models (DS) that account separately for each movement direction, may reduce statistical power by increasing the amount of data groupings. We propose a clustered-by-similarity (CS) in which movement directions with similar kinematic features are grouped together, maximizing model fit by decreasing unexplained variance while also decreasing uninformative sub-groupings. We tested model quality in measuring change over time in 10 kinematic features extracted from 72 chronic stroke patients participating in the VA-ROBOTICS trial, performing a targeted reaching task over 16 movement directions (8 targets, back- and forth from center) in the horizontal plane. Across 49 participants surviving a quality control sieve, 4.3 ± 1.1 (min: 3; max: 7) clusters were found among the 16 movement directions; clusters varied between participants. Among 49 participants, and averaged across 10 features, the better-fitting model for predicting change in features was found to be CS assessed by the Akaike Information criterion (61.6 ± 7.3%), versus DS (31.0 ± 7.8%) and DI (7.1 ± 7.1%). Confirmatory analysis via Extra Sum of Squares F-test showed the DS and CS models out-performed the DI model in head-to-head (pairwise) comparison in >85% of all specimens. Thus, we find overwhelming evidence that it is necessary to adjust for direction in the models of multi-directional movements, and that clustering kinematic data by feature similarly may yield the optimal configuration for this co-variate.

Introduction

There is overwhelming evidence that the kinematics of voluntary movements vary depending on the direction of movement (Georgopoulos et al., 1988; Mushiake et al., 1991; Schwartz, 1993; Hewitt et al., 2011; Turnham et al., 2012). And yet, paradoxically, there is minimal precedence for accounting for this variability when analyzing data collected from multi-directional movement tasks. Here, we propose that including a factor term to account for movement direction can improve the goodness of model fit to multi-directional data sets; we test this paradigm in the setting of a ubiquitous line of inquiry, i.e., change in performance over time.

The cells comprising the primary motor cortex activate selectively in association with trajectory variables: each cell’s tuning centers on a preferred movement direction (Mushiake et al., 1991; Scott et al., 1997). This tuning varies from cell to cell (Georgopoulos et al., 1982, 1988), and the total population can co-vary with the hand trajectory (Schwartz, 1992, 1993; Scott et al., 1997). Many studies incorporating a multi-directional movement task via planar robots present data which evidence a substantial difference in the character of trajectories by direction (Yamamoto et al., 2006; Richardson et al., 2008; Brown et al., 2010; Hewitt et al., 2011; Turnham et al., 2012; Berniker et al., 2014; Muceli et al., 2014), an effect which persists in brain-injured patients (Shadmehr et al., 1998; Patton and Mussa-Ivaldi, 2004) and especially for stroke patients (Yarosh et al., 2004; Coderre et al., 2010; Scott and Dukelow, 2011). This observation is congruous with the notion that each person has preferred sectors of space through which arm movement is more proficient (Georgopoulos et al., 1986; Caminiti et al., 1990).

Given the thoroughness with which the neuromotor system has been shown to be directionally sensitive, established through studies of both cellular physiology as well as motor behavior, there is great need to develop robust strategies for adjusting for movement direction in analysis of kinematic data. The preponderance of studies make no mention of inclusion of any factor for direction whatsoever (Shadmehr et al., 1998; Rohrer et al., 2002; Daly et al., 2005; Finley et al., 2005; Hogan et al., 2006; Scheidt and Stoeckmann, 2007; Bosecker et al., 2010; Ganesh et al., 2014; Muceli et al., 2014; Sakamoto and Kondo, 2015); thus, we regard this approach as the “directionally insensitive” (DI) approach, where data from all movement directions are included in a single dataset with no directional co-variate term. Following the presumptive kinematic heterogeneity across multiple movement directions, we speculate that analytic approaches using the DI are more prone to error because of their substantial unexplained variance. A few studies take into account the different movement directions, typically by including a factor term for each direction (Scott et al., 1997; Yarosh et al., 2004; Maschke et al., 2004). While “directionally specific” (DS) models reduce the unexplained variance, it increases the number of groups in an ANOVA test. Such approaches require the estimation of more parameters than necessary, thereby reducing the amount of data available for hypothesis testing.

We conjecture that for a given patient, there may be some movement directions for which the kinematics are sufficiently similar so as to be combined with each other with minimal dilution. In this “clustered by similarity” (CS) model, we anticipate that by merging directions based on kinematic homogeneity, the unexplained variance will remain low, but the number of factors will decrease, thus increasing cell size, providing an intermediate alternative between the DI and the DS approaches. In this paper, we propose a method for merging movement directions based on features extracted from the kinematic record, and we compare the three approaches (DI, DS, and CS) by model fit.

Methods

Overview

In this study we examine improvement in the analysis of multi-directional kinematic data when the directions are included as factors in the analysis. Specifically, we test fit of the DI model versus two alternative approaches: (1) inclusion of a factor term for each distinct movement direction (DS), and (2) grouping of movement directions according to similarity among a set of standard kinematical features (CS). We perform this inquiry over kinematic data obtained from a cohort of stroke patients rehearsing targeted pointing movements with their affected limb.

Study Population and Protocol

Data were obtained from participants in the VA-ROBOTICS study, a multi-center randomized and controlled clinical trial conducted at four VA medical centers between November 2006 and October 2008 (Lo et al., 2010). In the VA-ROBOTICS trial, patients with moderate-to-severe upper-limb impairment, at least 6 months removed from index stroke (i.e., the most recent stroke event, in case of recurrent stroke events), were randomly assigned to one of three treatment groups: Usual Care (UC), Robot-assisted Therapy (RT), or Intensive Comparison Therapy (ICT). Participants assigned to RT and ICT completed a comprehensive training regimen comprising 36 sessions over 12 weeks, driven by an adaptive robot (RT) or trained therapist (ICT); those assigned to UC did not receive prescribed training through the study, but were allowed to pursue collateral training independently from the study. Participants receiving ICT or UC were allowed to engage in robot-assisted rehabilitation after the completion of their follow-up; kinematic data were collected from the robot regardless of treatment assignment. Kinematic data were recorded from all participants interacting with the robot (Table 1).

TABLE 1

Table 1. Baseline Characteristics of Study Participants.

The robot-assisted therapy was delivered via the In-Motion rehabilitation robotic system (Interactive Motion Technologies) (Aisen et al., 1997), which is specifically designed for clinical neurological applications. It delivers goal-directed, assisted upper-extremity movement, with an interactive computer-generated video program providing visual feedback to the patient. Participants engaged in a point-to-point movement task, comprising a series of up to 5 unassisted clockwise rotations through a circular target; there were 8 directional compass targets (N, NE, E, SE, S, SW, W, and NW) with both outward (from center to periphery, “toward”: “t”), and return (from periphery to center, “back”: “b”), yielding 16 total movement directions. North and South are oriented in the position of maximum and minimum outward reach, respectively, i.e., the positions farthest- and nearest position relative to the trunk. For an individual interfacing with the robot with their right-hand, East and West reflect the maximally lateral- and medial hand position, respectively (Figure 1, Top).

FIGURE 1

Figure 1. Visualization participant interacting with robot (A) and demonstration of target task (B). Red dot shows current cursor location, dashed red line illustrates a trajectory toward the target location (yellow dot). Trajectory is shown in figure for clarity, but was not shown to participants as part of their visual feedback: only the targets and the cursor were shown.

The robot pre-positioned the participant’s arm at the center target, and participants were expected to move toward the target in a single, smooth motion. When the participant was unable to complete a movement, the attending clinician could override the task by pressing the space bar, progressing to the next sequential target. Participants were instructed: “Your goal is to reach toward each of the red targets. If you are able to reach the respective targets, the robot will prompt you to move toward the next one. You will complete five cycles around the circle in a clockwise fashion. In the event you are unable to reach the target, I will pause the device and move your arm passively to the next start position. Tell me when you’re ready.”

Participants sat upright in a comfortable chair, and interacted with the robot through a manipulandum (Figure 1, Bottom). Pointing movements were made in the transverse plane, with a monitor at eye level displaying a cursor against a target field. Data consisted of instantaneous velocity in two-dimensions (left-right and fore-aft), sampled at 200 Hz. We smoothed these data bi-directionally with a 1st order low-pass Butterworth’s filter with 20 Hz cutoff; filter was re-applied after each differentiation.

The VA-ROBOTICS study was approved by the Institutional Review Board at each medical center and all study participants provided written informed consent.

Data Conditioning

In order to avoid error associated with the robot’s automated data processing routine, serial files collected from the robot were concatenated and repartitioned according to movement reversals, yielding complete movement cycles reflecting a single point-to-point movement from start to finish (Beed et al., 2017). Movement cycles that were unusually short in duration were presumed to be measurement artefact due to clinician override because of participant inability to complete the movement. Accordingly, we censored movement cycles with duration 20 time samples (0.1 s) or with variance in velocity = 0, corresponding with a spuriously brief or stationary data sample, respectively. All data conditioning activities were performed on the comprehensive dataset, i.e., all participants, all days at once.

Feature Extraction

In order to adequately capture the kinematic character of these data, we extracted 10 features that have been shown to be labile to movement direction (de Rugy et al., 2013; Hamel-Paquet, 2006; Hewitt et al., 2011; Maschke, 2003; Muceli et al., 2014). In addition to basic movement parameters of amplitude, duration, average- and peak velocity, average- and peak acceleration, we calculated four descriptors related to movement smoothness including the number of peaks in the velocity trace (Fetters and Todd, 1987; Michaelsen et al., 2006; Thelen et al., 1993), the tent metric i.e., the ratio of the area under the convex hull fitted to the velocity trace versus the area under the velocity trace itself (Rohrer et al., 2002), and two formulations of jerk (Hogan and Sternad, 2009; Rohrer et al., 2002). Where jerk-based metrics have been shown to be especially sensitive to basic movement parameters (Wininger et al., 2009, 2012), we tested one formulation witha widely analyzed behavior ${\displaystyle {\mathrm{j}}_1=\frac{1}{v_{max}\cdot \left(t_2-t_1\right)}{\int }_{t_1}^{t_2}\left|\stackrel{⃛}{x}\right|dt}$, where t₁ and t₂ are the start-time and cessation time of the movement, v_max is the maximum velocity, and $\stackrel{⃛}{x}$ is the third time derivative of position; j₁ has dimensions T⁻² (Rohrer et al., 2002), and one recently-proposed measure designed to be robust to these parameters ${\displaystyle {\mathrm{j}}_2=\frac{{\left(t_2-t_1\right)}^5}{A^2}{\int }_{t_1}^{t_2}{\stackrel{⃛}{x}}^{2}dt}$, where A is the total movement distance; j₂ is non-dimensionalized (Hogan and Sternad, 2009).

Features were tested for normality via Shapiro-Wilk test; features found to be non-normal were log-transformed. Within each feature, outliers were identified as observations more than 3 standard deviations from the mean across all participants; movement cycles yielding one or more outlier value were discarded.

Our main objective in this study was to establish patterns of similarity or difference between movement cycles, based on their performance, i.e., characteristics of their trajectory. Because these features are expected to be highly collinear (Wininger et al., 2009, 2012), feature set dimensionality was reduced via Principal Components Analysis (PCA) over the entire dataset, i.e., all participants all days. A minimum threshold was established a priori at 80%: the minimum number of PCs accounting for greater than 80% of the variance would be retained for analysis.

Cluster Identification

In order to test for directions with kinematic homology, participant-level data were assessed via a greedy pair-wise optimization of the Mahalanobis Distance

${\displaystyle {\mathrm{D}}_{\mathrm{i}\mathrm{j}}={\left({\overrightarrow{\mathrm{x}}}_{\mathrm{i}}-{\overrightarrow{\mathrm{x}}}_{\mathrm{j}}\right)}^{\mathrm{T}}{\left(\frac{{\mathrm{S}}_{\mathrm{i}}}{{\mathrm{n}}_{\mathrm{i}}}+\frac{{\mathrm{S}}_{\mathrm{j}}}{{\mathrm{n}}_{\mathrm{j}}}\right)}^{-1}\left({\overrightarrow{\mathrm{x}}}_{\mathrm{i}}-{\overrightarrow{\mathrm{x}}}_{\mathrm{j}}\right)}$ (1)

where x_i and x_j are the means of the PC-transformed features of direction i and j in k dimensions, where k = the number of re-combined features surviving PCA threshold; n_i and n_j are the sample sizes (number of movements cycles in each direction), and S_i and S_j the covariance matrices of the PC-transformed features for each direction. We obtain the p-value of each D_ij due to the fact that D_ij is approximately Chi-square distributed with degree of freedom k, under the null hypothesis that the PC-transformed features in direction i are sampled from the same multi-variate Gaussian distribution as the PC-transformed features from direction j, equal to the number of dimensions mentioned above. The null hypothesis is: D_ij is not significantly different from zero, indicating the two directions are similar to each other based on Mahalanobis Distance and should be merged together. In the first iteration, this distance was computed for all 1 ≤ i ≤ 15 and i + 1 < j ≤ 16 (120 total comparisons); if the two directions with the minimum D_ij yielded a non-significant difference (p_max(Dij) > 0.05), the directions were merged into a new cluster (yielding 15 clusters total), and the process re-iterated for all 1 ≤ i ≤ 14 and i + 1 < j ≤ 15 (105 total remaining comparisons), and so-on. For each round of mergers, the Mahalanobis threshold was adjusted for multiple comparisons via the Bonferroni correction. This process continued until no pair-wise D yielded significance.

The Mahalanobis Distance requires S_i and S_j to be invertible, meaning that covariance matrix should be full-rank. Thus, the number of observations for each task must be equal to or greater than the number of PCs. To ensure viability, it was decided a priori to merge data over all days; participants with insufficient data were removed from the analysis stream.

Statistical Models

For each participant, three models were tested:

Directionally Insensitive model (DI)

${\displaystyle fi={\beta }_0+{\beta }_1\ast Day}$ (2) 

Directionally Specific model (DS)

${\displaystyle fi={\beta }_0+{\beta }_1\ast Day+\sum _{j=2}^{16}{\beta }_{2j}\ast Di{r}_j+\sum _{j=2}^{16}{\beta }_{3j}\ast \left(Day\ast Di{r}_j\right)}$ (3)

Clustered by Similarity model (CS)

${\displaystyle fi={\beta }_0+{\beta }_1\ast Day+\sum _{j=2}^J{\beta }_{2j}\ast Clus{t}_j+\sum _{j=2}^J{\beta }_{3j}\ast \left(Day\ast Clus{t}_j\right)}$ (4)

where f_i is the i^th feature, 1 ≤ i ≤ 10, and J is the total number of clusters for the participant; β₀ is the intercept and β_1,2,3 are the regression slope coefficients. Day is a factor variable related to the date of the visit (for assessing change across multiple visits), Dir indicates one of 16 directions, and Clust indicates one of J clusters. The DI represents the prevalent statistical model for showing change in a variable over a period of observation. The DS is designed to increase the proportion of explained variance of the DI by accounting for the 16 movement directions, but at the cost of reduced sample size in the model comparisons. The CS putatively decreases the complexity of the DS by consolidating directions where the Mahalanobis Distance is minimized, thus decreasing the number of comparisons thereby increasing comparison sample size. We hypothesized that both the DS and CS will improve model fit versus the DI. These equations are intended to variously extend the models presented in the study’s previous outcomes assessments (Lo et al., 2010; Wu et al., 2016). We are unaware of these specific models, particularly DS or CS, having been used elsewhere previously.

Model Quality

Firstly, a hierarchy of model quality was established via Akaike Information Criterion: AIC = 2 k – 2 ln [L], where L is the maximized likelihood, and k is the number of terms in the model. There is no rigorous test of model optimality for comparing heterogenous multi-level models, but the AIC is a useful heuristic (Akaike, 1981; Sin and White, 1996). For each participant, for each feature, we measured the frequency of AIC minimization by the three models.

Additionally, for each of the three possible outcomes (AIC minimized through DI, DS, or CS), the difference between the AIC values for the best-performing model was calculated for the two remaining models. Through this, we seek to measure the relative improvement of the optimum model over the two alternatives.

We observe that the DI can be considered nested within the CS, which in turn can be considered nested within the DS. Thus, model quality can be tested in these two pair-wise comparisons via the Extra Sum of Squares F test

${\displaystyle F=\frac{\frac{SS{E}_{Model1}-SS{E}_{Model2}}{df(SS{E}_{Model1})-df(SS{E}_{Model2})}}{\frac{SS{E}_{Model2}}{df(SS{E}_{Model2})}}}$(5)

where this test statistic has a F[df(SSE_Model1)−df(SSE_Model2), SSE_Model2] distribution, and Model1 corresponds to the model nested within Model2. Null hypothesis is the simpler model (Model1) fits adequately on the data (p-value > 0.05). For each participant, for each feature, we measured the frequency of significant test results in three comparisons: DI versus CS, DI versus DS, and CS versus DS.

Results

Descriptive Statistics

Data were collected from all participants engaging with the robot. In total 13,977 movement cycles were captured over 202 participant days. All movement cycles were found to be of adequate data quality to support analysis; 996 movement cycles were removed as outliers and 19 participants (1,134 movement cycles) were discarded due to a lack of a second recorded session with the robot. Among the surviving movement cycles, we observe that their basic kinematic parameters appear to commensurate with a moderately-impaired cohort (Table 2).

TABLE 2

Table 2. Basic kinematic parameters; F = Feature.

We note that in a two-way ANOVA, movement direction and the interaction term between movement direction and participant ID yielded significant p-values for all 10 features (results not shown). This supports the assumption of a sensitivity of kinematic performance to movement direction.

Principal Components

Most participants (35 of 49) required two PCs to reach the 80% threshold; the remainder (14 of 49) needed three. Subsequently ten participants (677 movement cycles) yielded datasets with one or more movement directions containing insufficient data to execute clustering, and were removed from the analysis stream. The final dataset comprised 11,119 movement cycles performed by 49 participants captured over 154 participant days. A correlation matrix for these features is shown in Table 3.

TABLE 3

Table 3. Correlation matrix of extracted features; F = Feature.

Clustering

Across the remaining participants, the data clustered into 4.3 ± 1.1 groupings (min: 3; max: 7). Among the single-most populated cluster, the number of movement directions contained therein was variable with a mean ± SD deviation of 7.5 ± 2.2 directions (min: 4; max: 14). Figure 2 shows two exemplars of participant data with clustering results.

FIGURE 2

Figure 2. Kinematic data collected from two stroke patients reflecting movement of the robot manipulandum in transverse plane. Data shown ensemble as recorded from robot (Left Panels), as well as rotated and grouped by cluster (Right Panels).

In Figure 2, the star-shaped portraits show kinematic data super-imposed in the plane of movement in toward movements (top) and back (bottom) directions. The columns showing break-out boxes with horizontal traces are the kinematic data, grouped by cluster, and rotated to left-to-right orientation for visualization. Participant #49 exhibits severely impaired movement; Participant #57 shows proficient movement. While space limitations preclude an extended analysis of the differences and similarities of the clusters across participants, we report that there is no one movement direction consistently grouped within the largest cluster (Table 4).

TABLE 4

Table 4. Representation of Directions among participants' largest clusters.

Model Quality

Table 5 presents the frequency of model optimization and the margin of model improvement for each feature at the p = 0.05 threshold, and for all features combined for all three thresholds.

TABLE 5

Table 5. Model Quality Results: Frequency of model superiority, and margin of improvement; F = Feature.

In particular, we note that the Directionally Insensitive model is best only 7% of the time (average across all features), and when the DI is best, it provides only a modest improvement over the clustered model, decreasing AIC by only 7.0 points versus CS (IQR 3.2–40.4 points decrease). The DS was the optimum model 31% of the time, although we note that in those cases where DS was the best model, its impact on AIC minimization was substantial: 62.8 points reduction in AIC versus CS (IQR: 36.7–94.3 points decrease). The CS yielded the superior model 62% of the time, with moderate improvement over the other two models: −33.1 (−61.0–−17.1) versus DI and −20.2 (−61.0–−11.0) versus DS.

Comparison via the F-test showed similar results (Table 6). The DS was consistently superior to the DI: range across 10 features, the F-test yielded p < 0.05 in at least 75% of samples, and as much as 96% (average: 86.7 ± 7.3%). Similarly, the CS proved superior to DI: range 78 to 100% (average: 91.2 ± 7.9%). We note that the margin of improvement of using DS versus CS was moderate: 36 to 67% (average: 49.6 ± 7.5%).

TABLE 6

Table 6. F test results between different model comparisons for each feature; within table: F = Feature.

Discussion

Study Validity

In this study, we test whether inclusion of a factor for movement direction increases model fit in the context of a simplistic, but common hypothesis test, i.e., change in a kinematical feature over time. Testing this method with data collected from an impaired cohort over an extended time period (stroke patients from the VA-ROBOTICS study, with 217 ± 172 days between first- and final measurements) is strategic in that there is a greater likelihood that some participants will show non-trivial change in these features over time, providing greater opportunity to capture an association between response and predictor.

While we considered assessing kinematic similarity based on the raw trajectory, we decided in the end to predicate our clustering on features extracted from the trajectory. While recognize that some will consider clustering-by-trajectory as a more direct approach, we believe clustering within the feature space to be a reasonable surrogate for trajectory, and caution that comparison of trajectories is a complex enterprise: patients with neuromotor dysfunction can move slowly or rapidly, and sometimes both slowly and rapidly in the same data sample. In moderately-impaired populations, movement arrest can be particularly prevalent (Beppu et al., 1984; Rohrer et al., 2002), and can confound attempts to obtain meaningful waveform correlations (Wininger et al., 2009).

And while not an explicit objective of this study, we note that this study partly replicates the findings of others, published in this journal who demonstrated that there is substantial efficiency to be gained –with minimal loss in accuracy– in reducing model complexity in kinematic analysis of human upper-limb movement in robotic systems (Averta et al., 2017).

Interpretations

The paradigm of adding a new factor to an analysis changes both statistical power, and model interpretation. Given the vastness of existing studies where planar robots are used to capture kinematics in multi-dimensional movement, it is impractical to attempt to address the implications of a DI approach with specificity. Broadly speaking, one prominent impact of the new CS approach is the opportunity to re-analyze data presented in previous studies: non-significant results may become significant due to the identification of homogeneity across multiple movement directions, and significant results could benefit from increased effect size. Equivalently, confirming negative findings through the incorporation of CS would add value.

There is inherent complexity in utilizing two assays (minimization of AIC versus frequency of significant F-tests). We assert the value in this pluralistic approach: two distinct analyses provide mutual validation in the absence of a gold-standard approach. At the same time, we recognize that the results, while largely consistent, diverge somewhat in the assessment of DS versus CS. Per AIC, CS is the better model 61% of the time (Table 5), per F-test, CS fails to yield significantly better results in more than half the cases (50.4%, Table 6). We view these as compatible results: the discrepancy is small (approximately 10% absolute difference), and is likely attributable to the disparate nature of the measures: AIC minimization is arithmetically based, where F-test involves transformation and thresholding. Taken together, we believe that CS does provide non-trivial improvement versus DS, but additional work may be helpful in illuminating the true magnitude of benefit.

Application

There is progressively greater evidence that rehabilitation robots are most efficacious when employing a cooperative control strategy versus passive support (Israel et al., 2006; Ziherl et al., 2010; Balasubramanian et al., 2012; Casadio and Sanguineti, 2012; Van der Loos et al., 2016). And it is generally preferred that a rehabilitation robot offer assistance in a way that reflects the level of impairment (Novak and Riener, 2015), which can vary across the workspace: movements in some directions are more proficient than others. Physical therapists are experts at identifying areas of weakness, and prioritizing those areas for focused rehabilitation (Sahrmann, 1988). Rehabilitation robots provide obvious advantages: precision placement and support, durability, and analytical sophistication. However, because the performance measures available for quantifying motor skill are many and varied, and collinearities are common (Table 3), there remain substantial opportunities to streamlining robot measurement and control strategies so that they are efficient and more human-like. We note that the interest in bridging the gap between rehabilitation robots and their human counterparts extends to identification of kinematical factors underlying the clinical functional scores (Bosecker et al., 2010; Scott and Dukelow, 2011; van Dokkum et al., 2014). We anticipate that our approach may provide new inroads towards understanding the subtle, intuitive approach of the human therapist in assessing motor skill.

Limitations

In terms of study design, this work can be considered generalizable in the sense that we analyzed data from patients in active rehabilitation programs, as well as those with no ongoing training. This study is limited by participant demography: we included primarily male patients, with mild- to moderate impairment due to stroke, an older population with substantial co-morbid burden. Furthermore, space limitations necessarily narrow the scope of this paper to purely technical matters; there is limited opportunity for meaningful clinical or physiological inquiry.

Methodologically, this approach faces some limitations. Primarily is that due to sample size. In order to prevent singularity, there must be at least N + 1 samples in each movement direction, where N is the dimensionality of the feature space. We accommodated this restriction by reducing feature space dimensionality through principal components analysis. PCs reflect a weighted sum of the variables, accounting for a progressively smaller amount of variability among the observed data. Increasing to a higher dimensionality, while preferable for clustering robustness, would have meant the removal of more participants’ data. We note that while many participants had ample data in most movement directions, many participants had difficulty in just one or two movement directions; this approach requires adequate samples in all movement directions. As a consequence of the use of principal components, the clustering is performed in a 2-D or 3-D space described by a linear re-combination of features, as opposed to the full 10-D space created by the raw features. On the other hand, our results (first two PCs yielded >80% of the variance) suggests the sufficiency of a low-dimensional transformed feature space, and the additional sample size gained through use of PCA evidences its value. Moreover, PCA provides excellent protection against feature collinearity: there is a reduced burden in feature design and selection when the features are re-combined in a way that maximizes their combined proportional variance.

Analytically, we note that this study is limited by the lack of clear tools for concise measurement of model goodness of fit. In order to convincingly test the three models, we analyzed in two ways: comparison via AIC, and via the Extra Sum of Squares F-test. We anticipate that in many cases, both the AIC and the F-test will be supported, but we acknowledge that the AIC is not universally considered a robust approach for model comparison, and the F-test does not easily avail to comparison beyond the binary assessment of above- or below the pre-defined significance threshold.

The Mahalanbois distance was chosen as the preferred approach for clustering out of expeditiousness. So-called greedy methods such as serial re-grouping based on minimum-distance criterion applied to a small number of candidate models, are not comprehensive, and therefore we cannot be certain that the results reported here reflect the true utility of the CS. A fully complete test of the CS would require assessment of model fit in all possible regression models. However, this poses an intractably large calculation: for the 1-cluster regression, there is a single candidate model (this is the DI); for the 2-cluster regression, there are ${\displaystyle \sum _{i=0}^{15}{C}_i^{15}}$ = 32,768 ways to divide 16 movement directions into two clusters; for the 3-cluster regression, there are ${\displaystyle C_2^{16}\cdot \sum _{i=0}^{14}\left(C_i^{14}\cdot \left(\sum _{j=0}^{14-i}{C}_j^{14-i}\right)\right)}$ = 573,956,280 ways to divide 16 movement directions into 3 clusters, and so on. A complete search across all candidate models would require many billions of calculations per participant-feature; given our sample size (n = 72 participants x 10 features), a full search is simply not feasible. Nevertheless, while our greedy model cannot be interpreted as exhaustive, the CS model showed consistent superiority versus the DI and DS, and could only be enhanced by more extensive search. Summarily, we assert that our incorporation of the Mahalanobis Distance as a heuristic for clustering directional data –while not exhaustive– was efficient and effective, and yielded results which strongly evidence the superiority of the CS.

Extension and Future Work

The present study is the first, known to the authors, to directly test the impact of movement direction as a co-variate in the kinematical analysis. Furthermore, ours is the first study to identify the tradeoff between explained variance and cell size in the two extant approaches (DI versus DS), and to propose an alternative, i.e., CS. As a result, this study has generated new knowledge, not only for methodological rigor, but for the study of human motor behavior: we find that the 16 movement directions cluster naturally into approximately 4 groups of kinematic similarity in mild- to moderately impaired chronic stroke patients (p < 0.05 threshold).

In the interest of brevity and clarity we defer focused analysis of the trajectory clustering. However, our preliminary inquiry reveals an intriguing cluster patterning. Movement directions were noted for number of times they clustered with each of the 15 other movement directions (Table 7). We highlight those cells with a large number of shared clusters (n ≥ 20), and note the highly diagonal orientation of these co-clusters within the matrix. From this, we identify the following approximate clusters in the all-participants view: {NWt, Nt, NEt}, {NWb, Nb, NEb}, {SEb, Sb, SWb}. In particular, we observe: (1) none of the backwards directions are associated with towards directions, (2) both towards and back directions have distinct groupings of North and South.

TABLE 7

Table 7. Co-clustered movement directions across participants.

Replication will be required in order to measure the clustering in other populations, including normative data on healthy controls. Additional exploration will give insight into whether these clusters are consistent across patients, whether the movement directions in some clusters are more sensitive to training than others, and whether a cluster-based approach will provide badly needed leverage in the ongoing research into treatment response (Lo et al., 2010), and the association between kinematics and clinical-based performance measures (Bosecker et al., 2010).

Disclaimer

Opinions herein are those of the individual authors and the contents do not represent the views of the Department of Veterans Affairs of the United States Government.

Ethics Statement

The study protocol was approved by the institutional review boards at each participating site and the human rights committee at the coordinating center, listed as follows: VA Maryland Health Care System (Baltimore), North Florida/South Georgia Veterans Health System (Gainesville), VA Puget Sound Health Care System (Seattle), and VA Connecticut Healthcare System (West Haven). All participants gave written informed consent prior to study participation. The study was registered on ClinicalTrials.gov (ClinicalTrials.gov identifier, NCT 00372411).

Author Contributions

This secondary study was initially conceived and conducted by LL, JH, XW, and MW. Original trial design and operations was overseen by PG and PP. Data conditioning and code review was performed by AB and MW. All authors were involved in the drafting and approving of this manuscript.

Funding

This study was conducted by the United States (U.S.) Department of Veterans Affairs, Office of Research and Development, Cooperative Studies Program. Additional support for this analysis was through the Veterans Affairs Connecticut Research and Education Foundation.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Aisen, M. L., Krebs, H. I., Hogan, N., McDowell, F., and Volpe, B. T. (1997). The effect of robot-assisted therapy and rehabilitative training on motor recovery following stroke. Arch. Neurol. 54 (4), 443–446. doi: 10.1001/archneur.1997.00550160075019

PubMed Abstract | CrossRef Full Text | Google Scholar

Akaike, H. (1981). Likelihood of a model and information criteria. J. Econom. 16 (1), 3–14. doi: 10.1016/0304-4076(81)90071-3

CrossRef Full Text | Google Scholar

Averta, G., Della Santina, C., Battaglia, E., Felici, F., Bianchi, M., and Bicchi, A. (2017). Unvealing the principal modes of human upper limb movements through functional analysis. Front. Robot. AI 4. doi: 10.3389/frobt.2017.00037

CrossRef Full Text | Google Scholar

Balasubramanian, S., Colombo, R., Sterpi, I., Sanguineti, V., and Burdet, E. (2012). Robotic assessment of upper limb motor function after stroke. J. Phys. Med. Rehabil. 91, S255–S269. doi: 10.1097/PHM.0b013e31826bcdc1

CrossRef Full Text | Google Scholar

Beed, A. T., Peduzzi, P., Guarino, P., and Wininger, M. (2017). A partitioning algorithm for extracting movement epochs from robot-derived kinematic data. Front. Robot. AI 4. doi: 10.3389/frobt.2017.00057

CrossRef Full Text | Google Scholar

Beppu, H., Suda, M., and Tanaka, R. (1984). Analysis of cerebellar motor disorders by visually guided elbow tracking movement. Brain J. Neurol 107 (3), 787–809. doi: 10.1093/brain/107.3.787

CrossRef Full Text | Google Scholar

Berniker, M., Franklin, D. W., Flanagan, J. R., Wolpert, D. M., and Kording, K. (2014). Motor learning of novel dynamics is not represented in a single global coordinate system: evaluation of mixed coordinate representations and local learning. J. Neurophysiol. 111 (6), 1165–1182. doi: 10.1152/jn.00493.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Bosecker, C., Dipietro, L., Volpe, B., Krebs, H. I., and Igo Krebs, H. (2010). Kinematic robot-based evaluation scales and clinical counterparts to measure upper limb motor performance in patients with chronic stroke. Neurorehabil. Neural Repair 24 (1), 62–69. doi: 10.1177/1545968309343214

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, L. E., Wilson, E. T., Obhi, S. S., and Gribble, P. L. (2010). Effect of trial order and error magnitude on motor learning by observing. J. Neurophysiol. 104 (3), 1409–1416. doi: 10.1152/jn.01047.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Caminiti, R., Johnson, P. B., and Urbano, A. (1990). Making arm movements within different parts of space: dynamic aspects in the primate motor cortex. J. Neurosci. 10 (7), 2039–2058. doi: 10.1523/JNEUROSCI.10-07-02039.1990

PubMed Abstract | CrossRef Full Text | Google Scholar

Casadio, M., and Sanguineti, V. (2012). Learning, retention, and slacking: a model of the dynamics of recovery in robot therapy. IEEE Trans. Neural Syst. Rehabil. Eng. 20 (3), 286–296. doi: 10.1109/TNSRE.2012.2190827

PubMed Abstract | CrossRef Full Text | Google Scholar

Coderre, A. M., Zeid, A. A., Dukelow, S. P., Demmer, M. J., Moore, K. D., Demers, M. J., et al. (2010). Assessment of upper-limb sensorimotor function of subacute stroke patients using visually guided reaching. Neurorehabil. Neural Repair 24 (6), 528–541. doi: 10.1177/1545968309356091

PubMed Abstract | CrossRef Full Text | Google Scholar

Daly, J. J., Hogan, N., Perepezko, E. M., Krebs, H. I., Rogers, J. M., Goyal, K. S., et al. (2005). Response to upper-limb robotics and functional neuromuscular stimulation following stroke. J. Rehabil. Res. Dev. 42 (6), 723–736.

PubMed Abstract | Google Scholar

de Rugy, A., Loeb, G. E., and Carroll, T. J. (2013). Are muscle synergies useful for neural control? Front. Comput. Neurosci. 7:19. doi: 10.3389/fncom.2013.00019

PubMed Abstract | CrossRef Full Text | Google Scholar

Fetters, L., and Todd, J. (1987). Quantitative assessment of infant reaching movements. J. Mot. Behav. 19 (2), 147–166. doi: 10.1080/00222895.1987.10735405

PubMed Abstract | CrossRef Full Text | Google Scholar

Finley, M. A., Fasoli, S. E., Dipietro, L., Ohlhoff, J., Macclellan, L., Meister, C., et al. (2005). Short-duration robotic therapy in stroke patients with severe upper-limb motor impairment. J. Rehabil. Res. Dev. 42 (5), 683–692. doi: 10.1682/JRRD.2004.12.0153

PubMed Abstract | CrossRef Full Text | Google Scholar

Ganesh, G., Yoshioka, T., Osu, R., and Ikegami, T. (2014). Immediate tool incorporation processes determine human motor planning with tools. Nat. Commun. 5:4524. doi: 10.1038/ncomms5524

PubMed Abstract | CrossRef Full Text | Google Scholar

Georgopoulos, A. P., Kalaska, J. F., Caminiti, R., and Massey, J. T. (1982). On the relations between the direction of two-dimensional arm movements and cell discharge in primate motor cortex. J. Neurosci. 2 (11), 1527–1537. doi: 10.1523/JNEUROSCI.02-11-01527.1982

PubMed Abstract | CrossRef Full Text | Google Scholar

Georgopoulos, A. P., Kettner, R. E., and Schwartz, A. B. (1988). Primate motor cortex and free arm movements to visual targets in three-dimensional space. II. Coding of the direction of movement by a neuronal population. J. Neurosci. 8 (8), 2928–2937. doi: 10.1523/JNEUROSCI.08-08-02928.1988

PubMed Abstract | CrossRef Full Text | Google Scholar

Georgopoulos, A. P., Schwartz, A. B., and Kettner, R. E. (1986). Neuronal population coding of movement direction. Science 233 (4771), 1416–1419. doi: 10.1126/science.3749885

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamel-Pâquet, C., Sergio, L. E., and Kalaska, J. F. (2006). Parietal area 5 activity does not reflect the differential time-course of motor output kinetics during arm-reaching and isometric-force tasks. J. Neurophysiol. 95 (6), 3353–3370. doi: 10.1152/jn.00789.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Hewitt, A. L., Popa, L. S., Pasalar, S., Hendrix, C. M., and Ebner, T. J. (2011). Representation of limb kinematics in Purkinje cell simple spike discharge is conserved across multiple tasks. J. Neurophysiol. 106 (5), 2232–2247. doi: 10.1152/jn.00886.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Hogan, N., Krebs, H. I., Rohrer, B., Palazzolo, J. J., Dipietro, L., Fasoli, S. E., et al. (2006). Motions or muscles? Some behavioral factors underlying robotic assistance of motor recovery. J. Rehabil. Res. Dev. 43 (5), 605–618. doi: 10.1682/JRRD.2005.06.0103

PubMed Abstract | CrossRef Full Text | Google Scholar

Hogan, N., and Sternad, D. (2009). Sensitivity of smoothness measures to movement duration, amplitude, and arrests. J. Mot. Behav. 41 (6), 529–534. doi: 10.3200/35-09-004-RC

PubMed Abstract | CrossRef Full Text | Google Scholar

Israel, J. F., Campbell, D. D., Kahn, J. H., and Hornby, T. G. (2006). Metabolic costs and muscle activity patterns during robotic- and therapist-assisted treadmill walking in individuals with incomplete spinal cord injury. Phys. Ther. 86 (11), 1466–1478. doi: 10.2522/ptj.20050266

PubMed Abstract | CrossRef Full Text | Google Scholar

Lo, A. C., Guarino, P. D., Richards, L. G., Haselkorn, J. K., Wittenberg, G. F., Federman, D. G., et al. (2010). Robot-assisted therapy for long-term upper-limb impairment after stroke. N. Engl. J. Med. 362 (19), 1772–1783. doi: 10.1056/NEJMoa0911341

PubMed Abstract | CrossRef Full Text | Google Scholar

Maschke, M., Gomez, C. M., Ebner, T. J., and Konczak, J. (2004). Hereditary cerebellar ataxia progressively impairs force adaptation during goal-directed arm movements. J. Neurophysiol. 91 (1), 230–238. doi: 10.1152/jn.00557.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Michaelsen, S. M., Dannenbaum, R., and Levin, M. F. (2006). Task-specific training with trunk restraint on arm recovery in stroke: randomized control trial. Stroke 37 (1), 186–192. doi: 10.1161/01.STR.0000196940.20446.c9

PubMed Abstract | CrossRef Full Text | Google Scholar

Muceli, S., Falla, D., and Farina, D. (2014). Reorganization of muscle synergies during multidirectional reaching in the horizontal plane with experimental muscle pain. J. Neurophysiol. 111 (8), 1615–1630. doi: 10.1152/jn.00147. 2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Mushiake, H., Inase, M., and Tanji, J. (1991). Neuronal activity in the primate premotor, supplementary, and precentral motor cortex during visually guided and internally determined sequential movements. J. Neurophysiol. 66 (3), 705–718. doi: 10.1152/jn.1991.66.3.705

PubMed Abstract | CrossRef Full Text | Google Scholar

Novak, D., and Riener, R. (2015). Control strategies and artificial intelligence in rehabilitation robotics. AIMag 36 (4), 23. doi: 10.1609/aimag.v36i4.2614

CrossRef Full Text | Google Scholar

Patton, J. L., and Mussa-Ivaldi, F. A. (2004). Robot-assisted adaptive training: custom force fields for teaching movement patterns. IEEE Trans. Biomed. Eng. 51 (4), 636–646. doi: 10.1109/TBME.2003.821035

PubMed Abstract | CrossRef Full Text | Google Scholar

Richardson, A. G., Lassi-Tucci, G., Padoa-Schioppa, C., and Bizzi, E. (2008). Neuronal activity in the cingulate motor areas during adaptation to a new dynamic environment. J. Neurophysiol. 99 (3), 1253–1266. doi: 10.1152/jn.01096. 2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Rohrer, B., Fasoli, S., Krebs, H. I., Hughes, R., Volpe, B., Frontera, W. R., et al. (2002). Movement smoothness changes during stroke recovery. J. Neurosci. 22 (18), 8297–8304. doi: 10.1523/JNEUROSCI.22-18-08297.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Sahrmann, S. A. (1988). Diagnosis by the physical therapist—a prerequisite for treatment. Phys. Ther. 68 (11), 1703–1706. doi: 10.1093/ptj/68.11.1703

CrossRef Full Text | Google Scholar

Sakamoto, T., and Kondo, T. (2015). Visuomotor learning by passive motor experience. Front. Hum. Neurosci. 9:279. doi: 10.3389/fnhum.2015. 00279

PubMed Abstract | CrossRef Full Text | Google Scholar

Scheidt, R. A., and Stoeckmann, T. (2007). Reach adaptation and final position control amid environmental uncertainty after stroke. J. Neurophysiol. 97 (4), 2824–2836. doi: 10.1152/jn.00870.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwartz, A. B. (1992). Motor cortical activity during drawing movements: single-unit activity during sinusoid tracing. J. Neurophysiol. 68 (2), 528–541. doi: 10.1152/jn.1992.68.2.528

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwartz, A. B. (1993). Motor cortical activity during drawing movements: population representation during sinusoid tracing. J. Neurophysiol. 70 (1), 28–36. doi: 10.1152/jn.1993.70.1.28

PubMed Abstract | CrossRef Full Text | Google Scholar

Scott, S. H., and Dukelow, S. P. (2011). Potential of robots as next-generation technology for clinical assessment of neurological disorders and upper-limb therapy. J. Rehabil. Res. Dev. 48 (4), 335. doi: 10.1682/JRRD.2010.04.0057

PubMed Abstract | CrossRef Full Text | Google Scholar

Scott, S. H., Sergio, L. E., and Kalaska, J. F. (1997). Reaching movements with similar hand paths but different arm orientations. II. Activity of individual cells in dorsal premotor cortex and parietal area 5. J. Neurophysiol. 78 (5), 2413–2426. doi: 10.1152/jn.1997.78.5.2413

PubMed Abstract | CrossRef Full Text | Google Scholar

Shadmehr, R., Brandt, J., and Corkin, S. (1998). Time-dependent motor memory processes in amnesic subjects. J. Neurophysiol. 80 (3), 1590–1597. doi: 10.1152/jn.1998.80.3.1590

PubMed Abstract | CrossRef Full Text | Google Scholar

Sin, C. -Y., and White, H. (1996). Information criteria for selecting possibly misspecified parametric models. J. Econom. 71 (1-2), 207–225. doi: 10.1016/0304-4076(94)01701-8

CrossRef Full Text | Google Scholar

Thelen, E., Corbetta, D., Kamm, K., Spencer, J. P., Schneider, K., and Zernicke, R. F. (1993). The transition to reaching: mapping intention and intrinsic dynamics. Child Dev. 64 (4), 1058–1098. doi: 10.2307/1131327

PubMed Abstract | CrossRef Full Text | Google Scholar

Turnham, E. J., Braun, D. A., and Wolpert, D. M. (2012). Facilitation of learning induced by both random and gradual visuomotor task variation. J. Neurophysiol. 107 (4), 1111–1122. doi: 10.1152/jn.00635.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

van Dokkum, L., Hauret, I., Mottet, D., Froger, J., Métrot, J., and Laffont, I. (2014). The contribution of kinematics in the assessment of upper limb motor recovery early after stroke. Neurorehabil. Neural Repair 28 (1), 4–12. doi: 10.1177/1545968313498514

PubMed Abstract | CrossRef Full Text | Google Scholar

Van der Loos, H. F. M., Reinkensmeyer, D. J., and Guglielmelli, E. (2016). “Rehabilitation and Health Care Robotics,” in Springer Handbook of Robotics, eds B. Siciliano, and O. Khatib (Cham, Switzerland: Springer International Publishing), 1685–1728.

Wininger, M., Kim, N. H., and Craelius, W. (2009). Spatial resolution of spontaneous accelerations in reaching tasks. J. Biomech. 42 (1), 29–34. doi: 10.1016/j.jbiomech.2008.10.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Wininger, M., Kim, N. H., and Craelius, W. (2012). Reformulation in the phase plane enhances smoothness rater accuracy in stroke. J. Mot. Behav. 44 (3), 149–159. doi: 10.1080/00222895.2012.663012

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, X., Guarino, P., Lo, A. C., Peduzzi, P., and Wininger, M. (2016). Long-term effectiveness of intensive therapy in chronic stroke. Neurorehabil. Neural Repair 30 (6), 583–590. doi: 10.1177/1545968315608448

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamamoto, K., Hoffman, D. S., and Strick, P. L. (2006). Rapid and long-lasting plasticity of input-output mapping. J. Neurophysiol. 96 (5), 2797–2801. doi: 10.1152/jn.00209.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Yarosh, C. A., Hoffman, D. S., and Strick, P. L. (2004). Deficits in movements of the wrist ipsilateral to a stroke in hemiparetic subjects. J. Neurophysiol. 92 (6), 3276–3285. doi: 10.1152/jn.00549.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Ziherl, J., Novak, D., Olenšek, A., Mihelj, M., and Munih, M. (2010). Evaluation of upper extremity robot-assistances in subacute and chronic stroke subjects. J. Neuroeng. Rehabil. 7:52. doi: 10.1186/1743-0003-7-52

PubMed Abstract | CrossRef Full Text | Google Scholar