# A Nonparametric Bayesian Item Response Modeling Approach for Clustering
Items and Individuals Simultaneously

Guanyu Hu Zhihua Ma Insu Paek Department of Statistics, University of
Connecticut; Department of Statistics, University of Missouri
ColumbiaCorresponding Author, Email<EMAIL_ADDRESS>College of
Economics, Shenzhen UniversityDepartment of Educational Psychology & Learning
Systems, Florida State University

###### Abstract

Item response theory (IRT) is a popular modeling paradigm for measuring
subject latent traits and item properties according to discrete responses in
tests or questionnaires. There are very limited discussions on heterogeneity
pattern detection for both items and individuals. In this paper, we introduce
a nonparametric Bayesian approach for clustering items and individuals
simultaneously under the Rasch model. Specifically, our proposed method is
based on the mixture of finite mixtures (MFM) model. MFM obtains the number of
clusters and the clustering configurations for both items and individuals
simultaneously. The performance of parameters estimation and parameters
clustering under the MFM Rasch model is evaluated by simulation studies, and a
real date set is applied to illustrate the MFM Rasch modeling.

Keywords: Bayesian Inference, Heterogeneity Pursuit, Mixture of Finite
Mixtures, Nonparametric Bayesian IRT, Rasch Model

##  1 Introduction

In psychometrics, item response theory (IRT) is a popular modeling paradigm
for measuring subject latent traits and item properties according to discrete
responses in tests or questionnaires. In IRT models, the discrete responses
are explained as a function of respondents’ latent traits and item properties.
The Rasch model (Rasch,, 1980), one of the most widely used IRT models, can be
used to infer the latent abilities of the test-takers as well as the item
difficulties. Generally, an assumption of most IRT models including the Rasch
model is measurement invariance, which means all items measure the latent
trait in the same way for all subjects, that is, the item response function
for an item is invariant across all subjects. However, this assumption may be
violated in some cases by poorly-constructed items or items with varying
difficulty levels for different subgroups of persons.

Heterogeneity in IRT analysis is typically addressed by latent class analysis
(LCA) approaches and mixture IRT models. The LCA model assumes that the sample
of respondents is drawn from a population divided into several latent classes,
with individuals in the same class sharing the same distribution of the
response variables (Lazarsfeld and Henry,, 1968). LCA model attempts to
identify distinct classes of individuals and classify the test-takers. For an
overview of LCA methods and applications, see e.g. Hagenaars and McCutcheon,
(2002). The mixture IRT models are regarded as a combination of IRT models and
LCA approaches, whose basic idea is that the same IRT model holds within
subjects in the same class, whereas different IRT models hold among different
classes (Gnaldi et al.,, 2016). Thus, with mixture IRT models, it is possible
to identify distinct subpopulations within a larger population, each of which
responds differently to a set of items (Rost,, 1990; Rost and Langeheine,,
1997).

Under both the LCA and mixture IRT frameworks, however, the number of clusters
should be pre-specified by the user before obtaining the estimates. In some
applications, the number of clusters is chosen by prior knowledge, which is
not always available. Another approach for selecting the number of clusters is
based on information criteria such as the Akaike information criterion (AIC;
Akaike,, 1998; Bartolucci,, 2007), and Bayesian information criterion (BIC;
Schwarz,, 1978), which are computationally ineffective and time consuming (Pan
and Huang,, 2014). Moreover, the two-stage approach that inferences
conditional on the specification of the number of clusters ignores the
uncertainty in the selection process (Yang et al.,, 2011).

This problem is circumvented by nonparametric Bayesian approaches.
Nonparametric Bayesian approaches allow joint inferences on the clustering
information, which includes both the number of clusters, clustering
configurations, and parameters on individuals and items simultaneously. For
example, Miyazaki and Hoshino, (2009) proposed a Bayesian nonparametric
mixture IRT model using a Dirichlet process prior. Liu, (2019) developed a
Bayesian nonparametric ordered latent class model with a modified Chinese
restaurant process prior. Although Bartolucci et al., (2017) used a Bayesian
nonparametric method on IRT model with a LCA formulation, the number of
clusters was still determined via model selection criteria. Research on
Bayesian nonparametric approaches for dealing with clustering in IRT models
are still limited, and the existing approaches mainly focus on Dirichlet
process prior.

The aim of this paper is to propose a new hierarchical Rasch model for
simultaneously clustering items and individuals under a Bayesian framework.
Bayesian nonparametric approaches offer choices to allow uncertainty in the
number of clusters, and provides an integrated probabilistic framework under
which the number of clusters, the clustering configuration, subject latent
abilities and item difficulties are simultaneously estimated. For the
clustering problem, Bayesian inference provides a probabilistic framework for
simultaneous inference of the number of clusters and the clustering
configurations based on a Dirichlet process mixture model (DPM; Ferguson,,
1973). Miller and Harrison, (2013) points out that DPM will produce extremely
small clusters and make the estimation for number of clusters inconsistent.
The mixture of finite mixtures (MFM) model proposed by Miller and Harrison,
(2018) provides a remedy to over-clustering problem for Bayesian nonparametric
methods. MFM is widely applied in different areas such as regional economics
(Hu et al.,, 2020), social science (Geng et al., 2019a, ) and environmental
science (Geng et al., 2019b, ). Thus, our key idea is to formulate the basic
Rasch model into a hierarchical form, and assign the MFM priors for the latent
ability and item difficulty parameters. By assigning the MFM priors, our
proposed approach allows for the uncertainty in the number of clusters, and
through a Bayesian framework, estimation is obtained for both the number of
clusters, clustering configurations and parameters in terms of person latent
traits and item difficulty parameters. Hence, the proposed method can be
effectively used to identify item and individual clustering. Test developers
may use the clustering results of item difficulties for their test
construction. In addition, the information on the clustering of subject latent
abilities could be useful when grouping of individuals is needed for placement
or other diagnostic assessment purposes.

The remainder of the paper is organized as follows. In Section 2, we present
the proposed hierarchical Rasch model with MFM. In Section 3, the Bayesian
inference procedure including prior and posterior distributions and model
comparison criterion are discussed. Simulation studies and the results are
described in Section 4. For the purpose of illustration, our proposed
methodology is applied to an English test data in Section 5. Finally, we
conclude this paper with a discussion in Section 6.

##  2 Methodology

###  2.1 Item Response Model

Suppose there are a total of $N$ subjects and $J$ items, and
$\bm{y}_{i}=(y_{i1},\cdots,y_{iJ})^{\prime}$ represents the binary response
vector for the $i$th individual at $J$ items. $y_{ij}=1$ means correct answer
at the $j$th item for individual $i$, while $y_{ij}=0$ represents wrong
answer. Denoting $P(y_{ij}=1|\theta_{i})=p_{ij}$, the Rasch model (Rasch,,
1980) is given as

$p_{ij}=\frac{\text{exp}(\theta_{i}-b_{j})}{1+\text{exp}(\theta_{i}-b_{j})},$
(1)

for $i=1,\cdots,N$, $j=1,\cdots,J$, where
$\bm{\theta}=(\theta_{1},\cdots,\theta_{N})$ is the vector of latent
abilities, and $\bm{b}=(b_{1},\cdots,b_{J})$ denotes the vector of item
difficulty parameters.

Under a Bayesian framework, the Rasch model in (1) is written as

$\begin{split}y_{ij}&\sim\text{Bernoulli}(p_{ij}),\\\
\text{log}(p_{ij}/(1-p_{ij}))&=\text{logit}(p_{ij})=\theta_{i}-b_{j}.\end{split}$
(2)

By giving the prior distributions of the unknown parameters and through the
Markov chain Mote Carlo (MCMC) sampling algorithm, the posterior estimates of
parameters for $i=1,\cdots,N$, $j=1,\cdots,J$ can be obtained.

###  2.2 Dirichlet Process

Next, we introduce nonparametric Bayesian methods for clustering the
individuals and items. Dirichlet process (DP) is an universally used Bayesian
nonparametric method for capturing heterogeneity effects of the data. DP is
applied to Rasch model for capturing heterogeneity on both item difficulties
and person abilities. In this part for simplicity, we discuss the grouping
methods for abilities $\theta_{1},\cdots,\theta_{N}$. Assume we have Rasch
model as (1), the prior distribution for $\theta_{i},i=1,\cdots,N$ is the
Dirichlet process with concentration parameter $\alpha$ and base distribution
$G_{0}$ as below:

$\begin{split}&\theta_{i}\sim G\\\ &G\sim\text{DP}(\alpha,G_{0}).\end{split}$
(3)

Since the DP yields distributions $G$ that are almost surely discrete
(Ferguson,, 1973), the distinguished values among
$\theta_{1},\cdots,\theta_{N}$ induce a grouping among the abilities. The
observations $i$ and $i^{\prime}$ are in the same group when
$\theta_{i}=\theta_{i^{\prime}}$. For observation $i$, we have its grouping
allocation label $z_{i}\in\\{1,\cdots,K\\}$, and the ability of observation
$i$ is denoted by $\theta_{z_{i}}$, where $K$ is the total number of groups.
In this way, we have an equivalent model as

