Maximum likelihood based analysis of equally spaced longitudinal count data with first-order antedependence and overdispersion

This manuscript implements a maximum likelihood based approach that is appropriate for equally spaced longitudinal count data with over-dispersion, so that the variance of the outcome variable is larger than expected for the assumed Poisson distribution. We implement the proposed method in the analysis of seizure data and a subset of German Socio-Economic Panel data. To demonstrate the importance of correctly modeling the over-dispersion, we make comparisons with the semi-parametric generalized estimating equations approach that incorrectly ignores any over-dispersion in the data. Our simulations demonstrate that accounting for over-dispersion results in improved small-sample efficiency and appropriate coverage probabilities. We also provide code in R so that readers can implement our approach in their own analyses.

linearity of the expectations for the conditional distributions, which are assumed to be Poisson. In addition, we assume that the correlation between adjacent measurements on a subject is constant.
The assumptions of the first-order Markov property, linearity in the conditional expectations, and constant adjacent correlations have been shown to induce a first-order autoregressive AR(1) correlation structure for the repeated outcomes on each subject (Guerra and Shults 2014). The AR(1) structure forces a decline in the intra-subject correlations with increasing separation in time. Our method is therefore most appropriate for analysis of equally spaced longitudinal count data with over-dispersion.
Other approaches for analysis of over-dispersed longitudinal count data include semiparametric approaches such as generalized estimating equations (GEE) (Liang and Zeger 1986). Vinod (2002) described econometric applications of GEE. Ghisletta and Spini (2004) provided a concise summary of GEE for the social sciences. GEE is widely used because it does not require specification of the full likelihood that can be quite complex for longitudinal discrete data. However, GEE does not account for over-dispersion. In addition, the relative ease of application of GEE for discrete data can also be a potential limitation for the approach. When only the first two moments of the distribution of the outcome variable are estimated, as they are for GEE, it is possible to obtain estimates that are not compatible with any valid parent distribution. As cautioned by Molenberghs and Kenward (2010), "the parent provides a natural description of the framework into which the semi-parametrically specified parameters fit. The implication is that such semi-parametric methods as GEE1, GEE2, ALR, etc. can always be applied because there is always a valid parent, and hence a probabilistic basis. " We make comparisons with GEE to evaluate the impact of incorrectly ignoring overdispersion when the models for the marginal mean and correlation structure are correct. We conduct simulations for moderately sized samples to demonstrate that when the likelihood is correctly specified, we have improved efficiency in estimation of the regression and correlation parameters for our approach relative to GEE that incorrectly ignores the over-dispersion.
Another model for longitudinal count data is the class of generalized linear mixedeffects models that incorporate random effects in the linear predictor. However, the implementation of likelihood based methods that involve random effects can be computationally challenging (p. 75, Fitzmaurice et al. 2008). In addition, in contrast to GEE, for mixed models it is not straightforward to specify a particular working correlation structure for the repeated measurements on subjects. For example, the AR(1) correlation structure is not among the covariance models that were suggested by Thall and Vail (1990). Mixed-effects models are typically employed when the goal is to estimate effects that are subject specific, because the analysis results are conditional on the random effects (Gardiner et al. 2009).
In general, likelihood based approaches like the one we implement in this paper enjoy several general advantages. Unlike semi-parametric approaches, they yield an estimated likelihood that can be used to conduct likelihood ratio tests and to compare the fit of models using criteria such as the Akaike information criterion (AIC) (Akaike 1974) and Bayesian information criterion (BIC) (Schwarz 1978). Maximum likelihood estimators are also most (asymptotically) efficient among a wide class of estimators (Serfling 2011) when the distribution is correctly specified. Our method in particular, allows for specification of the usual model for the marginal mean for Poisson data, while also accounting for over-dispersion and serial correlation in the data via the induced AR(1) correlation structure.
In "Methods" section we discuss the notation, model assumptions, the likelihood and likelihood equations. In 'Application" section we discuss an application of the methods followed by simulation studies in "Simulation studies" section. We conclude with a discussion in Conclusion" section.

