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

**|**HHS Author Manuscripts**|**PMC2875790

Formats

Article sections

- Summary
- 1 Introduction
- 2 Models for Longitudinal Data
- 3 Inference under Missing Data
- 4 Application
- 5 Discussion
- References

Authors

Related links

Biom J. Author manuscript; available in PMC 2010 August 1.

Published in final edited form as:

PMCID: PMC2875790

NIHMSID: NIHMS188151

Models for longitudinal data are employed in a wide range of behavioral, biomedical, psychosocial, and health-care related research. One popular model for continuous response is the linear mixed-effects model (LMM). Although simulations by recent studies show that LMM provides reliable estimates under departures from the normality assumption for complete data, the invariable occurrence of missing data in practical studies renders such robustness results less useful when applied to real study data. In this paper, we show by simulated studies that in the presence of missing data estimates of the fixed-effect of LMM are biased under departures from normality. We discuss two robust alternatives, the weighted generalized estimating equations (WGEE) and the augmented WGEE (AWGEE), and compare their performances with LMM using real as well as simulated data. Our simulation results show that both WGEE and AWGEE provide valid inference for skewed non-normal data when missing data follows the missing at random (MAR), the most popular missing data mechanism for real study data.

Models for longitudinal data are employed in a wide range of behavioral, biomedical, psychosocial, and health-care related research studies. One popular modeling paradigm is the latent variable, or random-effect, based approach for addressing correlated responses arising from such study data (Demidenko, 2004; Fitzmaurice et al. 2004; Raudenbush and Bryk, 2002). For example, for continuous response, the linear mixed-effects model (LMM) is one of the most widely used methods for modeling longitudinal data. By using random-effects to account for correlations across multiple within-subject outcomes, the LMM extends the classic multiple linear regression for modeling cross-sectional data to a general longitudinal data setting.

One major drawback of LMM as well as other random-effects based models is their dependence on distribution assumptions for inference. In recent years, many software packages have implemented robust methods to improve the validity of inference in the presence of departures from assumed parametric models. For example, simulation studies of Maas and Hox (2004; 2005) have found that estimates of the fixed-effects, or population parameters, of LMM are quite robust when model errors depart from the assumed normal distribution. However, these simulation studies have all been carried out under complete data. In most real studies, missing data invariably occurs and the robustness of estimates from LMM requires reassessment in the presence of missing data.

In this paper, we show by simulated data that the robustness of fixed-effect estimates of LMM is compromised under non-random missing data when the error term of the model deviates from the normality assumption. To address this key limitation, we investigate the performance of two robust alternatives, the weighted generalized estimating equations (WGEE) and augmented WGEE (AWGEE), for inference for non-normal longitudinal data. In Section 2, we briefly review the parametric LMM and the distribution-free formulation of this model. In Section 3, we discuss inference by the two dueling paradigms under missing data. In Section 4, we use simulated data to compare the performance of the various approaches for fitting non-normal longitudinal data under non-random missing data. In Section 5, we give our concluding remarks.

Over the past three decades, studies in biomedical and behavioral sciences have evolved from simple cross-sectional study designs to modern day longitudinal trials. As longitudinal study designs use subjects as their own controls, they provide a unique opportunity to study changes of outcomes of interest over time, causal effects and disease progression, in addition to providing more power for assessing treatment differences. In this section, we briefly review the two most popular approaches for modeling longitudinal data with continuous response.

First, consider a relatively simple longitudinal study design with only two assessments, i.e., the so-called pre-post design. Let *n* denote the number of subjects and *y _{it}* some continuous outcome of interest from the ith subject at time

$${y}_{it}={\beta}_{0}+{\beta}_{1}{I}_{\{t=2\}}+{b}_{i}+{\epsilon}_{it},\phantom{\rule{0.2em}{0ex}}{b}_{i}\sim N(0,{\sigma}_{b}^{2}),\phantom{\rule{0.3em}{0ex}}{\epsilon}_{it}\sim N(0,{\sigma}^{2}),\phantom{\rule{0.3em}{0ex}}1\le i\le n,\phantom{\rule{0.3em}{0ex}}1\le t\le 2,$$

(1)

where *I*_{{·}} denotes a set indicator, and *N* (*μ, σ*^{2}) denotes a normal with mean *μ* and variance *σ*^{2}. In the LMM above, since E (*y _{it}*) =

