# Automated Stroke Rehabilitation Assessment using Wearable Accelerometers in
Free-Living Environments

Xi Chen  School of Mathematics, Statistics & Physics, Newcastle
UniversityNewcastleUK ,  Yu Guan  Open Lab, School of Computing, Newcastle
UniversityNewcastleUK ,  Jian Qing Shi  School of Mathematics, Statistics &
Physics, Newcastle UniversityNewcastleUK ,  Xiu-Li Du  School of Mathematial
Sciences, Nanjing Normal University, China, and School of Mathematics,
Statistics & Physics, Newcastle University, UK and  Janet Eyre  Institute of
Neurosciences, Newcastle UniversityNewcastleUK

###### Abstract.

Stroke is known as a major global health problem, and for stroke survivors it
is key to monitor the recovery levels. However, traditional stroke
rehabilitation assessment methods (such as the popular clinical assessment)
can be subjective and expensive, and it is also less convenient for patients
to visit clinics in a high frequency. To address this issue, in this work
based on wearable sensing and machine learning techniques, we developed an
automated system that can predict the assessment score in an objective and
continues manner. With wrist-worn sensors, accelerometer data was collected
from 59 stroke survivors in free-living environments for a duration of 8
weeks, and we aim to map the week-wise accelerometer data (3 days per week) to
the assessment score by developing signal processing and predictive model
pipeline. To achieve this, we proposed two new features, which can encode the
rehabilitation information from both paralysed/non-paralysed sides while
suppressing the high-level noises such as irrelevant daily activities. We
further developed the longitudinal mixed-effects model with Gaussian process
prior (LMGP), which can model the random effects caused by different subjects
and time slots (during the 8 weeks). Comprehensive experiments were conducted
to evaluate our system on both acute and chronic patients, and the results
suggested its effectiveness.

Functional ability of upper limbs, wrist-worn accelerometer sensor, stroke
rehabilitation, CAHAI score, wavelet features, longitudinal data analysis

††copyright: usgovmixed

##  1\. Introduction

It is widely known that stroke is a worldwide health problem causing
disability and death (Donnan et al., 2008), and it occurs when a blood clot
cuts off oxygen supply to a region of the brain. Hemiparesis is a very common
symptom of post-stroke that is the fractional or intact paralysis of one side
of the body, i.e., the opposite side to where the blood clot occurred, and it
results in difficulties in performing activities, e.g., reduced arm movement.
Patients can recover some of their capabilities with intense therapeutic
input, so it is important to assess their recovery levels in time. There are
many approaches to assess patients’ recovery levels including brain imaging
(Wintermark et al., 2005), questionnaire-based (Ferrari et al., 2007), and
lab-based clinical assessment (Barreca et al., 2005b).

The brain imaging technique, is deemed as one of the most reliable approach,
which can provide the information of brain hemodynamics (Wintermark et al.,
2005). However, this approach requires special equipment and is very expensive
in cost. Questionnaire-based approaches investigate the functional ability
during a period using questionnaires, and it can be categorised into two
types: patient-completed and caregiver-completed (Ferrari et al., 2007).
Although it is much cheaper than brain imaging approaches, it may contain
high-level of bias. For instance, patients may not remember their daily
activities (i.e.,recall bias); the caregivers may not be able to observe the
patient all the time. These biases make questionnaire-based approaches less
precise. Lab-based clinical assessment approaches (Barreca et al.,
2005b)(Barreca et al., 2005a), on the other hand, provide an alternative
solution. The patients’ upper limb functionality will be assessed by
clinicians, e.g., by observing patients’ capabilities of finishing certain
pre-defined activities (Barreca et al., 2005b). Compared with braining imaging
or questionnaire-based approaches, the cost of lab-based clinical assessment
approaches is reasonable with high accuracy. However, this assessment is
normally taken in clinics/hospitals, which is not convenient for the patients,
making continuous monitoring less feasible.

In this work, we aim to build an automated stroke rehabilitation assessment
system using wearable sensing and machine learning techniques. Different from
the aforementioned approaches, our system can measure the patients objectively
and continuously in free-living environments. We collected accelerometer data
using wrist-worn accelerometer sensors, and designed compact features that can
capture rehabilitation-related movements, before mapping these features to
clinical assessment scores (i.e., the model training process). The trained
model can be used to infer recovery-level for other unknown patients. In free
living environments, there are different types of movements which may be
related to different frequencies. For example, activities such as running or
jumping may correspond to high-frequency signal, while sedentary or eating may
be low-frequency signal. In this study, instead of recognising the daily
activities, which is hard to achieve given limited annotation (i.e., no
frame/sample-wise annotation), we transformed the raw accelerometer data to
the frequency domain, where we design features that can encode the
rehabilitation-related movements. Specifically, wavelet transform (Donald B.
Percival, 2000) was used, and the wavelet coefficients can represent the
particular frequency information at certain decomposition scales. In (Preece*
et al., 2009), Preece et al. provided some commonly used wavelet features
extracted from accelerometer data. However, to capture stroke rehabilitation-
related activities, some domain knowledge should be taken into account to
design better features. After stroke, patients have difficulties in moving one
side (i.e., paralysed side) due to the brain injury, and data from paralysed
side tends to describe more about the upper limb functional ability, than the
non-paralysed side (i.e., normal side). However, such signals can be
significantly affected by personal behaviours or irrelevant daily activities,
and such noises should be suppressed before developing the predictive models.
Various wavelet features were studied, and we proposed two new types of
features that can encode information from both paralysed/non-paralysed sides,
before developing predictive models for stroke rehabilitation assessment.
Specifically, in this work our contributions can be summarised as follows:

Data and System:

The data was collected from 59 stroke survivors (26 acute patients, and 33
chronic patients) for a duration of 8 weeks. Both wrist-worn accelerometer
data (3 days per week in free living environment) and the clinically
assessment scores were collected, based on which we developed a system that
can provide automated rehabilitation assessment mechanism in an objective and
continuous manner.

Features:

We proposed two new types of compact wavelet-based features that can encode
information from both paralysed and non-paralysed sides to represent upper
limb functional abilities for stroke rehabilitation assessment. It can also
significantly suppress the influences of personal behaviours or irrelevant
daily activities for data collected in the noisy free-living environment.

Predictive Models:

We developed the predictive models using two approaches: linear fixed-effects
model(LM) and longitudinal mixed-effects model with Gaussian process prior
(LMGP). Compared with LM, LMGP can model the random effects due to the
heterogeneity nature among subjects in this 8-week longitudinal study.

Experiments:

Comprehensive experiments were designed to study the effectiveness of our
system. We empirically demonstrated the results with respect to feature subset
on modelling the mixed-effects of LMGP. Our system was evaluated on both acute
and chronic patients, and the results suggested the effectiveness of our
approach.

##  2\. Background and Related Work

As described in Sec.1, lab-based clinical assessment was one of the most
effective stroke rehabilitation assessment methods. In this section, we
introduce of the the lab-based approach named Chedoke Arm and Hand Activity
Inventory (CAHAI) scoring, based on which our automated system can be
developed. Some sensing and machine learning techniques for automated health
assessment were also introduced.

###  2.1. Chedoke Arm and Hand Activity Inventory (CAHAI)

Figure 1. Left: Clinical behaviour assessment for CAHAI scoring; Right: The
CAHAI score form (Barreca et al., 2006).

CAHAI scoring is a clinical assessment method for stroke rehabilitation, and
it is a fully validated measure (Barreca et al., 2006) of upper limb
functional ability with 9 tasks which are scored by using a 7-point
quantitative scale. In the assessment, the patient will be asked to perform 9
tasks, including opening a jar of coffee, drawing a line with a ruler, calling
911, etc. and the clinician will score these behaviours based on patient’s
performance at a scale from 1 (total assist weak) to 7 (complete independence
i.e., timely, safely) (Barreca et al., 2006). A task example ”call 911” is
shown in Fig. 1. Thus the minimum and maximum summation scores are 7 and 63
respectively.

###  2.2. Automated Behaviour Assessment using Wearables

Recently, wearable sensing and machine learning (ML) techniques were
comprehensively studied for automated health assessment. Compared with the
traditional assessment approaches (e.g., via self-reporting, clinical
assessment, etc.) which are normally subjective and expensive, the automated
systems may provide an objective, low-cost alternative, which can also be used
for continuous monitoring/assessment. Some automated systems were developed to
assess the behaviours of diseases such as Parkinson’s disease (Rehman et al.,
2019) (Hammerla et al., 2015), autism (Plötz et al., 2012), depression (Little
et al., 2020); or to monitor the health status such as sleep (Zhai et al.,
2020), fatigue (Bai et al., 2020), etc.

After collecting behaviour or physiological signals (e.g., accelerometers,
ECG, audio, etc.), assessment/monitoring models can be developed. For better
performance, feature engineering can be considered as an important step before
ML model development, which can reduce data redundancy while improving the
discrimination capacity (Plötz et al., 2012). However, it tends to be time-
consuming and sometimes may require domain knowledge (Zhai et al., 2020) (Bai
et al., 2020) (Rehman et al., 2019) for interpretable/clinically-relevant
feature extraction. On the other hand, when interpretability is less required,
deep learning methods can be directly applied to the raw signal for predictive
model development (Hammerla et al., 2015) (Zhai et al., 2020) (Bai et al.,
2020) (Little et al., 2020), yet it normally requires adequate data annotation
for better model generalisation.

