Package 'HMMpa'

Analysing Accelerometer Data Using Hidden Markov Models

Description

This package provides functions for analyzing accelerometer outpout data (known as a time-series of (impulse)-counts) to quantify length and intensity of physical activity.

Usually, so called activity ranges are used to classify an activity as “sedentary”, “moderate” and so on. Activity ranges are separated by certain thresholds (cut-off points). The choice of these cut-off points depends on different components like the subjects' age or the type of accelerometer device.

Cut-off point values and defined activity ranges are important input values of the following analyzing tools provided by this package:

1. Cut-off point method (assigns an activity range to a count given its total magintude). This traditional approach assigns an activity range to each count of the time-series independently of each other given its total magnitude.

2. HMM-based method (assigns an activity range to a count given its underlying PA-level). This approach uses a stochastic model (the hidden Markov model or HMM) to identify the (Markov dependent) time-series of physical activity states underlying the given time-series of accelerometer counts. In contrast to the cut-off point method, this approach assigns activity ranges to the estimated PA-levels corresponding to the hidden activity states, and not directly to the accelerometer counts.

Details

Package:	HMMpa
Type:	Package
Version:	1.0.1
Date:	2018-04-20
License:	GPL-3

The new procedure for analyzing accelerometer data can be roughly described as follows:
First, a hidden Markov model (HMM) is trained to estimate the number m of hidden physical activity states and the model specific parameters (delta, gamma, distribution_theta). Then, a user-sepcified decoding algorithm decodes the trainded HMM to classify each accelerometer count into the m hidden physical activity states. Finally, the estimated distribution mean values (PA-levels) corresponding to the hidden physical activity states are extracted and the accelerometer counts are assigned by the total magnitudes of their corresponding PA-levels to given physical activity ranges (e.g. "sedentary", "light", "moderate" and "vigorous") by the traditional cut-off point method.

Note

We thank Moritz Hanke for his help in realizing this package.

Author(s)

Vitali Witowski,
Ronja Foraita, Leibniz Institute for Prevention Research and Epidemiology (BIPS)

Maintainer: Ronja Foraita <[email protected]>

References

Baum, L., Petrie, T., Soules, G., Weiss, N. (1970). A maximization technique occurring in the statistical analysis of probabilistic functions of markov chains. The annals of mathematical statistics, vol. 41(1), 164–171.

Brachmann, B. (2011). Hidden-Markov-Modelle fuer Akzelerometerdaten. Diploma Thesis, University Bremen - Bremen Institute for Prevention Research and Social Medicine (BIPS).

Dempster, A., Laird, N., Rubin, D. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society. Series B (Methodological), vol. 39(1), 1–38.

Forney, G.D. (1973). The Viterbi algorithm. Proceeding of the IEE, vol. 61(3), 268–278.

Joe, H., Zhu, R. (2005). Generalized poisson distribution: the property of mixture of poisson and comparison with negative binomial distribution. Biometrical Journal, vol. 47(2), 219–229.

MacDonald, I. L., Zucchini, W. (2009) Hidden Markov Models for Time Series: An Introduction Using R, Boca Raton: Chapman & Hall.

Viterbi, A.J. (1967). Error Bounds for concolutional codes and an asymptotically optimal decoding algorithm. Information Theory, IEEE Transactions on, vol. 13(2), 260–269.

Witowski, V., Foraita, R., Pitsiladis, Y., Pigeot, I., Wirsik, N. (2014) Using hidden Markov models to improve quantifying physical activity in accelerometer data - A simulation study. PLOS ONE. 9(12), e114089. http://dx.doi.org/10.1371/journal.pone.0114089

Examples

################################################################
###     Example 1  (traditional approach)    ################### 
###     Solution of the cut-off point method ################### 
################################################################

################################################################
### Fictitious activity counts #################################
################################################################

x <- c(1,16,19,34,22,6,3,5,6,3,4,1,4,3,5,7,9,8,11,11,
  14,16,13,11,11,10,12,19,23,25,24,23,20,21,22,22,18,7,
  5,3,4,3,2,3,4,5,4,2,1,3,4,5,4,5,3,5,6,4,3,6,4,8,9,12,
  9,14,17,15,25,23,25,35,29,36,34,36,29,41,42,39,40,43,
  37,36,20,20,21,22,23,26,27,28,25,28,24,21,25,21,20,21,
  11,18,19,20,21,13,19,18,20,7,18,8,15,17,16,13,10,4,9,
  7,8,10,9,11,9,11,10,12,12,5,13,4,6,6,13,8,9,10,13,13,
  11,10,5,3,3,4,9,6,8,3,5,3,2,2,1,3,5,11,2,3,5,6,9,8,5,
  2,5,3,4,6,4,8,15,12,16,20,18,23,18,19,24,23,24,21,26,
  36,38,37,39,45,42,41,37,38,38,35,37,35,31,32,30,20,39,
  40,33,32,35,34,36,34,32,33,27,28,25,22,17,18,16,10,9,
  5,12,7,8,8,9,19,21,24,20,23,19,17,18,17,22,11,12,3,9,
  10,4,5,13,3,5,6,3,5,4,2,5,1,2,4,4,3,2,1)  
    
traditional_cut_off_point_method <- cut_off_point_method(x=x, 
 	cut_points = c(5,15,23), 
 	names_activity_ranges = c("SED","LIG","MOD","VIG"), 
 	bout_lengths = c(1,1,2,4,5,10,11,20,21,60,61,260),
 	plotting = 1)



################################################################ 
################################################################ 
###      Examples 2,3 and 4  (new approach)             ########
###      Solution of the HMM based cut-off point method ######## 
################################################################
###      Demonstrated both in three steps (Example 2)     ###### 
###      and condensed in one function (Examples 3 and 4) ######
################################################################

################################################################
### Example 2) Manually in three steps    ######################
################################################################

################################################################
### Fictitious activity counts #################################
################################################################
                                                                
x <- c(1,16,19,34,22,6,3,5,6,3,4,1,4,3,5,7,9,8,11,11,           
  14,16,13,11,11,10,12,19,23,25,24,23,20,21,22,22,18,7,         
  5,3,4,3,2,3,4,5,4,2,1,3,4,5,4,5,3,5,6,4,3,6,4,8,9,12,         
  9,14,17,15,25,23,25,35,29,36,34,36,29,41,42,39,40,43,         
  37,36,20,20,21,22,23,26,27,28,25,28,24,21,25,21,20,21,        
  11,18,19,20,21,13,19,18,20,7,18,8,15,17,16,13,10,4,9,         
  7,8,10,9,11,9,11,10,12,12,5,13,4,6,6,13,8,9,10,13,13,         
  11,10,5,3,3,4,9,6,8,3,5,3,2,2,1,3,5,11,2,3,5,6,9,8,5,         
  2,5,3,4,6,4,8,15,12,16,20,18,23,18,19,24,23,24,21,26,         
  36,38,37,39,45,42,41,37,38,38,35,37,35,31,32,30,20,39,        
  40,33,32,35,34,36,34,32,33,27,28,25,22,17,18,16,10,9,         
  5,12,7,8,8,9,19,21,24,20,23,19,17,18,17,22,11,12,3,9,         
  10,4,5,13,3,5,6,3,5,4,2,5,1,2,4,4,3,2,1)                      
     
################################################################
## Step 1: Training of a HMM ##################################
##         for the given time-series of counts ################ 
##         ####################################################
##         More precisely: training of a poisson ##############
##         distribution based hidden Markov model for #########
##         number of states m=2,...,6 #########################
##         and selection of the model with the most ########### 
##         plausible m ########################################
################################################################

m_trained_HMM <- HMM_training(x = x, 
 min_m = 2, 
 max_m = 6, 
 distribution_class = "pois")$trained_HMM_with_selected_m  
 	 
