ASSOCIATION BETWEEN EMERGENCY DEPARTMENT VISITS DURING PREGNANCY AND ENHANCED PRENATAL HEALTH CARE PROGRAM PARTICIPATION – AN APPLICATION OF FIXED-EFFECTS COUNT MODELS By Zhiheng Chen A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Biostatistics – Master of Science 2017 ABSTRACT ASSOCIATION BETWEEN EMERGENCY DEPARTMENT VISITS DURING PREGNANCY AND ENHANCED PRENATAL HEALTH CARE PROGRAM PARTICIPATION – AN APPLICATION OF FIXED-EFFECTS COUNT MODELS By Zhiheng Chen We study the number of Emergency Department (ED) visits during pregnancy for low income women in one county in southwest Michigan. The outcome distribution has a large proportion of zeros. This study examines the association between ED visits during pregnancy and the enhanced prenatal health care programs (Maternal Infant Health Program and Strong Beginning program). These enhanced prenatal health care programs enrolled high-risk women who might use ED more often than those with low-risk pregnancy. Due to selfselection and potential measurement errors in measuring the enrollment in enhanced prenatal health programs, fixed-effects models were used to fit the data. Results show that having enhanced prenatal health care programs helps reducing the expected number of ED visits during pregnancy. Keywords: Count Data; Endogenous variables; Fixed-effects models; Enhanced Prenatal Health Care Programs; Emergency Department TABLE OF CONTENTS LIST OF TABLES …………………………………………………………………………..iv LIST OF FIGURES ………………………………………………………………………….v Chapter 1 Introduction ……………………………………………………………………...1 Chapter 2 Background ………………………………………………………………………2 2.1 Count Data Models ………………………………………………………………….......2 2.1.1 Poisson Model ……………………………………………………………………….2 2.1.2 Negative Binomial Model …………………………………………………………….2 2.1.3 Hurdle Model ………………………………………………………………………..3 2.1.4 Zero-Inflated Model ………………………………………………………………….4 2.2 Causal Inference Methods ………………………………………………………..……...5 2.2.1 Causal Model ………………………………………………………………….……..5 2.2.2 Identification of Causal Effects ………………………………………………………..5 Chapter 3 Model and Data ………………………………………………….....…………….7 3.1 Model Specification ……………………………………………………………...…........7 3.1.1 Fixed-Effects Models …………………………………………………………………7 3.1.2 Fixed-Effects Count Models …………………………………………………………..7 3.1.3 Estimations …………………………………………………………………..………9 3.2 Data Collection ………………………………………………………………………...10 3.2.1 Sample …………………………………………………………………………......10 3.2.2 Missing Value ……………………………………………………………………….11 3.2.3 Outcome ……………………………………………………………………………11 3.2.4 Exposure …………………………………………………………………………...14 3.2.5 Covariates ………………………………………………………………………….14 A - Person Level ……………………………………………………………………….14 B - Census Tract Level ………………………………………………………………….14 Chapter 4 Results …………………………………………………………………………...18 Chapter 5 Discussion ……………………………………………………………………….24 APPENDICES ……………………………………………………………………………...25 APPENDIX A ………………………………………………………………………….......26 APPENDIX B ……………………………………………………………………………...27 APPENDIX C ……………………………………………………………………………...28 APPENDIX D ……………………………………………………………………………...29 REFERENCES ……………………………………………………………………………..30 iii LIST OF TABLES Table 1. Risk characteristics of full Medicaid-insured African American women with singleton birth in Kent County, MI 2009-2015………………………………………………13 Table 2. Distribution of the outcomes …………………………………………………...........15 Table 3. Test for non-nested models ………………………………………………………….16 Table 4. Results from Poisson model …………………………………………………............17 Table 5. Results from Poisson, Negative Binomial and ZIP model ………………………....19 Table 6. Results from FE-NB, FE-Hurdle and FE-ZIP model ………………………..............20 Table 7. Results for Fixed-Effect Hurdle model ………………………………………...........22 Table A1. Number of Missing Values ………………………………………………………28 iv LIST OF FIGURES Figure 1. Flow Chart…………………………………………………………………………26 Figure 2a. Distribution of the overall outcome in each treatment group. .…………………….27 Figure 2b. Distribution of the outcome (NE, PCT and PA using NYU algorithm) in each treatment group. ………………………………………………………………………….…..27 v Chapter 1 Introduction Adverse perinatal outcomes such as preterm birth and low birth weight represent a substantial clinical and public health problem. Getting early and regular prenatal care improves the chance of a healthy pregnancy. Many pregnant women continue to use the emergency department (ED) possibly due to lack of appropriate prenatal care. If unnecessary ED visits can be prevented, significant resources can be redirected toward more urgent needs. Maternal Infant Health Program (MIHP) and Strong Beginnings (SB) are two enhanced prenatal health care programs for pregnant women and infants in Michigan (SB is based in Kent county, Michigan). MIHP is Michigan’s largest home visiting program for Medicaideligible pregnant women and infants. It provides home visitation supports and care coordination to supplement regular prenatal care and promote healthy pregnancies and positive birth outcomes. The SB program is focused on high-risk pregnant African American women and infants in the Greater Grand Rapids area. It provides intensive outreach and case management in addition to MIHP. Our research question is whether there is an association between the ED visits during pregnancy and participation in different types of enhanced prenatal health care programs. We hypothesize that women in MIHP and SB will have fewer ED visits during pregnancy compared to similar women in neither program. Because of the potential measurement error in measuring the enrollment in MIHP and self-selection, which will be discussed in Chapter 3, the results from usual count models might be biased. For this reason we will use fixed-effects model to fit our data. 1 Chapter 2 Background 2.1 Count Data Models 2.1.1 Poisson Model As the outcome of this study, the number of ED visits is a count variable, the benchmark model is the Poisson model. Suppose there is a random variable mean , > 0. The density of is ( ) , ! ( = )= and ( )= under Poisson distribution with = 0,1,2, … ( )= . Similar to the linear regression models, Poisson regression models connect = (1, , ,…, with covariates, ) with unit as its first element. The Poisson model takes the form log ( | ) = ! 2.1.2 Negative Binomial Model There are 13 different derivations for the negative binomial (NB) distribution identified by Patil and Boswell in 1970 (Patil and Boswell, 1970). The most often cited one is the mixture of a Poisson and a gamma distribution. Suppose the gramma distribution has mean 1 and variance ", then the unconditional NB distribution is #$%( , "), with density Γ( + " ) " ( = )= ( ) Γ( + 1)Γ(" ) " + * +, where " > 0. (1 − " " + ) , = 0,1,2, … " is called the overdispersion parameter and Γ() is the gamma function: 4 Γ( ) = / 0 1 5 2 30 . One of the advantages of using a NB model instead of a Poisson model is that NB models can 2 deal with overdispersion. If ~#$%( , "), then ( )= ( ) = (1 + " ). Because " is always greater than 0, the NB model is not able to model underdispersion. This is why researchers do not simply fit count data using a NB model instead of a Poisson model (Hilbe, 2013). Hilbe (2013) suggests researchers to fit the data using Poisson regression first, and if apparent overdispersion is found, then use NB models instead. According to Hardin and Hilbe’s definition, apparent overdispersion – indicates the lack of fit, possibly due to: (a) the model omits important explanatory predictors; (b) the data include outliers; (c) the model fails to include interaction terms; (d) a predictor needs to be transformed to another scale; or (e) the assumed linear relationship between the link-function transformed response and predictors is misspecified. 2.1.3 Hurdle Model When the outcome contains a large proportion of zeros, the usual Poisson or NB models may not adequately fit the data. For example, in our study, there was a large proportion of pregnant women (39%) who had not gone to the ED during pregnancy. The decision to use ED or not might be governed by a different process from the process that determines how many times ED is used. A hurdle model (Cragg, 1971) can be used to deal with the different processes. Let the first process for using the ED at all be governed by density 8 ( ; : ) and the second process for the number of times using ED by density 8 ( ; : ). Then the density of a hurdle model is ( = )=; 8 (0; : ), 1 − 8 (0; : ) 3 8 ( ;: ) , 1 − 8 (0; : ) <8 <8 =0 >0 Specifically, take a logit-Poisson hurdle model as an example, suppose =5 = ( = 0) is modeled by a logistic regression, and is the mean of the Poisson part, then the density of the model is ( = )=; (1 − = 5 ) =5 , > ! (1 − >) <8 , <8 =0 >0 2.1.4 Zero-inflated Model Another way to fix the problem of excess zeros is the zero-inflated models. The zero-inflated models (Lambert, 1992) were introduced to allow the count data with a large proportion of zeros. Suppose there is a base count density 8 ( ), but it doesn’t predict too many zeros. One solution is to add another component that inflated the probability of zero. Then the density function for the zero-inflated model is ( = | )=? = + (1 − =)8 (0), (1 − =)8 ( ), <8 = 0 <8 > 0 where = is the proportion of the part that generate structural zeros. It can be a constant or a random probability depending on covariates via a binary outcome model such as logit. Specifically, for a zero-inflated Poisson (ZIP) model, the base distribution is a Poisson distribution with mean (Mullahy et al, 1986). The density is = + (1 − =) > , > ( = | )=; (1 − =) , ! The mean of ZIP is: and its variance is: <8 <8 =0 >0 ( | ) = (1 − =) , ( | ) = (1 − =)( + = ). Moreover, for a zero-inflated Negative Binomial (ZINB) variable , the base distribution is a #$( , "), so the ZINB density is 4 ( = | )=@ = + (1 − =)(1 + ") (1 − =)ℎ( , , "), where ℎ( , , ") is the NB density. The mean of ZINB is: * +, , <8 = 0 <8 > 0 ( | ) = (1 − = ) , and its variance is: ( | ) = (1 − =)B (1 + " ) + = C. 2.2 Causal Inference Methods 2.2.1 Causal Model A causal model contains BD, , C , where D is a set of exogenous variables, endogenous variables and is a set of is a set of structural equations (Pearl et al, 2000). Exogenous variables are independent variables that have some effects on a model but are determined by some factors outside the model. Instead, an endogenous variable is dependent on other variables in the model. 2.2.2 Identification of Causal Effects For a randomize control trial (RCT), the assignment of treatment is random, thus exogenous. If the experiment is conducted properly, the effect of the treatment can be directly estimated by comparing the outcome between the treated and control groups. However, for observation studies, take this study as an example, the participation in either MIHP or SB programs is voluntary, likely determined by many factors. Women in different programs may be systematically different in terms of risks for poor birth outcomes. One way to solve this problem is matching using propensity scores (PS). Although PS matching is not the same as randomization, it can balance observed baseline characteristics for study population under different exposures (Mnatzaganian et al, 2015). However, if there are some endogenous variables in the model, the PS matching method is no 5 longer valid. In this study, the enrollment of MIHP cannot be determined exactly by all the information we have. Measurement errors and self-selection might exist. Usually instrumental variables method is used to deal with these problems. A fixed-effects model can also be used. Despite the great usefulness of these methods, there are lots of assumptions for these methods. As Imbens and Rubin stated: “The potential outcomes for any unit do not vary with the treatments assigned to other units, and, for each unit, there are no different forms or versions of each treatment level, which lead to different potential outcomes” (Imbens & Rubin, 2015). We describe the fixed-effects models we use to mitigate the endogeneity problem with count data. 6 Chapter 3 Model and Data 3.1 Model Specification 3.1.1 Fixed-Effects Models Fixed-effects models are often used with panel data (Grafarend et al, 2006). Unlike the random- effects models, the distribution of “effect” parameter EF in fixed-effects models may have any distributions. Whether EF are correlated with other covariates or not is the important difference between random-effects and fixed-effect models (Wooldridge, 2002; Cameron, 2013). Linear fixed-effects models have the following form: FG = FG ! + EF + HFG where the index < represents individuals and I represents the I0ℎ observation for individual <. The error terms are assumed to be independent with other covariates and EF . Fixed-effects models for count outcomes can be written as: FG J = ( FG = FG | FG , EF ) FG ! + EF where g is the log link. If the conditional distribution for FG given by FG and EF is assumed to be Poisson distribution, the coefficients can be estimated using conditional maximum likelihood estimator (CMLE) or MLE with dummy variables. In other cases, if the conditional distribution is assumed to be NB distribution, then only MLE with dummy variables can be used (Gardiner et al, 2009). 3.1.2 Fixed-Effects Count Models Suppose there are treatment groups and a control group, each individual < chooses a group from ( + 1) choices. Like the linear fixed-effects models, we have exogenous covariates with associated coefficients, and an “effect” term which is allowed to be depend on those exogenous 7 covariates. Suppose binary variables 35 , 3 , … , 3 represent the observed treatment choices, and we assume the probability of treatment has a mixed multinomial logit structure. Let EF = (EF , EF , … , EF ) denote the individual fixed-effects term and KF denotes exogenous covariates, then the density function is defined as 3FG = 1LKF , EF = for I = 1, … , . MN *O PQNO 1 + ∑2T MN *S PQNS For a fixed-effects count data, the conditional outcome equation for each subject < is logB ( F |3F , where F F , EF )C F! = +U VG 3FG + U GT EFG GT are the exogenous covariates with associated parameters !. For a fixed-effects hurdle model, we assume that the second process which determines the counts is governed by the Poisson distribution, so that we can write the density of the hurdle model 8( F |3F , F , EF ) = ; 1 − WF , WF >N F F ! (1 − where log( F ) = ( F |3F , F , EF ) and XYJ<0(WF ) = F! Z N >N ) <8 , <8 F F =0 >0 + ∑GT VGZ 3FG + ∑GT EFG . The joint distribution of treatment and outcome variables can be written as the product of these two marginal densities ( F , 3F | F , KF , EF ) = MN *O PQNO ^ (1 − WF ) \ 1 + ∑` ] \WF [ >N F ! (1 − F N aT >N ) MN *_ PQN_ , MN *O PQNO 1 + ∑`aT MN *_ PQN_ <8 , <8 =0 F F >0 Similar to the hurdle model, the fixed-effect Zero-inflated Poisson model has density 8 ( F , 3F | F , KF , EF ) = ^B=F + (1 − =F ) \ ] \ (1 − =F ) [ >N F! >N F N C MN *O PQNO 1 + ∑aT ` MN *_ PQN_ MN *O PQNO 1 + ∑aT ` MN *_ PQN_ , , <8 <8 F F =0 >0 3.1.3 Estimations Because the fixed-effects terms are unknown, which incorporate the unobserved characteristics common to individual <′c treatment choice and outcome (Deb et al, 2006), we cannot estimate the model using the maximum likelihood estimators. We turn to the simulated maximum likelihood method. A STATA command “mtreatnb” was used to fit the Fixed-effect Negative Binomial model (Deb et al, 2006). All estimators for the FE-Hurdle model and FE-ZIP model were computed with Stata version 14.2 on a 2.4-GHz PC running macOS Sierra. There is no existing command for FE-Hurdle model and FE-ZIP model, we wrote our own programs using Stata’s ml routine to run the regression. First we assume that EFG are independently and identically distributed draws from the standard normal distribution and J(EF ) is the joint distribution of EFG . J(EF ) will be integrated out of the joint conditional (on the fixed effects) density ( F , 3F | F , KF , EF ) to derive the joint unconditional (on the fixed effects) density ( F , 3F | F , KF ). For the FE-Hurdle model and the FE-ZIP model, ( F , 3F | F , KF ) = / W( F , 3F | F , KF , EF )J(EF )3 EF Using simulation-based estimation (Gouriéroux et al, 1996), we can get the simulated joint density by drawing EFG # times from the density J(EF ), let EdFe be the f_0ℎ draw. The simulated density can be write by 9 i 1 ( F , 3F | F , KF ) ≈ U W( F , 3F | F , KF , EdFe ) # eT In our case, we can write our simulated log-likelihood function for the FE-Hurdle models as ln i ∑i eT {B(1 − WF ) n ln ( F , 3F | F , KF ) ≈ m o pq l N O NO r n P∑st, l mN o_ pqN_ C( N T5) + BWF n m o pq l N O NO l +uN >N vN N !( l +uN ) P∑r l mnN o_ pqN_ st, C( N w5) }. The simulated log-likelihood function for FE-ZIP model is ln i ∑i eT {{B=F + (1 − =F ) >N C l r ln ( F , 3F | F , KF ) ≈ mnN oOpqNO P∑st, n l mN o_ pqN_ }( T5) + B(1 − =F ) l +uN >N vN N l r mnN oOpqNO n P∑st, l mN o_ pqN_ C( N w5) }. 3.2 Data Collection As mentioned above, many pregnant women still have been visiting the ED possibly due to lack of appropriate prenatal care. While MIHP and SB program aim to provide enhanced prenatal care to pregnant women, we would like to see if these programs can reduce the number of ED visits during pregnancy. 3.2.1 Sample African American women in Medicaid in Kent county with singleton birth without congenital anomalies are our study population. From the vital record data, there were totally 445,331 births from 2009 to 2015 in Michigan. Among them, 31,613 were in Kent County. There were 6,752 births by African American women, of which 6,171 were by women who had full Medicaid four months before birth and 5,920 were singleton. 12 births were excluded because of fetal death. Four were excluded because of congenital anomalies. Congenital anomalies were defined if a birth had any of the following records: anencephaly, meningomyelocele/spina bifida, congenital heart disease, cyanotic congenital heart disease, omphalocele, gastroschisis, limb reduction, cleft lip with or without palate, cleft palate alone, Down syndrome – karyotype 10 confirmed, Down syndrome – karyotype pending, suspected chromosomal disorder – confirmed, suspected chromosomal disorder – pending, and hypospadias. After that, 71 births were excluded because their gestation ages were either shorter than 20 weeks or greater than 46 weeks, and another seven births were excluded because there were no claims records for the mothers. In the end, 5,826 births were included in this study. 533 women enrolled in MIHP and SB, 2,490 enrolled in MIHP only and the other 2,803 enrolled in neither of the two programs. The detailed numbers of births in each study group were given in Appendix - Figure 1. 3.2.2 Missing Value There were a few missing values in the study data set. The number of missing values for each variable was shown in Appendix - Table A1. To deal with the problem of missing values, the most frequent category of the variable was imputed for missing values of categorical variables. For continuous variables, the missing values were imputed with the mean of the non-missing values of the variable. Appendix - Table A1 also shows the imputed values for each variable with missing values. After solving the problem of missing values, the characteristics of the study population were shown in Table 1. 3.2.3 Outcome The outcome of this study is the number of ED visits during pregnancy. An ED visit was defined by using administrative claims data. Out-patient claims with revenue codes between 0450 and 0459 were defined as ED claims. Each unique date of those ED claims was counted as a unique ED visit, and only those visits during pregnancy were used to generate the outcome variable. The number of ED visits during pregnancy was counted as an overall outcome in our study. The distribution of the total number of ED visits during pregnancy was shown in Appendix – Figure 2a. It is well known that not all ED visits are truly emergent. Ambulatory care sensitive conditions 11 (ACSCs) are “a set of medical conditions, such as asthma, diabetes, gastroenteritis, or hypertension, for which timely, appropriate primary care can prevent or reduce the likelihood of so-called ‘avoidable’ hospitalization or emergency room visits.” (Falik et al, 2001) The number of ED visits is the outcome of this study, and it is important to distinguish ED visits for ACSCs from the truly unavoidable ones. The Emergency Department Algorithm (EDA) is developed by the New York University (NYU) Center for Health and Public Service Research for identifying ACSCs (NYU ED website). The EDA describes the likelihood of an ED visit in each of the following four categories (Jones et al, 2013): “(1) Non-emergent (NE): The patient’s initial complaint, presenting symptoms, vital signs, medical history, and age indicated that immediate medical care was not required within 12 hours. (2) Emergent/Primary care treatable (EPCT). (3) Emergent/ED care needed/Preventable or Avoidable (EPCA). (4) Emergent/ED care needed/Not preventable or Avoidable (ENPA).” The key to this classification is the principal diagnosis code, which is identified by the ICD9 codes (change to ICD-10 codes after 10/1/2015). This algorithm was developed based on a sample of nearly 6,000 ED records in a general hospital in New York, including the initial complaint, presenting symptoms, medical history, vital signs, patient characteristics, diagnoses, procedures performed and resources used in ED. The researchers estimated the probabilities for each diagnosis code being in different categories based on these ED records. To determine the category of each diagnosis code, we can either using the one with highest probability or the one with the probability over 75%. Those cases under Emergent/ED care needed/Not Preventable were not supposed to be reduced by prenatal care programs. The first three types of ED visits were combined as another outcome in which we were more interested. The distribution of the number of the combination of the 12 first three types of ED visits was shown in Appendix – Figure 2b. Table 1. Risk characteristics of full Medicaid-insured African American women with singleton birth in Kent County, MI 2009-2015. (2) (3) (1) Characteristics SB MIHP not No MIHP Total (% or mean) (N=533) SB (N=2.803) (N=5,826) (N=2,490) Average mother’s age, n (%) b <18 39 (6) 141 (6) 118 (4) 298 (5) 18-24 364 (55) 1223 (52) 1344 (48) 2931 (50) 25-29 147 (22) 559 (24) 764 (27) 1470 (25) 30-34 69 (10) 308 (13) 383 (14) 760 (13) 35+ 44 (7) 138 (6) 194 (7) 376 (6) ab Marital status, n (%) Married 56 (8) 343 (14) 454 (16) 853 (15) Unmarried, paternity acknowledged 274 (41) 970 (41) 1190 (42) 2434 (42) Mother only on birth certificate 333 (50) 1056 (45) 1159 (41) 2548 (44) Education, n (%) b HS 176 (27) 668 (28) 1065 (38) 1909 (33) 4 (1) 42 (2) 49 (2) 95 (2) Hispanic ancestry, n (%) ab 484 (73) 1760 (74) 1886 (67) 4130 (71) WIC program, n (%) b 453 (68) 1512 (64) 1625 (58) 3590 (62) Full Medicaid before conception, n (%) ab 100 (15) 323 (14) 310 (11) 733 (13) Smoking, n (%) b 32 (5) 103 (4) 94 (3) 229 (4) Quit smoking, n (%) b 122 (18) 399 (17) 366 (13) 887 (15) Others in household smoked, n (%) 6 (1) 27 (1) 27 (1) 60 (1) Drinking, n (%) 8 (1) 28 (1) 26 (1) 62 (1) Diabetes prior to pregnancy, n (%) 31 (5) 88 (4) 60 (2) 179 (3) Hypertension prior to pregnancy, n (%) b 42 (6) 163 (7) 195 (7) 400 (7) Previous preterm birth, n (%) Rapid repeat pregnancies <18 months from prior birth to current 161 (24) 653 (28) 831 (30) 1645 (28) conception >= 18 months from prior birth to current 272 (41) 878 (37) 1080 (39) 2230 (38) conception No prior deliveries 192 (29) 697 (29) 741 (26) 1630 (28) Unknown 38 (6) 141 (6) 151 (5) 330 (6) 29.6 (8.3) 28.9 (8.2) 28.6 (7.7) 28.9 (8.0) Pre-pregnancy BMI, mean (SD) b Abbreviations: MIHP, Maternal Infant Health Program; SB, Strong Beginnings. a : statistical significant difference (p<0.05) when comparing (1) vs (2). b : statistical significant difference (p<0.05) when comparing (1) vs (3). For simplicity, we use the category with highest probability to determine the visit type, and let %Fy denote the overall ED visits for woman < during pregnancy and %Fz denotes the ED visits after using NYU classify algorithm for woman < during pregnancy (combine NE, PCT and PA together). 13 3.2.4 Exposure There were three groups of women from the cohorts: those enrolled in SB; those enrolled in MIHP only; those not enrolled in either program. We used the procedure codes in Medicaid claims data to define the enrollment of MIHP. In this study, a woman was defined to be in MIHP if she has used MIHP for professional visits at least once during pregnancy. The professional visits include preventive counseling, prenatal care risk assessment, comprehensive multidiscipline evaluation and program intake assessment. Except for the claims for professional visits, there are also other claims for MIHP such as Nonemergency transportation – taxi, but we did not count these claims. Since we did not have a specific indicator for the enrollment of MIHP, the way we defined MIHP could only give us an approximation, which may lead to some measurement errors in the exposure. 3.2.5 Covariates A – Person Level Administrative claims and vital records for African American women who gave births between 2009 and 2015 in Kent county, Michigan were used in this study. Fifteen personal-level covariates were included in this study: moms’ age group, marital status, educational level, Hispanic or not, having received WIC program or not, full Medicaid before conception or not, smoking, quit smoking or not, having others in household who smoked, drinking, diabetes prior to pregnancy, Hypertension prior to pregnancy, previous preterm birth, rapid repeat pregnancies and pre-pregnancy BMI. B – Census Tract Level Data from American Community Surveys (ACS) in these years were used to create Census tract level social-demographic characteristics. Two neighborhood deprivation indices were created in the UK and were used for the UK populations – Townsend Material Deprivation Index (Nagaraja, 2015) and Jarman Underprivileged Area Score (Nagaraja, 2015). The 14 Table 2. Distribution of the outcomes Outcomes – Overall NE mean (min-max) SB 2.17 (01.17 17) (0-11) MIHP, no SB 1.87 (01.00 35) (0-25) None 1.31 (00.69 30) (0-18) PCT PA NPA AIDP UNCLS 0.32 (05) 0.28 (06) 0.21 (05) 0.03 (0-2) 0.02 (0-3) 0.01 (0-3) 0.03 (0-1) 0.03 (0-3) 0.02 (0-2) 0.06 (0-2) 0.07 (0-7) 0.04 (0-3) 0.55 (0-6) 0.48 (0-7) 0.34 (0-11) Abbreviations: MIHP, Maternal Infant Health Program; SB, Strong Beginnings; NE, Non-emergent; PCT, Emergent/Primary Care Treatable; PA, Emergent/Preventable or Avoidable; NPA, Emergent/Non-Preventable or Non-Avoidable; AIDP, alcohol, injury, drug or mental health; UNCLS, Unclassified. Townsend index is created based on four variables: 1. % of households with no vehicle 2. % of households with more than one occupant per room (i.e., overcrowding) 3. % of dwellings renter-occupied (i.e., housing tenure) 4. % of people above 16 years who are unemployed Jarman index is created based on eight variables: 1. % elderly living alone 2. % under 5 3. % persons living in a single parent household 4. % unemployment, above 16 5. % persons living in a household with more than 1 person per room 6. % moved within the last year 7. % born overseas from a non-English speaking county 8. % in social class 5, unskilled (UK) In 2015, these two indices were adapted to the US populations and fit for Census tract level data from ACS in New York and Bronx Counties by Nagaraja (Nagaraja, 2015). He found that both indices can identify deprived regions. However, the validity of these indices is still questionable in the wider US context (Nagaraja, 2015). In 2006, a standardized neighborhood deprivation index was created using Census tract level data from ACS (Messer et al, 2006). 15 Messer first included 20 variables and then used a principal component analysis to reduce them to eight variables: 1. % with less than a high school education 2. % unemployment rate, above 16 3. % crowded 4. % males in management and in professional occupations 5. % families in poverty 6. % female headed households on public assistance 7. % female headed households with dependent children 8. % households earning under $25,000/year A Vuong test was used to test three non-nested Poisson models (Vuong, 1989) where each model included all the personal covariates and one of the three indices. The results in Table 3 show that Messer’s index did a better job than Townsend’s and Jarman’s indices in fitting the data. Thus, Messer’s Neighborhood Deprivation Index was used in the following analyses. Table 3. Test for non-nested models Test Statistics 1.489 -7.731 -14.062 Townsend’s vs. Jarman’s Townsend’s vs. Messer’s Jarman’s vs. Messer’s p-value 0.136 0.000 0.000 All covariates mentioned above were included in the mixed multinomial logit part for the treatment group. For the count part of the FE-Hurdle model and FE-ZIP model, only some of the covariates were included. To decide which covariates to be included into the count part of the model, we ran a Poisson model with all the covariates and dropped those covariates with a high p-value (using %F as the outcome, p-value > 0.5). The results were shown below in z Table 4. According to the results, we selected the following covariates for the Hurdle and ZIP part: treatment groups, mom age groups, marital status, mom education levels, Hispanic, WIC program, Full Medicaid Before Conception, smoking, drinking, hypertension before 16 conception, previous preterm birth, rapid repeat pregnancies, pre-pregnancy BMI and Messer index. Table 4. Results from Poisson model Coef. Treatment group (ref: none enhanced prenatal health care program) SB program MIHP only Mother’s age group (ref: < 18) 18-24 25-29 30-34 35+ Marital status (ref: married) Unmarried, paternity acknowledged Mother only on birth certificate Education level (ref: < HS) HS >HS Hispanic ancestry WIC program Full Medicaid before conception Smoking Quit smoking Others in household smoked Drinking Diabetes prior to pregnancy Hypertension prior to pregnancy Previous preterm birth Rapid repeat pregnancies (ref: < 18 months from prior birth to current conception) >= 18 months from prior birth to current conception No prior deliveries Unknown Pre-pregnancy BMI Messer Index 17 Std. Err. p-value 0.39 0.30 0.06 0.04 0.000 0.000 0.48 0.43 0.26 0.15 0.10 0.11 0.12 0.15 0.000 0.000 0.028 0.318 0.29 0.31 0.07 0.06 0.000 0.000 -0.10 -0.13 -0.34 -0.05 0.56 0.09 -0.03 -0.01 0.25 0.00 0.14 0.22 0.04 0.05 0.14 0.04 0.04 0.05 0.09 0.05 0.26 0.14 0.11 0.06 0.022 0.008 0.017 0.195 0.000 0.104 0.703 0.897 0.342 0.980 0.183 0.000 0.07 0.04 0.112 0.01 -0.00 0.01 0.00 0.05 0.07 0.00 0.00 0.789 0.989 0.001 0.003 Chapter 4 Results Basically, in this study, two things were supposed to be checked. The first one is whether the estimations would change after using NYU algorithm to classify our ED visits for pregnant women. Two different outcomes were used in our models, one is the overall ED visits for each woman during pregnancy and the other one is the number of Non-emergent, Emergent/PC treatable and Emergent/Preventable ED visits based on the NYU algorithm. The second thing, which is the most important, is the validity of using the fixed-effects (FE) Hurdle model and FE-ZIP model when we have some potential selection bias or non-differential measurement error. Our data were used to fit six different models, including Poisson model, Negative Binomial model, ZIP model, FE-Negative Binomial model, FE-Hurdle model and FE-ZIP model. In Table 5, we gave the results from Poisson, Negative Binomial and Zero-Inflated Poisson model, using both %Fy and %F as the outcomes. z When we used %Fy as the outcome, the results we got were similar to the results if we used %F as the outcome. This might be because after use NYU algorithm, there were still a lot of z unclassified ED visits (nearly 30%). Among these unclassified ED visits, we could not know the distribution of non-emergent ED visits and “true” emergent ED visits, which made the results using %F similar to the results using %Fy as the outcome. z For now, let’s focus on the results using %F as the outcome. When we used Poisson and NB z models to fit the data, the coefficients for treatment groups are positive, which means that both SB group and MIHP only group have more ED visits during pregnancy than the no-MIHP group. For the ZIP model, the coefficients for the treatment groups in the count part are positive, while they are negative for the inflated part. Remember in the inflated part, the probability is the proportion for the part that generate structural zeros. So, in addition to the positive coefficients in the count part, the negative coefficients in the inflated part also mean that the 18 treatment groups (both SB and MIHP only) have more ED visits during pregnancy than the noMIHP group. Table 5. Results from Poisson, Negative Binomial and ZIP model SB Poisson NB ZIP, count part ZIP, inflation part Poisson NB ZIP, count part ZIP, inflation part Coefficients MIHP only 0.39* (0.28-0.51) 0.42* (0.30-0.53) 0.26* (0.14-0.37) -0.56* (-0.85 - -0.28) 0.41* (0.28-0.54) 0.44* (0.30-0.57) 0.25* (0.11-0.40) -0.51* (-0.81 - -0.21) {|~} 0.30* (0.23-0.38) 0.31* (0.24-0.39) 0.20* (0.12-0.28) -0.35* (-0.50 - -0.19) • {|} 0.31* (0.22-0.39) 0.32* (0.23-0.40) 0.20* (0.10-0.30) -0.30* (-0.48 - -0.12) Predictions – Rate Ratios SB vs. None MIHP only vs None 1.48 (1.33-1.66) 1.52 (1.35-1.70) 1.35 (1.25-1.45) 1.37 (1.27-1.47) 1.51 (1.34-1.68) 1.35 (1.25-1.45) 1.51 (1.33-1.71) 1.55 (1.36-1.76) 1.36 (1.25-1.48) 1.37 (1.26-1.49) 1.54 (1.34-1.74) 1.36 (1.25-1.48) Abbreviations: NB, Negative Binomial; ZIP, Zero-Inflated Poisson model; MIHP, Maternal Infant Health Program; SB, Strong Beginnings. 95% Confidence intervals in parentheses. *: statistical significance (p<0.05). Using the estimated coefficients, we held all the covariates except the treatment groups for each person constant and assumed that they were in SB group and calculated the predicted values of the outcome for SB group. Similarly, we calculated the predicted values of the outcome for MIHP only group and the reference group. Then we calculated the predicted rate ratio (RR) between SB group and the reference group and the predicted RR between MIHP only group and the reference group. The predicted RR for each model were listed in Table 5. Not surprisingly, for Poisson, NB and ZIP model, the RR between SB group and reference group and the RR between MIHP only and reference group are both greater than one. As discussed above, because of the measurement errors and self-selection, the results from Poisson, NB and ZIP model might be biased. We next used fixed-effects models to fit our data. 19 Table 6. Results from FE-NB, FE-Hurdle and FE-ZIP model Coefficients Predictions – Rate Ratios SB MIHP only SB vs. None MIHP only vs. None ~ {|} FE-NB -0.10 0.27* 0.92 1.22 (-0.29 – 0.09) (0.09 – 0.45) (0.75 – 1.10) (1.09 – 1.57) FE-Hurdle, -0.47* -0.51* count part (-0.64 - -0.30) (-0.60 - -0.41) 0.80 0.74 (0.69 – 0.92) (0.67 – 0.79) FE-Hurdle, -0.12 -0.31* hurdle part (-0.40 – 0.15) (-0.47 - -0.16) FE-ZIP, count -0.46* -0.47* part (-0.61 - -0.30) (-0.58 - -0.36) 0.67 0.63 (0.53 – 0.80) (0.49 – 0.76) FE-ZIP, -1.54 -0.60 inflation part (-4.50 – 1.37) (-1.82 – 0.59) • {|} FE-NB -0.11 0.26* 0.87 1.31 (-0.34 – 0.12) (0.12 – 0.40) (0.71-1.12) (1.13-1.50) FE-Hurdle, -0.52* -0.51* count part (-0.75 - -0.27) (-0.62 - -0.39) 0.81 0.75 (0.70-0.93) (0.69-0.81) FE-Hurdle, -0.16 -0.34* hurdle part (-0.43 – 0.11) (-0.49 - -0.19) FE-ZIP, count -0.52* -0.42* part (-0.70 - -0.35) (-0.59 - -0.26) --FE-ZIP, -12.94 -0.22 inflation part (-57.50 – 32.58) (-1.55 – 1.11) Abbreviations: FE-NB, Fixed-effect Negative Binomial; FE-Hurdle, Fixed-effect Hurdle model, FE-ZIP, Fixedeffect Zero-Inflated Poisson model; MIHP, Maternal Infant Health Program; SB, Strong Beginnings. *: statistical significance (p<0.05) Table 6 shows the estimates of the parameters from three Fixed-effect models. The estimates for FE-NB model came from the Stata command “mtreatnb” (Deb et al, 2006), and we programmed the MLE routine to get the estimates for FE-hurdle model and FE-ZIP model using simulated maximum likelihood estimation method (the codes were shown in Appendix). For the FE-Hurdle model and FE-ZIP model, as explained in Chapter 3, they generate two separate models. For the Hurdle model, we used FE truncated Poisson model for those with at least one visit and used an FE logit model to explain whether the number of visits is zero of not. For the FE-ZIP model, an FE Poisson model was used for those who had some probabilities to have an ED visit during pregnancy and an FE logit model was used to explain the inflated part. People in the inflated part was supposed to have no probability to visit Emergency 20 Department during pregnancy, this might due to a “very” healthy pregnancy. Same as before, the results using different outcomes were very closed to each other. Let’s look at the results (using %F as outcome). In the FE-Hurdle model, the estimated coefficients for z the treatment groups in the “Count” part are -0.52 and -0.51, which means that both SB women and MIHP only women had fewer ED visits during pregnancy. For the logit part, the estimated coefficients for the treatment groups are -0.16 and -0.34. In our FE-Hurdle model, the probability for the logit part in the left-hand side is the probability of having positive outcomes. Negative coefficients for the treatment groups mean that both SB program and MIHP only have fewer ED visits during pregnancy. The prediction values for the RR’s using the estimated coefficients in FE-NB and FE-Hurdle models were less than one except for the predicted RR for MIHP only vs. None in the FE-NB model, which means that women in both SB group and MIHP only group would have fewer ED visits during pregnancy. We used bootstrap to get the 95% confident interval (CI) for the predicted RR for FE-Hurdle and FE-ZIP models. We bootstrap 800 times for each model and found the 2.5 percentile as the lower bound of the 95% CI and the 97.5 percentile as the upper bound. The 95% CI showed that the predicted rate ratios for FE-Hurdle model and FE-ZIP model were statistical significant different from one. For the FE-ZIP model, the estimated coefficients for the treatment groups in the inflated part are -12.94 (SB vs. None) and -0.22 (MIHP only vs. None) which meant that those people who enrolled in SB and MIHP only group had lower probabilities to create structural zeros. This may be because those women who enrolled in these enhanced prenatal health programs were supposed to have more risk in their pregnancies. But the predictions based on these estimated coefficients meaningless, which may require further studies. Although the RR from FE-NB model was not statistically significant, it also showed that the results from regular count models such as Poisson model were biased. 21 Table 7. Results for Fixed-Effect Hurdle model Parameters Treatment group (ref: no prenatal program) SB a MIHP only ab Mom age group (ref: < 18) 18 – 24 ab 25 - 29 ab 30 - 34 a 35 + Marital status (ref: married) Paternity ab Mom only ab Education level (ref: < HS) HS ab > HS ab Hispanic a WIC program Full Medicaid before conception ab Smoking Drinking Hypertension prior to pregnancy Previous preterm birth ab Rapid repeat pregnancies (ref: < 18 months from prior birth to current conception) >= 18 months from prior birth to current conception Count Coef. Std. Err. Logit Coef. Std. Err. -0.52 -0.51 0.12 0.06 -0.16 -0.34 0.14 0.08 0.65 0.59 0.38 0.41 0.14 0.16 0.18 0.24 0.91 0.65 0.38 -0.16 0.18 0.20 0.22 0.25 0.24 0.27 0.12 0.12 0.28 0.40 0.12 0.12 -0.17 -0.24 -0.61 -0.05 0.68 0.11 -0.30 0.08 0.29 0.07 0.08 0.25 0.07 0.07 0.08 0.23 0.14 0.10 -0.14 -0.29 -0.52 0.07 0.98 0.18 0.27 0.38 0.47 0.09 0.10 0.30 0.08 0.08 0.11 0.37 0.21 0.15 0.07 0.08 0.25 0.09 0.13 0.14 0.01 0.01 0.09 0.12 0.00 0.00 -0.04 0.06 0.01 0.01 0.11 0.17 0.00 0.00 b No prior deliveries Unknown Pre-pregnancy BMI b Messer Index ab Abbreviations: MIHP, Maternal Infant Health Program; SB, Strong Beginnings. a : statistical significant difference (p<0.05) for Count part (versus reference group). b : statistical significant difference (p<0.05) for Logit part (versus reference group). For the other covariates, such as mom’s age and education level, the estimated coefficients and standard errors were listed in Table 7. For mothers’ age groups, comparing with those who gave births before 18 years old, those moms whose age was between 18 and 29 years old were more likely to have more ED visits during pregnancy. For marital status, those unmarried women were more likely to have more ED visits during pregnancy compared with those married women. Higher education level for mom seems to help reducing the number of ED visits during pregnancy. Those women with full Medicaid before conception had more ED visits during pregnancy than those women without full Medicaid before conception. Another important 22 covariate was the previous preterm birth, those women who had preterm birth previously had more ED visits during pregnancy significantly than those without preterm birth previously. 23 Chapter 5 Discussion Firstly, as discussed in the section of results, due to the large proportion of unclassified ED visits, there wasn’t any significant difference when we used %F as the outcome. For this study, z we would like to know if the enhanced prenatal health care could reduce the number of ED visits or not. The ED visits here should exclude those non-preventable emergent situations, because the enhanced prenatal health cares are not designed to help with these situations. That’s the reason why we wanted to use the NYU algorithm to classify our ED visits. However, there are lots of diagnosis codes which cannot be classified based on the NYU algorithm. That makes our outcome %F unclear. Further research is needed to refine the NYU ED algorithm. This is z also a limitation of using claims data. Nevertheless, according to Table 5, we can see if we used the usual count models, such as Poisson, NB or ZIP model, the predictions of the RR’s (SB vs. None, MIHP only vs. None) were greater than one, which means that the MIHP and SB program increased the number of ED visits during pregnancy. On the other hand, those women who enrolled in MIHP and SB program were supposed to have higher risk about their pregnancy, and this makes them to more likely to get ED visits during pregnancy. Because of the self-selection bias, it was not a surprise that the predictions of the RR’s (SB vs. None, MIHP only vs. None) changed when we turn to FE models. The results from FE models showed that the both MIHP and SB programs would decrease the number of ED visits during pregnancy. Since the emergence department resources are very valuable, it is better to reduce those unnecessary ED visits. The results show that the enhanced prenatal health care programs such as MIHP and SB can efficiently reduce the number of ED visits during pregnancy for African American women. 24 APPENDICES 25 APPENDIX A Figure 1. Flow Chart 26 APPENDIX B Figure 2a. Distribution of the overall outcome in each treatment group. Figure 2b. Distribution of the outcome (NE, PCT and PA using NYU algorithm) in each treatment group. 27 APPENDIX C Table A1. Number of Missing Values Variables Mother’s age group Marital status Paternity Education WIC program Smoking Others in household smoked Drinking Diabetes prior to pregnancy Hypertension prior to pregnancy Previous preterm birth Rapid repeat pregnancies Pre-pregnancy BMI # of missing values, N (%) 2 (0.03) 28 Imputed values 3 834 (14.32) 12 (0.21) 11 (0.19) 4 (0.07) 9 (0.15) 1 (0.02) 1 (0.02) 1 (0.02) 0 2 0 0 0 0 0 0 1 (0.02) 330 (5.66) 118 (2.03) 0 4 28.86 APPENDIX D ****************************Fixed-effect Hurdle Model*************************** cap program drop myfehurdle program myfehurdle args lnfj lnmu theta ome2 ome3 tempvar pi1 pi2 pi3 L1 mu p l pi qui { gen double `pi1' = 0 gen double `pi2' = 0 gen double `pi3' = 0 gen double `pi' = 0 gen double `mu' = 0 gen double `p' = 0 gen double `L1' = 0 gen double `l' = 0 } foreach x of num 1/500 { qui { replace `pi1' = 1 / (1+exp(`ome2'+l2`x')+exp(`ome3'+l3`x')) replace `pi2' = exp(`ome2'+l2`x') * `pi1' replace `pi3' = exp(`ome3'+l3`x') * `pi1' replace `pi' = `pi1'*d0+`pi2'*$ML_y2+`pi3'*$ML_y3 replace `mu' = exp(`lnmu'+l2`x'+l3`x') replace `p' = exp(`theta'+l2`x'+l3`x') / (1+exp(`theta'+l2`x'+l3`x')) replace `l' = (1-`p') * `pi' if $ML_y1 == 0 replace `l' = `p' * exp(-`mu') * `mu'^$ML_y1 / (1-exp(-`mu')) / round(exp(lnfactorial($ML_y1 )), 1) * `pi' if $ML_y1 > 0 replace `L1' = `L1' + `l' } } quietly replace `lnfj' = ln(`L1'/500) end ********************************Fixed-effect ZIP Model************************** cap program drop myfezip program myfezip args lnfj lnmu theta ome2 ome3 tempvar p1 p2 p3 L1 mu p l pi qui { gen double `p1' = 0 gen double `p2' = 0 gen double `p3' = 0 gen double `p' = 0 gen double `mu' = 0 gen double `pi' = 0 gen double `L1' = 0 gen double `l' = 0 } foreach x of num 1/500 { qui { replace `p1' = 1 / (1+exp(`ome2'+l2`x')+exp(`ome3'+l3`x')) replace `p2' = exp(`ome2'+l2`x') * `p1' replace `p3' = exp(`ome3'+l3`x') * `p1' replace `p' = `p1'*d0+`p2'*$ML_y2+`p3'*$ML_y3 replace `mu' = exp(`lnmu'+l2`x'+l3`x') replace `pi' = exp(`theta'+l2`x'+l3`x') / (1+exp(`theta'+l2`x'+l3`x')) replace `l' = (`pi'+(1-`pi')*exp(-`mu')) * `p' if $ML_y1 == 0 replace `l' = (1-`pi') * exp(-`mu') * `mu'^$ML_y1 / round(exp(lnfactorial($ML_y1 )), 1) * `p' if $ML_y1 > 0 replace `L1' = `L1' + `l' } } quietly replace `lnfj' = ln(`L1'/500) end 29 REFERENCES 30 REFERENCES Cameron, A., and Trivedi, P. (2005). Microeconomics: Methods and Applications. 1st ed. Cambridge: Cambridge University Press. Cameron, A. and Trivedi, P. (2013). Regression analysis of count data. Cambridge: Cambridge University Press. Cragg, J. (1971). Some Statistical Models for Limited Dependent Variables with Application to the Demand for Durable Goods. Econometrica, 39(5), p.829. Falik, M., Needleman, J., Wells, B. and Korb, J. (2001). Ambulatory Care Sensitive Hospitalizations and Emergency Visits: Experiences of Medicaid Patients Using Federally Qualified Health Centers. Medical Care, 39(6), pp.551-561. Freedman, D. (2006). Statistical Models for Causation: What Inferential Leverage Do They Provide?. Evaluation Review, 30(6), pp.691-713. Gardiner, J., Luo, Z. and Roman, L. (2009). Fixed effects, random effects and GEE: What are the differences?. Statistics in Medicine, 28(2), pp.221-239. Gourieroux, C. and Monfort, A. (1996). Simulation-based econometric methods. 1st ed. [Oxford]: Oxford Scholarship Online. Grafarend, E.W. (2006). Linear and Nonlinear Models: Fixed Effects, Random Effects, and Mixed Models. Germany: Walter de Gruyter. Hilbe, J. (2013). Negative binomial regression. 1st ed. Cambridge: Cambridge University Press. Hsiao, C. (2003). Analysis of Panel Data. 2nd ed. Cambridge: Cambridge University Press. Hu, M., Pavlicova, M. and Nunes, E. (2011). Zero-Inflated and Hurdle Models of Count Data With Extra Zeros: Examples from an HIV-Risk Reduction Intervention Trial. The American Journal of Drug and Alcohol Abuse, 37(5), pp.367-375. Imbens, G. and Rubin, D. (2015). Causal inference for statistics, social, and biomedical sciences. 1st ed. Jones, K., Paxton, H., Hagtvedt, R. and Etchason, J. (2013). An Analysis of the New York University Emergency Department Algorithm’s Suitability for Use in Gauging Changes in ED Usage Patterns. Medical Care, 51(7), pp.e41-e50. Lambert, D. (1992). Zero-inflated Poisson regression, with an application to defects in manufacturing. Technometrics. Volume 34, pp.1-14. Majo, M.C., van Soest, A.H.O. (2011). The Fixed-Effects Zero-Inflated Poisson Model with an Application to Health Care Utilization. (CentER Discussion Paper; Vol. 2011-083). Tilburg: Econometrics. 31 Messer, L., Laraia, B., Kaufman, J., Eyster, J., Holzman, C., Culhane, J., Elo, I., Burke, J. and O’Campo, P. (2006). The Development of a Standardized Neighborhood Deprivation Index. Journal of Urban Health, 83(6), pp.1041-1062. Mnatzaganian, G., Davidson, D.C., Hiller, J.E., and Ryan, P. (2015). Propensity score matching and randomization. Journal of Clinical Epidemiology, 68(7), pp.760-768. Mullahy, J. (1986). Specification and testing of some modified count data models. Journal of Econometrics, 33(3), pp.341-365. Nagaraja, C.H. (2015). Deprivation Indices for Census Tracts in Bronx and New York Counties. Partha, D., Pravin K.T. (2006). Maximum simulated likelihood estimation of a negative binomial regression model with multinomial endogenous treatment. The Stata Journal, Volume 6, Number 2, pp.246-255. Patil, G. and Boswell, M. (1970). A Characteristic Property of the Multivariate Normal Density Function and Some of its Applications. The Annals of Mathematical Statistics, 41(6), pp.19701977. Pearl, J. (2000). Causality: Models, Reasoning, and inference. Cambridge: Cambridge University Press. Peter, H., Arne, U. (2006). Estimation of multinomial logit models with unobserved heterogeneity using maximum simulated likelihood. The Stata Journal, Volume 6, Number 2, pp.229-245. Vuong, Q. (1989). Likelihood Ratio Tests for Model Selection and Non-Nested Hypotheses. Econometrica, 57(2), p.307. Wooldridge, J. (2002). Econometric Analysis of Cross Section and Panel Data. 1st ed. Cambridge, Mass.: MIT Press. 32