Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Biostatistics. Author manuscript; available in PMC 2009 December 10.
Published in final edited form as:
PMCID: PMC2791405

Longitudinal profiling of health care units based on continuous and discrete patient outcomes


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

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.

Table 1
Characteristics of service network outcomes in the VHA, 1993–1998

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

2.1 Profiling VHA mental health treatment quality

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.

2.2 Within-network model

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.

Within-network model for binary outcomes

Assume for binary-valued variables Zitjk ~ N(ξitjk, 1) such that Yitjk = I(Zitjk > 0), where


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 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.

Within-network model for continuous outcomes

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 [set membership] {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


where Σt(U) denotes the covariance of the continuous variables at time t. See Table 1 for an ordered list of the eight patient-level measures.

2.3 Between-network model

The K network effects are assumed to be represented by L latent variables that are a priori independent,


where λk = (λ1k,..., λLk)′ is an L × 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, ψtk2, represents heterogeneity among outcomes in a given year not explained by the L latent variables. Rl denotes a T × T correlation matrix that permits serial dependence in the L latent variables.

2.4 Serial dependence model for 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 ρl,t12+σlt2=1 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


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 σlt2=σl2 . 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.

2.5 Model identifiability

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.12.4), with


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


where A is a KT × KT matrix with K (T × T) blocks on the main diagonal, the kth of which has the form l=1Lλlk2Rl+diag(ψk12,,ψkT2) and with off-diagonal blocks, Ojk, of the form, lλljRlλlk.

For each t,


where ηi·t = (ηi1t,..., ηiKt)′, the kth row of K × 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,


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.

2.6 Characterizing networks 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


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


The trend may be characterized more formally by estimating and ranking the slope, Bθil. For example, we can compute


where Bθil = (X′X)−1 X′θil, and X is a 2 × 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.


for ll′. An overall summary may be represented as θit=lwlθilt for some selected weights, wl.

3. Inference

3.1 Prior distributions

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 (π, c2 ILK), where


We discuss the choices of π[open star] > 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, 1τt2 and 1ψtk2, with parameters (0.001, 0.001). A Wishart prior is placed on Σt(U)1 with degrees of freedom equal to dim(Σt(U)) 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.

3.2 Sampling and model fit

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 × 1 vectors βtk, for k = 1 and k = 5; the (3P) × 1 vectors of (βt2, βt3, βt4) and (βt6, βt7, βt8); the matrices Σt(U) for U and t; and in the between-network model, the L × 1 vector of λk for k = 1,..., K, and the T × 1 vectors θil for l = 1,..., L. The other parameters were sampled component-wise.

Model comparison and fit

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 × 1 dimensional response vector and ηit is a K × 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


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 O(it1nit), 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, ηi(r), 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 η^itkr. 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 (η^iηi(r))Vi1(η^iηi(r)). 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

4.1 Details on the VHA National Mental Health Performance Monitoring System data

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).

4.2 Sampling algorithm

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.

4.3 Prior for discrimination matrix, Λ

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, π[open star] 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 π[open star] = 0.25 and c2 = 0.01 (c = 0.1) and then examined sensitivity to fixing π[open star]/c and varying π[open star]; inferences were basically insensitive to this variation. We then fixed π[open star] and varied the ratio. Values of π[open star]/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 π[open star] = 0.25 and π[open star]/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.

Table 2
Posterior means and 95% credible intervals for the discrimination parameters, λk for the inpatient and outpatient latent variables for the final model given in Section 4.3

4.4 Serial dependence

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.

Table 3
Posterior mean of unstructured correlation matrices for inpatient (R1) and outpatient (R2) latent variables

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).

Fig. 1
Posterior distributions of several of the components of the Rl matrices under the unstructured R models.

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.

Table 4
DIC. Markov(1) corresponds to the first-order Markov model discussed in Section 2.4

4.5 Longitudinal patterns of quality

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.

Fig. 2
Posterior means for inpatient and outpatient latent variables for each year.

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.

Fig. 3
Posterior means and 95% pointwise credible intervals for selected inpatient and outpatient latent variables over time.

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.

4.6 Posterior predictive checks

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

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.


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.

Contributor Information

MICHAEL J. DANIELS, Department of Statistics, University of Florida, Gainesville, FL 32611, USA ; ude.lfu.tats@sleinadm.

