Original article| Volume 93, ISSUE 9, P1541-1552, September 01, 2012

# Longitudinal Patterns of Functional Recovery in Patients With Incomplete Spinal Cord Injury Receiving Activity-Based Rehabilitation

## Abstract

Lorenz DJ, Datta S, Harkema SJ. Longitudinal patterns of functional recovery in patients with incomplete spinal cord injury receiving activity-based rehabilitation.

### Objective

To model the progression of 3 functional outcome measures from patients with incomplete spinal cord injury (SCI) receiving standardized locomotor training.

### Design

Observational cohort.

### Setting

The NeuroRecovery Network (NRN), a specialized network of treatment centers providing standardized, activity-based therapy for SCI patients.

### Participants

Patients (N=337) with incomplete SCI (grade C or D on the International Standards for Neurological Classification of Spinal Cord Injury scale) who were enrolled in the NRN between February 2008 and March 2011.

### Intervention

All enrolled patients received standardized locomotor training sessions, as established by NRN protocol, and were evaluated monthly for progress.

### Main Outcome Measures

Berg Balance Scale, 6-minute walk test, and 10-meter walk test. Progression over time was analyzed via the fitting of linear mixed effects models.

### Results

There was significant improvement on each outcome measure and significant attenuation of improvement over time. Patients varied significantly across groups defined by recovery status and American Spinal Injury Association Impairment Scale (AIS) grade at enrollment with respect to baseline performance and rates of change over time. Time since SCI was a significant determinant of the rate of recovery for all measures.

### Conclusions

Locomotor training, as implemented in the NRN, results in significant improvement in functional outcome measures as treatment sessions accumulate. Variability in patterns of recovery over time suggest that time since SCI and patient functional status at enrollment, as measured by the Neuromuscular Recovery Scale, are important predictors of performance and recovery as measured by the targeted outcome measures.

## Key Words

#### List of Abbreviations:

AIC (Akaike information criterion), AIS (American Spinal Injury Association Impairment Scale), CI (confidence interval), ISNCSCI (International Standards for Neurological Classification of Spinal Cord Injury), NRN (NeuroRecovery Network), NRS (Neuromuscular Recovery Scale), SCI (spinal cord injury)
ACTIVITY-BASED interventions are emerging as a more successful approach for functional recovery after neurologic injury.
• Behrman A.L.
• Bowden M.G.
• Nair P.M.
Neuroplasticity after spinal cord injury and training: an emerging paradigm shift in rehabilitation and walking recovery.
• Harkema S.J.
• Behrman A.L.
• Bratta A.
• Sisto S.A.
• Edgerton V.R.
Establishing the NeuroRecovery Network: multisite rehabilitation centers that provide activity-based therapies and assessments for neurologic disorders.
Studies of locomotor training have demonstrated that recovery of walking and balance function can occur for individuals months and years after incomplete spinal cord injury (SCI).
• Harkema S.J.
• Lorenz D.J.
• Edgerton V.R.
• Behrman A.L.
Balance and ambulation improvements in individuals with chronic incomplete spinal cord injury using locomotor training-based rehabilitation.
However, observed magnitudes and rates of improvement in balance and walking function can be highly variable. Intensity of treatment is a critical determinant of recovery for locomotor training, and its cost per treatment session is higher than that for therapies that only require a single physical therapist. Therefore, knowledge of factors that can influence rates and magnitudes of recovery would directly aid prediction of the potential benefit for an individual and help guide treatment decisions.
Many evaluations of therapeutic interventions focus strictly on the endpoints of a study (enrollment and discharge), a strategy that provides information only on how much individuals recovered and not how that recovery occurred over time. This can leave critical questions unanswered, such as whether patients continued to receive benefit from rehabilitative therapy at discharge. A longitudinal examination of progression through an intervention can provide valuable information, which can be used to better target therapy for future patients. For example, a longitudinal examination can identify when the average patient will receive maximum impact (recover most rapidly) from a given intervention, determine when plateaus in recovery may be reached, or identify subsets of patients receiving the greatest/least average benefit from the intervention. This information can be used clinically to determine when, during the course of therapy, to specifically focus treatment efforts, when to discharge patients, and which subgroups of patients may benefit most from a given mode of therapy.
Such an examination has been conducted for stroke patients,
• Tilling K.
• Sterne J.A.
• Rudd A.G.
• Glass T.A.
• Wityk R.J.
• Wolfe C.D.
A new method for predicting recovery after stroke.
in which investigators modeled the Barthel Index
• Mahoney F.I.
• Barthel D.
Functional evaluation: the Barthel Index.
over time, as a function of various patient characteristics through the fitting of multilevel models. Results demonstrated that several patient characteristics—presence of prestroke disability, urinary incontinence, dysarthria, and sex—were associated with lower Barthel Index scores after stroke. Further, their models demonstrated that patients over 80 years of age, with dysphasia, or with a limb deficit exhibited poorer functional recovery over time.
The goal of our analysis was to identify characteristics that influenced the rate and extent of recovery in individuals who received standardized locomotor training
• Dobkin B.
• Apple D.
• Barbeau H.
• et al.
Weight-supported treadmill vs over-ground training for walking after acute incomplete SCI.
• Visintin M.
• Barbeau H.
The effects of body weight support on the locomotor pattern of spastic paretic patients.
• Wernig A.
• Nanassy A.
• Müller S.
Maintenance of locomotor abilities following Laufband (treadmill) therapy in para- and tetraplegic persons: follow-up studies.
• Dietz V.
• Colombo G.
• Jensen L.
• Baumgartner L.
Locomotor capacity of spinal cord in paraplegic patients.
• Behrman A.L.
• Harkema S.J.
Physical rehabilitation as an agent for recovery after spinal cord injury.
• Harkema S.J.
Plasticity of interneuronal networks of the functionally isolated human spinal cord.
• Behrman A.K.
• Lawless-Dixon A.R.
• Davis S.B.
• et al.
Locomotor training progression and outcomes after incomplete spinal cord injury.
within the NeuroRecovery Network (NRN), a multicenter rehabilitative program.
• Harkema S.J.
• Behrman A.L.
• Bratta A.
• Sisto S.A.
• Edgerton V.R.
Establishing the NeuroRecovery Network: multisite rehabilitation centers that provide activity-based therapies and assessments for neurologic disorders.
We hypothesized that severity of injury, time since injury, and age of the individual would influence the rate and magnitude of improvement of standard walking and balance measures.

## Methods

### Participants

