The modeling of medical expenditure data from a longitudinal survey using the Generalized Method of Moments (GMM) approach

1 Introduction

The rising medical expenditures have figured in the news rather prominently in the last several years and have prompted a substantial interest in the analysis of healthcare related data. Rising health care expenditures for older adults have provoked a lot of concern recently. The goal of our research is to investigate whether unmet need for help with disabilities in activities of daily living (ADL) is associated with higher medical expenditures. If this is the case, such a connection would inform policy makers about resource planning for older adults with unmet ADL need. Unmet ADL need is a serious problem among older adults. 15% of older adults need help from others to complete basic activities of daily living (ADL) [1] such as bathing, dressing, eating, toileting, and getting around inside. However, nearly 20% of older adults who need ADL help report unmet need for assistance with their ADL [2, 3, 4]. Unmet need for ADL assistance is associated with increased healthcare utilization including hospitalization [5], re-hospitalization [6], and nursing home placement [7]. Unmet need for ADL assistance is also associated with increased risk for death [8] and, for many older adults, medical care expenditures increase significantly in the months before death [9].

The medical expenditure data are commonly available from a number of longitudinal surveys such as the NLTCS (National Long Term Care Survey). NLTCS is one of the longest running longitudinal surveys in the USA that has been ongoing for more than 30 years and it is linked to Medicare claims data. Since the detailed description of NLTCS has been given already in [10], we only offer a very brief one. Let the starting point of our observation period be Sept 1st, 2004. Beginning with this date, subjects begin entering the study starting from the day of the interview. The interview day can fall on any day between Sept 1st, 2004 and Dec 31st, 2005. As a result, we observe the straggled entry of subjects. Some subjects die during this follow-up period while others survive until Dec 31st, 2005 at which point the follow-up period stops for all subjects. Note that medical expenditures are taken from linked Medicare claims data; the reason follow-up was only conducted until the end of 2005 had to do with availability of the pre-linked Medicare records. The total medical expenditures are subdivided into a number of categories, e.g. durable medical equipment, hospital expenditures, skilled nursing facility, home health agency etc. For each subject, a number of covariates are available, such as unmet need for ADL disabilities, the number of ADL disabilities, and age. All of the covariates are binary. In addition to these two, some others were the diabetes status, the heart disease status etc. as well as demographic covariates such as age and gender. The total number of subjects was 2400 of which 467 died during the period of study.

The medical expenditure data tend to be rather complicated data and they present numerous statistical analysis challenges. First, they tend to be highly skewed to the right, with a relatively small proportion of patients incurring very high medical costs while the rest of patients hardly incurring any. Second, these data usually have a lot of missing observations either due to the lack of follow-up or death related dropout or both factors. Since the probability of death is related to expenditures due to significantly higher expenditures for many in the last months of life, the missing observations due to death cannot be viewed as MCAR (missing completely at random) which makes the analysis even harder. Third, the common simple random sampling (SRS) assumption cannot be used in their analysis since each observation point has the survey-related weight. We now discuss these issues in some extra detail.

The presence of skewness in the data, together with a large number of zeros, implies that the choice of the modeling distribution may not be very straightforward. In general, such data cannot be viewed as generated by any particular continuous distribution. One of the approaches used to treat this problem is the use of OLS (ordinary least squares) regression with a positive shift at zero. For an overview of this approach, see, e.g. [11]. This approach has two significant shortcomings. First, the choice of the constant used to shift all of the observations away from zero tends to be rather arbitrary. Moreover, a retransformation back to the original scale is required in this case after the model has been fitted. Another possibility is to use the Tobit model which is, effectively, a censored normal regression that is based on the concept of latent variables. The genesis of this idea also goes back to econometric research; for details and a good overview, see [12] and [13]. The use of the Tobit model is problematic because it is very sensitive to violation of normality and heteroscedasticity assumptions (see [14] for more details). Moreover, Tobit model assumes that there is an underlying normal random variable that is censored due to some random mechanism; this implies, effectively, that zeros are not viewed as a valid response which is typically wrong in the medical cost data context.