$\displaystyle\begin{split}z_{i}|\bm{\pi}&\sim\text{Discrete}(\pi_{1},\cdots,\pi_{K}),\\\
\theta_{k}&\sim G_{0},\quad k=1,\cdots,K,\\\
\bm{\pi}&\sim\text{Dirichlet}(\alpha/K,\cdots,\alpha/K),\end{split}$ (4)

where $\bm{\pi}=(\pi_{1},\cdots,\pi_{K})$, $K\rightarrow\infty$ (Ishwaran and
James,, 2001).

The distribution of $z_{i}$ is marginally given by a stick-breaking
construction (Sethuraman,, 1991) of DP as

$\begin{split}z_{i}&\sim\sum_{h=1}^{\infty}\pi_{h}\delta_{h},\\\
\pi_{h}&=\nu_{h}\prod_{\ell\leq h}(1-\nu_{\ell}),\\\
\nu_{h}&\sim\text{Beta}(1,\alpha),\end{split}$ (5)

where $\delta_{h}$ is the Dirac function with mass at $h$. Based on the stick-
breaking construction of DP, a MCMC algorithm (Ishwaran and James,, 2001) can
be used for parameter estimation. However, DP results in producing extraneous
clusters in the posterior, leading to inconsistent estimation on $K$, the
number of clusters, even when the sample size goes to infinity.

###  2.3 Mixture of Finite Mixtures Model

Miller and Harrison, (2013) showed that the posterior distribution on the
number of clusters does not converge to the true number of components under
DP. Miller and Harrison, (2018) proposed a modification of the Chinese
restaurant process (CRP) (Pitman,, 1995), called a mixture of finite mixtures
(MFM) model, to circumvent this issue:

$\displaystyle\begin{split}K&\sim p(\cdot),\\\ (\pi_{1},\ldots,\pi_{K})\mid
K&\sim\text{Dirichlet}(\gamma,\ldots,\gamma),\\\ z_{i}\mid
K,\bm{\pi}&\sim\sum_{k=1}^{K}\pi_{k}\delta_{k},\quad i=1,\ldots,N,\end{split}$
(6)

where $p(\cdot)$ is a proper probability mass function (p.m.f) on
$\\{1,2,\ldots,+\infty\\}$ and $\delta_{k}$ is a point-mass at $k$. A default
choice of $p(\cdot)$ is a $\mbox{Poisson}(1)$ distribution truncated to be
positive (Miller and Harrison,, 2018), which is assumed through the rest of
the paper. Like the stick-breaking representation in (5) of Dirichlet process,
the MFM also has a similar construction. If we choose
$(K-1)\sim\mbox{Poisson}(\lambda)$ and $\gamma=1$ in (6), the mixture weights
$\pi_{1},\cdots,\pi_{K}$ can be constructed as follows:

  1. 1.

Generate
$\eta_{1},\eta_{2},\cdots\overset{\text{iid}}{\sim}\text{Exp}(\lambda)$,

  2. 2.

$K=\min\\{j:\sum_{k=1}^{j}\eta_{k}\geq 1\\}$,

  3. 3.

$\pi_{k}=\eta_{k}$, for $k=1,\cdots,K-1$,

  4. 4.

$\pi_{K}=1-\sum_{k=1}^{K-1}\pi_{k}$.

The Gibbs samplers are easily constructed in stick-breaking framework
(Ishwaran and James,, 2001) for MFM. For ease of exposition, we refer MFM as
$\text{MFM}(\gamma,\lambda)$ with default choice of $p(\cdot)$ being
$\mbox{Poisson}(\lambda)$.

###  2.4 Hierarchical Rasch Model with MFM

In order to allow for simultaneously heterogeneity detection of person and
item parameters, the MFM process in Section 2.3 is introduced for the
parameters $\bm{b}$ and $\bm{\theta}$ in the item response model. The
hierarchical Rasch model with MFM is given as follows, for $i=1,\cdots,N$,
$j=1,\cdots,J$:

$\begin{split}y_{ij}&\sim\text{Bernoulli}(p_{ij}),\\\
\text{log}(p_{ij}/(1-p_{ij}))&=\text{logit}(p_{ij})=\theta_{z_{i}}-b_{g_{j}},\\\
b_{g_{j}}&\sim N(0,\psi_{b}^{-1}),\\\ \theta_{z_{i}}&\sim
N(0,\psi_{\theta}^{-1}),\\\ \psi_{\theta}&\sim\text{Gamma}(100,1),\\\
g_{j}|K_{b},\pi_{b},\lambda_{b}&\sim\text{MFM}(\gamma,\lambda_{b}),\\\
z_{i}|K_{\theta},\pi_{\theta},\lambda_{\theta}&\sim\text{MFM}(\gamma,\lambda_{\theta}),\end{split}$
(7)

where $\psi_{b}^{-1}=0.01$, $\bm{g}=(g_{1},\cdots,g_{J})$ and
$\bm{z}=(z_{1},\cdots,z_{N})$ denote the vectors of cluster configurations of
parameters $\bm{b}$ and $\bm{\theta}$, respectively. The value of $\psi_{b}$
and the prior for $\psi_{\theta}$ are set for mean controlling purpose. MFM
denotes the clustering method introduced in Section 2.3.

##  3 Bayesian Inference

###  3.1 Prior Specification and Posterior Distribution

For the hierarchical Rasch model with MFM introduced in Section 2.4, the set
of parameters is denoted as
$\Theta=\\{(\theta_{z_{i}},b_{g_{j}},K_{b},K_{\theta},\pi_{b},\pi_{\theta},\lambda_{b},\lambda_{\theta}):i=1,\cdots,N,j=1,\cdots,J\\}$.
Priors for the hyperparameters are $\lambda_{b}\sim\text{Log-Normal}(0,1)$ and
$\lambda_{\theta}\sim\text{Log-Normal}(0,1)$. With the prior distributions
specified above, the posterior distribution of these parameters based on the
data $D=\\{(y_{ij}):i=1,\cdots,N,j=1,\cdots,J\\}$ is given by

$\displaystyle\pi(\Theta|D)$ $\displaystyle\propto L(\Theta|D)\pi(\Theta)$
$\displaystyle\prod_{i=1}^{n}\prod_{j=1}^{J}f(y_{ij}|b_{g_{j}},\theta_{z_{i}})f(b_{g_{j}})f(\theta_{z_{i}})f(g_{j})f(z_{i})\pi(\Theta).$

The analytical form of the posterior distribution of $\Theta$ is unavailable.
Therefore, we carry out the posterior inference using the MCMC sampling
algorithm to sample from the posterior distribution and then obtain the
posterior estimates of the unknown parameters. The estimated clustering
allocations $g_{j}$’s and $z_{i}$’s are obtained from the posterior modes of
the MCMC samples, while the other parameter estimations are obtained through
posterior means of the MCMC samples. Computation is facilitated by the
nimble(de Valpine et al.,, 2017) package in R (R Core Team,, 2013), which uses
syntax similar to WinBUGS (Lunn et al.,, 2000) and JAGS(Plummer et al.,,
2003), but generates C++ code for faster computation.

###  3.2 Bayesian Model Comparison

In the hierarchical Rasch model with MFM, the prior distributions of
$\lambda_{b}$ and $\lambda_{\theta}$ have several choices, including the Gamma
distribution, Uniform distribution and Log-normal distribution. In order to
choose the most suitable priors for these parameters, model comparison
criteria under the Bayesian framework are used. One of the most commonly used
criterion is the Deviance Information Criteria (DIC; Spiegelhalter et al.,,
2002), which measures the fitness as well as the complexity of the model. The
DIC is defined as

$\displaystyle\text{DIC}=\text{Dev}(\bar{\Theta})+2p_{D},$ (8)

where $\text{Dev}(\bar{\Theta})$ is the deviance function,
$p_{D}=\overline{\text{Dev}}({\Theta})-\text{Dev}(\bar{\Theta})$ is the
effective number of model parameters, $\bar{\Theta}$ is the posterior mean of
$\Theta$, and $\overline{\text{Dev}}({\Theta})$ is the posterior mean of
$\text{Dev}(\Theta)$. A model with a smaller DIC indicates a preferred model.

Another model comparison criterion is the logarithm of the pseudo-Marginal
likelihood (LPML; Ibrahim et al.,, 2013), which is obtained through the
conditional predictive ordinate (CPO) values. With
$\bm{y}^{*}_{(-i)}=(\bm{y}_{1},\ldots,\bm{y}_{i-1},\bm{y}_{i+1},\ldots,\bm{y}_{N})$
denoting the binary response vector with the $i$th subject response deleted,
CPO is regarded as leave-one-out-cross-validation under Bayesian framework,
and it estimates the probability of observing $\bm{y}_{i}$ in the future after
having already observed $\bm{y}^{*}_{(-i)}$. The CPO for the $i$th subject is
calculated as:

$\text{CPO}_{i}=\int
f(\bm{y}_{i}|\bm{b},\bm{\theta})\pi(\bm{b},\bm{\theta}|\bm{y}^{*}_{(-i)})d(\bm{b},\bm{\theta}),$
(9)

where

