Hidden Markov Model Based on Logistic Regression

Abstract

A hidden Markov model (HMM) is a useful tool for modeling dependent heterogeneous phenomena. It can be used to find factors that affect real-world events, even when those factors cannot be directly observed. HMMs differ from traditional methods by using state variables and mixture distributions to model the hidden states. This allows HMMs to find relationships between variables even when the variables cannot be directly observed. HMM can be extended, allowing the transition probabilities to depend on covariates. This makes HMMs more flexible and powerful, as they can be used to model a wider range of sequential data. Modeling covariates in a hidden Markov model is particularly difficult when the dimension of the state variable is large. To avoid these difficulties, Markovian properties are achieved by implanting the previous state variables to the logistic regression model. We apply the proposed method to find the factors that affect the hidden state of matsutake mushroom growth, in which it is hard to find covariates that directly affect matsutake mushroom growth in Korea. We believe that this method can be used to identify factors that are difficult to find using traditional methods.

1. Introduction

The hidden Markov model (HMM), introduced by Baum and his colleagues in a series of papers (Baum and Petrie [1]), is a useful way to represent dependent heterogeneous phenomena. It is often used to model stochastic processes and time-dependent sequences and has been applied to a wide range of applications such as medical diagnosis, financial forecasting, and natural language processing. In many applications, the marginal distribution of observations is clearly a multi-model, meaning that the observations come from a mixture of different distributions associated with different regimes. This behavior is a key attribute of HMMs. HMMs can be used to detect these regimes and find the hidden states that correspond to them. Once the hidden states have been found, a mixture of different distributions is identified on them (Titterington et al. [2]; Zucchini and MacDonald [3]). Inference in HMMs is typically based on maximum likelihood or Bayesian approaches. However, the dependence structure of HMMs can be more difficult to compute than that of regular mixtures. Robert et al. [4] provided an efficient Bayesian estimate of HMMs via Gibbs sampling, and Chib [5] introduced a widely used state simulation procedure (see also Campbell et al. [6]). Li et al. [7] developed an approach to multivariate time series anomaly detection in which an HMM is engaged here to detect anomalies in multivariate time series. Nguyen [8] used an HMM to predict a daily stock price of three active trading stocks, Apple, Google, and Facebook, based on their historical data.

HMMs also have a number of disadvantages, such as assumption of independence and limitation of state space, while they have a number of advantages, including flexibility, interpretability, and efficiency. Another limitation of HMMs lies in directly incorporating additional explanatory variables. Covariate-dependent hidden Markov models (CDHMMs) are a class of statistical models used to model sequential data where the transition probabilities between hidden states depend on observed covariates. This makes CDHMMs well-suited for a wide range of applications, such as medical diagnosis, financial forecasting, and natural language processing. CDHMMs can be used to identify different market states, such as bull markets, bear markets, and sideways markets, or used to model the progression of Alzheimer’s disease and to identify different stages of Alzheimer’s disease. CHMMs can be also used to segment customers into different groups based on their purchase behavior. For example, CHMMs can be used to identify customers who are likely to churn and to identify customers who are likely to be high-value customers.

One of the most common approaches to modeling covariates in HMMs is to use a mixture transition distribution model. In an MTD model, the transition probabilities between hidden states are modeled as a mixture of distributions, where the weights of the mixture components are determined by the covariates. This approach is relatively simple to implement and can be used to model a wide range of covariate effects. Another approach to modeling covariates in HMMs is to use a covariate-dependent transition matrix model. In this approach, the transition probabilities between hidden states are parameterized directly using the covariates. This approach is more flexible than the MTD approach, but it can also be more difficult to implement and estimate.

Rabiner and Juang [9] applied the CHMM to speech recognition and showed that it outperformed HMMs without covariate-dependent transition probabilities. Marshall and Jones [10] discussed the application of a multi-state model to diabetic retinopathy under the assumption that a continuous time Markov process determines the transition times between disease stages. Altman [11] presented a class of mixed HMMS where covariates and random effects captured differences between an observed process and an underlying hidden process. Chamroukhi et al. [12] proposed a finite mixture model of hidden Markov regression with covariate dependence. They applied the CHMM to financial time series analysis and showed that it outperformed other methods for financial time series modeling. Maruotti [13] provided a general review for statistical methods that combine HMMs and random effects models in a longitudinal setting.

Rubin et al. [14] proposed a joint logistic regression and Markov chain model to describe a binary cross-sectional response, where the unobserved transition rates of a two-state continuous-time Markov chain are included as covariates. Sarkar and Zhu [15] proposed a novel approach for detecting clusters and regimes in time series data in the presence of random covariates, which is based on a CHMM with covariate-dependent transition probabilities. In a covariate-dependent transition matrix model, HMMs are extended by allowing the transition probabilities to depend on covariates through regression models. In this case, each transition probability requires a different regression model, which can be especially problematic especially when the dimension of state variable is large. To avoid this problem, we here introduce the previous hidden state as another explanatory variable to HMMs through a logistic regression model. That is, Markovian properties are achieved by implanting the previous state variables to the logistic regression model. Note that the proposed model can be applied regardless of the dimension of the state variable in the hidden Markov model.

In some cases, it can be difficult to find the factors that affect a particular phenomenon. This is because the factors may be hidden or may interact in complex ways. HMMs can be used to analyze hidden states without using observed values as the dependent variable. This can be an effective way to find significant variables if they cannot be found directly. In Korea, matsutake mushrooms have not been cultivated on a large scale. Instead, they are harvested from natural forests. The effect of climatic factors on matsutake mushroom yield has not been studied in detail. Some studies have been conducted to investigate the relationship between matsutake mushroom occurrence and weather patterns, but they have not been able to identify the significant meteorological factors. The proposed model is used to identify the factors that indirectly affect matsutake mushroom production. This could help to improve the cultivation of matsutake mushrooms in Korea.

The remaining sections are organized as follows. In Section 2, a brief summary of the hidden Markov model is given. In Section 3, we introduce a hidden Markov model based on logistic regression. Then, we present the hierarchical Bayesian procedure of the proposed model in Section 4. In Section 5 and Section 6, we demonstrate the hidden Markov model based on logistic regression by way of a simulation study and a case study example. The last section is the conclusion of this paper and mentions further study.

