Goodness of fit test for higher order binary Markov chain models

Abstract: When the interest is in making statements about change based on repeated measurements of discrete data, one way to do so is using Markov chain models. Goodness of fit test to find a good model is very important in analyzing the underlying patterns and relationships in the repeated measures data. To test for the various associations in the models, the likelihood ratio and Wald tests are used. However, it has been observed that the efficient score tests can provide equally good tests and can provide an easier alternative. In this paper, we provide an extension of Tsiatis method for goodness of fit test on higher order Markov chains. In our method, we follow the approach of Tsiatis goodness of fit test in logistic regression models. New method provided in this paper is applied to real-life data to examine the suitability of the techniques.


Introduction
Markov chain models are used in various applied fields, such as time series analysis, longitudinal studies, life data, environmental problems. The behavior of a Markov chain depends on the transition

PUBLIC INTEREST STATEMENT
In this paper, an important test procedure arising from repeated measurements of discrete data have been introduced. In real-life situations, the change in the status of a disease or other outcome variables need to be examined. These changes need to be studied to understand the underlying factors influencing transitions in the status during specified time intervals. We can employ Markov models in order to find the relationships between the risk factors and outcome variables. The models can be of first or higher orders. One formidable challenge in employing these models is to confirm goodness of the model for first or higher order. This paper provides a simple test that can be used to test goodness of fit of higher order Markov models with covariate dependence for binary data. The proposed test procedure appears to be very useful technique in providing the goodness of fit of higher order binary Markov chains. matrix, which contains transitional probabilities. In most practical studies, the transition matrix is unknown and needs to be estimated. Several methods are available for the estimation and test procedure of transition probabilities. However, most researchers have worked primarily on the estimation of parameters and only a few reports on test procedures. One of the most important tests on Markov chain models is the stationarity of transition probabilities and the goodness of fit of Markov chain models. This section presents a brief summary of the tests performed on Markov chain. Anderson and Goodman (1957) obtained the maximum likelihood estimates and their asymptotic distribution for the transition probabilities in a Markov chain of arbitrary order with repeated observations of the chain. The likelihood ratio tests and chi-square tests used in contingency tables were obtained for testing these hypotheses. Billingsley (1961) used Whittle's formula, chi-square, and maximum likelihood methods to test for stationarity and order of the higher order Markov chain. Mcqueen and Thorley (1991) used Markov chain to analyze annual stock returns. Albert (1994) proposed a class of Markov models for analyzing sequences of ordinal data from a relapsing-remitting disease, where the state space was expanded to include information about the ordinal severity score as well as the relapsing-remitting status. He proposed a parameterization that can reduce the number of parameters. It is noteworthy that most of these research works have been conducted for estimating parameters based on the first-order Markov chain. Recently, several new methods for higher order Markov chains have been reported, where the estimation and test procedures became quite complex due to the increased order of the models (Chowdhury, Islam, Shah, & Al-Enezi, 2005;Islam & Chowdhury, 2006;Islam, Chowdhury, & Briollais, 2012;Rahman & Islam, 2007).
However, less effort is given towards studying the field of covariate-dependent Markov models (Muenz & Rubinstein, 1985;Yi, He, & Liang, 2009). In this paper, a test procedure for the goodness of fit of a binary Markov chain model is proposed by extending Tsiatis' procedure (Tsiatis, 1980). The proposed test was extended for the second-and higher order of the Markov chain model. The efficient score test was used for testing null hypotheses, which only required the estimate of parameters under true null hypothesis. The proposed model and test procedures were thoroughly examined using a set of data for the elderly population and employing simulations. Sirdari, Islam, and Awang (2013) proposed the goodness of fit test for higher order binary Markov chain models based on marginal distribution. The problem with this proposal was that marginal distribution has limited assumptions because of the correlations between variables, which are not easy to estimate. Thus, the proposed model in this study was based on the conditional transition probabilities, which means that there is no correlation between variables. Tsiatis (1980) proposed a goodness of fit test for the logistic regression model. In terms of binary data analysis, this model relates the probability of a response to a set of covariates 1 , … , p according to Equation (2.1):

A brief overview of the test proposed by Tsiatis (1980)
where p denotes the conditional probability of response given by the vector, = 1 , … , p , χ 0 = 1, and � = 0 , … , p denotes the regression coefficients. The space of covariates 1 , … , p is partitioned into k distinct region in p-dimensional space, denoted by R 1 , … , R k . The indicator functions, I (j) j = 1, … , k , are defined by I (j) = 1 if 1 , … , p ∈ R j and I (j) = 0.
Tsiatis considered the following model in Equation (2.2): where I � = I (1) , … , I (k) and � = 1 , … , k . The goodness of fit test consists of testing the This test is based on the efficient score test, as represented by Equation 2.3: where Z ′ is the k-dimensional vector l∕ 1 , … , l∕ k and l denotes the log-likelihood. The k × k matrix, V is equal to: where All previous terms were evaluated at = 0 and j =̂ j , where ̂ j is the maximum likelihood estimate of the parameters when is true.