###############################################################
## Step 2: Decoding of the trained HMM for the given ##########
##         time-series of accelerometer counts to extract #####
##         hidden PA-levels ###################################
###############################################################
hidden_PA_levels <- HMM_decoding(x = x, 
 	 m = m_trained_HMM$m, 
 	 delta = m_trained_HMM$delta, 
 	 gamma = m_trained_HMM$gamma, 
 	 distribution_class = m_trained_HMM
 	 $distribution_class, 
 	 distribution_theta = m_trained_HMM$
 	 distribution_theta)$decoding_distr_means

############################################################### 
## Step 3: Assigning of user-sepcified activity ranges ########
##         to the accelerometer counts via the total ##########
##         magnitudes of their corresponding ##################
##         hidden PA-levels ###################################
##         ####################################################
##         In this example four activity ranges ###############
##         (named as "sedentary", "light", "moderate" #########
##         and "vigorous" physical activity) are ##############
##         separated by the three cut-points 5, 15 and 23) ####
################################################################
HMM_based_cut_off_point_method <- cut_off_point_method(x = x, 
 	 hidden_PA_levels = hidden_PA_levels, 
 	 cut_points = c(5,15,23), 
 	 names_activity_ranges = c("SED","LIG","MOD","VIG"), 
 	 bout_lengths = c(1,1,2,4,5,10,11,20,21,60,61,260),
 	 plotting=1)      				 




###############################################################
## Example 3) In a single function (Step 1-3 of Example 2  ####
##            combined in one function)                    ####
###############################################################
	
###############################################################
## Fictitious activity counts #################################
###############################################################
                                                                
x <- c(1,16,19,34,22,6,3,5,6,3,4,1,4,3,5,7,9,8,11,11,           
  14,16,13,11,11,10,12,19,23,25,24,23,20,21,22,22,18,7,         
  5,3,4,3,2,3,4,5,4,2,1,3,4,5,4,5,3,5,6,4,3,6,4,8,9,12,         
  9,14,17,15,25,23,25,35,29,36,34,36,29,41,42,39,40,43,         
  37,36,20,20,21,22,23,26,27,28,25,28,24,21,25,21,20,21,        
  11,18,19,20,21,13,19,18,20,7,18,8,15,17,16,13,10,4,9,         
  7,8,10,9,11,9,11,10,12,12,5,13,4,6,6,13,8,9,10,13,13,         
  11,10,5,3,3,4,9,6,8,3,5,3,2,2,1,3,5,11,2,3,5,6,9,8,5,         
  2,5,3,4,6,4,8,15,12,16,20,18,23,18,19,24,23,24,21,26,         
  36,38,37,39,45,42,41,37,38,38,35,37,35,31,32,30,20,39,        
  40,33,32,35,34,36,34,32,33,27,28,25,22,17,18,16,10,9,         
  5,12,7,8,8,9,19,21,24,20,23,19,17,18,17,22,11,12,3,9,         
  10,4,5,13,3,5,6,3,5,4,2,5,1,2,4,4,3,2,1)                      

###############################################################
## Use a (m=4 state) hidden Markov model based on the #########
## generalized poisson distribution to assign an      #########
## activity range to the counts                       #########
###############################################################
## In this example three activity ranges                   ####
## (named as "light", "moderate" and "vigorous" physical   ####
## activity) are separated by the two cut-points 15 and 23 ####
###############################################################

HMM_based_cut_off_point_method <- HMM_based_method(x = x, 
 	 cut_points = c(15,23), 
 	 min_m = 4, 
 	 max_m = 4, 
 	 names_activity_ranges = c("LIG","MOD","VIG"), 
 	 distribution_class = "genpois", 
 	 training_method = "numerical",
 	 DNM_limit_accuracy = 0.05,
 	 DNM_max_iter = 10,
 	 bout_lengths = c(1,1,2,4,5,10,11,20,21,60,61,260),
 	 plotting = 1)
 	 

###############################################################
## Example 4) In a single function (Step 1-3 of Example 2  ####
##            combined in one function)                    ####
##            (large and highly scatterd time-series)      ####
###############################################################
	 
################################################################
### Generate a large time-series of highly scattered counts ####
################################################################

x <- HMM_simulation(
 size = 1500, 
 m = 10,
 gamma = 0.93 * diag(10) + rep(0.07 / 10, times = 10),
 distribution_class = "norm", 
 distribution_theta = list(mean = c(10, 100, 200, 300, 450, 
 600, 700, 900, 1100, 1300, 1500), 
 sd=c(rep(100,times=10))), 
 obs_round=TRUE, 
 obs_non_neg=TRUE,
 plotting=5)$observations

################################################################
### Compare results of the tradional cut-point method ##########
### and the (6-state-normal-)HMM based method ##################
################################################################
		
traditional_cut_off_point_method <- cut_off_point_method(x=x, 
 	 cut_points = c(200,500,1000), 
 	 names_activity_ranges = c("SED","LIG","MOD","VIG"), 
 	 bout_lengths = c(1,1,2,4,5,10,11,20,21,60,61,260),
 	 plotting = 1)

HMM_based_cut_off_point_method <- HMM_based_method(x=x,
	 max_scaled_x = 200, 
 	 cut_points = c(200,500,1000), 
 	 min_m = 6, 
 	 max_m = 6,
 	 BW_limit_accuracy = 0.5, 
 	 BW_max_iter = 10,
 	 names_activity_ranges = c("SED","LIG","MOD","VIG"), 
 	 distribution_class = "norm", 
 	 bout_lengths = c(1,1,2,4,5,10,11,20,21,60,61,260),
 	 plotting = 1)

---

AIC and BIC Value for a Discrete Time Hidden Markov Model

Description

Computes the Aikaike's information criterion and the Bayesian information criterion for a discrete time hidden Markov model, given a time-series of observations.

Usage

AIC_HMM(logL, m, k)
BIC_HMM(size, m, k, logL)

Arguments

size	length of the time-series of observations x (also T).
m	number of states in the Markov chain of the model.
k	single numeric value representing the number of parameters of the underlying distribution of the observation process (e.g. k=2 for the normal distribution (mean and standard deviation)).
logL	logarithmized likelihood of the model.

Details

For a discrete-time hidden Markov model, AIC and BIC are as follows (MacDonald & Zucchini (2009, Paragraph 6.1 and A.2.3)):

$\mbox{AIC} = -2 logL + 2p$

$\mbox{BIC} = -2 logL + p \log T$

where T indicates the length/size of the observation time-series and p denotes the number of independent parameters of the model. In case of a HMM as provided by this package, where k and m are defined as in the arguments above, p can be calculated as follows:

$p = m^2+km-1.$

Value

The AIC or BIC value of the fitted hidden Markov model.

Author(s)

Based on MacDonald & Zucchini (2009, Paragraph 6.1 and A.2.3). Implementation by Vitali Witowski (2013).

References

MacDonald, I. L., Zucchini, W. (2009) Hidden Markov Models for Time Series: An Introduction Using R, Boca Raton: Chapman & Hall.

See Also

HMM_training

Examples

################################################################
### Fictitious observations ####################################
################################################################

x <- c(1,16,19,34,22,6,3,5,6,3,4,1,4,3,5,7,9,8,11,11,
  14,16,13,11,11,10,12,19,23,25,24,23,20,21,22,22,18,7,
  5,3,4,3,2,3,4,5,4,2,1,3,4,5,4,5,3,5,6,4,3,6,4,8,9,12,
  9,14,17,15,25,23,25,35,29,36,34,36,29,41,42,39,40,43,
  37,36,20,20,21,22,23,26,27,28,25,28,24,21,25,21,20,21,
  11,18,19,20,21,13,19,18,20,7,18,8,15,17,16,13,10,4,9,
  7,8,10,9,11,9,11,10,12,12,5,13,4,6,6,13,8,9,10,13,13,
  11,10,5,3,3,4,9,6,8,3,5,3,2,2,1,3,5,11,2,3,5,6,9,8,5,
  2,5,3,4,6,4,8,15,12,16,20,18,23,18,19,24,23,24,21,26,
  36,38,37,39,45,42,41,37,38,38,35,37,35,31,32,30,20,39,
  40,33,32,35,34,36,34,32,33,27,28,25,22,17,18,16,10,9,
  5,12,7,8,8,9,19,21,24,20,23,19,17,18,17,22,11,12,3,9,
  10,4,5,13,3,5,6,3,5,4,2,5,1,2,4,4,3,2,1) 