Notation and model assumptions
The data comprise realizations y ij of ordered discrete random variables Y ij that are measured on subject i at time t ij (i = 1, . . . , m and j = 1, . . . , n i ). Associated with each y ij is a vector of explanatory variables (covariates) The expected value of measurement Y ij on subject i is given by and the variance by var (Y ij ) = σ 2 ij . We assume that observations on different subjects are independent. Further, the measurements within subjects are correlated with a structure that depends on parameter α . Let cov (Y ij , Y ik ) represent the covariance and corr (Y ij , Y ik ) represent the correlation between Y ij and Y ik .
We make three assumptions. First, we assume first-order antedependence, such that each Y ij , given the immediate antecedent Y ij−1 , is independent of all further preceding variables (Gabriel 1962). The joint probability mass function of Y i1 , . . . , Y in i can then be expressed as First-order antedependence is also referred to as the first-order Markov property in the literature (Feller 1968, p. 419).
Second, we assume that the correlation between adjacent measurements on a subject is constant, implying that where i = 1 . . . , m and j = 2, . . . , n i . Third, we assume that the conditional expectation of Y ij given Y ij−1 is a linear function of Y ij−1 , such that for i = 1 . . . , m and j = 2, . . . , n i .
These three assumptions imply the following results. From Theorem 2.1 of Guerra and Shults (2014), the conditional expectation is given by where i = 1, . . . , m and j = 2, . . . , n i . Next, from Section 2.5 of Zimmerman and Nuñez-Antón (2010), the correlation corr (Y ij , Y ij+t ) between Y ij and Y ij+t for t > 0 can be expressed as The induced correlation structure for Y i1 , . . . , Y in i ′ is therefore an AR(1) structure.
This AR(1) structure is plausible for longitudinal data because it requires the correlation between measurements on a subject to decline with increasing separation in time. For example, if α = 0.5, then the correlation between the 1st and 2nd measurements is 0.5, while the correlation between 1st and 3rd measurements is (0.5) 2 = 0.25.

Poisson likelihood
We assume Poisson distributions for the marginal and conditional distributions in Eq. 2. For each i = 1, . . . , m, the distribution of Y i1 is Poisson with µ i1 = i1 = exp x ′ i1 β and σ i1 2 = i1 , where β is a p × 1 vector of regression parameters. Then, for j = 2, . . . , n i , the conditional distribution of Y ij given Y ij−1 is Poisson with conditional mean E Y ij |Y ij−1 = ij * given by Eq. 3, with and for j = 2, . . . , n i and i = 1, . . . , m. The Y ij are over-dispersed relative to the Poisson distribution if j ≥ 2 and α � = 0, because in this case σ ij 2 = φ ij , where φ > 1.