###  2.3. Sensing Techniques for Automated Stroke Rehabilitation Monitoring

With the rapid development of the sensing/ML techniques, researchers also
started to apply various sensors for stroke rehabilitation monitoring. In
(Dolatabadi et al., 2017), Kinect sensor was used in a home-like environments
to detect the key joints such that stroke patients’ behaviour can be assessed.
In (Ganesh et al., 2018), a wireless surface Electromyography (sEMG) device
was used to monitor the muscle recruitment of the post-stroke patients to see
the effect of orthotic intervention. In clinical environments, five wearable
sensors were placed on the trunk, upper and forearm of the two upper limbs to
measure the reaching behaviours of the stroke survivors (Jung et al., 2018).
To monitor motor functions of stroke patients during rehabilitation sessions
at clinics, an ecosystem including a jack and a cube for hand grasping
monitoring, as well as a smart watch for arm dynamic monitoring was designed
(Bobin et al., 2019). These techniques can objectively assess/measure the
behaviours of the stroke patients, yet they are either limited to clinical
environments (Bobin et al., 2019)(Jung et al., 2018) (Ganesh et al., 2018) or
constrained environments (e.g., in front of a camera (Dolatabadi et al.,
2017)). Most recently, wrist-worn sensors were used for stroke rehabilitation
monitoring for patients in free-living environment (Halloran et al., 2019)
(Tang et al., 2020). However, the collected data included high-level of data
redundancy, and it is challenging to extract compact and discriminant features
for reliable assessment.

##  3\. Methodology

In this section, we introduced our method from data collection, data pre-
processing, feature design to predictive models. Our aim is to develop an
automated model which can map the free-living 3-day accelerometer data into
the CAHAI score. With the trained model, we can infer the CAHAI score in an
objective and continuous manner.

To achieve this, we first reduced the data redundancy via preprocessing and
design compact and discriminant features. Given the proposed features, a
longitudinal mixed-effects model with Gaussian Process prior (LMGP) was
developed, which can further reduce the impact of large variability (caused by
different subjects and time slots) for higher prediction results. The pipeline
of our system is demonstrated in Fig. 2.

###  3.1. Data Acquisition

##### Participants

Data was collected as part of a bigger research study which aims to use a
bespoke, professionally-written video game as a therapeutic tool for stroke
rehabilitation (Shi et al., 2013). Ethical approval was obtained from the
National Research Ethics Committee and all work undertaken was in accordance
with the Declaration of Helsinki. Written, informed consent from all the
subjects was obtained. A cohort of 59 stroke survivors, without significant
cognitive or visual impairment, were recruited for the study. Patients were
divided into two groups, i.e.,

  * •

Group 1: the acute patient group, consisting of 26 participants who enrolled
into the study within 6 months after stroke;

  * •

Group 2: the chronic patient group, was formed by 33 participants who were 6
months or more post onset of stroke.

It is also worth noting that these 59 patients have different paralysed sides.
They visited the clinic for the CAHAI scoring every week (a random day in
weekdays) for a duration of 8 weeks. In the 8 weeks, they were asked to wear
two wrist-worn sensors for 3 full days (including night time) a week. They
were also advised to remove the device during shower or swimming. Since some
patients needed time to get familiar to this data collection procedure, for
better data quality we did not use the first week’s accelerometer data. The
first week’s CAHAI scores were used as medical history information.

Figure 2. The pipeline of our automated stroke rehabilitation assessment
system

Figure 3. The accelerometer data collected from AX3 device.

##### Data collection

In contrast to other afore-mentioned sensing techniques (Jung et al.,
2018)(Bobin et al., 2019)(Ganesh et al., 2018)(Dolatabadi et al., 2017), in
this study we collected the accelerometer data from wrist-worn sensors in
free-living environments. The sensor used for this study, i.e., AX3 (Axivity
Ltd, [n.d.]), is a triaxial accelerometer logger that was designed for
physical activity/behaviour monitoring, and it has been widely used in the
medical community (e.g., for the UK Biobank physical activity study (Doherty
et al., 2017)). The wrist bands were also designed such that the users can
comfortably wear it without affecting their behaviours. The data was collected
at 100Hz sampling rate, which can well preserve the daily activities of human
being (Bouten et al., 1997). Different from human activity recognition which
requires sample-wise or frame-wise annotation (Guan and Plötz, 2017) (Plötz
and Guan, 2018), the data collection in this study is relatively straight-
forward. The patients put on both wrist-worn sensors 3 full days a week,
before visiting clinicians for CAHAI scoring (i.e., week-wise annotation). In
other words, we aim to use accelerometer data captured in free-living
environments to represent the stroke survivors’ upper limb activities to
measure the degree of paresis (Jørgensen] et al., 1999) (i.e., CAHAI score).

One problem with most commercial sensors is that only summary data (e.g., step
count from fitbit), instead of raw data, are available. The algorithms of
producing summary data are normally non-open source, and may vary from vendor
to vendor – making the data collection and analysis device-dependent, and thus
less practical in terms of generalisation and scalability. The AX3 device used
in this study, on the other hand, outputs the raw acceleration information in
x, y, z directions. In Fig. 3, we visualise the raw acceleration data
(measured in g units $1\text{g}=9.8m/s^{2}$) collected from both wrists. It is
simple and transparent, making the collected data re-usable, which is crucial
for research communities.

###  3.2. Data pre-processing

For accelerometer data, signal vector magnitude (VM) (Karantonis et al., 2006)
is a popular representation, which is simply the magnitude of the triaxial
acceleration data defined as follows:

$a(t)=\sqrt{a_{x}^{2}(t)+a_{y}^{2}(t)+a_{z}^{2}(t)},$

where $a_{x}(t)$,$a_{y}(t)$,$a_{z}(t)$ are the acceleration along the x, y, z
axes at timestamp $t$. The gravity effect can be removed by

$VM(t)=\left|a(t)-1\right|.$

Because its simplicity and effectiveness, VM has been widely used in health
monitoring tasks, such as fall detection (Karantonis et al., 2006), physical
activity monitoring (Doherty et al., 2017), perinatal stroke assessment (Gao
et al., 2019), etc. To further reduce the data volume, we used second-wise VM,
i.e., the mean VM over each second (including 100 samples per second) will be
used as new representation. Some second-wise VM examples can be found in Fig.
4.

###  3.3. Feature Design

####  3.3.1. Challenges

We aim to build a model that can map the 3-day time-series data to the CAHAI
score. The original input data, however, is very big. Even if we used the
second-wise VM data, in each trial it still included roughly $3$ days $\times$
$24$h/day $\times$ $3600$s/h $=259200$ samples with one annotation (i.e.,
CAHAI score). Compared with the popular human activity recognition tasks,
which tend to have rich annotations (in frame-wise or sample-wise level), the
lack of annotation makes this task extremely challenging. Moreover, since the
data was collected in free-living environments, and the 3 full days (per week)
can be taken in weekdays or weekends, which may increase the data variability
significantly. With data redundancy, limited annotations and large
variability, it is unrealistic to use the popular machine learning or deep
learning approaches directly, and in this case feature engineering may play a
major role in extracting compact and discriminant signatures.

Figure 4.  VM data of two patients (paralysed side).  Figure 5. $\mathbf{SAD}$
representation of two patients (paralysed side)

####  3.3.2. Wavelet Features

For time-series analysis, wavelet analysis is a powerful tool to represent
various aspects of non-stationary signals such as trends, discontinuities, and
repeated patterns (Ayachi et al., 2016) (Donald B. Percival, 2000) (Preece* et
al., 2009), which is especially useful in signal compression or noise
reduction. Given its properties, wavelet features have been widely used in
accelerometer-based daily living activity analytics (Ayachi et al., 2016). In
this work, we used discrete wavelet transform (DWT) and discrete wavelet
packet transform (DWPT) as feature extractors, based on which new features
were designed to preserve the stroke rehabilitation-related information. More
details of DWT and DWPT can be found at Appendix 5.1. After applying the DWT
and DWPT, VM signals can be transformed to the wavelet coefficients at
different decomposition scales. Specifically, we used $\mathbf{SAD}$
representation, and its entry $SAD_{j}$ (normalised Sum of Absolute value of
DWT coefficients at scale $j$) can preserve the energy of the daily activities
at different decomposition scale $j$, where
$j\in\\{1.1,1.2,1.3,1.4,2,3,4,5,6,7\\}$ in this work. The technical details of
$\mathbf{SAD}$ representation can be found in Appendix 5.2. Through wavelet
transformation, the long sequence (e.g., VM data in Fig. 4) can be transformed
into compact representation (e.g., $\mathbf{SAD}$ features in Fig. 5), based
on which we further design features for CAHAI score regression.

####  3.3.3. Proposed Features

Based on the compact $\mathbf{SAD}$ representation, we aim to further design
effective features for reliable CAHAI score regression. To have a better
understanding of the behaviour patterns, we selected two different patients
(in terms of recovery level) for visualisation.