### Assummptions (probability vector, transition matrix, 
### and distribution parameters)

delta <- c(0.25,0.25,0.25,0.25)

gamma <- 0.7 * diag(length(delta)) + rep(0.3 / length(delta))

distribution_class <- "pois"

distribution_theta <- list(lambda = c(4,9,17,25))


### log-likelihood

logL <- forward_backward_algorithm (x = x, 
                delta = delta, gamma=gamma,  
                distribution_class= distribution_class, 
                distribution_theta=distribution_theta)$logL


### the Poisson distribution has one paramter, hence k=1

AIC_HMM(logL = logL, m = length(delta), k = 1)

BIC_HMM(size = length(x) , logL = logL, m = length(delta), k = 1)

---

Estimation Using the Baum-Welch Algorithm

Description

Estimates the parameters of a (non-stationary) discrete-time hidden Markov model. The Baum-Welch algorithm is a version of the EM (Estimation/Maximization) algorithm. See MacDonald & Zucchini (2009, Paragraph 4.2) for further details.

Usage

Baum_Welch_algorithm(x, m, delta, gamma, distribution_class, 
                     distribution_theta, discr_logL = FALSE, 
                     discr_logL_eps = 0.5,
                     BW_max_iter = 50, BW_limit_accuracy = 0.001,
                     BW_print = TRUE, Mstep_numerical = FALSE, 
                     DNM_limit_accuracy = 0.001, DNM_max_iter = 50, 
                     DNM_print = 2)

Arguments

x	a vector object containing the time-series of observations that are assumed to be realizations of the (hidden Markov state dependent) observation process of the model.
m	a (finite) number of states in the hidden Markov chain.
delta	a vector object containing starting values for the marginal probability distribution of the m states of the Markov chain at the time point t=1 for the Baum-Welch algorithm.
gamma	a matrix (ncol=nrow=m) containing starting values for the transition matrix of the hidden Markov chain.
distribution_class	a single character string object with the abbreviated name of the m observation distributions of the Markov dependent observation process. The following distributions are supported: Poisson (pois); generalized Poisson (genpois, parameter estimation via the Baum-Welch algorithm is only supported if the M-step is performed numerically, i.e. if Mstep_numerical = TRUE); normal (norm).
distribution_theta	a list object containing starting values for the parameters of the m observation distributions of the observation process that are dependent on the hidden Markov state.
discr_logL	a logical object indicating whether the discrete log-likelihood should be used (for distribution_class="norm") for estimating the model specific parameters instead of the general log-likelihood. See MacDonald & Zucchini (2009, Paragraph 1.2.3) for further details. Default value is FALSE.
discr_logL_eps	a single numerical value to approximately determine the discrete likelihood for a hidden Markov model based on nomal distributions (for "norm"). Default value is 0.5. See MacDonald & Zucchini (2009, Paragraph 1.2.3) for further details.
BW_max_iter	a single numerical value representing the maximum number of iterations in the Baum-Welch algorithm. Default value is 50.
BW_limit_accuracy	a single numerical value representing the convergence criterion of the Baum-Welch algorithm. Default value is 0.001.
BW_print	a logical object indicating whether the log-likelihood at each iteration-step shall be printed. Default value is TRUE.
Mstep_numerical	a logical object indicating whether the Maximization Step of the Baum-Welch algorithm shall be performed by numerical maximization using the nlm-function. Default value is FALSE.
DNM_limit_accuracy	a single numerical value representing the convergence criterion of the numerical maximization algorithm using the nlm-function (used to perform the M-step of the Baum-Welch-algorithm). Default value is 0.001.
DNM_max_iter	a single numerical value representing the maximum number of iterations of the numerical maximization using the nlm-function (used to perform the M-step of the Baum-Welch-algorithm). Default value is 50.
DNM_print	a single numerical value to determine the level of printing of the nlm-function. See nlm-function for further informations. The value 0 suppresses, that no printing will be outputted. Default value is 2 for full printing.

Value

Baum_Welch_algorithm returns a list containing the estimated parameters of the hidden Markov model and other components. See MacDonald & Zucchini (2009, Paragraph 4.2) for further details on the calculated objects within this algorithm.

x	input time-series of observations.
m	input number of hidden states in the Markov chain.
zeta	a (T,m)-matrix (when T indicates the length/size of the observation time-series and m the number of states of the HMM) containing probabilities (estimates of the conditional expectations of the missing data given the observations and the estimated model specific parameters) calculated by the algorithm. See MacDonald & Zucchini (2009, Paragraph 4.2.2) for further details.
eta	a (T,m,m)-dimensional-array (when T indicates the length of the observation time-series and m the number of states of the HMM) containing probabilities (estimates of the conditional expectations of the missing data given the observations and the estimated model specific parameters) calculated by the algorithm. See MacDonald & Zucchini (2009, Paragraph 4.2.2) for further details.
logL	a numerical value representing the logarithmized likelihood calculated by the forward_backward_algorithm.
iter	number of performed iterations.
BIC	a numerical value representing the Bayesian information criterion for the hidden Markov model with estimated parameters.
delta	a vector object containing the estimates for the marginal probability distribution of the m states of the Markov chain at time-point point t=1.
gamma	a matrix containing the estimates for the transition matrix of the hidden Markov chain.
...	other input values (as arguments above). In the case that the algorithm stops before the targeted accuracy or the maximum number of iterations has been reached, further values are displayed and the estimates from the last successful iteration step are saved.

Author(s)

The basic algorithm for a Poisson-HMM is provided by MacDonald & Zucchini (2009, Paragraph 4.2, Paragraph A.2.3). Extension and implementation by Vitali Witowski (2013).

References

Baum, L., Petrie, T., Soules, G., Weiss, N. (1970). A maximization technique occurring in the statistical analysis of probabilistic functions of markov chains. The annals of mathematical statistics, vol. 41(1), 164–171.

Dempster, A., Laird, N., Rubin, D. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society. Series B (Methodological), vol. 39(1), 1–38.

MacDonald, I. L., Zucchini, W. (2009) Hidden Markov Models for Time Series: An Introduction Using R, Boca Raton: Chapman & Hall.

See Also

HMM_based_method, HMM_training, direct_numerical_maximization, forward_backward_algorithm, initial_parameter_training

Examples

################################################################
### Fictitious observations ####################################
################################################################

x <- c(1,16,19,34,22,6,3,5,6,3,4,1,4,3,5,7,9,8,11,11,
  14,16,13,11,11,10,12,19,23,25,24,23,20,21,22,22,18,7,
  5,3,4,3,2,3,4,5,4,2,1,3,4,5,4,5,3,5,6,4,3,6,4,8,9,12,
  9,14,17,15,25,23,25,35,29,36,34,36,29,41,42,39,40,43,
  37,36,20,20,21,22,23,26,27,28,25,28,24,21,25,21,20,21,
  11,18,19,20,21,13,19,18,20,7,18,8,15,17,16,13,10,4,9,
  7,8,10,9,11,9,11,10,12,12,5,13,4,6,6,13,8,9,10,13,13,
  11,10,5,3,3,4,9,6,8,3,5,3,2,2,1,3,5,11,2,3,5,6,9,8,5,
  2,5,3,4,6,4,8,15,12,16,20,18,23,18,19,24,23,24,21,26,
  36,38,37,39,45,42,41,37,38,38,35,37,35,31,32,30,20,39,
  40,33,32,35,34,36,34,32,33,27,28,25,22,17,18,16,10,9,
  5,12,7,8,8,9,19,21,24,20,23,19,17,18,17,22,11,12,3,9,
  10,4,5,13,3,5,6,3,5,4,2,5,1,2,4,4,3,2,1)  