Data from 337 patients enrolled in the NRN were examined (table 1). Those eligible for enrollment in the NRN had incomplete nonprogressive spinal cord lesions at level T10 or above, were classified as American Spinal Injury Association Impairment Scale (AIS) grade C or D on the International Standards for Neurological Classification of Spinal Cord Injury (ISNCSCI) scale, were not participants in an inpatient rehabilitation program, and met additional previously reported NRN eligibility criteria.
• Harkema S.J.
• Behrman A.L.
• Bratta A.
• Sisto S.A.
• Edgerton V.R.
Establishing the NeuroRecovery Network: multisite rehabilitation centers that provide activity-based therapies and assessments for neurologic disorders.
The patients considered in this analysis were recruited from 7 NRN centers—Boston Medical Center, Boston, MA (28 patients); Frazier Rehab Institute, Louisville, KY (23 patients); Kessler Institute for Rehabilitation, West Orange, NJ (60 patients); Magee Rehabilitation Hospital, Philadelphia, PA (74 patients); The Ohio State University Medical Center, Columbus, OH (25 patients); Shepherd Center, Atlanta, GA (82 patients); and The Institute for Rehabilitation and Research, Houston, TX (45 patients) —and were enrolled in the NRN between February 1, 2008 and March 31, 2010. All individuals who enrolled in the NRN between these dates and completed at least 1 functional evaluation were considered in this analysis. The submission of demographic, clinical, and functional outcome data to a centralized NRN database, from which the data for this article were gathered, was approved by each center's institutional review board. Each patient enrolled in the NRN provided a signed informed consent form prior to data collection.
Table 1Characteristics
CharacteristicsN=337
Sex
Female82 (24)
Male255 (76)
Age (y)40±17
C99 (29)
D238 (71)
Injury level
Cervical249 (74)
Thoracic88 (26)
Mechanism of injury
MVC107 (32)
Fall75 (22)
Sporting accident65 (19)
Medical/surgical29 (9)
Nontrauma26 (8)
Violence24 (7)
Other11 (3)
Assistive walking device
Nonambulatory137 (41)
Walker98 (29)
Cane(s)/crutch(es)63 (18)
None40 (12)
Time since SCI (y)1.0 (0.1, 52.3)
Phase at enrollment
Phase 1—1A, 1B, 1C26 (8), 59 (18), 69 (20)
Phase 2—2A, 2B, 2C59 (18), 37 (11), 30 (9)
Phase 3—3A, 3B, 3C29 (9), 20 (6), 8 (2)
Treatment and enrollment characteristics
Time of NRN enrollment (mo)
Significantly differed among the 3 numeric NRS phases (Kruskal-Wallis test, P<.001).
Thirty-two patients were enrolled in the NRN just before the cutoff point for this analysis (March 2011) or exited the NRN after their initial evaluation, and subsequently had zero enrollment time, received zero treatment sessions, and had 1 functional evaluation.
3.2 (0, 52.5)
Significantly differed among the 3 numeric NRS phases (Kruskal-Wallis test, P<.001).
Thirty-two patients were enrolled in the NRN just before the cutoff point for this analysis (March 2011) or exited the NRN after their initial evaluation, and subsequently had zero enrollment time, received zero treatment sessions, and had 1 functional evaluation.
40 (0, 353)
Cumulative no. of evaluations
Significantly differed among the 3 numeric NRS phases (Kruskal-Wallis test, P<.001).
Thirty-two patients were enrolled in the NRN just before the cutoff point for this analysis (March 2011) or exited the NRN after their initial evaluation, and subsequently had zero enrollment time, received zero treatment sessions, and had 1 functional evaluation.
3 (1, 18)
Treatment intensity (Tx/evaluation)18±5
NOTE. Values are mean ± SD, median (minimum, maximum), or counts (%).
Abbreviations: MVC, motor vehicle collision; Tx, treatments.
Significantly differed among the 3 numeric NRS phases (Kruskal-Wallis test, P<.001).
Thirty-two patients were enrolled in the NRN just before the cutoff point for this analysis (March 2011) or exited the NRN after their initial evaluation, and subsequently had zero enrollment time, received zero treatment sessions, and had 1 functional evaluation.

### Intervention

All patients enrolled in the NRN, and in particular those under consideration here, received standardized locomotor training and were periodically evaluated for progress on health, quality of life, and functional outcome measures. Functional evaluations, which were scheduled for every 20 treatment sessions, occurred on average every 18±5 sessions. The protocol for NRN locomotor training sessions and the standardization of the functional evaluations are detailed in Harkema et al
• Harkema S.J.
• Behrman A.L.
• Bratta A.
• Sisto S.A.
• Edgerton V.R.
Establishing the NeuroRecovery Network: multisite rehabilitation centers that provide activity-based therapies and assessments for neurologic disorders.
in this issue.

### Data Analysis

Linear mixed effects models were used to provide longitudinal models of recovery for Berg Balance Scale
• Berg K.O.
• Wood-Dauphinee S.L.
• Williams J.I.
• Maki B.
Measuring balance in the elderly: validation of an instrument.
• Berg K.O.
• Maki B.E.
• Williams J.I.
• Holliday P.J.
• Wood-Dauphinee S.L.
Clinical and laboratory measures of postural balance in an elderly population.
scores, 6-minute walk test distances, and 10-meter walk test speeds.
• van Hedel H.J.
• Wirz M.
• Dietz V.
Assessing walking ability in subjects with spinal cord injury: validity and reliability of 3 walking tests.
For longitudinal data, the mixed effects model predicts average performance on an outcome over time based on a set of covariates, the fixed effects, while allowing for patient-to-patient variation in the shape of the recovery curves through the random effects. Other than time, parameterized as the cumulative number of NRN treatment sessions received, the fixed effects included in our model were age, time since SCI, patient phase of recovery as defined by the Neuromuscular Recovery Scale (NRS),
• Behrman A.L.
• Ardolino E.
• VanHiel L.R.
• et al.
Assessment of functional improvement without compensation reduces variability of outcome measures after human spinal cord injury.
and AIS grade (C or D) on the ISNCSCI scale
American Spinal Injury Association
Reference manual for the International Standards for Neurological Classification of Spinal Cord Injury.
• Marino R.J.
• Barros T.
• Biering-Sorensen F.
• et al.
International standards for neurological classification of spinal cord injury.
; each was measured at enrollment. Note that we employed the alphanumeric specification of NRS phase (1A, 1B, 1C, 2A, etc) rather than the numeric NRS phase (1, 2, 3) considered in the source work for the NRS contained in this issue.
• Behrman A.L.
• Ardolino E.
• VanHiel L.R.
• et al.
Assessment of functional improvement without compensation reduces variability of outcome measures after human spinal cord injury.
Models of the 6-minute walk and 10-meter walk tests included an additional fixed effect denoting assistive device use. This term was parameterized as a 2-level factor, 1 level for the device the patient used at enrollment, termed the initial device, and 1 level for the patient's current use device, termed the current device. At every NRN evaluation, walk tests are assessed up to 2 times, once with the initial device and once with the current device (when available).
To avoid model overspecification and to circumvent natural covariation among several of the predictors (see Results section), we placed restrictions on the structure of the fixed effects, and in particular the interactions between the fixed effects permitted in the model. These restrictions are further discussed in the Model Specification subsection below.
The primary purpose of our analysis was to identify the best fitting longitudinal recovery model based on our set of covariates. To achieve this, we followed a previously detailed
• Verbecke G.
• Molenberghs G.
Linear mixed models for longitudinal data.
procedure in which a saturated or overparameterized fixed effects structure was fit and parameters that varied with time were selected as random effects. Homogeneity of variance was graphically assessed through diagnostic plotting after the fitting of this initial model. Variance heterogeneity identified in this diagnostic step was modeled through the introduction of residual variance functions, as detailed in mixed effects model literature.
• Pinheiro J.
• Bates D.
Mixed effects models in S and S-Plus.
Fixed effects in the final model were tested using sequential F tests from analyses of variance in a hierarchical fashion, with main effects being tested first, followed by polynomial terms, interactions, and polynomial-interaction terms.
A secondary goal of our analysis was to compare AIS grade and NRS phase as predictors of recovery. To accomplish this, we fit 3 different models for each outcome: (1) an NRS-only model, in which AIS was excluded as a predictor; (2) an AIS-only model, in which NRS was excluded as a predictor; and (3) an NRS + AIS model, in which both were included. The fixed effects specification for each of these models is detailed in the Model Specification section below. The NRS-only and AIS-only models were nested within the NRS + AIS model (because their fixed effects were subsets of the fixed effects in the NRS + AIS model), and comparisons of these 2 models against the more general NRS + AIS model were conducted through likelihood ratio tests from analyses of variance.
• Pinheiro J.
• Bates D.
Mixed effects models in S and S-Plus.
The NRS-only and AIS-only models, being nonnested, were compared informally using the Akaike information criterion (AIC), under which models with lower AIC are favored.
• Akaike H.
A new look at the statistical model identification.
In the Results section below, we provide comparisons of the NRS-only, AIS-only, and NRS + AIS models to determine which provided the best fit for each outcome measure. After identifying the best fitting models for each outcome, each model is further examined by considering the results of significance tests of model terms through sequential analyses of variance. The fixed effects estimates from the significant terms in the model are then detailed, and their practical interpretation by the plotting of average recovery curves for a variety of NRN patients defined by differing covariate values (NRS phases, time since SCI, etc) are discussed. Consideration of technical model details, such as model equations, covariance for the random effects, and residual variance functions is given in appendix 1. Demographic and clinical characteristics of our sample were summarized descriptively with means, SDs, medians, and extrema for continuous data and counts and percentages for categorical data. All analyses were conducted using the open-source R software package,
R Development Core Team
R: a language and environment for statistical computing.
and in particular the nlme package,
• Pinheiro J.
• Bates D.
• Debroy S.
• Sarkar D.
The R Core Team
nlme: Linear and nonlinear mixed effects models.
offered as an R extension. All significance tests were conducted at the .05 level.