Finally, the so-called two part model envisions a logit or probit GLM to model the probability of zero occurring while using another OLS or GLM to model the actual level of positive cost. This approach effectively models the fact that excessive zeros may be generated by a mechanism different from that of positive expenditures. Note that this amounts to the use of a degenerate mixture model where one of the components is concentrated at just one point. The two-part model has a long and distinguished history in various applications. A version of this model was used in 1970’s by meteorologists for rainfall; see, e.g. [15], [16] and [17]. The first ever example of its use in economic context was probably [18]. Later, this model was widely used in health economics as a result of the well known Health Insurance Experiment conducted by RAND Corporation; in that context, the two-part model was introduced in [19] and [20]. [21] provided a good overall review of the widespread use of the two-part model for health care cost data. Note that in the cross-sectional context the two parts of the model may be fitted separately. An excellent recent work on the practical implementation of the two-part model in the cross-sectional context is [22]. This approach was later extended to the longitudinal data context; the first occurrence was, probably, in [23].

The presence of missing data because of death also creates a significant problem in the data analysis. Some options considered so far in the literature include using an estimated probability of survival, obtained using, for example, Kaplan-Meier approach, within a short subinterval as a weight and then summing up the mean total cost weighted by it over all of the intervals. Such an approach was first suggested in [24]. [25] managed to extend this approach to develop an estimator whose asymptotic distribution is independent of the choice of partition. A different approach models the hazard function of the terminal event (i.e. death) based on subject specific covariates as a part of the joint model. For more detailed discussion of this approach, see [26] and [10].

Finally, the fact that the data come from a longitudinal survey also needs to be taken into account. It has been long known that ignoring sampling weights can lead to severely biased parameter estimates with underestimated standard errors (see [27] for a detailed discussion of this issue in the medical context). There has been relatively little research on how to account for sampling weights in biomedical modeling.

Our main goal in this manuscript is to estimate the extent of the influence that unmet ADL needs have on the amount of incurred expenditures within the follow-up interval. We also attempted to estimate the extent of this influence in the most unbiased way possible. The model is constructed to estimate net expenditures and also provide unbiased estimates of parameters (see the Table (4)), especially the one that reflects the influence that unmet ADL needs have on the probability of incurring expenditures as well as on the amount of expenditures. We constructed a model that allowed estimation of the influence of unmet ADL needs in the context of known contributors to older adults medical expenditures. These known contributors also include respondents baseline characteristics (e.g. age, ADL status). A number of subjects have died during the period of study and we incorporate the knowledge of their death in the form of Inverse Probability Weighting (IPW) procedure. Since observations missing due to death are clearly not Missing Completely At Random (MCAR), this information has to be incorporated to avoid possibly biased estimates of the cost [28]. Such a bias is typically a serious problem whenever the complete multivariate distribution of the data cannot be fully specified. It is possible to think of such a bias as resulting from two sources: one is the lack of information about the unmeasured disease severity between patients with different levels of ADL and the other is due to unbalanced nature of the data when some of the patients die during the period of study. Moreover, [29] (see p. 490) noted earlier that, ignoring the data that are not missing completely at random, such as MAR (Missing At Random), “…can potentially introduce bias in the estimates of regression parameters”. Due to these two concerns, we are introducing inverse probability weighting into our model to account for the missingness pattern due to the dropout.

Our secondary goal in conducting this research is to propose a longitudinal model that can relate highly skewed medical expenditures data with substantial missing data proportion to unmet need for ADL and some additional covariates. The ultimate hope is that the conclusions obtained will be of use in public policy. The current manuscript is structured as follows. Section 2 introduces the population averaged model we use to describe total medical expenditures across all categories and a novel method we introduce in order to fit it. Section 3 is dedicated to illustrating how the model works with simulated data. Section 4 shows how the model performs with real data. Finally, Section 5 describes possible directions of future research.

2 Model

Suppose a group of n subjects is followed with the medical expenditures being observed on a monthly basis. For each subject i, i = 1, …, n, we have m_i monthly expenditure recorded where m_i ≤ m = 13 observations. The reason there may be less than m observations for some subjects is because some subjects die before the end of study. The total number of observations is $N={\sum }_{i=1}^n{m}_i$. The observed expenditures of each subject are Y_i = (Y_i1, …, Y_i,mi)′ with Y_ij being expenditures in jth month for the ith subject. Also, denote $\mathbf{Y}={({\mathbf{Y}}_1^{\prime },\dots ,{\mathbf{Y}}_n^{\prime })}^{\prime }$. Of course, Y_ij ≥ 0; moreover, for some of i = 1, …, n, j = 1, …, m_i, Y_ij = 0; that happens when a subject doesn’t incur expenditures in any of the categories. All the subjects that did not die during the period of study are not observed any longer after Dec 31st, 2005 which is a fixed date that does not depend on a subject and also does not depend on the medical costs. In other words, the experiment ends at a pre-specified time. Also note that the time a subject enters the study is not fixed to be Sept 1st, 2004 but can occur on any day after that date and until the end of the study period. Technically, the time of entry into the study for a subject is his/her initial interview day. However, we will consider the 1st day of the month that follows the interview month as the true time of entry since this is more in line with how the Medicare claims are processed. In such a setting, the medical expenditures data should be viewed as truncated on the left.