### Assummptions (number of states, probability vector, 
### transition matrix, and distribution parameters)

m <- 4

delta <- c(0.25,0.25,0.25,0.25)

gamma <- 0.7 * diag(m) + rep(0.3 / m)

distribution_class <- "pois"

distribution_theta <- list(lambda = c(4,9,17,25))


### Estimation of a HMM using the Baum-Welch algorithm

trained_HMM_with_m_hidden_states <- 
    Baum_Welch_algorithm(x = x, 
        m = m, 
        delta = delta, 
        gamma = gamma,
        distribution_class = distribution_class, 
        distribution_theta = distribution_theta)

print(trained_HMM_with_m_hidden_states)

---

Cut-Off Point Method for Assigning Physical Activity Patterns

Description

This function assigns an activity range to each observation of a time-series, such as for a sequence of impulse counts recorded by an accelerometer. The activity ranges are defined by thresholds called “cut-off points”. Furthermore, bout periods are analysed (see Details for further informations).

Usage

cut_off_point_method (x, cut_points, 
        names_activity_ranges = NA, hidden_PA_levels = NA, 
        bout_lengths = NULL, plotting = 0)

Arguments

x	a vector object of length T containing non-negative observations of a time-series, such as a sequence of accelerometer impulse counts.
cut_points	a vector object containing cut-off points to separate activity ranges. For instance, the vector c(7,15,23) separates the four activity ranges [0,7);[7,15);[15,23);[23,Inf).
names_activity_ranges	an optional character string vector to name the activity ranges induced by the cut-points. This vector must contain one element more than the vector cut_points.
bout_lengths	a vector object (with even number of elemets) to define the range of the bout intervals (see Details for the definition of bouts). For instance, bout_lengths=c(1,1,2,2,3,10,11,20,1,20) defines the five bout intervals [1,1] (1 count); [2,2] (2 counts); [3,10] (3-10 counts); [11,20] (11-20 counts); [1,20] (1-20 counts - overlapping with other bout intervalls is possible). Default value is bout_lengths=NULL.
hidden_PA_levels	an optional vector object of length T containing a sequence of the estimated hidden physical activity levels (i.e. means) underlying the time-series of accelerometer counts. Such a sequence can be extracted by decoding a trained hidden Markov model. The cut-point method classifies then each count by its level in the hidden Markov chain that generates the physical activity counts, and does not use the observed count value (see HMM_based_method for further details). Default is NA (for the traditional cut-point method).
plotting	a numeric value between 0 and 5 (generates different outputs). NA suppresses graphical output. Default value is 0. 0: output 1-5 1: summary of all results 2: time series of activity counts, classified into activity ranges 3: time series of bouts (and, if available, the sequence of the estimated hidden physical activity levels, extracted by decoding a trained HMM, in green colour) 4: barplots of absolute and relative frequencies of time spent in different activity ranges 5: barplots of absolute frequencies of different bout intervals (overall and by activity ranges )

Details

A bout is defined as a period of time spending a defined intensity of physical activities in a specified physical activity range, without switching to activity intensities in a different activity range.

Value

cut_off_point_method returns a list containing the extracted sequence of activity ranges and plots key figures.

activity_ranges	an array object containing the cut-off intervals that indicate the activity ranges.
classification	an integer vector containing the sequence of activity ranges that were assigned to the observed time-series of accelerometer counts. If hidden_PA_levels=NA, then classification is the output of the traditional cut-point method, meaning that an activity range has been assigned to each accelerometer count over its observed value actual position. In case when hidden_PA_levels is available, classification is the output of the extendend cut-point method using hidden Markov models (see HMM_based_method for further details).
classification_per_activity_range	a pairlist object containing the classification of the observed counts by the assigned activity range.
freq_acitvity_range	table object containing the absolute frequencies of classifications into activity ranges.
rel_freq_acitvity_range	table object containing the relative frequencies of classifications into activity ranges.
quantity_of_bouts	overall number of bouts.
bout_periods	an array including the bout length assigned to acitiy ranges.
abs_freq_bouts_el	a pairlist object containing the absolute frequency of bout length per epoch length (aggregated).

Author(s)

Vitali Witowski (2013).

See Also

HMM_based_method

Examples

################################################################
### Fictitious activity counts #################################
################################################################

x <- c(1,16,19,34,22,6,3,5,6,3,4,1,4,3,5,7,9,8,11,11,
  14,16,13,11,11,10,12,19,23,25,24,23,20,21,22,22,18,7,
  5,3,4,3,2,3,4,5,4,2,1,3,4,5,4,5,3,5,6,4,3,6,4,8,9,12,
  9,14,17,15,25,23,25,35,29,36,34,36,29,41,42,39,40,43,
  37,36,20,20,21,22,23,26,27,28,25,28,24,21,25,21,20,21,
  11,18,19,20,21,13,19,18,20,7,18,8,15,17,16,13,10,4,9,
  7,8,10,9,11,9,11,10,12,12,5,13,4,6,6,13,8,9,10,13,13,
  11,10,5,3,3,4,9,6,8,3,5,3,2,2,1,3,5,11,2,3,5,6,9,8,5,
  2,5,3,4,6,4,8,15,12,16,20,18,23,18,19,24,23,24,21,26,
  36,38,37,39,45,42,41,37,38,38,35,37,35,31,32,30,20,39,
  40,33,32,35,34,36,34,32,33,27,28,25,22,17,18,16,10,9,
  5,12,7,8,8,9,19,21,24,20,23,19,17,18,17,22,11,12,3,9,
  10,4,5,13,3,5,6,3,5,4,2,5,1,2,4,4,3,2,1)  


################################################################
### i.) Tradionional cut_point method ##########################
################################################################

### Assigning activity ranges to activity counts using 
### fictitious cut-off points that produce the four activity 
### ranges "sedentary"", "light"", "moderate"", and "vigorous". 

solution_of_traditional_cut_off_point_method <- 
    cut_off_point_method(x = x, 
      cut_points = c(5,15,23), 
      names_activity_ranges = c("SED","LIG","MOD","VIG"), 
      bout_lengths = c(1,1,2,2,3,3,4,4,5,5,6,12,
      13,40,41,265,1,265), 
      plotting = 0)

print(solution_of_traditional_cut_off_point_method)


###############################################################
### ii.) Extension of the tradionional cut_point method #######
###      using HMMs      ######################################
###############################################################

## The following three steps define an extension of the 
## traditional cut-off method by first extracting the hidden 
## physical activity pattern behind the accelerometer counts 
## using a HMM (those three steps are basically combined in 
## the function HMM_based_method, see HMM_based_method for 
## further details and references): 


### Step 1 ##################################################### 

## Train hidden Markov model for different number of 
## states m=2,...,6 and select the model with the most 
## plausible m

m_trained_HMM <- 
    HMM_training(x = x, 
      min_m = 2, 
      max_m = 6, BW_print=FALSE,
      distribution_class = "pois")$trained_HMM_with_selected_m


### Step 2 ##################################################### 

## Decode the trained HMM (by using the 
## Viterbi algorithm (global decoding)) to get the estimated 
## sequence of hidden physical activity levels 
## underlying the the accelerometer counts 

## You have to compute 'm_trained_HMM' first (see Step 1)

global_decoding <- 
    HMM_decoding(x = x, 
      m = m_trained_HMM$m, 
      delta = m_trained_HMM$delta, 
      gamma = m_trained_HMM$gamma, 
      distribution_class = m_trained_HMM$distribution_class, 
      distribution_theta = m_trained_HMM$distribution_theta,
      decoding_method = "global")
        
