Novel approach to modeling high-frequency activity data to assess therapeutic effects of analgesics in chronic pain conditions

Osteoarthritis (OA) is a chronic condition often associated with pain, affecting approximately fourteen percent of the population, and increasing in prevalence. A globally aging population have made treating OA-associated pain as well as maintaining mobility and activity a public health priority. OA affects all mammals, and the use of spontaneous animal models is one promising approach for improving translational pain research and the development of effective treatment strategies. Accelerometers are a common tool for collecting high-frequency activity data on animals to study the effects of treatment on pain related activity patterns. There has recently been increasing interest in their use to understand treatment effects in human pain conditions. However, activity patterns vary widely across subjects; furthermore, the effects of treatment may manifest in higher or lower activity counts or in subtler ways like changes in the frequency of certain types of activities. We use a zero inflated Poisson hidden semi-Markov model to characterize activity patterns and subsequently derive estimators of the treatment effect in terms of changes in activity levels or frequency of activity type. We demonstrate the application of our model, and its advance over traditional analysis methods, using data from a naturally occurring feline OA-associated pain model.


Theoretical properties
Identifiability of the zero-inflated hidden Markov model Proof. According to Proposition 3.1 in Gassiat et al (2013), if the state-dependent probability measures are linearly independent and the weight matrix has full rank, then the parameters in the hidden Markov model are identifiable up to label switching of the hidden states as soon as the transition matrix has full rank. In our case, define the M state-dependent probability measures µ 1 , ..., µ M by where φ 1 is a point mass at zero, φ l (l=2,...,M+1) are the M different Poisson distributions, and ψ = [ψ jl ] is the weight matrix. Denote the structural zero proportion for the zero-inflated Poisson distribution in state 1 by p. Then the M × (M + 1) weight matrix ψ is where each row represents a latent state; column 1 represents the probability measure of the structural zeros; column 2 -column M+1 represent the M distinct Poisson distributions. It can be seen that ψ has rank M.

Proof of Lemma
Define the conditional log likelihood function for the subject-specific hidden semi-Markov model given initial state S 0 i as for subject i = 1, ..., n. Under (A1), this is equivalent to the conditional log likelihood for a correspondingly reparameterized hidden Markov model with an expanded state space. Under (A0) -(A4), the uniform convergence of the normalized log likelihood is proved in Proposition 2 in Douc et al (2004), i.e., sup where l i (θ i ) is a deterministic asymptotic function such that To show that l i (θ i ) is maximized at θ i = θ i for i = 1, ..., n, where the inner conditional expectation is the Kullback-Leibler divergence and is maximized at θ i = θ i . By the identifiability assumption (A5), θ i is the unique maximizer of l i (θ i ). Therefore, it follows that arg max with probability 1 as T → ∞.

Proof of Theorem
It suffices to proof the theorem for Ω Ω Ω m,n , which can then be naturally extended to Γ Γ Γ m,n and Λ Λ Λ m, ,n using the exact same approach, where m, = 1, ..., M. By Assumptions (A0) -(A5), we have Lemma 1, which implies that the MLEs for subject-specific treatment effects b b b 1,m,T from the first-stage fitting are strongly consistent; this implies that Recall that the linear regression model we fit in the second stage is of the form ; e m,i and is independent across all indices with mean 0 and variance σ 2 ; u m,i and furthermore is independent across all indices with mean 0 and variance Var[vec( b 1,m,T ) i ]. In the corresponding matrix notation, Consider the true linear model for We also assume that lim )) converges in distribution to a multivariate Gaussian distribution with mean 0 and covariance By the consistency of the sandwich variance estimator, Thus, by Slutsky's Theorem, the following quantity converges in distribution to a multivariate Gaussian distribution with mean zero and identity covariance,

