Biostatistics Advance Access originally published online on May 25, 2005
Biostatistics 2006 7(1):1-15; doi:10.1093/biostatistics/kxi036
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Longitudinal profiling of health care units based on continuous and discrete patient outcomes
Department of Statistics, University of Florida, Gainesville, FL 32611, USA mdaniels{at}stat.ufl.edu
Department of Health Care Policy, Harvard Medical School, Boston, MA 02115, USA and Department of Biostatistics, Harvard School of Public Health, Boston, MA 02115, USA
* To whom correspondence should be addressed.
 |
SUMMARY |
|---|
 TOP SUMMARY 1. INTRODUCTION2. JOINT MODELS FOR...3. INFERENCE4. APPLICATION TO THE...5. CONCLUDING REMARKSREFERENCES |
|---|
Monitoring health care quality involves combining continuous and discrete outcomes measured on subjects across health care units over time. This article describes a Bayesian approach to jointly modeling multilevel multidimensional continuous and discrete outcomes with serial dependence. The overall goal is to characterize trajectories of traits of each unit. Underlying normal regression models for each outcome are used and dependence among different outcomes is induced through latent variables. Serial dependence is accommodated through modeling the pairwise correlations of the latent variables. Methods are illustrated to assess trends in quality of health care units using continuous and discrete outcomes from a sample of adult veterans discharged from 1 of 22 Veterans Integrated Service Networks with a psychiatric diagnosis between 1993 and 1998.
Keywords: Bayesian hierarchical model; Correlation matrix; Informative priors; Latent variable; Mental health
|
1. INTRODUCTION |
|---|
|
TOPSUMMARY1. INTRODUCTION2. JOINT MODELS FOR...3. INFERENCE4. APPLICATION TO THE...5. CONCLUDING REMARKSREFERENCES |
|---|
Profiling health care providers involves combining information collected on subjects treated by health care providers and making inferences about the quality of care associated with each provider. The collected information has typically been one of three types: cross-sectional univariate outcomes (Goldstein and Spiegelhalter, 1996
; Normand et al., 1997); longitudinal univariate outcomes (Aguilar and West, 1998; Bronskill et al., 2002); or cross-sectional multivariate outcomes (Landrum et al., 2003). In all cases, the outcomes are clustered within health care units. It has become increasingly common to collect both continuous and discrete outcomes longitudinally on providers in order to infer trends in quality of care.
The Veterans Health Administration (VHA) is a case in point. Some have argued that care in the VHA, the largest vertically integrated health care system in the United States, is poorer than in non-VHA institutions (United States Senate, 1999). In 1995, the VHA initiated the National Mental Health Performance Monitoring System and began to publish systematic performance data for the 22 geographically defined service networks that provide mental health services to its veterans. Table 1 summarizes outcomes over a 6-year period for all adult veterans who had an inpatient admission in a general psychiatric program located in Veterans Affairs Medical Centers during the first 6 months of each fiscal year. Large values of the inpatient measures, with the exception of days to readmission, and small values of the outpatient measures, with the exception of days to first visit, indicate poor quality of care (Rosenheck and DeLilla, 1999). The goal for the VHA is to characterize trends in quality of individual service networks on the basis of the outcomes and to identify problems in quality for specific networks.
|
The observed data consist of multilevel continuous and discrete outcomes collected repeatedly within clusters over time. Heterogeneity is introduced from several sources: the correlation among the multivariate outcomes on the same patient, correlation among the multivariate outcomes on the same health care unit, the repeated measurement of health units over time, and the heterogeneity across health units. Several methods have been proposed to handle joint modeling of multilevel continuous and discrete outcomes in the cross-sectional setting (Dunson, 2000
; Dunson et al., 2003; Gueorguieva and Agresti, 2001; Landrum et al., 2003; Lee and Shi, 2001). The key idea is the introduction of a set of latent variables that reduce the dimensionality of the problem for outcomes having different measurement scales. Fewer methods are available for repeated measurements of multilevel discrete and continuous outcomes. Aguilar and West (1998) proposed a multivariate random-effects time series model for repeated binary outcomes through an auto-regressive structure for the cluster random effects. However, this model is restricted to joint modeling of binary outcomes. Outside the context of multilevel data, methodology for longitudinal data for the the latent variable (Dunson, 2003; Roy and Lin, 2000; Oort, 2001), latent class (Lin et al., 2002; Duncan et al., 2002; Muthén et al., 2001; Douglas et al., 1999; Douglas, 1999), and structural equation model (Lee et al., 1992; Mueller, 1996) settings have been developed.
In this article, we consider an analysis of I units measured repeatedly over T time points on the basis of K -dimensional response vectors of continuous and discrete data from a total of N individuals. A model-based approach that distinguishes these features has the advantage of providing estimates that filter the sampling and measurement noise from the underlying quality signals, thus providing more reliable estimates of quality (McClellan and Staiger, 1999). We propose multilevel multidimensional latent variable models for the multivariate data that permit inclusion of covariates and serial dependence among the latent variables. Underlying normal regression models for each outcome are used and dependence among different outcomes is induced through latent variables. Serial dependence is accommodated through modeling pairwise correlation of the latent variables. We use Bayesian methods for inference. We describe the substantive problem more completely in Section 2 as well as introduce the model and discuss identifiability issues. Section 3 describes our approach to inference, including the prior distributions and posterior computations. In Section 4, we apply our model to the VHA data and make concluding remarks in Section 5.
|
2. JOINT MODELS FOR LONGITUDINAL MULTILEVEL CONTINUOUS AND DISCRETE DATA |
|---|
|
TOPSUMMARY1. INTRODUCTION2. JOINT MODELS FOR...3. INFERENCE4. APPLICATION TO THE...5. CONCLUDING REMARKSREFERENCES |
|---|
The VHA provides health care to more than 4 million veterans annually and is the nation's largest provider of behavioral health services. Like the broader U.S. mental health care system, in the mid-1990s the VHA undertook several initiatives to reduce costs and to increase the number of veterans served (Kizer et al., 2000). In particular, the agency adopted a budget allocation system that gave service networks strong incentives to serve more patients and to reduce per patient costs. Typically, inpatient programs lose resources as a consequence of such health care reorganizations, as there is a shift in the type of treatment provided. A study by Chen et al. (2003) found that the shift from inpatient to outpatient mental health care in the VHA over the period 1995–2001 resulted in a 21% decrease in inpatient spending and a 63% increase in outpatient spending.
With such dramatic shifts in care settings, a natural concern relates to whether the quality of mental health services has declined. On the medical side, a study of the impact of VHA reorganization on the quality of care in long-term care units found an increase in pressure ulcer development (Berlowitz et al., 2001). Because regional service networks have strong financial incentives and large administrative control, the need to carefully monitor the quality of care in programs that may be losing resources as a result of the reorganization is especially elevated in a population of patients who may have difficulty navigating the health system.
Let Yitjk denote the value of the kth response for the jth subject in the ith unit (network) at time t. Subjects are clustered within units, and units are measured repeatedly.
Assume for binary-valued variables Zitjk
N(
itjk, 1) such that Yitjk = I(Zitjk > 0), where
 | (2.1) |
 | (2.2) |