### Model Specification

Based on inspection of individual level recovery curves (not shown), an attenuated recovery model, in which patients improve from enrollment at a steadily declining rate, was selected as most appropriate. Thus, the basic fixed effects specification included an intercept, linear rate term, and quadratic rate term, reflective of performance at enrollment, the rate of improvement, and the rate of attenuation of improvement, respectively. The additional predictors of recovery were age, time since SCI, NRS phase at enrollment, and AIS grade at enrollment. As noted in the Data Analysis section above, restrictions were placed on the fixed effects specification to avoid overfitting and predictor covariation. This included the following: (1) all predictors were included as main effects; (2) time since SCI and age were not permitted to interact; (3) NRS and AIS were not permitted to interact in the NRS + AIS model; and (4) all 2- and 3-way interaction terms were constrained to involve either time (the linear rate parameter) or time2 (the quadratic rate parameter).
These restrictions effectively defined a 9-parameter quadratic recovery curve for each phase of the NRS for the NRS-only model, both levels of AIS for the AIS-only model, and each level of both for the NRS + AIS model. The 9 parameters for each NRS phase or AIS level were (1) the intercept, modifications to the intercept by (2) time since SCI and (3) age, (4) the linear rate, modifications to the linear rate by (5) time since SCI and (6) age, (7) the quadratic rate, and modifications to the quadratic rate by (8) time since SCI and (9) age. Hence, age and time since SCI were permitted to independently impact the 3 parameters (intercept, linear rate, quadratic rate) of an attenuated recovery curve, and this impact was permitted to vary the NRS phase and/or AIS grade. We briefly note that 32 patients received only 1 assessment at enrollment, with no follow-up. While these patients do not contribute information about functional recovery, their at-enrollment assessments do inform the mixed effects model estimates of the model intercepts. Finally, the study site was included as a 7-level factor in each model in order to control for intersite variation in the outcome measures.

## Results

### Demographic and Clinical Characteristics

The 337 patients in our sample exhibited demographic characteristics corresponding to other NRN datasets
• Harkema S.J.
• Lorenz D.J.
• Edgerton V.R.
• Behrman A.L.
Balance and ambulation improvements in individuals with chronic incomplete spinal cord injury using locomotor training-based rehabilitation.
• Behrman A.L.
• Ardolino E.
• VanHiel L.R.
• et al.
Assessment of functional improvement without compensation reduces variability of outcome measures after human spinal cord injury.
and to the incomplete SCI population at large (see table 1). Our sample was composed largely of patients with AIS grade D injuries (71%) at the cervical level (74%), and represented a diverse set of injury mechanisms. Forty-one percent of our patients were nonambulatory at enrollment, and most were classified in NRS phase 1.
• Behrman A.L.
• Ardolino E.
• VanHiel L.R.
• et al.
Assessment of functional improvement without compensation reduces variability of outcome measures after human spinal cord injury.
Patients were well-distributed across the 3 NRS phases of recovery at enrollment. Time since SCI varied from just over 1 month to 52.3 years, with a median of just under 1 year. Treatment and enrollment characteristics in the NRN were highly variable, with enrollment time and cumulative number of locomotor training sessions received ranging as high as 52.5 months and 353 sessions, respectively. Patients were evaluated a median of 3 times, and this ranged from 1 to 18 evaluations. NRS phase 1 patients remained enrolled in the NRN longer and received more locomotor training sessions and evaluations than phase 2 and 3 patients (Kruskal-Wallis tests, P<.001), indicative of patients with lower function at enrollment requiring more treatment and corresponding with previous studies on NRN patients.
• Harkema S.J.
• Lorenz D.J.
• Edgerton V.R.
• Behrman A.L.
Balance and ambulation improvements in individuals with chronic incomplete spinal cord injury using locomotor training-based rehabilitation.
• Behrman A.L.
• Ardolino E.
• VanHiel L.R.
• et al.
Assessment of functional improvement without compensation reduces variability of outcome measures after human spinal cord injury.
These treatment characteristics did not significantly differ by AIS grade (P>.35).
NRS phase and AIS grade at enrollment were strongly associated (Fisher exact test, P<.001), because no AIS grade C patient was classified higher than phase 2B. Age was significantly associated with time since SCI (Spearman ρ=.13, P=.01), and AIS grade D patients were significantly older than AIS grade C patients (mean, 41 vs 36, P=.007, Wilcoxon rank-sum test). Time since SCI significantly differed over NRS phases, with NRS phase 1 patients having been enrolled farthest from injury and NRS phase 3 patients being enrolled closest to injury (Kruskal-Wallis test, P<.001). AIS grade C patients tended to enroll later after injury than AIS grade D patients (P=.004). These associations among predictors informed the specification of the fixed effects of our recovery models, which are discussed in the next subsection. We note that time since SCI was heavily skewed in our sample, ranging from 0.1 to 52.3 years with a mean of 2.9 and a median of 0.9. Hence, time since SCI was log-transformed when included in our recovery models.
We note that for each outcome measure, the study site was a significant model term (P<.001). Rates of recovery and attenuation, however, did not significantly vary by study site for any outcome measure (P>.09). Therefore, only the main effect for study site was included in the model for each outcome. The study site term was included only to adjust for intersite variation, and the fixed effects estimates for this term were not of interest and are not presented here. In what follows, fixed effects estimates of model intercepts are necessarily presented as averages over the 7 sites.

### Comparison of Phase and AIS Models

The NRS-only model for the Berg Balance Scale was superior to the AIS-only model (AIC=9352.3 vs 9616.2) in predicting recovery. The NRS + AIS model provided a significantly better fit to our data than both the NRS-only (χ211=42.5, P<.001) and AIS-only (χ272=432.4, P<.001) models. Thus, while the NRS-only model provided a better fit to the Berg Balance Scale data than the AIS-only model, the addition of AIS to the NRS-only model significantly contributed to the model fit, and the NRS + AIS model was selected for further consideration of the fixed effects.
In addition, for the 6-minute walk test, the NRS-only model was superior to the AIS-only model (AIC=–1867 vs –1575), and the NRS + AIS model was significantly better than the NRS-only (χ29=22.3, P=.01) and AIS-only (χ273=613.3, P<.001) models. For the 10-meter walk test, the NRS-only model provided a better fit than the AIS-only model (AIC=–1223 vs –988), and the NRS + AIS was significantly better than both (NRS: χ210=105.6, P<.001; AIS: χ273=466.8, P<.001). Thus, the NRS + AIS models were selected for further consideration for both walk tests.

### Berg Balance Scale Model