$\pi(\bm{b},\bm{\theta}|\bm{y}^{*}_{(-i)})=\frac{\prod_{q\neq
i}f(\bm{y}_{q}|\bm{b},\bm{\theta})\pi(\bm{b},\bm{\theta}|\bm{y}^{*}_{(-i)})}{c\big{(}\bm{y}^{*}_{(-i)}\big{)}},$

and $c\big{(}\bm{y}^{*}_{(-i)}\big{)}$ is the normalizing constant. Within the
Bayesian framework, a Monte Carlo estimate of the CPO is obtained as:

$\widehat{\text{CPO}}_{i}^{-1}=\frac{1}{M}\sum_{t=1}^{M}\frac{1}{f(\bm{y}_{i}|\bm{b}_{t},\bm{\theta}_{t})},$
(10)

An estimate of the LPML subsequently is calculated as:

$\widehat{\text{LPML}}=\sum_{i=1}^{n}\log(\widehat{\text{CPO}}_{i}).$ (11)

A model with a larger LPML value will be selected.

In addition, the effective number of parameters $p_{D}$ is applied to compare
the Rasch model and the proposed one from a model complexity point of view. A
larger number of $p_{D}$ indicates higher complexity of the model. Besides,
the area under the curve (AUC) value is calculated to measure the prediction
precision of the model, which is computed via the R package pROC (Robin et
al.,, 2011). A higher AUC value indicates the better performance under the
prediction.

##  4 Simulation

In order to evaluate the performance of the proposed method for clustering and
parameter estimation, simulation studies are conducted in this section. Two
situations with different number of items and subjects were considered, and
under each situation, two designs were set up. For both designs, 100 datasets
were generated with $N$ subjects and $J$ items in each dataset. The binary
response $y_{ij}$ with the probability $p_{ij}$ was generated according to
(1). Three clusters were assumed for both the subject parameters
$\theta_{i}(i=1,\cdots,N)$ and the item parameters $b_{j}(j=1,\cdots,J)$. In
the first simulation design, the values of $b_{j}(j=1,\cdots,J)$ and
$\theta_{i}(i=1,\cdots,N)$ were both randomly chosen from $\\{-2,0,2\\}$. Take
$b_{j}$’s for an example. The $j$th difficulty parameter was randomly assigned
one of the values of $\\{-2,0,2\\}$ as its true value, through which, three
clusters for the difficulty parameters were formed. In the second simulation
design, errors following a normal distribution $N(0,0.5^{2})$ were added to
the values of $b_{j}(j=1,\cdots,J)$ and $\theta_{i}(i=1,\cdots,N)$. Note that
the second condition is a more challenging situation because some disturbance
of the clustering of items and persons are added in the data by the error
$N(0,0.5^{2})$.

The performance of the posterior estimates were evaluated by the mean absolute
bias (MAB), the mean standard deviation (MSD) and the mean of mean squared
error (MMSE) in the following ways, take $\theta_{i}$ as an example:

$\displaystyle\text{MAB}=\frac{1}{100}\sum_{r=1}^{100}\frac{1}{N}\sum_{i=1}^{N}\left|\hat{\theta}_{ir}-\theta_{i}\right|,$
$\displaystyle\text{MSD}=\sqrt{\frac{1}{100}\sum_{r=1}^{100}\frac{1}{N}\sum_{i=1}^{N}\left(\hat{\theta}_{ir}-\bar{\hat{\theta}}_{i}\right)^{2}},$
$\displaystyle\text{MMSE}=\frac{1}{100}\sum_{r=1}^{100}\frac{1}{N}\sum_{i=1}^{N}\left(\hat{\theta}_{ir}-\theta_{i}\right)^{2}.$

To evaluate the accuracy of clustering, apart from showing the frequency of
cluster numbers in the simulations, the Precision, Recall and Rand Index (RI;
Rand,, 1971) were also applied. Let $TP$, $FP$, $TN$ and $FN$ denote the true
positive, false positive, true negative and false negative scores,
respectively. The Precision was calculated via $TP/(TP+FP)$, the Recall was
calculated via $TP/(TP+FN)$, and the RI was obtained by
$(TP+TN)/(TP+FP+TN+FN)$ (Ma et al.,, 2019). The average RI (MRI) was
calculated as the mean of RIs of the 100 simulations, and similarly the
average Precision (MP) and average Recall (MR). For both these three measures,
a higher value represents higher accuracy of clustering. In addition, the
average effective number (MpD), the average AUC (MAUC), and the average
computation time in seconds (MTime) were also calculated to compare the
complexity of the Bayesian Rasch model and the proposed model. Under different
settings, with a thinning interval of 2, 10000 samples were kept for
calculating the posterior estimates after a burn-in of 20000 samples.
Convergence was checked via the traceplots and autocorrelation plots of the
chains.

We firstly consider a situation with $N=200$ subjects and $J=60$ items. The
MpD, MAUC and MTime values of the two models under two designs are shown in
Table 1.

Table 1: Model performance results of the proposed model and the Bayesian Rasch model under two designs with $N=200$ and $J=60$ Design | The proposed model |  | Bayesian Rasch model
---|---|---|---
MpD | MAUC | MTime |  | MpD | MAUC | MTime
The first design | 27.698 | 0.831 | 4126.742 |  | 243.029 | 0.839 | 507.783
The second design | 85.518 | 0.834 | 4225.178 |  | 242.823 | 0.852 | 515.304

From Table 1, we see that under two different designs with $N=200$ and $J=60$,
the MpD values of the Bayesian Rasch model are larger than those of the
proposed model, meaning that the Bayesian Rasch model has greater complexity
and has more parameters to be estimated than the proposed model. On model
prediction performance, the MAUC values of the two models are merely the same,
which indicates that the prediction precision of the two models under
different simulation designs are similar. On the average computation time, the
proposed model requires much more time than the Bayesian Rasch model.

The frequencies of the number of clusters for $\bm{\theta}$ and $\bm{b}$ under
these two designs with $N=200$ and $J=60$ are shown in Figure 1. Under design
1, the percentages of clustering the $\theta_{i}$’s and the $b_{j}$’s into 3
clusters are nearly 60% and 100%, respectively. These results mean that for
$\theta_{i}$’s, about 60 out of 100 replicates, the number of clusters is
correctly inferenced, and for $b_{j}$’s, the proportion is 100%. In addition,
the MRI values in Table 2 show the good performance of clustering allocation,
which reach to 0.961 and 1.000 for these parameters. The MRI values indicate
that about 96% and 100% of the time, any two persons or items that are in the
same cluster are correctly grouped together. The MP and MR values are
consistent with the MRI values, which verify our conclusions. Under design 2,
the percentages of correct clustering and both the MRI, MP and MR values of
the two parameters are lower compared to design 1, but the most frequently
chosen cluster numbers are also 3 clusters and the MRI values are both larger
than 0.85, meaning that when the true values of $\theta_{i}$’s and $b_{j}$’s
are generated with errors, the precision of clustering are still high.

(a) Number of clusters of $\theta$ under design 1 (b) Number of clusters of
$b$ under design 1 (c) Number of clusters of $\theta$ under design 2 (d)
Number of clusters of $b$ under design 2

Figure 1: Frequency of number of clusters of $\theta$ and $b$ under different
designs with $N=200$ and $J=60$

The results of posterior estimates under two designs with $N=200$ and $J=60$
are also shown in Table 2. The MABs and MMSEs of the parameters are both
within a reasonable range. In general, performance of posterior estimates
under design 1 is better than those under design 2.

Table 2: Simulation results of the proposed model in two simulation designs with $N=200$ and $J=60$ Design | Parameter | MAB | MSD | MMSE | MRI | MP | MR
---|---|---|---|---|---|---|---
The first design | $\theta$ | 0.582 | 0.046 | 0.490 | 0.961 | 0.960 | 0.922
$b$ | 0.590 | 0.041 | 0.509 | 1.000 | 1.000 | 1.000
The second design | $\theta$ | 0.729 | 0.047 | 0.767 | 0.862 | 0.787 | 0.803
$b$ | 0.723 | 0.060 | 0.778 | 0.899 | 0.862 | 0.833

Another situation with $N=200$ subjects and $J=100$ items are considered.
Similarly, the MpD, MAUC and MTime values of these two models under two
designs are shown in Table 3. It is shown that the MpD values of the proposed
model under two different settings are much lower than those of the Bayesian
Rasch model, especially under the first design. The similar MAUC values show
that these two models have merely the same prediction precision, and the
proposed model still has high precision in addition to clustering. However,
the proposed model spends longer time for computation than the Bayesian Rasch
model. These conclusions are similar to those when $N=200$ and $J=60$.

Table 3: Model performance results of the proposed model and the Bayesian Rasch model under two designs with $N=200$ and $J=100$ Design | The proposed model |  | Bayesian Rasch model
---|---|---|---
MpD | MAUC | MTime |  | MpD | MAUC | MTime
The first design | 30.819 | 0.832 | 7158.913 |  | 289.058 | 0.840 | 999.524
The second design | 90.807 | 0.834 | 7487.103 |  | 288.721 | 0.853 | 971.540