Here, 
itj is a subject-specific random effect, 
itk, a network random effect, and xpitj (p = 1, ..., P), a subject-specific covariate. Let k = 1 and k = 5 correspond to the binary random variable for the inpatient and outpatient measures, respectively. The itj appears in the regression for the two binary responses for each subject inducing a patient-level correlation between the two responses. We specify the distribution of the itk in Section 2.3.
The continuous variables, Y, are grouped into three inpatient (k = 2, 3, 4) and three outpatient (k = 6, 7, 8) measures, and trivariate normal distributions are assumed for each. The two sets of continuous variables are conditional on a positive value for the corresponding binary-valued variables, that is, having an inpatient event and an outpatient event. Let U
{in, out} with U = in corresponding to k = 2, 3, 4 and U = out corresponding to k = 6, 7, 8. The model for continuous outcomes is
 | (2.3) |
 | (2.4) |
where 
denotes the covariance of the continuous variables at time t. See Table 1 for an ordered list of the eight patient-level measures.
The K network effects are assumed to be represented by L latent variables that are a priori independent,
 | (2.5) |
 | (2.6) |
where 
k = (1k, ..., Lk)' is an L x 1 vector of discrimination parameters assumed constant over time. The latent variables, 
ilt, represent unobserved traits of an individual unit, such as inpatient and outpatient quality, at time t. The variance of the unit random effects, 
represents heterogeneity among outcomes in a given year not explained by the L latent variables. Rl denotes a T x T correlation matrix that permits serial dependence in the L latent variables.
We model the within-unit dependency of the outcomes measured at different times through the correlation matrix, Rl. Our motivation for specifying some structure on Rl is two-fold: to predict latent quality at future times and to increase power through estimation of fewer correlation parameters. The latter consideration is important in our example as there are only 22 networks for estimating a six-dimensional correlation matrix. The complications in developing models for a correlation matrix arise for three reasons: (1) the marginal variances are fixed at one; (2) the matrix needs to be positive definite; and (3) in our setting, the parameters should be interpretable in a longitudinal context.
We consider first-order Markov models. Recall that il NT(0, Rl), which implies Var(ilt) = 1, t = 1, ..., T. For a first-order Markov model, we have
 |