Berg Balance Scale scores significantly improved from enrollment at an attenuated rate (table 2, Berg column). Performance at enrollment varied by NRS phase and AIS grade, but did not vary with time since SCI or age. NRS phase, AIS grade, and time since SCI impacted the rates of improvement and attenuation, as judged by their significant interactions with the linear (time) and quadratic (time2) rates. Age was not a significant determinant of rates of improvement or attenuation. The impact of time since SCI on the rate of improvement varied by NRS phase and AIS grade, as judged by 3-way interactions involving these terms. The impact of time since SCI on the attenuation rate also varied by NRS phase. The 3-way interaction of age, NRS phase, and the linear rate was significant but not of particular interest, because the lower order, 2-way interaction of age and the linear rate was nonsignificant.
Table 2P Values From Tests of the Fixed Effects Included in the Mixed Effects Model of the Berg Balance Scale, 6-Minute Walk Test, and 10-Meter Walk Test (intercept term not included)
Fixed EffectThis Term Tests Whether . . .Berg6MWT10MWT
TimeRate of improvement was significant<.001
Significant.
<.001
Significant.
< .001
Significant.
Time2Rate of attenuation was significant<.001
Significant.
<.001
Significant.
< .001
Significant.
NRS phaseEnrollment performance varied by NRS phase<.001
Significant.
<.001
Significant.
< .001
Significant.
DeviceEnrollment performance varied by deviceNA.160.020
Significant.
Time since SCIEnrollment performance impacted by time since SCI.240.690.190
AISEnrollment performance varied by AIS<.001
Significant.
.004
Significant.
.010
Significant.
AgeEnrollment performance impacted by age.500.020
Significant.
.020
Significant.
NRS phase×timeRate of improvement varied by NRS phase<.001
Significant.
<.001
Significant.
<.001
Significant.
Time since SCI×timeRate of improvement modified by time since SCI<.001
Significant.
<.001
Significant.
<.001
Significant.
Age×timeRate of improvement modified by age.570.140.110
AIS×timeRate of improvement varied by AIS.004
Significant.
.150.010
Significant.
Device×timeRate of improvement varied by deviceNA.540.680
NRS phase×time2Rate of attenuation varied by NRS phase<.001
Significant.
<.001
Significant.
<.001
Significant.
Time since SCI×time2Rate of attenuation modified by time since SCI<.001
Significant.
<.001
Significant.
<.001
Significant.
Age×time2Rate of attenuation modified by age.060.500.270
AIS×time2Rate of attenuation varied by AIS.001
Significant.
.110.090
Device×time2Rate of attenuation varied by deviceNA.540.110
NRS phase×time since SCI×timeImpact of time since SCI on rate of improvement varied by NRS phase<.001
Significant.
<.001
Significant.
<.001
Significant.
NRS phase×age×timeImpact of age on rate of improvement varied by NRS phase.001
Significant.
.040
Significant.
.420
AIS×time since SCI×timeImpact of time since SCI on rate of improvement varied by AIS.008
Significant.
.180.160
AIS×age×timeImpact of age on rate of improvement varied by AIS.430.850.920
Device×time since SCI×timeImpact of time since SCI on rate of improvement varied by deviceNA.910.990
Device×age×timeImpact of age on rate of improvement varied by deviceNA.690.260
NRS phase×time since SCI×time2Impact of time since SCI on rate of attenuation varied by NRS phase<.001
Significant.
<.001
Significant.
<.001
Significant.
NRS phase×age×time2Impact of age on rate of attenuation varied by NRS phase.220.020
Significant.
.260
AIS×time since SCI×time2Impact of time since SCI on rate of attenuation varied by AIS.850.790.430
AIS×age×time2Impact of age on rate of attenuation varied by AIS.070.880.830
Device×time since SCI×time2Impact of time since SCI on rate of attenuation varied by deviceNA.040
Significant.
.700
Device×age×time2Impact of age on rate of attenuation varied by deviceNA.010
Significant.
.100
NOTE. P values calculated from sequential F tests of model terms. Fixed effects are grouped according to type—main effects (top section), 2-way interactions, 2-way polynomial interactions, 3-way interactions, and 3-way polynomial interactions (bottom section).
Abbreviations: Berg, Berg Balance Scale; NA, no assistive device used for Berg Balance Scale assessments, and therefore a device term was not necessary; 6MWT, 6-minute walk test; 10MWT, 10-meter walk test.
Significant.
The estimates of significant model fixed effects (table 3) and plots of average recovery curves (fig 1) detail the characteristics of the Berg Balance Scale model. NRS phases significantly differed and were well-ordered with respect to performance at enrollment, with NRS phase 1A patients exhibiting the lowest average scores (1.1) and NRS phase 3C patients exhibiting the highest average scores (46.0).
Table 3Coefficient Estimates of Significant Fixed Effects With 95% CIs From the Linear Mixed Effects Model of the NRN Berg Balance Scale Data
NRS Phase at EnrollmentEnrollment Performance (intercept)Recovery Rate (linear term)Attenuation Rate (quadratic term)Time SCI- Recovery Rate InteractionTime SCI- Attenuation Rate Interaction
1A1.1 (–2.6 to 4.9)0.2 (–1.5 to 1.9)0 (–0.2 to 0.2)–0.7 (–2.0 to 0.6)0.1 (–0.1 to 0.2)
1B3.5 (0.8 to 6.1)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
1.8 (0.6 to 2.9)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.1 (–0.2 to 0)–0.7 (–1.6 to 0.1)0.1 (0 to 0.1)
1C6.0 (3.3 to 8.7)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
1.9 (0.5 to 3.3)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.1 (–0.2 to 0.1)–1 (–2 to 0)0.1 (–0.1 to 0.2)
2A15.4 (11.9 to 18.9)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
4.5 (2.7 to 6.2)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.3 (–0.5 to –0.1)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–3.6 (–4.7 to –2.5)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.3 (0.1 to 0.4)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
2B24.8 (20.6 to 28.9)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
3.8 (1.8 to 5.9)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.3 (–0.6 to –0.1)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–5.0 (–6.5 to –3.5)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.5 (0.2 to 0.8)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
2C30.9 (26.6 to 35.2)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
4.4 (2.2 to 6.7)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.5 (–0.8 to –0.2)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–2.9 (–4.3 to –1.5)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.3 (0.1 to 0.5)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
3A38.3 (33.8 to 42.8)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
4.2 (1.5 to 6.9)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–1.1 (–1.7 to –0.5)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–4.5 (–6.5 to –2.6)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
1.0 (0.3 to 1.7)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
3B40.2 (35.0 to 45.3)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.5 (–4.4 to 5.4)0.3 (–1.7 to 2.3)–4.3 (–8.3 to –0.3)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.8 (–1.1 to 2.8)
3C46.0 (38.4 to 53.6)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–2.2 (–20.8 to 16.4)0.9 (–16.7 to 18.4)–4 (–20.8 to 12.9)3.2 (–15.6 to 22)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
3.2 (1.8 to 4.5)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.2 (–0.4 to –0.1)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.6 (0.0 to 1.2)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
NS
NOTE. All rate parameters are given as rates per 20 sessions, approximately corresponding to rates per NRN evaluation. Time since SCI interaction effects are given in log-years. Intercepts are averages over study sites. The last row of the table provides estimates of differences between AIS grade C and AIS grade D patients for the given parameters.
Abbreviation: NS, not a significant model term.
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
Significant recovery was seen for patients in NRS phases 1B–3A, with patients in NRS phases 2A–3A improving most rapidly. Significant attenuation of recovery was observed for NRS Phases 2A–3A (see fig 1A). This attenuation increased with NRS phase (NRS phase 3A has the largest attenuation), a byproduct of the scale boundary (max score=56) and the fact that higher NRS phases performed better at enrollment and were in closer proximity to this boundary. Time since SCI had a significant negative impact on recovery rates for patients in NRS phases 2A–3B, as noted by the negative coefficient estimates for the time since SCI-linear rate interaction. Time since SCI also significantly positively impacted the attenuation rate for NRS phase 2A–3A patients. The positive coefficients indicate that recovery curves for patients farther removed from their injuries were less severely attenuated. Thus, patients farther removed from injury recovered balance less rapidly, but had flatter recovery curves because of lower attenuation (see fig 1B). We also note that there was no significant longitudinal progress for patients in NRS phases 1A, 3B, and 3C, because no rate parameters were marginally significantly different from 0. NRS phase 1A patients on average did not respond and NRS phase 3B and 3C patients possessed very high balance function at enrollment with little room for improvement (average scores=40.2 and 46.0). Injury severity (AIS grade) also significantly impacted performance at enrollment and recovery over time (see fig 1C). Patients with AIS grade D injuries tended to perform better at enrollment and exhibit a greater rate of improvement that was more significantly attenuated than AIS grade C patients. Further, the adverse impact of time since SCI on the recovery rate was significantly less for AIS grade D patients. We note that these conclusions regarding the difference between AIS grades C and D patients apply only to those NRS phases that AIS grades C and D patients occupied, 1A through 2B.