Some earlier work, e.g. [10], has only modeled the dependence of just one category - hospital expenditures - on a number of subject specific covariates. Our interest lies in modeling the total expenditures and, in particular, the influence of the number of ADL disabilities and unmet need for help with ADL disability on these expenditures. All of the covariates involved in this study were constant over time and so we denote X_i the vector of covariates for ith subject. We follow the approach similar to that used by [26] and [10] with several notable differences. First, we model jointly the probability of incurring total medical expenditures (that is, P(Y > 0)), and the amount of positive expenditures. Second, instead of the subject level approach used in [26] and [10], we use the population-averaged approach. In other words, instead of using subject-wise random effects to induce the necessary autocorrelation structure, the autocorrelation matrix for each subject is modeled directly. Also, unlike [26] and [10], we are using the GMM (Generalized Method of Moments) approach that is a generalization of the classical GEE (Generalized Estimating Equations) approach first pioneered in [30]. GMM was first proposed in econometric context in [31]. The classical GEE approach is very robust to misspecification of the subject-wise autocorrelation structure and can easily handle unbalanced designs. Also, the GEE choice is rather sensible when the number of subjects is large relative to the length of time period involved. However, GEE is known to be consistent only when the missing data are assumed to be MCAR which is clearly not the case for us. Moreover, our data are generated from what is effectively a mixture model that does not belong to an exponential family. The last two issues force us to consider a somewhat more general GMM which has been used extensively in econometrics for a long time.

Note also that our data have been obtained using a survey (NLTCS) with complicated weigtht structure, and, moreover, the death of certain subjects produces the dropout effect. There are a number of ways to compensate for the dropout effect. One of them is to use only the so-called “complete cases” (subjects that didn’t die until the end of the study period) with appropriate weights so that they account for “incomplete cases” as well. This approach is rather inefficient as it forces the researcher not to use a substantial proportion of the data. Therefore, we prefer to be able to weight the contribution of each subject at each month of observation to the total pseudo-likelihood of the model explicitly. This can be achieved using the Inverse Probability Weighting-GMM (IPW-GMM) approach; for detailed historical introduction see e.g. [29], Chapter 18. One of the few references in the literature using the GEE approach for data with excessive zeros is [32]; they consider just clustered (rather than specifically longitudinal) data that are not weighted in any way. [32] also uses an independent working correlation assumption and only takes into account the intrasubject correlation when computing “sandwich” estimators of standard errors of estimated parameters.

We assume that Y_ij may be equal to zero with a non-zero probability that depends on the covariate vector X_i = (X_i1, …, X_ip)′. Let γ be the p × 1 vector of parameters and denote p(X_i, γ) the probability of expenditures Y_ij being equal to zero. A number of choices are available to model this probability, e.g. the logistic or a probit model. To make our discussion simpler, we use a logistic model whereby

which implies that $\text{logit}P(Y_{ij}>0\mid {\mathbf{X}}_i)={\mathbf{X}}_i^{\prime }\gamma$. Since the positive expenditures are highly skewed to the right, a right skewed distribution should be used to model it. A gamma family represents a convenient choice, including exponential and χ² distributions as special cases. We parameterize gamma density as $f(y;\lambda ,\nu )=\frac{\lambda }{\mathrm{\Gamma }(\nu )}{(\lambda y)}^{\nu -1}{e}^{-\lambda y}$ defined for y ≥ 0. This parameterization (see, for example, [33]) is commonly used for purposes of generalized linear modeling. It is rather convenient since the mean is, then, $\frac{\nu }{\lambda }=\mu$ and the variance is $\frac{\nu }{{\lambda }^2}={\mu }^2\frac{1}{\nu }$ with $\frac{1}{\nu }$ being the dispersion parameter. There is some empirical evidence from various applications that the gamma distribution with the constant dispersion parameter may be insufficiently heavy tailed for healthcare modeling purposes; thus, we assume that, in general, the dispersion parameter may depend on a set of covariates through the log link. We use the log link to model the mean expenditure of ith subject in jth month as well. For simplicity, we assume that the covariate vector used to model the mean and the dispersion parameter is the same as that in the logistic model above though this need not always be true. Thus, letting δ and ρ be respective parameter vectors and $\eta =\frac{1}{\nu }$ the dispersion parameter, we define $log{\mu }_{ij}={\mathbf{X}}_i^{\prime }\delta$ and $log{\eta }_{ij}={\mathbf{X}}_i^{\prime }\rho$.