Using iterated expectations, it is easy to show that Var (ilt) = 1 implies that 
for l = 1, ..., L and t = 1, ..., T, where 
l0 = 0.
Models for l,t–1 can be developed subject to the constraint that l,t–1 lies in the open interval (–1, 1). In particular, we consider
 | (2.7) |
where z(·) is Fisher's z-transformation and gl(t – 1; 
) is a smooth function of time parameterized by for each of the L latent variables; this is a nonstationary, first-order Markov model. The form of these nonstationary models facilitates prediction by providing a mechanism to extrapolate the dependence structure past T. That is, we can use gl(·) to fill in l,t–1 for t > T. A simplified stationary version of this model involves setting gl(t – 1; ) = z(l) so that l,t–1 = l which implies 
The correlation between observations one time unit apart in this model is l. Markov models of higher order involve severely constrained parameters due to the restrictions on the marginal variance and are difficult to apply in this setting. The models proposed here are similar to the structured antedependence models proposed in Zimmerman and Nunez-Anton (1997) for a covariance matrix.
To examine identifiability of parameters of the between-unit model, we approximate the marginal posterior distribution of itk given only the within-unit model, defined in (2.1–2.4), with
 | (2.8) |
where Vi is assumed known and i = (i11, i12, ..., i1T, ..., iKT)'. The priors on i are given by (2.5)–(2.6), the between-unit model. Integrating over the latent variables, ilt, we obtain
 | (2.9) |
where A is a KT x KT matrix with K(T x T) blocks on the main diagonal, the kth of which has the form 
and with off-diagonal blocks, Ojk, of the form, 
 | (2.10) |
where i·t = (i1t, ..., iKt)', the kth row of K x L matrix 
is (k1, ..., kL)', and 
·t = diag(
1t, ..., Kt). This resembles a standard factor analytic model. Because Vi is known, and ·t will be identified under standard conditions for factor analysis models. We now need to show that Rl is identified. Consider,
 | (2.11) |
where ik· = (ik1, ..., ikT)'. As the components of kl and k· = diag(k1, ..., kT) are identified from (2.10), Rl will be identified from the correlation among ik over time.
We consider several quantities for evaluating networks longitudinally using the latent variables. Assuming small values of ilt correspond to poorer quality of care, networks worsening over time may be identified through
 | (2.12) |
for each i and l. Here, Tq(lt) is the qth quantile of the distribution of {ilt; i = 1, ..., I} and Y denotes the observed data. Similarly, units improving over time may be identified through calculation of
 | (2.13) |
The trend may be characterized more formally by estimating and ranking the slope, 
For example, we can compute
 | (2.14) |
where 
and X is a 2 x T matrix with tth row, (1, t), to provide inference about the rate of change. Alternatively, Spearman's correlation coefficient may be estimated as a nonparametric alternative.
Changes in the relationship between the latent variables may be assessed through the correlation, e.g.
 | (2.15) |
for l 
l'. An overall summary may be represented as 
for some selected weights, wl.
|
3. INFERENCE |
|---|
|
TOPSUMMARY1. INTRODUCTION2. JOINT MODELS FOR...3. INFERENCE4. APPLICATION TO THE...5. CONCLUDING REMARKSREFERENCES |
|---|
We adopt a Bayesian approach to inference and specify prior distributions for the hyperparameters. We assume the presence of two latent variables, one representing the quality of inpatient care and one the quality of outpatient care. This assumption was based on a variety of considerations. First, several studies have shown that the type of health care reorganization operating at the time of our study has generally resulted in different quality changes in the inpatient and outpatient settings. This shifting of inpatient to outpatient treatment settings has been observed in managed care organizations (Mark and Coffey, 2000). Second, subjective measures of quality of mental health services in the VHA, such as patient satisfaction, have demonstrated differences in inpatient and outpatient quality for some subgroups (Hoff et al., 1998). Third, in a prior cross-sectional analysis of the 1998 data, we found empirical evidence of two latent variables (Landrum et al., 2003).
We place informative priors on the discrimination parameters N(
, c2ILK), where
 | (3.1) |