### Six-Minute Walk Test Model

Distances for the 6-minute walk test significantly improved from enrollment at an attenuated rate (see table 2, 6MW column). Performance at enrollment varied with NRS phase, AIS grade, and age, but not by assistive device type or time since SCI. NRS phase and time since SCI significantly impacted the rates of improvement and attenuation, as judged by their significant interactions with the linear (time) and quadratic (time2) rates. Age, device type, and AIS grade were not significant determinants of rates of improvement or attenuation. The impact of time since SCI on the rates of improvement and attenuation significantly varied by NRS phase, as judged by 3-way interactions involving these terms. Three-way interactions between NRS phase, age, and recovery rate, NRS phase, age and attenuation rate, and device, age, and attenuation rate were significant, but lower order, 2-way interactions involving age, device, and recovery and attenuation rate were nonsignificant, and therefore these 3-way interactions were not of particular interest.
Estimates of significant model fixed effects (table 4) and plots of average recovery curves (fig 2) for the 6-minute walk test distances showed that NRS phases significantly differed and were well-ordered with respect to performance at enrollment. NRS phase 1A patients walked the shortest average distance (13m) and NRS phase 3C patients walked the greatest average distance (310m). Patients with AIS grade D injuries walked on average 32m farther (95% confidence interval [CI], 10–55m) than AIS C grade patients at enrollment. Age significantly impacted enrollment performance, as every 1-year increase in age at enrollment corresponded with a 0.6 m (95% CI, –1.1 to –0.1m) decrease in enrollment distance.
Table 4Coefficient Estimates of Significant Fixed Effects With 95% CIs From Linear Mixed Effects Model of NRN 6-Minute Walk Test Data
NRS Phase at EnrollmentEnrollment Performance (intercept)Recovery Rate (linear)Attenuation Rate (quadratic)Time SCI- Recovery Rate InteractionTime SCI- Attenuation Rate Interaction
1A13 (–22 to 48)1 (–13 to 16)0 (–1 to 1)0 (–11 to 11)0 (–1 to 1)
1B20 (–9 to 50)14 (4 to 24)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–1 (–2 to 0)0 (–7 to 7)0 (–1 to 1)
1C29 (1 to 57)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
13 (4 to 23)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–1 (–1 to 0)–8 (–16 to 0)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0 (0 to 1)
2A73 (41 to 106)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
42 (31 to 53)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–3 (–4 to –2)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–23 (–31 to –15)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
2 (1 to 3)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
2B133 (95 to 171)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
61 (46 to 76)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–6 (–7 to –4)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–42 (–54 to –30)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
5 (2 to 7)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
2C142 (104 to 179)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
68 (52 to 83)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–6 (–8 to –4)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–42 (–53 to –31)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
4 (2 to 5)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
3A222 (182 to 263)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
99 (73 to 125)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–19 (–26 to –12)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–71 (–92 to –51)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
16 (9 to 24)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
3B267 (224 to 311)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
63 (29 to 97)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–7 (–19 to 5)–39 (–67 to –10)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
12 (–1 to 24)
3C310 (237 to 383)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
29 (–142 to 199)24 (–95 to 142)–194 (–335 to –52)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
169 (38 to 299)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
NOTE. All rate parameters are given as rates per 20 sessions, approximately corresponding to rates per NRN evaluation. Time since SCI interaction effects are given in log-years. Intercepts are averages over study sites.
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
Significant recovery was seen for patients in NRS phases 1B–3B, and the rate of recovery increased with NRS phase up to NRS phase 3A. Significant attenuation of recovery was observed for NRS phases 2A–3A. This attenuation increased with NRS phase (NRS phase 3A had the largest attenuation), a likely effect of higher enrollment performances and rates of recovery for higher NRS phases (see fig 2A). Time since SCI had a significant negative impact on recovery rates for patients in NRS phases 1C–3B, as noted by the negative coefficient estimates for the time since SCI-linear rate interaction. Time since SCI also significantly positively impacted the attenuation rate for NRS phase 2A–3A patients. The positive coefficients indicated that recovery curves for patients farther removed from their injuries were less severely attenuated. Thus, patients farther removed from injury recovered walking endurance less rapidly, but had flatter recovery curves because of lower attenuation (see fig 2B). We also note that there was no significant longitudinal progress for patients in NRS phases 1A and 3C. No rate parameters were significant in NRS phase 1A, and rate parameters in NRS phase 3C were highly variable and imprecise, despite statistically significant time since SCI/NRS phase/rate interaction terms. NRS phase 1A patients, on average, did not respond and NRS phase 3C patients possessed very high average walking endurance at enrollment (average distance, 310m) with little room for improvement.

### Ten-Meter Walk Test Model

Speeds for the 10-meter walk significantly improved from enrollment at an attenuated rate (see table 2, 10MW column). Performance at enrollment significantly varied by NRS phase, AIS grade, assistive device used, and age, but not by time since SCI. NRS phase, time since SCI, and AIS grade significantly impacted the rate of improvement. NRS phase and time since SCI also significantly impacted the rate of attenuation, but not AIS grade. Age and device type were not significant determinants of rates of improvement or attenuation. The impact of time since SCI on the rates of improvement and attenuation significantly varied by NRS phase, as judged by 3-way interactions involving these terms. No other 3-way interaction terms were of significance.
Estimates of significant model fixed effects (table 5) and plots of average recovery curves (fig 3) for 10-meter walk test speeds showed that NRS phases significantly differed and were generally well-ordered with respect to performance at enrollment. NRS phase 1A patients walked slowest on average (.05m/s) and NRS phase 3C patients walked fastest (1.06m/s). Patients with AIS grade D injuries walked on average .10m/s (95% CI, 10–55m) faster than AIS grade C patients at enrollment. Age significantly impacted enrollment performance, because every 1-year increase in age at enrollment corresponded with a .002m/s (95% CI, –.004 to .000m/s) decrease in enrollment speed. Patients walked .04m/s (95% CI, –.07 to –.02m/s) slower with their current assistive walking device than their initial device.
Table 5Coefficient Estimates of Significant Fixed Effects With 95% CIs From the Linear Mixed Effects Model of NRN 10-Meter Walk Test Data
NRS Phase at EnrollmentEnrollment Performance (intercept)Recovery Rate (linear)Attenuation Rate (quadratic)Time SCI-Recovery Rate InteractionTime SCI-Attenuation Rate Interaction
1A0.05 (–0.07 to 0.18)–0.02 (–0.07 to 0.03)0.01 (0 to 0.01)0.01 (–0.03 to 0.05)0 (–0.01 to 0)
1B0.08 (–0.03 to 0.18)0.03 (–0.01 to 0.06)0 (–0.01 to 0)0 (–0.02 to 0.03)0 (0 to 0)
1C0.11 (0.01 to 0.21)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.04 (0.01 to 0.08)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.01 (–0.01 to 0)–0.03 (–0.06 to 0)0 (0 to 0.01)
2A0.25 (0.14 to 0.37)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.13 (0.09 to 0.18)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.01 (–0.02 to –0.01)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.08 (–0.11 to –0.05)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.008 (0.003 to 0.013)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
2B0.50 (0.36 to 0.64)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.18 (0.12 to 0.24)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.02 (–0.03 to –0.01)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.16 (–0.21 to –0.11)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.03 (0.02 to 0.04)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
2C0.45 (0.32 to 0.59)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.22 (0.16 to 0.28)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.02 (–0.03 to –0.01)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.10 (–0.14 to –0.06)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.010 (0.001 to 0.018)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
3A0.87 (0.72 to 1.01)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.20 (0.10 to 0.30)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.04 (–0.06 to –0.01)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.19 (–0.27 to –0.11)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.033 (0.004 to 0.061)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
3B0.92 (0.76 to 1.09)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.12 (–0.04 to 0.29)–0.01 (–0.08 to 0.05)–0.10 (–0.24 to 0.03)0 (–0.06 to 0.06)
3C1.06 (0.82 to 1.29)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
–0.22 (–0.68 to 0.25)0.25 (–0.07 to 0.57)–0.54 (–0.94 to –0.15)0.48 (0.12 to 0.84)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
0.03 (0.01 to 0.05)
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
NSNSNS
NOTE. All rate parameters are given as rates per 20 sessions, approximately corresponding to rates per NRN evaluation. Time since SCI interaction effects are given in log-years. Intercepts are averages over study sites.
Abbreviation: NS, not a significant model term.
Statistical significance (P<.05) of the marginal test of the coefficients against null value of zero.
Significant recovery was seen for patients in NRS phases 1C–3A, and the rate of recovery increased with NRS phase up to NRS phase 2C. Significant attenuation of recovery was observed for NRS phases 2A–3A. This attenuation increased with NRS phase (NRS phase 3A had the largest attenuation), an effect of higher enrollment performances and rates of recovery for higher NRS phases (see fig 3A). Time since SCI had a significant negative impact on recovery rates for patients in NRS phases 2A–3A, as noted by the negative coefficient estimates for the time since SCI-linear rate interaction. Time since SCI also significantly positively impacted the attenuation rate for NRS phases 2A–3A patients. Patients farther removed from injury thus recovered walking speed less rapidly, but had flatter recovery curves because of lower attenuation (see fig 3B). We also note that there was no significant longitudinal progress for patients in NRS phases 1A, 1B, and 3C. No rate parameters were significant in NRS phases 1A and 1B, and rate parameters in NRS phase 3C were highly variable and imprecise, despite statistically significant time since SCI/NRS phase/rate interaction terms. NRS phases 1A and 1B patients, on average, did not respond, and NRS phase 3C patients possessed very high average walking speed at enrollment (average speed, 1.06m/s), with little room for improvement. Injury severity (AIS grade) also significantly impacted the rate of recovery (see fig 3C). The rate of improvement for AIS grade D patients was significantly higher (.03m/s per session; 95% CI, .01–.05) than that for AIS grade C patients.