For an observation Y_ij, denote ξ_ij = 1 if Y_ij > 0 and ξ_ij = 0 if Y_ij = 0. Denote ${\text{logit}}^{-1}(x)=\frac{e^x}{1+e^x}$. Under the assumption of independence working model, we can write the full pseudolikelihood of our model as

where

The above is only true if there are no missing observations and all observations come from a simple random sample. In practice, some subjects die before the end of the study period and this has to be reflected in a weight assigned to each observation. More specifically, we suggest treating each death as a terminal event and use the inverse probability weighting (IPW) procedure to take that into account. For ith subject, define R_ij = 1 if the ith subject has a recorded expenditure amount (whether 0 or positive) in jth month and 0 otherwise. Denote the vector of response indicators for ith subject R_i = (R_i1, …, R_im)′. Then, the occasion when the subject experiences the terminal event (death) is $D_i=1+{\sum }_{i=1}^m{R}_i$. For a complete case, where m_i = m and the entire expenditure vector Y_i = (Y_i1, …, Y_im)′ is observed, we have D_i = m + 1. For an individual with an incomplete vector of m_i < m responses, we only observe Y_i = (Y_i1, …, Y_imi)′ and D_i = m_i + 1.

We assume that the dropout due to death in jth month can be thought of as occuring “at random” in the classical sense - the probability of death only depends on the subject specific covariates as well as prior expenditures (observations) in months up to j − 1st. In any month j, the probability of survival through the jth month for the ith subject is p_ij = P(D_i > j|D_i ≥ j) = P(R_ij = 1|R_i1 = ··· = R_i,j−1 = 1). It is commonly assumed that R_i1 = 1 for any subject i and, therefore, p_i1 = 1. Due to the MAR (Missing At Random) assumption, it seems reasonable to assume that the probability of survival for ith subject in jth month depends on both subject-wise covariate vector X_i and on the expenditures in all of the prior months. Let Y_ij = (Y_i1, …, Y_i,j−1)′. Define the joint covariate vector ${\mathbf{Z}}_{ij}={({\mathbf{Y}}_{ij}^{\prime },{\mathbf{X}}_i^{\prime })}^{\prime }$- a vector that typically consists of all observed expenditures prior to jth month for ith subject as well as the subject specific covariates X_i. The probability p_ij can then be modeled using, for example, logistic regression, as

In order to compute the weight assigned to subject i in jth month, we take the inverse of unconditional probability of survival occurring in jth month which is ${\prod }_{k=1}^{j-1}{p}_{ik}$ and multiply it by the sample weight under the independence assumption. If we denote the sampling weight of the ith subject ω_i, the weights to be used are defined as ${\omega }_{ij}={\omega }_i{\left[{\prod }_{k=1}^{j-1}{p}_{ik}\right]}^{-1}$.

Thus, if the ith individual had observed expenditures Y_ij in the jth month (whether zero or positive), that subject will receive weight ω_ij in that month. If, on the contrary, the subject has died before or in jth month, he will receive the weight of zero. For the ease of notation, let us introduce the joint parameter vector θ = (γ,δ,ρ). Then, the log-pseudolikelihood will become

where each

Note that, in our case, the marginal distribution of each observations Y_ij is a mixture that does not belong to an exponential family. Traditional GEE generally assumes that marginal distribution of observations Y_ij belongs to an exponential family; however, a slightly more general GMM (Generalized Method of Moments) (see, e.g., [31]) approach allows us to dispose of this assumption. GMM - based approach allows us, under weak assumptions on the true dependence within the sample, to obtain consistent and asymptotically normal estimators of the model parameters. As in the classical GEE method, this is done without taking into account the unknown dependence structure when formulating estimation equations.