SHARON-LISE T. NORMAND, 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.


  • Aguilar O, West M. Analysis of hospital quality monitors using hierarchical time series models. In: Gatsonis C, Kass RE, et al., editors. Case Studies in Bayesian Statistics. Vol. 4. Springer; New York, NY: 1998. pp. 287–302.
  • Barnard J, McCulloch R, Meng XL. Modeling covariance matrices in terms of standard deviations and correlations, with application to shrinkage. Statistica Sinica. 2000;10:1281–1311.
  • Berlowitz DR, Young GJ, Brandeis GH, Boris K, Anderson JJ. Health care reorganization and quality of care: unintended effects on pressure ulcer prevention. Medical Care. 2001;39:138–146. [PubMed]
  • Bronskill S, Normand S-L, Landrum MB, Rosenheck R. Longitudinal profiles of health care providers. Statistics in Medicine. 2002;21:1067–1088. [PubMed]
  • Chen S, Smith MW, Wagner TD, Barnett PG. Spending for specialized mental health treatment in the VA: 1995–2000. Health Affairs. 2003;22:256–263. [PubMed]
  • Daniels MJ. 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. 2005 (in press)
  • Daniels MJ, Kass RE. A note on first stage approximation in two stage hierarchical models. Sankhya, Series B. 1998;60:19–30.
  • Douglas JA. Item response models for longitudinal quality of life data in clinical trials. Statistics in Medicine. 1999;18:2917–2931. [PubMed]
  • Douglas JA, Kosorok MR, Chewning BA. A latent variable model or discrete multivariate psychometric waiting times. Psychometrika. 1999;64:69–82.
  • Duncan TE, Duncan SC, Okut H, Strycker LA, Li F. An extension of the general latent variable growth modeling framework to four levels of the hierarchy. Structural Equation Modeling. 2002;9:303–326.
  • Dunson DB. Bayesian latent variable models for clustered mixed outcomes. Journal of the Royal Statistical Society, B. 2000;62:355–366.
  • Dunson DB. Dynamic latent trail models for multidimensional longitudinal data. Journal of the American Statistical Association. 2003;98:555–563.
  • Dunson DB, Chen Z, Harry J. A Bayesian approach for joint modeling of cluster size and subunit-specific outcomes. Biometrics. 2003;59:521–530. [PubMed]
  • Gelman A, Carlin JB, Stern HS, Rubin DB. Bayesian Data Analysis. 2nd edition Chapman and Hall; 2003.
  • Goldstein H, Spiegelhalter DJ. Statistical aspects of institutional performance: league tables and their limitations (with discussion). Journal of the Royal Statistical Society, A. 1996;159:385–444.
  • Gueorguieva RV, Agresti A. A correlated probit model for joint modeling of clustered binary and continuous responses. Journal of the American Statistical Association. 2001;96:1102–1112.
  • Hoff RA, Rosenheck RA, Materko M, Wilson NJ. The quality of VA mental health services. Administration and Policy in Mental Health. 1998;26:45–56. [PubMed]
  • Kizer K, Demakis J, Feussner J. Reinventing VA health care: systematizing quality improvement and quality innovation. Medical Care. 2000;38(Suppl 1):I7–I16. [PubMed]
  • Landrum MB, Normand SLT, Rosenheck RA. Selection of related multivariate means: monitoring psychiatric care in the Department of Veterans Affairs. Journal of the American Statistical Association. 2003;98:7–16.
  • Lee SY, Poon WY, Bentler PM. Structural equation models with continuous and polytomous variables. Psychometrika. 1992;57:89–105.
  • Lee SY, Shi JQ. Maximum likelihood estimation of two-level latent variable models with mixed continuous and polytomous data. Biometrics. 2001;57:767–794. [PubMed]
  • Lin H, Turnbull B, McCulloch C, Slate E. 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. 2002;97:53–65.
  • Lopes HF, West M. Bayesian model assessment in factor analysis. Statistica Sinica. 2004;14:41–67.
  • Mark T, Coffey R. Spending on mental health and substance abuse treatment, 1987–1997. Health Affairs. 2000:108–120. July/August. [PubMed]
  • McClellan M, Staiger D. NBER Working Paper 7327. National Bureau of Economic Research; Cambridge, MA: 1999. The quality of health care providers.
  • Mueller RO. Basic Principles of Structural Equation Modeling: An Introduction to LISREL and EQS. Springer; New York: 1996.
  • Muthén B, Brown CH, Masyn K, Jo B, Khoo S-T, Yang C-C, Wang C-P, Kellam SG, Carlin JB, Liao J. General growth mixture modeling for randomized preventive interventions. Biostatistics. 2002;3:459–475. [PubMed]
  • Normand S-L, Glickman ME, Gatsonis CA. Statistical methods for profiling providers of medical care: issues and applications. Journal of the American Statistical Association. 1997;92:803–814.
  • Oort FJ. Three-mode models for multivariate longitudinal data. The British Journal of Mathematical and Statistical Psychology. 2001;54:49–78. [PubMed]
  • Rosenheck R, Delilla D. Department of Veterans Affairs National Mental Health Program Performance Monitoring System: Fiscal Year 1998 Report. VA Connecticut Healthcare System; West Haven, CT: 1999. Northeast Program Evaluation Center, VA Health Services Research & Development Service.
  • Roy J, Lin X. Latent variable models for longitudinal data with multiple continuous outcomes. Biometrics. 2000;56:1047–1054. [PubMed]
  • Spiegelhalter DJ, Best NG, Carlin BP, van der Linde A. Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society, Series B. 2002;64:583–616.
  • United States Senate 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. 1999. p. 62.
  • Wong F, Carter CK, Kohn R. Efficient estimation of covariance selection models. Biometrika. 2003;90:809–830.
  • Zimmerman D, Nunez-Anton V. Structured antedependence models for longitudinal data. In: Gregoire TG, Brillinger DR, Diggle PJ, Russek-Cohen E, Warren WG, Wolfinger RD, editors. Modelling Longitudinal Data and Spatially Correlated Data: Methods, Applications and Future Directions. Springer; 1997. pp. 63–76.