### Model Fit

We empirically examined the quality of fit for the final model for each outcome measure (fig 4). Residual SDs were 3.8 for the Berg Balance Scale model, 39.6m for the 6-minute walk test model, and .15m/s for the 10-meter walk test model, representing 6.8%, 5.9%, and 5.4% of the observed ranges for the outcome measures, respectively. There were scattered outlying observations (see fig 4), but overall the models exhibited fairly strong predictive value. Seventy-nine percent of predicted Berg Balance Scale values were within 3 points of observed values, 89% of predicted 6-minute walk test speeds were within 40m of observed speeds, and 84% percent of predicted 10-meter walk test speeds were within 0.1m/s of observed speeds. These statistics represent the quality of fit as measured on the existing data; wider margins would occur for the prediction of future observations. We briefly note that the distribution of the residuals for each model (not shown here), which are assumed to be normal in a mixed effects model, was heavy-tailed relative to the normative distribution, but also symmetric. This phenomenon has been shown to not substantially impact the estimation of the fixed effects but to inflate SEs of the fixed effects estimates, thereby rendering hypothesis tests more conservative than under a normal model.
• Pinheiro J.
• Bates D.
Mixed effects models in S and S-Plus.

## Discussion

The recovery models we have presented provide specific, previously unavailable details about longitudinal functional recovery of clinically incomplete SCI patients, including that time since SCI and NRS phase significantly impacted patterns of functional recovery. Patients who were further removed from their SCI recovered at slower rates. Further, patient age was not a significant predictor of recovery rates and only weakly correlated (ρ=.13) with time since SCI, and therefore the impact of time since SCI on recovery rates was independent of patient age. Patients in lower NRS phase groups exhibited lower balance and walking function, a phenomenon which has previously been demonstrated,
• Behrman A.L.
• Ardolino E.
• VanHiel L.R.
• et al.
Assessment of functional improvement without compensation reduces variability of outcome measures after human spinal cord injury.
but also recovered over time at slower rates. NRS phase was a substantially better predictor of recovery than AIS grade, but AIS grade did possess some predictive value. The addition of AIS grade to models including NRS phase significantly improved model fit for the Berg Balance Scale and 10-meter walk test, in that patients with AIS grade C injuries performed at a lower level and recovered along a slower trajectory than those with AIS grade D injuries.
The physiologic state of the spinal circuitry may have contributed to the rates of recovery both in regard to time since injury and the extent of recovery. In the case of time since injury, neuromuscular plasticity conceivably continues to occur over time, including deleterious changes
• Dietz V.
Neuronal plasticity after a human spinal cord injury: positive and negative effects.
and restoring the functional reorganization for behavioral changes in response to task-specific training, and thus would require more training the longer the intervention was delayed. For those within 1 year of injury, natural or spontaneous recovery should also be considered to contribute to the early recovery trajectory. However, for those several years postinjury, this would not likely be a contributing factor. For those individuals who initially could not maintain trunk stability or stand, the time to recovery standing balance and walking would also be longer, because those preceding functions would first need to be restored.
These models provide prognostic value to individuals receiving locomotor training based on simple baseline patient characteristics, that is, NRS phase, AIS grade, and time since SCI. For example, while NRS phase 2 patients can expect to see the largest gains in balance function, and NRS phase 3 patients the largest gains in walking ability, both can expect to see significant improvements in balance and walking function. These results correspond with previous results about the functional status of NRS phase groups.
• Behrman A.L.
• Ardolino E.
• VanHiel L.R.
• et al.
Assessment of functional improvement without compensation reduces variability of outcome measures after human spinal cord injury.
NRS phase 2 patients generally have standing capability and are in the early stages of relearning to walk, and therefore the immediate improvements after enrollment in the NRN were in regaining balance as a precursor to gaining walking function. NRS phase 3 patients had high balance function (ie, high Berg Balance Scale scores) and were generally ambulatory, although they exhibited gait deviations and/or were limited in speed and endurance. Hence, there was little room for improvement in balance capability, as measured by the Berg Balance Scale, whereas there was significant room for improvement in walking function. NRS phases 1B and 1C patients exhibited significant improvement in balance and walking ability as well, albeit at much lower rates than those observed for NRS phases 2 and 3 patients. At each end of the NRS spectrum—NRS phase 1A and phase 3C—there were no significant longitudinal changes in balance and walking ability. This may have been reflective of the high functional status of NRS phase 3C patients at enrollment, who had little room for improvement in these measures, and the limited gains made by NRS phase 1A patients because of the magnitudes of their deficits.
These results also suggest that recovery is a multifaceted process and a single outcome measure is not sufficient to adequately capture recovery after incomplete SCI,
• Forrest G.F.
• Lorenz D.J.
• Hutchinson K.
• et al.
Ambulation and balance outcomes measure different aspects of recovery in individuals with chronic incomplete spinal cord injury.
an important point to consider when designing rehabilitation clinical trials. We observed dramatic variability in recovery patterns for all patients, and in particular even after stratifying by AIS grade. Such high variability in recovery can predictably result in false negatives when a novel therapeutic intervention is compared with a control. As shown by our results, a substantial amount of this variability is explained by NRS phase and time since SCI. Thus, NRS phase and time since SCI should be considered when setting expectations for recovery because of therapy and in determining the number of therapy sessions needed to achieve a clinical effect.
Our models give a general indication of the expected recovery patterns for patients receiving locomotor training and can serve as a benchmark for functional recovery in locomotor training programs. Although these models can be used as an aid in planning rehabilitation for patients with SCI, including making decisions about treatment, and monitoring progress in a locomotor program, it is important to appreciate that such models are based on averages on a large number of patients and thus will not precisely predict the same outcomes for every individual. To our knowledge, this is the first examination of recovery over time by fitting longitudinal models for individuals with SCI receiving a standardized therapy intervention. Evidence-based practice is highly supported by the American Physical Therapy Association and rehabilitation specialists; however, randomized controlled trials are very difficult to conduct and are very expensive. We suggest that these models can be considered as a benchmark for other studies that conduct standardized rehabilitation therapies in large cohorts and may be used for clinical decision-making in the future.