Simulation experiments
We study the finite sample performance of the proposed two-stage estimator using a suite of simulation experiments. Across these experiments, we consider three different pathways for the treatment effect: (Scenario I) a change in the state-dependent means and the proportion of zeros; (Scenario II) a change in the latent state durations; and (Scenario III) a change in the state-transition probabilities. In each scenario, we simulate minute-by-minute activity counts during a two week period for n = 10 and n = 20 subjects. Each simulated trajectory is divided into an initial one-week baseline phase followed by a one-week treatment period; thus, each patient has 10,080 observations per period. The activity counts are generated using a three state zero-inflated Poisson hidden semi-Markov model, with the maximum latent state duration set to be 20. The initial state probabilities are distributed as follows: let U 1 , U 2 ∼ i.i.d. Uniform(0, 1) and set δ 1 = min(U 1 , U 2 ), δ 2 = max(U 1 , U 2 ) − δ 1 , and δ 3 = 1 − δ 1 − δ 2 . We generate subject covariates so that half of the subjects are male, and set the baseline Poisson means for the male subjects to be 5% less than the female subjects while the baseline zero proportion is 5% larger than the female subjects. Furthermore, for all subjects (regardless of sex) the Poisson means are 5% less at night than during the day and the zero proportions are 5% larger at night than during the day. These proportions are consistent with the accelerometer data collected in the study of the treatment effect of meloxicam. In Scenario I, treatment increases the log mean activity in all the latent states by 10% and decreases the log odds of zero in state 1 by 10%, so that (300), 0.1 2 , and b 0,3,i iid ∼ N log(700), 0.1 2 . Let N [a,b] (µ, σ 2 ) denote a normal distribution with mean µ and variance σ 2 truncated to [a, b]. In this scenario, we assume that both the transition probabilities and the latent state durations are unaffected by the treatment and that Q i (1, 2) iid ∼ N [0,1] (0.5, 0.01), This setup mimics the real accelerometer data where state 1 (lowest activity level) dominates the other latent states. We do not allow transition back to the same state as the latent state durations modeled separatetly; in Scenario 1, the latent state durations are uniformly distributed on {1, 2, . . . , 20} for each state.
In Scenario II, we assume that the treatment increases the log mean activity in state 3 by 10% but decreases the log mean 4/6 activity in state 1 and 2 by 10%. Further, the treatment increases the log odds of zero in state 1 by 10%, such that (700), 0.01). In this scenario, treatment will increase activity in the highest level state but make the subject less agitated in lower level activity states. In this scenario, both the transition probabilities and the latent state durations are unaffected by the treatment and match the settings of Scenario I.
In Scenario III, we assume that treatment accelerates the duration time by 10% in each latent state, so that subjects are more likely to actively switch between states when they are on treatment, i.e., for m = 1, 2, 3, where f 1,i , f 2,i , f 3,i are the exponential densities with rate parameters ξ 1,i , ξ 2,i , ξ 3,i , such that ξ 1,i iid ∼ N(5, 0.1 2 ), ξ 2,i iid ∼ N(3, 0.1 2 ), and ξ 3,i iid ∼ N(2, 0.1 2 ). This setup also mimics the study of the treatment effect of the meloxicam in that the expected duration in state 1 is longer than in other states. We assume that the treatment does not affect the state-dependent Poisson means and zero proportion so that (700), 0.01 . The transition probabilities are the same as in Scenario I.
In Scenario IV, we assume that treatment affects the state transition probabilities by increasing the odds of moving to the next most active state by 10% so that We compare the proposed estimator with the baseline treatment effect estimator constructed using a generalized estimating equation (GEE). This is equivalent to the repeated-measures ANOVA model as the data are complete and balanced. In the GEE model, denote the new response variable in the model byỸ τ i , which is the average activity count for subject i on day τ = 1, . . . , 14.
; β β β = [β 0 , β 1 , β 2 ] so that β 1 represents the treatment effect; the working covariance V i is assumed to be compound symmetric. Table 1 shows the mean and standard error for the estimated treatment effect estimators in the proposed method using 500 Monte Carlo replications. The estimators are approximately unbiased in all four scenarios and the standard error is small even with a moderate number of samples owing to within-subject pooling across observations. Table 2 shows the mean and standard error for the estimated treatment effect in the baseline GEE model. The GEE method manages to detect the treatment effect in Scenarios I and IV, where there is an increase in the marginal mean activity counts under treatment; however, it fails to detect the nonzero conditional treatment effect in different latent states using an α = 0.05 level test in Scenario II and III because there is only a small difference in the marginal mean activity counts.