Daniels M.J., Hogan J.W. Missing Data in Longitudinal Studies: Strategies for Bayesian Modeling and Sensitivity Analysis

Подождите немного. Документ загружается.

140 INFERENCE UNDER MAR

obs

)toobtain

DIC

=DIC

full

obs

)

= E

mis

{DIC

full

obs

, y

mis

)}

= −4E

mis

{(θ | y

obs

, y

mis

)}

+2E

mis

(

θ(y

mis

) | y

obs

, y

mis

)

, (6.16)

where E

mis

(·)istheexpectation with respect to p(y

mis

| y

obs

). Note that

we took a similar expectation with respect to p(y

mis

| y

obs

)toderivethe

observed data posterior from the full data posterior in (6.3).

DIC

cannot be computed directly in WinBUGS. To compute (6.16), we

use the samples from the data augmented Gibbs sampling algorithm; de-

note this sample as {y

(m)

mis

, θ

(m)

: m =1,...,M}.Theﬁrst term in (6.16)

is computed by averaging (θ | y

obs

, y

mis

)overthesample{y

(m)

mis

, θ

(m)

: m =

1,...,M}.Forthesecond term, we need to compute

θ(y

(m)

mis

)=E(θ | y

obs

, y

(m)

mis

)

at each iteration m,andthenaverage the entire term, {

θ(y

mis

) | y

obs

, y

mis

over the sample {y

(m)

mis

: m =1,...,M}.

Acomputational problem with this approach is that

θ(y

(m)

mis

)will typically

not be available in closed form. A straightforward, but computationally im-

practical approach would be to run a new Gibbs sampler for every value of

mis

sampled and compute θ(y

mis

)fromthese draws.

We recommend two alternative approaches that are more practical compu-

tationally. First, recall we have a sample {y

(m)

mis

, θ

(m)

: m =1,...,M} from

p(y

mis

, θ | y

obs

). The ﬁrst approach is to reweight this sample to estimate

E(θ | y

obs

, y

(m)

mis

)forall values y

(m)

mis

in the sample. In particular, we can use

the weightedaverage

∗

θ(y

(m)

mis

)=E(θ | y

obs

, y

(m)

mis

) ≈



a=1

(m)

(a)



a=1

(m)

, (6.17)

with weights w

(m)

given by

(m)

p(θ

(a)

, y

(m)

mis

, y

obs

)/p(y

(m)

mis

| y

obs

)

p(θ

(a)

, y

(a)

mis

, y

obs

)

and p(θ

(l)

, y

(l)

mis

, y

obs

)=p(y

obs

, y

(l)

mis

| θ)p(θ). The term p(y

(m)

mis

| y

obs

)canbe

∗

The approximate equalities in this section, signiﬁed by ≈,aredue only to Monte Carlo

error.

MODEL SELECTION AND MODEL FIT 141

estimated by

p(y

(m)

mis

| y

obs

) ≈



j=1

p(y

(m)

mis

| y

obs

, θ

(j)

where {θ

(j)

: j =1,...,M} is the sample from p(θ | y

obs

). This approach does

not require re-running the sampling algorithm and we expect these weights

to be fairly stable.

Asecond computational approach would be to run one additional Gibbs

sampler for a ﬁxed value of y

mis

,sayy

mis

= E(y

mis

| y

obs

); denote this sample

as {θ

(l)

: l =1,...,L}.Wecanthenuse the simpler weights

(m)

p(θ

(l)

, y

(m)

mis

, y

obs

)

p(θ

(l)

, y

mis

, y

obs

)

in (6.17), which are available in closed form. The sum is taken over the second

sample, {θ

(l)

: l =1,...,L}

Parts of DIC

can be computed in WinBUGS while other parts need to be

computed in other software (like R).

Summary

Both DIC constructions need further exploration and comparison to deter-

mine their behavior in the setting of ignorable missingness. DIC

removes

the missing data by averaging over the predictive distribution conditional

on θ,

p(y

obs

| θ)=



p(y

obs

, y

mis

| θ)dy

mis

to obtain the observed data likelihood with which the DIC is then constructed.

Using the observed-data log likelihood in DIC

is very similar to the AIC, a

frequentist model selection criterion.

On the other hand, DIC

removes the missing data by averaging the DIC

based on the full-data model, DIC

full

obs

, y

mis

)withrespect to the posterior

predictive distribution of y

mis

, p(y

mis

| y

obs

). So, essentially, this form of

the DIC is a weighted average of the DIC based on the complete data with

weights equal to how likely these ‘completed’ datasets are under the model

and the observed data. That the full-data response model and its parameters

are typically of primary interest provides some support for this approach. For

the data examples in Chapter 7, we use DIC

6.7.2 Posterior predictive checks

In Section 3.5.3, we discussed posterior predictive checks based on complete

data. Assessing model ﬁt for incomplete data using posterior predictive checks

142 INFERENCE UNDER MAR

can be based on statistics computed from replications of the observed data

or replications of completed datasets. A recent paper by Gelman et al. (2005)

advocates doing checks using completed (full-) data. The idea is to use the

value of y

mis

at each iteration of the data-augmented Gibbs sampler to create

an ‘observed complete’ dataset to then compare to a ‘replicated complete’

dataset. We provide more details below. In using this approach, we must keep

in mind that the ﬁt of the model is still only assessed to the observed data

since the missing data components of the observed datasets are ﬁlled in using

data augmentation (conditional on the model).

Basing the checks on completed datasets oﬀers several advantages. In prin-

ciple, as discussed throughout this chapter, interest most often is on the full-

data response model p(y | θ). Therefore, diagnostics on this model (as opposed

to the implied observed data model) will typically be of primary interest. As

such, the choices of test statistics discussed in the complete data setting in

Chapter 3 would be appropriate here. Second, model building under an as-

sumption of ignorability does not require explicit speciﬁcation of the missing

data mechanism, p(r | y

obs

, x, ψ). By using diagnostics based on completed

datasets, the replicated realizations of the full datavector y do not need to

diﬀerentiate the missing and observed components as speciﬁed by the response

indicator vector r;hence there is no need to specify or check p(r | y

obs

, x,ψ).

Forthese reasons, we conduct posterior predictive checks on the completed

datasets.

To assess the ﬁt of certain features of the model to the observed data,

appropriate test statistics need to be constructed. These test statistics have

the full data as their argument. For each MCMC sample, the test statistics

evaluated at the ‘observed’ data will have as their argument the observed

data y

obs

combined with the current realization of the posterior predictive

distribution of the missing data sampled in the data augmentation step. For

example, at iteration l,thefulldata would be (y

obs

, y

(l)

mis

), where within the

data augmented Gibbs sampler, we sample y

(l)

mis

from p(y

(l)

mis

| y

obs

, θ

(l)

). The

test statistics for the replicated data y

(l)

rep

will be sampled from the full-data

response model p(y

(l)

rep

| θ

(l)

)giventhecurrent value of the parameter at iter-

ation l, θ

(l)

(cf. posterior predictive distribution in (3.15)). This approach has

been implemented recently in several papers, including Ilk and Daniels (2007).

Choice of test statistics would be similar to those suggested in Chapter 3.

To be more speciﬁc, for some data summary T (·; θ)ofinterest, we compute

T (y

obs

, y

(l)

mis

; θ)andT (y

(l)

rep

; θ)atiteration l and compare the two realizations

at each iteration. The corresponding posterior predictive probability is de-

ﬁned as

p =



H(y

obs

, y

mis

, y

rep

,c; θ) p(y

rep

, y

mis

, θ | y

obs

) dθ dy

mis

rep

(6.18)

FURTHER READING 143

where

p(y

rep

, y

mis

, θ | y

obs

)=p(y

mis

| y

obs

, θ)p(y

rep

| θ)p(θ | y

obs

)

and

H(y

obs

, y

mis

, y

rep

,c; θ)=I



h(T (y

obs

, y

mis

; θ),T(y

rep

; θ)) >c



In Sections 7.2 and 7.5, we conduct checks for the growth hormone data and

CTQ II data, respectively.

For nonignorable missing data mechanisms, p(r | y, ω)isspeciﬁed and is

therefore part of the full-data model. We discuss posterior predictive checks

for these models in Chapter 8.

6.8 Further reading

Model speciﬁcation: dependence

Daniels and Pourahmadi (2002) propose prior distributions to shrink toward

structured GARP/IV models. Structured antedependence models (SAD) are

an alternative class of models that can be used to introduce structure into

acovariancematrix(N´u˜nez-Ant´on and Zimmerman, 2000) and are similar

to structured GARP/IV models. The log matrix parameterization is an al-

ternative unconstrained reparameterization of a covariance matrix that can

be used to introduce covariates (Leonard and Hsu, 1992; Chiu, Leonard, and

Tsui, 1996). For correlation matrices, Daniels and Pourahmadi (2007) have

been exploring alternative parameterizations for both structure and covari-

ates.

Another line of research on parsimoniously modeling covariance matrices

lets thedata itself ﬁnd a parsimonious structure. Wong et al. (2003) model de-

pendence parsimoniously (within the multivariate normal models in Example

2.3) by using Bayesian model averaging with priors that can ‘zero-out’ partial

correlations. Liu and Daniels (2007) use a similar approach in the context of

thejoint model in Section 6.6.3. Liechty et al. (2004) develop Bayesian mod-

els to parsimoniously estimate a correlation matrix through the construction

of appropriate priors. Chen and Dunson (2003) re-parametrize the random

eﬀects covariance matrix in order to ‘zero-out’ random eﬀects variances. This

parametrization is related to the modiﬁed Cholesky decomposition (Pourah-

madi, 2007).

Formisaligned measurement times in binary data, Su and Hogan (2007)

have developed continuous-time MTMs forlongitudinal binary data that allow

the dependence to depend on covariates.

144 INFERENCE UNDER MAR

Joint models

For the covariance structure for multivariate continuous longitudinal data

with misaligned times, additional structures can be found in Sy et al. (1997),

Sammeletal. (1999), Henderson et al. (2000), and Ferrer and McArdle (2003).

Directly speciﬁed models for binary data that do not use an underlying

multivariate normal or multivariate t latent structure, but provide marginal

logistic regressions for each component of Y , Y

ijk

,havebeen proposed. Fitz-

mauriceand Laird (1993) introduce a very general model for multivariate

(longitudinal) data based on the canonical log-linear parameterization. How-

ever, given the parameterization of this model, it is diﬃcult to parsimoniously

exploit the longitudinal correlation using serial correlation structures. Ilk and

Daniels (2007) constructed a model speciﬁcally for multivariate longitudinal

data, buildingonearlier work by Heagerty (1999, 2002) mixing both transition

(Markov) and random eﬀects structures; computations are complex though an

eﬃcient MCMC algorithm is proposed and implemented.

Indirectly speciﬁed conditional models (via random eﬀects or latent classes)

are another approach for modeling (multivariate) longitudinal binary data.

Relevant work includes Bandeen-Roche et al. (1997), Ribaudo and Thompson

(2002), Dunson (2003), and Miglioretti (2003).

There are many other approaches for general mixed longitudinal data that

have been proposed in the literature. Models for a set of mixed longitudinal

outcomes (more than just one continuous and binary longitudinal process)

have been developed without using an underlying normal latent structure, at

the cost of more complex computations. Some ofthesemodels use latent vari-

ables or latent classes to connect the processes, e.g., see Sammel et al. (1997),

Dunson (2003), and Miglioretti (2003). Instead of using latent variables, Lam-

bert et al. (2002) developed models for mixed outcomes using copulas (Nelson,

1999). Other authors have used the general location model for mixtures of

categorical and continuous outcomes (Fitzmaurice and Laird, 1997; Liu and

Rubin, 1998), but not speciﬁcally in the setting of multivariate longitudinal

data.

We did not discuss joint modeling of longitudinal and time to event data.

However, there is an extensive literatureonthistopic,bothinthecontext of

surrogate markers and for missing data where dropout is modeled as a time

to event process. See DeGruttola and Tu (1994), Faucett and Thomas (1996),

Wulfsohn and Tsiatis (1997), Wang and Taylor (2001), and Xu and Zeger

(2001), among others. For an early review of these methods, see Hogan and

Laird (1997b), and more recently, Tsiatis and Davidian (2004).

CHAPTER 7

Case Studies: Ignorable Missingness

7.1 Overview

In this chapter, we re-analyze several examples from Chapter 4 using all avail-

able data, under an ignorability assumption for the missing data. We also

analyze two additional examples. The ﬁrst uses the CTQ II smoking cessation

data (Section 1.4) to illustrate joint modeling and implementation of the aux-

iliary variable MAR assumption. The second uses the CD4 data from HERS

(Section 1.5) to illustrate modeling irregularly spaced longitudinal data with

missingness. In all these analyses we focusontheimportance of correctly

modeling/specifying the dependence under MAR.

7.2 Structured covariance matrices: Growth Hormone study

7.2.1 Models

The notation and model used here are the same as in Section 4.2. The distri-

bution of quadriceps strength (QS), Y

=(Y

)

,fortreatment group

= k is

| Z

= k ∼ N(µ

, Σ

), (7.1)

with µ

=(µ

,µ

)

and Σ

= Σ(η

), with η

containing the six

nonredundant parameters in the covariance matrix.

Aconvenientparameterization of the covariance matrix is based on the

GARP and IV parameters. To ﬁt these models inWinBUGS,were-parameterize

the multivariatenormalmodel as

| Z

= k ∼ N(µ

,σ

)

| y

= k ∼ N(β

+ φ

21,k

,σ

)

| y

= k ∼ N(β

+ φ

31,k

+ φ

32,k

,σ

where, for k =1,...,4,

=(φ

21,k

,φ

31,k

,φ

32,k

)

is the set of GARP for each treatment and σ

=(σ

,σ

)aretheIV

for each treatment. We set η

=(φ

, σ

). Using this re-parameterized model,

145

146 CASE STUDIES: IGNORABLE MISSINGNESS

the marginal means µ

can be computed recursively as

= β

j−2,k

j−1



l=1

jl,k

,j=2, 3.

7.2.2 Priors

Conditionally conjugate priors as described in Daniels and Pourahmadi (2002)

were used fortheGARP,

∼ N(0, 10

I),k=1,...,4.

We used diﬀuse normal priors on the intercepts (β

,β

)andtruncated uni-

forms on the square root of the IV parameters,

∼ N(0, 10

),j=0, 1; k =1,...,4,

∼ U(0, 100),j=1,...,3; k =1,...,4.

7.2.3 MCMC details

For all the mo dels, we ran three chains, each with 55, 000 iterations with a

burn-in of 5000 iterations. The chains mixed very well with minimal autocor-

relation.

Table 7.1 Growth hormone trial: posterior mean of the GARP/IV parameters of the

covariance matrices by treatment. The entries on the main diagonal are the innova-

tion variances (IV) and below the main diagonal are the generalized autoregressive

parameters (GARP).

EG G







697

.98 563

.45 .65 241













552

.90 176

.26 .61 91







EP P







741

.89 199

.21 .59 82













622

.74 203

−.01 .78 154







GROWTH HORMONE STUDY 147

7.2.4 Model selection and ﬁt

We ﬁrst ﬁt a model with Σ

distinct and unstructured for each treatment (k),

which we will call covariance model (1).Theposterior means of the GARP/IV

parameters η

for covariance model (1) are given inTable7.1. Examination of

this table suggestsmostoftheGARP/IV parameters do not vary substantially

across the treatments. As a result, we ﬁt two more parsimonious models:

covariance model (2) {η

= η : k =1,...,4} (common Σ across treatments)

and covariance model (3),

(σ

,φ

31,k

,φ

32,k

)=(σ

,φ

) k =1, 2, 3, 4

(σ

,σ

)=(σ

,σ

) k =2, 3, 4

21,k

= φ

k =1, 2, 3,

which allowed individual GARP/IV parameters to vary across (subsets) of

the treatments. Covariance model(3)waschosenbasedon our examination

of the estimates from covariance model (1) in Table 7.1. We use DIC

with

θ =(µ

, β

, φ

, {1/σ

: j =1,...,3; k =1,...,4})

to compare the ﬁt of the models. These results, which appear in Table 7.2,

support covariance model (3).

We also conducted two posterior predictive checks based on the residuals

to assess the ﬁt of covariance model (3). The posterior predictive probability

based on the multivariate version of Pearson’s χ

given in (3.19) was .65. The

one based on largest deviation in the empirical cdf’s, given in (3.22), was .40.

Both measures suggest the model ﬁts well.

Table 7.2 Growth hormone trial: DIC

for the threecovariance models.

Model DIC Dev(

θ) p

(1) 2674 2599 37

(2) 2667 2631 18

(3) 2652 2609 21

7.2.5 Results and comparison with completers-only analysis

Posterior means and 95% credible intervals for the mean parameters in co-

variance model (3) are given in Table 7.3; the posterior means and credible

intervals were similar under the three covariance models (not shown). In com-

paring these results to the complete case results in Table 4.3, we point out

148 CASE STUDIES: IGNORABLE MISSINGNESS

Table 7.3 Growth hormone trial: posterior means and 95% credible intervals of the

mean parameters for eachtreatmentassuming covariance model (3).

Month

Treatment 0 6 12

EG 69 (62, 77) 82 (71, 94) 81 (70, 92)

G68(61, 76) 66 (58, 74) 65 (58, 73)

EP 66 (58, 74) 81 (73, 90) 73 (65, 80)

P65(58, 73) 62 (55, 70) 63 (56, 70)

Table 7.4 Growth hormone trial: posterior probabilities that each of the pairwise

diﬀerences at month 12 is greater than zero under covariance model (3) and the

independence model.

Treatment Group Contrast

Model EG-G EG-EP EG-P G-EP G-P EP-P

(3) .99 .90 1.00 .08 .68 .97

Indep. 1.00 .98 1.00 .07 .52 .94

afewinteresting diﬀerences. For example, the credible intervals are narrower

under ignorability; this is expected as we are now using all the data. The diﬀer-

ence in the posterior mean of the 12-month mean in treatment EG illustrates

the potentially strong bias from the completers-only analysis; the posterior

mean was88(74, 103) in the completers-only analysis vs. 81 (70, 92) in the

MAR analysis). This is related to the fact that at baseline, there were large

diﬀerences between the mean of completers and the mean of all the subjects

(cf. Tables 4.3and7.3).

Posterior probabilities that each of the pairwise diﬀerences of the month 12

across the treatments are greater than zero are given in Table 7.4. The poste-

rior probabilities under ignorability showed diﬀerences from the completers-

only analysis (cf. the results in the caption of Figure 4.1). For example, the pos-

terior probabilities for the diﬀerence between treatment EG and EP changed

from .97 for the completers-only analysis to .90 for the MARanalysis. There

were also important diﬀerences between the three covariance models. For ex-

SCHIZOPHRENIA TRIAL 149

ample, in comparing covariance models (1) and (2) under ignorability, the

posterior probability for EG vs. EP was .90 and .80, respectively.

To further illustrate the importance of modeling dependence under MAR,

we ﬁt an independence model, Σ(η

)=σ

.TheDICforthis model was

2956 and the eﬀective number of parameters was 16. Based on the DIC, this

model ﬁt considerably worse than the dependence models (cf. Table 7.2).

The largest diﬀerence from the dependence models was seen in the month 12

mean for treatment EG. Under the independence model, this mean was 7 to 9

units higher than under the dependence models. This large diﬀerence was due

to the data augmentation step under the independence model not using the

observed values at month 0 and month 6 to ‘ﬁll-in’ the month 12 values for the

dropouts. This diﬀerence also led to more extreme (but incorrect) posterior

probabilities for the diﬀerences between month 12 means for treatment EG

and the other three treatments (see Table 7.4). It is also interesting to note

that the posterior distribution of the month 12 means under the independence

model are essentially equivalent to the posterior distribution under models

that allow dependence based on the completers only data.

7.2.6 Conclusions

The covariance model that allowed the individual GARP/IV parameters to

vary by treatment, covariance model (3), was the preferred model for this

analysis (as measured by the DIC) and also seemed to provide a good ﬁt

itself as measured by the posterior predictive checks. Mean QS at month 12

was signiﬁcantly higher on EG vs. the other three treatments, with posterior

probabilities ranging from .90 to 1.00. EP was signiﬁcantly higher than P

with a diﬀerence of 10 (ft.-lbs. of torque) and a posterior probability of .97.

An analysis of this data considering MNAR and sensitivity analysis can be

found in Chapter 10.

7.3 Normal random eﬀects model: Schizophrenia trial

7.3.1 Models and priors

The notation, models, and priors used here are the same as in Section 4.3.

Recall the goal of this 6-week trial is to compare the mean change in BPRS

scores (measure of schizophrenia severity) from baseline to week 6 between

the four treatment groups. For the full data Y

=(Y

,...,Y

)

,weassume

| b

= k ∼ N(x

(β

+ b

),σ

| Ω(φ

),Z

= k ∼ N(0, Ω(φ

)),

where the jth row of x

, x

is an orthogonal quadratic polynomial and Z

∈

{1, 2, 3, 4} is treatment group.