First, we visualised the paralysed side of patient la012 (with CAHAI score
55), and la040 (with CAHAI 26) using VM representation (Fig. 4) and
$\mathbf{SAD}$ representation (Fig. 5). From both figures, we can see the
limitations of both representations. Although VM can demonstrate distinct
patterns from both patients, it may be also related to the large intra-class
variability (e.g., personalised behaviour patterns). Moreover, the redundancy
as well as the high-dimensionality make it hard for modelling. On the other
hand, SAD has low dimensionality, yet both patients exhibited high-level of
similarity, indicating $\mathbf{SAD}$ of the paralysed side alone is not
enough for distinguishing patients with different recovery levels.

Figure 6. $\mathbf{SAD}$ representation of two patients (with both
paralysed/non-paralysed sides)

Figure 7. Two proposed $\mathbf{PNP}$ representations for two patients

Given the observations, we further visualised SAD features from both
paralysed/non-paralysed sides for both patients in Fig.6. We can see patient
la012 (with high recovery level) uses both hands (almost) equally while
patient la040 (with low recovery level) tends to use the non-paralysed side
more. These observations motivated us to design new features from both sides,
instead of the paralysed side alone. Specifically, we proposed two types of
features, which encode the ratio information between paralysed side and non-
paralysed side, namely, $PNP^{1}_{k}$ and $PNP^{2}_{k}$:

$\begin{split}&PNP^{1}_{k}=\frac{SAD_{k}^{p}}{SAD_{k}^{np}}\\\
&PNP^{2}_{k}=\frac{SAD_{k}^{np}-SAD_{k}^{p}}{SAD_{k}^{np}+SAD_{k}^{p}},\end{split}$

where $k=1.1,1.2,1.3,1.4,2,3,4,5,6,7$ and $p$ and $np$ referred to the
paralysed side and non-paralysed side respectively. We also visualised patient
la012 and patient la040 using the new proposed features $\mathbf{PNP^{1}}$ and
$\mathbf{PNP^{2}}$ in Fig. 7, from which we can see the proposed features can
well distinguish these two patients, in contrast to basic wavelet features SAD
(Fig. 5). Based on these distinguishable low-dimensional features, it is
feasible to build CAHAI regression models.

We listed 4 types of features, i.e., the original wavelet features extracted
from paralysed ($\mathbf{SAD^{p}}$) and non-paralysed sides
($\mathbf{SAD^{np}}$) separately, as well as the two new proposed features
($\mathbf{PNP^{1}}$ and $\mathbf{PNP^{2}}$). Based on 10 scales, we can form
40-dimensional feature vector, as shown in Table 1. However, there exist
certain level of noises and redundancy (especially on $\mathbf{SAD^{p}}$, and
$\mathbf{SAD^{np}}$ ), so it is crucial to develop feature selection mechanism
or powerful prediction model for higher performance.

Feature type | Dimension | feature entries for each type
---|---|---
$\mathbf{SAD^{p}}$ | 10 |  $SAD^{p}_{1.1}$, $SAD^{p}_{1.2}$, $SAD^{p}_{1.3}$, $SAD^{p}_{1.4}$, $SAD^{p}_{2}$, $SAD^{p}_{3}$, … , $SAD^{p}_{7}$
$\mathbf{SAD^{np}}$ | 10 |  $SAD^{np}_{1.1}$, $SAD^{np}_{1.2}$, $SAD^{np}_{1.3}$, $SAD^{np}_{1.4}$, $SAD^{np}_{2}$, $SAD^{np}_{3}$, … , $SAD^{np}_{7}$
$\mathbf{PNP^{1}}$ | 10 |  $PNP^{1}_{1.1}$, $PNP^{1}_{1.2}$, $PNP^{1}_{1.3}$ ,$PNP^{1}_{1.4}$, $PNP^{1}_{2}$, $PNP^{1}_{3}$, … , $PNP^{1}_{7}$
$\mathbf{PNP^{2}}$ | 10 |  $PNP^{2}_{1.1}$, $PNP^{2}_{1.2}$, $PNP^{2}_{1.3}$ ,$PNP^{2}_{1.4}$, $PNP^{2}_{2}$, $PNP^{2}_{3}$, … , $PNP^{2}_{7}$
Table 1. The wavelet features at 10 scales.

###  3.4. Predictive models

Based on the proposed wavelet features, we aim to develop predictive models
that can map features to the CAHAI score. Although we reduced the data
redundancy significantly, there still exists data noises, which may encode
irrelevant information. It is crucial to develop robust mechanism to select
the most relevant features, and here we used a popular feature selection
linear model (LASSO). To model the nonlinear random effects in the
longitudinal study, we also proposed to use a Gaussian Process (GP) regression
model.

It is worth noting that our model will also take advantage of the medical
history information (i.e., CAHAI score during the first visit) to predict
CAHAI scores for the rest 7 weeks (i.e., week 2 - week 8). From the
perspective of practical application, CAHAI score from the initial week
(referred to as $ini$) may be used as an important normalisation factor for
different individuals.

####  3.4.1. The linear fixed-effects model

Since there may exist some redundant or irrelevant features for the prediction
task, first we proposed to use LASSO (Least Absolute Shrinkage and Selection
Operator) for feature selection.

Given the $41$-dimensional input variables (40 wavelet features and 1 CAHAI
score from the initial week), first we standardised the data using z-norm, and
each feature entry $x_{k}$ will be normalised as
$x_{k}^{new}=(x_{k}-\overline{x})/s_{k}$, where $\overline{x}$ and $s_{k}$ are
the mean and standard deviation of the $k^{th}$ feature. Based on the
aforementioned model, namely LASSO, useful features can be selected, based on
which prediction model can be developed. For simplicity, we first used linear
model to predict the target CAHAI score ${y_{i}}$:

(1)
${y_{i}}=\mbox{\boldmath${x}$}_{i}^{\mathrm{T}}\mbox{\boldmath${\beta}$}+\epsilon_{i},\
\epsilon_{i}\sim N(0,\sigma^{2}),$

where $i$ stands for the $i^{th}$ trial from patients (out of all patients
during week 2 - week 8); $\mbox{\boldmath${x}$}_{i}$ represents the selected
feature vector; ${\beta}$ are the model parameter vector to be estimated, and
$\epsilon_{i}$ is the random noise term.

####  3.4.2. Longitudinal mixed-effects model with Gaussian process prior
(LMGP)

It is simple to use linear model for CAHAI score prediction. However, it
ignores the heterogeneity nature among subjects in this longitudinal study. To
model the heterogeneity, we proposed to use a nonlinear mixed-effects model
(Shi et al., 2012), which consists of the fixed-effects part and random-
effects part. Specifically, the random-effects part contributes mainly on
modelling the heterogeneity, making the the prediction process subject/time-
adaptive for longitudinal studies. The longitudinal mixed-effects model with
Gaussian Process prior (LMGP) is defined as follows:

(2)
${y_{i,j}}=\mbox{\boldmath${x}$}_{i,j}^{\mathrm{T}}\mbox{\boldmath${\beta}$}+g(\mbox{\boldmath${\phi}$}_{i,j})+\epsilon_{i,j},\
\epsilon_{i,j}\sim N(0,\sigma^{2}),$

where $i$,$j$ stand for the $i^{th}$ patient at the $j^{th}$ visit (from week
2 to week 8); $\epsilon_{i,j}$ refers to as independent random error and
$\sigma^{2}$ is its variance; In Eq(2),
$\mbox{\boldmath${x}$}_{i,j}^{T}\mbox{\boldmath${\beta}$}$ is the fixed-
effects part and $g(\mbox{\boldmath${\phi}$}_{i,j})$ represents the nonlinear
random-effects part, and the latter can be modelled using a non-parametric
Bayesian approach with a GP prior (Shi et al., 2012).

It is worth noting that in LMGP the fixed-effects part
$\mbox{\boldmath${x}$}_{i,j}^{\mathrm{T}}\mbox{\boldmath${\beta}$}$ explains a
linear relationship between input features and CAHAI, while the random-effects
part $g(\mbox{\boldmath${\phi}$}_{i,j})$ is used to explain the variability
caused by differences among individuals or time slots during different weeks.
By considering both parts, LMGP provides a solution of personalised modelling
for this longitudinal data analysis. In LMGP, it is important to select input
features to model both parts, and we referred them to as fixed-effects
features and random-effects features, respectively. The effect of the fixed-
effects features will be studied in the experimental evaluation section.

For LMGP training, we first ignored the random-effects part, and only
optimised the parameters $\hat{\mbox{\boldmath${\beta}$}}$ of the fixed-
effects part (via ordinary least squares OLS); With estimated parameters
$\hat{\mbox{\boldmath${\beta}$}}$, the residual
$r_{ij}=y_{ij}-\mbox{\boldmath${x}$}_{i,j}^{\mathrm{T}}\hat{\mbox{\boldmath${\beta}$}}=g(\mbox{\boldmath${\phi}$}_{i,j})+\epsilon_{i,j}$
can be calculated, from which we can model the random-effects

$g(\mbox{\boldmath${\phi}$}_{i,j})\sim
GP(0,K(\cdot,\cdot;\mbox{\boldmath${\theta}$})).$

In this paper we chose $K(\cdot,\cdot;\mbox{\boldmath${\theta}$})$ as the
following squared exponential (covariance) kernel function:

$K\left(\cdot,\cdot;\mbox{\boldmath${\theta}$}\right)=v_{0}\exp\left\\{-\frac{1}{2}\sum_{q=1}^{Q}w_{q}\left({\phi_{i,j,q}}-{\phi_{i,j,q}^{\prime}}\right)^{2}\right\\},$

and the model parameters $\mbox{\boldmath${\theta}$}=(v_{0},w_{1},...,w_{Q})$
can be estimated by using an empirical Bayesian method (and in this work we
used the R package GPFDA). It is worth noting that other covariance kernel
functions can also be applied and more information can be found in (Rasmussen
and Williams, 2006) (Shi and Choi, 2011).

Let $\mathbf{C}\in\mathbb{R}^{J\times J}$ be covariance matrix of
$g(\mbox{\boldmath${\phi}$}_{i,j})$, where $J$ referred to the total visit
number for each subject, and $\mathbf{C}$ can be calculated based on the
kernel function and its parameters ${\theta}$. Then the fitted value in
random-effects part can be calculated by

$\hat{g}(\mbox{\boldmath${\phi}$}_{i,j})=\mathbf{C}\left(\mathbf{C}+\sigma^{2}\mathbf{I}\right)^{-1}r_{ij},$

where $\mathbf{I}$ is the identity matrix. The variance of the random-effects
part is given by:

$\operatorname{Var}(\hat{g}(\mbox{\boldmath${\phi}$}_{i,j}))=\sigma^{2}\left(\mathbf{C}+\sigma^{2}\mathbf{I}\right)^{-1}\mathbf{C}.$

Finally, the fitted value of LMGP can be given by:

$\hat{y_{ij}}=\mbox{\boldmath${x}$}^{T}_{i,j}\hat{\mbox{\boldmath${\beta}$}}+\mathbf{C}\left(\mathbf{C}+\sigma^{2}\mathbf{I}\right)^{-1}r_{ij}.$

Given that, the $95\%$ predictive intervals can be written as:

$[\hat{y_{ij}}-1.96\sqrt{\operatorname{Var}(\hat{g}(\mbox{\boldmath${\phi}$}_{i,j}))},\
\hat{y_{ij}}+1.96\sqrt{\operatorname{Var}(\hat{g}(\mbox{\boldmath${\phi}$}_{i,j}))}].$

##  4\. Experimental Evaluation

In this section, several experiments were designed to evaluate the proposed
features as well as the proposed prediction system. The patients were splitted
into two groups according to the disease nature, i.e., the acute patient group
(26 subjects) and the chronic patient group (33 subjects), and experiments
were conducted on both group separately. Since CAHAI score prediction is a
typical regression problem, we used the root mean square error (RMSE) as the
evaluation metric, and lower mean RMSE values indicate better performance.

###  4.1. Evaluation of the Proposed Feature $\mathbf{PNP}$

In this subsection, we evaluated the effectiveness of the proposed
$\mathbf{PNP}$ features. One most straight-forward approach is to calculate
the correlation coefficients against the target CAHAI scores. In Table 2 we
reported the corresponding correlation coefficients ($PNP^{1}_{k}$, and
$PNP^{2}_{k}$ in 10 scales) for acute/chronic patients group. The correlation
coefficients of the original wavelet features (with paralysed side
$SAD^{p}_{k}$, and non-paralysed side $SAD^{p}_{k}$ in 10 scales) against
CAHAI score were also reported for comparison. From Table 2, we can see:

  * •

$\mathbf{PNP}$ features generally have higher correlation coefficients (than
$\mathbf{SAD}$) against the CAHAI scores.

  * •

for $\mathbf{PNP}$ features, from Scale $k=1.1$ to $k=5$ there are higher
correlations against the CAHAI scores.

  * •

for chronic patients, $\mathbf{SAD}$ features (on the non-paralysed side)
exhibit comparable correlation scores with $\mathbf{PNP}$ features.

These observations indicate the necessities of selecting useful features on
building the prediction system. Although $\mathbf{PNP}$ demonstrates more
powerful prediction capacity, in some cases, $\mathbf{SAD}$ (e.g., extracted
from the non-paralysed side) may also provide important information for a
certain population (e.g., chronic patients).

- |  | Acute | Patients |  |  | Chronic | Patients |
---|---|---|---|---|---|---|---|---
Scale (k) | $SAD^{p}_{k}$ | $SAD^{np}_{k}$ | $PNP^{1}_{k}$ | $PNP^{2}_{k}$ | $SAD^{p}_{k}$ | $SAD^{np}_{k}$ | $PNP^{1}_{k}$ | $PNP^{2}_{k}$
k=1.1 | -0.41 | 0.32 | 0.68 | -0.70 | 0.22 | 0.49 | 0.56 | -0.56
k=1.2 | -0.42 | 0.33 | 0.69 | -0.71 | 0.24 | 0.50 | 0.57 | -0.56
k=1.3 | -0.43 | 0.32 | 0.70 | -0.72 | 0.23 | 0.51 | 0.58 | -0.57
k=1.4 | -0.42 | 0.33 | 0.69 | -0.71 | 0.24 | 0.51 | 0.57 | -0.57
k=2 | -0.42 | 0.31 | 0.69 | -0.71 | 0.23 | 0.50 | 0.56 | -0.55
k=3 | -0.42 | 0.27 | 0.67 | -0.68 | 0.25 | 0.50 | 0.53 | -0.52
k=4 | -0.43 | 0.20 | 0.60 | -0.63 | 0.26 | 0.50 | 0.48 | -0.47
k=5 | -0.42 | 0.10 | 0.49 | -0.52 | 0.27 | 0.50 | 0.43 | -0.42
k=6 | -0.37 | -0.01 | 0.35 | -0.38 | 0.27 | 0.48 | 0.35 | -0.34
k=7 | -0.30 | -0.10 | 0.19 | -0.20 | 0.28 | 0.45 | 0.25 | -0.24
Table 2. Correlation coefficients of the wavelet features and CAHAI score.

(a) Acute patients.

(b) Chronic patients.

Figure 8. Cross-correlation of the candidate features

For better understanding the relationship between these features, we also
reported the cross-correlation between each feature pairs. Noting we also
included the medical history feature, i.e., the initial week-1 CAHAI score.
From Fig. 8, and we have the following observations:

  * •

For both patient groups, the $\mathbf{PNP}$ features are highly correlated.
$\mathbf{PNP}$ features within the same type ($\mathbf{PNP}^{1}$ or
$\mathbf{PNP}^{2}$) tend to be positively correlated, while $\mathbf{PNP}$
features from different types tend to be negatively correlated.

  * •

For acute patients, $\mathbf{SAD}$ features for each side (paralysed side
$\mathbf{SAD^{p}}$ or non-paralysed side $\mathbf{SAD^{np}}$) are highly
(positively) correlated, yet the $\mathbf{SAD}$ features from different sides
are less correlated. For chronic patients, however, $\mathbf{SAD}$ features
from both sides are highly (positively) correlated.

  * •

In general, $\mathbf{PNP}$ features, $\mathbf{SAD}$ features and the medical
history information $ini$ are less correlated, indicating them as potentially
complementary information to be fused.

Based on the above findings, it is clear that within each feature types, there
may exist high-level of feature redundancy, and it is necessary to select the
most relevant feature subsets. For acute and chronic patient groups, the
optimal feature subset may vary due to the different movement patterns (e.g.,
on paralysed/non-paralysed sides). Although the proposed $\mathbf{PNP}$
features can alleviate this problem to some extent, it is beneficial to
combine the less correlated features (i.e.,$\mathbf{PNP}$, $\mathbf{SAD}$, and
$ini$), and due to the feature redundancy, it is crucial to extract compact
representation for the prediction model development.

###  4.2. Evaluation of the Predictive Models

####  4.2.1. Feature Selection

Based on the feature correlation analysis in Sec. 4.1, it is important we
select the most relevant features from various sources (i.e., $\mathbf{PNP}$,
$\mathbf{SAD}$, and $ini$). Different from the correlation-based approach
which can select each feature independently (by the correlation coefficient),
LASSO can select the feature by solving a linear optimisation problem with
sparsity constraint, and it takes the relationship of the features into
consideration. Based on LASSO we selected the most important features for both
acute/chronic patients, as shown in Table 3.

| Acute Patients
---
Chronic Patients
|  $PNP^{2}_{3}$, $PNP^{1}_{6}$, $SAD_{2}^{np}$, $SAD_{1.2}^{p}$
---
$SAD_{6}^{np}$, $ini$
|  $PNP^{1}_{1.4}$, $SAD_{4}^{p}$, $SAD_{2}^{np}$, $PNP^{2}_{1.3}$
---
$PNP^{1}_{4}$, $PNP^{2}_{1.1}$, $ini$, $PNP^{1}_{6}$
$SAD_{1.4}^{np}$, $SAD_{6}^{np}$
Table 3. Selected features using LASSO

####  4.2.2. Performance of linear fixed-effects model