Let the score function be $U(\mathbf{Y}\mid \theta )=\frac{\partial l(\mathbf{Y}\mid \theta )}{\partial \theta }={\sum }_{i=1}^n{U}_i({\mathbf{Y}}_i\mid \theta )$ where $U_i({\mathbf{Y}}_i\mid \theta )=\frac{\partial }{\partial \theta }{l}_i({\mathbf{Y}}_i\mid \theta )$. Then, the sequence of estimating equations is defined as

At this point, we need to invoke some results on convergence of GMM estimators in cases where the sample is not the random draw from the population of interest. For example, [34] showed that the asymptotic normality is still valid for a variety of situations where weights are used as propensity scores and are dependent both on covariates and on other observations (which corresponds to the missing at random case). By this result, the normalized difference n^−1/2(θ̂−θ) converges to a normal distribution with mean zero and variance Σ_θ that needs to be estimated. The usual sandwich estimator has to be modified due to the presence of estimated weights. The usual “sandwich” estimator that would be ordinarily used to estimate Σ_θ is

In order to estimate the weights p_ij (and, therefore, ω_ij), we need to run a logistic regression analysis of the “stacked” dataset that includes observations for all subjects. This involves solving the following system of equations:

The system of equations above represents the sum of subject-wise score functions for the logistic regression model (2.2) used to compute propensity weights. Note that only when R_i,j−1 = 1 the ith subject’s contribution is not equal to zero. In other words, only when in a given month j − 1 the ith subject had recorded (possibly zero) expenditures, this subject contributes a nonzero term to its score function S_i(α). One can also say here that, after the death of ith subject he/she does not contribute anything to its score S_i(α). In order to adjust for estimated weights, we need to change the sandwich estimator. More specifically, we replace U_i(Y_i|θ̂) in the middle part of the “sandwich” estimator with the residual of the multivariate regression of U_i(Y_i|θ̂) on S_i(α). Such a residual is

Taking the above into account, the modified estimator of Σ̂ becomes

This amounts to effectively using a residual from a multivariate regression of U_i(Y_i|θ) on S_i(α) in order to reduce variability of estimated θ̂. For details of this approach, see e.g. [29], Chapter 18.

3 Simulation study

In order to assess model performance a simulation study was performed. Data characteristics of interest are sampling weights, positive right skewed response with a point mass on zero, correlated response over time, dropout due to death, and the staggered entry into the study due to differing interview dates.

In order to mimic the analysis of real data, 5 independent binary covariates X_ik, k = 1, …, 5, were generated from a Bernoulli distribution with the probability p = 0.5. Out of many ways of generating multivariate distribution with gamma marginals, we select the so-called Clayton copula. In brief, for a random vector X = (X₁, …, X_d)′ with continuous marginal distributions F_i(x) = P(X_i ≤ x), applying a probability integral transform results in a random vector U = (U₁, …, U_d)′ = (F₁(X₁), …, F_d(X_d))′ with uniform marginals. Then, the joint cumulative distribution of (U₁, …, U_d) is said to be a copula. The Clayton copula belongs to a specific class of copulas commonly called Archimedean copulas that allow modeling of dependence for an arbitrary high number of dimensions using only one parameter. For more details and an introduction to copulas, see e.g. [35]. In our cases, expenditures for a finite population of 50, 000 individuals were generated from a 13 dimensional Clayton Copula with the parameter θ = 2. In order to mimic expenditures Y_ij, we back transform individual univariate uniform random variables using the appropriate inverse gamma CDF and then multiply by the necessary scale factor. Let Γ⁻¹(u; α) be the value of inverse Gamma CDF at u with the shape parameter α. Now, the resulting “synthetic” expenditures are correlated with the correlation coefficient ρ ≈ 0.5 and are represented as

where u_ij is the random variate from the copula corresponding to X_ij, Y_ij is the expenditure of ith subject in jth month, and ν is the dispersion parameter,

In order to introduce a point mass on zero, the probability of no expenditure was calculated for each month of each subject. An indicator of non-zero expenditure was generated from a Bernoulli distribution with probability $p_{ij}^{\mathit{exp}}$ that was obtained as

