Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2791405

Formats

Article sections

- Summary
- 1. Introduction
- 2. Joint models for longitudinal multilevel continuous and discrete data
- 3. Inference
- 4. Application to the VHA study
- 5. Concluding remarks
- References

Authors

Related links

Biostatistics. Author manuscript; available in PMC 2009 December 10.

Published in final edited form as:

Published online 2005 May 25. doi: 10.1093/biostatistics/kxi036

PMCID: PMC2791405

NIHMSID: NIHMS143243

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

The publisher's final edited version of this article is available at Biostatistics

See other articles in PMC that cite the published article.

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.

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.

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 *Y _{itjk}* denote the value of the

Assume for binary-valued variables *Z _{itjk}* ~

$${\xi}_{itjk}={\beta}_{0tk}+{\beta}_{1tk}({x}_{1itj}-{\stackrel{\u2012}{x}}_{1t})+\cdots +{\beta}_{Ptk}({x}_{Pitj}-{\stackrel{\u2012}{x}}_{Pt})+{\varphi}_{itj}+{\eta}_{itk},$$

(2.1)

$${\varphi}_{itj}~N(0,{\tau}_{t}^{2}).$$

(2.2)

Here, *ϕ _{itj}* is a subject-specific random effect,

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

$${Y}_{itj}^{\left(U\right)}\mid {x}_{itj}^{\left(U\right)}=1~{N}_{3}({\xi}_{itj}^{\left(U\right)},{\Sigma}_{t}^{\left(U\right)}),$$

(2.3)

$${\xi}_{itjk}={\beta}_{0tk}+{\beta}_{1tk}({x}_{1itj}-{\stackrel{\u2012}{x}}_{1t})+\cdots +{\beta}_{Ptk}({x}_{Pitj}-{\stackrel{\u2012}{x}}_{Pt})+{\eta}_{itk},$$

(2.4)

where ${\Sigma}_{t}^{\left(U\right)}$ 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,

$${\eta}_{itk}\mid {\theta}_{il},{\lambda}_{k},{\psi}_{tk}^{2}~N({\lambda}_{1k}{\theta}_{i1t}+\cdots +{\lambda}_{Lk}{\theta}_{iLt},{\psi}_{tk}^{2}),$$

(2.5)

$${\theta}_{il}={({\theta}_{il1},\dots ,{\theta}_{ilT})}^{\prime}~{N}_{T}(0,{\mathbf{R}}_{l}),\phantom{\rule{1em}{0ex}}l=1,2,\dots ,L,$$

(2.6)