Based on the selected features, we further performed leave-one-patient-out
cross validation on these two patient groups respectively using the linear
fixed-effects model. As shown in Fig. 10, the prediction results of the
chronic patients (with mean RMSE 3.29) tend to be much better than the ones of
the acute group (with mean RMSE 7.24). One of the main reasons might be the
nature of the patient group. In Fig. 9, we plotted the clinical CAHAI
distribution (i.e., the ground truth CAHAI) from week 2 to week 8, and we can
see the clinical CAHAI scores are very stable for chronic patients. On the
other hand, for acute patients who suffered from stroke in the past 6 months,
their health statuses were less stable and affected significantly by various
factors, and in this case the simple linear fixed-effected model yields less
promising results.

Figure 9. Clinical assessed CAHAI distribution with respect to visit

Figure 10. Linear model prediction vs clinical CAHAI; Left: Acute patients
(RMSE 7.24); Right: Chronic patients (RMSE 3.29).

Figure 11. LMGP prediction vs clinical CAHAI; Left: Acute patients (RMSE
5.75); Right: Chronic patients (RMSE 3.12)

####  4.2.3. Performance of Longitudinal mixed-effects Model with Gaussian
Process prior (LMGP)

Based on the selected features, we also developed LMGP for both patient
groups. Here, we used the selected features (from Table 3) as the fixed-
effects features and random-effects features. Similar to the linear fixed-
effects model, we evaluated the performance based on leave-one-patient-out
cross validation, and the mean RMSE values were reported in Fig. 11, from
which can see LMGP can further reduced the errors when compared with the
fixed-effects linear model, with mean RMSE 5.75 for acute patients and 3.12
for chronic patients, respectively.

Based on LMGP, we also performed ”continuous monitoring” on 4 patients (two
for each patient group) from week 2 to week 8 by predicting the week-wise
CAHAI scores (with mean and $95\%$ confidence interval) in Fig. 12, which is
extremely helpful when uncertainty measurement is required.

Figure 12. Continues monitoring using LMGP for 4 patients (top: two chronic
patients; bottom: two acute patients)

####  4.2.4. On the fixed-effects part of LMGP

LMGP includes two key parts, i.e., the linear fixed-effects and the non-linear
random-effects part, and it is important to choose the key features for
modelling. Since the fixed-effects part measures the main (linear)
relationship between the input features and the predicted CAHAI, we studied
the corresponding feature subsets. For random-effects part, we used the full
LASSO features (as shown in Table 3).

To select the most important feature subset for the fixed-effects part
modelling, we ranked the features (from Table 3) based on two criteria: LASSO
coefficients, and correlation coefficients (between features and CAHAI, as
described in Sec.4.1). Table 4 demonstrates ranked features, and here only the
top $50\%$ features (i.e., top 3 features for acute patients and top 5
features for chronic patients) were used to model the fixed-effects part, and
the settings as well as the results were reported in Table 5.

Feature Ranking Criterion |  | Acute Patients
---
Chronic Patients
| LASSO Coefficients
---
(absolute value)
|  $PNP^{2}_{3}$, $PNP^{1}_{6}$, $SAD_{2}^{np}$, $SAD_{1.2}^{p}$
---
$SAD_{6}^{np}$, $ini$
|  $PNP^{1}_{1.4}$, $SAD_{4}^{p}$, $SAD_{2}^{np}$, $PNP^{2}_{1.3}$
---
$PNP^{1}_{4}$, $PNP^{2}_{1.1}$, $ini$, $PNP^{1}_{6}$
$SAD_{1.4}^{np}$, $SAD_{6}^{np}$
| Correlation Coefficients
---
(absolute value)
|  $PNP^{2}_{3}$, $ini$, $SAD_{1.2}^{p}$, $PNP^{1}_{6}$
---
$SAD_{2}^{np}$, $SAD_{6}^{np}$
|  $ini$, $PNP^{1}_{1.4}$, $PNP^{2}_{1.3}$, $PNP^{2}_{1.1}$
---
$SAD_{1.4}^{np}$, $SAD_{2}^{np}$, $PNP^{1}_{4}$, $SAD_{6}^{np}$
$PNP^{1}_{6}$, $SAD_{4}^{p}$
Table 4.  Feature importance ranking for acute/chronic patients.  Acute Patients |  | Fixed-effects features
---
| Random-effects features
---
RMSE
| full 6 features in Table 3
---
| full 6 features in Table 3
---
5.75
| top 3 features (Corr criterion in Table 4):
---
$PNP^{2}_{3}$, $ini$, $SAD_{1.2}^{p}$
| full 6 features in Table 3
---
5.37
| top 3 features (LASSO criterion in Table 4):
---
$PNP^{2}_{3}$, $PNP^{1}_{6}$, $SAD_{2}^{np}$
| full 6 features in Table 3
---
5.51
Chronic Patients |  | Fixed-effects features
---
| Random-effects features
---
RMSE
| full 10 features in Table 3
---
| full 10 features in Table 3
---
3.12
| top 5 features (Corr criterion in Table 4):
---
$ini$, $PNP^{1}_{1.4}$, $PNP^{2}_{1.3}$ $PNP^{2}_{1.1}$, $SAD_{1.4}^{np}$
| full 10 features in Table 3
---
3.20
| top 5 features (LASSO criterion in Table 4):
---
$PNP^{1}_{1.4}$, $SAD_{4}^{p}$, $SAD_{2}^{np}$ $PNP^{2}_{1.3}$, $PNP^{1}_{4}$
| full 10 features in Table 3
---
5.12
Table 5. LMGP’s fixed-effects part modelling results (RMSE) based on different
feature subsets

It is interesting to observe the performance may change substantially based on
different settings. Specifically, with the top feature subsets, modelling the
LMGP’s fixed-effects part can further reduce the errors for acute patients, in
contrast to chronic patients with increased errors. The top 5 features
selected via the LASSO criterion yields the worst performance for chronic
patients, and one possible explanation could be due to the lack of feature
$ini$ —–the initial health condition—–a major attribute for chronic patient
modelling (see Fig. 9).

####  4.2.5. Model comparison

Models |  | RMSE (Acute)
---
RMSE (Chronic)
Support vector regression (linear) | 7.47 | 3.25
Support vector regression (rbf) | 9.67 | 4.92
Random forest regression | 8.19 | 3.93
Linear fixed-effects model | 7.24 | 3.29
LMGP (ours) | 5.75 | 3.12
Table 6. Model comparison

Base on the selected features in Table 3, we also compared our model with
other baselines such as support vector regression (SVR) and random forest
regression(RF) for acute/chronic patient groups. Leave one subject out cross
validation and mean RMSE was used to measure the performance. We observed that
linear models yielded reasonable results (linear SVR and linear fixed-effects
model), while the non-linear baselines (SVR with rbf kernel, and RF) had the
worse performance. One of the explanation is the over-fitting effect, where
the trained non-linear models do not generalise well to the unseen
patients/environments in this longitudinal study setting. RF’s performance may
also suffer significantly from the low-dimensionality of the selected features
(6 features for acute patients and 10 features for chronic patients). Given
the simplicity of the linear models and the designed low-dimensional features,
linear models tend to suffer less from the over-fitting effect, with
reasonable results in these challenging environments. Compared with the
baselines, our LMGP can further model the longitudinal mixed-effects (i.e.,
with linear fixed-effect part and non-linear random-effects part), making the
system adaptive to different subjects/environments, with the lowest errors.

##  5\. Conclusions

In this work, we developed an automated stroke rehabilitation assessment
system using wearable sensing and machine learning techniques. We collected
accelerometer data using wrist-worn sensors, based on which we built models
for CAHAI score prediction, which can provide objective and continuous
rehabilitation assessment. To map the long time-series (i.e., 3-day
accelerometer data) to the CAHAI score, we proposed a pipeline which performed
data cleaning, feature design, to predictive model development. Specifically,
we proposed two compact features which can well capture the rehabilitation
characteristics while suppressing the irrelevant daily activities, which is
crucial on analysing the data collected in free-living environments. We
further employed LMGP, which can make the model adaptive to different subjects
and different time slots (across different weeks). Comprehensive experiments
were conducted on both acute/chronic patients, and very promising results were
achieved, especially on the chronic patient group. We also studied different
feature subsets on modelling the fixed-effects part in LMGP, and experiments
suggested the errors can be further reduced for the challenging acute patient
population.