hidden_PA_levels <- 
    global_decoding$decoding_distr_means


### Step 3 #####################################################

## Assigning activity ranges to activity counts using the 
## information extracted by decoding the HMM for the counts 
## (PA-levels) and fictitious cut-off points that produce 
## four so-called activity ranges:"sedentary", "light", 
## "moderate" and "vigorous":

## You have to compute 'm_trained_HMM' and 
## 'hidden_PA_levels' first (see above)

solution_of_HMM_based_cut_off_point_method <- 
 cut_off_point_method(x = x, 
  hidden_PA_levels = hidden_PA_levels, 
  cut_points = c(5,15,23), 
  names_activity_ranges = c("SED","LIG","MOD","VIG"), 
  bout_lengths = c(1,1,2,2,3,3,4,4,5,5,6,12,13,40,41,265,1,265), 
  plotting=1)

---

The Generalized Poisson Distribution

Description

Density, distribution function and random generation function for the generalized Poisson distribution.

Usage

dgenpois(x, lambda1, lambda2)
pgenpois(q, lambda1, lambda2)
rgenpois(n, lambda1, lambda2)

Arguments

x	a vector object of (non-negative integer) quantiles.
q	a numeric value.
n	number of random values to return.
lambda1	a single numeric value for parameter lambda1 with $lambda1 > 0$.
lambda2	a single numeric value for parameter lambda2 with $0 \le lamdba2 < 1$. When lambda2=0, the generalized Poisson distribution reduces to the Poisson distribution.

Details

The generalized Poisson distribution has the density

$p(x) = \lambda_1 (\lambda_1 + \lambda_2 \cdot x)^{x-1} \frac{ \exp(-\lambda_1-\lambda_2 \cdot x) )}{x!}$

for $x = 0,1,2,\ldots$,b

with $\mbox{E}(X)=\frac{\lambda_1}{1-\lambda_2}$ and variance $\mbox{var}(X)=\frac{\lambda_1}{(1-\lambda_2)^3}$.

Value

dgenpois gives the density, pgenpois gives the distribution function and rgenpois generates random deviates.

Author(s)

Based on Joe and Zhu (2005). Implementation by Vitali Witowski (2013).

References

Joe, H., Zhu, R. (2005). Generalized poisson distribution: the property of mixture of poisson and comparison with negative binomial distribution. Biometrical Journal 47(2):219–229.

See Also

Distributions for other standard distributions, including dpois for the Poisson distribution.

Examples

dgenpois(x = seq(0,20), lambda1 = 10, lambda2 = 0.5) 

pgenpois(q = 5, lambda1 = 10, lambda2 = 0.5) 

hist(rgenpois(n = 1000, lambda1 = 10, lambda2 = 0.5) )

---

Estimation by Directly Maximizing the log-Likelihood

Description

Estimates the parameters of a (stationary) discrete-time hidden Markov model by directly maximizing the log-likelihood of the model using the nlm-function. See MacDonald & Zucchini (2009, Paragraph 3) for further details.

Usage

direct_numerical_maximization(x, m, delta, gamma, 
     distribution_class, distribution_theta, 
     DNM_limit_accuracy = 0.001, DNM_max_iter = 50, 
     DNM_print = 2)

Arguments

x	a vector object containing the time-series of observations that are assumed to be realizations of the (hidden Markov state dependent) observation process of the model.
m	a (finite) number of states in the hidden Markov chain.
delta	a vector object containing starting values for the marginal probability distribution of the m states of the Markov chain at the time point t=1. This implementation of the algorithm uses the stationary distribution as delta.
gamma	a matrix (nrow=ncol=m) containing starting values for the transition matrix of the hidden Markov chain.
distribution_class	a single character string object with the abbreviated name of the m observation distributions of the Markov dependent observation process. The following distributions are supported by this algorithm: Poisson (pois); generalized Poisson (genpois); normal (norm, discrete log-Likelihood not applicable by this algorithm).
distribution_theta	a list object containing starting values for the parameters of the m observation distributions that are dependent on the hidden Markov state.
DNM_limit_accuracy	a single numerical value representing the convergence criterion of the direct numerical maximization algorithm using the nlm-function. Default value is 0.001.
DNM_max_iter	a single numerical value representing the maximum number of iterations of the direct numerical maximization using the nlm-function. Default value is 50.
DNM_print	a single numerical value to determine the level of printing of the nlm-function. See nlm-function for further informations. The value 0 suppresses, that no printing will be outputted. Default value is 2 for full printing.

Value

direct_numerical_maximization returns a list containing the estimated parameters of the hidden Markov model and other components.

x	input time-series of observations.
m	input number of hidden states in the Markov chain.
logL	a numerical value representing the logarithmized likelihood calculated by the forward_backward_algorithm.
AIC	a numerical value representing Akaike's information criterion for the hidden Markov model with estimated parameters.
BIC	a numerical value representing the Bayesian information criterion for the hidden Markov model with estimated parameters.
delta	a vector object containing the estimates for the marginal probability distribution of the m states of the Markov chain at time-point point t=1.
gamma	a matrix containing the estimates for the transition matrix of the hidden Markov chain.
distribution_theta	a list object containing estimates for the parameters of the m observation distributions that are dependent on the hidden Markov state.
distribution_class	input distribution class.

Author(s)

The basic algorithm of a Poisson-HMM is provided by MacDonald & Zucchini (2009, Paragraph A.1). Extension and implementation by Vitali Witowski (2013).

References

MacDonald, I. L., Zucchini, W. (2009) Hidden Markov Models for Time Series: An Introduction Using R, Boca Raton: Chapman & Hall.

See Also

HMM_based_method, HMM_training, Baum_Welch_algorithm, forward_backward_algorithm,

initial_parameter_training

Examples

################################################################
### Fictitious observations ####################################
################################################################

x <- c(1,16,19,34,22,6,3,5,6,3,4,1,4,3,5,7,9,8,11,11,
  14,16,13,11,11,10,12,19,23,25,24,23,20,21,22,22,18,7,
  5,3,4,3,2,3,4,5,4,2,1,3,4,5,4,5,3,5,6,4,3,6,4,8,9,12,
  9,14,17,15,25,23,25,35,29,36,34,36,29,41,42,39,40,43,
  37,36,20,20,21,22,23,26,27,28,25,28,24,21,25,21,20,21,
  11,18,19,20,21,13,19,18,20,7,18,8,15,17,16,13,10,4,9,
  7,8,10,9,11,9,11,10,12,12,5,13,4,6,6,13,8,9,10,13,13,
  11,10,5,3,3,4,9,6,8,3,5,3,2,2,1,3,5,11,2,3,5,6,9,8,5,
  2,5,3,4,6,4,8,15,12,16,20,18,23,18,19,24,23,24,21,26,
  36,38,37,39,45,42,41,37,38,38,35,37,35,31,32,30,20,39,
  40,33,32,35,34,36,34,32,33,27,28,25,22,17,18,16,10,9,
  5,12,7,8,8,9,19,21,24,20,23,19,17,18,17,22,11,12,3,9,
  10,4,5,13,3,5,6,3,5,4,2,5,1,2,4,4,3,2,1) 


### Assummptions (number of states, probability vector, 
### transition matrix, and distribution parameters)
    m <-4
delta <- c(0.25,0.25,0.25,0.25)
gamma <- 0.7 * diag(m) + rep(0.3 / m)
distribution_class <- "pois"
distribution_theta <- list(lambda = c(4,9,17,25))

### Estimation of a HMM using the method of 
### direct numerical maximization

trained_HMM_with_m_hidden_states <- 
		direct_numerical_maximization(x = x, 
      m = m, 
      delta = delta, 
      gamma = gamma, 
      distribution_class = distribution_class,
      DNM_max_iter=100,
      distribution_theta = distribution_theta)