We discuss the choices of 
> 0 and c2 in Section 4.2 and also examine sensitivity to these choices. This is the prior for the discrimination parameters in a two-latent variable model (L = 2) and is consistent with an inpatient and an outpatient quality variable. The 1's and –1's in the prior specification correspond to that response's (of the K) contribution to the inpatient and outpatient latent variable. For example, favorable inpatient characteristics include many days to first readmission but few number of bed days; in (3.1), these get the multipliers 1 and –1, respectively. The components of the prior means for were selected in order to place equal weights on the inpatient and outpatient components as this was deemed most reasonable by the investigator who collected the data. Informative priors on the discrimination parameters or restrictions on the form of the discrimination matrix, , are necessary for identifiability of the latent variables (Lopes and West, 2004).
Diffuse proper priors are specified for the regression coefficients, ßtk, i.e. normal priors with mean 0 and variance 10 000 and gamma priors on the inverse of the variance components, 
and 
with parameters (0.001, 0.001). A Wishart prior is placed on 
with degrees of freedom equal to 
and scale matrix equal to the identity matrix. Because the continuous responses were standardized, marginal variances of one for the scale matrix in the Wishart prior seemed reasonable. For the unstructured correlation matrix, Rl, similar to Barnard et al. (2000) and Daniels (2005), we specify a uniform prior over the compact subspace of the T(T – 1)/2 dimensional cubic, [–1, 1]T(T–1)/2, such that Rl is positive definite. A uniform prior on the interval (–1, 1) is specified for the correlation parameter l in the Markov model given in Section 2.3.
A Gibbs sampler was designed to sample from the posterior distribution of the parameters. The full conditional distributions are all known forms (gamma distributions for the inverse of the variance components, normal distributions for the regression parameters, latent variables, and discrimination parameters, and Wishart distributions for the covariance matrices) except for the correlation matrices, Rl, l = 1, ..., L. We implement an adaptation of the component-wise approach of Barnard et al. (2000) to sample Rl when Rl is assumed to be unstructured. In the case of Markov dependence, we use a random walk Metropolis–Hastings algorithm to sample the correlation parameter.
We sampled the parameters in blocks to increase the efficiency of the sampling algorithm. In particular, we sampled the following parameters as blocks in the within-network model: the P x 1 vectors ßtk, for k = 1 and k = 5; the (3P) x 1 vectors of (ßt2, ßt3, ßt4) and (ßt6, ßt7, ßt8); the matrices for U and t; and in the between-network model, the L x 1 vector of k for k = 1, ..., K, and the T x 1 vectors il for l = 1, ..., L. The other parameters were sampled component-wise.
We use the deviance information criterion (DIC) of Spiegelhalter et al. (2001)
as an overall measure of model fit. Because the main unit of interest is the network, we utilize the following integrated likelihood as the basis for the deviance in the DIC
 |
