- Email: [email protected]
A computer program to estimate the parameters of covariate dependent higher order Markov model
Computer Methods and Programs in Biomedicine (2005) 77, 175—181
A computer program to estimate the parameters of covariate dependent higher order Markov model Rafiqul Islam Chowdhurya,∗, M. Ataharul Islamb, Makhdoom Ali Shahc, Naser Al-Enezid a
Department of Health Information Administration, Kuwait University, P.O. Box 31470, Sulaibekhat 90805, Kuwait b Department of Statistics and OR, College of Science, King Saud University, Saudi Arabia c Department of Health Information Administration, Kuwait University, Kuwait d Department of Health Information Administration, Kuwait University, Kuwait Received 10 October 2003 ; received in revised form 21 August 2004; accepted 5 October 2004 KEYWORDS Longitudinal data; Covariate dependent; Higher order; Markov model; S-plus; Computer program
Summary This paper presents a computer program developed in S-plus to estimate the parameters of covariate dependent higher order Markov Chain and related tests. The program can be applied for two states Markov Chain with any order and any number of covariates depending on the PC capabilities. The program provides the maximum likelihood estimates of the parameters, together with their estimated standard error, t-value and significance level. It also produces the test results for likelihood ratio and model chi-square. To illustrate the program we have used a longitudinal data set on maternal morbidity of rural women in Bangladesh. The occurrences of haemorrhage, convulsion, or fits at different follow-ups were used as outcome variable. Economic status, wanted pregnancy, ages at marriage, and education of women were used as covariates. © 2004 Elsevier Ireland Ltd. All rights reserved.
1. Introduction Different statistical models are available to analyze the longitudinal data emerging from biological, epidemiological, demographic, sociological, econometric, and many other fields of study. Markov chain * Corresponding author. Tel.: +965 9503456; fax: +965 4830937. E-mail address: rafi[email protected] (R.I. Chowdhury).
model is one of the statistical techniques that play an important role to analyze such types of data [1—3]. In recent years, considerable interest has been witnessed in the development of multivariate models based on Markov chain and these models are employed to analyze the data generated from meteorology, reliability, survival analysis, etc. Muenz and Rubinstein introduced a first-order discrete time Markov chain for expressing the transition probabilities in terms of function of covari-
0169-2607/$ — see front matter © 2004 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2004.10.003
176
R.I. Chowdhury et al.
ates for a binary sequence of presence or absence of disease and suggested that this model can be extended for higher order Markov Chain Model [4]. They employed logistic regression models to analyze the transition probabilities from one state to another. The estimation of first-order Markov model is quite straightforward [5]. As the classical Markov chain model is a first order that is a process without memory, there are growing interests for higher order Markov models [6—9]. Islam et al. generalized the Muenz and Rubinstein models for covariate dependent higher order Markov model for binary outcomes [10]. They also discussed issues related to extension to more than two intercommunicating states, heterogeneity, and appropriate estimation and test procedures. Available statistical software do not always include procedure to estimate any newly developed model, which prevents researchers to apply the models for complexity in computational procedures. However, there are some softwares which have options to estimate transition matrix and some programs allow to estimate the Markov chain models with different specifications [11—13]. This paper describes a computer program designed to estimate the parameters of covariate dependent higher order Markov model for binary outcomes proposed by Islam et al., and the related test procedures [10]. The program was written using S-plus for windows and a follow-up data set on maternal morbidity of pregnant women was used to demonstrate the program.
2.1. Covariate dependent higher order model Following is the covariate dependent higher order Markov model proposed by Islam et al. [10] extending the Muenz and Rubinstein [4] model for first-order Markov chain. To illustrate the extension of the model for higher order, let us consider a second-order Markov model. The second-order Markov model for time points tj−2 , tj−1 , and tj with corresponding outcomes Yj−2 , Yj−1 , and Yj , respectively, is shown below: Y j−2
π01 (Y j = 1/Y j−1 = 0, X) =
eβ0 X
(1)
1 + eβ0 X
and
π11 (Y j = 1/Y j−1 = 1, X) =
eβ1 X 1+e
β1 X
.
(2)
where Xi = [1, X i1 , . . . , X ip ] is the vector of covariates for the ith person; 0 = [β00 , β01 , . . . , β0p ] the vector of parameters for the transition from 0; and 1 = [β10 , β11 , . . . , β1p ] is the vector of parameters for the transition from 1.
|
Yj
0
1
0
0
π000
π001
0
1
π010
π011
1
0
π100
π101
1
1
π110
π111
.
Using the definitions of vectors of covariates and parameters mentioned in the previous section, we can define the transition probabilities:
πuvw (Y j = 1/Y j−2 = 0, Y j−1 = 0, X) =
eβ001 X
1 + eβ001 X
u = 0, 1; v = 0, 1; w = 0, 1.
, (3)
It may be noted here that π000 + π001 = 1, π010 + π011 = 1, π100 + π101 = 1, and π110 + π111 = 1. Then the likelihood function, as a generalization can be expressed as follows:
2. Computational methods and theory The transition probabilities for the first order can be defined in terms of function of the covariates as follows:
Y j−1
L=
ni n
[{π000 }δ000ij {π001 }δ001ij ][{π010 }δ010ij
i=1 j=1
× {π011 }δ011ij ][{π100 }δ100ij {π101 }δ1011ij ][{π110 }δ110ij × {π111 }δ111ij ],
(4)
where L=L1 L 2 L3 L 4 . The first derivatives with respect to parameters are n
n
i ∂ ln L 1 = ∂βuvwq
i=1 j=1
× X iq δ000ij − (δ000ij + δ001ij )
eβ001 X
1 + eβ001 X
(5)
u = 0,1; v = 0,1; w = 0,1; q = 0, 1, 2, . . ., p. Similarly, we can differentiate ln L2 , ln L3 , ln L 4 , with respect to the corresponding parameters.
A computer program to estimate the parameters of covariate dependent Markov model We can solve for the sets of parameters by equating first derivatives to zero. Each set consists of (p + 1) parameters and hence the total number of parameters to be estimated here is 4(p + 1). The second derivatives for L1 are: ∂2 ln L 1 ∂βuvwq ∂βuvwl =−
ni n
vector: β = [β1 , β2 , . . . , β2k ] = [β , β , . . . , β k where βm mp ], m = 1, 2, . . ., 2 . m0 m1 To test the null hypothesis H0 :β = 0, we can employ the usual likelihood ratio test
−2[ln L(β0 ) − ln L(β)] ≈ χ22k p .
[X iq X il {(δuv0ij + δuv1ij )πuv0 (X)πuv1 (X)}],
i=1 j=1
(6)
where u = 0, 1; v = 0, 1; w = 0, 1. We can differentiate similarly for ln L2 , ln L3 , and ln L4 . Inverse of the (−1) (second derivative) provide estimates of the variance—covariance for the respective estimates of the parameters. Then the likelihood function for a Markov chain of order k, as a generalization, can be expressed as follows:
To test the significance of the qth parameter of the mth set of parameters, the null hypothesis is H0 :βmq = 0 and the corresponding Wald test is W=
ˆmq β . ˆmq ) se(β
Another alternative test procedure is illustrated here on the basis of the suggestions made by Billingsley [14] and Raftery and Tavarey [6]. The test statistic for the k-order Markov model is defined as:
k
L=
ni 2 n
[{πm }δmij {1 − πm }1−δmij ].
χ2 =
The log-likelihood function can be shown as ln L = ln L 1 + ln L 2 + · · · + ln L m + · · · + ln L 2k .
(8)
Here, Lm =
{n(j − k, j − k + 1, . . . j)
(7)
i=1 j=1 m=1
ni n
177
[{πm }δmij {1 − πm }1−δmij ]
(9)
i=1 j=1
The first derivatives with respect to parameters are ni n ∂ ln L m eβm X i X iq (1 − δmij ) − , = ∂βmq 1 + eβm X i i=1 j=1
m = 1, 2, . . . , 2k
and
q = 0, 1, 2, . . . , p.
(10)
The second derivatives for ln Lm are: n
n
i ∂2 ln L m =− [X iq X il {πm (1 − πm )}]. ∂βmq ∂βml
(11)
i=1 j=1
2.2. Testing for the significance of parameters The tests are discussed in Islam et al. [10]. Let us give an overview of the tests employed for testing the significance of the parameters. The vectors of 2k sets of parameters for the kth order Markov model can be represented by the following
− e(j − k, j − k + 1, . . . , j)}2 e(j − k, j − k + 1, . . . , j)
the sum is over all n(j − k, j − k + 1, . . . , j) for n(j − k, j − k + 1, . . . , +). The estimated frequencies can be obtained from the following relationship: e(j − k, j − k + 1, . . . , j) = n(j − k, j − k + 1, . . . , +)ˆ π(j − k, j − k + 1, . . . , j)
2.3. Fitting algorithm The following is an algorithm for fitting higher order Markov model from Eq. (7) using the Gauss—Newton approach: (1) Input data file, number of variables, order of Markov chain, and maximum number of iteration. (2) Read data file, initialize variables, calculate sets of models to estimate (e.g., for second order there are 22 = 4 sets of model, for third order 23 = 8 set of models and so on). (3) Compute the combination of outcomes needed to check different types of transitions for different sets of model. For example, for the second order we have to check the following combinations for outcomes for four sets of model:
178
R.I. Chowdhury et al. Y j−2
(4) (5) (6) (7) (8) (9)
Y j−1
0
0
| 0 π000
Yj 1
0
1
π010
1
0
π100
π011 π101
1
1
π110
π111
(11) Estimate  using the following: ˆj = β ˆj−1 − I −1 (β ˆj−1 )U(β ˆj−1 ) β
π001
(12) Update β with the latest value. (13) Calculate test statistics. (14) Iterate until convergence.
We need first row of the above matrix to check the ith individual for δ000ij = 1 if a transition type 0 0 0 is observed during jth follow-up for the ith individual, δ000ij = 0, otherwise; second row will need to check δ011ij = 1 if a transition type 0 1 1 is observed during jth follow-up for the ith individual, δ011ij = 0, otherwise; third row will need to check δ101ij = 1 if a transition type 1 0 1 is observed during jth follow-up for the ith individual, δ101ij = 0, otherwise; and last row will need to check δ111ij = 1 if a transition type 1 1 1 is observed during jth follow-up for the ith individual, δ111ij = 0, otherwise. Start the loop for iteration. Set initial value for parameter to 0 (β = 0). Start the loop for sets of models. Check the transition types for each set of model defined in (3). Compute the first derivatives and second derivatives for each set of model and store separately for each model. To get the score vector, combine the first derivative for all the models in a single vector. For example, if we have three covariates and a second-order Markov Chain then we will have a score vector of 16 elements constant + 3 covariates = 4 and 4 sets of model). If U1 , U2 , U3 , and U4 are the score vectors for four sets of model each with four elements, then we will get the following combined score vector: U(β) = [U 1 U 2 U 3 U 4 ]1×16
(10) Combine the second derivative off all the models as diagonal terms of the matrix, which will be the information matrix (for example, if we have 3 covariates and a second-order Markov Chain, the information matrix will be 16 × 16. Constant + 3 covariates = 4 and 4 sets of model). If I1 , I2 , I3 , and I4 are the 4 × 4 information matrix for four sets of model then we will get the following 16 × 16-combined information matrix: I1 0 0 0 0 I 0 0 2 I(β) = 0 0 I3 0 0
0
0
I4
16×16
3. Program description The function (markov.gen) can be run under S-plus command environment. It can be used for equal and unequal number of follow-ups. The function is written in an ASCII text file and could be copied and pasted in the S-plus command prompt to run [15]. To use the function five arguments are needed. Arguments are data file name including path (alphanumeric), total columns (numeric) in the data file, order of the Markov Chain (numeric, e.g., for first-order 1, for second-order 2, and so on), form of the transition probabilities (numeric, currently it is 1 for logistic form), and maximum number of iteration (numeric). The structure of the data file is ASCII free format. The first column in the data file should be case identification number; the second column is the follow-up number for each case. The third column is the outcome variable and the fourth and subsequent columns are the covariates. In addition, the first row in the data file should contain the entire variable name. Any software (e.g., SPSS, SAS, S-plus, etc.) can be used to prepare the data file in the required format to use our function for Markov Chain. To test the accuracy of the estimates of our program we used a hypothetical data set. The data set contains three observations with 12, 8, and 5 follow-ups and three covariates. Using the hypothetical data set the necessary calculations to estimate the parameter of our model were done manually. These results were compared with the results from the sample run of our S-plus program.
4. Sample run The application is based on the data from the Follow-up Survey on Maternal Morbidity in Bangladesh, conducted by the Bangladesh Institute of Research for Promotion of Essential and Reproductive Health and Technologies (BIRPERTH). The data were collected from November 1992 to December 1993. A total of 1020 pregnant women were included in the study (pregnancy less than 6 months). These subjects were followed-up with
A computer program to estimate the parameters of covariate dependent Markov model
Table 1
Distribution of follow-ups by occurrence of morbidity
Follow-up number
Occurrence of morbidities No n
1 2 3 4 5 6 7 8 Total
565 623 534 418 245 101 23 — 2509
Total Yes
% 57.0 67.9 69.3 70.4 66.2 68.2 67.6 — —
an interval on an average of 1 month, through full-term pregnancy, delivery and till 90 days postpartum period or 90 days after any other pregnancy outcome. The information on socio-economic background, pregnancy related care and practice, extent of morbidity during the index pregnancy, delivery and postpartum period, or abortions were collected. Out of 1020 women, 993 had at least one antenatal follow-up, and 1005 had information on pregnancy termination. Finally, 1006 had at least one postpartum follow-up. However, for our application, we used data for the antenatal follow-ups only. Table 1 shows the occurrence of morbidity at different follow-ups. The data comprise a total of 3827 follow-ups. There were maximum eight follow-ups, which was only for one woman. Dependent variable: Y it is defined as follows: Y it
179
= 1 if there was occurrence of haemorrhage, or convulsion or fits; = 0, otherwise
n
%
427 294 237 176 125 47 11 1
43.0 32.1 30.7 29.6 33.8 31.8 32.4 100.0
992 917 771 594 370 148 34 1
1318
—
3827
Table 2
Transition matrix
Transition type
→0
→1
0→0→0 1→0→0 0→1→0 0→0→1 1→1→0 1→0→1 0→1→1 1→1→1
459 91 31 49 40 26 8 44
72 30 16 27 28 20 23 132
On the S-plus command prompt the following command was used to invoke our function to estimate the parameters of third order Markov Chain with logistic link function: • markov.gen(‘‘c:/markov/marapp.dat’’,7,3,1,12) The model was converged after five iterations. The program produced the transition matrix and estimate of the parameter for each set of equations, which is presented in Table 2. The sets of equations follow the same order as the sequence of different transition types in transition matrix. It also computes the model chi-square, likelihood ratio, AIC, BIC, and Billingsley statistic, which are presented in Table 3.
the covariates were categorized and used as follows: economic status: good = 1, bad = 0; wanted pregnancy: yes = 1, and no = 0; education: no education = 0, primary or higher education = 1; and age at marriage: less than or equal to 15 years = 0, more than 15 years = 1.
5. Hardware and software specifications
We have already conducted several studies on maternal morbidity in Bangladesh using the same set of data for examining different aspects of maternal morbidity. Some of the relevant papers are Islam et al., Chakraborty et al., Chowdhury et al., and Bari et al. [16—21].
The function for Markov Chain is written in S-plus 6.0 and has approximately 325 lines of code. The program requires intensive computer resources, mainly to evaluate the likelihood function. The contribution to the likelihood function of two outcomes, in each case depending on the order of the
180
Table 3
R.I. Chowdhury et al.
Estimates of the parameters of third order Markov Chain
Characteristics and different transition types
Estimate
Standard error
t-value
p-value
0 → 0 → 0 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
−1.6804 0.2478 −0.1146 −0.1914 −0.1528
0.2548 0.3369 0.2774 0.2760 0.2602
−6.5945 0.7354 −0.4133 −0.6933 −0.5872
0.0000 0.4621 0.6794 0.4881 0.5570
1 → 0 → 0 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
−0.6063 −0.2463 −0.8580 0.6512 −0.5059
0.3613 0.6402 0.4387 0.4844 0.4598
−1.6781 −0.3847 −1.9556 1.3444 −1.1003
0.0933 0.7005 0.0505 0.1788 0.2712
0 → 1 → 0 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
−0.8476 −1.1475 0.2453 1.0430 −0.3283
0.5932 0.9642 0.6454 0.7268 0.6915
−1.4289 −1.1901 0.3801 1.4350 −0.4747
0.1530 0.2340 0.7038 0.1513 0.6350
0 → 0 → 1 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
0.2030 −0.0303 −0.6704 −0.5273 −0.5190
0.4478 0.6957 0.5028 0.5987 0.5298
0.4534 −0.0435 −1.3334 −0.8808 −0.9796
0.6503 0.9653 0.1824 0.3784 0.3273
1 → 1 → 0 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
−0.2380 1.2853 −0.4223 −0.1862 −0.2244
0.4623 0.6418 0.5361 0.5513 0.5527
−0.5147 2.0026 −0.7876 −0.3378 −0.4061
0.6068 0.0452 0.4309 0.7355 0.6847
1 → 0 → 1 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
0.0544 1.5457 −0.1858 −0.0901 −1.0149
0.5271 0.9456 0.6468 0.6757 0.7107
0.1032 1.6346 −0.2873 −0.1334 −1.4280
0.9178 0.1021 0.7738 0.8939 0.1533
0 → 1 → 1 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
1.4462 −0.1631 −0.8268 1.0574 0.5812
1.1843 1.0135 1.2022 1.1807 1.2298
1.2211 −0.1609 −0.6878 0.8956 0.4726
0.2221 0.8722 0.4916 0.3705 0.6365
1 → 1 → 1 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
1.4786 0.2558 −0.5342 −0.4100 0.1508
0.3299 0.4584 0.3743 0.3734 0.4103
4.4820 0.5581 −1.4274 −1.0980 0.3675
0.0000 0.5768 0.1535 0.2722 0.7133
Model chi-square (40 d.f.) Likelihood ratio AIC BIC BILLINGSLEY
403.08 445.05 1154.32 1350.35 1410.95
0.000 0.000 0.000 0.000 0.000
A computer program to estimate the parameters of covariate dependent Markov model Markov Chain, requires the evaluation of the transition probability matrix for a given set of covariates. The CPU time needed to run this function depends on the order of the Markov Chain, number of observations, the number of covariates, and the speed of PC processor and the amount of PC memory.
6. Mode of availability The full program ‘kern0pt‘markov.gen’’ is available on request. Any one, who is interested, can request the corresponding author for the complete program. The program file and the necessary instruction for users will be sent as e-mail attachment.
7. Conclusions In this paper we explained the development of S-plus program for application to real life data, which contains necessary instructions to use our S-plus function to estimate the parameter of covariate dependent higher order Markov Chain, using a follow-up data set. Despite the computational complexities, with the use of current state of high quality software like S-plus (StatSci Division, MathSoft Inc.), the estimates of parameters pose no problem. Work is in progress to refine and add new link functions and some other tests. We hope other researchers who are interested in using Markov Chain will use the function developed by us to estimate the parameter of the model.
Acknowledgements We gratefully acknowledge the permission of Dr. Halida Hanum Akhter, Director, BIRPERHT, for using the data for the application in this paper. The authors would like to thank Mahbub E. Elahi, K. Chowdhury, and Arindom Sen of the Bangladesh Institute of Research for Promotion of Essential and Reproductive Health and Technologies (BIRPERHT) for their assistance during different phases of this work. The authors are greatly indebted to the Ford Foundation for funding the data collection of the maternal morbidity study.
References [1] P.S. Albert, A Markov model for sequence of ordinal data from a relapsing-remitting disease, Biometrics 50 (1994) 51—60.
181
[2] P.S. Albert, A.W. Myron, A two state Markov Chain for heterogeneous transitional data: a quasilikelihood approach, Stat. Med. 17 (1998) 1481—1493. [3] R. Kay, A Markov model for analyzing cancer markers and disease states in survival studies, Biometrics 42 (1986) 855—865. [4] L.R. Muenz, L.V. Rubinstein, Markov models for covariate dependence of binary sequences, Biometrics 41 (1985) 91—101. [5] N. Harun, Higher order Markov Chain model for analyzing covariate dependence of longitudinal data. A Thesis Submitted in Partial Fulfillment of the Degree of M.Sc., Department of Statistics, University of Dhaka, 2002. [6] A.E. Raftery, S. Tavare, Estimating and modeling repeated patterns in higher order Markov Chains with mixture transition distribution, Appl. Stat. 43 (1) (1994) 179— 199. [7] A.E. Raftery, A model for higher order Markov Chains, J. R. Stat. Soc. 47 (1985) 1528—1539. [8] F.M. Bass, M.M. Givon, M.U. Kalwani, D. Reibstein, G.P. Wright, AN investigation into the order of the brand choice process, Marketing Sci. 3 (1984) 267—287. [9] J.A. Logan, A structural model of the higher-order Markov process incorporating reversion effects, J. Math. Sociol. 8 (1981) 75—89. [10] M.A. Islam, R.I. Chowdhury, A. Baharum, A higher Markov Model for analyzing covariate dependence, in: Proceedings of the National Conference of Management Science and Operations Research, Melaka, Malaysia, 24—25 June, 2003, pp. 25—33. [11] SAS Institute Inc., SAS/STAT User’s Guide, version 6, vol. 2, fourth edition, SAS Institute Inc., Cary, NC, 1989. [12] MARCH 1.1 User Guide, Markov chains Computation for Homogeneous and Non-Homogeneous data. Universite de Lausanne, Switzerland, 2001. [13] G. Marshall, W. Guo, R.H. Jones, MARKOV: a computer program for multi-state Markov models with covariables, Comput. Methods Program Biomed. 47 (1995) 147—156. [14] P. Billingsley, Statistical Inference for Markov Process, The University of Chicago Press, Chicago, USA, 1961. [15] S-plus 6 for Windows, Programmers guide Insightful Corporation, Seattle, Washington, USA, 2001. [16] M.A. Islam, R.I. Chowdhury, N. Chakraborty, W. Bari, A multistage model for maternal morbidity during antenatal, delivery and postpartum periods, Stat. Med. 23 (2004) 137—158. [17] N. Chakraborty, M.A. Islam, R.I. Chowdhury, W. Bari, Utilisation of postnatal care in Bangladesh: evidence from a longitudinal study, Health Soc. Commun. 10 (6) (2002) 492— 502. [18] N. Chakraborty, M.A. Islam, R.I. Chowdhury, W. Bari, Analysis of ante-partum maternal morbidity in rural Bangladesh, Aust. J. Rural Health 11 (2003) 22—27. [19] N. Chakraborty, M.A. Islam, R.I. Chowdhury, W. Bari, H.H. Akhter, Determinants of the use of maternal health services in rural Bangladesh, Health Promot. Int. 18 (4) (2003) 327—337. [20] R.I. Chowdhury, M.A. Islam, H.H. Akhter, N. Chakraborty, Analysis of postpartum complications in relation to selected delivery characteristics in rural Bangladesh, Health and Population in Developing Countries, 2004, URL: http://www.jhpdc.unc.edu/. [21] W. Bari, R.I. Chowdhury, M.A. Islam, N. Chakraborty, The differentials and determinants of perinatal mortality in rural Bangladesh, Eur. J. Contracept. Reproduct. Health Care 7 (2002) 216—222.
A computer program to estimate the parameters of covariate dependent higher order Markov model Rafiqul Islam Chowdhurya,∗, M. Ataharul Islamb, Makhdoom Ali Shahc, Naser Al-Enezid a
Department of Health Information Administration, Kuwait University, P.O. Box 31470, Sulaibekhat 90805, Kuwait b Department of Statistics and OR, College of Science, King Saud University, Saudi Arabia c Department of Health Information Administration, Kuwait University, Kuwait d Department of Health Information Administration, Kuwait University, Kuwait Received 10 October 2003 ; received in revised form 21 August 2004; accepted 5 October 2004 KEYWORDS Longitudinal data; Covariate dependent; Higher order; Markov model; S-plus; Computer program
Summary This paper presents a computer program developed in S-plus to estimate the parameters of covariate dependent higher order Markov Chain and related tests. The program can be applied for two states Markov Chain with any order and any number of covariates depending on the PC capabilities. The program provides the maximum likelihood estimates of the parameters, together with their estimated standard error, t-value and significance level. It also produces the test results for likelihood ratio and model chi-square. To illustrate the program we have used a longitudinal data set on maternal morbidity of rural women in Bangladesh. The occurrences of haemorrhage, convulsion, or fits at different follow-ups were used as outcome variable. Economic status, wanted pregnancy, ages at marriage, and education of women were used as covariates. © 2004 Elsevier Ireland Ltd. All rights reserved.
1. Introduction Different statistical models are available to analyze the longitudinal data emerging from biological, epidemiological, demographic, sociological, econometric, and many other fields of study. Markov chain * Corresponding author. Tel.: +965 9503456; fax: +965 4830937. E-mail address: rafi[email protected] (R.I. Chowdhury).
model is one of the statistical techniques that play an important role to analyze such types of data [1—3]. In recent years, considerable interest has been witnessed in the development of multivariate models based on Markov chain and these models are employed to analyze the data generated from meteorology, reliability, survival analysis, etc. Muenz and Rubinstein introduced a first-order discrete time Markov chain for expressing the transition probabilities in terms of function of covari-
0169-2607/$ — see front matter © 2004 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2004.10.003
176
R.I. Chowdhury et al.
ates for a binary sequence of presence or absence of disease and suggested that this model can be extended for higher order Markov Chain Model [4]. They employed logistic regression models to analyze the transition probabilities from one state to another. The estimation of first-order Markov model is quite straightforward [5]. As the classical Markov chain model is a first order that is a process without memory, there are growing interests for higher order Markov models [6—9]. Islam et al. generalized the Muenz and Rubinstein models for covariate dependent higher order Markov model for binary outcomes [10]. They also discussed issues related to extension to more than two intercommunicating states, heterogeneity, and appropriate estimation and test procedures. Available statistical software do not always include procedure to estimate any newly developed model, which prevents researchers to apply the models for complexity in computational procedures. However, there are some softwares which have options to estimate transition matrix and some programs allow to estimate the Markov chain models with different specifications [11—13]. This paper describes a computer program designed to estimate the parameters of covariate dependent higher order Markov model for binary outcomes proposed by Islam et al., and the related test procedures [10]. The program was written using S-plus for windows and a follow-up data set on maternal morbidity of pregnant women was used to demonstrate the program.
2.1. Covariate dependent higher order model Following is the covariate dependent higher order Markov model proposed by Islam et al. [10] extending the Muenz and Rubinstein [4] model for first-order Markov chain. To illustrate the extension of the model for higher order, let us consider a second-order Markov model. The second-order Markov model for time points tj−2 , tj−1 , and tj with corresponding outcomes Yj−2 , Yj−1 , and Yj , respectively, is shown below: Y j−2
π01 (Y j = 1/Y j−1 = 0, X) =
eβ0 X
(1)
1 + eβ0 X
and
π11 (Y j = 1/Y j−1 = 1, X) =
eβ1 X 1+e
β1 X
.
(2)
where Xi = [1, X i1 , . . . , X ip ] is the vector of covariates for the ith person; 0 = [β00 , β01 , . . . , β0p ] the vector of parameters for the transition from 0; and 1 = [β10 , β11 , . . . , β1p ] is the vector of parameters for the transition from 1.
|
Yj
0
1
0
0
π000
π001
0
1
π010
π011
1
0
π100
π101
1
1
π110
π111
.
Using the definitions of vectors of covariates and parameters mentioned in the previous section, we can define the transition probabilities:
πuvw (Y j = 1/Y j−2 = 0, Y j−1 = 0, X) =
eβ001 X
1 + eβ001 X
u = 0, 1; v = 0, 1; w = 0, 1.
, (3)
It may be noted here that π000 + π001 = 1, π010 + π011 = 1, π100 + π101 = 1, and π110 + π111 = 1. Then the likelihood function, as a generalization can be expressed as follows:
2. Computational methods and theory The transition probabilities for the first order can be defined in terms of function of the covariates as follows:
Y j−1
L=
ni n
[{π000 }δ000ij {π001 }δ001ij ][{π010 }δ010ij
i=1 j=1
× {π011 }δ011ij ][{π100 }δ100ij {π101 }δ1011ij ][{π110 }δ110ij × {π111 }δ111ij ],
(4)
where L=L1 L 2 L3 L 4 . The first derivatives with respect to parameters are n
n
i ∂ ln L 1 = ∂βuvwq
i=1 j=1
× X iq δ000ij − (δ000ij + δ001ij )
eβ001 X
1 + eβ001 X
(5)
u = 0,1; v = 0,1; w = 0,1; q = 0, 1, 2, . . ., p. Similarly, we can differentiate ln L2 , ln L3 , ln L 4 , with respect to the corresponding parameters.
A computer program to estimate the parameters of covariate dependent Markov model We can solve for the sets of parameters by equating first derivatives to zero. Each set consists of (p + 1) parameters and hence the total number of parameters to be estimated here is 4(p + 1). The second derivatives for L1 are: ∂2 ln L 1 ∂βuvwq ∂βuvwl =−
ni n
vector: β = [β1 , β2 , . . . , β2k ] = [β , β , . . . , β k where βm mp ], m = 1, 2, . . ., 2 . m0 m1 To test the null hypothesis H0 :β = 0, we can employ the usual likelihood ratio test
−2[ln L(β0 ) − ln L(β)] ≈ χ22k p .
[X iq X il {(δuv0ij + δuv1ij )πuv0 (X)πuv1 (X)}],
i=1 j=1
(6)
where u = 0, 1; v = 0, 1; w = 0, 1. We can differentiate similarly for ln L2 , ln L3 , and ln L4 . Inverse of the (−1) (second derivative) provide estimates of the variance—covariance for the respective estimates of the parameters. Then the likelihood function for a Markov chain of order k, as a generalization, can be expressed as follows:
To test the significance of the qth parameter of the mth set of parameters, the null hypothesis is H0 :βmq = 0 and the corresponding Wald test is W=
ˆmq β . ˆmq ) se(β
Another alternative test procedure is illustrated here on the basis of the suggestions made by Billingsley [14] and Raftery and Tavarey [6]. The test statistic for the k-order Markov model is defined as:
k
L=
ni 2 n
[{πm }δmij {1 − πm }1−δmij ].
χ2 =
The log-likelihood function can be shown as ln L = ln L 1 + ln L 2 + · · · + ln L m + · · · + ln L 2k .
(8)
Here, Lm =
{n(j − k, j − k + 1, . . . j)
(7)
i=1 j=1 m=1
ni n
177
[{πm }δmij {1 − πm }1−δmij ]
(9)
i=1 j=1
The first derivatives with respect to parameters are ni n ∂ ln L m eβm X i X iq (1 − δmij ) − , = ∂βmq 1 + eβm X i i=1 j=1
m = 1, 2, . . . , 2k
and
q = 0, 1, 2, . . . , p.
(10)
The second derivatives for ln Lm are: n
n
i ∂2 ln L m =− [X iq X il {πm (1 − πm )}]. ∂βmq ∂βml
(11)
i=1 j=1
2.2. Testing for the significance of parameters The tests are discussed in Islam et al. [10]. Let us give an overview of the tests employed for testing the significance of the parameters. The vectors of 2k sets of parameters for the kth order Markov model can be represented by the following
− e(j − k, j − k + 1, . . . , j)}2 e(j − k, j − k + 1, . . . , j)
the sum is over all n(j − k, j − k + 1, . . . , j) for n(j − k, j − k + 1, . . . , +). The estimated frequencies can be obtained from the following relationship: e(j − k, j − k + 1, . . . , j) = n(j − k, j − k + 1, . . . , +)ˆ π(j − k, j − k + 1, . . . , j)
2.3. Fitting algorithm The following is an algorithm for fitting higher order Markov model from Eq. (7) using the Gauss—Newton approach: (1) Input data file, number of variables, order of Markov chain, and maximum number of iteration. (2) Read data file, initialize variables, calculate sets of models to estimate (e.g., for second order there are 22 = 4 sets of model, for third order 23 = 8 set of models and so on). (3) Compute the combination of outcomes needed to check different types of transitions for different sets of model. For example, for the second order we have to check the following combinations for outcomes for four sets of model:
178
R.I. Chowdhury et al. Y j−2
(4) (5) (6) (7) (8) (9)
Y j−1
0
0
| 0 π000
Yj 1
0
1
π010
1
0
π100
π011 π101
1
1
π110
π111
(11) Estimate  using the following: ˆj = β ˆj−1 − I −1 (β ˆj−1 )U(β ˆj−1 ) β
π001
(12) Update β with the latest value. (13) Calculate test statistics. (14) Iterate until convergence.
We need first row of the above matrix to check the ith individual for δ000ij = 1 if a transition type 0 0 0 is observed during jth follow-up for the ith individual, δ000ij = 0, otherwise; second row will need to check δ011ij = 1 if a transition type 0 1 1 is observed during jth follow-up for the ith individual, δ011ij = 0, otherwise; third row will need to check δ101ij = 1 if a transition type 1 0 1 is observed during jth follow-up for the ith individual, δ101ij = 0, otherwise; and last row will need to check δ111ij = 1 if a transition type 1 1 1 is observed during jth follow-up for the ith individual, δ111ij = 0, otherwise. Start the loop for iteration. Set initial value for parameter to 0 (β = 0). Start the loop for sets of models. Check the transition types for each set of model defined in (3). Compute the first derivatives and second derivatives for each set of model and store separately for each model. To get the score vector, combine the first derivative for all the models in a single vector. For example, if we have three covariates and a second-order Markov Chain then we will have a score vector of 16 elements constant + 3 covariates = 4 and 4 sets of model). If U1 , U2 , U3 , and U4 are the score vectors for four sets of model each with four elements, then we will get the following combined score vector: U(β) = [U 1 U 2 U 3 U 4 ]1×16
(10) Combine the second derivative off all the models as diagonal terms of the matrix, which will be the information matrix (for example, if we have 3 covariates and a second-order Markov Chain, the information matrix will be 16 × 16. Constant + 3 covariates = 4 and 4 sets of model). If I1 , I2 , I3 , and I4 are the 4 × 4 information matrix for four sets of model then we will get the following 16 × 16-combined information matrix: I1 0 0 0 0 I 0 0 2 I(β) = 0 0 I3 0 0
0
0
I4
16×16
3. Program description The function (markov.gen) can be run under S-plus command environment. It can be used for equal and unequal number of follow-ups. The function is written in an ASCII text file and could be copied and pasted in the S-plus command prompt to run [15]. To use the function five arguments are needed. Arguments are data file name including path (alphanumeric), total columns (numeric) in the data file, order of the Markov Chain (numeric, e.g., for first-order 1, for second-order 2, and so on), form of the transition probabilities (numeric, currently it is 1 for logistic form), and maximum number of iteration (numeric). The structure of the data file is ASCII free format. The first column in the data file should be case identification number; the second column is the follow-up number for each case. The third column is the outcome variable and the fourth and subsequent columns are the covariates. In addition, the first row in the data file should contain the entire variable name. Any software (e.g., SPSS, SAS, S-plus, etc.) can be used to prepare the data file in the required format to use our function for Markov Chain. To test the accuracy of the estimates of our program we used a hypothetical data set. The data set contains three observations with 12, 8, and 5 follow-ups and three covariates. Using the hypothetical data set the necessary calculations to estimate the parameter of our model were done manually. These results were compared with the results from the sample run of our S-plus program.
4. Sample run The application is based on the data from the Follow-up Survey on Maternal Morbidity in Bangladesh, conducted by the Bangladesh Institute of Research for Promotion of Essential and Reproductive Health and Technologies (BIRPERTH). The data were collected from November 1992 to December 1993. A total of 1020 pregnant women were included in the study (pregnancy less than 6 months). These subjects were followed-up with
A computer program to estimate the parameters of covariate dependent Markov model
Table 1
Distribution of follow-ups by occurrence of morbidity
Follow-up number
Occurrence of morbidities No n
1 2 3 4 5 6 7 8 Total
565 623 534 418 245 101 23 — 2509
Total Yes
% 57.0 67.9 69.3 70.4 66.2 68.2 67.6 — —
an interval on an average of 1 month, through full-term pregnancy, delivery and till 90 days postpartum period or 90 days after any other pregnancy outcome. The information on socio-economic background, pregnancy related care and practice, extent of morbidity during the index pregnancy, delivery and postpartum period, or abortions were collected. Out of 1020 women, 993 had at least one antenatal follow-up, and 1005 had information on pregnancy termination. Finally, 1006 had at least one postpartum follow-up. However, for our application, we used data for the antenatal follow-ups only. Table 1 shows the occurrence of morbidity at different follow-ups. The data comprise a total of 3827 follow-ups. There were maximum eight follow-ups, which was only for one woman. Dependent variable: Y it is defined as follows: Y it
179
= 1 if there was occurrence of haemorrhage, or convulsion or fits; = 0, otherwise
n
%
427 294 237 176 125 47 11 1
43.0 32.1 30.7 29.6 33.8 31.8 32.4 100.0
992 917 771 594 370 148 34 1
1318
—
3827
Table 2
Transition matrix
Transition type
→0
→1
0→0→0 1→0→0 0→1→0 0→0→1 1→1→0 1→0→1 0→1→1 1→1→1
459 91 31 49 40 26 8 44
72 30 16 27 28 20 23 132
On the S-plus command prompt the following command was used to invoke our function to estimate the parameters of third order Markov Chain with logistic link function: • markov.gen(‘‘c:/markov/marapp.dat’’,7,3,1,12) The model was converged after five iterations. The program produced the transition matrix and estimate of the parameter for each set of equations, which is presented in Table 2. The sets of equations follow the same order as the sequence of different transition types in transition matrix. It also computes the model chi-square, likelihood ratio, AIC, BIC, and Billingsley statistic, which are presented in Table 3.
the covariates were categorized and used as follows: economic status: good = 1, bad = 0; wanted pregnancy: yes = 1, and no = 0; education: no education = 0, primary or higher education = 1; and age at marriage: less than or equal to 15 years = 0, more than 15 years = 1.
5. Hardware and software specifications
We have already conducted several studies on maternal morbidity in Bangladesh using the same set of data for examining different aspects of maternal morbidity. Some of the relevant papers are Islam et al., Chakraborty et al., Chowdhury et al., and Bari et al. [16—21].
The function for Markov Chain is written in S-plus 6.0 and has approximately 325 lines of code. The program requires intensive computer resources, mainly to evaluate the likelihood function. The contribution to the likelihood function of two outcomes, in each case depending on the order of the
180
Table 3
R.I. Chowdhury et al.
Estimates of the parameters of third order Markov Chain
Characteristics and different transition types
Estimate
Standard error
t-value
p-value
0 → 0 → 0 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
−1.6804 0.2478 −0.1146 −0.1914 −0.1528
0.2548 0.3369 0.2774 0.2760 0.2602
−6.5945 0.7354 −0.4133 −0.6933 −0.5872
0.0000 0.4621 0.6794 0.4881 0.5570
1 → 0 → 0 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
−0.6063 −0.2463 −0.8580 0.6512 −0.5059
0.3613 0.6402 0.4387 0.4844 0.4598
−1.6781 −0.3847 −1.9556 1.3444 −1.1003
0.0933 0.7005 0.0505 0.1788 0.2712
0 → 1 → 0 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
−0.8476 −1.1475 0.2453 1.0430 −0.3283
0.5932 0.9642 0.6454 0.7268 0.6915
−1.4289 −1.1901 0.3801 1.4350 −0.4747
0.1530 0.2340 0.7038 0.1513 0.6350
0 → 0 → 1 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
0.2030 −0.0303 −0.6704 −0.5273 −0.5190
0.4478 0.6957 0.5028 0.5987 0.5298
0.4534 −0.0435 −1.3334 −0.8808 −0.9796
0.6503 0.9653 0.1824 0.3784 0.3273
1 → 1 → 0 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
−0.2380 1.2853 −0.4223 −0.1862 −0.2244
0.4623 0.6418 0.5361 0.5513 0.5527
−0.5147 2.0026 −0.7876 −0.3378 −0.4061
0.6068 0.0452 0.4309 0.7355 0.6847
1 → 0 → 1 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
0.0544 1.5457 −0.1858 −0.0901 −1.0149
0.5271 0.9456 0.6468 0.6757 0.7107
0.1032 1.6346 −0.2873 −0.1334 −1.4280
0.9178 0.1021 0.7738 0.8939 0.1533
0 → 1 → 1 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
1.4462 −0.1631 −0.8268 1.0574 0.5812
1.1843 1.0135 1.2022 1.1807 1.2298
1.2211 −0.1609 −0.6878 0.8956 0.4726
0.2221 0.8722 0.4916 0.3705 0.6365
1 → 1 → 1 → 0 or 1 Constant Economic status Wanted pregnancy Age at marriage Education of women
1.4786 0.2558 −0.5342 −0.4100 0.1508
0.3299 0.4584 0.3743 0.3734 0.4103
4.4820 0.5581 −1.4274 −1.0980 0.3675
0.0000 0.5768 0.1535 0.2722 0.7133
Model chi-square (40 d.f.) Likelihood ratio AIC BIC BILLINGSLEY
403.08 445.05 1154.32 1350.35 1410.95
0.000 0.000 0.000 0.000 0.000
A computer program to estimate the parameters of covariate dependent Markov model Markov Chain, requires the evaluation of the transition probability matrix for a given set of covariates. The CPU time needed to run this function depends on the order of the Markov Chain, number of observations, the number of covariates, and the speed of PC processor and the amount of PC memory.
6. Mode of availability The full program ‘kern0pt‘markov.gen’’ is available on request. Any one, who is interested, can request the corresponding author for the complete program. The program file and the necessary instruction for users will be sent as e-mail attachment.
7. Conclusions In this paper we explained the development of S-plus program for application to real life data, which contains necessary instructions to use our S-plus function to estimate the parameter of covariate dependent higher order Markov Chain, using a follow-up data set. Despite the computational complexities, with the use of current state of high quality software like S-plus (StatSci Division, MathSoft Inc.), the estimates of parameters pose no problem. Work is in progress to refine and add new link functions and some other tests. We hope other researchers who are interested in using Markov Chain will use the function developed by us to estimate the parameter of the model.
Acknowledgements We gratefully acknowledge the permission of Dr. Halida Hanum Akhter, Director, BIRPERHT, for using the data for the application in this paper. The authors would like to thank Mahbub E. Elahi, K. Chowdhury, and Arindom Sen of the Bangladesh Institute of Research for Promotion of Essential and Reproductive Health and Technologies (BIRPERHT) for their assistance during different phases of this work. The authors are greatly indebted to the Ford Foundation for funding the data collection of the maternal morbidity study.
References [1] P.S. Albert, A Markov model for sequence of ordinal data from a relapsing-remitting disease, Biometrics 50 (1994) 51—60.
181
[2] P.S. Albert, A.W. Myron, A two state Markov Chain for heterogeneous transitional data: a quasilikelihood approach, Stat. Med. 17 (1998) 1481—1493. [3] R. Kay, A Markov model for analyzing cancer markers and disease states in survival studies, Biometrics 42 (1986) 855—865. [4] L.R. Muenz, L.V. Rubinstein, Markov models for covariate dependence of binary sequences, Biometrics 41 (1985) 91—101. [5] N. Harun, Higher order Markov Chain model for analyzing covariate dependence of longitudinal data. A Thesis Submitted in Partial Fulfillment of the Degree of M.Sc., Department of Statistics, University of Dhaka, 2002. [6] A.E. Raftery, S. Tavare, Estimating and modeling repeated patterns in higher order Markov Chains with mixture transition distribution, Appl. Stat. 43 (1) (1994) 179— 199. [7] A.E. Raftery, A model for higher order Markov Chains, J. R. Stat. Soc. 47 (1985) 1528—1539. [8] F.M. Bass, M.M. Givon, M.U. Kalwani, D. Reibstein, G.P. Wright, AN investigation into the order of the brand choice process, Marketing Sci. 3 (1984) 267—287. [9] J.A. Logan, A structural model of the higher-order Markov process incorporating reversion effects, J. Math. Sociol. 8 (1981) 75—89. [10] M.A. Islam, R.I. Chowdhury, A. Baharum, A higher Markov Model for analyzing covariate dependence, in: Proceedings of the National Conference of Management Science and Operations Research, Melaka, Malaysia, 24—25 June, 2003, pp. 25—33. [11] SAS Institute Inc., SAS/STAT User’s Guide, version 6, vol. 2, fourth edition, SAS Institute Inc., Cary, NC, 1989. [12] MARCH 1.1 User Guide, Markov chains Computation for Homogeneous and Non-Homogeneous data. Universite de Lausanne, Switzerland, 2001. [13] G. Marshall, W. Guo, R.H. Jones, MARKOV: a computer program for multi-state Markov models with covariables, Comput. Methods Program Biomed. 47 (1995) 147—156. [14] P. Billingsley, Statistical Inference for Markov Process, The University of Chicago Press, Chicago, USA, 1961. [15] S-plus 6 for Windows, Programmers guide Insightful Corporation, Seattle, Washington, USA, 2001. [16] M.A. Islam, R.I. Chowdhury, N. Chakraborty, W. Bari, A multistage model for maternal morbidity during antenatal, delivery and postpartum periods, Stat. Med. 23 (2004) 137—158. [17] N. Chakraborty, M.A. Islam, R.I. Chowdhury, W. Bari, Utilisation of postnatal care in Bangladesh: evidence from a longitudinal study, Health Soc. Commun. 10 (6) (2002) 492— 502. [18] N. Chakraborty, M.A. Islam, R.I. Chowdhury, W. Bari, Analysis of ante-partum maternal morbidity in rural Bangladesh, Aust. J. Rural Health 11 (2003) 22—27. [19] N. Chakraborty, M.A. Islam, R.I. Chowdhury, W. Bari, H.H. Akhter, Determinants of the use of maternal health services in rural Bangladesh, Health Promot. Int. 18 (4) (2003) 327—337. [20] R.I. Chowdhury, M.A. Islam, H.H. Akhter, N. Chakraborty, Analysis of postpartum complications in relation to selected delivery characteristics in rural Bangladesh, Health and Population in Developing Countries, 2004, URL: http://www.jhpdc.unc.edu/. [21] W. Bari, R.I. Chowdhury, M.A. Islam, N. Chakraborty, The differentials and determinants of perinatal mortality in rural Bangladesh, Eur. J. Contracept. Reproduct. Health Care 7 (2002) 216—222.