Goodness of fit test of first-order Markov chains
Consider the case of a single stationary process, Y 1 , … , Y T , generated by a binary Markov chain that uses values of 0 and 1. The transition matrix is defined by: The transition probabilities, p jt , can be modeled using logistic regression, as shown by Model (3.1): Vector t contains covariates and it is equal to t = 1, t1 , … , tp . j is the vector of parameters, j = j0 , j1 , … , jp . The likelihood function that corresponds to Model (3.1) is: where n 00t , n 01t , n 10t , and n 11t are the number of transitions of each type observed at time t. The log-likelihood is as shown by Equation (3.2): The log-likelihood can be shown as, l = ln L = ln L 0 + ln L 1 , where It was assumed that the space of covariate t1 , … , tp was partitioned into G distinct regions in Then, for a binary Markov chain, the following Model (3.3) was considered: t and � j = j1 , … , jG is an arbitrary covariate vector. This test was performed by testing the null hypothesis, H 0 : j1 = ⋯ = jG = 0. This hypothesis was proposed by partitioning the space of covariates into distinct regions and calculating a test statistic, which was a quadratic form of the observed counts, excluding the expected counts.
The efficient score test and the likelihood ratio test were also used. Both statistics have asymptotic chi-square distribution, with G degrees of freedom, as proven by Rao (1973). The current test used in this study was based on the efficient score test because it only requires an estimate of j under the null hypothesis, whereas the likelihood ratio statistics needs an estimate of j under the alternative model. The test statistics is defined by Equation (3.4): All previous terms were evaluated at j = 0 and =̂ , where ̂ is the maximum likelihood estimate of the parameters when H 0 is true. Using the standard likelihood theory, a test could be extracted in the quadratic form of observed counts minus expected counts, and whose large sample properties are easily established.
It is evident that the log-likelihood for Model (3.3) can be achieved by inserting Equation (3.2), as follows: The kth element of vector t used in the computation of Equation (3.4) is the partial derivative of where O jk and E jk are the observed and expected numbers of responses in the kth region. Therefore, Equation (3.4) is the quadratic form of the vector of observed counts minus expected counts.
Quantities necessary for computing the covariance matrix, V j , j = 0, 1 are as follows: where, ξ k denotes the set of indices t, such that

Extension of the model for higher-order Markov chains
Consider the nth-order Markov model for times, t − n, t − (n − 1), … , t − 1, and t, with transition matrix, P, and its components: The logistic regression model for p r⋯sjt is: where vector, t = 1, x t1 , … , x tp , contains covariates and l⋯sj = r⋯sj0 , r⋯sj1 , … , r⋯sjp is the vector of parameters.
The likelihood function corresponding to this model is as follows: The log-likelihood is: where for r, s, j = 0, 1, Then, Model (3.3) can be extended for order n, which can be written as shown by Equation (4.1): The related null hypothesis is H 0 : r⋯sj1 = … = r⋯sjG = 0, and the test statistic is shown by Equation where Z is a 2 n -dimensional vector with elements of The log-likelihood based on Model (4.1) is as follows: The partial derivative of l r⋯sj , l, s, j = 0, 1 with respect to r⋯sjk , r, s, j = 0, 1 at r⋯sj = 0 and r⋯sj =̂ r⋯sj ,