### Study Limitations

The principal limitation of our analyses is the imposition of a linear rather than nonlinear model to our data. Nonlinear models are selected to model the process generating the data, in this case the recovery process, while linear models are often selected to model the data, that is, to provide the best fit to the observed data.
• Pinheiro J.
• Bates D.
Mixed effects models in S and S-Plus.
This has several implications, particularly with regard to predictions generated from a linear model. In the present context, our models were quadratic in time, and therefore we inherently imposed the characteristics of quadratic recovery to our data—a high initial rate of recovery coupled with a plateau. The quadratic model effectively forces the highest rate of recovery to occur at the start of therapy, which may or may not be the case. Further, we tacitly ignore the behavior of the curve after the plateau, because a quadratic curve by definition declines after reaching its plateau and our data did not support this behavior. Thus, predictions given by our models after their plateaus are inherently unreliable.
While not beset with the aforementioned issues of linear models, nonlinear models require large amounts of data and, more importantly, a well-defined theoretical/functional form for the model. Our use of linear rather than nonlinear models was in large part motivated by the absence of any such a theoretical pattern of longitudinal recovery, because the clinical assumption has been that recovery rarely occurs in chronic SCI patients. Our recent study demonstrated recovery in these SCI patients,
• Harkema S.J.
• Lorenz D.J.
• Edgerton V.R.
• Behrman A.L.
Balance and ambulation improvements in individuals with chronic incomplete spinal cord injury using locomotor training-based rehabilitation.
and plots of the temporal patient profiles therein illustrated substantial variability in recovery, even within the NRS phase, that belied any mechanistic pattern. The selection of a quadratic recovery curve did seem theoretically reasonable, because patients do not interminably recover, and graphical examination of individual-level recovery curves supported attenuated, quadratic-shaped recovery. Despite the polynomial and high-order interaction fixed effects included in our model, the fixed effects were reasonably interpretable. Our temporally quadratic models were superior to linear models (P<.001), and models with higher order polynomial terms (eg, cubic) did not substantially improve model fit (P>.14). To assert that functional recovery was necessarily quadratic in time would be an overstatement of our results and would ignore the issues described in the preceding paragraph. It is reasonable to say, however, that among linear models, a quadratic curve makes most theoretical sense (attenuated recovery) and provides convincingly better fits than other candidate linear models.
In selecting explanatory variables for our model of functional recovery over time, we considered variables that were likely to impact the functional outcomes, namely, age, time since SCI, AIS grade, NRS phase, and assistive walking device (for the 6-minute walk and 10-meter walk tests). Other possible covariates, such as neurologic level of injury, AIS motor and sensory scores, and number of sessions per week (ie, treatment intensity) were not included and may warrant further investigation. Moreover, while we were interested in providing as accurate a recovery model as possible for our endpoints, we also sought to preserve some level of parsimony in model building to minimize the number of predictive factors and keep model parameters reasonably interpretable. It is possible that our restrictions on the fixed effects structure of our models missed characteristics of functional recovery. This needs further exploration on larger sets of data.

## Conclusions

We have fit longitudinal recovery models for patients with incomplete SCI receiving standardized locomotor training in the NRN for measures of balance and walking function. These models demonstrated that NRN patients experience significant improvement in balance and walking function and that NRS phase and time since SCI are important prognosticators of the rate of recovery over time. We have demonstrated that neuromuscular recovery score serves as a better predictor of recovery than AIS grade, although both can be used in conjunction. Our models demonstrated reasonable accuracy in prediction, but additional data are needed for further validation, as well as to explore additional model covariates and make attempts at fitting nonlinear models. These models can also be extended to other SCI outcomes including measures of cardiovascular function, health, quality of life, and muscle tone and strength.

## Acknowledgments

We thank Joe Canose, Susan Howley, and Michael Mangienello from the Christopher and Dana Reeve Foundation for their dedication and support. The extraordinary vision, compassion, and dedication of Christopher and Dana Reeve made the NeuroRecovery Network possible. We also thank the other current or past NeuroRecovery Network Center Directors: Steve Ahr, Linda Shelburne, PT, and Mark Sheridan, MSW (Frazier Rehab Institute); Mary Schmidt-Read, MS, DPT (Magee Rehabilitation); Steve Williams, MD (Boston Medical Center); Daniel Graves, PhD (Memorial Hermann/The Institute of Rehabilitation and Research); Sarah Morrison, PT and Keith Tansey, MD, PhD (Shepherd Center); Gail F. Forrest, PhD (Kessler Medical Rehabilitation Research and Education Corporation); and D. Michele Basso, PT, EdD (The Ohio State University Medical Center) plus all other current and previous Network members. We also thank the critical review and editorial support of Jessica Hillyer, PhD, and the leadership, foresight, and support of the NRN Advisory Board, Reggie Egerton, PhD, Moses Chao, PhD, Michael Fehlings, MD, PhD, Andrei Krassioukov, MD, PhD, and Shelly Sorani, MA.

## APPENDIX 1. Technical Model Specification

Here we provide more technical specification for the recovery models for our 3 outcomes. Specifically, we list the final model equations, define the random effects variance-covariance structures, and define any variance heterogeneity structures.

### Berg Balance Scale

The model equation for the Berg Balance Scale was
$yijk (t)=(α+ρj+λk+ai)+(β+δj+γk+ζXijk+(δζ)jXijk+(γζ)kXijk+bi)t+(κ+ηj+θXijk+(ηθ)jXijk+ci)t2+εijk,$

which states that at time t, the Berg Balance Scale score for patient i in NRS phase j with AIS grade k is a function of the following:
• Intercept terms
• α – overall intercept
• ρj – NRS phase j modifier
• λk – AIS grade k modifier
• ai – intercept random effect for patient i
• Linear rate terms
• β – overall linear rate
• δj – NRS phase j modifier
• γk – AIS grade k modifier
• ζ – time since SCI modifier
• xijk – time since SCI of patient i in NRS phase j and AIS grade k
• (δζ)j – interaction effect of time since SCI and NRS phase j
• (γζ)k – interaction effect of time since SCI and AIS grade k
• bi – linear rate random effect for patient i
• κ – overall quadratic rate
• ηj – NRS phase j modifier
• θ – time since SCI modifier
• (ηθ)j – interaction effect of time since SCI and NRS phase j
• ci – quadratic rate random effect for patient i
• εijk – residual error
There were 3 random effects defined at the patient level, and we left the variance-covariance matrix unspecified, leaving 3 variance and 3 covariance/correlation parameters to be estimated. Residual variance was found to be nonhomogeneous over the NRS phases, and therefore a variance function was defined to estimate residual variance within each NRS phase, resulting in the estimation of 8 additional parameters.

### Six-Minute Walk Test

The random effects for the 6-minute walk test model differed from the Berg Balance Scale model in 1 aspect: there was an added layer of random effects at the assistive device within patient level. This is seen in the model equation:
$yijk (t)=(α+ρj+ai+ai,k)+(β+δj+ζxijk+(δζ)jXijk+bi)t+(κ+ηj+θXijk+(ηθ)jXijk+ci)t2+εijk,$