## References

  * (1)
  * Axivity Ltd ([n.d.]) Axivity Ltd. [n.d.].  AX3, 3-Axis Logging Accelerometer.  https://axivity.com/product/ax3.  [Online; accessed July-2020].
  * Ayachi et al. (2016) Fouaz S Ayachi, Hung P Nguyen, Catherine Lavigne-Pelletier, Etienne Goubault, Patrick Boissy, and Christian Duval. 2016.  Wavelet-based algorithm for auto-detection of daily living activities of older adults captured by multiple inertial measurement units (IMUs).  _Physiological measurement_ 37, 3 (March 2016), 442—461.  https://doi.org/10.1088/0967-3334/37/3/442
  * Bai et al. (2020) Yang Bai, Yu Guan, and Wan-Fai Ng. 2020.  Fatigue Assessment Using ECG and Actigraphy Sensors. In _Proceedings of the 24rd International Symposium on Wearable Computers_ _(ISWC)_. Association for Computing Machinery, New York, NY, USA.  https://doi.org/10.1145/3410531.3414308
  * Barreca et al. (2005a) S. Barreca, P. Stratford, C. Lambert, L. Masters, and D. Streiner. 2005a.  Test-Retest Reliability, Validity, and Sensitivity of the Chedoke Arm and Hand Activity Inventory: a New Measure of Upper-Limb Function for Survivors of Stroke.  _Arch Phys Med Rehabil_ 86 (2005), 1616–1622.
  * Barreca et al. (2006) S. Barreca, P. Stratford, L. Masters, C. Lambert, J. Griffiths, and C. McBay. 2006\.  Validation of Three Shortened Versions of the Chedoke Arm and Hand Activity Inventory.  _Physiother. Can._ 58 (2006), 1–9.
  * Barreca et al. (2005b) Susan R Barreca, Paul W. Stratford, Cynthia L. Lambert, Lisa M. Masters, and David L Streiner. 2005b.  Test-retest reliability, validity, and sensitivity of the Chedoke arm and hand activity inventory: a new measure of upper-limb function for survivors of stroke.  _Archives of physical medicine and rehabilitation_ 86 8 (2005), 1616–22.
  * Bobin et al. (2019) Maxence Bobin, Franck Bimbard, Mehdi Boukallel, Margarita Anastassova, and Mehdi Ammi. 2019\.  SpECTRUM: Smart ECosystem for sTRoke patient’s Upper limbs Monitoring.  13 (02 2019).  https://doi.org/10.1016/j.smhl.2019.01.001
  * Bouten et al. (1997) Carlijn Bouten, Karel Koekkoek, Maarten Verduin, Rens Kodde, and Jan Janssen. 1997.  A Triaxial Accelerometer and Portable Data Processing Unit for the Assessment of Daily Physical Activity.  _IEEE transactions on bio-medical engineering_ 44 (04 1997), 136–47.  https://doi.org/10.1109/10.554760
  * Daubechies (2006) I. Daubechies. 2006\.  Orthonormal bases of compactly supported wavelets.  _Commun. Pure. Appl. Math_ (2006), 909–996.
  * Doherty et al. (2017) Aiden Doherty, Dan Jackson, Nils Hammerla, Thomas Plötz, Patrick Olivier, Malcolm H. Granat, Tom White, Vincent T. van Hees, Michael I. Trenell, Christoper G. Owen, Stephen J. Preece, Rob Gillions, Simon Sheard, Tim Peakman, Soren Brage, and Nicholas J. Wareham. 2017\.  Large Scale Population Assessment of Physical Activity Using Wrist Worn Accelerometers: The UK Biobank Study.  _PLOS ONE_ 12, 2 (02 2017), 1–14.  https://doi.org/10.1371/journal.pone.0169649
  * Dolatabadi et al. (2017) Elham Dolatabadi, Ying Zhi, Bing Ye, Marge Coahran, Giorgia Lupinacci, Alex Mihailidis, Rosalie Wang, and Babak Taati. 2017\.  The toronto rehab stroke pose dataset to detect compensation during stroke rehabilitation therapy. 375–381.  https://doi.org/10.1145/3154862.3154925
  * Donald B. Percival (2000) Andrew T. Walden Donald B. Percival. 2000\.  _Wavelet methods for time series analysis_ (1 ed.).  Cambridge University Press.  http://gen.lib.rus.ec/book/index.php?md5=6f8eb3e5445b7c12210cb5b1d50a40f3
  * Donnan et al. (2008) G. Donnan, M. Fisher, M. Macleod, and S. Davis. 2008\.  Stroke.  _Lancet_ 371, 2 (2008), 1612–1623.
  * Ferrari et al. (2007) Pietro Ferrari, Christine Friedenreich, and Charles Matthews. 2007\.  The Role of Measurement Error in Estimating Levels of Physical Activity.  _American journal of epidemiology_ 166 (11 2007), 832–40.  https://doi.org/10.1093/aje/kwm148
  * Ganesh et al. (2018) A. C. Ganesh, B. S. Renganathan, C. Rajakumaran, S. P. Preejith, K. Shubham, J. Jayaraj, and S. Mohanasankar. 2018.  Post-Stroke Rehabilitation Monitoring Using Wireless Surface Electromyography: A Case Study. In _2018 IEEE International Symposium on Medical Measurements and Applications (MeMeA)_. 1–6.
  * Gao et al. (2019) Yan Gao, Yang Long, Yu Guan, Anna Basu, Jessica Baggaley, and Thomas Ploetz. 2019\.  Towards Reliable, Automated General Movement Assessment for Perinatal Stroke Screening in Infants Using Wearable Accelerometers.  _Proc. ACM Interact. Mob. Wearable Ubiquitous Technol._ 3, 1, Article 12 (March 2019), 22 pages.  https://doi.org/10.1145/3314399
  * Guan and Plötz (2017) Yu Guan and Thomas Plötz. 2017.  Ensembles of Deep LSTM Learners for Activity Recognition Using Wearables.  _Proc. ACM Interact. Mob. Wearable Ubiquitous Technol._ 1, 2, Article 11 (June 2017), 28 pages.  https://doi.org/10.1145/3090076
  * Halloran et al. (2019) Shane Halloran, Lin Tang, Yu Guan, Jian Qing Shi, and Janet Eyre. 2019.  Remote Monitoring of Stroke Patients’ Rehabilitation Using Wearable Accelerometers. In _Proceedings of the 23rd International Symposium on Wearable Computers_ (London, United Kingdom) _(ISWC 19)_. Association for Computing Machinery, New York, NY, USA, 72–77.  https://doi.org/10.1145/3341163.3347731
  * Hammerla et al. (2015) Nils Y. Hammerla, James M. Fisher, Peter Andras, Lynn Rochester, Richard Walker, and Thomas Plotz. 2015\.  PD Disease State Assessment in Naturalistic Environments Using Deep Learning. In _Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence_ (Austin, Texas) _(AAAI)_. AAAI Press, 1742–1748.
  * Jørgensen] et al. (1999) Henrik [Stig Jørgensen], Hirofumi Nakayama, Hans Otto Raaschou, and Tom [Skyhøj Olsen]. 1999\.  Stroke: Neurologic and Functional Recovery The Copenhagen Stroke Study.  _Physical Medicine and Rehabilitation Clinics of North America_ 10, 4 (1999), 887 – 906.  https://doi.org/10.1016/S1047-9651(18)30169-4 A New Century Approach to Stroke Management and Rehabilitation.
  * Jung et al. (2018) H. Jung, J. Park, J. Jeong, T. Ryu, Y. Kim, and S. I. Lee. 2018\.  A wearable monitoring system for at-home stroke rehabilitation exercises: A preliminary study. In _2018 IEEE EMBS International Conference on Biomedical Health Informatics (BHI)_. 13–16.
  * Karantonis et al. (2006) D. M. Karantonis, M. R. Narayanan, M. Mathie, N. H. Lovell, and B. G. Celler. 2006\.  Implementation of a real-time human movement classifier using a triaxial accelerometer for ambulatory monitoring.  _IEEE Transactions on Information Technology in Biomedicine_ 10, 1 (2006), 156–167.
  * Little et al. (2020) Bethany Little, Ossama Alshabrawy, Daniel Stow, I. Ferrier, Roisin McNaney, Daniel Jackson, Karim Ladha, Cassim Ladha, Thomas Ploetz, Jaume Bacardit, Patrick Olivier, Peter Gallagher, and John O’Brien. 2020.  Deep learning-based automated speech detection as a marker of social functioning in late-life depression.  _Psychological Medicine_ (01 2020), 1–10.  https://doi.org/10.1017/S0033291719003994
  * Plötz et al. (2012) Thomas Plötz, Nils Y. Hammerla, Agata Rozga, Andrea Reavis, Nathan Call, and Gregory D. Abowd. 2012\.  Automatic Assessment of Problem Behavior in Individuals with Developmental Disabilities. In _Proceedings of the 2012 ACM Conference on Ubiquitous Computing_ (Pittsburgh, Pennsylvania) _(UbiComp)_. Association for Computing Machinery, New York, NY, USA, 391–400.  https://doi.org/10.1145/2370216.2370276
  * Plötz and Guan (2018) T. Plötz and Y. Guan. 2018.  Deep Learning for Human Activity Recognition in Mobile Computing.  _Computer_ 51, 5 (2018), 50–59.
  * Preece* et al. (2009) S. J. Preece*, J. Y. Goulermas, L. P. J. Kenney, and D. Howard. 2009.  A Comparison of Feature Extraction Methods for the Classification of Dynamic Activities From Accelerometer Data.  _IEEE Transactions on Biomedical Engineering_ 56, 3 (March 2009), 871–879.  https://doi.org/10.1109/TBME.2008.2006190
  * Rasmussen and Williams (2006) CE. Rasmussen and CKI. Williams. 2006.  _Gaussian Processes for Machine Learning_.  MIT Press, Cambridge, MA, USA. 248 pages.
  * Rehman et al. (2019) Rana zia ur Rehman, Silvia Din, Yu Guan, Alison Yarnall, Jian Shi, and Lynn Rochester. 2019\.  Selecting Clinically Relevant Gait Characteristics for Classification of Early Parkinson’s Disease: A Comprehensive Machine Learning Approach.  _Scientific Reports_ 9 (12 2019).  https://doi.org/10.1038/s41598-019-53656-7
  * Sekine et al. (1998) M. Sekine, T. Tamura, M. Ogawa, T. Togawa, and Y. Fukui. 1998.  Classification of acceleration waveform in a continuous walking record. In _Proceedings of the 20th Annual International Conference of the IEEE Engineering in Medicine and Biology Society. Vol.20 Biomedical Engineering Towards the Year 2000 and Beyond (Cat. No.98CH36286)_ , Vol. 3. 1523–1526 vol.3.  https://doi.org/10.1109/IEMBS.1998.747177
  * Shi et al. (2013) J.Q. Shi, Y. Cheng, J. Serradilla, G. Morgan, C. Lambden, G. Ford, C. Price, H. Rodgers, T. Cassidy, L. Rochester, and J.A. Eyre. 2013\.  Evaluating Functional Ability of Upper Limbs after Stroke Using Video Game Data. In _International Conference on Brain and Health Informatics_ _(Lecture Notes in Artificial Intelligence, Vol. 8211)_ , K. Imamura, S. Usui, T. Shirao, T. Kasamatsu, L. Schwabe, and N. Zhong (Eds.). Springer, 181–192.
  * Shi and Choi (2011) Jian Shi and Taeryon Choi. 2011.  _Gaussian Process Regression Analysis for Functional Data_.  London: Chapman and Hall/CRC.  https://doi.org/10.1201/b11038
  * Shi et al. (2012) J.Q. Shi, B. Wang, E.J. Will, and R.M. West. 2012\.  Mixed-effects Gaussian process functional regression models with application to dose response curve prediction.  _Statistics in Medicine_ 31, 26 (2012), 3165–3177.  https://doi.org/10.1002/sim.4502
  * Tang et al. (2020) Lin Tang, Shane Halloran, Jian Qing Shi, Yu Guan, Chunzheng Cao, and Janet Eyre. 2020\.  Evaluating upper limb function after stroke using the free-living accelerometer data.  _Statistical Methods in Medical Research_ (2020).  https://doi.org/10.1177/0962280220922259
  * Wintermark et al. (2005) Max Wintermark, Musa Sesay, Emmanuel Barbier, Katalin BorbÃ©ly, William P. Dillon, James D. Eastwood, Thomas C. Glenn, CÃ©cile B. Grandin, Salvador Pedraza, Jean-FranÃ§ois Soustiel, Tadashi Nariai, Greg Zaharchuk, Jean-Marie CaillÃ©, Vincent Dousset, and Howard Yonas. 2005.  Comparative Overview of Brain Perfusion Imaging Techniques.  _Stroke_ 36, 9 (2005), e83–e99.  https://doi.org/10.1161/01.str.0000177884.72657.8b
  * Zhai et al. (2020) Bing Zhai, Ignacio Perez-Pozuelo, Emma A. D. Clifton, Joao Palotti, and Yu Guan. 2020.  Making Sense of Sleep: Multimodal Sleep Stage Classification in a Large, Diverse Population Using Movement and Cardiac Sensing.  _Proc. ACM Interact. Mob. Wearable Ubiquitous Technol._ 4, 2, Article 67 (June 2020), 33 pages.  https://doi.org/10.1145/3397325