where *λ _{k}* = (

We model the within-unit dependency of the outcomes measured at different times through the correlation matrix, **R*** _{l}*. Our motivation for specifying some structure on

We consider first-order Markov models. Recall that *θ _{il}* ~

$$p({\theta}_{ilt}\mid {\theta}_{ilj},j<t)=p({\theta}_{ilt}\mid {\theta}_{ilt-1})=N({\rho}_{l,t-1}{\theta}_{ilt-1},{\sigma}_{lt}^{2}).$$

Using iterated expectations, it is easy to show that Var(*θ _{ilt}*) = 1 implies that ${\rho}_{l,t-1}^{2}+{\sigma}_{lt}^{2}=1$ for

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

$$z\left({\rho}_{l,t-1}\right)={g}_{l}(t-1;\alpha ),$$

(2.7)

where *z*(·) is Fisher's *z*-transformation and *g _{l}*(

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

$${\widehat{\eta}}_{i}~{N}_{KT}({\eta}_{i},{V}_{i}),\phantom{\rule{1em}{0ex}}i=1,\dots ,I,$$

(2.8)

where *V _{i}* is assumed known and

$${\eta}_{i}~{N}_{KT}(0,A),\phantom{\rule{1em}{0ex}}i=1,\dots ,I,$$

(2.9)

where *A* is a *KT* × *KT* matrix with *K* (*T* × *T*) blocks on the main diagonal, the *k*th of which has the form ${\sum}_{l=1}^{L}{\lambda}_{lk}^{2}{\mathbf{R}}_{l}+\text{diag}({\psi}_{k1}^{2},\dots ,{\psi}_{kT}^{2})$ and with off-diagonal blocks, *O*_{jk}, of the form, ${\sum}_{l}{\lambda}_{lj}{\mathbf{R}}_{l}{\lambda}_{lk}$.

For each *t*,

$${\eta}_{i\cdot t}~{N}_{K}(0,\Lambda {\Lambda}^{\prime}+{\Psi}_{\cdot t}),$$

(2.10)

where *η*_{i·t} = (*η*_{i1t},..., *η*_{iKt})′, the *k*th row of *K* × *L* matrix Λ is (*λ*_{k1},..., *λ*_{kL})′, and Ψ_{·t} = diag(*ψ*_{1t},..., *ψ*_{Kt}). This resembles a standard factor analytic model. Because *V _{i}* is known, Λ and Ψ

$${\eta}_{ik\cdot}~{N}_{T}\left(0,\sum _{l=1}^{L}{\lambda}_{lk}^{2}{\mathbf{R}}_{l}+{\Psi}_{k\cdot}\right),$$

(2.11)

where *η*_{ik·} = (*η*_{ik1},..., *η*_{ikT})′. As the components of *λ*_{kl} and Ψ_{k·} = diag(*ψ*_{k1},..., *ψ*_{kT}) are identified from (2.10), **R*** _{l}* will be identified from the correlation among

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

$${P}_{{\theta}_{il}}^{q}=P({\theta}_{il1}<{T}_{q}\left({\theta}_{l1}\right),\dots ,{\theta}_{ilT}<{T}_{q}\left({\theta}_{lT}\right)\mid \mathbf{Y}),$$

(2.12)

for each *i* and *l*. Here, *T _{q}*(

$${P}_{{\theta}_{il}}^{q}=P({\theta}_{il1}>{T}_{q}\left({\theta}_{l1}\right),\dots ,{\theta}_{ilT}>{T}_{q}\left({\theta}_{lT}\right)\mid \mathbf{Y}).$$

(2.13)

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

$$P({B}_{{\theta}_{il}}<0\mid \mathbf{Y}),$$

(2.14)

where *B _{θil}* = (

Changes in the relationship between the latent variables may be assessed through the correlation, e.g.

$$\text{Cor}({\theta}_{ilt},{\theta}_{i{l}^{\prime}t});\phantom{\rule{1em}{0ex}}t=1,2,\dots ,T,$$

(2.15)

for *l* ≠ *l′*. An overall summary may be represented as ${\stackrel{\u2012}{\theta}}_{it}={\sum}_{l}{w}_{l}{\theta}_{ilt}$ for some selected weights, *w _{l}*.

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 ** λ** ~

$$\pi ={\left(\begin{array}{c}\hfill {\pi}_{1k}\hfill \\ \hfill {\pi}_{2k}\hfill \end{array}\right)}^{\prime}={\pi}^{\star}{\left(\begin{array}{cccccccc}\hfill -1& \hfill -1& \hfill -1& \hfill 1\hfill & \hfill 0& \hfill 0\hfill & \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0\hfill & \hfill -1& \hfill 1\hfill & \hfill -1& \hfill -1\end{array}\right)}^{\prime}.$$

(3.1)

We discuss the choices of *π* > 0 and *c*^{2} 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\u2215{\tau}_{t}^{2}$ and $1\u2215{\psi}_{tk}^{2}$, with parameters (0.001, 0.001). A Wishart prior is placed on ${\Sigma}_{t}^{{\left(U\right)}_{-1}}$ with degrees of freedom equal to $\text{dim}\left({\Sigma}_{t}^{\left(U\right)}\right)$ 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,

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, **R*** _{l}*,

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

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

$$\prod _{t}\int \int \int \prod _{i}\prod _{j}\int p({y}_{itj}\mid {\eta}_{it},{\beta}_{t},{\varphi}_{itj},{\Sigma}_{t})p({\varphi}_{itj}\mid {\tau}_{t}^{2})p\left({\Sigma}_{t}\right)p\left({\beta}_{t}\right)p\left({\tau}_{t}^{2}\right)\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{\varphi}_{itj}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{\beta}_{t}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{\tau}_{t}^{2}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{\Sigma}_{t},$$

where *y _{itj}* is the

$${\widehat{\eta}}_{i}~{N}_{KT}({\eta}_{i},{V}_{i}),\phantom{\rule{1em}{0ex}}i=1,\dots ,I,$$

(3.2)

where ${\widehat{\eta}}_{i}$ is the maximum likelihood estimate (MLE) of *η _{i}* = {

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, ${\widehat{\eta}}_{itk}$, which are based on (3.2). Our goal is to assess the fit of the

We examine two correlations: (1) the correlation of the components of *η _{i}* across

At each iteration of the Gibbs sampler, we also calculate a confidence region given by ${({\widehat{\eta}}_{i}-{\eta}_{i}^{\left(r\right)})}^{\prime}{V}_{i}^{-1}({\widehat{\eta}}_{i}-{\eta}_{i}^{\left(r\right)})$. 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 ${\widehat{\eta}}_{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*.

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 *c*^{2}. 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 *c*^{2} = 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

To assess the importance of modeling the dependence of the networks over time, we first fit unstructured correlation matrices for both latent variables, **R*** _{l}*,

Posterior mean of unstructured correlation matrices for inpatient (**R**_{1}) *and outpatient* (**R**_{2}) *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 **R**_{1}. 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).

Posterior distributions of several of the components of the **R**_{l} 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, **R*** _{l}* =

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 **R*** _{l}* matrices is a result of more shrinkage of the

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.

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.

We first examined the correlation among the *η _{i}* corresponding to each of the eight responses, for each time,

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.

MICHAEL J. DANIELS, Department of Statistics, University of Florida, Gainesville, FL 32611, USA ; Email: 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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |