---
abstract: 'Previous research on automatic pain detection from facial expressions has focused primarily on ‘one-size-fits-all’ metrics (such as PSPI). In this work, we focus on directly estimating each individual’s self-reported visual-analog scale (VAS) pain metric, as this is considered the gold standard for pain measurement. However, the VAS pain score is highly subjective and context-dependent, and its range can vary significantly among different persons. To tackle these, we propose a novel two-stage personalized model, named DeepFaceLIFT, for automatic estimation of VAS. This model is based upon (1) Neural Networks and (2) Gaussian Process Regression learning frameworks and is used to personalize the estimation of self-reported pain via a set of hand-crafted personal features and multi-task learning. We show on the benchmark dataset for pain analysis (The UNBC-McMaster Shoulder Pain Expression Archive) that the proposed personalized model largely outperforms the traditional, unpersonalized models: the intra-class correlation improves from a baseline performance of 19% to a personalized performance of 35% while also providing confidence in the model’s estimates - in contrast to existing models for the target task. Also, DeepFaceLIFT automatically discovers the pain-relevant face regions for each person, allowing for an easy interpretation of the pain-related facial cues.'
author:
- |
    Dianbo Liu\*$^{2,3}$, Fengjiao Peng\*$^{1}$, Andrew Shea\*$^{3}$, Ognjen (Oggi) Rudovic$^{1}$ $\&$ Rosalind Picard$^{1}$\
    \
    \
    \
bibliography:
- 'Mendeley.bib'
- 'MendeleyOggie.bib'
- 'library.bib'
- 'bibfile.bib'
- 'ijcai17.bib'
nocite:
- '[@Lucey2011; @Ashraf2009]'
- '[@Hammal2012]'
- '[@Kaltwang2012]'
- '[@rudovic2013automatic]'
- '[@Rusk2016]'
- '[@Sikka2013]'
title: |
    DeepFaceLIFT: Interpretable Models\
    for Personalized Automatic Estimation of Self-reported Pain
---

Introduction
============

Research on automatic estimation of pain levels from facial expressions has focused primarily on estimating (arguably) objective metrics rather than self-reports. These objective metrics, such as Prkachin and Solomon Pain Intensity (PSPI) [@Prkachin2008], use facial action intensities to quantify facial expressions and produce pain scores. Using these metrics, two individuals who make the same facial expression should be given very similar pain scores. However, because pain is a fundamentally subjective experience [@textbookPain21] and because individuals express pain differently depending on qualitative factors such as their motor ability, gender, and age [@fillingim2009sex; @kunz2008impact; @scipio2011perception] these objective metrics may not sufficiently capture an individual’s true pain. Furthermore, many studies have found low correlations between facial expressions and self-reported pain levels [@Prkachin2008]. These observations pose the question of how well automatic estimates of objective metrics capture an individual’s actual experience. Self-reported pain scores are often considered the *gold standard* because they provide information about an individual’s subjective experience [@craig1992facial]. Therefore, an approach that automatically estimates self-reported pain scores rather than objective pain metrics may be more useful for individuals experiencing pain.

![Graphical representation of the proposed DeepFaceLIFT two-stage learning for estimation of VAS ($y_p$) from facial landmarks ($x$) extracted from each frame in the sequence. The neural network model is trained on per-frame data and outputs a per-frame estimation of VAS, which is concatenated into a sequence $\mathbf{V}$. The sequence $\mathbf{V}$ is next used to compute feature statistics $\mathbf{S}$, which is then used to train a Gaussian Process regression model to estimate per-sequence VAS.[]{data-label="fig:twostage"}](figures/two-stage-learning.png){width="50.00000%"}

Because individuals make different facial expressions depending on personal factors (e.g. their surroundings, age, and so on), a method for automatic pain estimation would likely improve by not only accounting for how individual facial expressions and specific facial actions vary, but also personalizing on qualitative characteristics by explicitly representing them as features. These features could then be incorporated into a model. This form of feature-level personalization differs from the model-level personalization (i.e., varying the model architecture) as a single model is consistently applied to all subjects. Using the personalized features may enable the model to learn additional personal information thus increasing its estimative power.

Another important aspect of personalized applications is ensuring that model outputs can be easily traced back to specific components of the input features, as doing so is crucial for explaining the specific characteristic and behavioral differences between individuals. This problem is broadly faced in deep learning research due to the lack of interpretability of learned representations with respect to the original input features [@zeiler2014visualizing; @shrikumar2016not]. For example, some methods have shown that personalizing models improve outputs [@zen2014unsupervised; @clifton2013gaussian]; however, the model outputs and learned internal representations cannot be easily nor intuitively linked to specific portions of the original inputs. As a result, while the improved performance achieved through personalization indicates that qualitative differences indeed exist between training subjects, exactly what these differences are and how they affect the final output are is not clear. Therefore, having the ability to examine how specific aspects of input features contribute to output estimations would provide insight into what differences exist amongst subjects and how these differences influence model estimations.

To address these challenges of subjectivity, personalization, and interpretability, we propose a two-stage hierarchical learning algorithm called DeepFaceLIFT (LIFT = Learning Important Features) [@shrikumar2016not] that estimates pain scores obtained using the self-reported Visual Analog Scale (VAS) [@carlsson1983assessment]. The first stage of the model comprises a fully-connected neural network (NN) that takes raw Active Appearance Model (AAM) landmark coordinates [@unbc] as input features. The second stage then processes the NN outputs and passes them as features into a Gaussian process (GP) model. In addition to producing accurate VAS estimates, the model also provides an interpretable learning framework by leveraging a model, called DeepLIFT [@shrikumar2016not], in the first stage to examine inputs and utilizing an RBF-ARD kernel in the second stage to examine the learned GP kernel weights. We show on the UNBC-McMaster Shoulder Pain Expression Archive Database [@unbc] that the proposed personalized approach for automatic VAS estimation outperforms traditional approaches. The model outline is shown in Fig.\[fig:twostage\].

Related Work
============

Although much research exists on automated recognition of affect from human facial expressions (for surveys, see [@Zeng2009; @tac_ref58]), until recently, only a handful of works have focused on automated pain estimation. Due to advances in computer vision and the recent release of the UNBC-McMaster dataset [@unbc], pain analysis from face images has seen significant advances. This dataset provides videos where each frame has been rated using the PSPI score. This commonly used pain score is considered to be an objective pain score, in contrast with subjective ones such as VAS. Although VAS is still the most commonly accepted pain score in clinical settings, all existing automatic methods for pain estimation from pain images have focused on estimating PSPI scores. We outline below the recently published works.

Face shape features have been used in NNs to classify images of faces showing typical mood versus pain expressions [@monwar2006pain]. Lucey *et al.*, 2011a, and Ashraf *et al.*, 2009, used AAM-based features combined with Support Vector Machine (SVM) classifiers to classify pain versus no-pain images. Rather than treating pain as binary classification task, Hammal and Cohn, 2012, attempted to estimate pain intensity on a 4 level scale using one-versus-all SVM classifiers, and Kaltwang *et al.*, 2012, performed estimation of the full 15 level PSPI scale using Relevance Vector Regression models. Likewise, Rudovic *et al.*, 2013, proposed a Conditional Random Field (CRF) model with heteroscedastic variance, showing that the changing model variance provides a better fit of individual pain sequences. More recently, several works using deep learning have been proposed for PSPI estimation [@egede2017fusing; @Rodriguez2017]. Other works have attempted to estimate the pain intensity of entire sequences. Ruiz *et al.*, 2016, proposed a HCRF framework and semi-supervised learning method to identify the peak of pain in image sequences. Sikka *et al.*, 2013, proposed a multi-instance learning framework to identify the most expressive (in terms of pain) segments within image sequences. However, these per sequence labels are derived using heuristics and do not relate to established pain ratings such as VAS.

Aside from the UNBC-McMaster dataset [@unbc], there have been a variety of pain recognition studies based on other datasets (e.g., see [@werner2012pain; @tac_ref66]). Pain detection has been attempted using physiological signals such as brain hemodynamic responses using NIRS [@Aasted2016; @Yucel2015]. Nevertheless, these works have neither attempted automatic estimation of VAS nor in a personalized manner.