Figure 2 shows that frequency of number of clusters under two designs with
$N=200$ and $J=100$. Under the first design, nearly 75% and $100\%$ of the
parameters $\theta_{i}$’s and $b_{j}$’s are grouped into three clusters.
Comparing these results with the results under design 1 with $J=60$, we see
that the clustering performance of parameters $\theta_{i}$’s improves, and the
percentage of clustering into 3 clusters increases from nearly 60% to 75%. The
higher MP, MR and MRI values shown in Table 4 also verifies this conclusion.
Under the second design, the percentages of clustering the parameters
$\theta_{i}$’s and $b_{j}$’s into 3 clusters are 66% and 63%, respectively.
Comparing with the results in Figure 1, the clustering performance of
$\theta_{i}$’s improves while the clustering performance of $b_{j}$’s
degenerates. This may due to the fact that a greater number of items increase
the difficulty of clustering $b_{j}$’s. The MRI values in Table 4 are both
larger than 0.88, meaning that the precision of clustering allocation under
this design is high.

(a) Number of clusters of $\theta$ under design 1 (b) Number of clusters of
$b$ under design 1 (c) Number of clusters of $\theta$ under design 2 (d)
Number of clusters of $b$ under design 2

Figure 2: Frequency of number of clusters of $\theta$ and $b$ under different
designs with $N=200$ and $J=100$

In Table 4, the results of posterior estimates under two designs with $N=200$
and $J=100$ are also shown. Comparing with the results in Table 2, posterior
estimation performance under the first setting with a smaller number of items
are relatively better than those under the one with $J=100$. However, the MABs
and MMSEs of the parameters are both within a reasonable range and in general,
performance of posterior estimates under design 1 is better than that under
design 2.

Table 4: Simulation results of the proposed model in two simulation designs with $N=200$ and $J=100$ Design | Parameter | MAB | MSD | MMSE | MRI | MP | MR
---|---|---|---|---|---|---|---
The first design | $\theta$ | 0.591 | 0.033 | 0.522 | 0.989 | 0.995 | 0.974
$b$ | 0.613 | 0.030 | 0.548 | 0.999 | 1.000 | 1.000
The second design | $\theta$ | 0.739 | 0.036 | 0.804 | 0.882 | 0.822 | 0.824
$b$ | 0.741 | 0.046 | 0.817 | 0.897 | 0.857 | 0.832

##  5 Application to English Test Data

###  5.1 Data Description and Preliminary Analysis

The dataset we analyzed is collected from the English exam results of the
2017-2018 academic year from a public middle school in China. This dataset was
previously analyzed in Liu et al., (2019). This exam is a midterm English exam
for grade 8 students in fall 2017. In our analysis, the total number of
students is 78 and the total number of items is 50. We used this dataset to
illustrate the MFM Rasch model estimation. For each student, the average
number of items of correct response is 23.09, while for each item, the average
number of students that answer correctly is 36.02, so the average proportion
of correct answer is 0.46. The boxplot of the proportions of correct answer
for the items and the students are shown in Figure 3.

Figure 3: Boxplot of the proportions of correct answer for the items and the
students in the English test data.

From Figure 3, we observe that the correct response proportions vary from
different items as well as different students, indicating that the items are
of different difficulty levels and the students are of different ability
levels. The Bayesian Rasch model in (1) was applied to analyze the data. For
the unknown parameters $b_{j}$ and $\theta_{i}$, the prior distributions were
normal distributions, i.e., $b_{j}\sim N(0,1)$ and $\theta_{i}\sim
N(0,\psi^{-1})$ with hyperprior $\psi\sim\text{Gamma}(0.001,0.001)$, for
$i=1,\cdots,78$, $j=1,\cdots,50$. With a thinning interval of 2, 10000 samples
were kept for calculating the posterior estimates after a burn-in of 20000
samples. The effective number of parameters $p_{D}$ of this model was also
calculated and turned out to be 107.676. The AUC value was 0.739. Posterior
estimates as well as the 68% highest posterior density (HPD) intervals of the
parameters $b$ and $\theta$ are shown in Figure 4.

(a) Posterior estimates of $\theta$ (b) Posterior estimates of $b$

Figure 4: Posterior estimates and 68% HPD intervals of parameters in the Rasch
model.

From Figure 4, we see that for the items’ difficulties and individuals’
abilities, the posterior means may be roughly divided into several groups.
Take parameter $\theta$’s as an example. Students with the posterior means
which are above 0.5 can be grouped together, indicating that the latent
abilities of students belonging to this group are close to each other and may
be nearly the same. For the ones whose posterior means are smaller than 0 may
also be grouped together. Due to the patterns we observed, the parameters of
different items/individuals may be clustered into different groups. A Rasch
model with clustering items and individuals simultaneously would capture the
heterogeneity patterns shown in the data.

###  5.2 MFM Rasch Clustering Analysis

In order to detect the clusters within the parameters, the proposed
methodology was applied to the English test data. We firstly applied the model
comparison criteria to select the most suitable priors for the parameters
$\lambda_{b}$ and $\lambda_{\theta}$. Three priors were compared, including
$\text{Gamma}(1,1)$, $\text{Unif}(0,1)$ and $\text{Log-Normal}(0,1)$. The
results of the model selection criteria are shown in Table 5.

Table 5: DIC and LPML values for different priors for $\lambda_{b}$ and $\lambda_{\theta}$ in the proposed model | $\text{Gamma}(1,1)$ | $\text{Unif}(0,1)$ | $\text{Log-Normal}(0,1)$
---|---|---|---
DIC | 4673.480 | 4681.919 | 4675.817
LPML | -2279.155 | -2282.574 | -2280.213

From Table 5, we find that the DIC and LPML criteria both prefer the prior
$\text{Gamma}(1,1)$. Thus, we used this prior for $\lambda_{b}$ and
$\lambda_{\theta}$ in our final analysis. With this prior, the effective
number of parameters was calculated and the value is 57.414. By comparing this
value with that of the Bayesian Rasch model (107.676), we see that the
complexity of our proposed model is lower and it fits less parameters compared
to the Rasch model. In addition, the AUC value of the proposed model is 0.738,
which is very close to that of the Bayesian Rasch model (0.739). This
indicates that the proposed model can not only cluster the parameters, but
also has a similar prediction performance with the Rasch model.

Finally, for the latent ability parameters $\theta_{i}$’s, two clusters are
identified, indicating that these 78 students are classified into 2 groups
according to different levels of latent ability. For the difficulty parameter
$b_{j}$’s, three clusters are identified, meaning that the 50 questions have 3
levels of difficulty. The clustering and estimation results are shown in Table
6. The clustering allocations are shown in Figure 5, where the x-axis
represents student index and item number, respectively, in (a) and (b).

Table 6: The clustering and estimation results for the real data analysis using the proposed model  Parameter | $\bm{\theta}$ | $\bm{b}$
---|---|---
Cluster 1 | Count | 22 | 18
Estimate | 0.612 (0.448, 0.788) | 0.643 (0.508, 0.775)
Cluster 2 | Count | 56 | 16
Estimate | -0.519 (-0.669, -0.390) | 0.005 (-0.144, 0.154)
Cluster 3 | Count | - | 16
Estimate | - | -0.728 (-0.857, -0.592)

(a) Posterior estimates of $\theta$ (b) Posterior estimates of $b$

Figure 5: Clustering allocations of individual abilities and item difficulties
in the proposed model

From Table 6 and Figure 5, we see that for latent ability, students are
clustered into 2 groups, and 56 among the 78 students (71.795% of the
students) are in Cluster 2, whose latent abilities are relatively lower. For
the students in Cluster 1, their latent abilities are higher. In order to
decide if these results are consistent to the original data, we calculated the
average proportion of correct answer of students within these two clusters
respectively. The results are shown in Table 7.

Table 7: Average proportion of correct answer in different clusters of latent ability Cluster | Cluster 1 | Cluster 2 | All
---|---|---|---
Proportion of correct answer | 68.45% | 37.43% | 46.18%

From Table 7, we notice that students in Cluster 1 have a proportion of
correct answer that is above the average level, while the ones in Cluster 2
have a lower proportion of correct answer, which is consistent with the
clustering results of latent ability.

For item difficulties, three clusters are identified. There are 18 items (36%
of the test length) are grouped into Cluster 1, while for Cluster 2 and
Cluster 3, there are 16 items (32% of the test length) in each cluster. It is
shown that questions belonging to Cluster 1 are the most difficult, while the
ones in Cluster 3 are the easiest. Similarly, we calculated the proportion of
correctly answer the questions in different clusters to see whether the
estimation results are consistent to the real data. The average proportions of
correct answer in three clusters of items are shown in Table 8.

Table 8: Average proportion of correct answer in different clusters of difficulty Cluster | Cluster 1 | Cluster 2 | Cluster 3 | All
---|---|---|---|---
Proportion of correct answer | 26.64% | 45.59% | 68.75% | 46.18%

From Table 8 we see that students have the highest proportion of correct
answer to items in Cluster 3, but the lowest proportion in Cluster 1,
indicating that items in Cluster 1 are the most difficult while items in
Cluster 3 are relatively easier. The difficulty level of items in Cluster 2
are at an average level since the proportion of correct answer in Cluster 2 is
similar to the average proportion of correct answer of all the items. These
results are also consistent with the conclusions we made according to the
clustering and estimation results.