print(trained_HMM_with_m_hidden_states)

---

Calculating Forward and Backward Probabilities and Likelihood

Description

The function calculates the logarithmized forward and backward probabilities and the logarithmized likelihood for a discrete time hidden Markov model, as defined in MacDonald & Zucchini (2009, Paragraph 3.1- Paragraph 3.3 and Paragraph 4.1).

Usage

forward_backward_algorithm(x, delta, gamma, distribution_class, 
distribution_theta, discr_logL=FALSE, discr_logL_eps=0.5) 

### Default for normal distributions 
### (calculate non-discrete likelihood)

Arguments

x	a vector object containing the time-series of observations that are assumed to be realizations of the (hidden Markov state dependent) observation process of the model.
delta	a vector object containing values for the marginal probability distribution of the m states of the Markov chain at the time point t=1.
gamma	a matrix (nrow=ncol=m) containing values for the transition matrix of the hidden Markov chain.
distribution_class	a single character string object with the abbreviated name of the m observation distributions of the Markov dependent observation process. The following distributions are supported: Poisson (pois); generalized Poisson (genpois); normal (norm); geometric (geom).
distribution_theta	a list object containing the parameter values for the m observation distributions of the observation process that are dependent on the hidden Markov state.
discr_logL	a logical object. It is TRUE if the discrete log-likelihood shall be calculated (for distribution_class="norm" ) instead of the general log-likelihood. See MacDonald & Zucchini (2009, Paragraph 1.2.3) for further details. Default is FALSE.
discr_logL_eps	a single numerical value to approximately determine the discrete log-likelihood for a hidden Markov model based on normal distributions (for "norm"). The default value is 0.5. See MacDonald & Zucchini (2009, Paragraph 1.2.3) for further details.

Value

forward_backward_algorithm returns a list containing the logarithmized forward and backward probabilities and the logarithmized likelihood.

log_alpha	a (T,m)-matrix (when T indicates the length/size of the observation time-series and m the number of states of the HMM) containing the logarithmized forward probabilities.
log_beta	a (T,m)-matrix (when T indicates the length/size of the observation time-series and m the number of states of the HMM) containing the logarithmized backward probabilities.
logL	a single numerical value representing the logarithmized likelihood.
logL_calculation	a single character string object which indicates how logL has been calculated (see Zucchini (2009) Paragrahp 3.1-3.4, 4.1, A.2.2, A.2.3 for further details).

Author(s)

The basic algorithm for a Poisson-HMM is provided by MacDonald & Zucchini (2009, Paragraph A.2.2). Extension and implementation by Vitali Witowski (2013).

References

MacDonald, I. L., Zucchini, W. (2009) Hidden Markov Models for Time Series: An Introduction Using R, Boca Raton: Chapman & Hall.

See Also

HMM_based_method, HMM_training, Baum_Welch_algorithm, direct_numerical_maximization, initial_parameter_training

Examples

################################################################
### Fictitious observations ####################################
################################################################

x <- c(1,16,19,34,22,6,3,5,6,3,4,1,4,3,5,7,9,8,11,11,
  14,16,13,11,11,10,12,19,23,25,24,23,20,21,22,22,18,7,
  5,3,4,3,2,3,4,5,4,2,1,3,4,5,4,5,3,5,6,4,3,6,4,8,9,12,
  9,14,17,15,25,23,25,35,29,36,34,36,29,41,42,39,40,43,
  37,36,20,20,21,22,23,26,27,28,25,28,24,21,25,21,20,21,
  11,18,19,20,21,13,19,18,20,7,18,8,15,17,16,13,10,4,9,
  7,8,10,9,11,9,11,10,12,12,5,13,4,6,6,13,8,9,10,13,13,
  11,10,5,3,3,4,9,6,8,3,5,3,2,2,1,3,5,11,2,3,5,6,9,8,5,
  2,5,3,4,6,4,8,15,12,16,20,18,23,18,19,24,23,24,21,26,
  36,38,37,39,45,42,41,37,38,38,35,37,35,31,32,30,20,39,
  40,33,32,35,34,36,34,32,33,27,28,25,22,17,18,16,10,9,
  5,12,7,8,8,9,19,21,24,20,23,19,17,18,17,22,11,12,3,9,
  10,4,5,13,3,5,6,3,5,4,2,5,1,2,4,4,3,2,1) 


### Assummptions (number of states, probability vector, 
### transition matrix, and distribution parameters)

m <-4

delta <- c(0.25,0.25,0.25,0.25)

gamma <- 0.7 * diag(m) + rep(0.3 / m)

distribution_class <- "pois"

distribution_theta <- list(lambda = c(4,9,17,25))


### Calculating logarithmized forward/backward probabilities 
### and logarithmized likelihood

forward_and_backward_probabilities_and_logL <- 
    forward_backward_algorithm (x = x, 
        delta = delta, 
        gamma = gamma, 
        distribution_class = distribution_class, 
        distribution_theta = distribution_theta)

print(forward_and_backward_probabilities_and_logL)

---

Hidden Markov Method for Predicting Physical Activity Patterns

Description

This function assigns a physical activity range to each observation of a time-series (such as a sequence of impulse counts recorded by an accelerometer) using hidden Markov models (HMM). The activity ranges are defined by thresholds called cut-off points. Basically, this function combines HMM_training, HMM_decoding and cut_off_point_method. See Details for further information.

Usage

HMM_based_method(x, cut_points, distribution_class, 
                 min_m = 2, max_m = 6, n = 100,
                 max_scaled_x = NA, names_activity_ranges = NA,  
                 discr_logL = FALSE, discr_logL_eps = 0.5, 
                 dynamical_selection = TRUE, training_method = "EM", 
                 Mstep_numerical = FALSE, BW_max_iter = 50, 
                 BW_limit_accuracy = 0.001, BW_print = TRUE,
                 DNM_max_iter = 50, DNM_limit_accuracy = 0.001, 
                 DNM_print = 2, decoding_method = 'global',
                 bout_lengths = NULL, plotting = 0)

Arguments