DeepFaceLIFT: The Method {#pmodel}
========================

For clarity, we will first briefly describe the dataset used for personalized modeling of the self-reported pain score, VAS. We used the UNBC-McMaster Shoulder Pain Expression Archive Database (UNBC-PAIN) [@unbc], a publicly available dataset containing facial video of patients who suffer from one-sided shoulder pain undergoing various arm range-of-motion exercises. The dataset contains 200 image frame sequences gathered from 25 subjects, totaling 48,398 image frames. Each frame is coded in terms of the facial action unit (AU) intensity on a six-point ordinal scale. The PSPI rating is computed as a non-linear combination of intensities of 6 AUs within each image as: $$PSPI = AU4 + max(AU6,AU7) + max(AU9,AU10) + AU43$$

For each image sequence, the dataset also includes observed pain index (OPI) scores, rated by experienced pain-observers on a Likert-like scale (0-5) [@unbc]. To evaluate the individual’s perception of pain, we do not use PSPI in the training or testing process. However, in certain settings, we do use the OPI scores. The dataset also includes for each image sequence the self-reported VAS ratings defined on an ordinal scale of 0-10. Note that although the UNBC-PAIN dataset contains a large amount of image frames, estimating VAS scores poses a challenge due to the limited number of sequence-level labels (only 200 sequences).

To personalize the models, we obtain three additional manually-labeled features derived from each subject’s appearance: complexion, age, and gender. Complexion is divided into three categories using a Fitzpatrick-like skin-tone scale: pale-fair, fair-olive, and olive-dark [@fitzpatrick1988validity]; age is also divided into three bins: young, middle-aged, and elderly; and gender is divided into two bins: male and female. These bins were independently assigned and then agreed upon by the primary authors. These personal features were selected based on their straightforward labeling and their relevance to pain expression. For example, individuals of different age and gender have been shown to express pain differently. Complexion is used because it is easily observed and thought to loosely contain information about an individual’s medical information, background, and/or daily habits.

AAM facial landmarks are used as generic input features instead of raw images. We do so for two reasons. The first is to reduce the feature space of each individual from 320x240 (pixels) to 132 (66 landmark coordinates). The second and primary reason is to ensure that the feature representations are distinct, easily-identifiable, and intuitive. Furthermore, doing so allows us to apply DeepLIFT to interpret how individual facial landmarks vary across persons and as a function of reported VAS scores. Note also that various methods have been reported that capture temporal information in a sequence of image frames, such as Recurrent Neural Networks [@mikolov2010recurrent] and Conditional Random Fields [@rudovic2013automatic]. Instead of using a temporal model, we consider the problem as a static one. This is accomplished by extracting sequence-level statistics from the outputs of the first stage, then passing these statistics into the GP regression framework in the second stage. In this way, the model encodes and captures both frame-level and sequence-level dependencies.

[**Notation.**]{} The estimation of VAS is formulated as a sequence-level regression problem, where we are given $N_i$ image sequences of person $i\in \{0,\dots,P\}$ with P being the number of target persons. The sequences of each person are annotated in terms of VAS as $\mathcal{V}=\{V_{1},\dots,V_{P}\}$, where $V_{i}=\{v_{i}^{1},\dots,v_{i}^{N_i}\}$ and individual per-sequence VAS scores are $v_i\in \{0,...,10\}$. For OPI scores we have: $\mathcal{O}=\{O_{1},\dots,O_{P}\}$, where $O_{i}=\{o_{i}^{1},\dots,o_{i}^{N_i}\}$ and $o_i\in \{0,...,5\}$. Additionally, let $f \in \{1, \dots, \mathcal{F}\}$ denote the frame number and $s \in \{1, \dots, \mathcal{S}\}$ denote the sequence number. Each sequence is represented by a set of input features $X_s=\{x_1, \dots,x_{L_s}\}$, where $L_s$ is duration or the number of frames in a given sequence $s$.

Weakly Supervised Neural Network {#stageone}
--------------------------------

In order to leverage all frame-level data before making the sequence-level estimations, the first stage of the DeepFaceLIFT model applies multi-instance learning to train a fully-connected NN (Fig.\[fig:twostage\]). This network uses frame-level AAM facial landmark features as input, where each feature is labeled with the VAS score of its corresponding sequence. Applying a NN is appropriate in this setting as the facial landmarks may be related in highly nonlinear fashions, and these relations can be learned automatically using NNs. A 4-layer NN with ReLU activation functions is used in the first stage of training to generate frame level VAS scores. The same hidden-layer architecture is used in all experiments, where each hidden layer in the NN contains 300, 100, 10 and 100 nodes, respectively. The NN is trained for 100 epochs with a batch size of 300.

The model generates frame-level VAS scores in 3 separate settings. Each setting differs in how personal information is used in the model (Not used, Appended to the 3rd layer of NN, Appended to Input Features). Table 1 shows these settings in the first column labeled with S1. Within each of these settings, another 2 separate labeling schemes are used to train the NN (Labeled with VAS, Labeled with VAS and OPI). Table 1 shows these settings in the second column labeled with S1. Therefore, for each frame $f \in \{1, \dots, M\}$, where M is the total number of frames in the dataset, the NN estimates either $V_f$ or both $V_f$ and $O_f$. These various settings are used to compare different personalization strategies.

GP with RBF-ARD Kernel {#stagetwo}
----------------------

Using the per frame VAS estimations, DeepFaceLIFT computes sequence level statistics (as described below), which are then fed as inputs into the GP framework to obtain the final VAS estimate. The GP framework is tested in 2 different settings. In the first setting, the input features only use the NN-estimated VAS. In the second, both NN-estimated VAS and OPI statistics are used to estimate VAS. Table 1 shows these settings the first column labeled with S2.

Because the target sequences vary in their length, for each sequence $s$, the model computes a set of statistics, resulting in a fixed-size feature vector ${\boldsymbol{S_s}}=$ {$\mu(V)=$ mean, $\eta(V)=$ median, $Min(V)=$ minimum, $Max(V)=$ maximum, $Var(V)=$ variance, $\mu^3(V)=$ 3rd moment, $\mu^4(V)=$ 4th moment, $\mu^5(V)=$ 5th moment, $Sum(V)=$ sum of all values, $IQR(V)=$ interquartile range}. These *sufficient statistics*  [@hogg1995introduction] are chosen because they capture diverse and important information that may be used to infer unknown parameters in many distribution families [@hogg1995introduction]. When using only the NN-estimated VAS to compute these statistics, the length of ${\boldsymbol{S_s}}$ is 10; when using NN-estimated VAS and OPI, the length of ${\boldsymbol{S_s}}$ is 20. ${\boldsymbol{S_s}}$ is then fed as input to the GP regression model [@rasmussen2006gaussian], which is reformulated as: $${\boldsymbol{y}} = f({\boldsymbol{S}}) + \epsilon$$ where $\epsilon \sim \mathcal{N}(0, \sigma^2_{v})$ is i.i.d. additive Gaussian noise. The objective is to infer the latent functions $f$ given the training dataset $\mathcal{D} = \{ {\boldsymbol{S}}, {\boldsymbol{Y}} \}$. Following the framework of GPs, we place a prior on the function $f$ so that the function values ${\boldsymbol{f}} = f({\boldsymbol{S}})$ follow a Gaussian distribution $p({\boldsymbol{F}}|{\boldsymbol{S}}) = \mathcal{N}({\boldsymbol{F}}|{\boldsymbol{0}}, {\boldsymbol{K}})$. Here, ${\boldsymbol{F}} = \{{\boldsymbol{f}}\}$, and ${\boldsymbol{K}} = k({\boldsymbol{S}},{\boldsymbol{S}}^T)$ is the kernel covariance function. Training of the GP consists of finding the hyper-parameters that maximize the log-marginal likelihood [$$\begin{aligned}
\small
\nonumber\log p({\boldsymbol{Y}}|{\boldsymbol{X}}, {\boldsymbol{\theta}}) =&-\textrm{tr}\left[({\boldsymbol{K}} + \sigma^2{\boldsymbol{I}})^{-1}{\boldsymbol{Y}}{{\boldsymbol{Y}}}^T\right] \\
\label{marg}  &- \log\vert {\boldsymbol{K}} + \sigma_v^2{\boldsymbol{I}} \vert + \textrm{const}.\end{aligned}$$ ]{} Given a test input ${\boldsymbol{s}}_\ast$ we obtain the GP estimative distribution by conditioning on the training data $\mathcal{D}$ as $p({\boldsymbol{f}}_\ast |{\boldsymbol{s}}_\ast, \mathcal{D} = \mathcal{N}(\mu({\boldsymbol{s}}_\ast, V({\boldsymbol{s}}_\ast)))$ with the mean and variance given by: $$\begin{aligned}
\small
\label{post_mu}\mu({\boldsymbol{s}}_\ast) &= {{\boldsymbol{k}}_\ast}^T ({\boldsymbol{K}} + \sigma_v^2{\boldsymbol{I}})^{-1}{\boldsymbol{Y}}\\ 
\label{post_s}V({\boldsymbol{s}}_\ast) &= k_{\ast\ast} - 
{{\boldsymbol{k}}_\ast}^T ({\boldsymbol{K}} + \sigma_v^2{\boldsymbol{I}})^{-1}
{\boldsymbol{k}}_\ast,\end{aligned}$$ where [${\boldsymbol{k}}_\ast = k({\boldsymbol{S}}, {\boldsymbol{s}}_\ast)$]{} and [$k_{\ast\ast} = k({\boldsymbol{s}}_\ast, {\boldsymbol{s}}_\ast)$]{}. For convenience, we denote ${\boldsymbol{\mu}}_\ast^{(v)} = \mu({\boldsymbol{s}}_\ast)$ and $V^{(v)}_{\ast\ast} = V({\boldsymbol{s}}_\ast)$.

The Radial Basis Function (RBF) kernel with Automatic Relevance Determination (ARD) [@rasmussen2006gaussian] is defined as: $$k({\boldsymbol{s}},{\boldsymbol{s}}^\prime) = \sigma_f^2\exp(-({\boldsymbol{s}} - {\boldsymbol{s}}^\prime)M({\boldsymbol{s}} - {\boldsymbol{s}}^\prime)),$$ where $M=diag(\ell_1,...,\ell_D)^{-2}$ is a diagonal matrix where each element reveals the importance of each statistic feature ($1,...,D$). The inverse of the length-scales ($\ell$) determines how relevant an input is to the output: at the extremes, if the length-scale is very large, then the covariance will be close to independent of that input, effectively removing that feature from the inference; if the ength scale is very small, then this indicates that the feature is of high relevance [@rasmussen2006gaussian]. The regression mapping is then fully defined by the set of hyper-parameters ${\boldsymbol{\theta}} = \{\ell_1,...,\ell_D, \sigma_f, \sigma_v\}$.

GP-ARD is used for two reasons. The first reason is that due to the small number of target sequences (=200), a nonparametric probabilistic modeling with GP is the best choice due to its ability to generalize from small data. The second is that GP-ARD provides uncertainty in its estimates, which is important for assessing the reliability of the estimated VAS.

To summarize, the learning and inference in this two-stage approach proceeds as follows: in the first stage, DeepFaceLIFT trains a NN multitask learning under different personalization strategies Sec.\[stageone\]. Next, the output of the trained NN is used to compute sequence-level statistics, which are then fed into a GP-ARD model to provide a personalized esimate of VAS (Sec.\[stagetwo\]). Inference in this approach is straightforward: features are propogated through the two-stage model, resulting in the VAS estimate for the test person, along with the estimated uncertainty (as can be seen in Fig.\[fig:Results\]). The per-subject estimation results of VAS and OPI are also shown in Fig.\[fig:Results\].

![ DeepFaceLIFT results for estimation of VAS ($y_p$) from AAM ($x_f$). MAE per subject for estimation of: A) VAS (average MAE=2.18), and B) OPI (average MAE=1.12). Note that the final MAE for DeepFaceLIFT are lower than the NN-MV (Mean Voting) baseline on 18/25 subjects. C-F) Estimated VAS vs. reported VAS of four randomly chosen individuals along with the uncertainty (shaded area) estimate produced by the GP framework. The model learns different uncertainty levels for different individuals. This may be due to various reasons, such as input-feature noise or subjectivity in self-reported VAS scores.[]{data-label="fig:Results"}](figures/ResultFigures.png){width="50.00000%"}

  -- ----- --------- ----------------- ------------- ------------ ------------- --
                                                                                
                                           NN-MV         RNN          HCRF      
            VAS,OPI     2.32(0.25)                                              
            **VAS**   **2.24 (0.27)**                                           
      VAS     VAS       2.30 (0.26)     2.82 (0.05)   3.01(0.05)   3.52 (0.12)  
            VAS,OPI     2.34 (0.27)                                             
            **VAS**   **2.18 (0.35)**                                           
      VAS     VAS       2.24 (0.25)     2.65 (0.17)   3.01(0.05)   3.67 (0.13)  
            VAS,OPI     2.41 (0.23)                                             
            **VAS**   **2.22 (0.27)**                                           
      VAS     VAS       2.22 (0.30)     2.53 (0.18)   3.01(0.05)   3.67 (0.13)  
  -- ----- --------- ----------------- ------------- ------------ ------------- --

Experiments
===========

[**Features and Evaluation Scores.**]{} 66 facial AAM landmarks provided by the database creators [@unbc] are used as input features. Each facial landmark has 2 dimensions (x,y), resulting in 132 AAM features when concatenated. For evaluation, Mean Absolute Error (MAE) is used since it is commonly used for ordinal estimation tasks [@corf]. Intra-Class Correlation ICC(3,1) [@Shrout1979] is also reported, as it is commonly used in behavioral sciences to measure agreement between annotators (in this case, the estimated and true VAS). Additionally, ICC prevents potentially deceptive MAE scores. For example, while MAE provides a good indication of accuracy, it does not capture the correlation between the model estimates and the true values. If a model estimates the mean VAS (3.25) for all the image sequences, then the resulting MAE is 2.44, which is lower than some of the baseline models. In this scenario, ICC would more strongly penalize these errors and capture the agreement between the estimated and true VAS, so ICC would be approximately zero.

[**Learning.**]{} The data was split into 5 folds, each fold containing the data for 5 persons. In each trial, 4 folds were used for training while 1 was held out for testing. The VAS estimates for each fold are concatenated and after all 5 trials, the estimated scores for all 200 sequences are used to compute the MAE and ICC. We compare the results of the two-stage DeepFaceLIFT to a baseline model RNN [@mikolov2010recurrent] comprising a 36-node LSTM layer and one dense layer and an HCRF model [@gunawardana2005hidden] comprising 5 hidden states. Both HCRF and LSTM are common choices for modeling sequential data. We also construct a NN baseline (NN-MV) by taking the mean voting of all the per-frame estimations of a sequence.[^1]

The test results are reported in Table \[my-label\]. When using VAS-derived statistics, DeepFaceLIFT achieves the lowest MAE, 2.18, compared to the other baseline result of the NN-MV (2.48), RNN (3.21), and HCRF (3.38). Regardless the setting, DeepFaceLIFT yields the best MAE and ICC scores compared to the baseline models. The second stage learning improves the ICC results from 0.19 (NN-MV) to 0.35, which demonstrates the effectiveness of the GP regression and choices of the feature statistics.

[**Personal features.**]{} Appending personal features to the NN yields better results than not including them (non-personalized model), as improvement is seen not only in DeepFaceLIFT, but also in RNN and HCRF. Personal features appended in the third layer in the NN reduces the MAE further, more than when appended to the inputs. This improvement suggests that the effect of personal features appended to the inputs are partially “washed out" by the larger (=132) input feature dimension. The discussion of the individual effects of the personal features are included in Sec.4.1 \[pf\].

[**Multi-task learning.**]{} When OPI is used as a second label in the naive NN-MV estimation, MAE and ICC improve upon using using VAS alone. This suggests that OPI can effectively control the personal-variant interpretations of pain to produce a better estimations of VAS.

DeepLIFT: Interpretable Learning {#exp}
--------------------------------

As mentioned earlier, NNs can be cumbersome due to their blackbox nature. There are several approaches available that help researchers infer how each input feature contributes to the final output [@zeiler2014visualizing; @sundararajan2016gradients]. However, most of these methods encounter issues, such as saturation and thresholding, and are computationally expensive. DeepLIFT is robust because it estimates input contributions while remaining computationally efficient and avoiding gradient discontinuities [@shrikumar2016not].

![image](figures/FigureDeepLIFT.png){width="70.00000%"}

\[fig:DeepLIFT\]

DeepLIFT estimates these contributions by performing backpropagation through each neuron [@shrikumar2016not], then calculating the difference between each neuron’s computed activation with its respective “reference activation" [@shrikumar2016not]. These differences taken between specific input and output activations allow DeepLIFT to explain variations in input feature contributions [@shrikumar2016not]. For clarity, let $V_{s}^{i}=\{v_{s}^{i},\dots,v_{s}^{i}\}$ and $X=\{x_1, \dots,x_T\}$ represent the output and input, respectively. For each person, features are first normalized using Z-scores. Next, DeepLIFT uses a zero-vector $X_{0}=\{x_{r_1}, \dots,x_{r_{T}}\}$ as a reference activation for the input. $X_0$ is then passed as an input into the NN to obtain the output $f(X_{0})$, which is then used as the reference activation for the output. Each input is assigned a contribution score $C_{\Delta x_{f} \Delta V_{s}^{i}}$ such that $\Delta V_{s}^{i}$=$\sum C_{\Delta x_{f} \Delta V_{s}^{i}}$. For NNs with multiple layers, $\Delta V_{s}^{i}$ is backpropagated to the input layer using the standard chain rule as follows: if $x$ is an input and $h$ is a neuron in a hidden layer, a multiplier that is used to calculate the contribution of $\Delta x$ to $\Delta y$ is defined as: $m_{\Delta x \Delta y} = \frac{C_{\Delta x \Delta h}}{\Delta x}$. Assuming $y$ is the target output and only $h$ is between $x$ and $y$, then $m_{\Delta x \Delta y} = \sum{m_{\Delta x \Delta h}m_{\Delta h \Delta y} }$.

The contribution of each landmark toward the frame-level VAS estimation is shown in Fig.\[fig:DeepLIFT\]A, where darker colors indicate higher contribution or importance towards the final output. In this figure, an average face is computed using all image frames to illustrate the locations of each AAM landmark. The geometric coordinates of the top 5 AAM landmarks are indicated as well. Note that the DeepLIFT model finds more importance in one subset of landmarks when estimating VAS and another subset when estimating OPI, with the minimum overlap of the top five landmarks. This difference may have occured for various reasons. One possibility is that the model may have simply model converged on different sets of parameters. Another possibility comes from the different ways in which VAS and OPI scores are obtained, i.e. with OPI, external observers may perceive a subject’s pain differently from the subject’s self-report while focusing on different facial regions when making their assessment. Additionally, the variation of importance of these landmarks among different subjects is shown alongside their respective figures. Fig. \[fig:DeepLIFT\]B. shows the standard deviations of importance scores across different subjects along with the boxplots for the landmarks with the top 5 standard deviations. Fig. \[fig:DeepLIFT\]C and D show the importance scores and their standards deviations for the final OPI estimates. Examples of four different subjects are given in Fig.\[fig:DeepLIFT\]E, F and G to illustrate the relative importance of VAS and OPI statistics and are accompanied by pictures of the subjects. These results suggest that landmarks contribute differently to VAS and OPI estimation and that the importance of each landmarks varies among subjects. Based on Fig.\[fig:DeepLIFT\], there appears to be a slight increase in the importance of landmarks on the left side of the face. This may be attributed to the possible redundancy of AAM landmarks due to facial symmetry, as well as to how the camera is aimed at the faces.

[**Personal Features.**]{} The importance of each of the personal features for the pain estimation is approximated by calculating the MAE of outputs in Stage 1 while excluding a single personal feature. The importance of OPI in Stage 1 is assessed by excluding it from the label and comparing the differences in the resulting MAEs obtained with and without the OPI label. As shown in Fig.\[fig:PFandWeights\]A, age contributes the most compared to the other personal features, while all of the features improve VAS estimation. These findings are consistent with the clinical studies which indicate that age, gender, and other personal features may affect pain perception and expression. Additionally, OPI labels appear to contain information that improves VAS estimates.

![ A) The contribution of personal features to the VAS estimation. The Y axis indicates the MAE increases (decrease in performance) when the target personal feature is excluded from training. B) The aggregated results of relevance of the feature statistics from Sec.\[stagetwo\] computed using the RBF-ARD kernel in the GP framework. A smaller or more negative value indicates that the statistic has more importance (In this figure, -ln($l_i)$, $i\in\{1\dots 10\}$, of the model parameters is shown). The above figure shows that the sum, mean, and IQR of the frame-level estimates are the three most important statistics.[]{data-label="fig:PFandWeights"}](figures/PFandWeights.png){width="48.00000%"}