The item characteristic curves using the item parameter estimates of the three
clusters from the proposed model are shown in Figure 6. Since the item
difficulty parameters are clustered into three groups, we have three item
characteristic curves for the proposed IRT model. For the three curves with
different item difficulty levels, the locations of the curves along the
ability scale are different, meaning that the difficulty of the items within
different groups vary. The curve of items within cluster 3 is at the left-hand
side, representing items in this cluster are the easiest since the probability
of correct response is high for low-ability students and approaches 1 for
high-ability students. Curve of cluster 2 is in the middle, meaning that items
within this cluster are of medium difficulty. The curve for items in cluster 1
is the right-hand curve, representing the items in this cluster are hard
items.

Figure 6: Item characteristic curves of the hierarchical Rasch model with MFM

##  6 Summary and Discussion

In the regular Rasch model, items are treated as invariant across test-takers,
thus an item’s difficulty is the same across all test-takers. Also, test-
takers are assumed to be from a homogeneous group. This homogeneity of item
and person parameters is not always feasible in reality. LCA or mixture IRT
model has been a major tool to investigate the heterogeneity or the violation
of the homogeneity assumption of item and person parameters. LCA could detect
heterogeneity in persons, permitting person clustering. Mixture IRT modeling
allows person and item heterogeneity, detecting different subgroups of test-
takers which have different constellations of item parameters. As pointed out,
these LCA and mixture IRT modeling requires a pre-specification of the number
of clusters by users before estimation. Also, a two-stage procedure, where
several models are estimated and compared to pick the most appropriate model,
is not efficient and disregards the uncertainties which take place in the
selection procedure.

In this paper, we proposed a Bayesian nonparametric model for capturing the
heterogeneity of items and individuals simultaneously under the Rasch model.
Our proposed model, MFM Rasch model, simultaneously estimate the number of
clusters and clustering configurations for both items and individuals. It does
not require a priori specification of the number of clusters by users. Because
of the simultaneous clustering of both the item and person parameters, it
avoids the problem of ignoring the uncertainties in the two-stage procedure
used in LA and mixture IRT modeling. A practical side of the MFM Rasch
modeling could include situations where users of the model are interested in
identifying clusters of items and persons, so that the information of clusters
may be used for their assessment purposes. For example, information on item
clusters may be useful in the test construction stage and information on
person clusters may also be utilized for student-placement purposes or
diagnostic purposes.

Our simulation results and real data analysis indicated that the MFM Rasch
estimation could achieve both accurate estimation and clustering performance
although the scope of the simulation is limited and the real data set may not
be the most typical application scenario regarding the sample size and test
length. Future endeavor should include more systematic studies through wider
scope of simulations (e.g., varying the number of clusters, sample sizes, and
test lengths) to provide a complete picture of the performance of the MFM
Rasch model. In addition, at least three topics beyond the scope of this paper
could be worth further investigation. First, extending the MFM Rasch model
into the 2-parameter or 3-parameter item response theory model seems to be an
interesting future work. Second, our proposed computing algorithm is based on
stick-breaking construction. Incorporating Pólya urn scheme in posterior
sample (Neal,, 2000) will speed up the MCMC algorithm. Third, extending the
current model to consider auxiliary information of individuals such as
previous exam performance for clustering detection may be pursued in future
research.