Application
We applied the proposed test on the Health and Retirement Study (HRS) (2009) data to demonstrate its application. This is a longitudinal household survey data-set for the study of retirement and health among the elderly in the United States. The RAND Centre collected these data to study aging, with funding and support from the National Institute on Aging (NIA) and the Social Security Administration (SSA). These data were collected from 1992 to 2006 in eight waves for 30,405 people. We considered individuals who attended the program in 1992 and then, followed up until 2006. The study was about depression among individuals (0 for no depression and 1 for depression), and age (yearly), gender (0 for male and 1 for female), body mass index (BMI), and drinking (0 for not drinking and 1 for drinking), which were considered as covariates. The space of covariate t1 , … , tp was partitioned into four distinct regions: (male and not drinking); (male and drinking); (female and not drinking); and (female and drinking). Some of these variables may contain missing values because the referenced person did not respond to the waves. Thus, we had to drop the ID of individuals from all waves if there were missing values for these covariates. There were 668 missing values in the covariates, which included 353 IDs, i.e. these individuals responded for the outcome variable, but not for the covariates. Thus, 353 IDs were dropped from the data in this work. Additionally, S-Plus functions modified by Chowdhury et al. (2005) were developed and used to estimate the parameters of the model. The Newton-Raphson method was used in this program for parameter estimation. Table 1 shows the different types of transition counts for the first-and second-order transitions. Meanwhile, Table 2 shows the estimated values for the covariate-dependent Markov models for different types of transitions. The results are for the first-and second-order Markov models.
Billingsley's chi-square statistics were computed using  Estimates of the parameters for the first-order transitions demonstrated negative association between the transitions from no depression to depression with age and drinking, while positive associations were obtained with BMI and sex (females have higher risks). Those who did not change their status from depression were found to be associated negatively with age and positively with drinking. Both the Billingsley and the proposed tests showed that the first-order models can be accepted. Similarly, the second-order models indicated that age and drinking were negatively associated, while sex was positively associated for the 0-0-1 type of transition. Sex and BMI were positively associated and drinking was negatively associated for the 1-0-1 transition type, while age was negatively associated, but drinking was positively associated for the 1-1-1 transition type. Interestingly, the second-order models also appeared to have good fit in favor of the null hypothesis. These results demonstrated that both the first-and second-order models could be employed for the given set of data.

Simulation
Data generated by the techniques provided by Ghosh and Mukerjee, and Leisch et al. was used to examine the suitability of the proposed models. In these techniques, bindata package in R were employed for generating correlated binary data. First, data were generated from the multivariate normal random variables, and then, they were transformed into binary data. In this study, two variables were generated as the outcome variables at time t and t − 1 for the first-order Markov model, with various combinations of probabilities of occurring 1 and 0 to obtain different correlations. These results were used to compare the models under independent and selected values of measure of association. For models 1, 2, and 3, the data were generated based on correlation of 0.4 between the outcome variables at time t − 1 and t for the first-order, and at time t − 2, t − 1, and t for the secondorder. Similarly, models 4, 5, and 6 were generated with correlation of 0, while models 7, 8, and 9 considered correlation of −0.4. For each model, four covariates were also generated, corresponding to the correlated response variables by considering different correlations with the outcome variables. These estimates and tests were repeated 500 times for all models, and for sample sizes of 250, 500, and 1,000 for different correlations between outcome variables. The models that were used for this simulation study were different applications of the conditional model verified in Model (3.3 was used to compare the obtained results, with and without covariates. In other words, the Markov models were estimated using covariates and were employed in this test. Table 3 shows the simulation results for the first-order model, which included frequencies by transition type, correlation between outcome variables in the bivariate Bernoulli population, average estimates of the parameters, and the number of rejected hypotheses in 500 times of simulation for these models using Billingsley's test and the proposed test as an extension of Tsiatis' test. Acceptance of the null hypothesis, H 0 : j1 = ⋯ = jG = 0, would indicate a good fit of Model (4.1) to the data. The percentage of rejection for Billingsley's test was 0 because the test procedure did not consider any covariate. However, in the proposed extension, models with covariate dependence were used. Hence, it can be concluded that the proposed test statistics depended on the covariate-dependant transition probabilities, where the selection of appropriate variables in the model may influence the goodness of fit, to a large extent. This observation implied that the proposed test may display deviations for a good fit in some instances. In other words, the goodness of fit test proposed by this study, as an extension of the Tsiatis' test, depended on the model's specifications in terms of the explanatory power of the selected variables. Based on the estimated covariates for the first-order transition from 0 to 1, we observed a positive association with variable 1, and a transition from 1 to 1, which has negative association with variable 1 for all models, except with model 4 (sample size of 250, with correlation of 0). The rejection percentage had varied for the first-order model, mainly in the range of 4.-6.6%. This result showed that the proposed test method was satisfactory for different sizes of samples, with different correlation of outcome variables based on the first-order Markov chain model.
The simulation results for the second-order model are given in Table 4. The table shows      for the models. Results for the second-order models showed that there was no association between covariates and outcome variables, which was expected because of the higher order of the underlying Markov chain. Only two models, 3 and 6, with sample size of 1000 and correlations of 0 and 0.4, have negative associations with variable 2 in different types of transitions. The range of rejection percentage of the null hypothesis for the extended Tsiatis' test was 3.6-6.2% for models 1-9 in the second-order. Thus, these models were acceptable for different sizes of samples and different correlations between outcome variables. Results from Billingsley's test were compared with results from the proposed extension of Tsiatis' test; the number of rejected null hypothesis was zero for Billingsley's test because it does not depend on covariates. The number of rejected null hypothesis for model 3 was the lowest.