## Appendix

###  5.1. Discrete wavelet transform and discrete wavelet packet transform

The DWT procedure includes two parts: decomposition and reconstruction.
Decomposition part will be the main focus in this project. We now consider
more details of the DWT using matrix algebra:

(3) $\mathbf{W}=\mathcal{W}\mathbf{X},$

where ${W}$ is the output of matrix of DWT coefficients in different scales.
${W}$ is the orthonormal matrix containing different orthonormal wavelet bases
(more details can be checked in (Daubechies, 2006) and (Donald B. Percival,
2000)) and it satisfies $\mathcal{W}^{T}\mbox{\boldmath${W}$}=\mathbf{I}_{n}$.
$\mathbf{X}$ is the raw signal. The signal $\mathbf{X}$ with length $N=2^{J}$,
the $N\times N$ orthonormal matrix ${W}$ can be separated into J+1
submatrices, each of which can produce a partitioning of the vector
$\mathbf{W}$ of DWT coefficients in each scale j, j = 1,2,…, J. To be more
specific, Eq(3) can be rewritten as follows:

(4) $\mathcal{W}\mathbf{X}=\begin{bmatrix}\mathcal{W}_{1}\\\
\mathcal{W}_{2}\\\ \vdots\\\ \mathcal{W}_{J}\\\ \mathcal{V}_{J}\\\
\end{bmatrix}\textbf{X}=\begin{bmatrix}\mathcal{W}_{1}\mathbf{X}\\\
\mathcal{W}_{2}\mathbf{X}\\\ \vdots\\\ \mathcal{W}_{J}\mathbf{X}\\\
\mathcal{V}_{J}\mathbf{X}\\\ \end{bmatrix}=\begin{bmatrix}\mathbf{W}_{1}\\\
\mathbf{W}_{2}\\\ \vdots\\\ \mathbf{W}_{J}\\\ \mathbf{V}_{J}\\\
\end{bmatrix}=\mathbf{W},$

where $\mathbf{W}_{j}$ is a column vector of length $N/{2^{j}}$ representing
the differences in adjacent weighted averages from scale 1 to scale J,
$\mathbf{V}_{J}$ is the last column contained in $\mathbf{W}$ which has the
same length with $\mathbf{W}_{J}$. $\mathbf{W}_{j}$ is defined as detailed
coefficients at scale j. $\mathbf{V}_{J}$ contains the approximated
coefficients at the J-th level. $\mathcal{W}_{j}$ has dimension
$N/{2^{j}}\times N$, where j = 1,2,…,J and $\mathbf{V}_{J}$ has the same
dimension with $\mathbf{W}_{J}$. Note that the rows of design orthonormal
matrix $\mathcal{W}$ depend on the decomposition level j-th. In other words,
the value of J depends on the DWT decomposition scale of the raw signal. The
maximum decomposition level j equals J since our signal $\mathbf{X}$ has
length $N=2^{J}$.

We now further consider wavelet packet transform DWPT. The DWPT is the
expansion of the discrete wavelet transformation. In DWT, each scale is
calculated by passing only the previous wavelet approximated coefficients
through discrete-time low and high pass quadrature mirror filters. However, in
the DWPT, both the detailed and approximation coefficients are decomposed to
create the full binary tree. More details can be found in (Donald B. Percival,
2000).

Finally, at 3 scales of decomposition as example, a very striking plot will be
given in Fig.13, where the difference between DWT and DWPT will be seen very
clearly.

Figure 13. Comparing with the DWT and DWPT at three scales decomposition.

###  5.2. Commonly used wavelet features

In the discrete wavelet transform (DWT), $\mathbf{W}_{j}$ represents DWT
coefficients in the j-th decomposition scale. DWT can be written as $\
\mathbf{W}=\mathcal{W}\mathbf{X}$, where $\mathbf{W}$ is a column vector with
length $2^{j}$ and
$\mathbf{W}=[\mathbf{W}_{1},\mathbf{W}_{2},...,\mathbf{W}_{J},\mathbf{V}_{J}]^{\mathrm{T}}$,
$\mathcal{W}$ is the orthonormal matrix which satisfies
$\mathcal{W}^{T}\mathcal{W}=\mathbf{I}_{n}$ and contains different filters.
Due to the orthonormality of DWT, which means that
$\mathbf{X}=\mathcal{W}^{\mathrm{T}}\mathbf{W}$ and
$\left\|\mathbf{X}\right\|{{}^{2}}=\left\|\mathbf{W}\right\|{{}^{2}}$,
$\left\|\mathbf{W}_{j}\right\|{{}^{2}}$ shows energy in the DWT coefficients
with decomposition level j.

Now the energy preserving condition can be written as:

(5)
$\left\|\mathbf{X}\right\|{{}^{2}}=\left\|\mathbf{W}\right\|{{}^{2}}=\sum_{j=1}^{J}\left\|\mathbf{W}_{j}\right\|{{}^{2}}+\left\|\mathbf{V}_{J}\right\|{{}^{2}},$

where $\mathbf{X}$ is our VM data (the signal vector magnitude of
accelerometer data; see Sec.3.2) with length N, $j=1,2,...,J$ is the discrete
wavelet transform decomposition level. $\mathbf{W}_{j}$ denotes the detailed
coefficient in scale j, and is a vector of length $N/{2^{j}}$ representing the
differences in adjacent weighted averages from scale 1 to scale J.
$\mathbf{V}_{J}$ denotes the approximated coefficients in the Jth level and
has the same length as $\mathbf{W}_{J}$. Based on the decomposition, each
$\left\|\mathbf{W}_{j}\right\|{{}^{2}}$ represents a special part of the
energy in our VM data which relates to the certain frequency domain (Preece*
et al., 2009) (Donald B. Percival, 2000).

Then the sample variance from (Donald B. Percival, 2000) can be decomposed as:

(6)
$\widehat{\sigma}_{\mathbf{X}}^{2}=\frac{1}{N}\left\|\mathbf{W}\right\|{{}^{2}}-\overline{X}=\sum_{j=1}^{J}\frac{\left\|\mathbf{W}_{j}\right\|{{}^{2}}}{N}.$