## References

  * Akaike, (1998) Akaike, H. (1998).  Information theory and an extension of the maximum likelihood principle.  In Selected papers of hirotugu akaike, pages 199–213. Springer.
  * Bartolucci, (2007) Bartolucci, F. (2007).  A class of multidimensional irt models for testing unidimensionality and clustering items.  Psychometrika, 72(2):141.
  * Bartolucci et al., (2017) Bartolucci, F., Farcomeni, A., and Scaccia, L. (2017).  A nonparametric multidimensional latent class irt model in a bayesian framework.  Psychometrika, 82(4):952–978.
  * de Valpine et al., (2017) de Valpine, P., Turek, D., Paciorek, C. J., Anderson-Bergman, C., Lang, D. T., and Bodik, R. (2017).  Programming with models: writing statistical algorithms for general model structures with NIMBLE.  Journal of Computational and Graphical Statistics, 26(2):403–413.
  * Ferguson, (1973) Ferguson, T. S. (1973).  A Bayesian analysis of some nonparametric problems.  The Annals of Statistics, 1(2):209–230.
  * (6) Geng, J., Bhattacharya, A., and Pati, D. (2019a).  Probabilistic community detection with unknown number of communities.  Journal of the American Statistical Association, 114(526):893–905.
  * (7) Geng, J., Shi, W., and Hu, G. (2019b).  Bayesian nonparametric nonhomogeneous poisson process with applications to usgs earthquake data.  arXiv preprint arXiv:1907.03186.
  * Gnaldi et al., (2016) Gnaldi, M., Bacci, S., and Bartolucci, F. (2016).  A multilevel finite mixture item response model to cluster examinees and schools.  Advances in Data Analysis and Classification, 10(1):53–70.
  * Hagenaars and McCutcheon, (2002) Hagenaars, J. A. and McCutcheon, A. L. (2002).  Applied latent class analysis.  Cambridge University Press.
  * Hu et al., (2020) Hu, G., Geng, J., Xue, Y., and Sang, H. (2020).  Bayesian spatial homogeneity pursuit of functional data: an application to the us income distribution.  arXiv preprint arXiv:2002.06663.
  * Ibrahim et al., (2013) Ibrahim, J. G., Chen, M.-H., and Sinha, D. (2013).  Bayesian survival analysis.  Springer Science & Business Media.
  * Ishwaran and James, (2001) Ishwaran, H. and James, L. F. (2001).  Gibbs sampling methods for stick-breaking priors.  Journal of the American Statistical Association, 96(453):161–173.
  * Lazarsfeld and Henry, (1968) Lazarsfeld, P. F. and Henry, N. W. (1968).  Latent structure analysis.  Houghton Mifflin Co.
  * Liu, (2019) Liu, X. (2019).  Three Contributions to Latent Variable Modeling.  PhD thesis, Columbia University.
  * Liu et al., (2019) Liu, Y., Hu, G., Cao, L., Wang, X., and Chen, M.-H. (2019).  A comparison of monte carlo methods for computing marginal likelihoods of item response theory models.  Journal of the Korean Statistical Society.
  * Lunn et al., (2000) Lunn, D. J., Thomas, A., Best, N., and Spiegelhalter, D. (2000).  Winbugs-a bayesian modelling framework: concepts, structure, and extensibility.  Statistics and computing, 10(4):325–337.
  * Ma et al., (2019) Ma, Z., Xue, Y., and Hu, G. (2019).  Heterogeneous regression models for clusters of spatial dependent data.  page arXiv:1907.02212.
  * Miller and Harrison, (2013) Miller, J. W. and Harrison, M. T. (2013).  A simple example of Dirichlet process mixture inconsistency for the number of components.  In Advances in Neural Information Processing Systems, pages 199–206.
  * Miller and Harrison, (2018) Miller, J. W. and Harrison, M. T. (2018).  Mixture models with a prior on the number of components.  Journal of the American Statistical Association, 113(521):340–356.
  * Miyazaki and Hoshino, (2009) Miyazaki, K. and Hoshino, T. (2009).  A Bayesian semiparametric item response model with dirichlet process priors.  Psychometrika, 74(3):375–393.
  * Neal, (2000) Neal, R. M. (2000).  Markov chain sampling methods for Dirichlet process mixture models.  Journal of Computational and Graphical Statistics, 9(2):249–265.
  * Pan and Huang, (2014) Pan, J.-C. and Huang, G.-H. (2014).  Bayesian inferences of latent class models with an unknown number of classes.  Psychometrika, 79(4):621–646.
  * Pitman, (1995) Pitman, J. (1995).  Exchangeable and partially exchangeable random partitions.  Probability Theory and Related Fields, 102(2):145–158.
  * Plummer et al., (2003) Plummer, M. et al. (2003).  Jags: A program for analysis of bayesian graphical models using gibbs sampling.  In Proceedings of the 3rd international workshop on distributed statistical computing, volume 124, page 10. Vienna, Austria.
  * R Core Team, (2013) R Core Team (2013).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Rand, (1971) Rand, W. M. (1971).  Objective criteria for the evaluation of clustering methods.  Journal of the American Statistical Association, 66(336):846–850.
  * Rasch, (1980) Rasch, G. (1960/1980).  Probabilistic models for some intelligence and attainment tests (ExpandEd Edition 1980).  Danish Institute for Education Research, The University of Chicago Press, Chicago: Copenhagen.
  * Robin et al., (2011) Robin, X., Turck, N., Hainard, A., Tiberti, N., Lisacek, F., Sanchez, J.-C., and Müller, M. (2011).  pROC: an open-source package for r and s+ to analyze and compare roc curves.  BMC Bioinformatics, 12(1):77.  R package version 0.3.0.
  * Rost, (1990) Rost, J. (1990).  Rasch models in latent classes: An integration of two approaches to item analysis.  Applied Psychological Measurement, 14(3):271–282.
  * Rost and Langeheine, (1997) Rost, J. E. and Langeheine, R. E. (1997).  Applications of latent trait and latent class models in the social sciences. Waxmann Publishing Co.
  * Schwarz, (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  The Annals of Statistics, 6(2):461–464.
  * Sethuraman, (1991) Sethuraman, J. (1991).  A constructive definition of Dirichlet priors.  Statistics Sinica, 4(2):639–650.
  * Spiegelhalter et al., (2002) Spiegelhalter, D. J., Best, N. G., Carlin, B. P., and Van Der Linde, A. (2002).  Bayesian measures of model complexity and fit.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64(4):583–639.
  * Yang et al., (2011) Yang, H., O’Brien, S., and Dunson, D. B. (2011).  Nonparametric bayes stochastically ordered latent class models.  Journal of the American Statistical Association, 106(495):807–817.


# Studying Accuracy of Machine Learning Models Trained on Lab Lifting Data in
Solving Real-World Problems Using Wearable Sensors for Workplace Safety

Joseph Bertrand, Nick Griffey, Ming-Lun Lu, Rashmi Jha

###### Abstract

Porting ML models trained on lab data to real-world situations has long been a
challenge. This paper discusses porting a lab-trained lifting identification
model to the real-world. With performance much lower than on training data, we
explored causes of the failure and proposed four potential solutions to
increase model performance.

##  I Introduction

Back pain is a leading cause of occupation-induced disability and results in
large amounts of lost productivity annually [1]. Most occupational back pain
is caused by roles requiring repetitive tasks such as heavy lifting of
objects. There are many ergonomic guidelines for lifting with the goal of
reducing the risk of back pain, but it is difficult for workers to
consistently follow these guidelines in all situations. Many may be unable to
lift in the required manner due to poor workstation designs or physical
limitations [2]. Monitoring a lifting workstation can ensure compliance with
safe lifting techniques as well as aid in determining if workers are
consistently lifting in an unsafe manner to perform the work the job requires.
However, this type of monitoring requires significant overhead and doesn’t
provide sufficient data to create an objective methodology to assess risk. An
automated system would allow users to classify risk with minimal overhead and
provide real-time feedback to workers in an attempt to reduce long-term risk.

Both lift classification and lift detection are known as human activity
recognition (HAR), which have been extensively studied in machine learning
under various sensor modalities [3]. HAR systems have been successful in
video-based deployments [4]. However, automatic lift assessment would require
cameras to be placed in every lift location, which is prohibitively expensive
and raises concerns for the privacy of the workers. Furthermore, it is
impractical for temporary workstations such as construction zones and contract
work. In contrast, Inertial Measurement Unit (IMU) sensors are wearable, can
detect motion regardless of workplace location, and don’t have the privacy
concerns associated with cameras. IMU sensors are widespread in common
consumer devices, such as smart watches and phones and therefore provide
significant advantage for HAR problems.

Previous work was done to develop a model that could identify a lifting event
from a subset of laboratory-gathered data with an F1 score of 97% [5]. F1 is
the harmonic mean of precision and recall and is commonly used to evaluate
binary classifiers. The evaluation data was randomly sampled using K-fold
Cross Validation. The model trained was able to reliably identify lifting and
non-lifting events from data gathered under the same conditions as its
training data. However, when this model was applied to a real-world
environment (data collection in lab and real world is discussed in section
II), the F1 score dropped to 32.8%, showing that while the model could
reliably identify lifting events within a dataset environment, it failed to
identify more general lifting events.

Poor general performance on a model is not uncommon, and there is a
significant amount of general guidance published[6]. However, most of this
work pertains to datasets where the capture device use does not significantly
affect the data. For example, creating an object recognition model from
images, the camera settings do not majorly change the underlying dataset,
since all reasonable cameras produce clear images, either with or without the
relevant recognition object. However, in HAR models, the location of the IMU
sensor majorly changes the data, and it is not currently possible to
standardize IMU data from sensors placed in different locations on the body.
This makes it unfeasible to use some common techniques for improving model
performance, such as transfer learning and generalizing training data from
public datasets. Therefore, adjustments to previous techniques were developed
to make them suitable for HAR model applications.

##  II Background

###  II-A Data Collection

To gather data, each subject had six IMU sensors attached to various locations
on their body. One sensor on each wrist, one on the upper-right thigh, one on
the upper back, one on the upper-right arm, and one on the waist (figure 1).
Each IMU sensor gathered accelerometer and gyroscopic data at a sampling rate
of 25hz.

Figure 1: Placement of IMU sensors on subjects

####  II-A1 Lab Data

Lab data was gathered in Phase 1 & 2 of the NIOSH research (IRB approved study
protocol: 16-DART-05XP). Therefore, it may be referred to as phase 1&2 data
throughout this paper. To gather the data, each subject started in the center
of a room. They would then walk forward and pick up a small, wired box (this
denoted the start of the lift). They would then turn and walk with the box
across the room and set the box down on a table. Setting the box down denoted
the end of the lifting event. After they set the box down, they would turn and
walk back to the center of the room, where the trial would stop.

NIOSH recorded videos of each trial and used a motion capture system (MoCap)
to accurately track the location of the box. The start and end lift times were
determined from the motion capture system. To identify the start and stop lift
time for training and testing the IMU sensor model, code was written to
synchronize the IMU data with the MoCap data. The time of the lift, in
hh:MM:ss:ms format, was stored, from MoCap data, as metadata for each lift
trial and each IMU sensor data point had a Unix Epoch timestamp associated
with it. To synchronize the two, code was written to extract the time of day
from the Unix timestamp, then using the time of day from the metadata, and the
knowledge that the sampling rate of the IMUs was 25HZ, we could calculate
which IMU frame the lift started and ended. The calculations and code were
manually verified.

####  II-A2 Real-World Data

Real-world data was gathered in Phase 3 of the NIOSH research (IRB approved
study protocol: 16-DART-05XP). Therefore, it may be referred to as phase 3
data throughout this paper. Phase 3 data was gathered in the field at a light-
industrial plant. Subjects wore the sensors in the same arrangements as they
did in the phase 1&2 data. Two Sony Handycam digital cameras were placed at
two angles of view (side and front view) of the subject to allow the NIOSH
team to review and manually label the beginning and ending of each lift. The
labeled lifts were used to assess the model in the evaluation stage (See
section V). There were some points where certain sensors were non-functional.
That data was removed from the dataset to maintain consistency with the
expected six-sensor data format.

To perform model evaluation, the team sampled a small segment of phase 3 data
and used subjective analysis to determine if the model was performing well. As
model development continued, the team decided to go back and label the phase 3
dataset so models could be compared objectively. After balancing the dataset
there were approximately 950 seconds of data used for evaluation.

###  II-B Identification Model Training

For the purposes of this paper the beginning of the lifting motion (BOL) was
defined by when the object that the subject was lifting began to move, and the
end of the lifting motion (EOL) was defined by when the subject had fully
lifted the object and had either placed it down, or brought it close to their
body. We defined lifting with these parameters because the majority of the
risk associated with the lifting motion is within these two time intervals [7]
The lifting labeling in Phase 1&2 training data differed slightly from this
definition. For the NIOSH labeling, two researchers reviewed the trials and
identified the frame number where the grid started to move (Beginning of Lift)
and where the grid was set down completely by two hands (End of Lift). If
there was any doubt, the researchers would discuss and reach an agreement on
the final frame [8]. In order to fit the NIOSH labeling to the definition of
lifting used in this project, the EOL was adjusted to be 1.52 seconds after
the BOL. This time was chosen by reviewing the lifting videos and finding the
average time it took for a subject to complete the core lifting motion. For
the evaluation dataset, the beginning of lift and end of lift labels provided
by NIOSH matched our lifting definition very closely, so no adjustments were
necessary.

To train the lift identification model, the phase 1&2 dataset was divided
equally into two categories, lifting events and non-lifting events. Since the
raw dataset contained more non-lifting events than lifting events, the dataset
was balanced randomly to include roughly the same number of both classes. The
lifting event was a 1.2 second window, starting from the beginning of lift
index provided by NIOSH. 1.2 seconds was the approximate average lifting time
based on examining videos of the subject’s lifts. Everything outside that 1.2
second window was considered to be a non-lifting event. There were also
several files provided which contained no lifting events, and instead
contained events such as walking, sitting, or jumping. Those files were also
considered to be non-lifting events.

Once a balanced dataset was generated, the model was trained using keras’s
.fit method. The hyperparamaters batch size, feature length, epoch, and
validation split, were varied in model in an attempt to find the combination
that provided the best performance. A total of 3456 models were trained, using
a combination of various predefined hyperparamater combinations. All models
had the same neuron architecture, consisting of an input layer, followed by a
128-layer deep LSTM layer, followed by two Dense functions of width 5, then
with an dense output layer of size 1 to provide the final result. This model
architecture was evaluated on phase 1&2 data and had performance parity with
previous work [5].

Each model was then evaluated on phase 3 data. The evaluation consisted of
obtaining the accuracy, F1, and confusion matrix from each model. The results
were cataloged to create the results data seen below in the results section.

##  III Challenge Points

###  III-A Data Investigation

####  III-A1 Data Offset

Results from the base lab-trained model were promising and produced a very
high accuracy [5]. When it was discovered the model performed poorly on real-
world data, the first step was to perform a detailed examination of the
training data (Phase 1&2). It was discovered that several subjects had
systematic offsets between the IMU time and the labeling time. This caused the
predictions to be delayed by up to several seconds (see figure 2 and 3). Upon
investigation, it appears that the IMU clock drifts out of sync with the
camera clock over time. The degree of offset varies from subject to subject,
but typically is consistent across trials that occurred back-to-back.

This error was less pronounced on the original lab-trained model, likely
because the size of the sliding window was 25 frames (1 second) [5]. However,
when shrinking the size of the sliding window, the issue becomes much more
pronounced. The incorrect offsets were adjusted manually and trials without
the necessary information to manually fix were removed from the dataset. The
model was retrained on Phase 1&2 and re-evaluated on Phase 3 data. The model
showed a noticeable improvement in evaluation and the optimal size of the
sliding window changed from 25 frames to 10 frames.

####  III-A2 Incorrect Sensor Placement

Further into the dataset investigation process it was discovered that the back
sensor had been incorrectly placed on the subjects’ back in several phase 1&2
trials. Fortunately, the incorrect orientation was consistent, meaning that we
could use data from other sensors to infer if the sensor was incorrectly
placed and restructure the data mathematically to emulate the correct sensor
placement. The resulting sensor data was compared against a correctly placed
sensor on the same subject performing a similar motion and was confirmed to be
effective.

Typically, it would be advisable to remove any abnormal data. However, given
our very small dataset we decided to invest time in fixing the data instead of
removing it.

Figure 2: Typical lifting predictions on Phase 1&2 data Figure 3: Offset
Predictions on Phase 1&2 data. Model is consistently and systematically
predicting several seconds to early, implying a fault in the dataset.

###  III-B Non-Representative Training Data

Another challenge we faced porting this model onto real-world data was going
from specific laboratory-controlled movement to general real-world movement.
The Phase 1&2 non-lifting dataset consisted of primarily of sitting, standing,
walking, and jumping. In the real world, workers are performing much more
complex non-lifting motions. This causes the model to struggle to
differentiate between these complex movements and a lifting event. For
example, the model frequently predicts a lifting event when a worker is
utilizing a control panel.

Machine learning models learn their classifiers by identifying features that
exist when the classifier is present, but don’t exist when the classifier is
not. If the training data is not fully representative of the deployment
environment, the model may assume certain feature combinations always
correlate with the classifier when in fact they do not. In our situation, it
seems the model highly correlates lifting with a subject having their arms
extended out in front of them. This is a reasonable and true correlation in
Phase 1&2 but is not true in the more general Phase 3 data.

###  III-C Poor Model Generalization

As mentioned before, due to our relatively small and controlled training
dataset, our model struggled to interpret the complex motions in phase 3 data.
One technique we used to understand _how_ the model interpreted the data was
saliency mapping. Typically, neural network models are thought of as ”black
boxes”, where input data is fed into the model and outputs are generated. The
user has no concept of how the model processes the input data to provide the
outputs. This makes it very challenging to debug poorly performing models.
Saliency mapping is a technique which uses the weights between neurons as well
as an input to create a heatmap of the input features that have the most
significance on the output [9].

Historically, saliency maps are applied to images where the heatmap can be
overlaid on the input image. However, saliency maps can be used on non-image
data as well. For this project, we applied a saliency map to our IMU data. See
figure 4. Each row is a ’frame’ in our dataset (representing 1/25th of a
second) and each column is a feature from an IMU sensor. The darker squares
represent less contribution, while the lighter squares represent more. As we
can see in figure 4, this model uses a single point in our input data as a
major indicator of whether a lifting event is occurring or not.

Figure 4: IMU Saliency Map. Lighter colors: higher saliency, darker colors:
lower saliency. Columns: IMU channels (gyro or accelerometer) and rows: time
steps (1/25th of a second)

Unless we expect a single feature to effectively represent our class, this is
typically a bad sign as it implies the model did not generalize well. Here, we
can see that the upper arm gyroscope feature is the especially salient feature
in our input data. Seeing this should raise a key question ”Why is this such a
salient feature in our training dataset?” Upon further investigation, we
discovered that the only time a subject’s hands are out in front of them is
when they are lifting. Therefore, using the arm gyro is very effective on the
lab data, but as we know, there are frequently times in real-world situations
where the arms are extended in front of the body that aren’t lifting movements
(sitting at a desk, using machinery, etc.).

There are several effective ways to resolve this issue. The first and most
intuitive is to simply gather additional training data to address this edge
case. Unfortunately, gathering additional test data was not financially
feasible for this project, so we had to resort to more creative methods.

##  IV Methodology

With these three challenge-areas identified we investigated four potential
solutions to resolve these issues.

###  IV-A Synthetic IMU Data

One idea was to create additional data using lifting data gathered in other
public studies. This was originally deemed not feasible as there is no
standard placement of IMU sensors. Therefore, the data from other studies
utilizing IMUs would not match up with the data we gathered for our study.

However, many studies utilize Motion Capture data (mocap) to capture the
movements of subjects. The AMASS dataset [10] is a unified database of 15
MoCap datasets containing over 40 hours of motion data. BABEL [11] builds on
the AMASS dataset by labeling frames of each sequence with the action
occurring. Additionally, [12] outlines a methodology to synthesize IMU data
from the AMASS dataset by placing virtual IMU sensors on the mesh surface. The
combination of frame-by-frame data and synthetic IMU data generation makes
this dataset a feasible avenue for incorporating additional data.

###  IV-B Transfer Learning

Another idea was to use transfer learning. Transfer learning is where a model
is trained on a large, generalized dataset relating to its final classifier,
then is specifically trained on a small amount of data. A common example of
this is a voice assistant that responds only to a specific person. The voice
assistant is trained on general voice recognition, then later uses a few
samples of the user’s specific voice to identify that voice uniquely.

This idea doesn’t work well for IMU data because there are no standard IMU
placement locations. Models and datasets that are publicly available do not
match our IMU placements, and therefore are not fit for transfer learning in
the traditional sense. However, we can still implement a small transfer
learning style procedure in this model.

Instead of using external data to augment our dataset, we use our lab dataset
as the ”general” dataset, then a few examples of lifting and non-lifting from
a real-world environment. This solution wouldn’t fix the lapses in the lab
dataset, but it would help the model adjust to any personal bias in lifting
technique.

###  IV-C Filtering

Filtering data (or preprocessing) is a common machine learning technique. The
idea here was that filtering the data would mask the gaps in our training set
and force the model to learn more generally. To do this, we tested two types
of filters: Extended Kalman filters and Mahony filters.

Extended Kalman filters (EKFs) are a modification of standard Kalman filters,
both of which have seen widespread use in IMUs. Both operate fundamentally the
same way. Kalman filters first “predicts” the current state of the sensor
using a state transition model, then “updates” using a sensor reading, derived
from an observation model. This predict-update cycle allows for handling of
sensor noise. In particular, it can assist in the handling of gyroscope drift.
EKFs differ from conventional Kalman filters by using non-linear,
differentiable functions as part of the model, rather than linear models.

An EKF filter is readily available in the AHRS Python library. Mahony filters
(available in the same library) were also tested, despite being less widely
used, for a more comprehensive evaluation of the performance of filtering.

###  IV-D Sensor Removal

Removing sensors from the dataset seems very un-intuitive. We’d expect that
providing the model with more data would allow it to identify correlations
that we humans may not see, and while that can be true, in our situation, it
actually ends up hurting performance more than helping. With a small dataset,
we found that removing features not deemed salient to lifting (by referencing
NIOSH equation [7]) improved performance significantly.

This process must be done carefully, as removing data limits the capabilities
of our model. In our situation, our baseline accuracy was not acceptable for
production deployment, so we were willing to sacrifice some sensitivity to
gain specificity. It’s most important to identify high-risk lifts and since
high-risk lifts often have pronounced trunk flexion, we created a set of
sensors we felt best captured the motions of high-risk lifting. Models were
retrained on their specific sets of sensors and reevaluated on the entire
real-world dataset (not just high-risk events).

##  V Results and Discussion

All results are from the model evaluated on Phase 3 data. The ’Evaluation
data’ is a subset of the Phase 3 data.

###  V-A Synthetic IMU Data

The authors of [12] have provided their tooling for generating synthetic IMU
data. However, while the tooling provides a simulation of accelerometer data,
gyroscopic data is provided as orientations, requiring modifications to
provide angular velocities as in the NIOSH data set. Gyroscopic data is vital
in capturing key lifting characteristics; therefore, reconciling the two
formats is essential for incorporating additional data. Unfortunately, the
necessary modifications proved too far outside our area of experience to
incorporate any synthetic data into our training set.

###  V-B Transfer Learning

Our transfer learning strategy did not produce a significant improvement in
accuracy or F1. Typically, transfer learning models are trained on much more
data ahead of time and have very good performance on general class
detection/identification ahead of the transfer learning stage. In our process,
the model did not have a general understanding of what a lifting event was, so
finishing training with a small sample of real-world data did not improve
performance because the model was already poorly generalized. Simply put,
transfer learning works best from taking a powerful generalized model and
applying it to a specific sub-set of its original purpose (i.e., written digit
recognition). Our model was not a well-generalized model, so our specific
transfer training was unsuccessful.

###  V-C Filtering

See Figure 5. It was found that both EKF and Mahony filters resulted in
slightly lower median accuracy and F1 scores on training data, and notably
lower median accuracy and F1 scores on evaluation data. Furthermore, Figure 6
shows that even the models with the best F1 scores with filtering applied are
outperformed by the median model with no filtering.

It was noted in [5] that filtering the NIOSH data set was considered but
ultimately not used due to both better performance with raw data and reduced
pre-processing time; although it is unclear what filtering methods were
evaluated, the results we have obtained seem to align with these preliminary
findings.

One additional source of reduced performance from the EKF filter may be from
un-optimized filter parameters. Kalman filters incorporate noise variances of
the sensors. The default noise levels of the AHRS library were used in testing
due to unfamiliarity with the details of the sensors in the NIOSH data set.
Similarly, default values were used for gain parameters of the Mahony filter.

Figure 5: Results from filtering. No filtering provides the highest accuracy
and F1 scores on the evaluation dataset Figure 6: Heatmap results from
filtering. There are no models that perform exceptionally better than the mean
for any given filter technique

###  V-D Sensor Removal

Removing sensor channels produced excellent results. The best model trained
with only wrist and back sensors had an F1 score 8% higher than the top
performing model trained with all sensors. See figure 7 for a breakdown of the
different core channel combinations and their associated metrics.

Hiding data from a Machine Learning model typically does not produce better
results, however in this situation it did because the sensors we removed were
the _least_ representative of real-world lifting motions. In real-world
lifting environments, subjects are constantly moving their arms around as part
of the normal working environment, however in our training data, non-lifting
working motions were not included in the dataset. Therefore, when the model
learned our lifting class, it used features in the other sensors that were not
unique to lifting generally, but were in our training set. Removing these
sensor channels forced the model to learn a more generalized definition of a
lifting event from the remaining sensors.

Figure 7: Results from removing sensors. Using only the wrist and back sensors
provided the best F1 and accuracy scores on the Phase 3 dataset, showing
significant improvement over training with all sensors

HAR IMU models have come a long way since their original development. We can
design models to identify repetitive consistent laboratory-like data reliably
using multiple different model architectures. However, the vision of utilizing
a laboratory trained lifting model on real-world data is still largely out of
reach. Non-standard IMU placements cause limited datasets which make it
challenging to incorporate all reasonable movements in training. The authors
of [12] have successfully emulated IMU sensor data from motion-capture data by
placing virtual sensors on models, which could significantly expand the reach
of future training data sets. Furthermore, the authors have provided their
tooling for generating synthetic IMU data, although modifications may be
necessary to reconcile differences between synthetic and real-world data sets.

However, even with limited training data, we are still able to achieve
significant insight with saliency mapping techniques. We were able to apply
saliency mapping to our IMU model and see that the model was fixating on a
feature not vitally significant to lifting. From this learning, we were able
to target our solutions to resolve this challenge and achieved significant
results for our efforts.

Our baseline balanced Phase 3 F1 was 32.8%. By perfecting training data,
resolving sensor issues, and performing a targeted removal of non-salient
sensors, we were able to increase our typical F1 scores up to 55%, with our
top model nearly doubling our baseline with a F1 of 64.66%

Unfortunately, even with these substantial increases to performance, we still
fall short of the 90+% necessary for deployment. However, with more work in
synthetic data future work could improve performance enough to make this
solution viable for deployment.

##  VI Conclusion

Outlined below are three key research considerations informed from this
project:

  1. 1.

Training Data is very important

     * •

Thoroughly ensure the accuracy of your training set and insure it is
representative of your deployment environment. If the dataset isn’t
representative, it will be very challenging to develop a successful model.

     * •

Researching ways to augment datasets: Larger training sets, as long as they’re
representative, will help your model generalize better. Augmenting data with
synthetic data or publicly available datasets is a powerful way to increase
your model’s robustness.

  2. 2.

Studying model’s saliency: Having a sense of how your model interprets data is
vital to understanding how it will behave in new environments. It can also
inform modifications to make ahead of time to improve the chances of success.

  3. 3.

Considering transfer learning: Transfer learning is a great way to adjust a
model slightly to remove user bias. While transfer learning was unsuccessful
on this particular project, transfer learning can provide a ”personal touch”
to a model that can significantly improve its reliability on a specific user.

##  VII Acknowledgments

This work was funded by CDC/NIOSH. The authors would like to thank Marie
Hayden, Menekse Barim and Dwight Werren for their support in creating the
datasets. Findings and conclusions in this report are those of the authors and
do not necessarily represent the official positions of NIOSH or the Centers
for Disease Control and Prevention (CDC). The mention of any company or
product does not constitute endorsement by NIOSH or the CDC.

## References

  * [1] V. Anderson, B. Bernard, S. Burt, L. Cole, C. Estill, L. Fine, K. Grant, C. Gjessing, L. Jenkins, J. Hurrell, N. Nelson, D. Pfirman, R. Roberts, D. Stetson, M. Haring-Sweeney, and S. Tanaka, _Musculoskeletal Disorders and Workplace Factors: A Critical Review of Epidemiologic Evidence for Work-Related Musculoskeletal Disorders of the Neck, Upper Extremity, and Low Back_ , 01 1997.
  * [2] _Work practices guide for manual lifting._ , ser. DHHS (NIOSH) publication ; no. 81-122. Cincinnati, Ohio: U.S. Dept. of Health and Human Services, Public Health Service, Centers for Disease Control, National Institute for Occupational Safety and Health, Division of Biomedical and Behavioral Science, 1981.
  * [3] D. Cook, K. D. Feuz, and N. C. Krishnan, “Transfer learning for activity recognition: a survey,” _Knowledge and Information Systems_ , vol. 36, no. 3, pp. 537–556, Sep 2013. [Online]. Available: https://doi.org/10.1007/s10115-013-0665-3
  * [4] I. Rodríguez-Moreno, J. M. Martínez-Otzeta, B. Sierra, I. Rodriguez, and E. Jauregi, “Video activity recognition: State-of-the-art,” _Sensors_ , vol. 19, no. 14, 2019. [Online]. Available: https://www.mdpi.com/1424-8220/19/14/3160
  * [5] B. Thomas, M.-L. Lu, R. Jha, and J. Bertrand, “Machine learning for detection and risk assessment of lifting action,” _IEEE Transactions on Human-Machine Systems_ , vol. 52, no. 6, pp. 1196–1204, 2022.
  * [6] X. Ying, “An overview of overfitting and its solutions,” _Journal of Physics: Conference Series_ , vol. 1168, no. 2, p. 022022, feb 2019. [Online]. Available: https://dx.doi.org/10.1088/1742-6596/1168/2/022022
  * [7] L. Donisi, G. Cesarelli, A. Coccia, M. Panigazzi, E. M. Capodaglio, and G. D’Addio, “Work-related risk assessment according to the revised niosh lifting equation: A preliminary study using a wearable inertial sensor and machine learning,” _Sensors_ , vol. 21, no. 8, 2021. [Online]. Available: https://www.mdpi.com/1424-8220/21/8/2593
  * [8] M. Barim, M.-L. Lu, S. Feng, G. Hughes, M. Hayden, and D. Werren, “Accuracy of an algorithm using motion data of five wearable imu sensors for estimating lifting duration and lifting risk factors,” _Proceedings of the Human Factors and Ergonomics Society Annual Meeting_ , vol. 63, no. 1, p. 1105–1111, 2019.
  * [9] K. Simonyan, A. Vedaldi, and A. Zisserman, “Deep inside convolutional networks: Visualising image classification models and saliency maps,” 2014.
  * [10] N. Mahmood, N. Ghorbani, N. F. Troje, G. Pons-Moll, and M. J. Black, “AMASS: Archive of motion capture as surface shapes,” in _2019 IEEE/CVF International Conference on Computer Vision (ICCV)_ , Oct. 2019, pp. 5441–5450.
  * [11] A. R. Punnakkal, A. Chandrasekaran, N. Athanasiou, A. Quiros-Ramirez, and M. J. Black, “BABEL: Bodies, action and behavior with english labels,” in _Proceedings IEEE/CVF Conf. on Computer Vision and Pattern Recognition (CVPR)_ , Jun. 2021, pp. 722–731.
  * [12] Y. Huang, M. Kaufmann, E. Aksan, M. J. Black, O. Hilliges, and G. Pons-Moll, “Deep inertial poser learning to reconstruct human pose from sparseinertial measurements in real time,” _ACM Transactions on Graphics, (Proc. SIGGRAPH Asia)_ , vol. 37, no. 6, pp. 185:1–185:15, Nov. 2018.