Let y* _{i}* = (

$$V=\text{Var}({\text{y}}_{i})=\left(\begin{array}{cc}{\sigma}^{2}+{\sigma}_{b}^{2}& {\sigma}_{b}^{2}\\ {\sigma}_{b}^{2}& {\sigma}^{2}+{\sigma}_{b}^{2}\end{array}\right)=({\sigma}^{2}+{\sigma}_{b}^{2}){C}_{2}(\rho ),$$

(2)

where
$\rho ={\sigma}_{b}^{2}/({\sigma}_{b}^{2}+{\sigma}^{2})$ is known as the intraclass correlation (ICC) and *C _{m}* (

More generally, the LMM for a longitudinal data with *n* assessments has the following form:

$${y}_{it}={\mathbf{\text{x}}}_{it}^{\top}\beta +{\mathbf{\text{z}}}_{it}^{\top}{\mathbf{\text{b}}}_{i}+{\epsilon}_{i},\phantom{\rule{0.3em}{0ex}}{\mathbf{\text{b}}}_{i}\sim N(\mathbf{0},{\mathbf{\sum}}_{b}),\phantom{\rule{0.3em}{0ex}}{\epsilon}_{i}\sim N(\mathbf{0},{\sigma}^{2}{\mathbf{\text{I}}}_{m}),\phantom{\rule{0.3em}{0ex}}1\le i\le n,\phantom{\rule{0.3em}{0ex}}1\le t\le m,$$

(3)

where x* _{it}* = (1,

$${V}_{i}=\text{Var}({\mathbf{\text{y}}}_{i}|{\mathbf{\text{x}}}_{i},{\mathbf{\text{z}}}_{i})={\mathbf{\text{Z}}}_{i}{\mathbf{\sum}}_{b}{\mathbf{\text{Z}}}_{i}^{\top}+{\sigma}^{2}{\mathbf{\text{I}}}_{m}.$$

(4)

By setting
${\mathbf{\text{x}}}_{it}^{\top}={\left(1,{I}_{\{t=2\}}\right)}^{\top}$ and **z*** _{it}* = 1 in (4), we immediately obtain (2).

Maximum likelihood (ML) is the most popular inference procedure for LMM (Demidenko, 2004; Raudenbush and Bryk, 2002). One major drawback of ML estimate is the dependence on the parametric assumptions; if data does not follow the normality assumptions in (1), model estimates may become biased or inconsistent. In recent years, many packages have adopted the sandwich variance estimate to address this issue (Goldstein, 1995; Rasbash et al., 2000; Raudenbush and Bryk, 2002). In this case, these procedures essentially yield variance estimates equivalent to those from the generalized estimating equations (GEE), which we review next.

Since the seminal work of Liang and Zeger (1986), the generalized estimating equations (GEE) approach has been widely used as an alternative for modeling longitudinal data. By modeling the marginal mean of the response at each assessment time, GEE eliminates both layers of distribution assumptions for the random effect and error term, thereby providing consistent estimates regardless of data distributions and the complexities of structure of correlated responses.

Within our context, consider the following linear model:

$$\text{E}({y}_{it}|{\mathbf{\text{x}}}_{it})={\mu}_{it}={\mathbf{\text{x}}}_{it}^{\top}\mathit{\beta},\phantom{\rule{0.3em}{0ex}}1\le i\le n,\phantom{\rule{0.3em}{0ex}}1\le t\le m.$$

(5)

In (5), only the marginal mean at each time *t* is specified, which models the fixed-effect of LMM in (1) at time *t*. Under GEE, inference is based on the following score-like vector equation:

$${\mathbf{\text{W}}}_{n}(\mathit{\beta})=\sum _{i=1}^{n}{\mathbf{\text{W}}}_{ni}(\mathit{\beta})=\sum _{i=1}^{n}{\mathbf{\text{G}}}_{i}({\mathbf{\text{x}}}_{i}){\mathbf{\text{S}}}_{i}(\mathit{\beta})=\sum _{i=1}^{n}{\mathbf{\text{G}}}_{i}({\mathbf{\text{x}}}_{i})({\mathbf{\text{y}}}_{i}-{\mathit{\mu}}_{i})=0,$$

(6)

where *μ** _{i}* = (

$${\mathbf{\text{G}}}_{i}({\mathbf{\text{x}}}_{i})={\mathbf{\text{D}}}_{i}{\mathbf{\text{V}}}^{-1}(\alpha ),\phantom{\rule{0.3em}{0ex}}{\mathbf{\text{D}}}_{i}=\frac{\partial}{\partial \beta}{\mu}_{i},\phantom{\rule{0.3em}{0ex}}1\le i\le n,$$

(7)

where **V** (** α**) denotes a working variance matrix parameterized by

The phrase working variance is used to emphasize the fact that **V** (** α**) is not necessarily the true variance matrix of

In addition to *β*, (6) also depends on ** α**, though we have suppressed this dependence to highlight the fact that (6) is the equation for estimating

$${\mathbf{\sum}}_{\beta}={\mathbf{\text{B}}}^{-1}\text{E}\left({\mathbf{\text{G}}}_{i}{\mathbf{\text{S}}}_{i}{\mathbf{\text{S}}}_{i}^{\top}{\mathbf{\text{G}}}_{i}^{\top}\right){\mathbf{\text{B}}}^{-\top},\phantom{\rule{0.3em}{0ex}}\mathbf{\text{B}}=\text{E}\left(\frac{{\partial}^{\top}}{\partial \mathit{\beta}}{\mathbf{\text{W}}}_{ni}\left(\mathit{\beta},\mathit{\alpha}\right)\right).$$

A consistent estimate of **Σ*** _{β}* is obtained by substituting consistent estimates in place of the respective quantities above. For example, we can estimate

In longitudinal studies, missing data are inevitable; subjects may simply quit or they may not show up at follow-up visits. We characterize the impact of missing data on model estimates through assumptions or missing data mechanisms, which allow us to ignore the multitude of reasons for missing data and focus on addressing their impact on estimation of model parameters. The missing completely at random (MCAR) assumption models a class of missing data that does not affect model estimates when completely ignored. For example, in a treatment study, missing data resulting from patient's relocation or scheduling conflict falls into this category. However, MCAR is not a plausible model when missing data are associated with treatment interventions such as patients' deteriorated or improved health conditions due to treatment. By modeling the occurrence of missing data as a function of observed responses prior to the assessment point, the missing at random (MAR) assumption addresses this class of treatment related or response-dependent missing data.

Within the longitudinal study setting in Section 2, define a missing (or rather, observed) data indicator as:

$${\mathit{\text{r}}}_{it}=\{\begin{array}{ll}1\hfill & \text{if}\phantom{\rule{0.2em}{0ex}}{y}_{it}\phantom{\rule{0.2em}{0ex}}\text{is observed}\hfill \\ 0\hfill & \text{if}\phantom{\rule{0.2em}{0ex}}{y}_{it}\phantom{\rule{0.2em}{0ex}}\text{is missing}\hfill \end{array},\phantom{\rule{0.3em}{0ex}}{\mathbf{\text{r}}}_{i}={\left({r}_{i1},\dots ,{r}_{im}\right)}^{\top},\phantom{\rule{0.3em}{0ex}}1\le t\le m,\phantom{\rule{0.3em}{0ex}}1\le i\le n.$$

We assume no missing data at baseline *t* = 1 such that *r _{i}*

Let **y*** _{i}* = (

The joint density function, *f* (y* _{i}*,

$$f({\mathbf{\text{y}}}_{i},{\mathbf{\text{r}}}_{i}|{\mathbf{\text{x}}}_{i})=f({\mathbf{\text{y}}}_{i}|{\mathbf{\text{x}}}_{i})\phantom{\rule{0.3em}{0ex}}f({\mathbf{\text{r}}}_{i}|{\mathbf{\text{y}}}_{i},{\mathbf{\text{x}}}_{i}).$$

(8)

Under MAR, the distribution of r* _{i}* depends only on the observed response,
${\mathbf{\text{y}}}_{i}^{o}$, and thus:

$$f({\mathbf{\text{r}}}_{i}|{\mathbf{\text{y}}}_{i},{\mathbf{\text{x}}}_{i})=f({\mathbf{\text{r}}}_{i}|{\mathbf{\text{y}}}_{i}^{o},{\mathbf{\text{y}}}_{i}^{m},{\mathbf{\text{x}}}_{i})=f({\mathbf{\text{r}}}_{i}|{\mathbf{\text{y}}}_{i}^{o},{\mathbf{\text{x}}}_{i}).$$

(9)

It follows from (8) and (9) that:

$$\begin{array}{ll}f({\mathbf{\text{y}}}_{i}^{o},{\mathbf{\text{r}}}_{i}|{\mathbf{\text{x}}}_{i})& =\int f({\mathbf{\text{y}}}_{i}^{o},{\mathbf{\text{y}}}_{i}^{m}|{\mathbf{\text{x}}}_{i})\phantom{\rule{0.2em}{0ex}}f({\mathbf{\text{r}}}_{i}|{\mathbf{\text{y}}}_{i}^{o},{\mathbf{\text{x}}}_{i})\phantom{\rule{0.2em}{0ex}}d{\mathbf{\text{y}}}_{i}^{m}\\ & =f({\mathbf{\text{r}}}_{i}|{\mathbf{\text{y}}}_{i}^{o},{\mathbf{\text{x}}}_{i},{\mathit{\text{\theta}}}_{\mathit{\text{y}}|\mathit{\text{r}}})\int f({\mathbf{\text{y}}}_{i}^{o},{\mathbf{\text{y}}}_{i}^{m}|{\mathbf{\text{x}}}_{i})\phantom{\rule{0.2em}{0ex}}d{\mathbf{\text{y}}}_{i}^{m}\\ & =f({\mathbf{\text{y}}}_{i}^{o}|{\mathbf{\text{x}}}_{i},{\mathit{\text{\theta}}}_{y})\phantom{\rule{0.2em}{0ex}}f({r}_{i}|{y}_{i}^{o},{\mathbf{\text{x}}}_{i},{\mathit{\text{\theta}}}_{y|r}).\end{array}$$

(10)

If *θ** _{y}* and

$$l(\mathit{\text{\theta}})=\sum _{i=1}^{n}log\phantom{\rule{0.2em}{0ex}}(f({\mathbf{\text{y}}}_{i}^{o}|{\mathbf{\text{x}}}_{i},{\mathit{\text{\theta}}}_{\mathit{\text{y}}}))+\sum _{i=1}^{n}log\phantom{\rule{0.2em}{0ex}}(f({\mathbf{\text{r}}}_{i}|{y}_{i}^{o}|{\mathbf{\text{x}}}_{i},{\mathit{\text{\theta}}}_{y|r}))={l}_{1}({\mathit{\text{\theta}}}_{y})+{l}_{2}({\mathit{\text{\theta}}}_{y|r}).$$

Thus, inference about *θ** _{y}* can simply be based on the log-likelihood

Most packages provide inference about the parameters of interest *θ** _{y}* based on maximizing the likelihood function

In the presence of missing data, we may apply the GEE in (6) to the observed responses, i.e.,

$${\mathbf{\text{W}}}_{n}(\mathit{\beta})=\sum _{i=1}^{n}{\mathbf{\text{W}}}_{ni}(\mathit{\beta})=\sum _{i=1}^{n}{\mathbf{\text{G}}}_{i}({\mathbf{\text{x}}}_{i}){\mathbf{\text{R}}}_{i}{\mathbf{\text{S}}}_{i}=\sum _{i=1}^{n}{\mathbf{\text{G}}}_{i}({\mathbf{\text{x}}}_{i}){\mathbf{\Delta}}_{i}\left({\mathbf{\text{y}}}_{i}-{\mathit{\mu}}_{i}\right)=0,$$

(11)

where **Δ_{i}** = diag

To illustrate the basic idea for modification, consider the relatively simple, pre-post design, with a homogeneous sample. We are interested in estimating the mean response at pre- and post-assessment, ** μ** = E (y

$${\mathbf{\text{W}}}_{n}(\mathit{\mu})=R{(\mathit{\alpha})}^{-1}\left(\begin{array}{ll}1\hfill & 0\hfill \\ 0\hfill & {r}_{i2}\hfill \end{array}\right)\phantom{\rule{0.3em}{0ex}}\left[\sum _{i=1}^{n}\left(\begin{array}{l}{y}_{i1}-{\mu}_{1}\hfill \\ {y}_{i2}-{\mu}_{2}\hfill \end{array}\right)\right]=\mathbf{0}.$$

(12)

Solving the equations above for ** μ** yields:

$$\widehat{\mathit{\mu}}={\left({\widehat{\mu}}_{1},{\widehat{\mu}}_{2}\right)}^{\top}={\left(\frac{1}{n}\sum _{i=1}^{n}{y}_{i1},\frac{1}{{n}_{2}}\sum _{i=1}^{n}{r}_{i2}{y}_{i2}\right)}^{\top},\phantom{\rule{0.3em}{0ex}}{n}_{2}=\sum _{i=1}^{n}{r}_{i2}.$$

(13)

If the missingness of *y _{i}*

Under MAR, the missingness of *y _{i}*

$${\pi}_{i2}=\text{P}({r}_{i2}=1|{y}_{i1},{y}_{i2})=\text{P}({r}_{i2}=1|{y}_{i1}).$$

This probability *π _{i}*

$${\mathbf{\text{W}}}_{n}(\mathit{\mu})=\mathbf{\text{R}}{(\mathit{\alpha})}^{-1}\left(\begin{array}{ll}1\hfill & 0\hfill \\ 0\hfill & \frac{{r}_{i2}}{{\pi}_{i2}}\hfill \end{array}\right)\phantom{\rule{0.2em}{0ex}}\left[\sum _{i=1}^{n}\phantom{\rule{0.2em}{0ex}}\left(\begin{array}{l}{y}_{i1}-{\mu}_{1}\hfill \\ {y}_{i2}-{\mu}_{2}\hfill \end{array}\right)\right]=\mathbf{0}.$$

(14)

It is readily checked that E(**W**_{n}(** μ**)) =

$${\widehat{\mu}}_{2}=\frac{1}{n}\sum _{i=1}^{n}\frac{{r}_{i2}}{{\pi}_{i2}}{y}_{i2}{\to}_{p}\text{E}\phantom{\rule{0.2em}{0ex}}\left(\frac{{r}_{i2}}{{\pi}_{i2}}{y}_{i2}\right)=\text{E}\phantom{\rule{0.2em}{0ex}}\left[{y}_{i2}\frac{1}{{\pi}_{i2}}E({r}_{i2}|{y}_{i1},{y}_{i2})\right]=\text{E}\phantom{\rule{0.2em}{0ex}}\left[{y}_{i2}\frac{1}{{\pi}_{i2}}E({r}_{i2}|{y}_{i1})\right]={\mu}_{2},\phantom{\rule{0.2em}{0ex}}$$

where →* _{p}* denotes convergence in probability.

By comparing (12) with (14), it is seen that the latter differs from the former only in the definition of **Δ*** _{i}*. By carrying this modification over to a general setting with

$${\pi}_{it}=\text{P}({r}_{it}=1|{\mathbf{\text{x}}}_{i},{\mathbf{\text{y}}}_{i}),\phantom{\rule{0.3em}{0ex}}{\mathbf{\Delta}}_{it}=\frac{{r}_{it}}{{\pi}_{it}},\phantom{\rule{0.3em}{0ex}}{\mathbf{\Delta}}_{i}={\text{diag}}_{t}({\mathbf{\Delta}}_{it}),\phantom{\rule{0.3em}{0ex}}1\le t\le m,\phantom{\rule{0.3em}{0ex}}1\le i\le n.$$

(15)

It is again readily checked that E [**G*** _{i}* (x

To use WGEE, we must know or have estimates of *π _{it}*. In some cases, subject dropout is created by study design and

The WGEE discussed in Section 3.2.1 depends on the model for missing data in (15). In most studies, *π _{it}* are unknown and must be modeled and estimated. If such a model is misspecified, the WGEE estimate may be inconsistent. In applications, reliable models may also exist for directly relating the missing response to the observed ones and other covariates. The augmented WGEE (AWGEE) is developed to take advantage of this additional source of modeling information to ensure valid inference when the model for

To illustrate, consider again the pre-post study design. Suppose that we can predict *y _{i}*

$$\text{E}({y}_{i2}|{y}_{i1})={\eta}_{0}+{\eta}_{1}{y}_{i1},\phantom{\rule{0.3em}{0ex}}1\le i\le n.$$

(16)

Then, we can estimate *μ _{2}* without using WGEE by
${\stackrel{\sim}{\mu}}_{2}={\sum}_{i=1}^{n}\left({\widehat{\eta}}_{0}+{\widehat{\eta}}_{1}{y}_{i1}\right)/n$. This new estimate is consistent if (16) is a correct model, since

$${\stackrel{\sim}{\mu}}_{2}={\widehat{\eta}}_{0}+{\widehat{\eta}}_{1}\frac{1}{n}\sum _{i=1}^{n}{y}_{i1}{\to}_{p}{\eta}_{0}+{\eta}_{1}\text{E}({y}_{i1})=\text{E}\phantom{\rule{0.2em}{0ex}}[\text{E}({y}_{i2}|{y}_{i1})]=E({y}_{i2})={\mu}_{2}.$$

By combining both the prediction model in (16) and the WGEE in (14), we obtain an augmented WGEE to estimate *μ* as follows:

$${\mathbf{\text{W}}}_{n}(\mathit{\mu})=\sum _{i=1}^{n}\mathbf{\text{R}}{(\mathit{\alpha})}^{-1}\left({\mathbf{\Delta}}_{i}{\mathbf{\text{S}}}_{i}-{\mathbf{\Delta}}_{i}^{c}{\stackrel{\sim}{S}}_{i}\right),\phantom{\rule{0.3em}{0ex}}{\mathbf{\Delta}}_{i}^{c}={\mathbf{\Delta}}_{i}-{\mathbf{\text{I}}}_{2},$$

(17)

$${\stackrel{\sim}{S}}_{i}={\left({y}_{i1}-{\mu}_{1},{\widehat{y}}_{i2}-{\mu}_{2}\right)}^{\top},{\widehat{y}}_{i2}={\eta}_{0}+{\eta}_{1}{y}_{i1}.$$

It is readily checked that
$\text{E}\left[{\mathbf{\text{G}}}_{\mathbf{\text{i}}}({\mathbf{\text{x}}}_{\mathbf{\text{i}}})\phantom{\rule{0.2em}{0ex}}\left({\mathbf{\Delta}}_{\mathbf{\text{i}}}{\mathbf{\text{S}}}_{i}-{\mathbf{\Delta}}_{i}^{c}\stackrel{\sim}{{\mathbf{\text{S}}}_{i}}\right)\right]=\mathbf{0}$ if either (15) or (18) or both are correct. Thus, the AWGEE above yields consistent estimates of ** μ** if at least one of these models is correct. Further, when both models are correct, the AWGEE estimate from (17) may also be more efficient than the WGEE estimate (Robins et al, 1995; Tsiatis, 2006).

The above is readily extended to a more general setting where the prediction model in (16) also involves other baseline covariates. Let u* _{i}* be a set of baseline variables including

$$\text{E}({y}_{i2}|{\mathbf{\text{u}}}_{i})={\eta}_{0}+{\mathbf{\text{u}}}_{i}^{\top}{\eta}_{1},\phantom{\rule{0.3em}{0ex}}1\le t\le 2,\phantom{\rule{0.3em}{0ex}}1\le i\le n.$$

(18)

The AWGEE is defined by

$${\mathbf{\text{W}}}_{n}(\mathit{\beta})=\sum _{i=1}^{n}{\mathbf{\text{W}}}_{ni}(\mathit{\beta})=\sum _{i=1}^{n}{\mathbf{\text{G}}}_{i}\left({\mathbf{\text{x}}}_{i}\right)\phantom{\rule{0.2em}{0ex}}\left({\mathbf{\Delta}}_{i}{\mathbf{\text{S}}}_{i}-{\mathbf{\Delta}}_{i}^{c}{\stackrel{\sim}{S}}_{i}\right)=\mathbf{0},$$

(19)

where
${\widehat{y}}_{i2}={\eta}_{0}+{\text{u}}_{i}^{\top}{\eta}_{1}$. To ensure consistent estimation for the regression model, we assume a surrogacy-type assumption, [*y _{i}*

Although feasible in principle, it is more complex to implement AWGEE for a general longitudinal study with more than two assessments. For example, for *m* = 3, we need to consider two missing data patterns when predicting missing *y _{i}*

Under MCAR, **r_{i}** are independent of x

As noted earlier, missing data in longitudinal trials often occur as the result of subject dropout due to deteriorated/improved health conditions and other related conditions, exhibiting the so-called monotone missing data pattern (MMDP). The structured patterns under MMDP make it possible to model *π _{it}* in most studies.

Under MMDP, if *y _{it}* is observed at time

$$\begin{array}{l}{\mathbf{\text{H}}}_{it}=\{{\stackrel{\sim}{\mathbf{\text{x}}}}_{is},{\stackrel{\sim}{\mathbf{\text{y}}}}_{is};1\le s\le t-1\},\phantom{\rule{0.3em}{0ex}}{\stackrel{\sim}{\mathbf{\text{x}}}}_{it}={({x}_{i1}^{\top},\dots ,{x}_{i1(t-1)}^{\top})}^{\top},\\ {\stackrel{\sim}{\mathbf{\text{y}}}}_{it}={({y}_{i1},\dots ,{y}_{i(t-1)})}^{\top},\phantom{\rule{0.3em}{0ex}}2\le t\le m.\end{array}$$

The subset **H_{it}** contains all observed data prior to time

$${\pi}_{it}=\text{P}({r}_{it}=1|{\mathbf{\text{x}}}_{i},{\mathbf{\text{y}}}_{i})=\text{P}({r}_{it}=1|{\mathbf{\text{H}}}_{it})=\text{P}({r}_{it}=1|{\stackrel{\sim}{x}}_{it},{\stackrel{\sim}{y}}_{it}).$$

(20)

Thus, under MMDP and MAR, *π _{it}* are a function of observed data only, making it possible to estimate these selection probabilities.

Let *p _{it}* = E (

$${\pi}_{it}=\prod _{s=2}^{t}{p}_{is},\phantom{\rule{0.3em}{0ex}}2\le t\le m,\phantom{\rule{0.3em}{0ex}}1\le i\le n.$$

(21)

Thus, we can estimate *π _{it}* by modeling the

$$\text{logit}\phantom{\rule{0.2em}{0ex}}({p}_{it}\phantom{\rule{0.2em}{0ex}}({\alpha}_{t}))=\text{logit}\phantom{\rule{0.2em}{0ex}}(\text{E}({r}_{it}=1|{r}_{i(t-1)}=1,{\mathbf{\text{H}}}_{it}))={\mathbf{\xi}}_{t}^{\top}{\stackrel{\sim}{w}}_{it},\phantom{\rule{0.3em}{0ex}}2\le t\le m,$$

(22)

where
${\stackrel{\sim}{\mathbf{\text{w}}}}_{it}={\left(1,{\stackrel{\sim}{\mathbf{\text{x}}}}_{it}^{\top},{\stackrel{\sim}{\mathbf{\text{y}}}}_{it}^{\top}\right)}^{\top}$ and
${\stackrel{\sim}{\phantom{\rule{0.2em}{0ex}}\mathbf{\text{w}}}}_{i}={\left({\stackrel{\sim}{\mathbf{\text{w}}}}_{i2}^{\top},\dots ,{\stackrel{\sim}{\mathbf{\text{w}}}}_{im}^{\top}\right)}^{\top}$. For each *t*, we can estimate **ξ**_{t} by maximum likelihood or GEE conditional on the observed **_{it}** at

For AWGEE, we again consider the pre-post study design. In this case, we can readily estimate ** η** in (18) using GEE, where
$\eta ={({\eta}_{0},{\eta}_{1}^{\top})}^{\top}$. With an estimate , we can predict

For inference based on WGEE, let **Δ_{i}**, (

$${\sum}_{i=1}^{n}{\mathbf{\text{q}}}_{i}(\mathit{\xi})={\sum}_{i=1}^{n}{\left({\mathbf{\text{q}}}_{i2}^{\top},\dots ,{\mathbf{\text{q}}}_{im}^{\top}\right)}^{\top}=0,{\mathbf{\text{q}}}_{it}=\frac{\partial}{\partial {\mathit{\xi}}_{t}}\left\{{r}_{i(t-1)}\phantom{\rule{0.2em}{0ex}}\left[{r}_{it}log({p}_{it})-(1-{r}_{it})log(1-{p}_{it})\right]\right\}.$$

(23)

Now, let

$$\mathbf{\text{B}}=\text{E}\left(\frac{{\partial}^{\top}}{\partial \mathit{\beta}}{\mathbf{\text{W}}}_{ni}\right),\phantom{\rule{0.3em}{0ex}}\mathbf{\text{C}}=\text{E}\left(\frac{{\partial}^{\top}}{\partial \mathit{\xi}}{\mathbf{\text{W}}}_{ni}\right),\phantom{\rule{0.3em}{0ex}}\mathbf{\text{F}}=\text{E}\left(\frac{{\partial}^{\top}}{\partial \mathit{\xi}}{\mathbf{\text{q}}}_{i}(\mathit{\xi})\right),$$

where **W*** _{ni}* is defined in (11). Then, as shown in Appendix A, under
$\sqrt{n}$-consistency of

$${\sum}_{\mathit{\beta}}={\mathbf{\text{B}}}^{-1}(\text{Var}({\mathbf{\text{W}}}_{ni})+\mathbf{\Phi}){\mathbf{\text{B}}}^{-\top},\mathbf{\Phi}={\mathbf{\text{CF}}}^{-1}\text{Var}({\mathbf{\text{q}}}_{i}){\mathbf{\text{F}}}^{-\top}{\mathbf{\text{C}}}^{\top}-\text{E}\left({\mathbf{\text{W}}}_{ni}{\mathbf{\text{q}}}_{i}^{\top}{\mathbf{\text{F}}}^{-\top}{\mathbf{\text{C}}}^{\top}\right)-{\left[\text{E}\left({\mathbf{\text{W}}}_{ni}{\mathbf{\text{q}}}_{i}^{\top}{\mathbf{\text{F}}}^{-\top}{\mathbf{\text{C}}}^{\top}\right)\right]}^{\top}.$$

(24)

In (24), **Φ** accounts for the variability of estimated **. We can estimate ****Σ*** _{β}* by substituting consistent estimates in place of the respective quantities.

For AWGEE inference, we also need to estimate *η* for the prediction model in (18). Using GEE, the vector estimating equation is given by:

$${\sum}_{i=1}^{n}{\mathbf{\text{g}}}_{i}(\eta )={\sum}_{i=1}^{n}\frac{\partial}{\partial \eta}\left\{{r}_{i2}\left[{y}_{i2}-\left({\eta}_{0}+{\mathbf{\text{u}}}_{i}^{\top}{\mathit{\eta}}_{2}\right)\right]\right\}=\mathbf{0}.$$

(25)

Let
$\widehat{\phi}={\left({\widehat{\phi}}_{1}^{\top},{\widehat{\phi}}_{2}^{\top}\right)}^{\top}$ with *ϕ*_{1} = ** ξ** and

$${\sum}_{i=1}^{n}{\mathbf{\text{s}}}_{i}={\sum}_{i=1}^{n}{\left({\mathbf{\text{q}}}_{i}^{\top}({\varphi}_{1}),{\mathbf{\text{g}}}_{i}^{\top}({\varphi}_{2})\right)}^{\top}=\mathbf{0},\phantom{\rule{0.3em}{0ex}}{\mathbf{\text{s}}}_{i}={\left({\mathbf{\text{q}}}_{i}^{\top}({\varphi}_{1}),{\mathbf{\text{g}}}_{i}^{\top}({\varphi}_{2})\right)}^{\top}.$$

The AWGEE estimate ** is also asymptotically normal under
$\sqrt{n}$-consistency of ****, with the asymptotic variance having the same form as in (24) except for substituting ****s*** _{i}* for

We illustrate our considerations with both real and simulated data. We first present an application to data from a longitudinal study in depression research and then investigate the performance of the approach with small to moderate sample sizes by simulation. In all the examples, we set the statistical significance at *α* = 0.05. All analyses are carried out using a code we have developed for implementing the proposed approach using the R software platform (Free Software Foundation, 1999). This code is available from the author upon request.

In a study on geriatric depression and associated medical comorbidities for old primary care patients, 744 subjects were enrolled from private practices and University-affiliated clinics in general internal medicine, geriatrics, and family medicine in Monroe County, New York (Lyness et al., 2007). All patients age 65 years and older who presented for care on selected days and were capable of giving informed consent. Enrolled subjects underwent semi-structured interviews, administered by trained raters in the subjects' homes or in research offices at the UR Medical Center. The raters' subject interviews included assessments of cognition, functional status, and psychopathology, the latter including the Structured Clinical Interview for DSM-IV (SCID) (Spitzer et al. 1986). Interviews and chart reviews were conducted at study intake, and again at one- and two-year follow-up time points.

In geriatric research, overall functional disability is of particular importance, as it reflects both the mental and physical health conditions of the individual. Primary measures of overall functional status include the Instrumental Activities of Daily Living (IADL), Physical Self-Maintenance Scales (PSMS), Global Assessment of Functioning (GAF), and the Karnofsky Performance Status Scale (KPSS) (Lawton MP and Brody, 1969, Karnofsky DA and Burchenal JH, 1949, Ware JE, Jr. and Sherbourne CD). For illustration purposes, we analyzed the change of IADL from baseline to one-year follow-up, as this measure assesses instrumental activities such as shopping or using the telephone and is particularly popular in geriatric research. Further, we only included the baseline value as a predictor when modeling the missingness of this outcome as well as the outcome itself at the follow-up using the respective logistic (20) and linear (18) models.

Of the 744 enrolled, 468 completed the IADL at the one-year follow-up. Shown in Table 1 are the estimates of the intercept and slope from the fitted logistic regression for modeling the missingness and the linear regression for modeling the outcome of IADL as a function of its baseline value at the follow-up. The baseline IADL was significant in both models, indicating that it did predict the occurrence of missing IADL as well as the outcome itself at the follow-up. Note that the negative sign of the estimate of the coefficient for baseline IADL in the logistic model indicates that the subjects with lower baseline IADL were more likely to come for assessment at the one-year follow-up. As lower IADL is associated with poorer functioning status, the observed sample at the follow-up visit seemed to be biased towards those with more severe overall functional disability at baseline.

Estimates of parameters of (1) logistic regression for modeling missingness at one-year follow-up, and (2) linear model for predicting IADL at one-year follow-up for the study on geriatric depression and associated medical comorbidities.

We fit the LMM in (1) and the distribution-free alternative in (5) to examine the change of IADL from baseline to the one-year follow-up. Shown in Table 2 are the estimates of the intercept *β*_{0} and slope *β*_{1} for the respective models under the different inference procedures. As the estimates of *β*_{1} were positive across the board, the mean IADL increased at the follow-up visit, indicating better functioning status for the old primary care patients in this observational study.

Estimates of parameters of (1) linear mixed-effects model (ML), and (2) distribution-free linear model (GEE, WGEE and AWGEE) for change of IADL from baseline to one-year follow-up for the study on geriatric depression and medical comorbidities.

However, the magnitude of the estimate of *β*_{1} did vary substantially — not only between the models, but also across the different procedures within the same distribution-free linear model. The WGEE and AWGEE yielded quite similar estimates, with the latter AWGEE also providing improved efficiency, as indicated by smaller asymptotic standard errors. The GEE performed poorly, with a whopping 50% downward bias, as compared to its counterparts WGEE and AWGEE estimates. For the between-model comparison, the ML estimate of *β*_{1} from the fitted LMM also incurred a downward bias, albeit with a much smaller magnitude relative to the GEE estimate. The downward bias in both cases was consistent with the fact that those assessed at the follow-up visit represented a subgroup with more severe overall functional disability at baseline.

Shown in the Figure is the normal-based Q-Q plot of the conditional residuals obtained from the estimated fixed and random effects for the fitted LMM model (Nobre and Singer, 2007). The plot indicates clearly that the residuals did not follow a normal distribution, which may explain the difference in the estimates of the fixed effects between the parametric (LMM) and distribution-free models (WGEE and AWGEE).

Given the discrepant estimates between LMM and WGEE (AWGEE) for the real study data in 4.2, we conducted a simulation study with a pre-post study design to investigate this issue further. We considered two non-normal distributions for the model error term of the LMM: a rescaled central chi-square with one degree of freedom and a uniform between −1 and 1. Since the results are quite similar, we only discuss and report the results from the chi-square-distributed error. To examine the performance of the models under small, moderate and large samples, we performed the simualation study with three sample sizes: *n* = 50, 100 and 2,000. All simulations were performed with a Monte Carlo sample of 1,000 using the R software (R Development Core Team, 2007).

We considered the pre- and post-treatment design and simulated the outcome according to the LMM in (1) by setting *β*_{0} = *β _{1}* = 1 and

$$\text{logit}\phantom{\rule{0.2em}{0ex}}({\pi}_{i2})=\text{logit}\phantom{\rule{0.2em}{0ex}}(\text{P}\phantom{\rule{0.2em}{0ex}}({r}_{i2}=1|{y}_{i1}))={\xi}_{0}+{\xi}_{1}{y}_{i1},\phantom{\rule{0.4em}{0ex}}1\le i\le n.$$

(26)

We set ξ_{0} = 0.5 and ξ_{1} = 1.2 to create about 25% missing response *y _{i}*

Under (1), it is readily checked that regardless of the distributions for *b*_{i} and * _{it}* (see Appendix)

$$\text{E}\phantom{\rule{0.2em}{0ex}}({y}_{i2}|{y}_{i1})={\eta}_{0}+{\eta}_{1}{y}_{i1},\phantom{\rule{0.4em}{0ex}}{\eta}_{0}={\beta}_{0}\phantom{\rule{0.2em}{0ex}}\left(1-\frac{{\sigma}_{b}^{2}}{{\sigma}_{b}^{2}+{\sigma}^{2}}\right)+{\beta}_{1},\phantom{\rule{0.4em}{0ex}}{\eta}_{1}=\rho =\frac{{\sigma}_{b}^{2}}{{\sigma}_{b}^{2}+{\sigma}^{2}},\phantom{\rule{0.4em}{0ex}}1\le i\le n.$$

(27)

The above was used to predict missing *y _{i}*

We considered the null *H*_{0} : *β*_{0} = *β*_{1}, i.e., the mean at post-treatment is twice that at pretreatment. We tested *H*_{0} using the Wald statistic,
${\mathit{\text{Q}}}_{n}^{2}=n{\widehat{\mathit{\beta}}}^{\top}{\mathbf{\text{K}}}^{\top}{(\mathbf{\text{K}}{\widehat{\sum}}_{\beta}{\mathbf{\text{K}}}^{\top})}^{-1}\mathbf{\text{K}}\widehat{\mathit{\beta}}$, which has an asymptotic central
${\chi}_{1}^{2}$ distribution. We estimated the type I error rate *α* based on the distribution of
${Q}_{n}^{2}$ from 1, 000 Monte Carlo (MC) replications,
$\widehat{\alpha}=({\sum}_{j=1}^{1000}{I}_{\{{Q}_{nj}\ge q0.95\}})/1000$, where
${Q}_{nj}^{2}$ denotes the value of
${Q}_{n}^{2}$ from the *j*th MC replication and *q*_{0.95} the 95th percentile of
${\chi}_{1}^{2}$. The maximum likelihood (ML) inference about *β* = (*β*_{0},*β*_{1})^{┬} for LMM was obtained from the LME procedure in R, while the WGEE and AWGEE estimates for the distribution-free alternative in (5) were computed based on the asymptotic results in Section 3.2.

Shown in Table 3 are the estimates of *β* and associated asymptotic standard errors averaged over 1,000 MC replications obtained from ML, GEE, WGEE and AWGEE for the respective models, with a within-subject correlation *ρ* = 0.1 (or
${\sigma}_{b}^{2}=2/9$ and *σ*^{2} = 2). The results confirmed that the baseline mean *β _{0}* is consistently estimated by all four procedures. For

Averaged estimates of *β* over 1,000 Monte Carlo replications along with asymptotic standard errors (s.e.) and type I error rates *α* for sample size 50, 100, 2000, with about 25% missing data at post-treatment based on ML for linear mixed-effects **...**

Under the correct weight function, WGEE performed well. When the incorrect constant weight function was used, WGEE yielded biased estimates, while the AWGEE estimates remained close to the true value of *β*_{1} under the correct prediction model. However, AWGEE did not show any significant gain in efficiency; in fact, the estimate of the slope *β*_{1} had a larger standard error under AWGEE than under WGEE across small sample sizes *n*.

To further investigate the relative efficiency between WGEE and AWGEE, we replicated the above analysis by increasing the within-subject correlation. For example, shown in Tables 4 and and55 are the estimates from one such replicated analysis with *ρ* = 0.3 (or
${\sigma}_{b}^{2}=1$ and *σ*^{2} = 2 and 0.6 (or
${\sigma}_{b}^{2}=1.5$ and *σ*^{2} = 1), respectively. It is seen that all conclusions above remain the same, except for the relative efficiency between WGEE and AWGEE under the correct weight function and prediction model. As *ρ* increased to 0.3, AWGEE have smaller standard errors than its counterpart WGEE when *n* = 2000. At *ρ* = 0.6, not only did AWGEE show smaller standard errors than WGEE across all sample sizes, their differences also widened as compared to those with *ρ* = 0.1.

Averaged estimates of *β* over 1,000 Monte Carlo replications along with asymptotic standard errors (s.e.) and type I error rates *α* for sample size 50, 100, 2000, with about 25% missing data at post-treatment based on ML for linear mixed-effects **...**

We investigated the two primary modeling strategies for longitudinal continuous response, the linear mixed-effects model (LMM) and the distribution-free linear model, with respect to their performance under missing data, and illustrated our considerations using real as well as simulated study data. Our results show that LMM and the GEE procedure for the distribution-free alternative generally yield biased estimates under MAR when the normality assumption for LMM is violated. Further, as indicated by the simulation results, the standard errors of these estimates do not increase to reflect model misspecification, making inference prone to misleading findings. Thus, when modeling longitudinal data, it is important to test the MCAR assumption as discussed in Section 3.2 before applying any of the models and inference procedures considered. If the null of MCAR is rejected, WGEE and/or AWGEE should be considered, unless there is strong evidence to support the use of the alternative linear mixed-effects model.

Our simulation results also indicate that the gain in efficiency by AWGEE over WGEE depends on the magnitude of the within-subject correlation *ρ*. Within the context of the particular simulation model considered, AWGEE is less efficient than WGEE for small sample size under small *ρ* such as 0.1. But, AWGEE edged out WGEE to be a more efficient procedure as *ρ* increased to 0.6.

This research was supported in part by NIH grants R01-DA012249, R21-AG023956, UL1 RR024160 and R24-MH071604. We thank Ms. Bliss-Clark for her careful proofreading of the manuscript, and an anonymous, an Associate Editor and Editor Prof. Leonhard Held for their constructive comments that led to a substantial improvement in the presentation of this research.

Under (1), we have

$${\rho}_{\mathit{\text{yb}}}=\text{Corr}\phantom{\rule{0.2em}{0ex}}({y}_{i1},{b}_{i})=\frac{\text{Cov}\phantom{\rule{0.2em}{0ex}}({y}_{i1},{b}_{i})}{\sqrt{\text{Var}\phantom{\rule{0.2em}{0ex}}({y}_{i1})\phantom{\rule{0.2em}{0ex}}\text{Var}\phantom{\rule{0.2em}{0ex}}({b}_{i})}}=\sqrt{\frac{{\sigma}_{b}^{2}}{{\sigma}_{b}^{2}+{\sigma}^{2}}}.$$

(28)

Further, the above holds regardless of the distributions for *b _{i}* and

$${b}_{i}={\varphi}_{0}+{\varphi}_{1}{y}_{i1}+{\delta}_{i},\phantom{\rule{0.4em}{0ex}}{\delta}_{i}\sim (0,{\sigma}_{\delta}^{2}),\phantom{\rule{0.4em}{0ex}}1\le i\le n,$$

Where $(0,{\sigma}_{\delta}^{2})$ denotes a distribution with mean 0 and variance ${\sigma}_{\delta}^{2}$. It follows from (28) and the relationship between linear regression coefficients and correlation that

$${\varphi}_{1}={\rho}_{\mathit{\text{yb}}}\sqrt{\frac{\text{Var}\phantom{\rule{0.2em}{0ex}}({b}_{i})}{\text{Var}\phantom{\rule{0.2em}{0ex}}({y}_{i1})}}=\sqrt{\frac{{\sigma}_{b}^{2}}{{\sigma}_{b}^{2}+{\sigma}^{2}}}\sqrt{\frac{{\sigma}_{b}^{2}}{{\sigma}_{b}^{2}+{\sigma}^{2}}}=\frac{{\sigma}_{b}^{2}}{{\sigma}_{b}^{2}+{\sigma}^{2}}.$$

Also, since

$$0=\text{E}\phantom{\rule{0.2em}{0ex}}({b}_{i})={\varphi}_{0}+{\varphi}_{1}\text{E}\phantom{\rule{0.2em}{0ex}}({y}_{i1})={\varphi}_{0}+{\varphi}_{1}{\beta}_{0},$$

it follows that ${\phi}_{0}=-{\beta}_{0}{\sigma}_{b}^{2}/({\sigma}_{b}^{2}+{\sigma}^{2})$.

Under (1) and regardless of the distributions for *b _{i}* and

$$\text{E}\phantom{\rule{0.2em}{0ex}}({y}_{i2}|{y}_{i1})={\beta}_{0}+{\beta}_{1}+\text{E}\phantom{\rule{0.2em}{0ex}}({b}_{i}|{y}_{i1})={\beta}_{0}+{\beta}_{1}+{\varphi}_{0}+{\varphi}_{1}{y}_{i1}.$$

By substituting the expressions of _{0} and _{0} into the above and combining the coefficients, we obtain (27).

1. Demidenko E. Mixed Models: Theory and Applications. New York: Wiley; 2004.

2. Diggle PJ, Heagerty P, Liang KY, Zeger SL. Analysis of Longitudinal Data. 2nd. Oxford University Press; 2002.

3. Fitzmaurice GM, Laird NM, Ware JH. Applied Longitudinal Analysis. New York: Wiley; 2004.

4. Goldstein H. Multilevel Statistical Models. Edward Arnold; London; Halsted, New York: 1995.

5. Karnofsky DA, Burchenal JH. The clinical evaluation of chemotherapeutic agents in cancer. In: MacLeod CM, editor. Evaluation of Chemotherapeutic Agents. New York: Columbia; 1949.

6. Kowalski J, Tu XM. A GEE approach to modeling longitudinal data with incompatible data formats and measurement error: application to HIV immune markers. Journal of the Royal Statistical Society, Series C. 2002;51:91–114.

7. Kowalski J, Tu XM. Modern Applied U Statistics. New York: Wiley; 2007.

8. Lawton MP, Brody EM, editors. Gerontologist. Vol. 9. 1969. Assessment of older people; self-maintaining and instrumental activities of daily living; pp. 179–186. [PubMed]

9. Lyness JM, Niculescu A, Tu XM, Reynolds CF, III, Caine ED. The relationship of medical comorbidity to depression in older primary care patients. Psychosomatics. 2007;47:435–439. [PubMed]

10. Maas CJM, Hox JJ. The influence of violations of assumptions on multilevel parameter estimates and their standard errors. Computational Statistics & Data Analysis. 2004;46:427–440.

11. Maas CJM, Hox JJ. Sufficient Sample Sizes for Multilevel Modeling. Methodology. 2005;1:86–92.

12. Nobre JS, Singer JM. Residual analysis for linear mixed models. Biometrical Journal. 2007;49:863–875. [PubMed]

13. Prentice RL. Surrogate endpoints in clinical trials: Definition and operational criteria. Statistics in Medicine. 1989;8:431–440. [PubMed]

14. Rasbash J, Browne W, Goldstein H, Yang M, Plewis I, Healy M, Woodhouse G, Draper D, Langford I, Lewis T. A User's Guide to MLwiN. Multilevel Models Project. University of London; London: 2000.

15. Raudenbush SW, Bryk AS. Hierarchical Linear Models. 2nd. Sage; Thousand Oaks, CA: 2002.

16. R Development Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing; Vienna, Austria: 2007. URL http://www.R-project.org.

17. Robins JM, Rotnitzky A, Zhao LP. Analysis of semiparametric regression models for repeated outcomes in the presence of missing data. Journal of the American Statistical Association. 1995;90:106–121.

18. Spitzer RL, Williams JB, Gibbon M. Structured Clinical Interview for DSM-III-R (SCID) New York: Biometrics Research; 1986.

19. Tsiatis AA. Semiparametric Theory and Missing Data. New York: Springer; 2006.

20. Ware JE, Sherbourne CD. The MOS 36-item short-form health survey (SF-36). I. Conceptual framework and item selection. Medical Care. 1992;30:473–483. [PubMed]

21. Zeger SL, Liang KY. Longitudinal data analysis for discrete and continuous outcomes. Biometrics. 1986;42:121–130. [PubMed]

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