We used the same set of covariates here as the one used to model the distribution of expenditures Y_ij. If the indicator was one, the expenditure was left as it was, if the indicator was zero, the expense was set to zero, to simulate months with no expenditures. The introduction of zeros reduced correlation to levels consistent with the observed data (ρ ≈ 0.3).

In order to introduce staggered entry due to interview date, a random start month was drawn from a multinomial distribution with π = (0.1810, 0.3101, 0.3142, 0.1634, 0.0294, 0.0011) for months 1 to 6 to mimic the sample exactly. Unconditional probabilities of survival are calculated based on a logistic model that depends on the same 5 binary covariates as before and the expenditure of the previous month. We assume the dependence on the expenditure of ith subject in the previous month in order to simulate the MAR (Missing At Random) compliant dropout times due to death. Unconditional probabilities of not dropping out in the jth month for the ith subject are computed as a product of the probability of not dropping out j − 1st month multiplied by the probability of survival through jth month. It is assumed that the probability of survival in the first month is 1. Therefore, the probability that ith subject dies in jth month is given by the expression

Beginning with the month following the start month, a binary death indicator was sampled from a Bernoulli distribution with $P=p_{ij}^{\mathit{death}}$. For each subject, expenditures prior to the start month and after the month of death (if any) were dropped as being unobserved and not occurring respectively. In order to create sampling weights, data was structured to have 100 strata with 10 clusters each containing 50 subjects. A two stage sampling design was implemented. The first stage consisted of selecting a simple random sample of two clusters per each stratum. The probability of selection for each cluster was, therefore, equal to 0.2. At the second stage, a sample of four subjects per cluster was chosen with probability proportional to size so that we ended up with a sample of size 800. The size measure for ith subject was defined similarly to [10], namely

where δ_i is an indicator variable equal to one if the i^th individual died before the end of the study and zero otherwise; ȳ_i is the average monthly cost of ith subject. Note that this is an informative sampling scheme where the subjects who died during the period of study are oversampled. Sampling weights were calculated as the inverse of probability of selection. This was repeated 1000 times to create 1000 samples with 800 subjects each.

For each sample, a logistic model was fitted to estimate death weights, and then the IPW-GEE model was fitted and the Robust Sandwich Variance Estimator was used for standard errors. Coefficient estimates, variability of those estimates, mean of the Robust standard errors, and coverage probability of true parameter values are presented in Table (1). Note that in 3.3% of simulations, estimation procedure did not converge and the respective results were not included in the Table (1). It is clear that our approach seems to work rather well with simulated data in this informative sampling scheme. Note, in particular, an excellent coverage probability for all of the confidence intervals involved.

Since we estimate monthly expenditures in this setting, the presence of the dropout related weights does not change the values of estimated parameters greatly. We conducted a separate simulation where the weights used were just sampling weights, that is, ω_ij = ω_i for each subject i. The resulting Table is almost identical to the Table (1) and is omitted in the interest of brevity. What is different, however, is the numerical stability of the estimation procedure. Whereas before, as we mentioned earlier, the estimation procedure did not converge in only about 3.3% of all cases, when the dropout related weights are omitted, this number rises to approximately 15%. We believe this is yet another reason to include the dropout related weights in our estimation procedure.

For comparison purposes, we also fit the same model to the data with no use of weights; therefore, the presence of missing data due to death as well as sampling weights are both ignored. The results are given in the Table (2). Note that the fitted coefficients are now biased, with the particularly pronounced bias present in estimated intercept of the medical expenditure model and in the dispersion parameter of the gamma distribution. The informative sampling scheme we are using results in oversampling of subjects with higher average expenditures and those who died during the observation period. If this element of the sampling design is ignored, one would expect a larger intercept in the medical expenditure model as well as a larger estimated dispersion parameter. The latter would account for increased variability in medical expenditures that was not properly accounted for by the sampling design.

4 Analysis of the NLTCS data

Included in the analysis are the 2400 respondents of the 2004 NLTCS community survey of whom 467 died during the observation period. Each respondent was followed for up to 13 months after the interview. Observation stopped with either death or being censored at the end of 2005. On average, respondents were followed for 10.66 months with a standard deviation of 2.46 months. Monthly medical expenditures in the follow-up period were obtained from the Medicare claims data by accumulating recorded expenditures within the same month. 86.12% of respondents had a nonzero total medical expenditure during the follow-up period.