where yitj is the K x 1 dimensional response vector and it is a K x 1 dimensional latent variable representing K outcomes at the network level for each time, t. Because this is not available in closed form and the number of patients per network is quite large, we approximate the integrated likelihood using
|
| (3.2) |
where 
i is the maximum likelihood estimate (MLE) of i = {itk} (here, we use posterior mean) and Vi is the estimated variance (posterior variance). We computed this approximation by fitting the within-unit model, given by (2.1)–(2.4) using diffuse normal priors on the itk, i.e. normal priors with mean 0 and variance 10 000. Specifically, we do not use the between-unit model, (2.5)–(2.6), when constructing this approximation to the network likelihood. We also note that Vi is a block diagonal because of the assumed independence across years, and our data in this setting are the i. Results in Daniels and Kass (1998) can be used to show that the approximation in (3.2) has the same order of error as a Laplace approximation, here 
where nit is the number of patients in network i in year t, for estimating the between-network parameters; in addition, results given in that paper show that the order of the error is unchanged by replacing the MLE with the posterior mean. Finally, the network-level data have been adjusted for patient mix and the correlation among the multivariate outcomes for each patient.
We use the DIC to aid in the choice of an appropriate form for the correlation matrices although it could also be used to choose the number of latent variables in the between-unit model (2.5)–(2.6).
We also use this approximation to conduct posterior predictive checks (Gelman et al., 2003) for assessing model fit. In this context, we determine whether the features of the posterior distribution of itk from fitting the full model are consistent with the data, itk, which are based on (3.2). Our goal is to assess the fit of the between-network model. Because the network-level ‘data’ constructed using the approximation given by (3.2) did not use the between-network model, we will not obtain overly optimistic results from the posterior predictive checks. This approach is somewhat nonstandard as the are not observed data—the yitj are observed. To do this, we sample a set of r using (3.2) given the current values of the parameters, 
sampled from i|y using the Gibbs sampler.
We examine two correlations: (1) the correlation of the components of i across k given t, corr(ikt, ik't) and (2) the correlation of the components of i across t given k, corr(ikt, ikt'). For each, we compute a single summary based on the data, itk, and a distribution of the quantity from the posterior predictive distribution using 
. We summarize these checks via posterior predictive p-values. If the model fits well, we expect the p-values to not be extreme (e.g. close to 0 or 1) as this would indicate that the particular feature of the data is not captured well by the model.
At each iteration of the Gibbs sampler, we also calculate a confidence region given by 
. A 95% cutoff for this region is calculated based on a 
2 distribution which should be approximately correct at each iteration of the sampler, given the parameters values at that iteration. We evaluate an empirical coverage probability by averaging across iterations. Because i is high dimensional (48), we examine the coverage for both this 48-dimensional vector and also for the 8-dimensional vector of i·t at each time, t = 1, ..., T.
|
4. APPLICATION TO THE VHA STUDY |
|---|
|
TOPSUMMARY1. INTRODUCTION2. JOINT MODELS FOR...3. INFERENCE4. APPLICATION TO THE...5. CONCLUDING REMARKSREFERENCES |
|---|
The continuous inpatient and outpatient responses were standardized to have mean 0 and variance 1. Patient-level covariates included in the within-unit model [see (2.1) and (2.4)] included age (categorized as 
39, 40–59, or 
60), race (categorized as white, black, or other), and primary mental health diagnosis (schizophrenia, other psychoses, posttraumatic stress disorder, alcohol abuse, or other). Over the 6-year period, the majority (95%) of patients were male; one-quarter black; 60% aged between 40 and 59 years of age at the time of their index discharge; and 27% had a primary mental health diagnosis of schizophrenia (data not shown). Table 1 indicates a 3.4% absolute decline between 1993 and 1998 in 6-month readmission rates, and a corresponding 5.3% increase in the outpatient visit rate over the 22 networks. These changes are consistent with expected patterns following health care reorganization. In what follows, networks are labeled using roman letters (A–W).
We ran five parallel chains, each of length 1100. The chains appeared to have converged within 100 iterations when examining the trace plots of selected parameters within each chain; we therefore used the 100 iterations as the burn-in. To minimize autocorrelation within the chains, we sampled every 10th iteration. The total posterior sample size was 500.
We fit a two-latent variable model (L = 2) where, a priori, the latent variables were thought to represent inpatient and outpatient quality. The prior specification given in Section 3.1 required specification of two hyperparameters, and c2. We wanted to specify priors flexible enough to verify the a priori belief of an inpatient and outpatient latent variable, but informative enough to avoid identifiability problems. We first set = 0.25 and c2 = 0.01 (c = 0.1) and then examined sensitivity to fixing /c and varying ; inferences were basically insensitive to this variation. We then fixed and varied the ratio. Values of /c less than two were getting close to identifiability problems with components of the latent variables switching and signs changing on the discrimination parameters. Values of the ratio over three seemed to be too informative. Thus, for inferences in Sections 4.4–4.6, we set = 0.25 and /c = 2.5.
The posterior means of the discrimination parameters were consistent with the prior which differentiated inpatient and outpatient measures (Table 2), although they were considerably smaller than their prior means. The discriminating abilities of the components of the inpatient and outpatient latent variables were not very consistent with equality among the four inpatient and outpatient variables. For example, the discrimination parameters for the number of visits and periods of continuity were larger than visits within 180 days and days to first visit. In addition, the posterior variances of the discrimination parameters were much tighter than the prior variances, indicating that the data very much determined the posterior on . Finally, there was no evidence that the discrimination parameters varied over time implying no need to index k by t.
|
To assess the importance of modeling the dependence of the networks over time, we first fit unstructured correlation matrices for both latent variables, Rl, l = 1, 2. The posterior mean of the correlation matrices for the inpatient and outpatient latent variables in the unstructured model are given in Table 3.
|
In general, the posterior means of the lag 1 correlations were around 0.5–0.7. The correlation generally decreased as the lag increased. The posterior distribution of the correlations were quite skewed with the posterior medians and mode for the large positive correlations closer to one than the posterior means (Figure 1). A first-order Markov model seemed most reasonable for R1. Because the lag k correlations did not vary much over time (moving down each off-diagonal of the correlation matrix), we fit a first-order Markov model with
l,t–1 = l. In this first-order model, the lag 1 correlations was quite high with a posterior mean (95% credible interval) of 0.71 (0.63, 0.79).
|
Ultimately, in addition to the unstructured and Markov specification of the correlation matrix, we also considered independence over time, Rl = I. The R = I model was fit to assess the potential efficiency gains over modeling the data separately by year as this model did not seem feasible when examining the estimated correlation matrices in Table 3.
Using the DIC, models that employed serial dependence for the inpatient latent variable fit best (Table 4). Fewer effective number of parameters in the dependence models which have more parameters in the Rl matrices is a result of more shrinkage of the i in these models. We used the best-fitting model, R1 = Rm, R2 = R, for the basis of inference. An alternative to picking the ‘best’ model would involve model averaging over the various specifications of the correlation matrix. We discuss this in Section 5.
|
Figure 2 displays scatterplots of the posterior means of the inpatient and outpatient latent variables for each network. These variables summarize the inpatient and outpatient quality of the networks over time with higher values corresponding to better care. Network W had consistently high inpatient quality, network N had consistently high outpatient quality, and network U had consistently poor inpatient and outpatient quality (lower left quadrant) across all years.
|
Figure 3 displays the posterior means of selected inpatient and outpatient latent variables for a subset of the 22 networks. Some networks appear to be ‘significantly’ improving or worsening over the 6-year period. For example, network V appears to provide improving outpatient quality while network P appears to have worsening outpatient quality.
|
As a comparison, the probabilities in parentheses below are derived from assuming independence over time, the R = I model. There was no strong evidence that any of the networks provided consistently poor care over time using using (2.12). Setting q = 0.20 (20th percentile), the highest probability of being in the lower 20% for all 6 years was network U on inpatient quality with probability 0.49 (0.20 under R = I) and networks U and L on outpatient quality, with probabilities 0.37 (0.04 under R = I) and 0.18 (0.04 under R = I), respectively. In terms of excellent inpatient quality (2.13), network W was in the upper 20% for all 6 years for the inpatient latent variable with probability 0.89 (0.54 under R = I) while network D had probability 0.51 (0.15 under R = I). For excellent outpatient quality, network N had probability 0.52 (0.14 under R = I). From these results, it is clear that inferences change considerably under R = I model and stronger inferences are possible about the behavior of the networks over time by accounting for the temporal correlation.
As a less stringent criterion for longitudinal performance, we computed posterior probabilities of a positive (negative) slope for each latent variable for each network as given by (2.14). Networks C, J, and V had posterior probabilities of 0.85 or greater for improving inpatient quality over the 6-year period while network M had a similar probability associated with worsening inpatient quality. Using the same criterion, networks F, K, L, and V had large probabilities for outpatient quality improvement. Degrading outpatient quality trends were associated with networks C, O, and P.
We first examined the correlation among the i corresponding to each of the eight responses, for each time, t, corr(ikt, ik't), in order to ensure the model characterized the correlation among the network effects. This check resulted in 168 p-values with only 3 outside the interval (0.05, 0.95) and none outside the interval (0.01,0.99). We then examined the correlation among the i over time for each of the eight responses, corr(ikt, ikt'), to assess the appropriateness of temporal correlations. This check resulted in 120 p-values. The fit was not as good as the previous check, but not unreasonable—24 of the 120 p-values were outside the interval (0.05, 0.95), but only 5 were outside (0.01, 0.99). There was no consistent pattern in terms of particular correlations not being captured well. In addition, for the extreme p-values, the observed correlation was typically quite close to the observed extreme of the posterior predictive distribution. Finally, the average coverage for the confidence regions described in Section 3.2 for the entire i vector was 0.86 and for the lower dimensional i·k vectors, 0.92, which is quite good given the complexity of the data and the ‘parsimonious’ hierarchical model.
|
5. CONCLUDING REMARKS |
|---|
|
TOPSUMMARY1. INTRODUCTION2. JOINT MODELS FOR...3. INFERENCE4. APPLICATION TO THE...5. CONCLUDING REMARKSREFERENCES |
|---|
We have proposed a multilevel multidimensional model to characterize units measured repeatedly on the basis of within-unit continuous and discrete outcomes. The aggregation level of interest, the health care network, was fit fairly well by the model as assessed by posterior predictive checks. We also observed an increase in efficiency in terms of longitudinal profiling by exploiting the temporal correlation in our model for the data. In our problem, subjects were not measured repeatedly over time, rather the health care unit was measured repeatedly. We introduced temporal correlation at the network level through modeling the pairwise correlations of the latent variables. An alternative approach could involve reparameterization of the matrix from pairwise correlations to partial correlations (Wong et al., 2003
). We note however that the correlation structure of longitudinal data is difficult to model parsimoniously using partial correlations and additional constraints would be required. An important area where more research is needed relates to model fit. Given the complexity of the model, we used an approximate network-level likelihood to evaluate the DIC. Moreover, we focused on the fit of the model at the network level. In complex models, such as those we proposed here, there are several steps where choices are made: the number of latent variables underlying the observed data, the structure of the serial correlation, and the inclusion of covariates. While we did not include unit (network)-level covariates, these could be easily accommodated in the dependence model for the latent variables.
In this paper, we fixed the number of latent variables based on prior research, theory, and expert opinion. An extension would involve accounting for uncertainty in the number of latent variables through a model-averaging approach. However, given the goals of this paper, such an approach would greatly hinder interpretation. As pointed out by a referee, an interesting extension would involve model averaging over the different covariance structures for the latent variables. To accomplish this, methods to compute the marginal distribution of the data would need to be developed and ideally a broader class of covariance models considered.
|
ACKNOWLEDGMENTS |
|---|
This research was supported by Grants CA85295 (Daniels) from the National Cancer Institute and MH61434 (Normand) from the National Institute of Mental Health. The VHA data were generously provided through the efforts of Robert Rosenheck, MD, of West Haven, CT. The authors are grateful to Myung Joon Kim for programming expertise.
|
REFERENCES |
|---|
|
TOPSUMMARY1. INTRODUCTION2. JOINT MODELS FOR...3. INFERENCE4. APPLICATION TO THE...5. CONCLUDING REMARKSREFERENCES |
|---|
-
AGUILAR, O. AND WEST, M. (1998). Analysis of hospital quality monitors using hierarchical time series models. In Gatsonis, C., Kass, R. E. et al. (eds), Case Studies in Bayesian Statistics, Volume 4. New York, NY: Springer, pp. 287–302.
BARNARD, J., MCCULLOCH, R. AND MENG, X. L. (2000). Modeling covariance matrices in terms of standard deviations and correlations, with application to shrinkage. Statistica Sinica 10, 1281–1311.
BERLOWITZ, D. R., YOUNG, G. J., BRANDEIS, G. H., BORIS, K. AND ANDERSON, J. J. (2001). Health care reorganization and quality of care: unintended effects on pressure ulcer prevention. Medical Care 39, 138–146.[CrossRef][Medline]
BRONSKILL, S., NORMAND, S.-L., LANDRUM, M. B. AND ROSENHECK, R. (2002). Longitudinal profiles of health care providers. Statistics in Medicine 21, 1067–1088.[CrossRef][Web of Science][Medline]
CHEN, S., SMITH, M. W., WAGNER, T. D. AND BARNETT, P. G. (2003). Spending for specialized mental health treatment in the VA: 1995–2000. Health Affairs 22, 256–263.
DANIELS, M. J. (2005). Bayesian modelling of several covariance matrices and some results on the propriety of the posterior for linear regression with correlated and/or heterogeneous errors. Journal of Multivariate Analysis (in press).
DANIELS, M. J. AND KASS, R. E. (1998). A note on first stage approximation in two stage hierarchical models. Sankhya, Series B. 60, 19–30.
DOUGLAS, J. A. (1999). Item response models for longitudinal quality of life data in clinical trials. Statistics in Medicine 18, 2917–2931.[Medline]
DOUGLAS, J. A., KOSOROK, M. R. AND CHEWNING, B. A. (1999). A latent variable model or discrete multivariate psychometric waiting times. Psychometrika 64, 69–82.
DUNCAN, T. E., DUNCAN, S. C., OKUT, H., STRYCKER, L. A. AND LI, F. (2002). An extension of the general latent variable growth modeling framework to four levels of the hierarchy. Structural Equation Modeling 9 303–326.[CrossRef]
DUNSON, D. B. (2000). Bayesian latent variable models for clustered mixed outcomes. Journal of the Royal Statistical Society, B 62, 355–366.[CrossRef]
DUNSON, D. B. (2003). Dynamic latent trail models for multidimensional longitudinal data. Journal of the American Statistical Association 98, 555–563.[CrossRef]
DUNSON, D. B., CHEN, Z. AND HARRY, J. (2003). A Bayesian approach for joint modeling of cluster size and subunit-specific outcomes. Biometrics 59, 521–530.[Medline]
GELMAN, A., CARLIN, J. B., STERN, H. S. AND RUBIN, D. B. (2003). Bayesian Data Analysis, 2nd edition. Chapman and Hall.
GOLDSTEIN, H. AND SPIEGELHALTER, D. J. (1996). Statistical aspects of institutional performance: league tables and their limitations (with discussion). Journal of the Royal Statistical Society, A 159, 385–444.[CrossRef]
GUEORGUIEVA, R. V. AND AGRESTI, A. (2001). A correlated probit model for joint modeling of clustered binary and continuous responses. Journal of the American Statistical Association 96, 1102–1112.[CrossRef]
HOFF, R. A., ROSENHECK, R. A., MATERKO, M. AND WILSON, N. J. (1998). The quality of VA mental health services. Administration and Policy in Mental Health 26, 45–56.[CrossRef][Web of Science][Medline]
KIZER, K., DEMAKIS, J. AND FEUSSNER, J. (2000). Reinventing VA health care: systematizing quality improvement and quality innovation. Medical Care 38 (Suppl 1), I7–I16.[CrossRef][Medline]
LANDRUM, M. B., NORMAND, S. L. T. AND ROSENHECK, R. A. (2003). Selection of related multivariate means: monitoring psychiatric care in the Department of Veterans Affairs. Journal of the American Statistical Association 98, 7–16.[CrossRef]
LEE, S. Y., POON, W. Y. AND BENTLER, P. M. (1992). Structural equation models with continuous and polytomous variables. Psychometrika 57, 89–105.[CrossRef]
LEE, S. Y. AND SHI, J. Q. (2001). Maximum likelihood estimation of two-level latent variable models with mixed continuous and polytomous data. Biometrics 57, 767–794.
LIN, H., TURNBULL, B., MCCULLOCH, C. AND SLATE, E. (2002). Latent class models for joint analysis of longitudinal biomarker and event process data: application to longitudinal prostate-specific antigen readings and prostate cancer. Journal of the American Statistical Association 97, 53–65.[CrossRef]
LOPES, H. F. AND WEST, M. (2004). Bayesian model assessment in factor analysis. Statistica Sinica 14, 41–67.
MARK, T. AND COFFEY, R. (2000). Spending on mental health and substance abuse treatment, 1987–1997. Health Affairs July/August, 108–120.
MCCLELLAN, M. AND STAIGER, D. (1999). The quality of health care providers. NBER Working Paper 7327. Cambridge, MA: National Bureau of Economic Research.
MUELLER, R. O. (1996). Basic Principles of Structural Equation Modeling: An Introduction to LISREL and EQS. New York: Springer.
MUTHÉN, B., BROWN, C. H., MASYN, K., JO, B., KHOO, S.-T., YANG, C.-C., WANG, C.-P., KELLAM, S. G., CARLIN, J. B. AND LIAO, J. (2002). General growth mixture modeling for randomized preventive interventions. Biostatistics 3, 459–475.[Abstract]
NORMAND, S.-L., GLICKMAN, M. E. AND GATSONIS, C. A. (1997). Statistical methods for profiling providers of medical care: issues and applications. Journal of the American Statistical Association 92, 803–814.[CrossRef]
OORT, F. J. (2001). Three-mode models for multivariate longitudinal data. The British Journal of Mathematical and Statistical Psychology 54, 49–78.[Medline]
ROSENHECK, R. AND DELILLA, D. (1999). Department of Veterans Affairs National Mental Health Program Performance Monitoring System: Fiscal Year 1998 Report. Northeast Program Evaluation Center, VA Health Services Research & Development Service. West Haven, CT: VA Connecticut Healthcare System.
ROY, J. AND LIN, X. (2000). Latent variable models for longitudinal data with multiple continuous outcomes. Biometrics 56, 1047–1054.[CrossRef][Web of Science][Medline]
SPIEGELHALTER, D. J., BEST, N. G., CARLIN, B. P. AND VAN DER LINDE, A. (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society, Series B 64, 583–616.[CrossRef]
UNITED STATES SENATE (1999). Hearing of ‘quality of care in the Veterans Affairs health care system’, September 22, 1998. One Hundred Fifth Congress, Second Session. Washington, DC: Government Printing Office, p. 62.
WONG, F., CARTER, C. K. AND KOHN, R. (2003). Efficient estimation of covariance selection models. Biometrika 90, 809–830.
ZIMMERMAN, D. AND NUNEZ-ANTON, V. (1997). Structured antedependence models for longitudinal data. In Gregoire, T. G., Brillinger, D. R., Diggle, P. J., Russek-Cohen, E., Warren, W. G. and Wolfinger, R. D. (eds), Modelling Longitudinal Data and Spatially Correlated Data: Methods, Applications and Future Directions. Springer, pp. 63–76.
Received August 16, 2004; revised May 6, 2005; accepted for publication May 20, 2005.
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  M. Plummer Penalized loss functions for Bayesian model comparison Biostat., July 1, 2008; 9(3): 523 - 539. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 D. B. Dunson Bayesian methods for latent trait modelling of longitudinal data Statistical Methods in Medical Research, October 1, 2007; 16(5): 399 - 415. [Abstract] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||