which states that at time t, the 6-minute walk test distance for patient i in NRS phase j using assistive walking device k is a function of the following:
• Intercept terms
• α – overall intercept
• ρj – NRS phase j modifier
• ai – intercept random effect for patient i
• ai,k – intercept random effect for device k within patient i
• Linear rate terms
• β – overall linear rate
• δj – NRS phase j modifier
• ζ – time since SCI modifier
• xijk – time since SCI of patient i in NRS phase j and AIS grade k
• (δζ)j – interaction effect of time since SCI and NRS phase j
• bi – linear rate random effect for patient i
• κ – overall quadratic rate
• ηj – NRS phase j modifier
• θ – time since SCI modifier
• (ηθ)j – interaction effect of time since SCI and NRS phase j
• ci – quadratic rate random effect for patient i
• εijk – residual error
There were 3 random effects defined at the patient level, and we left the variance-covariance matrix unspecified, leaving 3 variance and 3 covariance/correlation parameters to be estimated. There was an additional nested random effect—a nested random intercept at the assistive walking device within patient level, adding another parameter (random effects variance) to be estimated. Models considering random linear and quadratic rates at this nested level failed to improve the model fit (P=.98). As with the Berg Balance Scale model, residual variance was found to be nonhomogeneous over the NRS phases, and therefore a variance function was defined to estimate residual variance within each NRS phase, resulting in the estimation of 8 additional parameters.

### Ten-Meter Walk Test

The 10-meter walk test model was similar in definition to the 6-minute walk test model, with a few slight differences:x
$yijk (t)=(α+ρj+ai+ai,k)+(β+δj+ζxijk+(δζ)jXijk+bi)t+(κ+ηj+θXijk+(ηθ)jXijk+ci)t2+εijk,$

which states that at time t, the 10-meter walk test speed for patient i in NRS phase j using assistive walking device k is a function of the following:
• Intercept terms
• α – overall intercept
• ρj – NRS phase j modifier
• ai – intercept random effect for patient i
• ai,k – intercept random effect for device k within patient i
• Linear rate terms
• β – overall linear rate
• δj – NRS phase j modifier
• ζ – time since SCI modifier
• xijk – time since SCI of patient i in NRS phase j and AIS grade k
• (δζ)j – interaction effect of time since SCI and NRS phase j
• bi – linear rate random effect for patient i
• κ – overall quadratic rate
• ηj – NRS phase j modifier
• θ – time since SCI modifier
• (ηθ)j – interaction effect of time since SCI and NRS phase j
• ci – quadratic rate random effect for patient i
• εijk – residual error
The random effects and variance function specification were identical so that for the 6-minute walk test model—3 random effects defined at the patient level with an unspecified variance-covariance matrix, a nested random intercept at the assistive walking device within patient level, and a variance function defined to estimate residual variance within each NRS phase. Thus, 7 random effects variance parameters and 8 variance function parameters were estimated. The addition of linear and quadratic rate random effects at the device within patient level did not improve the fit of the model (P=.49).

## References

• Behrman A.L.
• Bowden M.G.
• Nair P.M.
Neuroplasticity after spinal cord injury and training: an emerging paradigm shift in rehabilitation and walking recovery.
Phys Ther. 2006; 86: 1406-1425
• Harkema S.J.
• Behrman A.L.
• Bratta A.
• Sisto S.A.
• Edgerton V.R.
Establishing the NeuroRecovery Network: multisite rehabilitation centers that provide activity-based therapies and assessments for neurologic disorders.
Arch Phys Med Rehabil. 2012; 93: 1498-1507
• Harkema S.J.
• Lorenz D.J.
• Edgerton V.R.
• Behrman A.L.
Balance and ambulation improvements in individuals with chronic incomplete spinal cord injury using locomotor training-based rehabilitation.
Arch Phys Med Rehabil. 2012; 93: 1508-1517
• Tilling K.
• Sterne J.A.
• Rudd A.G.
• Glass T.A.
• Wityk R.J.
• Wolfe C.D.
A new method for predicting recovery after stroke.
Stroke. 2001; 32: 2867-2873
• Mahoney F.I.
• Barthel D.
Functional evaluation: the Barthel Index.
Md State Med J. 1965; 14: 61-65
• Dobkin B.
• Apple D.
• Barbeau H.
• et al.
Weight-supported treadmill vs over-ground training for walking after acute incomplete SCI.
Neurology. 2006; 66: 484-493
• Visintin M.
• Barbeau H.
The effects of body weight support on the locomotor pattern of spastic paretic patients.
Can J Neurol Sci. 1989; 16: 315-325
• Wernig A.
• Nanassy A.
• Müller S.
Maintenance of locomotor abilities following Laufband (treadmill) therapy in para- and tetraplegic persons: follow-up studies.
Spinal Cord. 1998; 36: 744-749
• Dietz V.
• Colombo G.
• Jensen L.
• Baumgartner L.
Locomotor capacity of spinal cord in paraplegic patients.
Ann Neurol. 1995; 37: 574-582
• Behrman A.L.
• Harkema S.J.
Physical rehabilitation as an agent for recovery after spinal cord injury.
Phys Med Rehabil Clin N Am. 2007; 18: 183-202
• Harkema S.J.
Plasticity of interneuronal networks of the functionally isolated human spinal cord.
Brain Res Rev. 2008; 57: 255-264
• Behrman A.K.
• Lawless-Dixon A.R.
• Davis S.B.
• et al.
Locomotor training progression and outcomes after incomplete spinal cord injury.
Phys Ther. 2005; 85: 1356-1371
• Berg K.O.
• Wood-Dauphinee S.L.
• Williams J.I.
• Maki B.
Measuring balance in the elderly: validation of an instrument.
Can J Public Health. 1992; 83: S7-S11
• Berg K.O.
• Maki B.E.
• Williams J.I.
• Holliday P.J.
• Wood-Dauphinee S.L.
Clinical and laboratory measures of postural balance in an elderly population.
Arch Phys Med Rehabil. 1992; 73: 1073-1080
• van Hedel H.J.
• Wirz M.
• Dietz V.
Assessing walking ability in subjects with spinal cord injury: validity and reliability of 3 walking tests.
Arch Phys Med Rehabil. 2005; 86: 190-196
• Behrman A.L.
• Ardolino E.
• VanHiel L.R.
• et al.
Assessment of functional improvement without compensation reduces variability of outcome measures after human spinal cord injury.
Arch Phys Med Rehabil. 2012; 93: 1518-1529
• American Spinal Injury Association
Reference manual for the International Standards for Neurological Classification of Spinal Cord Injury.
American Spinal Injury Association, Chicago2003
• Marino R.J.
• Barros T.
• Biering-Sorensen F.
• et al.
International standards for neurological classification of spinal cord injury.
J Spinal Cord Med. 2003; 26: S50-S56
• Verbecke G.
• Molenberghs G.
Linear mixed models for longitudinal data.
(Corrected edition) Springer-Verlag, New York2001
• Pinheiro J.
• Bates D.
Mixed effects models in S and S-Plus.
Springer-Verlag, New York2002
• Akaike H.
A new look at the statistical model identification.
IEEE Trans Automat Contr. 1974; 19: 716-723
• R Development Core Team
R: a language and environment for statistical computing.
R Foundation for Statistical Computing, Vienna2009
• Pinheiro J.
• Bates D.
• Debroy S.
• Sarkar D.
• The R Core Team
nlme: Linear and nonlinear mixed effects models.
Version R package version 3.1-96. 2009; ([computer program])
• Dietz V.
Neuronal plasticity after a human spinal cord injury: positive and negative effects.
Exp Neurol. 2012; 235: 110-115
• Forrest G.F.
• Lorenz D.J.
• Hutchinson K.
• et al.
Ambulation and balance outcomes measure different aspects of recovery in individuals with chronic incomplete spinal cord injury.
Arch Phys Med Rehabil. 2012; 93: 1553-1564