A descriptive summary of monthly expenditures and some weighted sample characteristics is given in the Table (3). Ordinarily, SAS PROC GENMOD is used for the GEE approach to estimation of longitudinal data-based models. However, PROC GENMOD requires an explicit distributional assumption out of the (fairly short) list; our data have been generated by a mixture of right-skewed continuous gamma distribution and the degenerate point-mass distribution. Due to that, we decided to use PROC NLMIXED instead since it can approximately solve the necessary non-linear estimating equation.

For comparison purposes, we fit gamma models with both constant and non-constant dispersion parameters. We examine goodness of fit of these distributions to the positive monthly healthcare expenditures. More specifically, 1000 data sets are simulated for each model; afterwards, the expected quantiles based on the model are computed and plotted against sample quantiles. These plots show that, in both cases, the fit is relatively good although the heteroscedastic version does offer certain improvement by providing a more heavy tailed distribution. We also tried alternative heavy tailed distributions, such as inverse gamma and Weibull. Respective quantile to quantile plots showed a much worse pattern and are not shown here in the interests of brevity.

The weights have been fitted using the same set of covariates as the main model with the addition of the prior months medical expenses. That set of covariates Z_ij consists of moderate to severe ADL disability (3–5 categories), unmet need for ADL disability, age (viewed as a continuous varibale), gender, race (white or nonwhite), and medical expenses in j − 1st month.

The results of our analysis are presented in (4) and (5). The first two parts of these models are very similar. All of the parameters that are statistically significant in one of them are significant in the other and vice versa. For convenience, we will refer to p-values from (4). They suggest that both increased ADL and the unmet need for medical assistance are strongly associated with the higher amount of incurred expenditures (p < 0.0001 and p < 0.0008, respectively). The increased ADL is also strongly predictive of the higher probability of incurring expenditures (p < 0.0072); however, the unmet need is not strongly predictive of that same probability (p < 0.2473). Both gender and race also seem to be strongly associated with the higher probability of incurring medical expenditures; females seem to incur costs less often than males (p < 0.0001) and whites seem to incur expenditures more often than people of other races (p < 0.0002). Finally, the older age is strongly associated with larger medical expenditures (p < 0.0007).

5 Discussion

Our work continues the recently observed trend of utilizing complex survey data in statistical healthcare research. A particularly beneficial feature of a survey such as NLTCS is its linkage to Medicaid/Medicare claims which provides an opportunity to study the relationship between demographic characteristics and medical expenditures. Modeling medical expenditures has long been known to be difficult due to the presence of a large number of zeroes as well as the highly skewed nature of the non-zero part of observations and commonly present missing observations. All of these features need to be taken into account when analyzing the data.

A distinctive feature of our analysis is the use of a marginal model that is estimated using the GEE/GMM approach with inverse probability weighting. This approach has been less popular in the literature than the competing mixed modeling approach; [32] is one of the few available examples but they only modeled general clustered (and not specifically longitudinal) data that didn’t come from a survey and were not subject to a dropout phenomenon. We found the GEE/GMM based approach particularly applicable in our case since our dataset had a large number of participants relative to the length of observations available; this latter fact makes GEE/GMM asymptotic results practically applicable. Its other appealing properties is robustness to misspecification of the data covariance structure as well as the ease of adapting it to the presence of dropout in the data under the MAR (Missing At Random) Assumption.

It is also of interest to note that, due to the “working independence” assumption that we have, in fact, been using in our GEE/GMM analysis, we can utilize any of the subject-wise covariates X_i, even though some of them may have missing values (see, e.g. the discussion in [29] on p. 529. Unlike in a classical IPW method that is only using the so-called “complete cases” (only subjects who didn’t die until the end of the study period), the IPW-GEE/GMM estimator that we are using in this research, does not disregard any of the available data.

The choice of a particular distributional model for the data analysis in this situation is rather complicated. A large number of different distributions have been used in practice before, e.g. the generalized gamma distribution [10], Pareto distribution [36] and several others. We conjecture that a combination of excessive zeros and a right-skewed heavy-tailed distribution may also be modeled as a member of a general exponential dispersion family, e.g. one of the so-called Tweedie family distributions (see e.g. [37]). Such an approach would avoid the necessity of using a mixture model to describe the data and may result in a simplified estimation procedure. Future research is needed to investigate such a possibility.