x	a vector object of length T containing non-negative observations of a time-series, such as a sequence of accelerometer impulse counts, which are assumed to be realizations of the (hidden Markov state dependent) observation process of a HMM.
cut_points	a vector object containing cut-off points to separate activity ranges. For instance, the vector c(7,15,23) separates the four activity ranges [0,7), [7,15), [15,23) and [23,Inf).
distribution_class	a single character string object with the abbreviated name of the m observation distributions of the Markov dependent observation process. The following distributions are supported: Poisson (pois); generalized Poisson (genpois); normal (norm)).
min_m	miminum number of hidden states in the hidden Markov chain. Default value is 2.
max_m	maximum number of hidden states in the hidden Markov chain. Default value is 6.
n	a single numerical value specifying the number of samples. Default value is 100.
max_scaled_x	an optional numerical value, to be used to scale the observations of the time-series x before the hidden Markov model is trained and decoded (see Details). Default value is NA.
names_activity_ranges	an optional character string vector to name the activity ranges induced by the cut-points. This vector must contain one element more than the vector cut_points.
discr_logL	a logical object indicating whether the discrete log-likelihood should be used (for "norm") for estimating the model specific parameters instead of the general log-likelihood. See MacDonald & Zucchini (2009, Paragraph 1.2.3) for further details. Default is FALSE.
discr_logL_eps	a single numerical value to approximate the discrete log-likelihood for a hidden Markov model based on nomal distributions (for distribution_class="norm"). The default value is 0.5.
dynamical_selection	a logical value indicating whether the method of dynamical initial parameter selection should be applied (see HMM_training for details). Default is TRUE.
training_method	a logical value indicating whether the Baum-Welch algorithm ("EM") or the method of direct numerical maximization ("numerical") should be applied for estimating the model specific parameters of the HMM. See Baum_Welch_algorithm and direct_numerical_maximization for further details. Default is training_method="EM".
Mstep_numerical	a logical object indicating whether the Maximization Step of the Baum-Welch algorithm shall be performed by numerical maximization. Default is FALSE.
BW_max_iter	a single numerical value representing the maximum number of iterations in the Baum-Welch algorithm. Default value is 50.
BW_limit_accuracy	a single numerical value representing the convergence criterion of the Baum-Welch algorithm. Default value is 0.001.
BW_print	a logical object indicating whether the log-likelihood at each iteration-step shall be printed. Default is TRUE.
DNM_max_iter	a single numerical value representing the maximum number of iterations of the numerical maximization using the nlm-function (used to perform the M-step of the Baum-Welch-algorithm). Default value is 50.
DNM_limit_accuracy	a single numerical value representing the convergence criterion of the numerical maximization algorithm using the nlm function (used to perform the M-step of the Baum-Welch-algorithm). Default value is 0.001.
DNM_print	a single numerical value to determine the level of printing of the nlm-function. See nlm-function for further informations. The value 0 suppresses, that no printing will be outputted. Default value is 2 for full printing.
decoding_method	a string object to choose the applied decoding-method to decode the HMM given the time-series of observations x. Possible values are "global" (for the use of the Viterbi_algorithm) and "local" (for the use of the local_decoding_algorithm). Default value is "global".
bout_lengths	a vector object (with even number of elemets) to define the range of the bout intervals (see Details for the definition of bouts). For instance, bout_lengths=c(1,1,2,2,3,10,11,20,1,20) defines the five bout intervals [1,1] (1 count); [2,2] (2 counts); [3,10] (3-10 counts); [11,20] (11-20 counts); [1,20] (1-20 counts - overlapping with other bout intervalls is possible). Default value is bout_lengths=NULL.
plotting	a numeric value between 0 and 5 (generates different outputs). NA suppresses graphical output. Default value is 0. 0: output 1-5 1: summary of all results 2: time series of activity counts, classified into activity ranges 3: time series of bouts (and, if available, the sequence of the estimated hidden physical activity levels, extracted by decoding a trained HMM, in green colour) 4: barplots of absolute and relative frequencies of time spent in different activity ranges 5: barplots of relative frequencies of the lenghts of bout intervals (overall and by activity ranges )

Details

The function combines HMM_training, HMM_decoding and cut_off_point_method as follows:

Step 1: HMM_training trains the most likely HMM for a given time-series of accelerometer counts.
Step 2: HMM_decoding decodes the trained HMM (Step 1) into the most likely sequence of hidden states corresponding to the given time-series of observations (respectively the most likely sequence of physical activity levels corresponding to the time-series of accelerometer counts).
Step 3. cut_off_point_method assigns an activity range to each accelerometer count by its hidden physical activity level (extracted in Step 2).

Value

HMM_based_method returns a list containing the output of the trained hidden Markov model, including the selected number of states m (i.e., number of physical activities) and plots key figures.

trained_HMM_with_selected_m	a list object containing the trained hidden Markov model including the selected number of states m (see HMM_training for further details).
decoding	a list object containing the output of the decoding (see HMM_decoding for further details)

.

extendend_cut_off_point_method

a list object containing the output of the cut-off point method. The counts x are classified into the activity ranges by the corresponding sequence of hidden PA-levels, which were decoded by the HMM (see cut_off_point_method for further details).

Note

The parameter max_scaled_x can be applied to scale the values of the observations. This might prevent the alogrithm from numerical instabilities. At the end, the results are internaly rescaled to the original scale. For instance, a value of max_scaled_x=200 shrinks the count values of the complete time-series x to a maximum of 200. Training and decoding of the HMM is carried out using the scaled time-series.
From our experience, especially time-series with observations values >1500, or where T > 1000, show numerical instabilities. We then advice to make use of max_scaled_x .

The extention of the cut-off point method using a Poisson based HMM has been provided and evaluated successfully on simulated data firstly by Barbara Brachmann in her diploma thesis (see References).

Author(s)

Vitali Witowski (2013).

References

Brachmann, B. (2011). Hidden-Markov-Modelle fuer Akzelerometerdaten. Diploma Thesis, University Bremen - Bremen Institute for Prevention Research and Social Medicine (BIPS).

MacDonald, I. L., Zucchini, W. (2009) Hidden Markov Models for Time Series: An Introduction Using R, Boca Raton: Chapman & Hall.

Witowski, V., Foraita, R., Pitsiladis, Y., Pigeot, I., Wirsik, N. (2014) Using hidden Markov models to improve quantifying physical activity in accelerometer data - A simulation study. PLOS ONE. 9(12), e114089. http://dx.doi.org/10.1371/journal.pone.0114089

See Also

initial_parameter_training, Baum_Welch_algorithm, direct_numerical_maximization, AIC_HMM, BIC_HMM, HMM_training, Viterbi_algorithm, local_decoding_algorithm, cut_off_point_method

Examples

################################################################
### Fictitious activity counts #################################
################################################################

x <- 100 * c(1,16,19,34,22,6,3,5,6,3,4,1,4,3,5,7,9,8,11,11,
  14,16,13,11,11,10,12,19,23,25,24,23,20,21,22,22,18,7,
  5,3,4,3,2,3,4,5,4,2,1,3,4,5,4,5,3,5,6,4,3,6,4,8,9,12,
  9,14,17,15,25,23,25,35,29,36,34,36,29,41,42,39,40,43,
  37,36,20,20,21,22,23,26,27,28,25,28,24,21,25,21,20,21,
  11,18,19,20,21,13,19,18,20,7,18,8,15,17,16,13,10,4,9,
  7,8,10,9,11,9,11,10,12,12,5,13,4,6,6,13,8,9,10,13,13,
  11,10,5,3,3,4,9,6,8,3,5,3,2,2,1,3,5,11,2,3,5,6,9,8,5,
  2,5,3,4,6,4,8,15,12,16,20,18,23,18,19,24,23,24,21,26,
  36,38,37,39,45,42,41,37,38,38,35,37,35,31,32,30,20,39,
  40,33,32,35,34,36,34,32,33,27,28,25,22,17,18,16,10,9,
  5,12,7,8,8,9,19,21,24,20,23,19,17,18,17,22,11,12,3,9,
  10,4,5,13,3,5,6,3,5,4,2,5,1,2,4,4,3,2,1)  

	   
### Fictitious cut-off points that produce four so-called 
### activity ranges "sedentary", "light", "moderate", 
### and "vigorous".
cut_points <- 100 * c(7,15,23)
names_activity_ranges <- c("SED","LIG","MOD","VIG")


### Plot fictitious activity counts
plot(x, main = "counts with high values", 
     xlab = "time/epoch", ylab = "counts")
abline(h = cut_points, col = "grey50", lty = "dashed")


################################################################
### Comparing the results of the traditional ################### 
### cut-off point method and the new HMM-based method ##########
################################################################

### Apply the traditional cut-off point method to assign 
### physical activity ranges to each observed count

solution_of_tradtionional_cut_off_point_method <-
   cut_off_point_method(x = x, 
       hidden_PA_levels = NA, 
       cut_points = cut_points, 
       names_activity_ranges = names_activity_ranges, 
       bout_lengths = c(1,1,2,2,3,3,4,4,5,5,6,12, 
       13,40,41,265,1,265), 
	     plotting = 1)

### Apply the HMM-based method to assign physical activity 
### ranges to the hidden physical activity level of each count

solution_of_HMM_based_method <- 
    HMM_based_method(x = x, 
      max_scaled_x = 50, 
      cut_points  =cut_points, 
    	min_m = 2, 
    	max_m = 6, 
    	names_activity_ranges = names_activity_ranges, 
      distribution_class = "pois", 
      training_method = "EM", 
      decoding_method = "global", 
      bout_lengths = c(1,1,2,2,3,3,4,4,5,5,6,12,
      13,40,41,265,1,265), 
      plotting = 1)

		