[**Statistics.**]{} We ran a 5-fold cross validation to aggregate the relevance of the feature statistics. The results are shown in Fig.\[fig:PFandWeights\]B. The smaller the numerical kernel value is, the more that statistic contributes to the final VAS estimate. According to the model, the top 3 most important statistics are $\Sigma (P),\mu (P),$ and $IQR(P)$, with $\Sigma (P)$ being the most important. This finding suggests that the average result of the NN and the length of the image sequence both contribute greatly to the estimation of VAS, as $\mu(P) = \frac{\Sigma (P)}{L}$, where L is the length of the sequence. The fact that $L$ increases the importance of the statistic indicates a correlation between the length of the pain-generating experience and VAS score.\

Discussion and Conclusions
==========================

We presented a two-stage interpretable machine learning model to automatically estimate subjective pain level from face images. We demonstrate that a weakly supervised learning model is able to estimate VAS in each video sequence accurately, and the inclusion of personal features and objective pain level (OPI) during training improves model performance. Applying a GP framework to estimate sequence-level VAS scores using statistics of frame-level embeddings further improves the performance. Another important contribution of this paper is that the importance of features in every stage of the learning can be tracked under the framework of DeepFaceLIFT. This makes it possible to explore the relations between input features and performance of the algorithm, improving the understanding of how humans perceive pain. Importance scores estimated by DeepFaceLIFT suggest that subjective pain and objective pain are perceived in different ways and the importance of facial landmarks for pain detection varies significantly among individuals. However, more analysis is needed to verify whether these differences reflect the model’s behavior solely, or are also a consequence of the raters’ differences (VAS vs OPI). In addition, important personal features and statistics from embeddings were identified. We showed that both of these are important for improving the models performance. We also showed that utilizing personal features at different levels in the model results in different performance gains, with the optimal performance achieved when they are embedded at the intermediate layers of the NN. Furthermore, by exploiting the benefits of multi-task learning, we obtain a more robust estimation of VAS [@caruana1998multitask]. We attribute this to the fact that OPI, coming from an external observer, acts as a strong reference point when estimating the VAS scores, which are highly subjective.

To summarize, we proposed a personalized machine learning approach for estimation of VAS. We illustrated the personalization process using several different settings, and in the context of NNs and GP modeling frameworks. The results show clear benefits of the personalization. In future works, we plan to further investigate the relationships between VAS and other pain scores (such as PSPI), and its relationships to facial action units. As context is another important aspect not explored in this paper, we plan to derive more advanced statistics that would capture additional information that can enhance estimation of target VAS. Finally, we hope that these findings will advance applications of pain estimation in clinical settings.

Acknowledgements
================

The work of O. Rudovic is funded by European Union H2020, Marie Curie Action - Individual Fellowship no. 701236 (EngageMe). The work of D. Liu is funded by the Vest Scholarship and Welcome Trust Fellowship. We also thank Daniel Lopez Martinez (Media Lab, MIT) and Kelly N Peterson (MIT) for their help in the early stages of the project.

[^1]: We have also attempted using other statistics, such as the median, max, etc., to construct the baseline neural network model. The mean voting provides the best estimates and is thus used as the baseline.