2. Hidden Markov Model

We here introduce a hidden Markov model. A brief definition of a hidden Markov model and the assumptions of the model are introduced. A hidden Markov model $\{y_t:t\in N\}$ is a particular kind of dependent mixture. In its most general form, the hidden Markov model is defined as follows:

$$\begin{array}{c}\hfill P\left(z_t\right|{\mathbf{Z}}^{(t-1)})=P\left(z_t\right|z_{t-1})fort=2,3,\dots ,T,\end{array}$$

(1)

$$\begin{array}{c}\hfill P\left(y_t\right|{\mathbf{Y}}^{(t-1)},{\mathbf{Z}}^{\left(t\right)})=P\left(y_t\right|z_t)fort=1,2,\dots ,T,\end{array}$$

(2)

where $z_t$ and $y_t$ represent an unobserved state and an observation at the time t, respectively. Note that ${\mathbf{Z}}^{\left(t\right)}$ denotes a vector of unobserved states from time 1 to time t, ${\mathbf{Y}}^{\left(t\right)}$ denotes a vector of observations from time 1 to time t, and T is the number of total observations. That is, ${\mathbf{Z}}^{\left(t\right)}=(z_1,z_2,\dots ,z_t)$ and ${\mathbf{Y}}^{\left(t\right)}=(y_1,y_2,\dots ,y_t)$.

Unlike typical mixture models, the HMM assumes that observed data are generated through a finite-valued, unobserved process. This unobserved process is assumed to be in one of a finite number of discrete states at each discrete time point, and given the previous state or states, it is assumed to transit stochastically in the Markov fashion. The data observed at each time point depend only on the value of the corresponding hidden state and are independent of others. The heterogeneity of data is represented by hidden Markov states. In other words, a pair $(z_t,y_t)$ with the state $z_t\in \{1,\dots ,k\}$ and its random value $y_t|z_t\sim f_{z_t}\left(y_t\right)$. Here, $y_t$ is considered conditionally independent given $z_t$. The HMM derives its name from the following two defining attributes. First, $z_t$ is distributed as a (finite-state) Markov chain. Second, $z_t$ is not observed.

3. Logistic Regression on Hidden State

A logistic regression is used to model the probability of a certain class or event existing such as pass or fail or to model a binary dependent variable. Consider a model with p predictors, $x_1,\dots ,x_p$, and one binary response variable Y with a probability $\pi =P(Y=1)$. A linear relationship between the predictor variables and the log-odds of the event that $Y=1$ is assumed. This linear relationship can be written in the following mathematical form, where $({\beta }_0,{\beta }_1,\dots ,{\beta }_p)$ are parameters of the model:

$$\begin{array}{c}\hfill log\left(\frac{\pi }{1-\pi }\right)={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_p{x}_p\end{array}$$

By simple algebraic manipulation, the probability that $Y=1$ is

$$\begin{array}{c}\hfill \pi =\frac{\exp ({\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_p{x}_p)}{1+\exp ({\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_p{x}_p)}.\end{array}$$

Therefore, predictors affecting the binary dependent variable Y can be found with logistic regression analysis by estimating parameters $({\beta }_0,{\beta }_1,\dots ,{\beta }_p)$.

We assume that $(x_{t1},x_{t2},\dots ,x_{tp})$ are the predictor variables and the space of hidden state $z_t$ of HMM is $\{0,1\}$ at the time point t. That is, $z_t\in \{0,1\}$ is a binary variable. Then, we model the following logistic regression:

$$\begin{array}{c}\hfill P(z_t=1)={\pi }_t,\end{array}$$

$$log\left(\frac{{\pi }_t}{1-{\pi }_t}\right)={\beta }_0+{\beta }_1{x}_{t1}+\cdots +{\beta }_p{x}_{tp}+\gamma z_{t-1}$$

for $t=1,\dots ,T$, and then the probability that $z_t=1$ is

$$\begin{array}{c}\hfill {\pi }_t=\frac{\exp ({\beta }_0+{\beta }_1{x}_{t1}+\cdots +{\beta }_p{x}_{tp}+\gamma z_{t-1})}{1+\exp ({\beta }_0+{\beta }_1{x}_{t1}+\cdots +{\beta }_p{x}_{tp}+\gamma z_{t-1})},\end{array}$$

where $\gamma z_{t-1}$ represents the Markovian property in Equation (1). That is, the term involving $z_{t-1}$ allows for a shift in the transition probability depending on the previous hidden state. Therefore, hidden state model $f_{z_t|\mathit{\beta },\gamma ,z_{t-1}}\left(z_t\right)$ can be expressed as follows:

$$\begin{array}{c}\hfill f_{z_t|\mathit{\beta },\gamma ,z_{t-1}}\left(z_t\right)\sim Ber\left({\pi }_t\right),{\pi }_t=\frac{\exp ({\mathit{\beta }}^{\prime }{\mathbf{X}}_t+\gamma z_{t-1})}{1+\exp ({\mathit{\beta }}^{\prime }{\mathbf{X}}_t+\gamma z_{t-1})},\end{array}$$

where ${\mathbf{X}}_t={(1,x_{t1},x_{t2},\dots ,x_{tp})}^{\prime }$ and $\mathit{\beta }={({\beta }_0,{\beta }_1,\dots ,{\beta }_p)}^{\prime }$. Additionally, in the function $f_{z_1|\mathit{\beta },\gamma ,z_0}\left(z_1\right)$, $z_0$ is a starting hidden state immediately before the hidden state $z_1$, and it is assumed as follows:

$$\begin{array}{c}\hfill P(z_0=1)={\pi }_0,\end{array}$$

$$\begin{array}{c}\hfill f_{z_0|{\pi }_0}\left(z_0\right)\sim Ber\left({\pi }_0\right),\end{array}$$

where ${\pi }_0$ is a parameter of the state $z_0$ distribution. Note that $Ber(·)$ denotes Bernoulli distribution.

4. Hierarchical Bayesian Approach to HMM Based on Logistic Regression

Now we introdue explanatory variables to HMMs through a logistic regression model in which we add the previous state as another explanatory variable in the logistic regression model to maintain the temporal dependence of the HMM. At first, hierarchical structures are described, then we check the full conditional distributions of each parameter for Gibbs sampling. In addition, prediction steps of the proposed model are shown.

4.1. HMM Based on Logistic Regression

In the transition probability distribution $f_{y_t|z_t}\left(y_t\right)$, it is assumed that

$$\begin{array}{c}\hfill y_t|z_t=1\sim f_1\left(y_t\right|{\theta }_1),\\ \hfill y_t|z_t=0\sim f_0\left(y_t\right|{\theta }_0)\end{array}$$

for $t=1,\dots ,T$, where ${\theta }_1$ and ${\theta }_0$ are parameters of the mixture distributions. Therefore, the data model $f_{{\mathbf{Y}}^{\left(T\right)}|{\mathbf{Z}}^{\left(T\right)},\mathit{\theta }}\left({\mathbf{Y}}^{\left(T\right)}\right)$ can be expressed as

$$\begin{array}{c}\hfill f_{{\mathbf{Y}}^{\left(T\right)}|{\mathbf{Z}}^{\left(T\right)},\mathit{\theta }}\left({\mathbf{Y}}^{\left(T\right)}\right)=\prod _{t=1}^T{f}_{y_t|z_t,\mathit{\theta }}\left(y_t\right),\end{array}$$

where $\mathit{\theta }=({\theta }_0,{\theta }_1)$. The process model $f_{{\mathbf{Z}}^{\left(T\right)},z_0,{\pi }_0|\mathit{\beta },\gamma }\left({\mathbf{Z}}^{\left(T\right)}\right)$ is given by

$$\begin{array}{c}\hfill f_{{\mathbf{Z}}^{\left(T\right)}|\mathit{\beta },\gamma ,z_0}\left({\mathbf{Z}}^{\left(T\right)}\right)=\prod _{t=1}^T{f}_{z_t|\mathit{\beta },\gamma ,z_{t-1}}\left(z_t\right).\end{array}$$

Finally, we consider prior distributions for all model parameters $\mathit{\beta },\gamma ,\mathit{\theta }$ and an initial state $z_0$. In summary, the hierarchical structure can be expressed as follows:

• Data model: $f_{{\mathbf{Y}}^{\left(T\right)}|{\mathbf{Z}}^{\left(T\right)},\mathit{\theta }}\left({\mathbf{Y}}^{\left(T\right)}\right)={\prod }_{t=1}^T{f}_{y_t|z_t,\mathit{\theta }}\left(y_t\right)$
• Process model: $f_{{\mathbf{Z}}^{\left(T\right)}|\mathit{\beta },\gamma ,z_0}\left({\mathbf{Z}}^{\left(T\right)}\right)={\prod }_{t=1}^T{f}_{z_t|\mathit{\beta },\gamma ,z_{t-1}}\left(z_t\right)$
• Prior model: $\pi (\mathit{\beta },\gamma )\pi \left(\mathit{\theta }\right)\pi (z_0,{\pi }_0)=\pi (\mathit{\beta },\gamma )\pi \left(\mathit{\theta }\right)\pi \left(z_0\right|{\pi }_0)\pi \left({\pi }_0\right)$

4.2. Bayesian Analysis

Since a direct sampling from a joint posterior distribution of model parameters is computationally difficult, simplified or approximate MCMC approaches such as the Metropolis–Hasting (M-H) approach within the Gibbs algorithm are needed.

The joint posterior distribution $\pi (\mathit{\beta },\gamma ,\mathit{\theta },{\pi }_0,z_0,{\mathbf{Z}}^{\left(T\right)}|{\mathbf{Y}}^{\left(T\right)})$ is proportional to

$$\begin{array}{ccc}& & f_{{\mathbf{Y}}^{\left(T\right)}|{\mathbf{Z}}^{\left(T\right)},\mathit{\theta }}\left({\mathbf{Y}}^{\left(T\right)}\right)f_{{\mathbf{Z}}^{\left(T\right)}|\mathit{\beta },\gamma ,z_0}\left({\mathbf{Z}}^{\left(T\right)}\right)\pi (\mathit{\beta },\gamma )\pi \left(\mathit{\theta }\right)\pi (z_0,{\pi }_0)\hfill \\ & & =\left[\prod _{t=1}^T{f}_{y_t|z_t,\mathit{\theta }}\left(y_t\right)f_{z_t|\mathit{\beta },\gamma ,z_{t-1}}\left(z_t\right)\right]\pi (\mathit{\beta },\gamma )\pi \left(\mathit{\theta }\right)\pi (z_0,{\pi }_0).\hfill \end{array}$$

Note that $\pi (·)$ and $\pi (·|{\mathbf{Y}}^{\left(T\right)})$ denote a prior distribution function and a posterior distribution function given ${\mathbf{Y}}^{\left(T\right)}$, respectively. Note that the full conditional distribution used in the Gibbs algorithm can be easily found through the joint posterior distribution.

The full conditional distribution of logistic regression parameters, $(\mathit{\beta },\gamma )$, is of the form

$$\begin{array}{c}\hfill \pi (\mathit{\beta },\gamma |\mathit{\theta },{\pi }_0,z_0,{\mathbf{Z}}^{\left(T\right)},{\mathbf{Y}}^{\left(T\right)})\propto \left[\prod _{t=1}^T{f}_{z_t|\mathit{\beta },\gamma ,z_{t-1}}\left(z_t\right)\right]\pi (\mathit{\beta },\gamma ),\end{array}$$

where the prior distribution $\pi (\mathit{\beta },\gamma )$ generally is assumed to be non-informative, which is a constant. Note that $\pi (·|·,{\mathbf{Y}}^{\left(T\right)})$ denotes a conditionl posterior distribution function.

The full conditional distribution of data model parameters, $\mathit{\theta }$, can be expressed as

$$\begin{array}{c}\hfill \pi \left(\mathit{\theta }\right|\mathit{\beta },\gamma ,{\pi }_0,z_0,{\mathbf{Z}}^{\left(T\right)},{\mathbf{Y}}^{\left(T\right)})\propto \left[\prod _{t=1}^T{f}_{y_t|z_t,\mathit{\theta }}\left(y_t\right)\right]\pi \left(\mathit{\theta }\right),\end{array}$$

where the prior distribution $\pi \left(\mathit{\theta }\right)$ is also assumed to be non-informative prior, which is a constant.

The full conditional distributions of parameters, ${\pi }_0$ and $z_0$, are given by

$$\begin{array}{c}\hfill \pi \left({\pi }_0\right|\mathit{\beta },\gamma ,\mathit{\theta },z_0,{\mathbf{Z}}^{\left(T\right)},{\mathbf{Y}}^{\left(T\right)})\propto f_{z_0|{\pi }_0}\left(z_0\right)\pi \left({\pi }_0\right)\propto {\pi }_0^{z_0}{(1-{\pi }_0)}^{1-z_0}\pi \left({\pi }_0\right)\end{array}$$

(3)

and

$$\begin{array}{c}\hfill \pi \left(z_0\right|\mathit{\beta },\gamma ,\mathit{\theta },{\pi }_0,{\mathbf{Z}}^{\left(T\right)},{\mathbf{Y}}^{\left(T\right)})\propto f_{z_0|{\pi }_0}\left(z_0\right)f_{z_1|\mathit{\beta },\gamma ,z_0}\left(z_1\right),\end{array}$$

respectively. Note that we here assume $U(0,1)$ as the prior distribution $\pi \left({\pi }_0\right)$. From a non-informative prior assumption, the full conditional distribution of parameter ${\pi }_0$ is a beta distribution of the form

$$\begin{array}{c}\hfill \pi \left({\pi }_0\right|\mathit{\beta },\gamma ,\mathit{\theta },z_0,{\mathbf{Z}}^{\left(T\right)},{\mathbf{Y}}^{\left(T\right)})\sim Beta(z_0+1,2-z_0).\end{array}$$

Note that $Beta(·,·)$ denotes a beta distribution. In addition, since the hidden state space is assumed to be $\{0,1\}$, the probability of the hidden state $z_0$ can be expressed as

$$\begin{array}{c}\hfill P\left[z_0=s|\mathit{\beta },\gamma ,\mathit{\theta },{\pi }_0,{\mathbf{Z}}^{\left(T\right)},{\mathbf{Y}}^{\left(T\right)}\right]=\frac{P(z_0=s|{\pi }_0)f_{z_1|\mathit{\beta },\gamma ,z_0=s}\left(z_1\right)}{{\sum }_{k=0}^{1}P(z_0=k|{\pi }_0)f_{z_1|\mathit{\beta },\gamma ,z_0=k}\left(z_1\right)}\end{array}$$

for $s=0,1$. As a result, the posterior distribution of $z_0$ can be defined as follows:

$$\begin{array}{c}\hfill \pi \left(z_0\right|\mathit{\beta },\gamma ,\mathit{\theta },{\pi }_0,{\mathbf{Z}}^{\left(T\right)},{\mathbf{Y}}^{\left(T\right)})\sim Ber\left[\frac{P(z_0=s|{\pi }_0)f_{z_1|\mathit{\beta },\gamma ,z_0=s}\left(z_1\right)}{{\sum }_{k=0}^{1}P(z_0=k|{\pi }_0)f_{z_1|\mathit{\beta },\gamma ,z_0=k}\left(z_1\right)}\right].\end{array}$$

The full conditional distribution of a hidden state $z_j$ can be expressed as

$$\begin{array}{c}\hfill \pi \left(z_j\right|\mathit{\beta },\gamma ,\mathit{\theta },{\pi }_0,z_0,{\mathbf{Z}}_{-j}^{\left(T\right)},{\mathbf{Y}}^{\left(T\right)})\propto f_{y_j|z_j,\mathit{\theta }}\left(y_j\right)f_{z_{j+1}|\mathit{\beta },\gamma ,z_j}\left(z_{j+1}\right)f_{z_j|\mathit{\beta },\gamma ,z_{j-1}}\left(z_j\right)\end{array}$$

for $j=1,\dots ,T-1$, and

$$\begin{array}{c}\hfill \pi \left(z_T\right|\mathit{\beta },\gamma ,\mathit{\theta },{\pi }_0,z_0,{\mathbf{Z}}^{(T-1)},{\mathbf{Y}}^{\left(T\right)})\propto f_{y_T|z_T,\mathit{\theta }}\left(y_T\right)f_{z_T|\mathit{\beta },\gamma ,z_{T-1}}\left(z_T\right)\end{array}$$

where ${\mathbf{Z}}_{-j}^{\left(T\right)}$ are the unobserved states from time 1 to time t except time j. That is, ${\mathbf{Z}}_{-j}^{\left(T\right)}=(z_1,\dots ,z_{j-1},z_{j+1},\dots ,z_T)$. In the same way of $z_0$, the posterior distribution of $z_j$ can be obtained as follows:

$$\begin{array}{ccc}& & \pi \left(z_j\right|\mathit{\beta },\gamma ,\mathit{\theta },{\pi }_0,z_0,{\mathbf{Z}}_{-j}^{\left(T\right)},{\mathbf{Y}}^{\left(T\right)})\hfill \\ & & \sim Ber\left[\frac{f_{y_j|z_j=1,\mathit{\theta }}\left(y_j\right)f_{z_{j+1}|\mathit{\beta },\gamma ,z_j=1}\left(z_{j+1}\right)P(z_j=1|\mathit{\beta },\gamma ,z_{j-1})}{{\sum }_{k=0}^1{f}_{y_j|z_j=k,\mathit{\theta }}\left(y_j\right)f_{z_{j+1}|\mathit{\beta },\gamma ,z_j=k}\left(z_{j+1}\right)P(z_j=k|\mathit{\beta },\gamma ,z_{j-1})}\right]\hfill \end{array}$$

for $j=1,\dots ,T-1$, and

$$\begin{array}{c}\hfill \pi \left(z_T\right|\mathit{\beta },\gamma ,\mathit{\theta },{\pi }_0,z_0,{\mathbf{Z}}^{(T-1)},{\mathbf{Y}}^{\left(T\right)})\sim Ber\left[\frac{f_{y_T|z_T=1,\mathit{\theta }}\left(y_T\right)P(z_T=1|\mathit{\beta },\gamma ,z_{T-1})}{{\sum }_{k=0}^1{f}_{y_T|z_T=k,\mathit{\theta }}\left(y_T\right)P(z_T=k|\mathit{\beta },\gamma ,z_{T-1})}\right].\end{array}$$

Each full conditional distribution obtained in this way can be used for sampling with the Gibbs algorithm. Note that the full conditional distributions of logistic parameters $(\mathit{\beta },\gamma )$, data model parameters $\mathit{\theta }$, and $z_0$ are not given as a closed form. Thus the M-H algorithm is required for sampling from these distributions.

4.3. Prediction

The posterior distribution sample can be generated in Section 4.2, and the observation model and process model are already assumed in (1) and (2), and then the Bayesian prediction process is carried out in the following steps:

• Generate model parameters $\mathit{\beta },\gamma ,\mathit{\theta },z_T$ from the $\pi (\mathit{\beta },\gamma ,\mathit{\theta },{\mathbf{Z}}^{\left(T\right)}|{\mathbf{Y}}^{\left(T\right)})$.
• Predict the state variable $z_{T+1}$ from the $f_{z_{T+1}|\mathit{\beta },\gamma ,z_T}\left(z_{T+1}\right)$.
• Finally predict $y_{T+1}$ from $f_{y_{T+1}|z_{T+1},\mathit{\theta }}\left(y_{T+1}\right)$.

Note that step 1 can be done through the Gibbs algorithm.

5. Simulation Study

We present the results of a small simulation study designed to investigate the proposed model-based approach in terms of the parameter estimation. We consider two models based on Gaussian mixture and exponential mixture, respectively. Because our application of interest is the effects of the previous hidden state and the degree of mixing of two distributions on the estimations, which is also considered in the simulation, for each model, we execute the following two steps:

• Generate hidden states sequentially based on transition probabilities through the logistic regression.
• Generate a sequence of observations from the mixture distribution corresponding the hidden states.

Let $Y_t$ be the observation at time t and let $Z_t$ be the associated hidden state. Let $x_t$ be the covariates from $N(0,1)$, respectively. We consider the following eight models:

• $y_t|z_t=0\sim N({\mu }_0,{\sigma }_0^2)$, $y_t|z_t=1\sim N({\mu }_1,{\sigma }_1^2)$ and $logit\left(P(z_t=1|z_{t-1})\right)={\beta }_0+{\beta }_1{x}_t+\gamma z_{t-1}$.
  • (${\mu }_0=5$, ${\sigma }_0^2=4$, ${\mu }_1=13$, ${\sigma }_1^2=9$) and (${\beta }_0=1$, ${\beta }_1=2$, $\gamma =0.9$).
  • (${\mu }_0=5$, ${\sigma }_0^2=4$, ${\mu }_1=13$, ${\sigma }_1^2=9$) and (${\beta }_0=1$, ${\beta }_1=2$, $\gamma =0.3$).
  • (${\mu }_0=5$, ${\sigma }_0^2=9$, ${\mu }_1=15$, ${\sigma }_1^2=16$) and (${\beta }_0=1$, ${\beta }_1=2$, $\gamma =0.9$).
  • (${\mu }_0=5$, ${\sigma }_0^2=9$, ${\mu }_1=15$, ${\sigma }_1^2=16$) and (${\beta }_0=1$, ${\beta }_1=2$, $\gamma =0.3$).
• $y_t|z_t=0\sim Exp\left({\theta }_0\right)$, $y_t|z_t=1\sim Exp\left({\theta }_1\right)$ and $logit\left(P(z_t=1|z_{t-1})\right)={\beta }_0+{\beta }_1{x}_t+\gamma z_{t-1}$.
  • (${\theta }_0=2$, ${\theta }_1=5$) and (${\beta }_0=1$, ${\beta }_1=2$, $\gamma =0.9$).
  • (${\theta }_0=2$, ${\theta }_1=5$) and (${\beta }_0=1$, ${\beta }_1=2$, $\gamma =0.4$).
  • (${\theta }_0=2$, ${\theta }_1=3$) and (${\beta }_0=1$, ${\beta }_1=2$, $\gamma =0.9$).
  • (${\theta }_0=2$, ${\theta }_1=3$) and (${\beta }_0=1$, ${\beta }_1=2$, $\gamma =0.4$).

Note that the settings of (${\mu }_0=5$, ${\sigma }_0^2=4$, ${\mu }_1=13$, ${\sigma }_1^2=9$) and (${\theta }_0=2$, ${\theta }_1=5$) mean a mild mixing of two distributions, and the settings of (${\mu }_0=5$, ${\sigma }_0^2=9$), (${\mu }_1=15$, ${\sigma }_1^2=16$), and (${\theta }_0=2$, ${\theta }_1=3$) mean a strong mixing of two distributions. In addition, $\gamma$ controls the effect of the previous hidden state. The eight models were simulated and averaged statistics calculated, with the simulation exercise repeated 200 times.

Table 1. Averaged posterior means (with posterior standard deviations) for parameters based on Gaussian mixtures.

True Parameters							Posterior Means (Posterior Standard Deviations)
${\mathit{\mu }}_{\mathbf{0}}$	${\mathit{\sigma }}_{\mathbf{0}}^{\mathbf{2}}$	${\mathit{\mu }}_{\mathbf{1}}$	${\mathit{\sigma }}_{\mathbf{1}}^{\mathbf{2}}$	${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$	$\mathit{\gamma }$	${\mathit{\mu }}_{\mathbf{0}}$	${\mathit{\sigma }}_{\mathbf{0}}^{\mathbf{2}}$	${\mathit{\mu }}_{\mathbf{1}}$	${\mathit{\sigma }}_{\mathbf{1}}^{\mathbf{2}}$	${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$	$\mathit{\gamma }$
5	4	13	9	1	2	1	5.0219	4.1093	13.1212	8.9132	1.0426	1.9873	0.8852
							(0.1983)	(0.5922)	(0.4342)	(1.4984)	(0.0897)	(0.1464)	(0.1742)
						0.3	4.9823	4.1232	12.8334	8.9342	1.0426	1.9873	0.2942
							(0.2134)	(0.6312)	(0.4724)	(1.6243)	(0.0903)	(0.1483)	(0.1383)
5	9	15	16	1	2	1	5.1423	9.2742	14.8976	16.5872	1.0634	1.9674	0.9853
							(0.2072)	(0.6001)	(0.4872)	(1.9874)	(0.0932)	(0.1512)	(0.1932)
						0.3	5.2014	9.4234	14.8773	16.9532	1.0721	1.9532	0.2912
							(0.2303)	(0.7423)	(0.5023)	(2.1532)	(0.1102)	(0.1772)	(0.1492)

and

Table 2. Averaged posterior means (with posterior standard deviations) for parameters based on exponential mixtures.

True Parameters					Posterior Means (Posterior Standard Deviations)
${\mathit{\theta }}_{\mathbf{0}}$	${\mathit{\theta }}_{\mathbf{1}}$	${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$	$\mathit{\gamma }$	${\mathit{\theta }}_{\mathbf{0}}$	${\mathit{\theta }}_{\mathbf{1}}$	${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$	$\mathit{\gamma }$
2	5	1	2	1	2.1032	4.9882	1.0723	1.991	0.8872
					(0.1072)	(0.2314)	(0.0932)	(0.1498)	(0.1821)
				0.4	2.1032	4.9882	1.0723	1.991	0.4123
					(0.1212)	(0.2532)	(0.0987)	(0.1534)	(0.1431)
2	3	1	2	1	2.1432	2.9432	1.0934	1.9845	0.8773
					(0.1523)	(0.2944)	(0.0963)	(0.1517)	(0.1889)
				0.4	2.1632	2.8982	1.0723	1.9741	0.3889
					(0.1873)	(0.3214)	(0.1083)	(0.1562)	(0.1589)

show the averaged posterior means and posterior standard deviations for parameters based on Gaussian mixtures and exponential mixtures, respectively. It turns out that the stronger the mixing and the weaker the influence of hidden states, the greater the uncertainty in the estimation. However, overall, it demonstrates a good performance.

6. Data Analysis

Matsutake mushrooms are a valuable forest product that can increase rural incomes by up to 100%. They are wild mushrooms that grow naturally in Pinus densiflora forests. In Korea, matsutake mushrooms have not been cultivated due to a lack of agricultural technology. However, the impact of climate on matsutake mushroom yield has not been well studied. Previous studies have attempted to relate matsutake mushroom occurrence to weather patterns in Korea but have not been able to identify the specific meteorological factors that are most influential. Here, we use hidden Markov models (HMMs) to identify hidden states in a specific area of Bonghwa-gun, Gyeongsangbuk-do. We then use these hidden states to identify the meteorological factors that indirectly affect matsutake mushroom yield.

6.1. Data Description

The data for this study are the annual production of matsutake mushrooms (kg) observed in a region of Bonghwa-gun, Gyeongsangbuk-do, from 1997 to 2016. Matsutake mushrooms are harvested in Korea from late August to late October, with the peak harvest period occurring for about 10 days from late September to early October. The production of matsutake mushrooms can vary significantly depending on the weather. Therefore, this study considered the meteorological factors in May, June, July, and August from 1997 to 2016. The meteorological factors were combined and analyzed in terms of time and space. The detailed variables are summarized in

Table 3. Variables and their descriptions.

Variables	Description
m_temp [5/6/7/8]	mean temperature during [May/June/July/August] (°C)
h_temp [5/6/7/8]	maximum temperature during [May/June/July/August] (°C)
l_temp [5/6/7/8]	minimum temperature during [May/June/July/August] (°C)
mdew_temp [5/6/7/8]	mean dew point temperature during [May/June/July/August] (°C)
ms_pre [5/6/7/8]	mean spot atmospheric pressure during [May/June/July/August] (hPa)
msea_pre [5/6/7/8]	mean sea-level pressure during [May/June/July/August] (hPa)
hsea_pre [5/6/7/8]	maximum sea-level pressure during [May/June/July/August] (hPa)
lsea_pre [5/6/7/8]	minimum sea-level pressure during [May/June/July/August] (hPa)
mwat_pre [5/6/7/8]	mean water vapor pressure during [May/June/July/August] (hPa)
hwat_pre [5/6/7/8]	maximum water vapor pressure during [May/June/July/August] (hPa)
lwat_pre [5/6/7/8]	minimum water vapor pressure during [May/June/July/August] (hPa)
sun [5/6/7/8]	amount of sunlight during [May/June/July/August] (hour)
daylight [5/6/7/8]	daylight ratio in [May/June/July/August] (%)
m_humid [5/6/7/8]	mean relative humidity in [May/June/July/August] (%)
l_humid [5/6/7/8]	minimum relative humidity in [May/June/July/August] (%)
precip [5/6/7/8]	total precipitation during [May/June/July/August] (mm)
day_precip [5/6/7/8]	the largest daily precipitation in [May/June/July/August] (mm)
hhour_precip [5/6/7/8]	the largest amount of precipitation in one hour in [May/June/July/August] (mm)
hmin_precip [5/6/7/8]	the largest amount of precipitation in ten minutes in [May/June/July/August] (mm)
m_wind [5/6/7/8]	mean wind speed during [May/June/July/August] (m/s)
h_wind [5/6/7/8]	maximum wind speed during [May/June/July/August] (m/s)
h_instwind [5/6/7/8]	the maximum instantaneous wind speed during [May/June/July/August] (m/s)
lgra_temp [5/6/7/8]	minimum grass temperature during [May/June/July/August] (°C)
mground_temp [5/6/7/8]	mean ground temperature during [May/June/July/August] (°C)

.

6.2. Hierarchical Modeling

We here assume that the production of matsutake mushroom is in one of two states: a lean year or a bumper year. Each hidden state has an independent distribution of total matsutake mushroom production, which is assumed to be a Gaussian distribution. From year to year, persistence varies from state to state because it is governed by state and predictor each time. Let $z_t$ denote the (unobserved) annual year state at time t (i.e., $z_t=1$ for a bumper year, $z_t=0$ for a lean year). Let $y_t$ be the (observed) total matsutake mushroom production amount at time t for $1\le t\le T$. The total matsutake mushroom production amount $y_t$ is conditionally independent given the current year state $z_t$. For a hierarchical modeling, we express the full conditional distribution more specifically.

We consider a Gaussian mixture for an observation $y_t$. That is, $f_1\left(y_t\right)\sim N({\theta }_1,{\sigma }_1^2)$, and $f_0\left(y_t\right)\sim N({\theta }_0,{\sigma }_0^2)$. We assume informative prior $\pi (\mathit{\beta },\gamma )\propto c$, where c is a constant. Then, the full conditional distribution of the logistic parameters $(\mathit{\beta },\gamma )$ is obtained by

$$\begin{array}{c}\hfill \pi (\mathit{\beta },\gamma |\mathit{\theta },{\pi }_0,z_0,{\mathbf{Z}}^{\left(T\right)},{\mathbf{Y}}^{\left(T\right)})\propto \prod _{t=1}^T{\left[\frac{\exp ({\mathit{\beta }}^{\prime }{\mathbf{X}}_t+\gamma z_{t-1})}{1+\exp ({\mathit{\beta }}^{\prime }{\mathbf{X}}_t+\gamma z_{t-1})}\right]}^{z_t}{\left[\frac{1}{1+\exp ({\mathit{\beta }}^{\prime }{\mathbf{X}}_t+\gamma z_{t-1})}\right]}^{1-z_t}.\end{array}$$

Based on the informative prior $\pi \left(\mathit{\theta }\right)\propto c$, where c is a constant, the full conditional distribution of data model parameters $\mathit{\theta }=({\theta }_0,{\theta }_1,{\sigma }_0^2,{\sigma }_1^2)$ is given by

$$\begin{array}{c}\hfill \pi \left(\mathit{\theta }\right|\mathit{\beta },\gamma ,{\pi }_0,z_0,{\mathbf{Z}}^{\left(T\right)},{\mathbf{Y}}^{\left(T\right)})\propto \prod _{t=1}^T\left[z_t{f}_1\left(y_t\right)+(1-z_t)f_0\left(y_t\right)\right].\end{array}$$

Assuming $\pi \left({\pi }_0\right)\sim U(0,1)$, the full conditional distribution of hyperparameter ${\pi }_0$ and initial state $z_0$ can be expressed as

$$\begin{array}{c}\hfill \pi \left({\pi }_0\right|\mathit{\beta },\gamma ,\mathit{\theta },z_0,{\mathbf{Z}}^{\left(T\right)},{\mathbf{Y}}^{\left(T\right)})\sim Beta(z_0+1,2-z_0)\end{array}$$

and

$$\begin{array}{c}\hfill \pi \left(z_0\right|\mathit{\beta },\gamma ,\mathit{\theta },{\pi }_0,{\mathbf{Z}}^{\left(T\right)},{\mathbf{Y}}^{\left(T\right)})\sim Ber\left[\frac{P(z_0=1|{\pi }_0)f_{z_1|\mathit{\beta },\gamma ,z_0=1}\left(z_1\right)}{{\sum }_{k=0}^{1}P(z_0=k|{\pi }_0)f_{z_1|\mathit{\beta },\gamma ,z_0=k}\left(z_1\right)}\right],\end{array}$$

where

$$\begin{array}{ccc}& & f_{z_0|{\pi }_0}\left(z_0\right)f_{z_1|\mathit{\beta },\gamma ,z_0}\left(z_1\right)\hfill \\ & & ={\pi }_0^{z_0}{(1-{\pi }_0)}^{1-z_0}{\left[\frac{\exp ({\mathit{\beta }}^{\prime }{\mathbf{X}}_1+\gamma z_0)}{1+\exp ({\mathit{\beta }}^{\prime }{\mathbf{X}}_1+\gamma z_0)}\right]}^{z_1}{\left[\frac{1}{1+\exp ({\mathit{\beta }}^{\prime }{\mathbf{X}}_1+\gamma z_0)}\right]}^{1-z_1}.\hfill \end{array}$$

Finally, the full conditional distribution of hidden state $z_j$ is of the form

$$\begin{array}{ccc}& & \pi \left(z_j\right|\mathit{\beta },\gamma ,\mathit{\theta },{\pi }_0,z_0,{\mathbf{Z}}_{-j}^{\left(T\right)},{\mathbf{Y}}^{\left(T\right)})\hfill \\ & & \sim Ber\left[\frac{f_{y_j|z_j=1,\mathit{\theta }}\left(y_j\right)f_{z_{j+1}|\mathit{\beta },\gamma ,z_j=1}\left(z_{j+1}\right)P(z_j=1|\mathit{\beta },\gamma ,z_{j-1})}{{\sum }_{k=0}^1{f}_{y_j|z_j=k,\mathit{\theta }}\left(y_j\right)f_{z_{j+1}|\mathit{\beta },\gamma ,z_j=k}\left(z_{j+1}\right)P(z_j=k|\mathit{\beta },\gamma ,z_{j-1})}\right],\hfill \end{array}$$

where

$$\begin{array}{ccc}& & f_{y_j|z_j,\mathit{\theta }}\left(y_j\right)f_{z_{j+1}|\mathit{\beta },\gamma ,z_j}\left(z_{j+1}\right)f_{z_j|\mathit{\beta },\gamma ,z_{j-1}}\left(z_j\right)\hfill \\ & & =\left[z_j{f}_1\left(y_j\right)+(1-z_j)f_0\left(y_j\right)\right]\hfill \\ & & \times \prod _{k=j}^{j+1}{\left[\frac{\exp ({\mathit{\beta }}^{\prime }{\mathbf{X}}_k+\gamma z_{k-1})}{1+\exp ({\mathit{\beta }}^{\prime }{\mathbf{X}}_k+\gamma z_{k-1})}\right]}^{z_k}{\left[\frac{1}{1+\exp ({\mathit{\beta }}^{\prime }{\mathbf{X}}_k+\gamma z_{k-1})}\right]}^{1-z_k}\hfill \end{array}$$

for $j=1,\dots ,T-1$, and

$$\begin{array}{c}\hfill \pi \left(z_T\right|\mathit{\beta },\gamma ,\mathit{\theta },{\pi }_0,z_0,{\mathbf{Z}}^{(T-1)},{\mathbf{Y}}^{\left(T\right)})\sim Ber\left[\frac{f_{y_T|z_T=1,\mathit{\theta }}\left(y_T\right)P(z_T=1|\mathit{\beta },\gamma ,z_{T-1})}{{\sum }_{k=0}^1{f}_{y_T|z_T=k,\mathit{\theta }}\left(y_T\right)P(z_T=k|\mathit{\beta },\gamma ,z_{T-1})}\right],\end{array}$$

where

$$\begin{array}{ccc}& & f_{y_T|z_T,\mathit{\theta }}\left(y_T\right)f_{z_T|\mathit{\beta },\gamma ,z_{T-1}}\left(z_T\right)\hfill \\ & & =\left[z_T{f}_1\left(y_T\right)+(1-z_T)f_0\left(y_T\right)\right]\hfill \\ & & \times {\left[\frac{\exp ({\mathit{\beta }}^{\prime }{\mathbf{X}}_T+\gamma z_{T-1})}{1+\exp ({\mathit{\beta }}^{\prime }{\mathbf{X}}_T+\gamma z_{T-1})}\right]}^{z_T}{\left[\frac{1}{1+\exp ({\mathit{\beta }}^{\prime }{\mathbf{X}}_T+\gamma z_{T-1})}\right]}^{1-z_T}.\hfill \end{array}$$

6.3. Analysis

Before examining the parameters from the posterior distribution, we first perform a regression analysis using the production of matsutake mushrooms as the dependent variable to identify the meteorological factors that explain the variation in matsutake mushroom production. Two weather variables are found to be significant (see

Table 4. Posterior means (Est.) and posterior standard deviations (S.E.) for parameters of regression model.

Variables	Intercept	h_temp6	lwat_pre6
Est.	−240,276.5654	153.3783	8741.5625
S.E.	110,388.0470	36.5569	3981.2952

).

In a comparison for our proposed method, the Altman’s approach can be considered. Altman [11] considered the following model:

• $logit\left(P(z_t=1)\right|P(z_{t-1}=k))={\gamma }_k$ at the time point t.
• Conditional on the hidden state $z_t$, observation $y_t$ is assumed to be normally distributed with mean ${\theta }_t$ and variance ${\sigma }^2$, where ${\theta }_t=\alpha z_t+{\mathit{\beta }}^{\prime }{\mathbf{X}}_t$.

It can be also extended by allowing the transition probabilities to depend on covariates or random effects through logistic regression models.

Table 5. Posterior means (Est.) and posterior standard deviations (S.E.) by the Altman’s approach.

Model for Transition Probability
Variables	${\mathit{\gamma }}_{\mathbf{0}}$	${\mathit{\gamma }}_{\mathbf{1}}$
Est.	0.4167	0.3750
S.E.	0.1263	0.1089
Model for Observation
Variables	Intercept	h_temp6	lwat_pre6
Est.	190,359.78	−6054.89	106.57
S.E.	115,617.48	2408.28	1085.01
Variables	precip8	mground_temp5	$Z_t$
Est.	36.61	2518.04	−51,275.40
S.E.	22.49	3561.73	8208.23

shows the estimated coefficients with their standard error for the model to be compared. In Altman’s approach, only hidden state and maximum temperature during June are found to be significant, which means that most meteorological variables do not directly affect the observed annual production of matsutake mushrooms.

For a Bayesian analysis for the proposed HMM based on logistic regression, two parallel chains are used to check the convergence of the chain. After generating 35,000 samples, the initial 15,000 samples are discarded to eliminate the influence on the initial value, then 1000 samples are extracted by selecting every 20th samples to eliminate the autocorrelation. Gelman–Rubin statistics (G-R) are also checked (Gilks et al. [16]). First, let us look at the sequence of hidden states, $z_t$, corresponding to the matsutake production (see Figure 1). The dashed line represents the annual matsutake mushroom production, and the solid line shows the year hidden state corresponding to the annual matsutake mushroom production. When $z_t=0$, it means a lean year, and when $z_t=1$, it means a bumper year. Second, let us take a look at the logistic parameters, $(\mathit{\beta },\gamma )$. In a logistic regression based on the hidden Markov model (HMM), we use the hidden state in each year as the dependent variable of the logistic regression model to identify the meteorological factors that indirectly explain the variation in matsutake mushroom production (see

Table 6. Posterior means (Est.) and posterior standard deviations (S.E.) for parameters of the proposed model.

Variables	${\mathit{\theta }}_0$	${\mathit{\theta }}_1$	${\mathit{\sigma }}_0$	${\mathit{\sigma }}_1$
Est.	11,136.43	58,555.87	8829.52	14,328.17
S.E.	284.76	395.82	149.51	263.08
Variables	Intercept	h_temp6	lwat_pre6
Est.	−13.6391	0.2367	−0.1378
S.E.	2.4270	0.0482	0.0451
Variables	precip8	mground_temp5	$Z_{t-1}$
Est.	0.0032	0.2520	0.3321
S.E.	0.0005	0.0576	0.1533

). As a result of the analysis, four variables are found to be significant. The analysis shows that the total precipitation in August and the mean ground temperature in May are additional variables that affect matsutake mushroom production, compared to the previous model. In addition, there is a significant shift in the transition probability depending on whether the previous year was a lean year or a bumper year.

7. Concluding Remarks

Covariates-dependent hidden Markov models extended by allowing the transition probabilities to depend on covariates are a powerful class of statistical models and are used to model a wider range of sequential data. In this paper, Markovian properties were achieved by implanting the previous state variables to the logistic regression model. The logistic regression based on the hidden Markov model (HMM) is applied to identify significant variables affecting the annual production of matsutake mushrooms (kg) data observed in the region of Bonghwa-gun from 1997 to 2016. The proposed method differs from the existing analysis methods by using the state variable in logistic regression and the mixture distribution of the state, rather than using the observed value directly in the analysis. It is particularly useful to identify the relevance of variables that are difficult to show in the existing models. As a result, we find additional meteorological factors affecting the annual production of matsutake mushrooms compared to the existing methods. The proposed covariates-dependent hidden Markov models can be useful tools for sequential data in the presence of covariates and used in a variety of applications, including financial time series analysis, medical diagnosis, and customer segmentation. In addition, it can be applied in modeling covariates in HMMs regardless of the dimension of the state variable in the hidden Markov model.