The term $\frac{\left\|\textbf{W}_{j}\right\|{{}^{2}}}{N}$ represents the
sample variance (corresponding to $j$ at different scales of DWT
decomposition) in our VM data $\mathbf{X}$.

There are many wavelet features (e.g., (Preece* et al., 2009)) for the
classification of dynamic activities from accelerometer data using DWT. On
this basis, we extract the features from the energy preserving condition and
sample variance mentioned previously.

We aim to look for the features which imply the recovery level among the
stroke patients (see Sec.3.3). Now, we define the features in the j-th level
discrete wavelet transform and discrete wavelet packet transform:

$\mathbf{SSD}_{j}=\frac{\left\|\textbf{W}_{j}\right\|^{2}}{N/{2^{j}}}=2^{j}\frac{\left\|\textbf{W}_{j}\right\|^{2}}{N}.$

For the detailed coefficients $\textbf{W}_{j}$ at decomposition level j,
$\left\|\textbf{W}_{j}\right\|^{2}$ presents its energy and the raw data with
length N. Hence the physical explanation of $\mathbf{SSD}_{j}$ is that it
stands for the point energy at the decomposition level j. Moreover, from the
Eq(6), $\frac{\left\|\textbf{W}_{j}\right\|{{}^{2}}}{N}$ represents the sample
variance at the decomposition level j, $\mathbf{SSD}_{j}$ also has properties
of both the energy preserving condition and the sample variance in wavelet
analysis with constant ${2^{j}}$ .

Comparing with $\mathbf{SSD}_{j}$ (sum of Square value of DWT coefficients at
scale $j$ (with normalisation)), we define other features call
$\mathbf{SAD}_{j}$, which is sum of Absolute value of DWT coefficients at
scale $j$ (with normalisation):

$\mathbf{SAD}_{j}=\frac{\left\|\mathbf{W}_{j}\right\|_{1}}{N/{2^{j}}}=2^{j}\frac{\left\|\mathbf{W}_{j}\right\|_{1}}{N}.$

After we check the correlation between the important wavelet feature
$\mathbf{PNP}$ ( Sec.3.3) and CAHAI score, the branch of features
$\mathbf{PNP}$ using $\mathbf{SAD}$ based perform better than those using
$\mathbf{SSD}$ based in Table 7. Hence we consider the commonly used feature
$\mathbf{SAD}_{j}$ in this paper.

- |  | Acute | Patients |  |  | Chronic | Patients |
---|---|---|---|---|---|---|---|---
Scale (k) |  | $PNP^{1}_{k}$
---
$(SSD)$
| $PNP^{2}_{k}$
---
$(SSD)$
| $PNP^{1}_{k}$
---
$(SAD)$
| $PNP^{2}_{k}$
---
$(SAD)$
| $PNP^{1}_{k}$
---
$(SSD)$
| $PNP^{2}_{k}$
---
$(SSD)$
| $PNP^{1}_{k}$
---
$(SAD)$
| $PNP^{1}_{k}$
---
$(SAD)$
k=1.1 | 0.60 | -0.65 | 0.68 | -0.70 | 0.45 | -0.45 | 0.56 | -0.56
k=1.2 | 0.60 | -0.66 | 0.69 | -0.71 | 0.46 | -0.45 | 0.57 | -0.56
k=1.3 | 0.63 | -0.69 | 0.70 | -0.72 | 0.49 | -0.48 | 0.58 | -0.57
k=1.4 | 0.62 | -0.68 | 0.69 | -0.71 | 0.47 | -0.47 | 0.57 | -0.57
k=2 | 0.65 | -0.69 | 0.69 | -0.71 | 0.45 | -0.45 | 0.56 | -0.55
k=3 | 0.63 | -0.67 | 0.67 | -0.68 | 0.39 | -0.38 | 0.53 | -0.52
k=4 | 0.59 | -0.63 | 0.60 | -0.63 | 0.31 | -0.30 | 0.48 | -0.47
k=5 | 0.46 | -0.50 | 0.49 | -0.52 | 0.29 | -0.27 | 0.43 | -0.42
k=6 | 0.32 | -0.38 | 0.35 | -0.38 | 0.20 | -0.16 | 0.35 | -0.34
k=7 | 0.16 | -0.19 | 0.19 | -0.20 | 0.13 | -0.10 | 0.25 | -0.24
Table 7. The correlation between SAD and SSD based wavelet features and CAHAI
score for acute and chronic patients .

In our analysis, we assume the discrete wavelet decomposition level $J=7$
which is the same level as in (Sekine et al., 1998) and contains enough low-
frequency component as the stroke patients’ movement. The frequency domain
with seven scales is shown in Table 8:

| Scale 7 | Scale 6 | Scale 5
---|---|---|---
Frequency | 0.0078hz-0.0156hz | 0.0156hz - 0.0312hz | 0.0312hz - 0.0625hz
| Scale 4 | Scale 3 | Scale 2
Frequency | 0.0625hz - 0.125hz | 0.125hz - 0.25hz | 0.25hz - 0.50h
| Scale 1 |  |
Frequency | 0.50hz - 1hz |  |
Table 8. The frequency domain from scale 1 to scale 7 by using DWT.

So far, we have decomposed the VM data X to get $\textbf{W}_{1}$,
$\textbf{W}_{2}$, … , $\textbf{W}_{7}$ using DWT. Since the frequency domain
at scale 1 is so wide (0.50hz - 1hz), it is better to divide it into smaller
one, then using DWPT in Appendix 5.1, we can further decompose
$\textbf{W}_{1}$ into $\textbf{W}_{3.4}$, $\textbf{W}_{3.5}$,
$\textbf{W}_{3.6}$ and $\textbf{W}_{3.7}$ which are the results of the $3$-rd
stage of DWPT, each coefficient vector with length $N/{2^{3}}$ has the same
dimension as the coefficients in the third level of DWT decomposition, that is

$\left\|\textbf{X}\right\|{{}^{2}}=\left\|\textbf{W}\right\|{{}^{2}}=\left\|\textbf{W}_{3.4}\right\|{{}^{2}}+\left\|\textbf{W}_{3.5}\right\|{{}^{2}}+\left\|\textbf{W}_{3.6}\right\|{{}^{2}}+\left\|\textbf{W}_{3.7}\right\|{{}^{2}}+\sum_{j=2}^{J}\left\|\textbf{W}_{j}\right\|{{}^{2}}+\left\|\textbf{V}_{J}\right\|{{}^{2}}.$

Now we have coefficients at 10 decomposition scales by using DWT and DWPT:
$\textbf{W}_{3.4}$, $\textbf{W}_{3.5}$, $\textbf{W}_{3.6}$,
$\textbf{W}_{3.7}$, $\textbf{W}_{2}$, $\textbf{W}_{3}$, $\textbf{W}_{4}$,
$\textbf{W}_{5}$, $\textbf{W}_{6}$ and $\textbf{W}_{7}$. Based on these
detailed coefficients, we define the commonly used wavelet features again:

$\begin{split}\textbf{Scale 1.1}:\
&SAD_{1.1}=\frac{\left\|\textbf{W}_{3.4}\right\|_{1}}{N/{2^{3}}}=2^{3}\frac{\left\|\textbf{W}_{3.4}\right\|_{1}}{N},\\\
\textbf{Scale 1.2}:\
&SAD_{1.2}=\frac{\left\|\textbf{W}_{3.5}\right\|_{1}}{N/{2^{3}}}=2^{3}\frac{\left\|\textbf{W}_{3.5}\right\|_{1}}{N},\\\
\textbf{Scale 1.3}:\
&SAD_{1.3}=\frac{\left\|\textbf{W}_{3.6}\right\|_{1}}{N/{2^{3}}}=2^{3}\frac{\left\|\textbf{W}_{3.6}\right\|_{1}}{N},\\\
\textbf{Scale 1.4}:\
&SAD_{1.4}=\frac{\left\|\textbf{W}_{3.7}\right\|_{1}}{N/{2^{3}}}=2^{3}\frac{\left\|\textbf{W}_{3.7}\right\|_{1}}{N},\\\
\textbf{Scale \ j}:\
&SAD_{j}=\frac{\left\|\textbf{W}_{j}\right\|_{1}}{N/{2^{j}}}=2^{j}\frac{\left\|\textbf{W}_{j}\right\|_{1}}{N},\qquad\qquad
j=2,3,4,5,6,7.\\\ \end{split}$

There are 10 features which provide reliable and valid information
(corresponding to more frequency domains) from different frequency domains.
The frequency domain of these features, among 10 scales, is listed in Table 9:

| Scale 1.1 | Scale 1.2 | Scale 1.3
---|---|---|---
Frequency | 0.5hz - 0.625hz | 0.625hz - 0.75hz | 0.75hz - 0.875hz
| Scale 1.4 | Scale 2 | Scale 3
Frequency | 0.875hz - 1hz | 0.25-0.50hz | 0.125hz - 0.25hz
| Scale 4 | Scale 5 | Scale 6
Frequency | 0.0625hz - 0.125hz | 0.0312hz - 0.0625hz | 0.0156hz - 0.0312hz
| Scale 7 |  |
Frequency | 0.0078hz - 0.0156hz |  |
Table 9. The frequency domain from scale 1.1 to scale 7 by using DWPT and DWT.