Likelihood equations
To obtain maximum likelihood estimates of β and α, we need to obtain simultaneous solutions to the following estimating equations for β and α, respectively: and The derivatives are provided in Appendix A of the longer, working version of this paper (Gamerman et al. 2016 at http://biostats.bepress.com/upennbiostat/art45). We maximized the likelihood using an adaptive barrier algorithm as implemented in the constrOptim function in R (R Core Team 2014). We applied the Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimization method by Broyden (1970), Fletcher (1970, Goldfarb (1970), and Shanno (1970), Shanno and Kettler (1970), which is implemented in constrOptim when the gradient is provided.
The following algorithm summarizes our estimation procedure for a particular model: 1. Choose initial estimates (starting values) of α and β. Starting values can be obtained using GEE to fit a Poisson model with an AR(1) correlation structure; however, we should check that the starting values satisfy the constraints ("Poisson Likelihood"). If the estimates violate the constraints, change the starting values by choosing a value for α that is closer to zero or by applying Poisson regression, which is equivalent to assuming that α = 0. 2. Obtain solutions to the likelihood Eqs. 7 and 8 using the adaptive barrier algorithm that is implemented in the R package constrOptim. R code for the log likelihood function and for the gradient function, both of which are implemented in the Application, are provided in Appendix B of Gamerman et al. 2016 (at http://biostats. bepress.com/upennbiostat/art45).

Asymptotic distribution of the estimators
If the model is correctly satisfied and standard regularity conditions are satisfied, the ML approach described here will yield estimates that are consistent and asymptotically normal. Let θ = (β, α) T and the maximum likelihood estimators θ = (β,α) T . We estimated the asymptotic covariance matrix of θ with the inverse of the observed information, (i(θ )) −1 , that we estimated using the inverse of the negative Hessian matrix, which is defined and implemented in Appendices A and B, respectively, of Gamerman et al.

Doctor visits data
Here we consider an analysis of a subset of data from the German Socio-Economic Panel data (Winkelmann 2004) that we obtained within Stata (http://www.stata-press.com/ data/r14/drvisits) and then exported for analysis in R (StataCorp 2013). Here we compare the results of an analysis using the proposed ML approach with the results obtained using Poisson regression and GEE. The goal of the analysis was to assess the impact of the 1997 health reform on the reduction of government expenditures. A sample of 1518 women who were employed full time in the year before or after the reform was implemented were evaluated. The outcome we considered was the self-reported number of doctor visits in the three months prior to the interview. The main covariate of interest was an indicator variable that took value 1 if the interview took place after the reform was implemented (or took value 0 otherwwise). Additional covariate information was available on each participant's age, education, marital status, self-reported health status, and the logarithm of household income. Of the 1518 women in the dataset, 709 were interviewed both before and after the reform was implemented; 391 were only interviewed before; and 418 were only interviewed after the reform went into effect. This resulted in a total of 2227 observations available for the analysis.
We assumed Eq. 5 with the following linear predictor: where x ij1 was the indicator variable for health care reform (1 if after implementation; 0 if before), x ij2 was age in years, x ij3 was education in years, x ij4 was marital status (1 if married; 0 if not married), x ij5 was self-reported health status (1 if bad; 0 if not bad), and x ij6 was the logarithm of household income. We first fit the above model using Poisson regression as implemented in the glm function in R; the results are provided in Table 1. Among women with the same household income, marital status, self-reported health, and education, there was a reduction in the log count of doctor visits of 0.140 after health care reform was implemented (p < 0.0001 ).
Next, we used the geeglm function in R to implement GEE with an assumed AR(1) working correlation structure; the results are shown in Table 1. As for Poisson regression, there was a significant reduction in the log count of doctor visits (β 1 = −0.123, p = 0.0202). The estimated correlation parameter was 0.213.
When we fit the GEE model we assumed that the scalar parameters φ = 1 ∀ i, j. After fitting GEE, we assessed the adequacy of this assumption by obtaining an estimate of φ based on the final GEE estimates of β: . The estimated φ was φ = 4.33, which is much greater than 1 and was therefore suggestive of over-dispersion in the data.
Lastly, we fit the proposed ML approach using the algorithm for estimation described in "Likelihood equations" section. We obtained starting values for our approach using GEE, after first confirming that α satisfied the necessary constraint to guarantee a valid parent distribution, which in this case was α < 0.4494. Table 1 shows the results for the ML approach. The estimated correlation parameter was 0.313 with a 95% confidence interval of (0.272, 0.354). After adjusting for the correlation among the counts of doctors visits, for over-dispersion, and for the other covariates, we again found that there was a significant impact of initiation of health care reform on the number of doctor visits (β 1 = −0.113, p < 0.0001).
Overall, the parameter estimates were similar for the proposed ML approach, GEE, and the Poisson regression. While the impact of age was similar across the approaches, it was significant in both the ML and Poisson approaches but not significant in the GEE model (ML p = 0.0005, GEE p = 0.1182, and Poisson p = 0.0008). Similarly, the logarithm of household income was significant in both the ML and Poisson approaches but not significant in the GEE model (ML p < 0.0001, GEE p = 0.0809, and Poisson p < 0.0001).
With estimates of the log-likelihood for Poisson regression and the proposed ML approach, it was possible to calculate the AIC and BIC criteria as measures of the relative quality of the models for this set of data. Both BIC and AIC incorporate a penalty term for the number of parameters used in the model because it is possible to increase the numerical value of the likelihood solely by including additional parameters in the model, which may result in over-fitting the model to the data. This penalty term is larger in the BIC as compared to the AIC. For the Poisson regression model, the AIC and BIC values were 11,899 and 11,939, which were both greater than the AIC and BIC values for the ML approach (AIC = 11,707 and BIC = 11,746), which indicates that the ML approach had improved model fit over Poisson regression.

Epilepsy seizure data
Here we implement the proposed ML method and GEE for analysis of the epilepsy seizure data (Thall and Vail 1990;Farewell and Farewell 2013) that is available as part of the MASS package in R (Venables and Ripley 2002). We assumed Eq. 5 with the following linear predictor: where x ij1 represents an indicator for treatment, x ij2 represents baseline seizure count (number of seizures in the 3 month time period prior to the start of the study), x ij3 represents subject age in years, and x ij4 represents two-week time period (coded as 1, 2, 3, 4). We initially included a time period by treatment interaction term, but the interaction term was not significant for the proposed approach or for GEE (both p-values > 0.05); we therefore initially focused on the simpler model 9 for this demonstration. Table 2 shows the sample mean and variance of seizure counts at baseline and the four subsequent two-week periods (denoted as Y1 through Y4) for the placebo and drug groups for the seizure counts; it also displays the sample mean and variance of age at baseline. From the table, the sample variance for the outcome variables, Y1 through Y4, were greater than their respective means, which suggested that there was over-dispersion in the seizure counts. Table 3 shows the results of the analysis. The estimates were similar for the proposed ML method and GEE. The estimate of treatment was negative for both approaches, which suggested that the number of seizures was lower for subjects in the treatment group. However, treatment only differed significantly from 0 for the proposed ML approach (p = 0.0124 for ML versus p = 0.3014 for GEE). In addition, time period only differed significantly from 0 for the proposed ML approach (p = 0.0032 for ML versus p = 0.0580 for GEE). The likelihood ratio test of the hypothesis that the regression parameter for time period is 0 also suggested that time period should be retained in the model for the proposed ML approach (p = 0.0030.) However, since the GEE analysis suggested that time period might not be important, we removed time period from the model for both GEE and the proposed ML approach. As shown in Table 4, treatment differed significantly from 0 for the proposed ML approach, but was not significant for GEE (p = 0.0121 for ML versus p = 0.2977 for GEE).
We next compared the AIC and BIC for the models that included and excluded time period. As shown in the Tables, both the AIC and BIC values were smaller for the larger model that included time period. The respective AIC and BIC values were 1566 and

Simulation studies
In the previous section we identified significant treatment effects for the proposed ML approach that were not observed for GEE. Since the results depended on choice of approach, it was of interest to compare the performance of the methods for finite samples. We therefore performed simulations to assess the properties of the estimators of α and β for the proposed ML approach and GEE.

Set-up
We compared the performance of the ML and GEE estimators for where the x ijk were defined in the previous section.
Covariates were simulated based on the observed epilepsy seizure data in the previous section. Treatment was specified as present (equal to 1) for one group and as absent (equal to 0) for the other group. Baseline seizure count was simulated from a Poisson distribution with a random seed and mean = 31.22 based on the mean baseline age from the epilepsy data. Similarly, age was simulated from a normal distribution based on the epilepsy data for which the minimum age was 18, the mean was 28.3, and the standard deviation was 6.261. Simulated age values below 18 were discarded and the next (10) simulated age value was assigned. Age was then rounded to a whole number, as it was recorded in the epilepsy data. The approach proposed by Guerra and Shults (2014) was used to simulate the correlated Poisson seizure counts with specified means, over-dispersion, and AR(1) correlation structure.

Assessments
We wrote code in R to evaluate mean square error (MSE), percent bias, small sample efficiency, and 95% coverage probabilities using the observed information matrix. The mean square error (MSE) for estimator θ is defined as where θ is the true value. The percent bias for estimator θ is defined as Lastly, to evaluate the coverage probabilities, a 95% confidence interval was computed for each parameter estimate within each simulation run. The coverage probabilities represent the proportion of the R simulation runs in which the true parameter fell within the 95% confidence bounds. GEE coverage probabilities were computed similarly using the naïve variance estimates obtained from geeglm in R.

Results
Table 5 displays the MSE and Table 6 displays the percent bias for the simulations. For the ML method, the MSE for β and α and the percent bias for α decreased as m increased. As compared to GEE, the ML approach had lower MSE and percent bias for all sample sizes for α. For β , the percent bias was similar for ML and GEE; however, the MSE was slightly smaller for ML than for GEE. For scenarios with high correlation (α = 0.6 or 1 R R i=1 θ −θ i /θ * 100. 0.7), the intercept and treatment estimates, β 0 and β 1 , had smaller MSE and percent bias for the proposed ML approach than for GEE, for all samples sizes. Table 7 then displays the estimated coverage probabilities. With respect to β , the coverage probabilities were similar for the ML and GEE approach and were close to the nominal 95% level. With respect to α, the ML approach model-based coverage probabilities were close to the nominal 95%, which outperformed the GEE approach, whose model-based coverage probabilities were below the nominal 95% level. Coverage probabilities for α were better for the ML based approach than GEE across all sample sizes and correlations (α = 0.2, 0.4, 0.6, 0.7).

Conclusion
We proposed an ML approach for analysis of equally spaced longitudinal count data that accounts for intra-subject correlation of measurements and over-dispersion. Our application of the ML approach and GEE demonstrated that the results of the analysis differed between approaches, with significant treatment differences observed for some models for the ML approach, but not for GEE. The availability of the AIC and BIC criteria for the ML approach was useful for selecting between nested models. The interested Table 7 Coverage probabilities for the ML and GEE approaches with the AR(1) correlation structure for varying values of α and sample size per group The true correlation structure is AR (1) There are equal sample sizes of m 2 per group and β = (β 0 , β drug , β baseline , β age ) ′ = (0.4467, −0.1659, 0.0232, 0.0258) ′