### Print details of the traditional cut-off point method 
### and the new HMM-based method
print(solution_of_tradtionional_cut_off_point_method)
print(solution_of_HMM_based_method)

---

Algorithm for Decoding Hidden Markov Models (local or global)

Description

The function decodes a hidden Markov model into a most likely sequence of hidden states. Furthermore this function provides estimated observation values along the most likely sequence of hidden states. See Details for more information.

Usage

HMM_decoding(x, m, delta, gamma, distribution_class, 
		     distribution_theta, decoding_method = "global", 
		     discr_logL = FALSE, discr_logL_eps = 0.5)

Arguments

x	a vector object containing the time-series of observations that are assumed to be realizations of the (hidden Markov state dependent) observation process of the model.
m	a (finite) number of states in the hidden Markov chain.
delta	a vector object containing values for the marginal probability distribution of the m states of the Markov chain at the time point t=1.
gamma	a matrix (ncol=nrow=m) containing values for the transition matrix of the hidden Markov chain.
distribution_class	a single character string object with the abbreviated name of the m observation distributions of the Markov dependent observation process. The following distributions are supported by this algorithm: Poisson (pois); generalized Poisson (genpois); normal (norm); geometric (geom).
distribution_theta	a list object containing the parameter values for the m observation distributions that are dependent on the hidden Markov state.
decoding_method	a string object to choose the applied decoding-method to decode the HMM given the time-series of observations x. Possible values are "global" (for the use of the Viterbi_algorithm) and "local" (for the use of the local_decoding_algorithm). Default value is "global".
discr_logL	a logical object. It is TRUE if the discrete log-likelihood shall be calculated (for distribution_class="norm" instead of the general log-likelihood). Default is FALSE.
discr_logL_eps	a single numerical value to approximately determine the discrete log-likelihood for a hidden Markov model based on nomal distributions (for "norm"). The default value is 0.5.

Details

More precisely, the function works as follows:

Step 1: In a first step, the algorithm decodes a HMM into the most likely sequence of hidden states, given a time-series of observations. The user can choose between a global and a local approch.
If decoding_method="global" is applied, the function calls Viterbi_algorithm to determine the sequence of most likely hidden states for all time points simultaneously.
If decoding_method="local" is applied, the function calls local_decoding_algorithm to determine the most likely hidden state for each time point seperately.

Step 2: In a second step, this function links each observation to the mean of the distribution, that corresponds to the decoded state at this point in time.

Value

HMM_decoding returns a list containing the following two components:

decoding_method	a string object indicating the applied decoding method.
decoding	a numerical vector containing the most likely sequence of hidden states as decoded by the Viterbi_algorithm (if "global" was applied) or by the local_decoding_algorithm (if "local" was applied).
decoding_distr_means	a numerical vector of estimated oberservation values along the most likely seuquence of hidden states (see decoding and Step 2).

Author(s)

Vitali Witowski (2013).

References

MacDonald, I. L., Zucchini, W. (2009) Hidden Markov Models for Time Series: An Introduction Using R, Boca Raton: Chapman & Hall.

See Also

local_decoding_algorithm, Viterbi_algorithm

Examples

################################################################
### i) HMM-training  ###########################################
################################################################

################################################################
### Fictitious observations ####################################
################################################################

x <- c(1,16,19,34,22,6,3,5,6,3,4,1,4,3,5,7,9,8,11,11,
  14,16,13,11,11,10,12,19,23,25,24,23,20,21,22,22,18,7,
  5,3,4,3,2,3,4,5,4,2,1,3,4,5,4,5,3,5,6,4,3,6,4,8,9,12,
  9,14,17,15,25,23,25,35,29,36,34,36,29,41,42,39,40,43,
  37,36,20,20,21,22,23,26,27,28,25,28,24,21,25,21,20,21,
  11,18,19,20,21,13,19,18,20,7,18,8,15,17,16,13,10,4,9,
  7,8,10,9,11,9,11,10,12,12,5,13,4,6,6,13,8,9,10,13,13,
  11,10,5,3,3,4,9,6,8,3,5,3,2,2,1,3,5,11,2,3,5,6,9,8,5,
  2,5,3,4,6,4,8,15,12,16,20,18,23,18,19,24,23,24,21,26,
  36,38,37,39,45,42,41,37,38,38,35,37,35,31,32,30,20,39,
  40,33,32,35,34,36,34,32,33,27,28,25,22,17,18,16,10,9,
  5,12,7,8,8,9,19,21,24,20,23,19,17,18,17,22,11,12,3,9,
  10,4,5,13,3,5,6,3,5,4,2,5,1,2,4,4,3,2,1)  

### Train hidden Markov model for different number of 
### states m=2,...,6 and select the optimal model

m_trained_HMM <- 
    HMM_training(x = x, 
       min_m = 2, 
       max_m = 6, 
       distribution_class="pois")$trained_HMM_with_selected_m
         
                 
################################################################
################################################################
### ii) Global decoding ########################################
################################################################
################################################################

### Decode the trained HMM using the Viterbi algorithm to get 
### the estimated sequence of hidden physical activity levels
global_decoding <- 
    HMM_decoding(x = x, 
        m = m_trained_HMM$m, 
        delta = m_trained_HMM$delta, 
        gamma = m_trained_HMM$gamma, 
        distribution_class = m_trained_HMM$distribution_class, 
        distribution_theta = m_trained_HMM$distribution_theta,
        decoding_method = "global")
          
### Globally most likely sequence of hidden states, 
### i.e. in this case sequence of activity levels
global_decoding$decoding

par(mfrow = c(1,1))
plot(global_decoding$decoding)


### Plot the observed impulse counts and the most likely 
### sequence (green) according to the Viterbi algorithm that 
### generated these observations
plot(x)
lines(global_decoding$decoding_distr_means, col = "green")

################################################################
################################################################ 
### iii) Local decoding ########################################
################################################################ 
################################################################

### Decode the trained HMM using the local decoding algorithm 
### to get the estimated sequence of hidden physical activity 
### levels
local_decoding <- 
   HMM_decoding(x = x, 
   m = m_trained_HMM$m, 
   delta = m_trained_HMM$delta, 
   gamma = m_trained_HMM$gamma, 
   distribution_class = m_trained_HMM$distribution_class, 
   distribution_theta = m_trained_HMM$distribution_theta,
   decoding_method = "local")
        

### Locally most likely sequence of hidden states, 
### i.e. in this case sequence of activity levels
local_decoding$decoding

par(mfrow=c(1,1))
plot(local_decoding$decoding)


### Plot the observed impulse counts and the most likely 
### sequence (green) according to the local decoding algorithm 
### that generated these observations
plot(x)
lines(local_decoding$decoding_distr_means, col = "red")

################################################################
################################################################
### iv) Comparison of global and local decoding ################
################################################################
################################################################

### Comparison of global decoding (green), local decoding (red) 
### and the connection to the closest mean (blue)
print(global_decoding$decoding)  
print(local_decoding$decoding)

### Plot comparison 
par(mfrow = c(2,2))
plot(global_decoding$decoding[seq(230,260)], col = "green", 
  ylab = "global decoding", main = "(zooming)")
  plot(x[seq(230,260)], ylab = "global decoding", 
  main = "(zooming x[seq(230,260)])")
lines(global_decoding$decoding_distr_means[seq(230,260)], 
  col = "green")
plot(local_decoding$decoding[seq(230,260)], col = "red", 
  ylab = "local decoding", main = "(zooming)")
plot(x[seq(230,260)], ylab = "local decoding", 
  main = "(zooming x[seq(230,260)])")
lines(local_decoding$decoding_distr_means[seq(230,260)], 
  col = "red")
par(mfrow = c(1,1))

---