Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Evaluating the associations between the baseline and other exposure variables with the longitudinal trajectory when responses are measured with error
(USC Thesis Other)
Evaluating the associations between the baseline and other exposure variables with the longitudinal trajectory when responses are measured with error
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
Evaluating the Associations between the Baseline and Other Exposure
Variables with the Longitudinal Trajectory When Responses are Measured
with Error
by Zhanghua Chen
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(Biostatistics)
December 2013
Copyright 2013 Zhanghua Chen
ii
Dedication
I dedicate my dissertation work to my family and many friends. I would like to
express my deepest gratitude to my loving parents, Jianguo Zhang and Deying Chen.
Their support, encouragement, and sacrifice have sustained me throughout my life.
I would also like to give special thanks to my wonderful husband, Xueming Yu,
and my two lovely kids, Tianzhang Yu and Yunjin Zhang Chen, for being there for me
throughout the entire doctorate program. All of you have been my best cheerleaders.
iii
Acknowledgements
I would like to express my deepest gratitude to my advisor, Dr. Richard M.
Watanabe, for having me as a student. He has provided excellent guidance and an
atmosphere to conduct independent research. I have been fortunate learning from him
throughout these years and I am greatly indebted for all of his time, patience, and
guidance.
I would also like to express my deepest gratitude to Dr. Anny H. Xiang for all her
help and support throughout these years in terms of mentorship and research
opportunities. Most importantly, she has taught me the importance of being an
independent researcher. I would like to thank her for giving me the chance to prove
myself as a research assistant and for providing invaluable funding opportunities.
I would like to thank Dr. Daniel O. Stram for being on my committee and for
guidance on my dissertation topic. He has been generous with his insightful ideas and
advice. I also owe a huge debt of gratitude to his funding support during these years on
the genetic research.
I would like to thank Dr. Wendy J. Mack for her encouragement and insightful
advice on my dissertation research and writing.
I would like to thank Dr. Thomas A. Buchanan for agreeing to be on my
committee and serving as the external member. He has provided invaluable advice on the
biological aspects of my work.
iv
Table of Contents
Acknowledgements iii
List of Tables vii
List of Figures ix
Abstract xi
Chapter 1: Introduction 1
1.1 Statistical methods for analyzing longitudinal data 1
1.2 Measurement error 10
Chapter 2: Assessing the baseline prediction of the longitudinal trajectory when
responses are measured with error 24
2.1 Literature review 24
2.2 Commonly used mixed effects models for the baseline association and the
longitudinal trajectory 29
2.3 Theoretical derivation of bias correction formula for the baseline association
with the longitudinal trajectory using mixed effects models 32
2.4 Simulation study 38
Chapter 3: Assessing the baseline-adjusted association between the exposure
variable and the longitudinal trajectory when responses are measured with error 63
3.1 Literature review 63
3.2 Commonly used mixed effects model for the baseline-adjusted association
between the exposure variable and the longitudinal trajectory 67
v
3.3 Theoretical derivation of bias correction formula for baseline-adjusted
association between the exposure variable and the longitudinal trajectory 69
3.4 Simulation study 74
Chapter 4: Obtaining bias corrected p-values and 95% confidence intervals 100
4.1 Theoretical derivations of test statistics 100
4.2 Simulation results 102
Chapter 5: GDM cohort study 109
5.1 Introduction 109
5.2 Study design and methods 110
5.3 Simulation study 114
5.4 Results 116
5.5 Conclusion 119
Chapter 6: The association between intrauterine GDM exposure and offspring’s
BMI growth from age 5 to 11 129
6.1 Introduction 129
6.2 Study design and methods 130
6.3 Results 133
6.4 Conclusion 136
Chapter 7: Summary and future directions 145
7.1 Overview 145
7.2 Correcting measurement error with no additional data 147
7.3 Analysis of non-linear growth trajectories 149
vi
7.4 Adjustment for additional confounders 150
7.5 Exposure variables also contain measurement error 152
7.6 Conclusions 153
Bibliography 155
vii
List of Tables
Table 2.1: Simulation results from two, four and six follow-up visits 47
Table 2.2: Comparison of simulation results from Model I(a), Model II(a) and Byth’s
model 49
Table 2.3: Validation results from simulated data 50
Table 2.4: Validation results from simulated data with 20% data missing at random 51
Table 2.5: Comparison of bias-corrected estimates
of mixed effects model and linear
regression 52
Table 3.1: Comparison between λ and λ' at Different Correlation Coefficients (r) 83
Table 3.2: Simulation results of
from Two Follow-up Visits
a
Using Model I(b) 84
Table 3.3: Simulation results of
from Two Follow-up Visits
a
Using Model I(b) 86
Table 3.4: Simulation results of
from Two Follow-up Visits
a
Using Model II(b) 88
Table 3.5: Simulation results of
from Two Follow-up Visits
a
Using Model II(b) 90
Table 3.6: Validation results from two follow-up visits
a
using bias-correction formula
for Model I(b) and Modle II(b) 92
Table 3.7: Comparison of bias-corrected estimates
of mixed effects model and linear
regression 93
Table 4.1: Comparison of test statistics for
and
using delta method and
bootstrapping method 104
Table 5.1: Validation results of data simulated to represent glucose and ivGTT traits 123
Table 5.2: Baseline characteristics and rates of change during the follow-up 125
Table 5.3: Baseline prediction of longitudinal changes in oGTT and ivGTT traits 126
Table 5.4: Baseline-adjusted associations between baseline caloric intake and
longitudinal changes in oGTT and ivGTT traits 127
viii
Table 6.1: Comparison of characteristics of offspring exposed and unexposed to GDM in
utero for children entering the study at five years of age 139
Table 6.2: The association between baseline BMI and the longitudinal change in BMI
between the ages of 5 to 11 140
Table 6.3: The association between GDM exposure in untero and the rate of increase in
BMI between the ages of 5 to 11 after adjusting for baseline BMI, race, sex,
maternal age, parity, education, gestational week at delivery and birth
weight 141
ix
List of Figures
Figure 2.1: Bias of naï ve estimates (
) using Model I(a) and Model II(a)
against the proportion of variance in the total variance of observed baseline
value due to measurement error 54
Figure 2.2: Bias of naï ve estimates (
) using Model I(a) and Model II(a)
against the true underlying parameter of the association between baseline
value and longitudinal slop (
) 55
Figure 2.3: Bias of naï ve estimates (
) using Model I(a) against the inverse of
number of follow-up visits 56
Figure 2.4: Bias of naï ve estimates (
) using Model I(a), Model II(a) and
Byth’s model against the true variance of model residual (
) 57
Figure 2.5: Comparisons of biases in the naï ve estimate (
) and bias-corrected
estimate (
) using Model I(a) and Model II(a) over ME parameter ( )
settings (0 - 0.4) 58
Figure 2.6: Comparisons of biases in the naï ve estimate (
) and bias-corrected
estimate (
) using Model I(a) and Model II(a) over follow-up visits from 3
to 10 59
Figure 2.7: Comparisons of biases in the naï ve estimate (
) and bias-corrected
estimate (
) using Model I(a) and Model II(a) over ME parameter ( )
settings (0 - 0.4) when follow-up interval is varied from 8-16 months 60
Figure 2.8: Comparisons of biases in the naï ve the estimate (
) and the bias-corrected
estimate (
) using Model I(a) and Model II(a) over follow-up visits from 3
to 10 when follow-up interval is varied from 8-16 months 61
Figure 2.9: Comparisons of percent biases in the bias-corrected estimate (
) using
Model I(a) and Model II(a) over ME parameter ( ) settings when 50% data is
missing at random 62
Figure 3.1: Comparisons of biases in the naï ve estimate (
) and the bias-corrected
estimate (
) using Model I
B
and Model II
B
over different correlations (r)
between baselineY
0
and exposure Z 95
x
Figure 3.2: Comparisons of biases in the naï ve estimate (
) and the bias-corrected
estimate (
) using Model I
A
and Model II
A
over ME parameter ( ) settings
(0 - 0.4) 96
Figure 3.3: Comparisons of biases in the naï ve estimate (
) and the bias-corrected
estimate (
) using Model I
A
and Model II
A
over ME parameter ( ) settings
(0 - 0.4) when follow-up interval is varied from 8-16 months 97
Figure 3.4: Comparisons of biases in the naï ve estimate (
) and the bias-corrected
estimate (
) using Model I(b) and Model II(b) over ME parameter ( )
settings (0 - 0.4) when the exposure Z is categorical variable 98
Figure 3.5: Comparisons of biases in the bias-corrected estimate (
) using Model I
B
,
Model II
B
and Harrison’s method over different residual variations (
)
settings (0.005, 0.01, 0.02, 0.04, 0.06, 0.08, 0.1) 99
Figure 4.1: Comparisons of type I error rate of z-statistics of naï ve estimates (
and
) and bias-corrected estimates (
and
) using Model I(b) and Model
II(b) over ME parameter ( ) settings (0 - 0.4) 106
Figure 4.2: Comparisons of power of z-statistics of naï ve estimates (
and
) and
bias-corrected estimates (
and
) using Model I(b) and Model II(b)
over ME parameter ( ) settings (0 - 0.4) 107
Figure 4.3: Comparisons of power of z-statistics of naï ve estimates (
and
) and
bias-corrected estimates (
, and
) using Model I(b) and Model II(b)
over follow-up visits n=3-10 108
Figure 5.1: The association between the baseline caloric intakes and rates of change in 2-
hr glucose, S
I
, and DI during the follow-up 128
Figure 6.1: Offspring BMI growth between the ages of 5 to 11. Means and standard
errors of BMI at each visit are presented. n is the number of sample size at
each visit 142
Figure 6.2: Offspring BMI growth curves between the ages of 5 to 11 stratified by
median baseline BMI (16.2kg/m
2
). 143
Figure 6.3: Offspring BMI growth curves between the ages of 5 to 11 stratified by
children exposed and unexposed to GMD in utero. 144
xi
Abstract
In clinical research, the impact of baseline value of the response variable on its
longitudinal growth trajectory has been frequently discussed. Researchers are not only
interested in how the baseline value directly affects the longitudinal trajectory, but also
how baseline value and other correlated exposures jointly alter the trajectory. In this
dissertation, we will focus on one topic in this field: analysis of the associations between
baseline and other exposure variables with the longitudinal trajectory when responses are
subject to intra-individual variations such as biological random variations and random
measurement errors. Conventionally, researchers combined different sources of random
variations as a general topic of random measurement error problem. The theoretical work
of this dissertation was motivated by a study investigating the pathogenesis of diabetes
among women with prior gestational diabetes mellitus (GDM). Disposition index (DI) is
a derived measure from intravenous glucose tolerance test (ivGTT) to assess β-cell
function for insulin resistance (Bergman, 1989). DI is an important outcome measuring
-cell health, and the deterioration in DI indicates progression to type 2 diabetes. The
intra-individual variations in DI cannot be ignored due to subject daily variations, glucose
injection variations, glucose and insulin measurement errors, etc. Because of the intra-
individual variations in observed DI, the evaluation of potential factors (including
baseline DI) associated with the deterioration of DI over time can be incorrect. In order to
solve this estimation problem, we need to first assess the influence of intra-individual
variations in an association study and then derive a bias-correction formula to calibrate
the naï ve association estimate. In the following contents, we will treat all kinds of intra-
xii
individual variations as the general random measurement error. Because we will only
focus on continuous responses, such as DI, we will use linear mixed effects models to
analyze the associations between exposure variables and the longitudinal trajectory.
In Chapter 1, we first introduce a general review of traditional statistical methods
for analyses of longitudinal data and assessments of the associations between baseline
and other exposure variables with the longitudinal trajectory (1.1). The measurement
error problem will be described in 1.2. A broad review of various methods developed to
correct measurement error in covariates and outcomes for linear regression and mixed
effects models are compared for strengths and weaknesses.
Chapter 2 starts with a literature review of bias correction methods for the
estimation of the baseline association with the longitudinal trajectory for both studies
with single follow-up visit and multiple follow-up visits. The limitations of previous
methods are discussed. After the literature review, we propose the mixed effects model
including an interaction term between baseline and time as a fixed effect to assess the
baseline prediction of the longitudinal trajectory. The bias of naï ve association estimates
will be assessed and bias-correction formulae will be derived for our proposed models.
Finally, a simulation study will be conducted to assess the bias of naï ve estimates for our
proposed models and to validate the bias-correction formulae for those naï ve estimates.
Moreover, the simulation study is used to compare our new method to previous methods.
In Chapter 3, we first present a literature review of previous methods to correct
biased estimates for the baseline-adjusted association between the exposure variable and
the longitudinal trajectory. After the discussion of limitations of previous methods,
xiii
Chapter 3 extends the method developed in Chapter 2 and investigates the bias in the
naï ve estimate of baseline-adjusted association between the exposure of interest and the
longitudinal trajectory in a mixed effects model. A bias-correction formula is developed
for the naï ve estimate of the baseline-adjusted association between the exposure variable
and the trajectory. A simulation study is completed to assess the bias in the naï ve
estimates, to validate our bias-correction formulae, and to compare our method with other
published methods.
In Chapter 4, we developed test statistics of our previously developed bias-
corrected estimates using two different methods. A simulation study is conducted to
assess the type I error and power of our proposed test statistics of the bias-corrected
estimates.
Chapter 5 and Chapter 6 present two real study examples to apply our bias-
correction formulae to naï ve association estimates. Chapter 5 assesses the association
between baseline caloric intake and longitudinal trajectories of type 2 diabetes-related
traits after adjusting for the baseline value of the trait in the GDM cohort study (Xiang et
al., 2010). Chapter 6 investigates the impact of intrauterine GDM exposure on the
offspring BMI growth trajectory after adjusting for offspring baseline BMI using a
sample dataset arising from Kaiser Permanente.
Chapter 7 summarizes findings from previous chapters and discusses advantages
and limitations of our methods. We also propose some future directions to extend our
methods to a more general application.
1
Chapter 1
Introduction
1.1 Statistical methods for analyzing longitudinal data
1.1.1 Comparison among ANOVA, MANOVA, summary measures, mixed effects
model and GEE
Longitudinal studies have a defined feature in that repeated measurements are
taken from the same individual over time. Therefore, longitudinal studies allow the
investigation of the change of responses over time and associations with the change. The
goal of longitudinal studies in biomedical research is to investigate the factors which
have causal relationships with the development and persistence of disease or disease-
related biomarkers. Longitudinal data have a corresponding feature where measurements
from the same individual are positively correlated and the correlations often decrease
with longer time between measurements. Therefore, the analysis of longitudinal data
should consider the correlations of repeated measurements within individuals. Ignoring
the correlations within individuals leads to an overestimation of the variance of the
change in responses over time, thus incorrectly inflating the p-values of test statistics.
For continuous longitudinal outcomes, historical approaches involve univariate
repeated measures analysis of variance (ANOVA) (Fisher, 1925; Greenhouse and Geisser,
1959; Winer, 1971; Yates, 1935); multivariate repeated measures analysis of variance
(MANOVA) (Box, 1950; Cole and Grizzle, 1966; Geisser, 1963; Hand and Taylor, 1987);
summary measures such as the area under the curve (AUC) and growth curve models
2
(Box, 1950; Everitt, 1995; Matthews et al., 1990; Rowell and Walters, 1976; Wishart,
1938); linear mixed effects models (Fisher, 1918; Laird and Ware, 1982) and generalized
estimating equations (GEE) (Burton, Gurrin and Sly, 1998). Garrett et al. compared the
pros and cons of each of these methods (Garrett, Nan and James, 2004).
ANOVA and MANOVA: Common limitations of ANOVA and MANOVA are
that they require the covariates to be time-invariant and discrete factors, the study design
to be regular time spaced, and the data to be complete. Moreover, ANOVA and
MANOVA only compare group means in the analysis, ignoring the information of
individual-specific growth curves. One advantage of MANOVA over ANOVA is that
MANOVA allows a more general variance-covariance structure for repeated
measurements. ANOVA restrictively assumes a compound-symmetric structure among
successive measurements, which is not appealing for longitudinal data because the
correlations between repeated measurements are usually expected to decay with
increasing time between measurements.
Summary measures: Compared with ANOVA and MANOVA, summary
measures allow unbalanced design and data missing at random. However, summary
measures are limited in focusing on only one summary aspect of the data, thereby losing
the information on heterogeneity among individuals. Different individual profiles of
repeated measures can be represented in the same summary measure of responses.
Summary measures also require the covariate to be time-invariant.
Mixed model and GEE: Modern analysis methods for longitudinal data,
including mixed effects models and general estimating equations (GEE), handle
3
continuous and time-varing covariates. Both methods can be applied to unbalanced
designs where measurements are collected at different times or include randomly missing
measurements for each individual. Although previous publications showed the magnitude
of regression coefficients and inference tests using mixed effects models and GEE are
almost the same (Gardiner, Luo and Roman, 2009; Twisk, 2004), there is an argument
that GEE are not favored when the outcome is continuous because widely used statistical
software such as SAS restricts the scale parameter used in estimating the variance of
responses to be time-invariant. This restricts the variance of repeated measurements to be
consistent during the follow-up (Garrett et al., 2004).
For discrete longitudinal responses, generalized mixed effects models and GEE
are widely used. The choice of these two methods depends on the study endpoint.
Generalized mixed effects models are applied when the inference focuses on individual
prediction and needs to capture the heterogeneity across individuals. The GEE approach
is used when the inference is the entire population and the investigation is more interested
in a population-averaged effect (Gardiner et al., 2009; Garrett et al., 2004; Twisk, 2004).
1.1.2 Mixed effects model
By definition, a mixed effects model contains both fixed and random effects.
Fixed effects are population-averaged characteristics that apply to all individuals.
Random effects are subject-specific effects that are unique to each individual in the
population. Random effects models were first introduced by Fisher to study traits
following Mendelian inheritance between relatives (Fisher, 1918). Henderson (Henderson,
1975; Henderson et al., 1959) developed a best linear unbiased estimator (BLUE) of
4
fixed effects and Searle, Harville et al. (Harville, 1988, 1990; McLean, Sanders and
Stroup, 1991; Robinson, 1991; Searle, 1971) developed a best linear unbiased predictor
(BLUP ) for random effects. Based on Herderson’s notation for a mixed effects model:
, where Y is an vector of responses, X is the design matrix. B is a
vector of unknown fixed effects. Z is design matrix. U is a vector of
unknown random effects with E(U)=0. is an vector with E( )= 0, and
. Thus, and
. For this model, estimates of
and are obtained through likelihood-based methods. Classical maximum likelihood
(ML) estimators of B and V are from the same sample data, which cause the denominator
of
to be overestimated by the rank of B, thereby the entire
is underestimated. It
can be shown that
. In order to address this problem in estimating V,
restricted maximum likelihood (REML) uses the residuals (
) to estimate V. Hence,
REML is widely applied in mixed models especially for small samples(Corbeil and
Searle, 1976; Garrett et al., 2004). The restricted likelihood for the mixed effects model
can be expressed as:
,
where
. To obtain the estimates of fixed and random effects,
mixed model equations (Henderson, 1984) are generally applied:
.
Henderson shows the corresponding solutions are:
5
.
Laird first introduced the mixed effects model to longitudinal data (Laird and
Ware, 1982). In this publication, the mixed effects model was separated into two stages.
Stage one incorporated the population parameters, individual effects, and within-person
variation:
for each individual i, where is the vector of
unknown population parameters;
is the
known design matrix;
is the
vector of unknown individual effect;
is the
design matrix linking
to
individual response
;
is the within-individual variation with distribution of
.
The covariance matrix
is a
positive-definite covariance matrix. Stage two
represents the between-individual variation, in which
are distributed as and
are independent of each other and of
. The covariance matrix of G is a positive-
definite covariance matrix.
1.1.3 Baseline impact on longitudinal slope under mixed effects model
In longitudinal studies, the baseline value of the response variable may impact the
growth curve of responses over time (Gill et al., 1985; Halverson and RE., 1981). In
clinical studies, findings of baseline predictors of the longitudinal trajectory are of great
interest (Ironson et al., 2005; Sen, Jen and Djuric, 2007; Xiang et al., 2010). For studies
with only the baseline and one follow-up measurement, a change model using the linear
regression approach is generally applied to investigate the association between the
baseline value of responses and the longitudinal change over time (the change between
the baseline value and the follow-up value regressed on the baseline value) (Blomqvist,
6
1977; Jin, 1992). When the follow-up intervals are varied across subjects, the outcome of
the regression model is often expressed as a rate of change, which is calculated by the
change between the baseline value and the follow-up value divided by the follow-time
(Xiang et al., 2010). The regression coefficient associated with the baseline value
estimates the change association with the baseline.
For studies with multiple follow-up visits and continuous responses, linear mixed
effects models are usually performed where responses are modeled with a random
intercept and random slope of time effects. Generally, there are two approaches for mixed
effects models to estimate the baseline association with the longitudinal change. First,
responses over time are regressed on the fixed effect of the baseline value, the random
effect of follow-up time, and the interaction between the baseline value and time
(Gardiner, Johnson and Demirel, 2012; Garrett et al., 2004; Josephs et al., 2011; Konstan
et al., 2007; Wattmo et al., 2011). The outcome in the model can include responses at
each time point or the changes from the baseline value to each follow-up value. The
regression coefficient of the interaction term estimates the baseline prediction of the
longitudinal trajectory. In practice, there are two different ways of constructing the
dependent variable in the model. Some studies have included the baseline value in the
dependent variable (Gardiner et al., 2012; Josephs et al., 2011), and others have excluded
the baseline value from the dependent variable (Wattmo et al., 2011). In the study
investigating baseline predictors of the functional decline (assessed by Clinical Dementia
Rating Sum of Boxes scale (CDR-SB)) of behavioral variant frontotemporal dementia,
CDR-SB scores at each visit, including the baseline visit, were modeled as the dependent
7
variable. The baseline CDR-SB is analyzed for its association with the rate of change in
CDR-SB using a mixed effects model by estimating the interaction between baseline
CDR-SB and the follow-up time (Josephs et al., 2011). A similar method was used in the
study of predictive factors of eye functional progression in early glaucoma patients
(Gardiner et al., 2012). The rate of change in mean deviation (MD), a measure of eye
function, is calculated using measurements at all visits in a mixed effects model, and then
the rate of change in MD is regressed on the baseline MD. This two-stage modeling
strategy can be treated the same as a mixed effects model including the interaction term
between the baseline and time, and the baseline value is included in the dependent
variable. On the other hand, an Alzheimer’s disease study investigates whether a
baseline cognitive score predicts the long-term changes in cognitive scores by testing the
interaction between the baseline cognitive score and the follow-up time as a fixed effect
in the mixed effects model. The dependent variable is composed of successive cognitive
scores starting from the first follow-up visit, and then all subsequent assessments of each
patient (Wattmo et al., 2011). In another study of lung function, researchers used the
mixed effects model with the interaction term between baseline forced expiratory volume
in one second (FEV1) and the follow-up time to describe the association between
baseline FEV1 and the decline of FEV1 among children with cystic fibrosis (Konstan et
al., 2007). However, these studies do not describe the details about whether they include
baseline FEV1 in the model outcome or not. Taken together, mixed effect models with
the interaction term between baseline value and time are widely used to evaluate the
baseline prediction of the subsequent changes in responses. However, little attention has
8
been paid to whether baseline value should be included as the dependent variable or not.
It should be noted that models with the baseline in the dependent variable is an incorrect
model because the baseline value appears in both the dependent and the independent
variables, which can cause biases in parameter estimations. In our paper, we will discuss
more details of two modeling methods in Chapter 2, and show that we should always use
the model excluding the baseline from the dependent variable for real data analyses.
The second modeling approach includes regressing longitudinal responses on a
random intercept and a random effect of time. The covariance matrix can be obtained
from the random intercept and the random slope of time. The random intercept is used to
estimate the individual baseline value. The association between the change and baseline
is estimated by the ratio of the covariance between the intercept and the slope of time to
the variance of the intercept (Byth and Cox, 2005; Kenward and Roger, 2010). More
details of this approach will be described in 1.2.6 and in Chapter 2.
1.1.4 Adjustment for baseline value of the response variable under a mixed effects
model
When the research question is not restricted to the baseline prediction of the
longitudinal change, and is more concerned with exploring the association between the
other exposure variable and the longitudinal change over time, the baseline value of the
response variable is often a confounder to be adjusted for in the study design and
statistical analysis. For randomized experiments, various randomization methods have
been applied to achieve the balance of baseline characteristics among the study groups
(Laird and Fong, 1990; Su, 2011; Yoshioka, 1998). For other clinical studies, particularly
9
the observational studies, complete randomization is difficult to achieve by the study
design. Therefore, the baseline value in the statistical analysis is usually adjusted for.
When there is only one baseline and one follow-up visit in the study, the baseline value is
adjusted for as a covariate in the change model (Yanez, Kronmal and Shemanski, 1998)
using linear regression. When the study is designed to include multiple follow-up visits,
there are two ways to adjust for baseline values in the mixed effect model. Both are
extensions of two mixed effects models used for the baseline prediction of the
longitudinal change described in 1.1.3. In the first approach, the mixed effects model
regresses successive responses against fixed effects of time, the exposure of interest, the
baseline value, and two interaction terms (the baseline × time and the exposure × time),
as well as the random effect of time (Garrett et al., 2004). The regression coefficient of
the interaction (the exposure × time) estimates the association between the exposure
variable and the growth trajectory after adjusting for the baseline. The second approach
uses a mixed effects model including a random intercept and a random slope of the time
effect, a cross-sectional fixed effect of the exposure, and the interaction between the
exposure and time. The association between the baseline and the longitudinal trajectory is
estimated in the same way as Byth’s method (Byth and Cox, 2005). The baseline adjusted
association between the exposure variable and the longitudinal trajectory is estimated by
subtracting the regression coefficient of the interaction term (the exposure × time) by the
product of baseline association with the longitudinal trajectory, and with the exposure
variable (Harrison et al., 2009). More details of these methods will be illustrated in 1.2.7
and Chapter 3.
10
1.2 Measurement error
1.2.1 Definition of measurement error
The accurate measurement of variables is required for correct estimates of
baseline prediction of the longitudinal trajectory or the association between the exposure
variable and the trajectory adjusting for baseline. When there is measurement error in
either covariates or responses, the estimates of associations between the baseline and
exposure variables with the longitudinal trajectory are likely to be biased. The general
measurement error is defined as the within subject variations of several measurements
completed on the same quantity (Bland and Altman, 1996). There are two categories of
measurement error: random error and systematic error. Random measurement error is the
statistical fluctuation around the true value among a series of measurements due to
limitations in the precision (Fluke., 1994; Taylor, 1997). It should be noted that we use a
more general definition of random measurement error in our work, which also includes
individual biological variations. Naï ve estimates or statistical inferences on correlations
or relationships using measurements subject to random measurement error can lead to
biased estimates and incorrect test properties. The second type of measurement error is
called systematic measurement error, which is the occurrence of consistent discrepancies
in the same direction between the observed measures and the true underlying value. Thus,
the inaccurate measurements are always too high or too low (Baird, 1995; Bevington and
Robinson, 1992; Fluke., 1994; ISO, 1993; Taylor, 1997). For categorical variable,
measurement error is also referred to as misclassification.
11
1.2.2 Statistical approaches to address measurement error
Most previous studies on the measurement error problem in longitudinal studies
have been focused on error-prone variables in covariates. When there is only one baseline
and one follow-up visit for each individual, the bias correction for naive estimates can be
developed based on measurement error models for linear regressions. Fuller, Carrol et al.
and John have provided a broad review of measurement error models for linear and
nonlinear models (Carroll et al., 2006; Fuller, 1987; John, 2010).
Least square estimates with known reliability ratio: Fuller studied the bias of
least square estimates for simple and multiple linear regression when the covariates in the
model are measured with error (Fuller, 1987). A bias-corrected estimate was developed
based on a least square estimate with known reliability ratio. The reliability ratio is the
variance of the true variable without error (
) divided by the variance of the observed
variable measured with error (
). For continuous variables, an additive measurement
error model is widely used, in which the variance of the true variable equals the variance
of the observed variable subtracted by the variance of the measurement error (
).
Although the variance of measurement error is unknown, it can be estimated using extra
data such as replications of all sample or a subsample:
, where
is the number of replicates for ith individual;
is the number of individuals with replicate data; and
is
the sample variance of replicated measurements of each individual in the replication
sample (Chan and Mak, 1979; Yanez et al., 1998).
12
Induced model approach: For nonlinear regression or generalized linear
regression, an induced model followed by the quasi-likelihood method or regression
calibration approach is usually performed to correct the bias of naï ve estimates (Carroll et
al., 2006; Fuller, 1987; John, 2010). Let Y be the outcome, X be the true covariates and
W be the observed covariates.
, where is the function regressing Y on X, and is
the density function of X given W. An induced model requires W to be a surrogate of X, ie.
. The density function involves a measurement error model,
which is usually a Berkson error model. The distinction between a classical measurement
error model and a Berkson error model is that the former model represents the observed
value given the true value (W|X), but the latter model represents the true value given the
observed value (X|W). Another important assumption for the induced model is that the
true value (X) must be a random variable rather than a fixed value.
The quasi-likelihood method and regression calibration approaches are both
based on the induced model (Carroll et al., 2006; John, 2010). For the quasi-likelihood
method, the induced model , where , is
the parameter in the function of expected X given W,
. Let
. There are three general steps for the quasi-likelihood method: First, find an
expression of ; Second, use the external data such as replication or
validation data to estimate the parameters η; Third, fit the regression model with observed
data using the expression in the first step and the estimated parameter . Because the
13
standard quasi-likelihood method usually requires complex integrations, it is
computationally unattractive.
Regression calibration is a special case within the quasi-likelihood method and
is simpler to apply (Carroll et al., 2006; John, 2010). In the regression calibration
approach, the imputed true covariate can be estimated from the observed covariates
through a calibration equation (linear regression etc.) using external information such as
validation data where both true and observed covariates are known,
. If a
validation sample is not available, a true covariate can also be estimated from the
replication sample under the assumption of the measurement error model. For example,
when the classical additive error is assumed,
where
is the mean estimate of sample ,
and
are estimated variance matrices of true
covariate , and measurement error from the replication sample, then the regression-
calibrated parameter (
) is estimated from the model of the outcome against the
expected true covariates given observed covariates. The power loss compared to using
the true X is due to the fact that the variance of expected true covariates given observed
covariates ( ) is always smaller than the variance of true covariates
(X): . Moreover, the imprecision of estimating
will further lead to imprecision in estimating
.
When the induced model assumptions are not applicable, other correction
methods can be employed such as the moment-based method, estimating equation
approaches (including score function methods) (Carroll et al., 2006; Nakamura, 1990),
simulation extrapolation (SIMEX) (Cook and Stefanski, 1994; Stefanski and Cook, 1995;
14
Wang et al., 1998), and Bayesian methods (Dominici, Zeger and Samet, 2000; Espino-
Hernandez, Gustafson and Burstyn, 2011; Fox and Glas, 2003; Gustafson, 2004).
The moment-based method derives the bias-corrected estimators based on the
summary statistics (mean, variance, covariance, etc.) obtained from the observed data
(Fuller, 1987; John, 2010). The moment-based method is applicable to cases when the
true covariates (X) are random variables or fixed values. In the normal structure model,
where measurement error is normal, distributed with constant variance, the moment-
based estimator is same as the least-square estimator (John, 2010).
Estimating equations: Like the moment-based method, estimating equations
allow the covariates X to be both random and fixed (Carroll et al., 2006; Nakamura,
1990). The classical measurement error model can be assumed in estimating equations.
Under maximum likelihood, estimating equations are derived from the first derivative of
the log-likelihood in terms of the parameter of interest which links the outcome and
covariates:
, where
is the outcome;
is true covariates; and
is the parameter of interest. When there is measurement error in the covariates, observed
is fitted in the estimating equations instead of
. Generally, there is no closed form
solution to
. One approach to address the problem is to use
the Taylor Series expansion where
. In practice, the expected values in this expression
always need further approximations because they are usually nonlinear functions of or
measurement error u (John, 2010).
15
SIMEX is more widely used due to its simple application. SIMEX is a
simulation-based method to correct bias due to measurement error (Cook and Stefanski,
1994; Stefanski and Cook, 1995; Wang et al., 1998). The general steps for SIMEX are:
First, simulate data based on observed covariate (W) with different settings of
measurement error parameter (
):
, where
is the simulated
covariates based on observed covariates (W) for each simulated sample b =1,…B, and
is the measurement error i.i.d. with mean 0 and covariance
. The covariance of
measurement error
is estimated based on extra data such as replications; Second,
establish a series of naï ve estimates (
) from the simulated data for different settings of
,
; Third, build extrapolation function
=
. The two choices
provided by STATA for extrapolation functions are rational extrapolant:
and quadratic extrapolant:
, where
are
constants to achieve the equations; Fourth, obtain SIMEX estimator (
)
based on the extrapolation function decided in the third step,. When the multiple linear
regression is analyzed with the rational extrapolant, SIMEX estimators are the same as
the moment-based estimator.
Bayesian methods: All the above methods are non-Bayesian approaches.
Gustafson (2004) provided an overview of Bayesian methods in measurement error
problems. Bayesian analysis requires specifying the likelihood of the regression model
and the prior distribution for the unknown parameters. The unknown parameters include
the parameters in the outcome models, measurement error models, and the distribution
16
parameters for the true covariates. The posterior distribution is then expressed as the joint
density of the conditional likelihood of the observed data given all the parameters and the
prior information of those parameters.
1.2.3 Measurement error in covariates under the mixed effects model
For longitudinal studies with multiple follow-up visits, the bias in the naï ve
estimate and its corresponding correction methods for both linear (Paterno and Amemiya,
1996; Tosteson, Buonaccorsi and Demidenko, 1998; Wang et al., 1998) and non-linear
mixed effects models (Wang and Davidian, 1996; Wang et al., 1998; Wang, Lin and
Guttierrez, 1999) have been studied. Many of the tools introduced in previous linear
regression frameworks can also be applied to the mixed effects model. However, the
additional complexity in the mixed effects model creates a computational challenge to the
traditional methods. For example, the induced model approach is restricted to apply
within the conditions that the mean structure keeps the same between the model using
observed covariates and the model using true covariates. However, in the mixed effects
model, the mean structure of fixed effects is only preserved when the product of
covariance matrix of true covariates and the inverse covariance matrix of observed
covariates is proportional to the identity matrix (
). Although the induced
model is not always applicable to mixed effect models, regression calibration can still be
applied by fitting the covariates using true estimated covariates given the observed
covariates. Likelihood-based methods and estimating equations can be very non-trivial to
compute for mixed effects models. SIMEX is relatively easier to apply with an additive
error and an estimated measurement error variance. Different from linear regressions,
17
some measurement error problems within the mixed effects model can be corrected
without additional data to estimate the measurement error variance (Tosteson,
Buonaccorsi et al. 1998; Buonaccorsi, Demidenko et al. 2000; Liang, Wu et al. 2003; Li,
Shao et al. 2005). For example, a mixed effects model for the unobserved time-varying
true covariate can be constructed:
, where
is the mean of X;
is a
vector for random effects with mean zero and covariance matrix
; and is the
design matrix of rank q. This mixed effects model and additional measurement error
model can jointly demonstrate the observed time-varying covariates (
). Hence, likelihood can be constructed based on the original outcome model inserted
with the model for observed time-varying covariates:
, where Z is the design matrix for random effects in the outcome
model;
is the vector of random effects in outcome model; and R is assumed to be
equal to Z for simplification.
1.2.4 Measurement error in the outcome under linear regression
In contrast to the error-prone covariates, few studies have been performed on
evaluating or correcting the bias in the naï ve estimate of the association between the
outcome and covariates due to the measurement error in the outcome. Under the linear
regression framework, Fuller assessed the bias of least square estimates based on the
variance and covariance between measurement error of the outcome and measurement
error of the covariate (Fuller, 1987; John, 2010). It has been shown that when the
measurement error of the outcome is uncorrelated with the measurement error of the
covariate, the bias in the association estimate of the outcome with the covariate does not
18
depend on the measurement error in the outcome. Wang et al. proposed a general solution
to several problems including measurement error, misclassification, and missing data
altogether based on the expected estimating equations (EEE) (Wang et al., 2008). In this
paper, additive measurement error model is assumed for both the outcome and the
covariate. Estimating equations using likelihood scores were constructed for the
parameter (β) of the true outcome (
) given the true covariate (
, the parameter (κ) of
the observed covariate (
) given the true covariate (
), and the parameter (η) of the
observed outcome (
) given the true outcome (
). Finally, three estimating equations
can be constructed:
However, solving these equations demands some
complicated integrations and approximations.
1.2.5 Regression to the mean and mathematical coupling
Much attention has been paid on the assessment of the baseline prediction of the
longitudinal trajectory or the adjustment for the baseline value as a covariate when the
outcome is prone to measurement errors. There are two main methodological concerns of
investigating the baseline prediction of the longitudinal trajectory: regression to the mean
(RTM) (Altman, 1982, 1991a; Altman, 1991b; Barnett, van der Pols and Dobson, 2005;
Blomqvist, 1987; Hayes, 1988; Kirkwood and JAC., 2003; Oldham, 1962) and
mathematical coupling (Archie, 1981; Moreno et al., 1986; Stratton, Feustel and Newell,
1987; Tu et al., 2004).
RTM is a statistical phenomenon in which extremely large or small
measurements tend to be closer to their mean value on subsequent measurements. This
19
occurs when observed values are subject to random intra-individual variations, generally
treated as random measurement error. It should be noted that RTM does not arise from
systematic measurement error (Barnett et al., 2005). For studies with only one baseline
and one follow-up visit, several methods have been developed to address RTM. Oldman
proposed to plot the change of two repeated measurements over the mean value of these
two measurements (Oldham, 1962). MacGregor et al. criticized the orthogonal indices
derived in Oldman’s plot as lack of biological interpretation (MacGregor, Sagnella and
K.D., 1985), and suggested plotting the final measurement against the baseline
measurement in both linear and logarithmic scales. Blomqvist assessed the bias in the
naï ve estimate of the association between the baseline value and the change, and derived
the bias-correction formula in a linear regression framework (Blomqvist, 1977, 1987;
Blomqvist et al., 1977). Several subsequent studies validated that the Blomqvist formula
is most preferred when the RTM is caused by random measurement error in responses
over time (Hayes, 1988; Tu and Gilthorpe, 2007).
Mathematical coupling is defined as two variables in a regression or correlation
analysis that share a common component, or when one variable is contained by the other
variable. The definition of mathematical coupling also implies that two variables are
correlated with each other. In longitudinal studies, the change from baseline is calculated
as follow-up measurement (
) minus the baseline measurement (
). It is obvious that
the baseline (
) is part of the change (
). Due to the internal correlations between
the change and the baseline value in terms of mathematical expressions, naive statistical
analysis of the correlation between the change and the baseline can be led to erroneous
20
conclusions (Anderson, 1990). Oldham’s method (Oldham, 1962) and the variance ratio
test (Green and van de Vijver, 1993; Guilford and Fruchter, 1973) were proposed to solve
mathematical coupling. Oldham’s method addressed the internal correlation between
and
by creating two orthogonal indices (
) and ½(
). Instead of analyzing
the correlations between (
) and
, Oldman suggested to analyze the correlation
between (
) and ½(
), which is interpreted as the correlation between the
change and the half-way value of the responses. The variance ratio test is similar to
Oldham’s method and provides a t statistic to test the correlation between (
) and
(
) with n-2 degrees of freedom expressed as:
, where
is the standard deviation of the baseline
;
is the
standard deviation of the follow-up
;
is the correlation between
and
; and n is
the sample size.
1.2.6 Measurement error problem for the baseline prediction of the longitudinal
trajectory using mixed effects models
For longitudinal studies with multiple follow-up visits, few studies have
investigated the bias in naï ve estimates of the baseline prediction of its subsequent
changes using mixed effects models for continuous responses. One study of Huntington’s
disease assessed the association between the baseline functional score and the
longitudinal trajectory of functional scores over around five years (Byth and Cox, 2005).
The method used to analyze the baseline prediction of the change in functional scores is
illustrated in 1.1.3 with more details included in Chapter 2. Briefly, the mixed effects
21
model with a random intercept and a random time effect was used to analyze the data.
This approach uses the model residual to estimate the measurement error in the response
variable. However, it cannot solve the measurement error problem completely because
measurement error is only a proportion of the model residual. Other factors such as
model misspecification and true biological variations can also influence the model
residual.
Another approach using mixed effects models with the interaction term between
the baseline value of the outcome and the follow-up time, described in 1.1.3, produces the
correct estimate of the baseline prediction of the longitudinal trajectory only when
responses are measured precisely. When responses are subject to measurement error,
naive estimates can be biased using this approach. In the study of eye functions (Gardiner
et al., 2012), authors discussed the results that the observed relationship between the
baseline MD and follow-up changes in MD could be attenuated due to the measurement
error in MD. However, no bias-correction method is provided in their study. In our
previous study of the predictive impact of the baseline DI on the longitudinal trajectory of
DI over time (Xiang et al., 2010), an empirical formula was derived from a simulation
study to correct the bias in the naï ve association estimate. In Chapter 2, we will present
simulation results to show that the naï ve estimate from the mixed effects model is biased
when successive responses are subject to measurement errors. More importantly, we will
provide the new method to correct the bias in the naï ve association estimate. The
comparison between our methods with the Byth’s method will also be performed in the
simulation study.
22
1.2.7 Measurement error problem for the baseline-adjusted association between the
exposure variable and the longitudinal trajectory of responses over time
The modeling strategies of analyzing the association between the exposure
variable and the longitudinal trajectory adjusting for the baseline value were reviewed in
1.1.4. One of the analysis problems occurs when the response variable is subject to
measurement errors. The measurement error in the response variable can cause a biased
estimate of the association between the exposure variable and the longitudinal trajectory
after adjusting for the baseline, if the exposure variable is cross-sectionally correlated
with the baseline value of the outcome, and the baseline value predicts the longitudinal
trajectory over time (Yanez et al., 1998). This is an example of residual confounding
(Fraser and Stram, 2001) where estimates for variables that are not measured with error
can also be biased. In clinical trials, it is suggested that the measurement error in the
response variable can be ignored because the baseline value is expected to be equal
among treatment groups if the trial is designed to be completely randomized and the
expected value of the measurement error is zero (Carroll, Senn and Carroll, 1990; Chan et
al., 2004; Senn, 1990). However, for observational studies, baseline measurement errors
cannot be overlooked because the study design does not have the advantage of
randomization (Chambless and Davis, 2003; Yanez et al., 1998). For studies with only
two time points, Yanez et al. compares the bias of ordinary least square (OLS) estimates
of the association between various risk factors of coronary heart disease and the change
in wall thickness of the carotid artery with or without adjusting for the baseline artery
thickness (Yanez et al., 1998). The results show that the measurement error in the
23
observed artery thickness caused the biased estimate of the association between other risk
factors and the trajectory of the artery thickness when the model was adjusted for the
baseline artery thickness. Based on this observation, Yanez developed the bias-corrected
estimate using the least square estimate of the multiple linear regression. For studies with
multiple visits, Harrison extends the Byth’s approach to investigate the difference of the
longitudinal trajectory of CD4 cell counts after highly active antiretroviral therapy
between two subtypes of HIV (Harrison et al., 2009). Although the model residual is used
to represent the measurement error in responses, Harrison’s approach has the same
limitation as the Byth’s approach that it does not specify the residual variability due to
measurement error or other variations. When the responses are subject to measurement
errors and the exposure variable is correlated with the baseline value, another approach
using two interaction terms (exposure × time + baseline × time) as fixed effects in one
mixed effects model (introduced in 1.1.4) also produces a biased estimate of the baseline-
adjusted association between the exposure variable and the longitudinal trajectory of
responses. Thus, we developed bias correction formulae for the naive estimate using the
mixed effects model with two aforementioned interaction terms in Chapter 3. A
simulation study has been conducted to validate our new methods and to compare our
methods to the Harrison’s method.
24
Chapter 2
Assessing the Baseline Prediction of the Longitudinal
Trajectory when Responses Are Measured with Error
2.1 Literature review
2.1.1 Linear regressions for studies with only one follow-up visit when responses are
measured with error
Blomqvist derived a bias correction formula under linear regression for the
simple case when studies only contain one baseline and one follow-up visit (Blomqvist,
1977). When there are only two time visits in the study, linear regression using the
change model can be fitted to analyze the association between the baseline and the
longitudinal change:
, where
; is the intercept;
;
represents the baseline
measurement;
represents the follow-up measurement; and
captures the association
between the baseline
and the longitudinal change
. When there is measurement
error in the responses and this measurement error follows the normal additive
measurement error model, observed responses
; ; and
. Therefore, we can set the change model using observed data as:
25
, where
is the intercept; and
. Blomqvist’s correction
formula of the naï ve estimate of parameter
was constructed as:
, where equals to one minus the reliability ratio, i.e
. Because measurement
error
is assumed to be independent of the true value
and
is independent and
identically distributed with distribution N(0,
). The additive measurement error model
can be expressed as
,
. The variance of measurement
error
can be estimated using the replication data of individuals from the subsample or
the external sample where the characteristics of subjects are similar to the current study:
(2.1)
, where
is the number of replicates for ith individual; is the number of individuals
with replication data; and
is the sample variance of
replicated measurements of each individual. Although Blomqvist does not discuss the
case when the follow-up time interval is not fixed but varying within a range, we propose
that Blomqvist’s formula is also applicable to the data with a varied time interval.
However, the slope of
is more appropriately used as the outcome of linear
regression compared to the absolute change
. We will compare the results of using
and
as the outcome in the simulation study.
26
2.1.2 Mixed effects models analyzing the change association with exposures for
studies with multiple follow-up visits
For longitudinal data with multiple follow-up visits, a mixed effects model is
generally used to describe the growth trajectory of successive responses considering the
correlation among successive measurements. The between-individual variations of
trajectories can be modeled by the random time effect in the model. The correlations
among measurements taken from the same individual over time can be modeled by
various correlation structures such as exchangeable, autoregressive, and unstructured.
For the ith individual ( ), let
be the measurement taken from the
ith individual at time
. When ,
, and
represents the baseline
measurement for the ith individual. The simplest mixed effects model describing the
linear trend of measurements over time can be modeled as:
(2.2)
, where
is the random intercept and is assumed to be normally distributed with mean
and variance
.
is the fixed effect of follow-up time
.
represents the random
variations among individual trajectories and is assumed to be normally distributed with
mean 0 and variance
; the covariance between
and
is denoted as
; and
is
the residual error and is assumed to be independently and normally distributed with mean
0 and variance
. Based on formula 2.2, we can see that the intercept
is the expected
baseline value
:
(2.3)
27
In real studies, researchers are interested in the exposure which is associated with
the longitudinal trajectory over time. A mixed effects model with the interaction term
between the exposure and time as fixed effect is widely applied to estimate the
association between the exposure and the longitudinal trajectory. Let
be a
vector of the exposure of interest. A mixed effects model used to evaluate the association
can be formulated as:
(2.4)
, where
is the time-invariant or time-variant exposure variable.
is a vector
of parameters representing the cross-sectional correlations between
and
.
is a
vector of parameters describing the magnitude of deviation in the slope of due to
the factor .
represents the magnitude of slope deviation conditioned on . Parameter
is the slope variation across individuals and is assumed to be normally distributed
with the mean 0 and variance
. From another aspect, the mixed model in formula 2.4
can also be expressed using a two-level model:
Level 1:
(2.5)
Level 2:
, where
is the individual slope of follow-up measures over time.
For analysis of the baseline prediction of the longitudinal trajectory, the mixed
effects model can be expressed as:
(2.6)
28
, where
represents the baseline value.
is the parameter of the cross-sectional
association between repeated measurements
and the baseline value
.
is the
parameter representing the magnitude of deviation in the slope of due to the difference
in baseline value
. Because the baseline value is already a covariate in the model 2.6,
the intercept
is not identifiable and can be excluded from the model. Hence, model
2.6 can be simplified as:
(2.7)
,where
represents the time for baseline visit. Model 2.7 can be fitted in a SAS
Proc Mixed procedure (SAS Institute Inc., Cary, NC, USA) or R nlme package (Version
3.1-101) as a no intercept model.
2.1.3 Byth’s random intercept and random slope model for baseline prediction of
the longitudinal trajectory when responses are measured with error
Byth and Cox developed the estimate of baseline prediction of the longitudinal
trajectory using a random intercept and random slope model (Byth and Cox, 2005). Let
the model for true measurements
to be:
, where for each individual. is the number of visits. The
and
. As shown
in formula (2.3), the intercept estimate
is the expected baseline value
. Based on
29
the OLS estimate, the association between baseline
and the longitudinal trajectory can
be estimated by
. When responses are measured with error,
, where
are the observed measurements corresponding to the underlying true value
,
where ;
;
;
and
.
Byth noted that if the residual
is equal to the measurement error
,
is the unbiased estimate of the true association parameter
. In reality,
only
presents a proportion of
. Other variability such as model misspecification or true
biological variation can also contribute to the residual
. This causes the limitation of
its application in real studies. In the simulation study, we will show the situations when
the Byth’s approach fails to produce the correct estimate of the baseline association with
the longitudinal trajectory.
2.2 Commonly used mixed effects models for the baseline association
with the longitudinal trajectory
For statistical analysis, there are various ways to set mixed effect models for the
association between the baseline and the longitudinal trajectory. We considered four
models that are commonly used by analysts and they are listed as follows:
30
Model I(A):
, where ,
,
.
Model II(A):
, where ,
,
.
Model III(A):
, where
, ,
,
.
Model IV(A):
, where
, ,
,
.
The difference between Model I(A) and Model II(A) is that Model I(A) excludes
the baseline measurement
from the dependent variable, but Model II(A) retains the
baseline measurement
in the dependent variable. Model I(A) is described in Garrett’s
book and is the strategy to analyze the association between the baseline value and the
longitudinal trajectory of responses (Garrett et al., 2004). Model II(A) was used in some
studies (Gardiner et al., 2012; Josephs et al., 2011) for association analysis, but it is an
incorrect model because the baseline value is included in both the dependent and
independent variable. Comparing Model III(A) with Model I(A), the difference is that the
dependent variable of Model III(A) is the change of post-baseline measurements from the
baseline measurement. Using absolute changes instead of successive measurements in the
dependent variable can be viewed as a shift of all the measurements downward by the
magnitude of baseline value. Therefore, Model III(A) yields exactly the same estimates
of
,
and
as Model I(A). The only difference is
from Model III(A) is equal
31
to the
from Model I(A) minus one. Model IV(A) does the same shift as Model II(A),
thereby yields identical estimates of
,
and
as Model II(A).
from model
IV(A) is equivalent to the
from Model II(A) minus one. Although Model III(A) and
Model IV(A) are treated as the same for the estimate of
as Model I(A) and Model
II(A), Garrett et al. (2004) suggests to use successive measurements as the dependent
variable because using the absolute changes as the dependent variable will make the
interpretation of coefficient estimates be more complicated. Hence, we will use Model
I(A) and Model II(A) as the main analysis model in the following work.
When there is measurement error in the response variable
, the parameter
estimate of the baseline association with the longitudinal trajectory using observed
measurements tends to be biased. First, we apply the additive measurement error model
to depict the relationship between the true responses and the observed responses
measured with error. Let
be the true responses and
be the corresponding observed
responses with measurement error. Let
be the measurement error for jth measurement
of ith individual.
is assumed to be identically, independently, and normally distributed
with mean 0 and variance
. The additive measurement error model can be expressed as:
, where
is independent of
. However, only observed data were
available in real studies, thus, Model I(A) and Model II(A) for data analysis can be
expressed as follows where we use to represent the naï ve coefficient estimates from
models using observed responses
, to distinguish from,
which are the true
parameters of models using the true responses
:
32
Model I(a):
, where ,
,
.
Model II(a):
, where ,
,
.
2.3 Theoretical derivation of the bias correction formula for the baseline
association with the longitudinal change using mixed effects models
2.3.1 Derivation of the bias correction formula for Model I(a) when the baseline is
excluded from the dependent variable
As shown in formula 2.5, a mixed effects model can be formulated as a two-level model.
Thus, Model I(A):
can be expressed as
a two-level model:
Level 1:
(2.8)
Level 2:
.
Model I(a) for observed responses:
can be expressed as a two-level model:
Level 1:
(2.9)
Level 2:
.
33
Based on level 1 of formula 2.8, let
and
. The least square
estimate
, where
.
Therefore,
is equal to
.
Similarly, for Model I(a) using observed responses,
(2.10)
and
.
Based on level 2 of formula 2.8, least square estimate of
can be obtained as:
(2.11).
Because
and
,
from result 2.11 can be manipulated as:
.
Here, we employ the same definition as the Blomqvist’s method that
can be estimated using the replication data as formula 2.1;
can be estimated
from the sample variance of observed baseline measurements:
; and
can be interpreted as the proportion of measurement error
34
variance in the entire variance of the baseline value of the response variable. Finally, the
bias-corrected estimate of
can be simplified as:
(2.12)
Result 2.12 shows that the bias-corrected estimate
depends on the naï ve estimate
and measurement error parameter .
Although we developed the correction formula 2.12 based on the mixed effects
model with a linear time variable, we can apply the same method to any polynomial time
effects. The only change is that we need to use the polynomial time effect instead of
linear
in the correction formula. For example, if the repeated measurements follow a
quadratic trend over time instead of a linear trend, the mixed effects model for the true
responses can be set as:
.
The analysis model for the observed responses is:
.
Finally, following the same derivation process as we did for linear time variable, the
correction formula
can be expressed as:
(2.13)
2.3.2 Derivation of the bias correction formula for Model II(a) when the baseline is
included in the dependent variable
Although Model II(a) is an incorrect analysis model, we derive the bias correction
formula for this model if researchers still use it for data analysis. The derivation of the
35
bias-corrected estimate for Model II (a) is similar to Model I(a). First, model II(A) and
model II(a) can be transformed into two level models using formula 2.8 and 2.9,
respectively. However, since the baseline value is included in the dependent variable,
and
. Based on the level one of formula 2.8, the least square
estimate of
(2.14)
, where
.
Because
, result 2.14 can be extended as :
. (2.15)
After this mathematical manipulation,
is separated into two components. The first two
summations in the parenthesis of 2.15 compose the component including all post-baseline
measurements and the third element in the parenthesis including the component of
baseline measurement. The correlation between
and the baseline
can be expressed
as:
(2.16)
, where
, representing the sample estimate of
.
Similarly, for Model II(a) using observed data,
and
.
36
The least square estimate of
(2.17).
Hence,
(2.18).
Based on level two of formula 2.8, the least square estimate of
can be expressed as:
(2.19).
Based on the assumption of additive measurement error model,
=0.
Additionally,
(2.20). After inserting the result 2.18 and 2.20 into 2.19,
(2.21).
Based on level two of formula 2.9, the least square estimate of
. Thus,
result 2.21 can be reformulated as
(2.22).
37
Let the estimate of the measurement error variance be
, the variance of the observed
baseline response be
, and
. The final bias correction formula for
can
be simplified as:
(2.23).
It can been seen that if the study is designed as fixed-length of follow-up intervals and
there is only one baseline and one follow-up visit, . In this case, formula 2.23 will
be exactly the same as the Blomqvist’s formula. Additionally, result 2.23 shows that the
estimate of the true association parameter
depends on the naï ve estimate of
, the
number of follow-up visits, and
.
When a polynomial time effect is used instead of a linear time variable, the bias-
corrected estimate
can be expressed as
(2.22)
,where
represents any mathematical transformation of time variable, such as linear,
quadratic, logarithm, etc.
and
=
.
38
2.4 Simulation study
2.4.1 Simulation design
For simulation studies, we generated sample size of 50, 100, and 500 per cohort to
mimic the sample size of usual clinical studies. For each set of underlying parameters, we
repeated simulation 1000 times and the results were summarized by the average of 1000
replications. The bias of the parameter estimate was calculated by subtracting the true
parameter from the parameter estimate obtained from the simulated data. Percent bias
was calculated by dividing the absolute bias by the true parameter. The true baseline
measurement
was generated following a standard normal distribution . In real
studies, standard normal distribution can be achieved by standardizing the baseline
measurement using the baseline mean and standard deviation. Baseline measurement
error
is assumed to be identically, independently, and normally distributed with mean
0 and variance
. The simulated observed baseline measurement was equal to the true
baseline measurement plus the measurement error, i.e.,
. The variance of
the observed baseline measurement was equal to (1+
). To be consistent with the
Blomqvist’s approach, was generated as the proportion of the total variance of the
observed baseline measurement due to measurement error, i.e.
.
In this simulation study, we set to be varied from 0 to 0.5 (0, 0.01, 0.05, 0.1, 0.2, 0.3,
0.4, and 0.5). Since we assumed
follows the standardized normal distribution with a
variance equal to one, the corresponding
was set to be varied from 0 to 1 (0, 0.1005,
0.2294, 0.3333, 0.5, 0.6547, 0.8165, and 1).
39
The growth trajectory of repeated measurements was assumed to be linearly
associated with the baseline value in a regression model as:
.
was the population-averaged time effect, which was set to be 0, -0.5, and 0.5.
Considering the variance of the baseline
was one,
was set to be varied among (-
0.50, -0.05, 0.0, 0.05, and 0.5) representing reasonable moderate to small associations in
various directions between the true baseline value and the longitudinal trajectory.
was
assumed to be normally distributed with N(0,
).
was set to be 0.2 and 0.4 in the
simulation study. We generated two settings for follow-up time intervals. First, each visit
was assumed to have an equal time interval corresponding to each annual visit. Second,
the time interval varied from 8 to 16 months with a scheduled visit at one year. Uniform
distribution was used to generate the varied time intervals by .
Datasets were generated with up to ten follow-up visits (j=1, 2, 3, …, and 10). Based on
the true baseline measurements and the individual trajectories we generated above, the
true follow-up measurements for the th individual at th time were simulated using the
following model:
, where
is randomly distributed with N(0,
).
was set to be 0.01 to 4. The
corresponding observed measurements at each follow-up visit was equal to the true
measurement plus measurement error
:
, where
is identically,
independently, and normally distributed with mean 0 and variance
.
In clinical studies, missing data always occurs during the follow-up. Hence, we
generated 20% to 50% data missing at random in the follow-up measurements following
40
a uniform distribution . Only the situation of completely random missing data
was studied. Situations other than completely random missing data were not investigated
in this study because non-random missing data can cause biased estimates itself.
We verified that the coefficient estimate
does not depend on
and
based
on our simulation results. Therefore, we only present the result using
and
to be 0
and 0.2. In the association analysis, Model I(a), Model II(a), and the Byth’ approach were
applied sequentially using the same dataset. Thus, fair comparison can be performed
among three methods. An unstructured within-individual correlation matrix was used in
the mixed effects models. The bias-correction formulae were applied to the naï ve
association estimates obtained from Model I(a) and Model II(a). The sample estimate of
the variance of measurement error
was calculated as:
(2.26).
Finally, the naive estimate
and the bias-corrected estimate
were
compared among Model I(a), Model II(a) and Byth’s approach. For the case that there
were only one baseline and one follow-up visit (j=1), bias-corrected estimates
using
Model I(a) and Model II(a) were also compared to the linear regression estimate
corrected by the Blomqvist’s formula. In this comparison, we used a fixed time effect for
Model I(a) and Model II(a) instead of a random time effect because the parameter
estimates of random effect did not converge. Moreover, when there was only one follow-
up visit and the follow-up interval was fixed among individuals, Model I(a) still has the
convergence problem even if the time effect was set to be fixed. Therefore, results were
41
not compared between the bias-corrected estimate of Model I(a) and the Blomqvist’s
approach when the simulated data had only two time-points and the time intervals were
fixed.
2.4.2 Simulation results
Table 2.1 reports the simulation results with two, four, and six follow-up visits
with fixed time interval of Model I(a) and Model II(a) with the sample size of 100. Model
I(a) shows the following: 1) when =0, i.e., there is no measurement error, the regression
coefficient estimating the baseline association with the longitudinal trajectory using the
observed data is essentially equal to the true association, thus no bias was observed, as
expected; 2) when there is no baseline association with the trajectory, the coefficient
estimate
is not biased even there is measurement error, i.e. ; 3) when there is
truly an association between the baseline value and the longitudinal trajectory, the bias of
compared to the true parameter increases when the measurement error increases; 4)
the magnitude of bias does not depend on the number of follow-up visits. Model II(a)
shows the following: 1) when =0, the estimate of the baseline association with the
longitudinal trajectory using the observed data is not biased; 2) the bias increases with the
increase of measurement error parameter conditioned on the same number of follow-up
visits; 3) bias decreases with the increase of follow-up visits if other parameters are fixed.
The relationships between the biases of naï ve estimates and the degree of random
measurement error, the number of follow-up visits, and the true underlying association
parameter using Model I(a) and Model II(a) are summarized in Figures 2.1-2.3. Figure
2.1 shows the scatter plot of biases (
) against the random measurement error
42
( ) by the number of follow-up visits when the true association (
) is zero (Model I(a):
Fig. 2.1a; Model II(a): Fig. 2.1b), -0.25 (Model I(a): Fig. 2.1c; Model II(a): Fig. 2.1d),
and 0.25 (Model I(a): Fig. 2.1e; Model II(a): Fig. 2.1f). For Model I(a), the bias of
has a linear relationship with measurement error ( ) and slopes of this relationship varies
by the true association parameter (
). However, the slope of bias over measurement
error does not vary by the number of follow-up visits. For Model II(a), the bias in
and measurement error ( ) also follows a linear relationship. Slopes of this linear
relationship vary by the number of follow-up visits and the true association parameter
(
). Larger bias is associated with larger measurement error, and fewer numbers of
follow-ups. Figure 2.2 is the scatter plot of biases (
) against the true underling
association parameter (
) by the number of follow-up visits when the degree of
measurement error ( ) equals to 10% (Model I(a): Fig. 2.2a; Model II(a): Fig. 2.2b), 30%
(Model I(a): Fig. 2.2c; Model II(a): Fig. 2.2d), and 50% (Model I(a): Fig. 2.2e; Model
II(a): Fig. 2.2f). For both Model I(a) and Model II(a), bias (
) is linearly related
with the true parameter (
). For Model I(a), slopes of the linear relationship vary by
the measurement error , but not by the number of follow-up visits. However, for Model
II(a), slopes of the linear relationship are not only varied by measurement error , but
also across the vertical axis with different intercepts for different numbers of follow-up
visits. Because the bias of
is not associated with the number of follow-up visits in
Model I(a), we only present the figures of bias against the inverse number of visits for
Model II(a). Figure 2.3 shows the bias against the inverse number of follow-up visits
43
across a different measurement error ( ) when the true association parameter (
) is 0
(Fig. 2.3a), -0.25 (Fig. 2.3b), and 0.25 (Fig. 2.3c). The biases in naï ve estimates
and
the inverse number of follow-up visits can be described as a linear relationship and slopes
of this relationship vary by the degree of random measurement error. Simulated data with
data missing at random or with varying follow-up time intervals for follow-up visits do
not change the above results.
2.4.3 Comparison of Model I(a), Model II(a) and Byth’s method
Table 2.2 compares the naï ve estimate
using Model I(a), Model II(a) and
Byth’s method for the robustness of the estimate of the baseline association with the
longitudinal change across various settings of residual variance (
). The results
presented are derived from the simulated data with a fixed time interval and sample size
of 100. This shows that the parameter estimates (
) using Model I(a) and Model II(a)
do not vary by the magnitude of residual variance (
). For the Byth’s model, when
is
close to 0,
is very close to the true relationship
. The bias of
using the
Byth’s model increases when the residual variance increases. In our simulation study,
when the residual variance increases to be the same as or larger than the variance of the
baseline value
, which is set to be one in the simulation, the estimate of baseline
variance (
) can be close to zero, and
of the Byth’s model turns out to be
tremendously large. Figure 2.4 draws the picture of the bias in
using Model I(a) (Fig.
2.4a, 2.4b), Model II(a) (Fig. 2.4c, 2.4d) and the Byth’s model (Fig. 2.4e, 2.4f) against
residual variance (
) with a fixed follow-up time interval and the number of visits from
44
two to ten. Both Model I(a) and Model II(a) have consistent estimates of
across the
different residual variance
in each setting of the true association parameter (
) and
the measurement error ( ). When there is no measurement error, i.e. =0,
derived
from Model I(a) and Model II(a) are not biased regardless of the residual variance and the
number of follow-up visits. For the Byth’s method,
is not biased when the residual
variance
is close to zero. However, when even there is no measurement error,
is
still biased when
is not equal to zero. For a small number of follow-up visits and
relatively large
,
turns out to be unexpectedly large as we explained in table 2.2.
Therefore, the Byth’s model is only appropriate when the underlying true trajectory of
successive measurements is perfectly linear over time. Otherwise, even there is no
measurement error, Byth’s method is still not able to provide correct estimates of the
association between the baseline value and the longitudinal trajectory.
2.4.3 Validation of proposed bias correction formulae
To evaluate the performance of our new bias-correction formulae for Model I(a)
(2.12) and Model II(a) (2.23), we simulated various datasets with different combinations
of parameter settings. Table 2.3, Figure 2.5, and Figure 2.6 present some selected
examples of the validation results from the simulated data with both fixed and varied time
intervals. Naively applying mixed effects models can cause very biased results. However,
applying Formulae (2.12) and (2.23) produce bias-corrected estimates (
) very close to
the true association parameter (
) for all cases, even for a cohort size as small as 50.
Results using varied time intervals for follow-up visits are similar to the results using
45
fixed-time intervals (Table 2.3, Figure 2.7, and Figure 2.8). Moreover, we assessed the
bias correction Formulae (2.12) and (2.23) when there were data missing at random.
Table 2.4 shows the simulation results using the same parameter settings as Table 2.3.
The results of naï ve estimates (
) and bias-corrected estimates (
) are very close to
the results in Table 2.3. However, the bias in the bias-corrected estimates (
) using
Model II(a) is greater than the bias using Model I(a) (Figure 2.9). In conclusion, our bias
correction formulae successfully correct the bias in the naï ve estimate of the baseline
association with the longitudinal trajectory for both Model I(a) and Model II(a). When
there are data missing at random, the bias in the bias-corrected estimate using Model I(a)
is even smaller than the bias using Model II(a).
2.4.4 Comparison of proposed bias correction formula with the Blomqvis’s formula
When there is only one baseline and one follow-up visit, both the mixed effects
model and linear regression can be used to assess the association between the baseline
value and the longitudinal trajectory. Table 2.5 compares the bias-corrected estimate
(
) using bias-correction formulae (2.12) and (2.23) under the mixed effects models
framework with the bias-corrected estimates using the Blomqvist’s approach under the
linear regression framework. The simulation results are based on a sample size of 100
with 1000 replications. Both data with fixed time intervals and varied time intervals were
analyzed. When the time interval is varied from 8 to 16 months for one year visit, linear
regressions with two different outcomes are additionally compared in two ways. First, the
outcome is raw change from the baseline to the follow-up; and Second, the outcome is
the rate of change (the raw change divided by the follow-up time). When the follow-up
46
time interval is fixed across subjects, regression coefficients using the raw change and the
rate of change as outcomes are identical. Hence, we only present the results from the
model using the raw change as the outcome. For a fixed time interval, bias-corrected
estimates (
) using correction formula (2.23) for Model II(a) and the Blomqvist’s
formula for linear regression are identical and very close to the true association parameter
for the baseline with the longitudinal trajectory (
) when the outcomes are prone to
measurement error. Model I(a) cannot be fitted for data with fixed time intervals due to
insufficient data variations for parameter estimations. For a varied time interval, the
Blomqvist formula provides similar results as the bias-correction formula (2.12) for
Model I(a). Using formula (2.23) based on Model II(a) presents more biased estimates.
Comparing two linear regression models using different outcomes, corrected estimates
using the Blomqvist’s formula based on model using the rate of change as the outcome is
closer to the true association parameter
than using the raw change as the outcome
when there is measurement error in the outcomes. The fact that bias-corrected estimates
(
) using Model II(a) and formula (2.23) are more biased when there is measurement
error in the outcomes can be because the baseline value is included as both the dependent
and independent variables in the model, thus the estimation of baseline variance is
incorrect using such model, which further influences the correctness of association
estimates. Even under the case of two time points, we still see correction formula (2.12)
for Model I(a) successfully corrects biases in the naï ve estimates when responses are
measured with error.
47
Table 2.1. Simulation Results from Two, Four and Six Follow-up Visits
a
True Parameters
b
Model I(a)
c
(
) Model II(a)
c
(
)
λ n=2 n=4 n=6 n=2 n=4 n=6
0 0 0.000 0.000 0.000 0.000 0.000 0.000
0 0.01 0.001 0.001 0.000 -0.006 -0.003 0.000
0 0.05 -0.001 0.000 0.001 -0.025 -0.011 -0.005
0 0.1 0.003 -0.001 0.000 -0.048 -0.020 -0.010
0 0.2 -0.001 0.001 -0.001 -0.099 -0.041 -0.021
0 0.3 -0.002 -0.001 0.000 -0.150 -0.060 -0.033
0 0.4 -0.001 0.001 0.000 -0.200 -0.080 -0.043
0 0.5 -0.004 0.000 0.000 -0.252 -0.100 -0.054
-0.25 0 -0.250 -0.249 -0.250 -0.250 -0.250 -0.250
-0.25 0.01 -0.249 -0.247 -0.248 -0.252 -0.25 -0.248
-0.25 0.05 -0.239 -0.237 -0.237 -0.263 -0.247 -0.242
-0.25 0.1 -0.226 -0.225 -0.225 -0.274 -0.245 -0.235
-0.25 0.2 -0.202 -0.202 -0.198 -0.300 -0.239 -0.222
-0.25 0.3 -0.177 -0.176 -0.175 -0.324 -0.235 -0.207
-0.25 0.4 -0.148 -0.151 -0.150 -0.350 -0.231 -0.192
-0.25 0.5 -0.123 -0.126 -0.125 -0.379 -0.226 -0.179
-0.5 0 -0.499 -0.500 -0.498 -0.501 -0.500 -0.500
-0.5 0.01 -0.495 -0.495 -0.495 -0.501 -0.496 -0.496
-0.5 0.05 -0.476 -0.476 -0.475 -0.502 -0.484 -0.481
-0.5 0.1 -0.451 -0.448 -0.45 -0.501 -0.47 -0.46
-0.5 0.2 -0.398 -0.399 -0.4 -0.499 -0.441 -0.421
-0.5 0.3 -0.348 -0.35 -0.352 -0.498 -0.411 -0.382
-0.5 0.4 -0.297 -0.3 -0.3 -0.5 -0.38 -0.343
-0.5 0.5 -0.251 -0.252 -0.249 -0.501 -0.35 -0.303
0.25 0 0.249 0.249 0.250 0.249 0.249 0.250
0.25 0.01 0.249 0.247 0.247 0.242 0.246 0.247
0.25 0.05 0.237 0.237 0.237 0.212 0.226 0.232
0.25 0.1 0.226 0.226 0.225 0.174 0.205 0.216
0.25 0.2 0.195 0.202 0.200 0.099 0.159 0.178
0.25 0.3 0.176 0.176 0.174 0.026 0.114 0.141
0.25 0.4 0.147 0.149 0.15 -0.05 0.07 0.107
0.25 0.5 0.132 0.123 0.124 -0.123 0.026 0.071
0.5 0 0.501 0.500 0.500 0.500 0.499 0.501
0.5 0.01 0.495 0.495 0.495 0.491 0.493 0.494
0.5 0.05 0.477 0.475 0.476 0.45 0.465 0.469
0.5 0.1 0.451 0.450 0.451 0.399 0.430 0.438
0.5 0.2 0.398 0.401 0.400 0.300 0.360 0.379
0.5 0.3 0.353 0.351 0.35 0.200 0.290 0.318
0.5 0.4 0.305 0.299 0.301 0.099 0.219 0.258
0.5 0.5 0.247 0.250 0.249 -0.003 0.152 0.199
a. Data shown from sample size of 100 with 1000 replications. Linear regression was used for one follow-up visit,
random effect mixed model was used for two and five follow-up visits.
48
b. Parameters used in generating the data in simulation study: n is the number of follow-up visits;
is the true
change per visit associated with per unit difference of baseline value; is the proportion of total observed variance in
outcome due to random measurement error; The residual variance
was set to be 0.01.
c. Mean of estimated change per visit associated with per unit difference of baseline value (
)) from 1000
replications.
Table 2.2. Comparison of Simulation Results from Model I(a), Model II(a) and Byth’s model
a
True Parameters
b
Model I(a) Estimated
Parameters
c
(
)
Model II(a) Estimated
Parameters
c
(
)
Byth’s model(a) Estimated
Parameters
c
(
)
λ
n=2 n=4 n=6 n=2 n=4 n=6 n=2 n=4 n=6
0 0 0.01 -0.001 0.000 0.000 0.000 0.000 0.000 0.003 0.001 0.000
0 0 0.25 -0.002 0.000 0.000 -0.001 0.000 0.000 0.091 0.027 0.011
0 0 1 -0.003 -0.004 0.000 -0.001 -0.002 0.000 0.522 0.139 0.050
0 0 4 0.003 0.002 -0.001 0.003 0.001 -0.001 2.78E+30 5.8E+28 1.06E+13
0 0.5 0.01 -0.005 -0.001 0.000 -0.251 -0.100 -0.053 0.000 0.009 0.005
0 0.5 0.25 0.001 0.001 0.000 -0.249 -0.100 -0.054 0.074 0.038 0.014
0 0.5 1 0.003 0.003 0.000 -0.248 -0.099 -0.053 0.948 0.167 0.062
0 0.5 4 -0.013 0.001 -0.001 -0.252 -0.100 -0.054 3.38E+15 1.7E+29 1.63E+13
-0.25 0.5 0.01 -0.127 -0.124 -0.125 -0.376 -0.224 -0.178 -0.245 -0.247 -0.250
-0.25 0.5 0.25 -0.121 -0.125 -0.124 -0.374 -0.225 -0.178 -0.198 -0.233 -0.251
-0.25 0.5 1 -0.120 -0.127 -0.125 -0.376 -0.226 -0.179 1.73E+15 -0.200 -0.254
-0.25 0.5 4
-0.112 -0.127 -0.123 -0.370 -0.226 -0.178 4.99E+45 -2.1E+62 -7.4E+47
0.5 0.5 0.01 0.248 0.250 0.250 -0.002 0.150 0.196 0.558 0.524 0.512
0.5 0.5 0.25 0.242 0.249 0.250 -0.002 0.148 0.196 0.808 0.597 0.549
0.5 0.5 1 0.241 0.250 0.250 -0.002 0.150 0.198 3.45E+13 1.007 0.712
0.5 0.5 4 0.238 0.250 0.250 -0.004 0.150 0.196 1.33E+16 3.69E+14 5.91E+13
a. Data shown from sample size of 100 with 1000 replications.
b. Parameters used in generating the data in simulation study: n is the number of follow-up visits;
is the true change per visit associated with per unit difference of
baseline value; is the proportion of total observed variance in outcome due to random measurement error.
is the residual variance of the mixed effects models;
c. Mean of estimated change per visit associated with per unit difference of baseline value with 1000 replications. Uncorrected
is the model based estimate from
random effect mixed model.
49
Table 2.3. Validation Results from Simulated Data
Simulation Parameters
a
Estimated Parameters
b
Model I(a) Model II(a)
Cohort Size n
λ* Uncorrected
Bias Corrected
Uncorrected
Bias Corrected
Fixed Follow-up Time Interval
50 3 -0.5 0.3 0.25
0.225 0.301 0.151 0.315
100 3 -0.5 0.3 0.25 0.224 0.298 0.150 0.305
500 3 -0.5 0.3 0.25
0.224 0.299 0.150 0.300
50 5 0.5 -0.4 0.2 -0.320 -0.400 -0.349 -0.405
100 5 0.5 -0.4 0.2
-0.321 -0.401 -0.348 -0.404
500 5 0.5 -0.4 0.2
-0.320 -0.400 -0.349 -0.401
Varied Follow-up Time Interval
50 3 -0.5 0.3 0.25
0.223 0.297 0.121 0.319
100 3 -0.5 0.3 0.25
0.221 0.295 0.121 0.310
500 3 -0.5 0.3 0.25
0.226 0.301 0.122 0.304
50 5 0.5 -0.4 0.2
-0.320 -0.400 -0.363 -0.403
100 5 0.5 -0.4 0.2 -0.319 -0.398 -0.362 -0.400
500 5 0.5 -0.4 0.2
-0.321 -0.401 -0.363 -0.400
a. Parameters used in generating the data in simulation study: n is the number of follow-up visits;
is the true change per visit not associated with baseline value;
is
the true change per visit associated with per unit difference of baseline value; is the proportion of total observed variance in outcome due to random measurement
error; The residual variance
was set to be 0.01.
b. Mean of estimated change per visit associated with per unit difference of baseline value with 1000 replications. Uncorrected
is the model based estimate from
random effect mixed model; bias corrected
is calculated using formula (2.12) for model I(a) and formula (2.23) for model II(a).
c. Fixed follow-up time interval assumed each visit having equal time interval and corresponding to annual visit, varied follow-up time interval varied from 8 to 16
months with a scheduled visit at one year
50
Table 2.4. Validation Results from Simulated Data with 20% data missing at random
Simulation Parameters
a
Estimated Parameters
b
Model I(a) Model II(a)
Cohort Size n
λ* Uncorrected
Bias Corrected
Uncorrected
Bias Corrected
Fixed Follow-up Time Interval
c
50 3 -0.5 0.3 0.25 0.224 0.298 0.142 0.306
100 3 -0.5 0.3 0.25 0.225 0.301 0.145 0.299
500
50
3
5
-0.5
0.5
0.3
-0.4
0.25
0.2
0.226 0.301 0.145 0.294
-0.321 -0.402 -0.353 -0.411
100 5 0.5 -0.4 0.2 -0.321 -0.402 -0.352 -0.407
500 5 0.5 -0.4 0.2 -0.320 -0.400 -0.352 -0.404
Varied Follow-up Time Interval
c
50 3 -0.5 0.3 0.25 0.228 0.303 0.120 0.305
100 3 -0.5 0.3 0.25 0.227 0.303 0.119 0.296
500
50
3
5
-0.5
0.5
0.3
-0.4
0.25
0.2
0.225 0.301 0.119 0.291
-0.323 -0.404 -0.369 -0.416
100 5 0.5 -0.4 0.2 -0.320 -0.401 -0.367 -0.411
500 5 0.5 -0.4 0.2 -0.320 -0.400 -0.366 -0.410
a. Parameters used in generating the data in simulation study: n is the number of follow-up visits;
is the true change per visit not associated with baseline value;
is
the true change per visit associated with per unit difference of baseline value; is the proportion of total observed variance in outcome due to random measurement
error; The residual variance
was set to be 0.01.
b. Mean of estimated change per visit associated with per unit difference of baseline value with 1000 replications. Uncorrected
is the model based estimate from
random effect mixed model; bias corrected
is calculated using formula (2.12) for model I(a) and formula (2.23) for model II(a).
c. Fixed follow-up time interval assumed each visit having equal time interval and corresponding to annual visit, varied follow-up time interval varied from 8 to 16
months with a scheduled visit at one year.
51
52
Table 2.5. Comparison of bias-corrected estimates
of mixed effects model and
linear regression
True Parameters
b
Bias-corrected Estimates
c
(
) – True parameter (
)
λ
Fixed Time Interval Varied Time Interval
ModelI(a) Blomqvist
1
ModelI(a) ModelII(a) Blomqvist
1
Blomqvist
2
0 0 -0.0010 -0.0010 0.0020 0.0010 0.0010 0.0000
0 0.01 0.0010 0.0010 -0.0070 0.0010 0.0010 0.0010
0 0.05 0.0020 0.0020 -0.0050 0.0060 0.0000 -0.0020
0 0.1 0.0060 0.0060 0.0020 0.0170 0.0050 0.0010
0 0.2 0.0080 0.0080 -0.0010 0.0420 0.0150 0.0050
0 0.3 0.0090 0.0090 -0.0080 0.0610 0.0160 -0.0010
0 0.4 0.0310 0.0310 0.0050 0.1080 0.0370 0.0100
0 0.5 0.0400 0.0400 -0.0070 0.1530 0.0410 -0.0010
-0.25 0 0.0000 0.0000 -0.0030 0.0000 -0.0010 0.0000
-0.25 0.01 0.0000 0.0000 0.0020 0.0010 0.0000 -0.0010
-0.25 0.05 0.0010 0.0010 0.0000 0.0060 0.0000 -0.0020
-0.25 0.1 0.0030 0.0030 0.0090 0.0140 0.0010 -0.0030
-0.25 0.2 0.0040 0.0040 0.0280 0.0380 0.0090 -0.0010
-0.25 0.3 0.0140 0.0140 -0.0020 0.0590 0.0120 -0.0050
-0.25 0.4 0.0160 0.0160 0.0030 0.0890 0.0170 -0.0110
-0.25 0.5 0.0150 0.0150 0.0200 0.1470 0.0310 -0.0040
-0.5 0 0.0010 0.0010 0.0020 0.0000 -0.0010 -0.0010
-0.5 0.01 0.0000 0.0000 0.0020 0.0010 0.0000 -0.0010
-0.5 0.05 0.0010 0.0010 -0.0060 0.0070 0.0020 0.0000
-0.5 0.1 0.0010 0.0010 -0.0110 0.0110 -0.0020 -0.0050
-0.5 0.2 0.0010 0.0010 0.0150 0.0340 0.0050 -0.0040
-0.5 0.3 0.0070 0.0070 -0.0130 0.0520 0.0090 -0.0100
-0.5 0.4 0.0120 0.0120 -0.0070 0.0790 0.0060 -0.0190
-0.5 0.5 0.0230 0.0230 0.0490 0.1370 0.0240 -0.0210
0.25 0 0.0010 0.0010 0.0000 0.0000 0.0000 0.0000
0.25 0.01 0.0010 0.0010 -0.0010 0.0010 0.0000 0.0000
0.25 0.05 0.0050 0.0050 -0.0110 0.0070 0.0020 0.0000
0.25 0.1 0.0050 0.0050 -0.0190 0.0140 0.0040 0.0000
0.25 0.2 0.0080 0.0080 0.0030 0.0400 0.0130 0.0030
0.25 0.3 0.0190 0.0190 -0.0080 0.0630 0.0170 -0.0010
0.25 0.4 0.0290 0.0290 0.0420 0.1050 0.0280 -0.0010
0.25 0.5 0.0660 0.0660 0.0070 0.1700 0.0550 0.0130
0.5 0 -0.0010 -0.0010 0.0050 0.0000 0.0000 0.0000
0.5 0.01 0.0010 0.0010 -0.0010 0.0010 0.0000 -0.0010
0.5 0.05 0.0020 0.0020 0.0010 0.0080 0.0020 0.0000
0.5 0.1 0.0020 0.0020 -0.0180 0.0160 0.0050 0.0010
0.5 0.2 0.0110 0.0110 0.0160 0.0270 -0.0020 -0.0130
0.5 0.3 0.0280 0.0280 -0.0100 0.0630 0.0180 0.0000
0.5 0.4 0.0380 0.0380 -0.0340 0.1010 0.0310 0.0040
0.5 0.5 0.0560 0.0560 0.0080 0.2090 0.0930 0.0500
53
a. Data shown from sample size of 100 with 1000 replications. Linear regression was used for one follow-up visit,
random effect mixed model was used for two and five follow-up visits.
b. Parameters used in generating the data in simulation study: n is the number of follow-up visits;
is the true
change per visit associated with per unit difference of baseline value; is the proportion of total observed
variance in outcome due to random measurement error; The residual variance
was set to be 0.01.
c. Mean of bias-corrected estimated change per visit associated with per unit difference of baseline value (
) from
1000 replications. Correction formula (2.20) was used for model I(a). Correction formula (2.24) was used for
model II(a). Blomqvist’s formula was used for linear regression of change model. Blomqvist
1
used raw change
as the outcome. Blomqvist
2
used slope of
as the outcome.
54
Figure 2.1. Bias of naï ve estimates (
) using Model I(a) and Model II(a)
against the proportion of variance in the total variance of observed baseline value due to
measurement error. The scatter plot are presented across different settings of true
parameter (
) of the association between baseline and longitudinal trend by number of
follow-up visits.
b) Model I(a):
d) Model I(a):
f) Model I(a):
a) Model II(a):
c) Model II(a):
e) Model II(a):
Bias (estimate – true) Bias (estimate – true) Bias (estimate – true)
Measurement error (λ) Measurement error (λ)
55
Figure 2.2. Bias of naï ve estimates (
) using Model I(a) and Model II(a)
against the true underlying parameter of the association between baseline value and
longitudinal slop (
). The scatter plots are presented across different settings of the
ratio of variance of measurement error to the total variance of observed baseline value ( )
by number of follow-up visits.
b) Model I(a): error = 10%
d) Model I(a): error = 30%
f) Model I(a): error = 50%
a) Model II(a): error = 10%
c) Model II(a): error = 30%
e) Model II(a): error = 50%
Bias (estimate – true) Bias (estimate – true) Bias (estimate – true)
Regression Parameter Regression Parameter
56
Figure 2.3. Bias of naï ve estimates (
) using Model I(a) against the inverse of
number of follow-up visits. The scatter plots are presented across different settings of
settings of true parameter (
) of the association between baseline and longitudinal
trend by number of follow-up visits.
Bias (estimate – true) Bias (estimate – true) Bias (estimate – true)
1 / Number of follow-up visits (n)
a) Model II(a):
b) Model II(a):
c) Model II(a):
57
Figure 2.4. Bias of naï ve estimates (
) using Model I(a), Model II(a) and
Byth’s model against the true variance of model residual (
). The scatter plots are
presented across different settings of true parameter (
) of the association between
baseline and longitudinal trend and the proportion of variance in the observed baseline
value due to measurement error (λ) by number of follow-up visits.
Bias (estimate – true) Bias (estimate – true) Bias (estimate – true)
a) Model I(a): b) Model I(a):
c) Model II(a):
d) Model II(a):
e) Byth’s Model: f) Byth’s Model:
Residual Variance Residual Variance
58
Figure 2.5. Comparisons of biases in naï ve estimate(
) and bias-corrected estimate
(
) using Model I(a) and Model II(a) over ME parameter ( ) settings (0 - 0.4). Panel a)
and b) are results using Model I(a) under true parameter (
) equals to -0.2 and 0.2,
respectively. Panel c) and d) are results using Model II(a) under true parameter (
)
equals to -0.2 and 0.2, respectively. Follow-up visits number n = 5.
Bias%
0 0.1 0.2 0.3 0.4
Model I(a) Model II(a)
β
*
y0t
= -0.2 β
*
y0t
= 0.2
a)
c)
b)
d)
0 0.1 0.2 0.3 0.4
Bias%
Naï ve estimate ( ) Bias-corrected estimate ( )
0
yt
0
*
yt
59
Figure 2.6. Comparisons of biases in the naï ve estimate (
) and bias-corrected estimate
(
) using Model I(a) and Model II(a) over follow-up visits from 3 to 10. Panel a) and b)
are results using Model I(a) under true parameter (
) equals to -0.2 and 0.2,
respectively. Panel c) and d) are results using Model II(a) under true parameter (
)
equals to -0.2 and 0.2, respectively. ME = 0.3.
Bias%
Model I(a) Model II(a)
β
*
y0t
= -0.2 β
*
y0t
= 0.2
a)
c)
b) d)
Bias%
Naï ve estimate ( ) Bias-corrected estimate ( )
0
yt
0
*
yt
Follow-up visits (n) Follow-up visits (n)
2 3 4 5 6 7 8 9 2 3 4 5 6 7 8 9
60
Figure 2.7. Comparisons of biases in the naï ve estimate (
) and bias-corrected estimate
(
) using Model I(a) and Model II(a) over ME parameter ( ) settings (0 - 0.4) when
follow-up interval is varied from 8-16 months. Panel a) and b) are results using Model I(a)
under true parameter (
) equals to -0.2 and 0.2, respectively. Panel c) and d) are results
using Model II(a) under true parameter (
) equals to -0.2 and 0.2, respectively.
Follow-up visits number n = 5.
Bias%
0 0.1 0.2 0.3 0.4
Model I(a) Model II(a)
β
*
y0t
= -0.2 β
*
y0t
= 0.2
a)
c)
b) d)
0 0.1 0.2 0.3 0.4
Bias%
Naï ve estimate ( ) Bias-corrected estimate ( )
0
yt
0
*
yt
61
Figure 2.8. Comparisons of biases in naï ve the estimate (
) and bias-corrected estimate
(
) using Model I(a) and Model II(a) over follow-up visits from 3 to 10 when follow-
up interval is varied from 8-16 months. Panel a) and b) are results using Model I(a) under
true parameter (
) equals to -0.2 and 0.2, respectively. Panel c) and d) are results using
Model II(a) under true parameter (
) equals to -0.2 and 0.2, respectively. ME = 0.3.
Bias%
Model I(a) Model II(a)
β
*
y0t
= -0.2 β
*
y0t
= 0.2
a)
c)
b) d)
Bias%
Naï ve estimate ( ) Bias-corrected estimate ( )
0
yt
0
*
yt
Follow-up visits (n) Follow-up visits (n)
2 3 4 5 6 7 8 9 2 3 4 5 6 7 8 9
62
Figure 2.9. Comparisons of percent biases in the bias-corrected estimate (
) using
Model I(a) and Model II(a) over ME parameter ( ) settings when 50% data is missing at
random. Panel a) and b) are results for follow-up visits n = 5 under true parameter (
)
equals to -0.2 and 0.2, respectively. Panel c) and d) are results for follow-up visits n = 10
under true parameter (
) equals to -0.2 and 0.2, respectively.
Bias%
β
*
y0t
= -0.2 β
*
y0t
= 0.2
a) c)
b)
d)
Bias% Bias%
Bias%
n=5 n=10
0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4
Model I(a) Model II(a)
63
Chapter 3
Assessing the Baseline-adjusted Association between the
Exposure Variable and the Longitudinal Trajectory when
Responses Are Measured with Error
3.1 Literature review
3.1.1 Linear regression with one follow-up visit when responses are measured with
error
For studies with only a baseline and one follow-up visit, linear regression with the
outcome variable as the change from baseline value to the follow-up measurement can be
applied to investigate the association between the change association with the exposure
variable (e.g. gene, environmental factors, and so on) adjusting for the baseline value of
the outcome. Yanez et al. assessed the bias of the naï ve estimate using a linear regression
when responses are measured with error, and proposed the bias correction formula for the
naive estimate of the baseline-adjusted association between the exposure variable and the
change over time (Yanez et al., 1998). As shown by Yanez, the linear change model with
error-prone responses can be set as
, where
is the intercept, and
. Yanez provided the bias-
correction formula for the naï ve estimate
as
64
(3.1)
, and the bias-correction formula for the naï ve estimate
as
(3.2)
, where is the regression coefficient of the exposure using a linear regression of the
baseline value
regressed on , i.e.
;
is the ratio of variance of
measurement error to the variance of the observed baseline value conditioned on , i.e.
.
can be estimated using formula (2.1), and
can be estimated by
the residual variance of
regressed on , i.e.
, where is the residual
and
.
3.1.2 Mixed effects models analyzing the baseline-adjusted association the exposure
variable and the longitudinal trajectory for the study with multiple follow-up visits
If the study is designed with multiple follow-up visits and the hypothesis is to
assess the baseline adjusted association between the exposure variable and the
longitudinal trajectory, the mixed effects model shown in model (2.4) can be applied for
the analysis. Here, we assume the exposure to be time-invariant fixed effect without
measurement error and we start with a simpler scenario in that there is only one exposure
of interest other than the baseline value. Questions for multiple exposures of interest can
be solved following the same method we propose using necessary expansions. If the
baseline value of the outcome is associated with the exposure of interest and the baseline
value is also predicting its subsequent change over time, it is appropriate to adjust for
65
baseline value in the analysis of the exposure association with the longitudinal trajectory.
For such analysis, we can set a mixed effects model as
(3.3)
, where
,
,
, and
share the same definitions as in formula (2.6). is the
exposure of interest for its association with the longitudinal change.
is the parameter
for the cross-sectional association between successive measurements and the fixed effect
of exposure .
is the parameter representing the magnitude of deviation in the slope
of due to the difference in the exposure . For the same reason as formula (2.7),
can be excluded from the model because it is unidentifiable when
appears in the
model as a fixed effect. Thus, formula (3.3) can be shortened as
(3.4).
The mixed effects model (3.4) can be fitted in SAS Proc Mixed procedure (SAS Institute
Inc., Cary, NC, USA) or the R nlme package (Version 3.1-101) as no intercept model.
3.1.3 Harrison’s random intercept and random slope model for the baseline-
adjusted association between the exposure and the longitudinal trajectory when
responses are measured with error
As an extension of the Byth’s model, Harrison developed a random intercept and
random slope mixed effects model to estimate the exposure association with the change
over time adjusting for the baseline value of the outcome (Harrison et al., 2009). Harrison
built the model for true successive measurements
as:
66
, where ,
,
, and
.
and
are parameters for fixed effects of the time-
invariant exposure and the interaction of . Harrison shows
is biased compared
to the true association parameter for the change association with the exposure if is
associated with the baseline
and the baseline is also associated with the longitudinal
trajectory. Similar to the Byth’s model, the OLS estimate for the baseline association
with the longitudinal trajectory is constructed as:
. The true estimate of the
baseline-adjusted association between the exposure variable and the trajectory is derived
as:
.
For observed responses (
) measured with error, and when the measurement
error follows a normal additive model, Harrison’s model can be constructed as:
, where ,
,
, and
. Harrison argues that if the residual
only represents
the measurement error, the estimates derived from
and
are close to the true parameters for the longitudinal change associations with the baseline
and the exposure Z. However, Harrison’s method suffers the same limitation as Byth’s
method in that there are always some random variations in the residual due to the model
67
misspecification. More evidence of the limitation of Harrison’s method will be presented
in the simulation study.
3.2 Commonly used mixed effects model for the baseline-adjusted
association between the exposure variable and longitudinal trajectory
Following the same constructions of modeling in Chapter 2, there are also various
mixed effects models to be fitted when assessing the association between the exposure
and the longitudinal trajectory adjusting for the baseline value. Four different models can
be set as follows:
Model I(B):
, where ,
,
.
Model II(B):
, where ,
,
.
Model III(B):
, where
, ,
,
.
Model IV(B):
, where
, ,
,
.
Similar to Chapter 2, the difference between Model I(B) and Model II(B) is that
Model I(B) excludes the baseline value from the dependent variable, but Model II(B)
68
incorrectly kept the baseline in the dependent variable. Compared with Model I(B),
Model III(B) obtains all the same estimates, except
. Similarly, Model IV(B) derives
the same estimates as Model II(B), except
. Thus, we will only focus on the
discussion of
and
using Model I(B) and Model II(B).
When there is measurement error in the outcome, the naï ve parameter estimates of
associations between the baseline and the exposure of interest with the longitudinal
trajectory tend to be biased. We assume the same additive measurement error model as
Chapter 2 and defined the mixed effects model with error-prone outcomes corresponding
to Model I(B) and Model II(B) as
Model I(b):
, where ,
,
.
Model II(b):
, where ,
,
.
In this chapter, we will assess the bias in naï ve estimates of the change associations with
the exposure Z (
) and the baseline value (
) using Model I(b) and Model II(b).
69
3.3 Theoretical derivation of the bias correction formula for baseline-
adjusted association between the exposure variable and the longitudinal
trajectory
3.3.1 The bias-correction formula for Model I(b) when the baseline is exclude from
the dependent variable
Based on formula (2.5), we can rewrite the mixed effects model into a two-level
model. Model I(B) can be transformed as
Level 1:
(3.5)
Level 2:
.
Model I(b) can be transformed as
Level 1:
(3.6)
Level 2:
.
Because the baseline value is excluded from the outcome and formula (2.10) shows
, the naï ve estimate of the longitudinal slope ( )
is not correlated with the baseline measurement error
i.e.
in Model I(b).
Carroll et. al. showed that naï ve estimates using level 2 of formula 3.6 are biased as
(Carroll, Gallo and Leon Jay, 1985)
(3.7).
70
, where
is the variance of the true baseline value
;
is the variance of measurement
error;
is the variance of exposure Z; and
is the covariance between the true
baseline value and the exposure . After some algebra, we can get
, (3.8)
(3.9).
After we apply the same and
as Yanez’s paper, formula (3.8) and (3.9) can be further
written as
(3.10)
(3.11).
Solving formula (3.10) and (3.11) gives us the bias correction formula for
as
(3.12),
and the bias correction formula for
as
(3.13).
Based on the correction formula (3.12) and (3.13), we conclude that the bias in naï ve
estimates
and
of Model II(b) do not depend on the follow-up time, i.e. the
number of visits, because there is no time component in the correction formulae. formula
(3.12) has a similar mathematical structure as formula (2.12), except is replaced by ’.
When the exposure Z is not correlated with the baseline value, the naï ve estimate (
) is
not biased, even the responses are subject to error. Because the correction formula for
naï ve estimates do not depend on the time variable (
), Formulae (3.12) and (3.13) can
71
be directly generalized to a non-linear trajectory analysis using polynomial time effects.
Let
be the power-transformation of time variables (e.g.
or log(
)), and we can
then write the general correction formula as
(3.14),
(3.15).
3.3.2 Bias-correction formula for Model II(b) when the baseline is included in the
dependent variable
Although Model II(b) is an incorrect analysis model, we derive the bias-correction
formula if researchers still use it in their study. Similar to Model I(b), we can rewrite the
mixed effects model into a two-level model as (3.6). Because baseline measurement is
included in the outcome
, formula (3.7) should be changed as
(3.16)
, where
is the covariance of the longitudinal slope
and the baseline
measurement error (
). By solving formula (3.17), we get
(3.18)
72
, and
(3.19)
.
As shown in formula (2.17),
. Thus,
, where
. Using the same definition in Yanez’s paper, let
, and the least square estimate
. Because we assume there
is no measurement error in the exposure Z,
and
. Also,
. After using
parameters defined above, formula (3.18) can be written as
(3.19).
Formula (3.19) can be transformed as
(
) (3.20).
By solving both formulas (3.19) and (3.20), we get the bias correction formula for
as
(3.21).
The bias correction formula for
is
(
) (3.22).
The bias-correction formula for the naï ve estimate of the baseline association with
the longitudinal trajectory is almost the same as formula (2.23), except is replaced by ’.
73
The estimate of the true baseline association with the longitudinal trajectory is
independent of the parameter estimate of the trajectory association with the exposure Z.
In contrast, the bias in the naï ve estimate of the trajectory association with the exposure Z
depends on the true association between the baseline and the longitudinal trajectory over
time (
); the correlation between the baseline value and the exposure Z ( ); the ratio of
measurement error variance to the conditional variance of observed baseline value
conditioned on the exposure Z (
); and the number of follow-up visits ( ). Obviously, if
there is no association between the exposure Z and the baseline value, i.e. , the
naï ve estimate of the association between the exposure Z and the longitudinal trajectory is
not biased even if responses are measured with error. When there are only two time
points in the study and the follow-up time interval is fixed among subjects, is equal to
one. In this situation, bias correction Formulae (3.21) and (3.22) using a mixed effects
model are the same as the correction formulae derived by Yanez using a linear regression.
Also, we can explore our bias-correction formulae to any polynomial time effects to
allow the analysis of non-linear trends. Assuming
is the power-transformation of the
time variable (e.g.
or log(
)), we can then write the general correction formulae as
(3.23),
and
(
) (3.24),
, where
and
=
.
74
3.4 Simulation study
3.4.1 Simulation design
For simulation studies, we generated sample sizes of 50, 100, and 500 per cohort
to mimic the sample size of usual clinical studies and epidemiological studies. For each
set of underlying parameters, we repeated the simulation 1000 times and the results were
summarized by the average of 1000 replications. The bias of parameter estimates were
calculated by subtracting the true parameter from the naï ve parameter estimates obtained
from the simulated data. The percent bias was calculated as the division of the absolute
bias by the true parameter times 100. The true baseline
and the exposure Z were
generated using a multivariate normal distribution with both means equal to 0 and
variances equal to 1. The correlation coefficient ( ) is given among (-0.1, -0.2, -0.3, -0.4,
-0.5, 0, 0.1, 0.2, 0.3, 0.4, and 0.5). We did not choose a correlation coefficient greater
than 0.6 or less than -0.6 because if
and Z are highly correlated with each other, they
should not be included as covariates in one model in order to avoid the problem of
collinearity. The problem of collinearity is that estimates of regression parameters can be
unreliable when there are strongly correlated covariates in one model. Random sampling
with multivariate normal distribution can be generated by the RANDNORMAL statement
in SAS/IML or the “mvtnorm” package in R. Besides continuous exposure Z, we also
generated categorical exposure Z’ based on the order of the continuous Z generated above
and the given frequencies of different levels of categorical Z’. The process includes the
following: First, sorting Z in ascending order; and Second, cutting the sorted Z into
different groups based on the given frequencies of each group. Dichotomized variable
75
(coded as 0 and 1) was firstly generated with frequencies of 50% for both 0 and 1
categories. Then, categorical variable
containing three levels was also generated. This
(coded as 0, 1, 2) was simulated to mimic alleles of a single-nucleotide polymorphism
(SNP) represented by aa, Aa, AA under the assumption of genetic additive model
following Hardy-Weinberg equilibrium. We assumed minor allele (a) frequency q to be
0.2 and 0.4. Then the major allele frequency p was 0.8 and 0.6 correspondently. The
frequencies of aa, Aa and AA were calculated by
, and
. For minor allele
frequency equaling to 0.2, the frequencies of aa, Aa and AA were set to be 0.04, 0.16,
and 0.64. For minor allele frequency identical to 0.4, the frequencies of aa, Aa and AA
were set to be 0.16, 0.48, and 0.36. It should be noted that the correlation between the
categorical exposures (
and
) and the baseline value
was not the same as the
given correlation between the continuous Z and the baseline
, but these two correlations
were close to each other. The additive measurement error model was applied in the same
way as the simulation study in Chapter 2. In the construction of the simulated data, we
still used to set the magnitude of measurement error. was defined in Chapter 2 as the
proportion of total variance of the observed baseline value due to measurement error, i.e.
. We set to be varied from 0 to 0.5 (0, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, and 0.5).
Since we assumed
followed the standardized normal distribution with variance equal
to one, the corresponding
was set to be varied from 0 to 1 (0, 0.1005, 0.2294, 0.3333,
0.5, 0.6547, 0.8165, and 1).
76
The slope of successive responses was assumed to be associated with both the
exposure Z and the baseline value
in a linear regression:
.
was the intercept of the longitudinal slope, which was set to be 0, 0.5 and -0.5.
was set to be varied among the combinations [(-0.5,-0.5), (-0.5,0), (-0.5, 0.5), (0.25,
-0.25). (0.25, 0.25), (0,0), (0.5, -0.25), (-0.5, 0.25), (0.25, -0.5), (-0.25, 0.5)].
was
assumed to be normally distributed with N(0,
).
was set to be 0.2 and 0.4. For the
purpose of validating the bias-correction formulae using mixed effects models with
categorical exposure Z’,
and
were also implemented instead of the continuous
variable Z in the regression to predict
. The same
settings as continuous Z were
applied for both
and
. In Chapter 2, we generated two settings for the follow-up
time. First, the follow-up visit was generated as equal time intervals and corresponding to
each annual visit. Second, time intervals varied from 8 to 16 months with a scheduled
visit at one year. Uniform distribution, , was used to generate the
follow-up time variable. The simulated data was generated with up to 10 follow-up visits
(j=1, 2, 3, …, 10). Based on the true baseline
and the longitudinal slope generated
above, the true follow-up measurement for th individual at th time was generated using
the following model:
, where
is randomly distributed with N(0,
).
was set to be 0.1, 0.5, 1, and 2. The corresponding observed measurement at each
visit was equal to the true measurement plus the measurement error
:
,
where
is identically, independently and normally distributed with mean 0 and
variance
. We additionally set 20% or 50% follow-up data missing at random as in
77
Chapter 2 using a uniform distribution at follow-up visits to assess the performance of
correction formulae (3.12), (3.13), (3.21) and (3.22) when there are data missing
completely at random.
Similar to Chapter 2, we verified that the estimation of
and
do not
depend on
and
. Therefore, we only presented the result using
and
to be 0
and 0.2 respectively. Model I(b), Model II(b) and the Harrison’s model were applied to
estimate the baseline-adjusted association between the exposure Z and the longitudinal
trajectory. An unstructured within-individual correlation matrix was used in the mixed
effects model. The bias correction Formulae (3.12) and (3.13) were applied to the
coefficients using Model I(b). Formulae (3.21) and (3.22) were applied to correct the bias
in naï ve estimates using Model II(b). The correlation ( ) between the baseline value and
the exposure Z was estimated by the coefficient of the linear regression of
against Z.
was estimated by the proportion of observed baseline variance conditioned on the
exposure Z due to the measurement error, i.e.
. The variance of the observed
baseline value conditioned on the exposure Z was estimated by calculating the variance
of the residual in the linear regression of the observed baseline value regressed on the
exposure Z. The sampling variance of measurement error was estimated by formula
(2.26). Besides the continuous variable Z, we also assessed correction formulae for
categorical exposures
and
. Using the simulated data, we compared the naive
estimates
and
with corresponding bias-corrected estimates
and
among
Model I(b), Model II(b) and the Harrison’s model. In the condition that there were only
78
two time points (j=1), we compared the true association estimates
and
using
bias-correction formulae combined with Model I(b) and Model II(b) to the linear
regression estimates corrected by the method of Yanez et. al. (Formulae (3.3) and (3.4)).
3.4.2 Simulation results
Because λ’, instead of λ, is used in the bias-correction formulae for the estimate of
the baseline-adjusted association between the exposure and the longitudinal trajectory,
we compare λ and λ’ to demonstrate their differences at different correlations between the
baseline
and the exposure Z. The comparison results are presented in Table 3.1. To be
consistent with the simulation settings, the variances of
and Z are both defined as one.
λ and λ’ are generally very close to each other. When λ is equal to 0, λ’ is also equal to 0
for all correlation coefficients (r). With the absolute value of the correlation coefficient (r)
increased, the discrepancy between λ and λ’ is also increased. The amounts of λ’ over λ
are the same for both positive and negative correlations as long as the absolute values of r
are the same.
Table 3.2 and Table 3.3 present some simulation results of naï ve estimates of the
change associations with the baseline value (
) and the continuous exposure Z (
)
using mixed effects Model I(b) for the simulated data of three successive measurements
with fixed time intervals. We elected five different settings of the true relationship
and
to represent combinations of different scenarios of the slope associations with the
baseline Y
0
and covariate Z. Table 3.2 shows that when there is no measurement error
( =0), the naï ve estimates of the association between the baseline and the longitudinal
trajectory are not biased. It also shows the bias in the naï ve estimate
compared to the
79
true parameter
increases when the measurement error ( ) increases, and when the
correlation between the baseline
and the exposure Z increases. Table 3.3 indicates that
the naï ve estimates
are not biased when there is no measurement error in the baseline
( =0). Additionally, it shows the observed
is not biased when there is no
correlation between the baseline
and the exposure Z (r=0), even the baseline
is
prone to measurement error. Furthermore, when
and Z are correlated, the increase of
either measurement error ( ) or the absolute value of the correlation coefficient (|r|)
between
and Z can enlarge the bias of
. Additionally, Table 3.2 and Table 3.3 both
show that when there is no relationship between the baseline and the longitudinal change
(
=0), the naï ve estimates of
and
are not biased for any correlations between
and Z even
is measured with error.
Table 3.4 and Table 3.5 present the simulation results of naï ve estimates
and
using the same settings as Table 3.2 and Table 3.3, but applying Model II(b) as the
analysis model.
and
from Model II(b) share most of the characteristics of the
corresponding estimates from Model I(b) explained in Table 3.2 and Table 3.3, except
when
=0, there are some biases in the naï ve estimates of
and
when
is
subject to measurement error.
3.4.3 Validation of proposed bias-correction formulae
Table 3.6, Figure 3.1, and Figure 3.2 present some examples of bias-corrected
estimates
and
using the bias-correction Formulae (3.12) and (3.13) for Model
I(b), and bias-correction Formulae (3.21) and (3.22) for Model II(b). The exposure Z was
80
treated as a continuous variable and the follow-up time interval is fixed at one year. The
naï ve estimates of the baseline adjusted association between the exposure Z and the
longitudinal trajectory (
) are all biased from the true parameters at different
combinations of the simulation parameters for both Model I(b) and Model II(b). The bias
in the naï ve estimate using Model II(b) is even larger than Model I(b) (Figure 3.1 and
Figure 3.2). In contrast, the proposed bias-correction formulae correct the naï ve estimates
to be very close to the true association parameter
.
Results using the simulated data with time-varied follow-up intervals are similar
to previous results using a fixed time interval (Figure 3.3). Our bias-correction formulae
also successfully corrected the bias in the naï ve estimate
for ordinal exposures
and
(Figure 3.4).
3.4.4 Comparisons of Model I(b), Model II(b) and Harrison’s method
Comparing Model I(b) and Model II(b) with Harrison’s method (Harrison et al.,
2009) (Figure 3.5), we observe bias-corrected estimates
using both Model I(b) and
Model II(b), and they are not varied by different residual variance
, and are both close
to the true parameter. When there is measurement error in the response variable, the bias
in the corrected estimate
using Model II(b) is slightly larger than the bias using
Model I(b). Harrison’s method is biased when there is residual variability beyond
measurement error, and the bias increases when the residual variance
becomes larger.
81
3.4.6. Comparisons of Model I(b), Model II(b) and Yanez’s method when there is
only one follow-up visit
When there is only one baseline and one follow-up visit, both the mixed effects
model and linear regression can be used to assess the baseline-adjusted association
between the exposure variable and the longitudinal trajectory. Table 3.5 compares the
bias-corrected estimate (
) using bias-correction Formulae (3.13) and (3.22) under a
mixed effects models framework with the bias-corrected estimates using Yanez’s method
under a linear regression framework. The simulation results are based on sample size of
100 with 1000 replications. Both data with fixed-time intervals and varied-time intervals
were analyzed. When the time interval varies from 8 to 16 months for one year visit,
linear regressions with two different outcomes were additionally compared: First, using
the raw change from the baseline to the follow-up as the outcome; Second, using the rate
of change (the raw change divided by the follow-up time) as the outcome. When the
follow-up time interval is fixed across subjects, regression coefficients using the raw
change and the rate of change as outcomes are identical. Therefore, we only present the
results from the model using the raw change as the outcome. For a fixed-time interval,
bias-corrected estimates (
) using correction formula (3.22) for Model II(b) and
Yanez’s method for linear regression are identical and very close to the true parameter
(
) for the association between the exposure Z and the longitudinal trajectory
conditioned on the baseline value (
) when the outcomes are prone to measurement error.
Model I(b) cannot be fitted for the data with fixed-time intervals due to insufficient data
variations for parameter estimations. For a varied-time interval, Yanez’s method provides
82
similar results as bias-correction formula (3.13) for Model I(b). Bias-corrected estimates
using formula (3.22) based on Model II(b) has larger bias compared to Model I(b) and
Yanez’s method. After using Yanez’s correction method, estimates from two linear
regression models using either the absolute change over time or the rate of change in the
outcome are both close to the true association parameter.
In conclusion, when the study has only a baseline and one follow-up visit, linear
regression is suggested to be used for the longitudinal analysis. When random
measurement error exists in the observed responses, Yanez’s bias-correction method
should be used to obtain the correct estimate of the true baseline-adjusted association
between the exposure variable and the longitudinal trajectory. Additionally, when the
follow-up time interval is varied across subjects, the mixed effects model excluding the
baseline value from the dependent variable (model I(b)) can also be used for the analysis,
and the bias-correction formula (3.13) should be used to correct the bias in the naï ve
estimate when response are prone to measurement error.
Table 3.1. Comparison between λ and λ' at different correlation coefficients (r)
a
λ
λ’
r=0 r= -0.1 r= -0.2 r= -0.3 r= -0.4 r= -0.5 r=0.1 r=0.2 r=0.3 r=0.4 r=0.5
0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
0.1 0.10 0.10 0.10 0.11 0.12 0.13 0.10 0.10 0.11 0.12 0.13
0.2 0.20 0.20 0.21 0.22 0.23 0.25 0.20 0.21 0.22 0.23 0.25
0.3 0.30 0.30 0.31 0.32 0.34 0.36 0.30 0.31 0.32 0.34 0.36
0.4 0.40 0.40 0.41 0.42 0.44 0.47 0.40 0.41 0.42 0.44 0.47
0.5 0.50 0.50 0.51 0.52 0.54 0.57 0.50 0.51 0.52 0.54 0.57
a. λ
is the proportion of total observed variance in outcome due to random measurement error. λ’ is the proportion of measurement error variance in the observed
baseline variance conditioned on factor Z. r is the correlation coefficient for the correlations between baseline measurement Y
0
and exposure Z.
b. λ’ is calculated for different λ values based on the simulation settings that
equals to 1 and
equals to 1.
83
84
Table 3.2. Simulation results of
from two follow-up visits
a
using Model I(b)
True Parameters
b
Model I(b) Parameter Estimates
c
(
)
λ r = 0 r = -0.1 r = -0.5 r = 0.1 r = 0.5
0 0 0 0.0005 -0.0002 0.0005 -0.0004 -0.0005
0 0 0.01 0.0000 -0.0011 0.0006 0.0004 0.0013
0 0 0.1 0.0004 -0.0001 -0.0005 0.0020 0.0000
0 0 0.2 -0.0016 -0.0001 0.0005 0.0015 0.0003
0 0 0.3 0.0008 -0.0006 -0.0010 0.0014 -0.0026
0 0 0.4 -0.0030 0.0011 0.0015 -0.0010 -0.0005
0 0 0.5 0.0007 0.0005 -0.0026 -0.0014 0.0004
0 -0.25 0 0.0006 0.0008 0.0007 -0.0009 0.0002
0 -0.25 0.01 -0.0010 0.0000 0.0001 0.0008 0.0005
0 -0.25 0.1 -0.0002 -0.0004 -0.0014 0.0008 -0.0013
0 -0.25 0.2 0.0013 -0.0019 0.0005 0.0007 0.0003
0 -0.25 0.3 0.0001 -0.0009 0.0011 0.0000 -0.0006
0 -0.25 0.4 -0.0004 -0.0002 -0.0009 0.0003 0.0005
0 -0.25 0.5 0.0003 0.0010 0.0000 -0.0002 -0.0011
0.5 0 0 0.4995 0.5001 0.5003 0.4989 0.4994
0.5 0 0.01 0.4942 0.4948 0.4942 0.4961 0.4933
0.5 0 0.1 0.4493 0.4498 0.4350 0.4511 0.4353
0.5 0 0.2 0.4000 0.3997 0.3753 0.3984 0.3749
0.5 0 0.3 0.3509 0.3491 0.3199 0.3483 0.3164
0.5 0 0.4 0.3013 0.3021 0.2653 0.3002 0.2652
0.5 0 0.5 0.2485 0.2490 0.2146 0.2516 0.2167
0.25 0.5 0 0.2511 0.2504 0.2493 0.2510 0.2501
0.25 0.5 0.01 0.2478 0.2473 0.2474 0.2479 0.2472
0.25 0.5 0.1 0.2258 0.2258 0.2166 0.2259 0.2179
0.25 0.5 0.2 0.1983 0.2009 0.1883 0.1981 0.1867
0.25 0.5 0.3 0.1754 0.1745 0.1579 0.1754 0.1582
0.25 0.5 0.4 0.1492 0.1497 0.1322 0.1484 0.1306
0.25 0.5 0.5 0.1253 0.1239 0.1078 0.1241 0.1085
-0.25 0.25 0 -0.2497 -0.2500 -0.2503 -0.2503 -0.2492
-0.25 0.25 0.01 -0.2473 -0.2473 -0.2459 -0.2477 -0.2459
-0.25 0.25 0.1 -0.2255 -0.2250 -0.2178 -0.2254 -0.2175
-0.25 0.25 0.2 -0.2000 -0.2005 -0.1888 -0.1978 -0.1855
-0.25 0.25 0.3 -0.1751 -0.1743 -0.1565 -0.1727 -0.1599
-0.25 0.25 0.4 -0.1510 -0.1485 -0.1320 -0.1484 -0.1322
-0.25 0.25 0.5 -0.1243 -0.1259 -0.1071 -0.1246 -0.1040
85
a. Data shown from sample size of 500 with 1000 replications. Follow-up intervals are defined as fixed for 1 year
and covariate Z is treated as continuous.
b. Parameters used in generating the data in simulation study:
is the true change per visit associated with per
unit difference of baseline value;
is the true change per visit associated with per unit difference in factor Z. λ
is
the proportion of total observed variance in outcome due to random measurement error; The residual variance
was set to be 0.01.
c. Mean of estimated change per visit associated with per unit difference of baseline value (
) from 1000
replications.
86
Table 3.3. Simulation results of
from two follow-up visits
a
using Model I(b)
True Parameters
b
Model I(b) Parameter Estimates
c
(
)
λ r = 0 r = -0.1 r = -0.5 r = 0.1 r = 0.5
0 0 0 -0.0006 0.0003 0.0004 0.0001 0.0004
0 0 0.01 0.0004 -0.0009 -0.0006 -0.0002 0.0001
0 0 0.1 -0.0021 0.0015 -0.0006 -0.0006 0.0000
0 0 0.2 -0.0010 -0.0006 0.0005 0.0002 -0.0016
0 0 0.3 0.0000 -0.0025 0.0023 -0.0014 -0.0006
0 0 0.4 0.0042 0.0001 -0.0005 -0.0017 0.0017
0 0 0.5 0.0006 0.0002 -0.0004 0.0003 0.0016
0 -0.25 0 -0.2508 -0.2490 -0.2488 -0.2498 -0.2504
0 -0.25 0.01 -0.2497 -0.2504 -0.2501 -0.2496 -0.2494
0 -0.25 0.1 -0.2506 -0.2503 -0.2490 -0.2490 -0.2478
0 -0.25 0.2 -0.2496 -0.2516 -0.2485 -0.2498 -0.2488
0 -0.25 0.3 -0.2498 -0.2504 -0.2485 -0.2505 -0.2490
0 -0.25 0.4 -0.2511 -0.2528 -0.2512 -0.2501 -0.2472
0 -0.25 0.5 -0.2484 -0.2504 -0.2473 -0.2511 -0.2487
0.5 0 0 0.0005 0.0003 -0.0002 0.0012 0.0011
0.5 0 0.01 -0.0004 0.0003 -0.0037 -0.0001 0.0033
0.5 0 0.1 0.0018 -0.0048 -0.0335 0.0053 0.0324
0.5 0 0.2 -0.0001 -0.0086 -0.0611 0.0092 0.0617
0.5 0 0.3 -0.0026 -0.0134 -0.0905 0.0142 0.0915
0.5 0 0.4 -0.0016 -0.0191 -0.1167 0.0172 0.1164
0.5 0 0.5 -0.0002 -0.0295 -0.1422 0.0256 0.1410
0.25 0.5 0 0.5004 0.5004 0.4997 0.4996 0.5005
0.25 0.5 0.01 0.4983 0.4986 0.4985 0.4995 0.5031
0.25 0.5 0.1 0.4997 0.4975 0.4834 0.5015 0.5155
0.25 0.5 0.2 0.4975 0.4929 0.4705 0.5034 0.5304
0.25 0.5 0.3 0.4997 0.4932 0.4541 0.5078 0.5467
0.25 0.5 0.4 0.4993 0.4894 0.4423 0.5102 0.5589
0.25 0.5 0.5 0.4991 0.4895 0.4320 0.5158 0.5699
-0.25 0.25 0 0.2500 0.2495 0.2486 0.2507 0.2494
-0.25 0.25 0.01 0.2501 0.2497 0.2520 0.2489 0.2474
-0.25 0.25 0.1 0.2490 0.2532 0.2665 0.2459 0.2344
-0.25 0.25 0.2 0.2515 0.2547 0.2815 0.2429 0.2169
-0.25 0.25 0.3 0.2506 0.2558 0.2970 0.2405 0.2042
-0.25 0.25 0.4 0.2501 0.2603 0.3084 0.2420 0.1892
-0.25 0.25 0.5 0.2501 0.2624 0.3207 0.2405 0.1767
87
a. Data shown from sample size of 500 with 1000 replications. Follow-up intervals are defined as fixed for 1 year
and covariate Z is treated as continuous.
b. Parameters used in generating the data in simulation study:
is the true change per visit associated with per
unit difference of baseline value;
is the true change per visit associated with per unit difference in factor Z. λ
is the proportion of total observed variance in outcome due to random measurement error; The residual variance
was set to be 0.01.
c. Mean of estimated change per visit associated with per unit difference of factor Z (
) from 1000 replications.
88
Table 3.4. Simulation results of
from two follow-up visits
a
using Model II(b)
True Parameters
b
Model II(b) Parameter Estimates
c
(
)
λ r = 0 r = -0.1 r = -0.5 r = 0.1 r = 0.5
0 0 0 0.0001 0.0002 0.0001 0.0000 -0.0003
0 0 0.01 -0.0050 -0.0053 -0.0065 -0.0047 -0.0064
0 0 0.1 -0.0499 -0.0511 -0.0649 -0.0498 -0.0641
0 0 0.2 -0.1005 -0.1009 -0.1249 -0.0998 -0.1247
0 0 0.3 -0.1495 -0.1508 -0.1827 -0.1503 -0.1821
0 0 0.4 -0.2010 -0.2021 -0.2357 -0.2018 -0.2359
0 0 0.5 -0.2496 -0.2510 -0.2862 -0.2521 -0.2857
0 -0.25 0 0.0003 0.0007 0.0003 -0.0005 0.0004
0 -0.25 0.01 -0.0052 -0.0048 -0.0065 -0.0049 -0.0063
0 -0.25 0.1 -0.0497 -0.0504 -0.0649 -0.0503 -0.0653
0 -0.25 0.2 -0.0997 -0.1018 -0.1246 -0.1005 -0.1246
0 -0.25 0.3 -0.1504 -0.1520 -0.1813 -0.1512 -0.1818
0 -0.25 0.4 -0.2006 -0.2013 -0.2358 -0.2019 -0.2354
0 -0.25 0.5 -0.2502 -0.2516 -0.2870 -0.2511 -0.2874
0.5 0 0 0.5000 0.5002 0.5000 0.4996 0.4999
0.5 0 0.01 0.4897 0.4898 0.4871 0.4906 0.4872
0.5 0 0.1 0.3999 0.3990 0.3709 0.3995 0.3700
0.5 0 0.2 0.3003 0.2982 0.2494 0.2989 0.2498
0.5 0 0.3 0.2003 0.1974 0.1381 0.1966 0.1355
0.5 0 0.4 0.0997 0.0989 0.0293 0.0984 0.0290
0.5 0 0.5 -0.0005 -0.0025 -0.0715 -0.0019 -0.0703
0.25 0.5 0 0.2503 0.2502 0.2495 0.2506 0.2498
0.25 0.5 0.01 0.2424 0.2423 0.2401 0.2423 0.2404
0.25 0.5 0.1 0.1760 0.1740 0.1534 0.1750 0.1528
0.25 0.5 0.2 0.0999 0.0997 0.0637 0.0977 0.0630
0.25 0.5 0.3 0.0246 0.0230 -0.0236 0.0230 -0.0227
0.25 0.5 0.4 -0.0502 -0.0505 -0.1031 -0.0515 -0.1028
0.25 0.5 0.5 -0.1249 -0.1278 -0.1780 -0.1275 -0.1791
-0.25 0.25 0 -0.2498 -0.2503 -0.2502 -0.2501 -0.2493
-0.25 0.25 0.01 -0.2523 -0.2524 -0.2528 -0.2524 -0.2531
-0.25 0.25 0.1 -0.2757 -0.2753 -0.2821 -0.2756 -0.2822
-0.25 0.25 0.2 -0.3004 -0.3005 -0.3125 -0.2991 -0.3119
-0.25 0.25 0.3 -0.3252 -0.3253 -0.3398 -0.3247 -0.3410
-0.25 0.25 0.4 -0.3511 -0.3505 -0.3670 -0.3492 -0.3676
-0.25 0.25 0.5 -0.3745 -0.3763 -0.3936 -0.3755 -0.3923
89
a. Data shown from sample size of 500 with 1000 replications. Follow-up intervals are defined as fixed for 1 year
and covariate Z is treated as continuous.
b. Parameters used in generating the data in simulation study:
is the true change per visit associated with per
unit difference of baseline value;
is the true change per visit associated with per unit difference in factor Z. λ
is the proportion of total observed variance in outcome due to random measurement error; The residual variance
was set to be 0.01.
c. Mean of estimated change per visit associated with per unit difference of baseline value (
) from 1000
replications.
90
Table 3.5. Simulation results of
from two follow-up visits
a
using Model II(b)
True Parameters
b
Model II(b) Parameter Estimates
c
(
)
λ r = 0 r = -0.1 r = -0.5 r = 0.1 r = 0.5
0 0 0 -0.0004 0.0002 0.0001 -0.0001 0.0002
0 0 0.01 0.0000 -0.0010 -0.0032 0.0002 0.0038
0 0 0.1 -0.0006 -0.0049 -0.0323 0.0049 0.0319
0 0 0.2 0.0002 -0.0096 -0.0619 0.0104 0.0616
0 0 0.3 -0.0006 -0.0158 -0.0902 0.0153 0.0907
0 0 0.4 0.0019 -0.0206 -0.1185 0.0200 0.1195
0 0 0.5 -0.0005 -0.0253 -0.1419 0.0252 0.1438
0 -0.25 0 -0.2502 -0.2495 -0.2496 -0.2496 -0.2503
0 -0.25 0.01 -0.2498 -0.2509 -0.2533 -0.2494 -0.2465
0 -0.25 0.1 -0.2499 -0.2542 -0.2820 -0.2446 -0.2169
0 -0.25 0.2 -0.2492 -0.2617 -0.3116 -0.2397 -0.1874
0 -0.25 0.3 -0.2501 -0.2649 -0.3396 -0.2342 -0.1592
0 -0.25 0.4 -0.2503 -0.2707 -0.3681 -0.2304 -0.1325
0 -0.25 0.5 -0.2497 -0.2748 -0.3922 -0.2241 -0.1068
0.5 0 0 0.0001 0.0004 -0.0002 0.0004 0.0003
0.5 0 0.01 -0.0002 -0.0008 -0.0067 0.0005 0.0065
0.5 0 0.1 0.0007 -0.0094 -0.0641 0.0105 0.0650
0.5 0 0.2 -0.0002 -0.0197 -0.1245 0.0198 0.1249
0.5 0 0.3 -0.0011 -0.0290 -0.1796 0.0301 0.1828
0.5 0 0.4 -0.0002 -0.0403 -0.2354 0.0399 0.2343
0.5 0 0.5 -0.0005 -0.0523 -0.2851 0.0511 0.2859
0.25 0.5 0 0.5000 0.5001 0.4998 0.4997 0.5002
0.25 0.5 0.01 0.4991 0.4988 0.4949 0.5008 0.5056
0.25 0.5 0.1 0.4999 0.4921 0.4523 0.5075 0.5485
0.25 0.5 0.2 0.4998 0.4840 0.4069 0.5142 0.5934
0.25 0.5 0.3 0.5002 0.4782 0.3630 0.5230 0.6361
0.25 0.5 0.4 0.5005 0.4699 0.3237 0.5303 0.6774
0.25 0.5 0.5 0.5008 0.4618 0.2863 0.5389 0.7149
-0.25 0.25 0 0.2499 0.2498 0.2492 0.2503 0.2498
-0.25 0.25 0.01 0.2500 0.2495 0.2487 0.2496 0.2513
-0.25 0.25 0.1 0.2497 0.2477 0.2336 0.2521 0.2666
-0.25 0.25 0.2 0.2513 0.2447 0.2188 0.2543 0.2809
-0.25 0.25 0.3 0.2501 0.2417 0.2048 0.2569 0.2952
-0.25 0.25 0.4 0.2504 0.2404 0.1921 0.2608 0.3082
-0.25 0.25 0.5 0.2506 0.2365 0.1781 0.2632 0.3219
91
a. Data shown from sample size of 500 with 1000 replications. Follow-up intervals are defined as fixed for 1 year
and covariate Z is treated as continuous.
b. Parameters used in generating the data in simulation study:
is the true change per visit associated with per
unit difference of baseline value;
is the true change per visit associated with per unit difference in factor Z. λ
is the proportion of total observed variance in outcome due to random measurement error; The residual variance
was set to be 0.01.
c. Mean of estimated change per visit associated with per unit difference of baseline value (
) from 1000
replications.
Table 3.6.Validation results from two follow-up visits
a
using bias-correction formula for Model I(b) and Modle II(b)
a
True Parameters
b
Estimated Parameters
c
Cohort Size
λ r Uncorrected (
) Bias-corrected (
) Uncorrected (
) Bias-corrected (
)
Model I(b)
500 0 -0.25 0.2 -0.5 0.0323 0.0010 -0.1850 -0.2477
500 -0.5 -0.25 0.2 -0.5 -0.4712 -0.5033 -0.1890 -0.2528
500 -0.25 0.5 0.2 -0.5 -0.3143 -0.2510 0.3745 0.5004
500 0 -0.25 0.4 0.5 -0.0593 0.0003 -0.1323 -0.2517
500 -0.5 -0.25 0.4 0.5 -0.5592 -0.5000 -0.1315 -0.2497
500 -0.25 0.5 0.4 0.5 -0.1322 -0.2507 0.2637 0.5017
Model II(b)
500 0 -0.25 0.2 -0.5 -0.0311 0.0007 -0.3122 -0.2486
500 -0.5 -0.25 0.2 -0.5 -0.5320 -0.5002 -0.3124 -0.2491
500 -0.25 0.5 0.2 -0.5 -0.3754 -0.2486 0.2495 0.5016
500 0 -0.25 0.4 0.5 0.0584 -0.0016 -0.3674 -0.2474
500 -0.5 -0.25 0.4 0.5 -0.4408 -0.5007 -0.3670 -0.2475
500 -0.25 0.5 0.4 0.5 -0.0146 -0.2522 0.0288 0.5061
a. Validation results are from the simulated data with fixed time intervals for baseline and two follow-up visits. Covariate Z is defined as continuous variable.
was set
to be 0 and residual variance
was set to be 0.01
b. Parameters used in generating the data in simulation study:
is the true change per visit associated with per unit difference of baseline value;
is the true change
per visit associated with per unit difference of factor Z; λ is the proportion of total observed variance in outcome due to random measurement error; r is the correlation
coefficient between covariate Z and baseline.
c.. Mean of estimated change per visit associated with per unit difference of baseline value or covariate Z with 1000 replications. Uncorrected
and
are the model
based estimates from random effect mixed model; bias corrected
is calculated using formula (3.12) for model I(b) and formula (3.21) for model II(b); bias corrected
is calculated using formula (3.13) for model I(b) and formula (3.22) for model II(b).
92
93
Table 3.7. Comparison of bias-corrected estimates
of mixed effects model and linear
regression
True Parameters
b
Bias-corrected Estimates
c
(
) – True Parameter(
)
λ
Fixed Time Interval Varied Time Interval
Model
II(b)
Yanez
1
Model
I(b)
Model
II(b)
Yanez
1
Yanez
2
0 0 0 0.0000 0.0000 0.0031 0.0003 0.0001 0.0000
0 0 0.05 0.0004 0.0004 -0.0038 -0.0035 -0.0014 -0.0005
0 0 0.1 0.0002 0.0002 0.0007 -0.0053 -0.0015 0.0000
0 0 0.2 -0.0013 -0.0013 -0.0165 -0.0110 -0.0008 0.0031
0 0 0.3 -0.0052 -0.0052 -0.0202 -0.0182 -0.0014 0.0048
0 0 0.4 -0.0022 -0.0022 -0.0011 -0.0345 -0.0104 -0.0018
0 0 0.5 -0.0007 -0.0007 -0.0065 -0.0428 -0.0055 0.0088
-0.4 0 0 -0.0010 -0.0010 0.0032 0.0006 0.0002 0.0003
-0.4 0 0.05 -0.0013 -0.0013 -0.0063 -0.0030 -0.0005 0.0003
-0.4 0 0.1 0.0000 0.0000 -0.0028 -0.0069 -0.0035 -0.0017
-0.4 0 0.2 0.0013 0.0013 -0.0252 -0.0124 -0.0016 0.0022
-0.4 0 0.3 -0.0010 -0.0010 -0.0021 -0.0168 -0.0016 0.0048
-0.4 0 0.4 -0.0044 -0.0044 0.0011 -0.0254 -0.0014 0.0074
-0.4 0 0.5 -0.0126 -0.0126 0.0031 -0.0438 -0.0069 0.0065
0.4 0 0 0.0005 0.0005 -0.0002 0.0002 0.0003 0.0002
0.4 0 0.05 0.0009 0.0009 -0.0029 -0.0028 -0.0006 0.0000
0.4 0 0.1 -0.0008 -0.0008 0.0015 -0.0038 -0.0003 0.0015
0.4 0 0.2 -0.0049 -0.0049 -0.0019 -0.0100 -0.0012 0.0027
0.4 0 0.3 -0.0017 -0.0017 -0.0031 -0.0174 -0.0004 0.0046
0.4 0 0.4 -0.0018 -0.0018 0.0029 -0.0305 -0.0066 0.0024
0.4 0 0.5 -0.0149 -0.0149 0.0251 -0.0449 -0.0100 0.0024
0.4 -0.2 0 0.0001 0.0001 0.0041 0.0002 0.0001 -0.0002
0.4 -0.2 0.05 0.0005 0.0005 0.0022 -0.0010 0.0004 0.0015
0.4 -0.2 0.1 -0.0003 -0.0003 0.0084 -0.0040 -0.0008 0.0006
0.4 -0.2 0.2 -0.0013 -0.0013 -0.0021 -0.0117 -0.0027 0.0006
0.4 -0.2 0.3 0.0022 0.0022 0.0106 -0.0159 -0.0014 0.0045
0.4 -0.2 0.4 -0.0061 -0.0061 -0.0125 -0.0269 -0.0027 0.0080
0.4 -0.2 0.5 -0.0119 -0.0119 0.0134 -0.0417 -0.0061 0.0073
0.4 0.2 0 0.0003 0.0003 0.0013 0.0000 -0.0005 -0.0002
0.4 0.2 0.05 -0.0021 -0.0021 -0.0032 -0.0028 -0.0003 0.0001
0.4 0.2 0.1 -0.0021 -0.0021 -0.0052 -0.0068 -0.0026 -0.0010
0.4 0.2 0.2 0.0014 0.0014 0.0214 -0.0099 -0.0027 -0.0001
0.4 0.2 0.3 -0.0024 -0.0024 -0.0066 -0.0159 0.0007 0.0064
0.4 0.2 0.4 -0.0009 -0.0009 -0.0092 -0.0307 -0.0068 0.0029
0.4 0.2 0.5 -0.0143 -0.0143 0.0037 -0.0476 -0.0111 0.0019
94
a. Data shown from sample size of 100 with 1000 replications. Linear regression was used for one follow-up visit,
random effect mixed model was used for two and five follow-up visits.
b. Parameters used in generating the data in simulation study: n is the number of follow-up visits;
is the true
parameter for the baseline-adjusted association between the exposure Z and the longitudinal change;
is the
true change per visit associated with per unit difference of baseline value; is the proportion of total observed
variance in outcome due to random measurement error; The residual variance
was set to be 0.01.
c. Mean of bias-corrected estimated change per visit associated with per unit difference of the exposure Z (
) from
1000 replications. Correction formula (3.12) was used for model I(b). Correction formula (3.23) was used for
model II(b). Yanez’s method was used for linear regression of change model. Yanez
1
used raw change
as the
outcome. Yanez
2
used slope of
as the outcome.
95
Figure 3.1. Comparisons of biases in the naï ve estimate (
) and the bias-corrected
estimate (
) using Model I
B
and Model II
B
over different correlations (r) between the
baselineY
0
and exposure Z. Panel a) is the result using Model I
B
, and panel b) is the result
using Model II
B
. Follow-up visits number n = 4.True parameter
= - 0.2,
= - 0.4,
ME = 0.3.
a) Model I(b)
b) Model II(b)
Bias % Bias %
Correlation (r)
Naï ve estimate ( ) Bias-corrected estimate ( )
zt
*
zt
96
Figure 3.2. Comparisons of biases in the naï ve estimate(
) and the bias-corrected
estimate (
) using Model I
A
and Model II
A
over measurement error parameter ( )
settings (0 - 0.4). Panel a) and b) are the results of using Model IA under true parameter
(
) equals to -0.4 and 0.4, respectively. Panel c) and d) demonstrate the results of using
Model IIA under true parameter (
) equals -0.4 and 0.4, respectively. Follow-up visits
number n = 4,
= - 0.2, correlation (r) between baselineY
0
and exposure Z equals 0.2.
Bias%
0 0.1 0.2 0.3 0.4
Model I(b) Model II(b)
β
*
zt
= -0.4 β
*
zt
= 0.4
a)
c)
b)
d)
0 0.1 0.2 0.3 0.4
Bias%
Naï ve estimate ( ) Bias-corrected estimate ( ) zt
*
zt
97
Figure 3.3. Comparisons of biases in the naï ve estimate(
) and the bias-corrected
estimate (
) using Model I
A
and Model II
A
over measurement error parameter ( )
settings (0 - 0.4) when follow-up interval is varied from 8-16 months. Panel a) and b) are
results of using Model IA under true parameter (
) equals -0.4 and 0.4, respectively.
Panel c) and d) are results of using Model IIA under true parameter (
) equals -0.4 and
0.4, respectively. Follow-up visits n = 4,
= - 0.2, correlation (r) between baselineY
0
and exposure Z equals 0.2.
Bias%
0 0.1 0.2 0.3 0.4
Model I(b) Model II(b)
β
*
zt
= -0.4 β
*
zt
= 0.4
a)
c)
b)
d)
0 0.1 0.2 0.3 0.4
Bias%
Naï ve estimate ( ) Bias-corrected estimate ( ) zt
*
zt
98
Figure 3.4. Comparisons of biases in the naï ve estimate (
) and the bias-corrected
estimate (
) using Model I(b) and Model II(b) over measurement error parameter ( )
settings (0 - 0.4) when the exposure Z is categorical variable. Z is categorized into 0 and 1
by frequencies of 0.5 and 0.5 for panel a) and c) (Z
D
), and is categorized into 0,1 and 2 by
frequencies of 0.04, 0.32 and 0.64 for panel b) and d) (Z
G
). Panel a) and b) are results
using Model I(b). Panel c) and d) are results of using Model II(b). Follow-up visits
number n = 5,
= - 0.2,
= - 0.4, correlation (r) between the baselineY
0
and exposure
Z equals 0.2.
Bias%
0 0.1 0.2 0.3 0.4
Model I(b) Model II(b)
a)
c)
b)
d)
0 0.1 0.2 0.3 0.4
Bias%
Naï ve estimate ( ) Bias-corrected estimate ( ) zt
*
zt
Z
D
Z
G
99
Figure 3.5 Comparisons of biases in the bias-corrected estimate (
) using Model I
B
,
Model II
B
and Harrison’s method over different residual variations (
) settings (0.005,
0.01, 0.02, 0.04, 0.06, 0.08, and 0.1). Panel a) is the result when there is no ME ( ),
and panel b) is the result when . Follow-up visits number n = 5,
= - 0.2,
=
- 0.4, correlation (r) between baselineY
0
and exposure Z equals 0.2.
a) λ = 0
b) λ = 0.3
0 0.02 0.04 0.06 0.08 0.1
Bias % Bias %
2
Model I(b) Model II(b) Harrison’s method
100
Chapter 4
Obtaining bias corrected p-values and 95% confidence
intervals
4.1. Theoretical derivations of test statistics
4.1.1. Bootstrapping method
The bootstrapping method (Efron, 1982) is proposed to obtain p-values and 95%
confidence intervals to test whether the bias-corrected estimates of the baseline
association with the longitudinal trajectory and the baseline-adjusted exposure
association with the longitudinal trajectory are significantly different from null
hypotheses of no associations. First, a pseudo cohort is generated with the sample size
equal to the original study cohort where each subject is randomly selected from the
original study cohort with replacement. Second, data analysis is performed on this pseudo
cohort using the mixed effects models as described in Chapter 2 and Chapter 3, and bias
correction formulae (2.12), (2.23), (3.14), (3.15), (3.21), and (3.22) are applied to obtain
bias corrected estimates according to various mixed effects model settings. This process
is usually repeated 100 to 10000 times. We used 10000 times in our simulation study and
real data applications. The standard deviation of the corrected estimates from the repeated
pseudo cohort is the standard error of the corrected estimate from the original cohort. The
p-value can be obtained from the standard z-test if the distribution of the estimate is
101
normal. The 95% confidence interval is equal to the 2.5
th
percentile and 97.5
th
percentile
of the distribution of the corrected estimates from the repeated pseudo cohort.
4.1.2. Delta method
For a large dataset, the bootstrapping method is not efficient for data analysis.
Thus, we developed closed forms of test statistics using the delta method (Casella and
Berger, 2002). However, estimates of variances of parameter estimates require some
approximations where estimates of measurement error parameter
or
, and association
estimates between the baseline value and the exposure of interest ( ) are treated as fixed
values without variations. When using Model I(a) or Model II(a) for the analysis of the
change association with the baseline,
[4.1].
When using Model I(b) or Model II(b) for the analysis of the baseline-adjusted
association between the exposure variable and the longitudinal trajectory,
[4.2],
[4.3].
Estimates of
,
, and
can be obtained from the variance-
covariance matrix of fixed-effect parameter estimates using the mixed effects model. The
SAS Proc Mixed procedure provides COVB option to present such matrix, then p-values
can be obtained from the standard z-test.
102
4.2 Simulation results
We used the data generated in Chapter 3, and compared variances of parameter
estimates using the bootstrapping method and the delta method. Overall, estimates of
and
variances using Model I(b) or Model II(b) were very similar between the
bootstrapping method and the delta method for data with a relatively large sample size
(Table 4.1). However, when the sample size is small (N=50), variance estimates of both
and
tended to be smaller comparing the delta method to the bootstrapping
method. Thus, the type I error rate using the delta method was inflated for the small
sample size. This happens even when there is no measurement error in the outcome.
However, when the sample size increased, the inflated type I error became around the
nominal value 0.05 for Model I(b) using delta method. In contrast, although the type I
error rate decreased slightly, it was still inflated above 0.05 for Model II(b) using the
delta method. The bootstrapping kept the type I error rate around 0.05 even for the small
sample size and data with measurement error.
Because the delta method is as effective as the bootstrapping method when the
sample size is relatively large, we compared the performance of test statistics using the
delta method between Model I(b) and Model II(b) for a sample size of 200. Z-statistics
for the naï ve estimates (
and
) or bias-corrected estimates (
and
) using
Model I(b)
successfully kept the type I error rate all around assigned significance level
0.05 across different measurement error parameter λ settings (Figure 4.1). The type I
error rate of z-statistics for the naï ve estimates
and
using Model II(b)
were
103
dramatically inflated when λ is greater than zero. After using bias correction formulae
(3.21) and (3.22), the type I error rate of z-statistics for
and
were closer to 0.05.
However, when λ was greater than 0.2, type I error rates were still above 0.06 but below
0.09. The power of z-statistics for both parameter estimates using Model II(b) was larger
than Model I(b) (Figure 4.2). However, this gain in power might be caused by an inflated
type I error of Model II(b). Another reason for lower power using Model I(b) was
because Model I(b) excluded the baseline from the dependent variable, thus containing
fewer data points than Model II(b) to be used in the association estimation. However, this
loss of information became smaller when the number of follow-up visits increased
(Figure 4.3).
Table 4.1. Comparison of test statistics for
and
using delta method and bootstrapping method
a
.
N
Var
(
)
c
%Sig
(
)
b,c
Var
(
)
d
%Sig
(
)
b,d
Var
(
)
c
%Sig
(
)
b,c
Var
(
)
d
%Sig
(
)
b,d
Model I(b)
0 0 0 50 0.00026 0.0014 8.1 0.00170 5.2 0.00112 0.00086 8.1 0.00103 5.3
0 0 0 200 0.00004 0.00038 5.6 0.00043 5.4 0.00003 0.00023 5.2 0.00022 5.4
0 0 0 500 -0.00058 0.00015 5.8 0.00018 5.5 -0.00006 0.00009 5.1 0.00010 5.1
0 0 0.3 50 -0.00196 0.00564 7.8 0.00992 5.1 0.00200 0.00205 8.0 0.00242 4.2
0 0 0.3 200 0.00166 0.00121 5.7 0.00144 5.5 0.00073 0.00050 5.6 0.00051 5.4
0 0 0.3 500 0.00065 0.00048 5.1 0.00044 5.3 0.00045 0.00020 5.2 0.00020 4.9
-0.02 -0.04 0 50 -0.02314 0.00133 10.7 0.00137 9.8 -0.04101 0.00086 31.3 0.00096 31.2
-0.02 -0.04 0 200 -0.02107 0.00032 22.5 0.00034 19.8 -0.04023 0.00019 80.9 0.00022 78.2
-0.02 -0.04 0 500 -0.02021 0.00013 44.4 0.00014 39.8 -0.04002 0.00008 99.6 0.00008 99.5
-0.02 -0.04 0.3 50 -0.02422 0.00564 8.8 0.00947 4.5 -0.04041 0.00207 17.3 0.00216 13.4
-0.02 -0.04 0.3 200 -0.01789 0.00122 9.9 0.00125 8.3 -0.04121 0.00050 44.9 0.00047 47.6
-0.02 -0.04 0.3 500 -0.02024 0.00055 10.2 0.00042 13.8 -0.03967 0.00023 77.1 0.00020 81.9
Model II(b)
0 0 0 50 0.00033 0.00060 8.2 0.00068 5.5 0.00057 0.00036 8.3 0.00042 6.5
0 0 0 200 0.00007 0.00015 5.8 0.00017 4.6 0.00006 0.00008 5.6 0.00009 5.1
0 0 0 500 -0.00031 0.00006 4.8 0.00007 5.4 -0.00009 0.00004 5.3 0.00004 5
0 0 0.3 50 0.00725 0.00240 9.7 0.01169 3.4 0.00252 0.00108 9.3 0.00154 4.6
0 0 0.3 200 0.00230 0.00052 7.3 0.00071 6.8 0.00078 0.00021 6.7 0.00024 4.9
0 0 0.3 500 0.00073 0.00020 6.7 0.00024 4.7 0.00012 0.00008 6.1 0.00009 4.7
-0.02 -0.04 0 50 -0.02201 0.00059 16.5 0.00060 15.4 -0.04066 0.00036 58.1 0.00040 56.8
-0.02 -0.04 0 200 -0.02075 0.00015 40.8 0.00014 41.5 -0.03991 0.00009 98.7 0.00009 98.6
-0.02 -0.04 0 500 -0.02013 0.00006 73.8 0.00006 73.5 -0.04003 0.00004 100 0.00003 100
-0.02 -0.04 0.3 50 -0.01365 0.00231 11.6 0.00726 8.9 -0.03974 0.00086 33.5 0.00122 21.9
-0.02 -0.04 0.3 200 -0.01649 0.00050 17.4 0.00061 12.9 -0.04066 0.00020 78.9 0.00021 77.7
-0.02 -0.04 0.3 500 -0.01940 0.00019 31.2 0.00010 36 -0.03986 0.00008 96.6 0.00007 97.9
104
105
a) Simulation results are averages from 1000 replications of sample size (N) from 50 to 500. The number
of follow-up visits n=5, and the data is complete.
b) %sig represents the proportions of significance tests in all 1000 tests. When true parameter
and
equal to zero, %sig represents type I error rate of test statistic. When true parameters
and
do not equal to zero, % sig represents the power of test statistic.
c) Variance estimates and % sig are estimated using delta method.
d) Variance estimates and % sig are estimated using bootstrapping method.
106
Figure 4.1. Comparisons of the type I error rate of z-statistics of naï ve estimates (
,
and
) and the bias-corrected estimates (
, and
) using Model I(b) and Model II(b)
over ME parameter ( ) settings (0 - 0.4). Panels a) and b) are comparisons of
and
using Model I(b) and II(b). Panels c) and d) are comparisons of
and
using
Model I(b) and II(b). Follow-up visits number n = 5,
= 0,
= 0, correlation (r)
between baselineY
0
and exposure Z equals 0.2.
Type I error (%) Type I error (%)
λ λ
0
yt
zt
a)
b)
c)
d)
Naï ve estimates Bias-corrected estimates
Model I (b) Model II (b)
107
Figure 4.2 Comparisons of the power of z-statistics of naï ve estimates (
, and
) and
the bias-corrected estimates (
, and
) using Model I(b) and Model II(b) over ME
parameter ( ) settings (0 - 0.4). Panels a) and b) are comparisons of
and
using
Model I(b) and II(b). Panels c) and d) are comparisons of
and
using Model I(b)
and II(b). Follow-up visits number n = 5,
= - 0.2,
= - 0.4, correlation (r) between
baselineY
0
and exposure Z equals 0.2.
Power (%)
Power (%)
λ λ
0
yt
zt
Naï ve estimates Bias-corrected estimates
a)
b)
c)
d)
Model I (b) Model II (b)
108
Figure 4.3 Comparisons of power of z-statistics of the naï ve estimates (
, and
) and
the bias-corrected estimates (
, and
) using Model I(b) and Model II(b) over the
follow-up visits n=3-10. Panels a) and b) are comparisons of
and
using Model
I(b) and II(b). Panels c) and d) are comparisons of
and
using Model I(b) and II(b).
λ = 0.3,
= - 0.2,
= - 0.4, correlation (r) between baselineY
0
and exposure Z equals
0.2.
Power (%)
Power (%)
0
yt
zt
Naï ve estimates Bias-corrected estimates
a)
b)
c)
d)
3 4 5 6 7 8 9 10 3 4 5 6 7 8 9 10
Follow-up visits (n) Follow-up visits (n)
Model I (b) Model II (b)
109
Chapter 5
GDM cohort study
5.1 Introduction
The development of type 2 diabetes is characterized by worsening glucose
tolerance and falling beta-cell compensation for insulin resistance (Cnop et al., 2007;
Festa et al., 2006; Weyer et al., 2001; Xiang et al., 2006). Previous studies have
demonstrated that caloric restriction can delay the onset of type 2 diabetes (Knowler WC
et al., 2002; Lindströ m et al., 2003). Caloric restriction can attenuate obesity and reduce
insulin resistance (Despres et al., 1991; Larson-Meyer et al., 2006). Caloric restriction
can also enhance glucose tolerance and ameliorate insulin resistance independent of
weight loss (HENRY, SCHEAFFER and OLEFSKY, 1985; Kelley et al., 1993).
Severe, short-term caloric restriction has been shown to improve insulin
sensitivity and enhance beta-cell function (Lim et al., 2011; Utzschneider et al., 2004).
However, little is known about the relationship between spontaneous caloric intake and
chronic changes in beta-cell function under free-living conditions. In this report, we
examine this relationship using data from a cohort of non-diabetic Hispanic women with
recent gestations diabetes mellitus (GDM) who were prospectively followed for up to 12
years in the USC Gestational Diabetes Cohort Study.
110
5.2 Study design and method
5.2.1 The Study Participants
Selection of the original cohort has been described in detail (Buchanan et al.,
1998; Xiang et al., 1999). Briefly, all Latino women referred to Los Angeles County
Women’s Hospital for management of GDM between August 1993 and March 1995 were
asked to participate if they met all of the following criteria: 1) gestational age between 28
and 34 weeks, 2) no current or prior insulin therapy, 3) all fasting serum glucose
concentrations <130 mg/dl (7.2 mmol/l) during pregnancy, 4) otherwise uncomplicated
singleton pregnancy, 5) negative for anti-islet cell antibodies, and 6) both parents and at
least three of four grandparents from Mexico, Guatemala, or El Salvador. All women had
detailed metabolic testing during the third trimester (Xiang et al., 1999). They were asked
to return for a 75-g oral glucose tolerance test (oGTT) six months postpartum and then
for more detailed testing at 12-15 month intervals starting at fifteen months postpartum.
This testing included measurement of body morphometry and composition, assessment of
diet and physical activity (PA), and performance of oGTTs and intravenous glucose
tolerance tests (ivGTTs) as described below. Women who had abnormal fasting or post-
challenge glucose levels at a follow-up visit were referred to a dietitian for advice on
nutrition and daily walking. Participants remained in follow-up until they withdrew
consent, were lost to follow-up, reached the final scheduled study visit 12 years
postpartum, or developed a fasting plasma glucose concentration>140 mg/dl (>7.8
mmol/l), at which time they were referred for pharmacological treatment. Women who
were pregnant at the time of a scheduled follow-up visit were studied at least 4 months
111
after pregnancy and at least 1month after completion of breastfeeding. Family history of
diabetes was not assessed for this cohort.
All participants gave written informed consent for participation in the study,
which was approved by the Institutional Review Board of the University of Southern
California and the Los Angeles County plus University of Southern California Medical
Center.
For the present report, we analyzed data from all women who 1) had a baseline
oGTT and ivGTT without diabetes within 60 months postpartum and 2) returned for at
least one additional oGTT and ivGTT by twelve years postpartum.
5.2.2 Testing Procedures and Assays
For each follow-up assessment, women came to the General Clinical Research
Center on two separate days, at least 48 h apart, after an 8- to 12-h overnight fast and at
least three days on an unrestricted diet. On day 1, bioelectrical impedance (BIA) and
oGTT were performed (Xiang et al., 2010). On day 2, an IVGTT was performed (Xiang
et al., 2010). Glucose was measured by glucose oxidase (Beckman Glucose Analyzer II;
Beckman Coulter, Brea, CA). Insulin was measured by a radioimmunoassay (RIA) (Novo
Pharmaceuticals, Danbury, CT) that measured insulin and proinsulin with intra-and inter-
observer coefficient of variations (CV) of < 2.3 and 4.4%.
5.2.3 Dietary and PA Assessment
Trained bi-lingual interviewers administered both dietary and PA questionnaires.
Dietary intake was assessed by the structured food frequency questionnaires developed
for Hispanics by the Hawaii-Los Angeles Multi-ethnic Cohort Study (MEC) (Kolonel et
112
al., 2000; Stram et al., 2000). The questionnaire consists of a list of food items including
Hispanic-specific foods such as tamales and tortillas. Three choices of portion sizes, eight
frequency categories for foods and nine frequency categories for beverages of each item
were used to record the food intake during the past year. Total caloric and nutrient intakes
were computed from the percent contribution of the individual food items by the MEC
team. The questionnaire was validated and calibrated in multiethnic populations by
comparing diet reported from single questionnaire with three 24-hour recalls using the
same questionnaires (Stram et al., 2000).
The amount and intensity of PA was assessed by the questionnaires developed by
MEC (Kolonel et al., 2000; Nö thlings et al., 2007). The questionnaire was comprised of a
list of structured questions describing various types of activity (sitting, strenuous sports,
vigorous work, and moderate activities including sports and work) during the past year.
Responses were then used to estimate the total minutes of moderate and vigorous activity
per week. The total minutes of all activities was calculated as minutes of moderate plus
minutes of vigorous activities. The questionnaire was validated against doubly-labeled
water and showed a correlation coefficient of 0.3 between hours of total activity and total
energy expenditure (unpublished data).
5.2.4 Data Analysis
Insulin sensitivity (S
I
), and the acute insulin response to glucose from the first
10-mins (AIRg) under ivGTT were calculated from the ivGTT glucose and insulin
concentrations using the MIMMOD program (Bergman, 1989). Disposition index (DI), a
measure of beta-cell compensation for insulin resistance, was computed as the product of
113
AIRg and S
I
. Percent body fat and fat mass were calculated by the formula of Kotler et al.
(Kotler et al., 1996) using height, weight, and BIA measurements. PA levels were
compared to the recommendation of the World Health Organization’s (WHO): at least 75
minutes per week of vigorous or 150 minutes per week of moderate activity for
Americans.
Characteristics of the cohort were presented as medians, 25
th
, and 75
th
percentiles.
Fasting and post-challenge glucose and insulin concentrations, incremental insulin
responses, S
I
, DI, S
G
, and total calories were log transformed to approximate normal
distribution prior to analysis. An additional quadratic time effect was tested for evidence
against linear trend. If no significant quadratic time effect was observed, association
analyses were conducted using linear rates of change in traits. Mixed effects models were
used to assess the baseline prediction of the longitudinal change in oGTT and ivGTT
traits, as well as associations between baseline caloric intake and rates of change in oGTT
and ivGTT traits. Baseline characteristics, rates of change in body composition, caloric
intake, PA levels, and additional pregnancies during the follow-up were included as
covariates to adjust for potential confounding. Due to biological variations and
measurement errors inherent in oGTT and ivGTT traits, the association between the
baseline caloric intake and rates of changes in OGTT and IVGTT traits can be biased
when baseline values of oGTT and IVGTT traits are adjusted for. We corrected the
potential bias by applying our bias-corrected formulae described in Chapter 2 and
Chapter 3 for the parameter estimates of baseline associations with longitudinal
trajectories of metabolic traits, and baseline-adjusted associations between the baseline
114
caloric intake and rates of change in metabolic traits using mixed effects models. Random
variations in glucose, insulin, S
I
, AIRg, and DI were estimated from 109 Hispanic women
with similar characteristics who participated in a separate study where ivGTTs were
repeated after three months without any treatment (Buchanan et al., 2002). Additionally,
we assume the amount of measurement error in 2-hr glucose and insulin is same as
fasting glucose and insulin.
In order to minimize the distortion of sample estimates due to the small sample
size, we report bootstrapping association estimates and corresponding p-values (Efron
and Tibshirani, 1993). Sixty-two subjects were randomly selected with replacement from
the original sample 10000 times. The reported association estimates and their standard
errors were derived as the means and standard deviations of regression coefficients of all
replications. A Z-test was used to calculate the p-values for the estimates based on the
fact that the empirical distributions of our estimates are close to normal.
All statistical tests were two-sided. SAS version 9.2 (SAS Institute Inc., Cary,
North Carolina) was used for data analysis.
5.3 Simulation Study
Although we did a comprehensive simulation study in Chapter 2 and Chapter 3,
we further compared various methods using simulated data to mimic our real GDM
cohort study. We simulated the outcome and the exposure variables using the method
described in 3.4 to mimic oGTT, ivGTT traits, and the baseline caloric intake. Model I(b),
Model II(b), and Harrison’s method were applied to analyze the data, and biases in
115
parameter estimates were compared across three methods. We also validated our bias-
correction formulae (3.12), (3.13), (3.21), and (3.22) using this simulated data.
For the baseline association with longitudinal trajectories of outcomes over time,
the biases in the naï ve association estimates (
) were from -32% to 27%. After using
the bias-correction formula (3.12) for Model I(b), the biases in corrected estimates (
)
shrank to -6% to 8%. Biases in both naï ve estimates (-55% to 208%) and bias-corrected
estimates (-101% to 13%) using model II(b) were both higher than biases using model
I(b), due to the fact that there are 50% data missing at random. Estimates using
Harrison’s method are very unstable, and can be tremendously biased for some
simulation runs.
For the baseline-adjusted association between the exposure variable and
longitudinal trajectories in outcomes, the biases in the naï ve association estimates (
)
were from -11% to 13%. After using the bias-correction formula (3.22) for Model I(b),
the biases in corrected estimates (
) were attenuated to -5% to 5%. Biases in both naï ve
estimates (-13% to 19%) and bias-corrected estimates (-7% to 7%) using Model II(b)
were both higher than biases using Model I(b). Because
using Harrison’s method is
very unstable, the parameter estimates of
is also not stable and can be extremely
biased in some cases. It should be noted that the influence of bias-corrections on the
naï ve estimate
was not obvious in this data largely due to small correlations between
the exposure variable and baseline values of the outcomes.
116
5.4 Results
5.4.1Baseline Characteristics and rates of change during the follow-up
Sixty-two women met the inclusion criteria of baseline oGTTs, ivGTTs, and at
least one follow-up visit. At baseline (Table 5.1), the median age was 32 years (range:
21-42 years) and the median BMI was 30.3 kg/m
2
. Thirty-six (58%) of the participants
had impaired glucose tolerance. The median caloric intake was 2091 kcal/day, 53% from
carbohydrates, 33% from fat, and 17% from protein. Total fat and saturated fat
consumption were above the upper limits recommended by the WHO guideline for
preventing diet-related chronic diseases (WHO, 2003). Carbohydrate consumption was
below the lower limit of the WHO recommendation. Very few subjects were physically
active; the median duration of total activity was 0 minute per week. Only six (<10%) of
the women met the WHO recommendation of more than 75 minutes per week vigorous
activity or 150 minutes per week moderate activity adults (116).
The median duration of the follow-up was 7.5 years (range 1.3 – 10.8 years) with
a median of five sets of oGTT, ivGTT, and BIA per participant. During the follow-up, 27
(44%) of the women developed diabetes defined by American Diabetes Association
(Diagnosis and Mellitus*, 2003). A total of 15 (24%) of the women experienced one or
more additional pregnancies.
No significant quadratic trends were observed for any metabolic trait during the
follow-up (p>0.14). Thus, only linear rates of change are presented in Table 1. BMI,
fasting glucose, 2-hr glucose, fasting insulin, and moderate and total activity amounts
increased significantly over time (all p<0.001), while S
I
, AIRg, and DI decreased
117
significantly over time (p=0.0003, 0.001, <0.0001, respectively). No other significant
changes were found during the follow-up (all p>0.10). Dietary intake and vigorous
activity did not change over time.
5.4.2 Baseline prediction of longitudinal trajectories of oGTT and ivGTT traits
Table 5.2 shows Model I(a), Model II(a) and Byth’s method provide different
results for selected associations between baseline values of traits and longitudinal
trajectories of traits. Before the bias correction, we found higher baseline 2-hr glucose
was associated with faster increase in 2-hr glucose and higher baseline AIRg was
associated with faster decrease in AIRg over time using Model I(a) (p=0.022 and 0.034,
respectively). Compared to Model I(a), more negative baseline associations with rates of
change in insulin levels, S
I
, DI over time (all p<0.025) were found using Model II(a) and
Byth’s method. However, these additional negative associations are hard to interpret for
their biological meanings. Instead, Model II(a) and Byth’s method did not find the
significant association between the baseline 2-hr glucose and the rate of increase in 2-hr
glucose over time (both p>0.11).
Because random variations were expected for all oGTT and ivGTT measurements,
we estimated these random variations from another study with similar characteristics as
the GDM cohort study (Buchanan et al., 2002). We assumed the true values of oGTT and
ivGTT measurements did not change within the three-month follow-up time. The
proportion of measurement error variance in the total variance (λ) of glucose levels,
insulin levels, SI, AIRg, and DI were 32%, 28%, 16%, 14%, and 18%, respectively. After
using bias-correction formulae (2.12 and 2.22) for Model I(a) and Model II(a),
118
magnitudes of association parameter estimates changed from 5% to 345%. Directions and
p-values of these association estimates were not changed using Model I(a) after
correcting biases. However, for Model II(a), bias correction changed significant baseline
associations with rates of change in 2-hr insulin and DI to be non-significant.
5.4.3 Associations between baseline caloric intake and rates of change in oGTT and
ivGTT traits after adjusting for baseline values of traits
Based on the naï ve estimates, we found higher baseline caloric intake was
associated with faster decreases in S
I
and DI after adjusting for baseline S
I
and DI using
Model I(b) (p=0.013 and 0.009, respectively). However, Model II(b) and Harrison’s
model did not find these significant associations. Additionally, we a found borderline
significant association between the baseline caloric intake and the rate of increase in 2-hr
glucose (p=0.095) using Model I(b). After bias-correction, magnitudes of parameter
estimates were changed from 0.1% to 6% for Model I(b), and from 2% to 31% for Model
II(b). However, directions of associations and the corresponding p-values were not
changed after correcting the bias.
Based on our previous simulation results, Model I(a) and Model I(b) provided
more accurate association estimates after bias-correction than Model II(a) and Model
II(b), as well as Byth’s method and Harrison’s method. Thus we selected Model I(a) and
Model I(b) as our final analysis model. We also assessed the relative contribution of body
fat to the observed associations between the baseline caloric intake and metabolic traits
by additionally adjusting for BMI, percent body fat, or waist/hip ratio using Model I(b)
combined with bias-correction formula (3.13). After further adjustment for the rate of
119
change in BMI, significant associations between the baseline caloric intake and rates of
change in S
I
and DI remained (p= 0.013 and 0.003, respectively), and regression
coefficients changed by less than 7%. Adjustment for changes in BMI increased the
predictive effect of baseline caloric intake on the rate of change in 2-hr glucose by 20%
and the association turned to be significant (p=0.037). Results were similar when models
were adjusted for percent body fat or waist/hip ratio instead of BMI. The baseline value-
and BMI- adjusted means for 2-hr glucose, S
I
, and DI over follow-up by baseline low
versus high caloric intake groups were presented in figures 5.1.
Adjustment for additional confounders including baseline age, baseline parity,
rates of change in caloric intake, additional pregnancies during follow-up, baseline, and
rates of change in PA levels did not change the significant associations between baseline
caloric intake and rates of change in S
I
, DI, and 2-hr glucose (changes in regression
coefficients <15%).
5.5 Conclusions
In this study, we observed that higher baseline daily caloric intake was predictive
of the faster increase of 2-hr glucose, and faster declines of insulin sensitivity and beta-
cell function among Hispanic women at high risk of type 2 diabetes. These predictive
effects were not modified significantly by the adjustment for baseline values, the rate of
change in adiposity, and other confounders including baseline age, parities, physical
activity levels, and additional pregnancies during the follow-up. Taken together, our
120
findings indicate that higher caloric intake in a free-living environment may have
contributed to rising insulin resistance and falling beta-cell function.
This study served as an application example using the mixed effects model and
bias-correction formulae for the baseline association with the longitudinal trajectory, and
baseline-adjusted associations between the exposure variable and the longitudinal
trajectory. We compared Model I(a), Model II(a) with Byth’s method for the baseline
associations with longitudinal trajectories in oGTT and ivGTT traits. Also, we compared
Model I(b) and Model II(b) with Harrison’s method for baseline-adjusted associations
between baseline caloric intake and rates of change in metabolic traits.
In the simulation study, we compared bias in association estimates using various
methods based on data simulated to mimic the real data in the GDM cohort study.
Because Model I(b), Model II(b), and Harrison’s method include the baseline association
with the longitudinal trajectory as described in Model I(a), Model II(a) and Byth’s
method, we did not perform an additional simulation study comparing Model I(a), Model
II(a) and Byth’s method. We found Model I(b) had the most correct association
estimates for all traits after bias-correction. Although our bias-correction formulae had
some impact on estimates of baseline associations with longitudinal trajectories of traits,
they did not significantly change the estimates of baseline-adjusted associations between
the baseline caloric intake and rates of change in outcomes due to small correlations
between the baseline caloric intake and metabolic traits. If those correlations become
larger, our bias-correction formulae will show more alterations of association estimates.
121
Results from the real data showed different methods could lead to different
conclusions of the baseline associations with longitudinal trajectories of metabolic traits,
and baseline-adjusted associations between the baseline caloric intake and trajectories of
metabolic traits. We observed that Model II(a) and Byth’s method failed to find a positive
association between the baseline 2-hr glucose and its longitudinal changes over time.
Instead, Model II(a) and Byth’s method found more significantly negative associations
between baseline insulin levels, S
I
, DI, and their rates of change over time. However,
these negative associations do not have biological meanings. As we showed in the
simulation study, association estimates using Model II(a) have more bias than Model I(a)
when there were data missing at random. However, using bias-correction formula (2.22)
combined with Model II(a) successfully detected a significant association for the 2-hr
glucose, and corrected unexpected findings of 2-hr insulin and DI. Because residual
variations existed beyond measurement error variations in our data, estimates using
Byth’s method were biased. Due to the incorrectness of association estimates of Model
II(b) and Harrison’s method, both methods missed significant associations between the
baseline caloric intake and rates of change in S
I
and DI during the follow-up. Because
correlations between the baseline caloric intake and baseline traits were small, we did not
observe significant changes in association estimates after using bias correction formula.
In conclusion, Model I(a) and Model I(b) are suggested to be used for real data
analysis, especially for incomplete data. When random measurement error exists in the
outcome variable, it is important to correct the bias in the naï ve estimate of baseline
association with its longitudinal change over time. When another exposure variable is
122
correlated with the baseline value of the outcome, and the baseline value needs to be
adjusted for in the association between the exposure variable and the longitudinal change
in the outcome, bias-corrections should also be applied to the naï ve association estimate.
However, the alteration to the naï ve estimates after bias-correction might not be
significant when the correlation between the exposure variable and the baseline value of
the outcome is small.
Table 5.1. Validation results of data simulated to represent glucose and ivGTT traits.
Representative traits λ Var(Y
0
) r V ar( ε )
%bias (
) %bias (
)
%bias (
) %bias (
)
Model I(b)
Fasting Glucose 0.32 0.005 0.14 0.009 0.22 -31.1 7.7 0.03 13.0 -4.8
2-hr Glucose 0.32 0.029 0.02 0.023 0.12 -32.2 4.1 0.04 0.5 -0.9
Fasting Insulin 0.28 0.223 0.09 0.183 -0.02 -30.4 -0.2 0.01 -10.5 -3.5
2-hr Insulin 0.28 0.348 0.02 0.227 -0.01 27.3 3.1 -0.01 2.4 1.8
SI 0.16 0.112 -0.14 0.055 -0.03 -16.6 0.6 -0.05 0.6 2.0
AIRg 0.14 0.218 -0.05 0.084 -0.05 -13.0 3.0 -0.04 3.9 4.8
DI 0.18 0.093 0.01 0.024 -0.05 -24.3 -6.0 -0.03 0.2 0.3
Model II(b)
Fasting Glucose 0.32 0.005 0.14 0.009 0.22 -43.6 13.3 0.03 19.0 -7.3
2-hr Glucose 0.32 0.029 0.02 0.023 0.12 -55.0 13.4 0.04 0.9 -2.0
Fasting Insulin 0.28 0.223 0.09 0.183 -0.02 50.2 -30.5 0.01 -12.8 6.5
2-hr Insulin 0.28 0.348 0.02 0.227 -0.01 207.5 -101.0 -0.01 -10.3 -4.2
SI 0.16 0.112 -0.14 0.055 -0.03 23.8 -12.8 -0.05 3.7 0.8
AIRg 0.14 0.218 -0.05 0.084 -0.05 12.8 -6.1 -0.04 3.4 2.4
DI 0.18 0.093 0.01 0.024 -0.05 8.7 -13.5 -0.03 -1.5 -1.4
Harrison's estimates
Fasting Glucose 0.32 0.005 0.14 0.009 0.22
3.8E+32 0.03 43.5 -2.1E+32
2-hr Glucose 0.32 0.029 0.02 0.023 0.12
3.3E+18 0.04 2.9 -5.5E+17
Fasting Insulin 0.28 0.223 0.09 0.183 -0.02
-1.1E+16 0.01 -28.1 -1.1E+16
2-hr Insulin 0.28 0.348 0.02 0.227 -0.01
-3.7E+16 -0.01 -8.0 -6.1+E15
SI 0.16 0.112 -0.14 0.055 -0.03
-6.9E+32 -0.05 -6.2 -4.0E+28
AIRg 0.14 0.218 -0.05 0.084 -0.05
-7.9E+16 -0.04 -2.9 -9.8E+15
DI
0.18 0.093 0.01 0.024 -0.05 -1.7E+32 -0.03 -1.8 1.9E+31 123
124
*Data was simulated as sample size of 60 with 10 visits. 50% of data was set to be missing at random. Parameters used
in generating the data in simulation study: variance of baseline caloric intake = 0.28; was the proportion of total
observed variance in outcome due to random measurement error; Var(Y
0
) was the variance of true baseline values; r
was the correlation coefficient between baseline traits and baseline caloric intake. Var(ε) was the residual variance.
The association estimates (
,
,
,
) were estimated by averaging the results from 1000 replications. %bias
in the estimates was calculated as (regression estimate – true association parameter)/true association parameter× 100.
125
Table 5.2. Baseline characteristics and rates of change during the follow-up
Baseline
Annualized Rates of
change
Median
(25
th
, 75
th
percentiles)
Median
(25
th
, 75
th
percentiles)
Anthropometrics
Age [years] 32 (28, 36)
BMI [kg/m
2
] 30.3 (28.0, 32.6) 0.30 (0.06, 0.59)‡
Body fat (%) 44.3 (40.6, 48.4) 0.36 (-0.06, 0.66)‡
Weight [kg] 70.0 (64.1, 78.9) 0.82 (0.20, 1.38)‡
Fat Mass [kg] 31.9 (26.8, 37.0) 0.59 (0.04, 1.08)‡
Waist Hip ratio × 100 0.83 (0.79, 0.85) 0.006 (0.0004, 0.011)‡
oGTT-derived Measures
Fasting glucose [mmol/l]* 5.3 (5.0, 5.7) 0.02 (0.003, 0.04)‡
2hr glucose [mmol/l]* 8.0 (6.8, 9.5) 0.03 (0.01, 0.09)‡
Fasting insulin [pmol/l]* 115 (87, 167) 0.06 (-0.03, 0.09)‡
30-min ∆ Insulin [pmol/l]* 639 (483, 906) -0.04 (-0.11, 0.02)
2hr Insulin [pmol/l]* 695 (514, 1285) 0.002 (-0.117, 0.063)
ivGTT-derived Measures
S
I
[× 10
-3
min
-1
per pmol/l]* 1.4 (1.0, 2.2) -0.03 (-0.10, 0.04)†
AIRg [pmol/l × 10 min]* 3167 (1854, 5372) -0.05 (-0.11, -0.01)‡
DI [S
I
× AIRg]* 5172 (2374, 8214) -0.04 (-0.09, -0.01)‡
Dietary Intake
Calorie Intake [kcal/day]* 2091 (1526, 3013) -0.02 (-0.07, 0.03)
Percentage of Total Energy intake
Total Fat 32.8 (28.9, 37.6) 0.04 (-0.64, 0.70)
Saturated Fatty Acids 10.9 (9.3, 12.7) -0.01 (-0.37, 0.23)
Polyunsaturated Fatty Acids 6.8 (6.0, 7.7) -0.01 (-0.15, 0.12)
Monounsaturated Fatty Acids 11.8 (10.4, 13.5) 0.02 (-0.29, 0.26)
Total Carbohydrate 53.1 (46.3, 57.8) -0.14 (-0.98, 0.74)
Protein 16.7 (15.4, 18.3) 0.15 (-0.13, 0.40)
Physical Activity
Vigorous Activity [min per week] 0 (0, 0) 0.0 (-3.1, 13.9)
Moderate Activity [min per
week]
0 (0, 44.9) 3.5 (0, 28.2)‡
Total activity [min per week] 0 (0, 44.9) 8.0 (0, 70.5)‡
Meets WHO Guidelines 6 (10%)
n = 62. Rates of change are presented as medians (25
th
, 75
th
percentiles) of units per year. * Log transformation was
applied for data analysis. † P < 0.05, ‡ P < 0.001.
126
Table 5.3. Baseline prediction of longitudinal changes in oGTT and ivGTT traits.
λ
P
1
P
2
%Diff*
Model I(a)
Fasting glucose (mmol/l per year)
0.32 0.25 0.228 0.37 0.228 47.1
2-hr glucose (mmol/l per year) 0.32 0.08 0.022 0.12 0.022 47.1
Fasting insulin (pmol/l per year) 0.28 -0.02 0.423 -0.03 0.423 38.9
2-hr insulin (pmol/l per year) 0.28 0.01 0.806 0.01 0.806 38.9
S
I
(min
-1
per pmol/l× 10
-3
)
0.16 -0.03 0.402 -0.04 0.402 19.0
AIRg (pmol/l × 10 min) 0.14 -0.05 0.034 -0.06 0.034 16.3
DI
0.18 -0.02 0.370 -0.03 0.370 22.0
Model II(a)
Fasting glucose (mmol/l per year) 0.32 0.35 0.163 0.59 0.113 68.2
2-hr glucose (mmol/l per year) 0.32 0.02 0.504 0.11 0.046 345.4
Fasting insulin (pmol/l per year) 0.28 -0.09 <0.0001 -0.06 0.004 -29.2
2-hr insulin (pmol/l per year) 0.28 -0.06 0.0002 -0.03 0.291 -58.6
S
I
(min
-1
per pmol/l× 10
-3
)
0.16 -0.13 0.002 -0.12 0.011 -4.5
AIRg (pmol/l × 10 min) 0.14 -0.08 0.004 -0.07 0.038 -16.0
DI 0.18 -0.06 0.010 -0.04 0.159 -33.7
Byth's Method
Fasting glucose (mmol/l per year) 0.32 NA NA 0.54 0.117
2-hr glucose (mmol/l per year)
0.32 NA NA 0.08 0.115
Fasting insulin (pmol/l per year) 0.28 NA NA -0.07 0.025
2-hr insulin (pmol/l per year) 0.28 NA NA -0.08 0.018
S
I
(min
-1
per pmol/l× 10
-3
)
0.16 NA NA -0.11 0.014
AIRg (pmol/l × 10 min) 0.14 NA NA -0.06 0.015
DI
0.18 NA NA -0.05 0.011
P
1
is the p-value of naï ve estimate
.
P
2
is the p-value of bias-corrected estimate
.
*%Diff equals to
.
127
Table 5.4. Baseline-adjusted associations between baseline caloric intake and
longitudinal changes in oGTT and ivGTT traits.
λ r
2
P
1
P
2
%Diff*
Model I(b)
Fasting glucose (mmol/l per year)
0.32 0.01990 0.021 0.152 0.020 0.209 -5.9
2-hr glucose (mmol/l per year)
0.32 0.00042 0.024 0.095 0.025 0.087 4.9
Fasting insulin (pmol/l per year)
0.28 0.00893 0.020 0.456 0.020 0.453 2.1
2-hr insulin (pmol/l per year)
0.28 0.00039 -0.026 0.362 -0.026 0.366 -0.1
S
I
(min
-1
per pmol/l× 10
-3
)
0.16 0.01955 -0.064 0.013 -0.065 0.012 1.5
AIRg (pmol/l × 10 min)
0.14 0.00233 -0.022 0.405 -0.021 0.420 -3.6
DI
0.18 0.00006 -0.040 0.009 -0.040 0.009 0.1
Model II(b)
Fasting glucose (mmol/l per year)
0.32 0.01990 0.021 0.181 0.019 0.331 -13.0
2-hr glucose (mmol/l per year)
0.32 0.00042 0.018 0.205 0.020 0.177 14.4
Fasting insulin (pmol/l per year)
0.28 0.00893 0.008 0.640 0.007 0.692 -10.6
2-hr insulin (pmol/l per year)
0.28 0.00039 -0.011 0.588 -0.007 0.711 -31.4
S
I
(min
-1
per pmol/l× 10
-3
)
0.16 0.01955 -0.021 0.329 -0.020 0.342 -2.9
AIRg (pmol/l × 10 min)
0.14 0.00233 -0.023 0.370 -0.025 0.350 5.7
DI
0.18 0.00006 -0.022 0.130 -0.022 0.141 2.4
Harrison's Method
Fasting glucose (mmol/l per year)
0.32 0.01990 NA NA 0.023 0.165 NA
2-hr glucose (mmol/l per year)
0.32 0.00042 NA NA 0.021 0.070 NA
Fasting insulin (pmol/l per year)
0.28 0.00893 NA NA 0.005 0.765 NA
2-hr insulin (pmol/l per year)
0.28 0.00039 NA NA 0.003 0.830 NA
S
I
(min
-1
per pmol/l× 10
-3
)
0.16 0.01955 NA NA -0.003 0.878 NA
AIRg (pmol/l × 10 min)
0.14 0.00233 NA NA -0.027 0.119 NA
DI
0.18 0.00006 NA NA -0.019 0.124 NA
P
1
is the p-value of naï ve estimate
.
P
2
is the p-value of bias-corrected estimate
.
*%Diff equals to
.
128
Figure 5.1. The association between baseline caloric intakes and rates of change in 2-hr
glucose, S
I
, and DI during the follow-up. Data are presented as means and standard
deviations by lower 50
th
percentile and higher 50
th
percentile of baseline caloric intake
(calories < 2055 kcal/day vs. calories ≥ 2055 kcal/day). β and P are regression
coefficients and p-values for associations between baseline caloric intakes and the rates
of change in 2-hr glucose, S
I
, and DI after adjusting for baseline values and the change in
BMI.
129
Chapter 6
The association between intrauterine GDM exposure and
offspring’s BMI growth from age 5 to 11
6.1 Introduction
Previous studies in Pima Indians have showed that offspring of diabetic women
were more likely to be predisposed to obesity in their later life than offspring of non-
diabetic women (Pettitt et al., 1983; Pettitt et al., 1987; Pettitt et al., 1993). The
mechanisms of such associations were not completely known. There is evidence that the
development in a diabetic intrauterine environment results in additional insulin reactions
by the fetal pancreas due to excess glucose load. Additional insulin can act as a fetal
growth hormone promoting growth and adiposity (Dabelea, 2007).
Gestational diabetes mellitus (GDM) women are hyperglycemic during the
pregnancy without pre-existing diabetes (Association, 2011; Buchanan et al., 2007).
Many studies have observed associations between intrauterine exposure to GDM and
increased offspring birth weight, body mass index (BMI), adiposity, and being
overweight at various ages (Baptiste-Roberts et al., 2012; Boerschmann et al., 2010;
Crume et al., 2011b; Franks et al., 2006; Gillman et al., 2003; Lawlor et al., 2010; Lawlor,
Lichtenstein and Langstrom, 2011; Tsadok et al., 2011; Vaarasmaki et al., 2009).
To our knowledge, only one small study of 504 children investigated the association
between GDM exposure and offspring longitudinal BMI growth trajectory from birth to the
130
age of 13 years (Crume et al., 2011a). Although this study includes non-Hispanic Whites,
African Americans, and Hispanics, it does not assess the association in the Asian population,
which has the highest prevalence of GDM (Berkowitz et al., 1992; Lawrence et al., 2008) .
Thus, there is a need to evaluate such association in a larger and ethnically diverse
population.
In this study, we used a large population based data from a health care plan with a
demographically and socioeconomically diverse population to investigate the association
between intrauterine GDM exposure and offspring BMI growth trajectory over an age
span of 5 to 11 years. For this analysis, eligible participants were GDM offspring and a
random sample of children unexposed to GDM who were matched with GDM offspring
on their age of entering the study, sex, race, birth year, and maternal overweight (yes/no)
defined by BMI higher or lower than 25 kg/m
2
(WHO, 2000).
6.2 Study Design and Methods
6.2.1 Population and Data Sources
This is a retrospective cohort study of children who were born as a singleton at
28-44 weeks gestation in Kaiser Permanente Southern California (KPSC) hospitals
between Jan.1, 1995 and Dec. 31, 2009 and had measured weight and height data from
electronic medical records (EMR) prior to Dec. 31, 2011. Only children > 2-years-old with
valid BMI and birth weight were included (n=274,174). KPSC is a large prepaid group-
practice managed health care organization with approximately 30,000 deliveries per year
and over 3 million members in 2011. Members receive their health care in KPSC-owned
131
facilities throughout the 7-county region. According to institutional standards, all
pregnant women receiving prenatal care at KPSC were requested to take either a 1-hour
50g oral glucose challenge (GCT) test followed by a 3-hour 100g oral glucose tolerance
test (OGTT) if the GCT test results was positive, or a single 2-hour 75g OGTT between
24 and 28 weeks gestation. GDM was identified based on the following hierarchy:1) at
least two 3-hour OGTT values meeting or exceeding the Carpenter and Coustan threshold
values (Association, 2010) (fasting > 95 mg/dl, 1-hr >180 mg/dl, 2-hr >155 mg/dl, and 3-
hr >140 mg/dl); or 2) at least two 2-hour OGTT values greater than or equal to the same
threshold values (Association, 2010); or 3) a single level of 1-hr 50-g GCT > 200 mg/dl.
Women with pre-existing diabetes prior to each pregnancy were excluded (n=6,172).
This study was approved by the KPSC Institutional Review Board.
Child sex, gestational week at birth, and birth weight were obtained from each
infant birth certificate or EMR. Maternal age at delivery, education, race/ethnicity, and
parity were also obtained from the infant birth certificate. Maternal education was
classified into < high school graduate, some college, > college graduate, or missing.
Parity was categorized into 0, 1, or > 2. Maternal race/ethnicity was categorized as
Hispanic (regardless of race), non-Hispanic White (NHW), Black, and Asian/Pacific
Islanders (API). Children with maternal race/ethnicity other than those mentioned above
or multiple race/ethnicity were excluded (n=3,250). We further excluded children with
missing mother’s parity (n=3,167).
Because left truncation on children’s BMI data occurred due to the gradual
implementation of BMI data collection since 2004 and the staged entry of the study
132
cohort, we only selected children who entered the study when they were five years of age
as examples for the association analysis using bias-correction formula. The choice of this
sample is based on the fact that the growth trajectory of BMI over the follow-up years is
linear from age 5 to 11, and the sample size is the largest when compared to other age
ranges with linear BMI growth trajectories.
6.2.2 Data Analysis
BMI was calculated as weight in kilograms divided by height in meters squared.
Gestational week at delivery, birth weight, sex, as well as maternal age, parity, education,
and race/ethnicity at delivery were compared between children exposed to GDM in utero
and those unexposed to GDM in utero. Mixed effects models were used to assess the
offspring baseline BMI association with their BMI growth trajectory over time, and the
association between intrauterine GDM exposure and offspring BMI growth trajectory.
Baseline BMI, race, sex, maternal age, parity, education, gestational week at delivery,
and birth weight were included as covariates to adjust for potential confounding.
Residual analysis was performed to adjust for these confounders. Due to measurement
errors inherent in BMI, the estimates of baseline BMI association with the longitudinal
trajectory of BMI, as well as the estimates of the baseline-BMI adjusted association
between intrauterine GDM exposure and offspring BMI growth trajectory could be
biased. We corrected the potential bias by applying our bias-correction formulae
described in Chapter 2 and Chapter 3 for the association parameter estimates. Variation
due to random measurement error in BMI was estimated from 1246 adults in the
BetaGene study who had repeated measurements of weight and height at one day apart
133
oGTT and ivGTT test (Black et al., 2008). We assumed that the true values of weight and
height did not change during one day period. The z-statistic for bias-corrected association
estimates were calculated using the delta method described in Chapter 5. Because we
expected the estimation of measurement error variations in BMI from the BetaGene study
could be smaller than measurement error variations in clinical measurements of BMI, we
performed an exploratory analysis assuming a larger proportion of the total variation of
BMI was accounted by random measurement error (5-40%). Conditioned on these unreal
measurement error estimates, we assessed the impact of our bias-correction formulae on
the naï ve association estimates for both model I(b) and model II(b).
All statistical tests were two-sided. SAS version 9.3 (SAS Institute Inc., Cary,
North Carolina) was used for data analysis.
6.3 Results
6.3.1 Characteristics of GDM exposed and unexposed offspring at baseline and
during the follow-up
Among 1220 pairs of GDM exposed and unexposed offspring who entered the
study at five years of age, 48.8% were females, and the distribution of race groups is 16.0%
non-Hispanic Whites (NHW), 9.9% African Americans, 57.8% Hispanics, and 16.3%
Asian/Pacific Islanders (API). Table 6.1 compared the cohort characteristics between the
GDM exposed and unexposed children. Overall, maternal characteristics at the time of
delivery differed between GDM exposed and unexposed, where mothers in the GDM
exposed group were on average approximately 3-years older and had higher parity. Children
134
in the GDM exposed group were delivered a half week earlier than children in the
unexposed group. Levels of maternal education were lower among children exposed to
GDM, and children in the GDM exposed group were 24g lower in birth weight. BMI at five
years of age was higher in the GDM exposed group compared to the GDM unexposed group.
The mean follow-up years was 3.6 years (range: 1-6 years). BMI significantly increased
during the follow-up (mean increasing rate of 0.7 kg/m
2
per year, p<0.0001, Figure 6.1).
6.3.2 The association between baseline BMI and the rate of increase in BMI over time
Table 6.2 compared results of baseline BMI association with the longitudinal change
in BMI using Model I(a) and Model II(a). Based on naï ve association estimates, both
models found higher baseline BMI significantly predicted a faster increase of BMI between
the age of 5 to 11 years (p<0.0001, Figure 6.2), and the association estimates from both
models were very close to each other. Because measurement error variation was small in
BMI measures (the estimated proportion of measurement error variation in the entire
BMI variation was 0.7% based on the BetaGene study), bias-correction formulae (2.12)
and (2.23) did not significantly change the naï ve estimates. Even when we presumably
added more measurement error in BMI, bias-corrected association estimates were still
significantly positive and the magnitude of the estimates was enlarged (Table 6.2).
Additionally, Byth’s method provided similar result that the association estimate was
0.119 (p<0.0001).
Because not all subjects had complete follow-up visits over time, there were data
missing at random. Based on our previous evaluations of biases in association estimates
using various mixed effects models, we selected Model I(a) and bias-correction formulae
135
(2.12) for final data analysis to obtain most accurate association estimate. After adjusting
for sex, race, maternal age, parity, education, gestational week at delivery, and birth
weight, the significant association between baseline BMI and the rate of increase in BMI
still remained (bias-corrected estimate =0.092, p<0.0001).
6.3.3 The association between intrauterine GDM exposure and offspring BMI growth
curve adjusting for baseline BMI
Univariate analysis showed children exposed to GDM had higher increase of BMI
between 5 to 11 years of age (regression coefficient=0.125, p<0.0001, Figure 6.3). This
association remained significant after additionally adjusting for race, sex, maternal age,
parity, education, gestational week at delivery, and birth weight (regression
coefficient=0.108, p=0.0002). Because higher baseline BMI was significantly associated
with intrauterine GDM exposure (p=0.0001) and faster increase in BMI over time
(p<0.0001), baseline BMI might be a potential confounder to be adjusted for. Table 6.3
shows the association between offspring GDM exposure status and the rate of increase in
BMI after adjusting for baseline BMI and addition aforementioned confounders. Results
are compared between Model I(b) and Model II(b). Before bias-correction, adjustment
for baseline BMI attenuated the regression coefficient of GDM exposure association with
the rate of increase in BMI by 48.8% using Model I(b) (p=0.13), and by 23.1% using
Model II(b) (p=0.007). Because the measurement error variation was small in BMI
measures, bias-correction formulae (3.13) and (3.22) did not significantly change the
naï ve association estimates after further adjustment for baseline BMI for both Model I(b)
and Model II(b). If the measurement error variation became larger, bias-correction
136
formulae changed naï ve estimates by -0.5% to -47.7% for Model I(b), and by -1.4% to -
143.3% for Model II(b) when the measurement error variation was enlarged to account
for 5% to 40% variation in BMI. The association estimate using Harrison’s method was
closer to the estimate using Model I(b) (
=0.061).
Same as before, we selected Model I(b) and bias-correction formulae (3.13) as the
most correct estimation of the association between intrauterine GDM exposure and the
rate of increase in BMI during the follow-up after adjusting for potential confounders
including baseline BMI. Thus, after adjustment for potential confounders, intrauterine
GDM exposure was not significantly associated with the rate of increase in BMI during
the follow-up (p=0.13). It should also be noted that higher baseline BMI still significantly
predicted faster increase of BMI between the ages of 5 to 11 years even after adjusting
for intrauterine GDM exposure, and additional founders mentioned before (p<0.0001).
6.4 Conclusion
We found that a higher baseline BMI significantly predicted the faster increase of
BMI between the ages of 5 to 11 years. However, intrauterine GDM exposure was not
significantly associated with the rate of increase in BMI after adjusting for baseline BMI,
race, sex, maternal age, parity, education, gestational week at delivery, and birth weight.
This non-significant finding could be because intrauterine GDM exposure might cause
the difference of offspring BMI growth earlier than five years old. Thus the impact of
intrauterine GDM exposure on the BMI growth can be reflected by baseline BMI at five
137
years of age. Adjustment for baseline BMI in this case may mask the association between
intrauterine GDM exposure and the rate of increase of BMI after five years old.
Because measurement error in BMI is very small, biases in naï ve estimates using
Model I(a), Model II(a), Model I(b) and Model II(b) were small, thus our bias-correction
formulae did not obviously change the naï ve association estimates. However, if more
measurement error variation existed in the outcome, bias-correction formulae would alter
the naï ve estimates more significantly.
Although Model I(a), Model II(a) and Byth’s method provided consistent results
for the association between baseline BMI and the rate of increase in BMI over time,
results of the association between intrauterine GDM exposure and the rate of increase in
BMI were different between using Model I(b) and Model II(b). Although Model II(b)
found the significant association, based on our previous simulation results, Model I(b)
produced a more accurate association estimate after bias-correction, especially for the
study with data missing at random. It should be noted that Harrison’s method provided a
similar association estimate as Model I(b) for the association between intrauterine GDM
exposure and the rate of increase in BMI. Thus the result from Model I(b) was used to
generate our final conclusion that no significant association was found between
intrauterine GDM exposure and the rate of increase in BMI after adjusting for baseline
BMI and other confounders.
In conclusion, higher baseline BMI at five years predicts a faster increase of BMI
between the ages of 5 to 11 years. The intrauterine GDM exposure was not associated
with offspring BMI growth trajectories after adjusting for baseline BMI and other
138
potential confounders. Adjustment for baseline BMI significantly influenced the
association between intrauterine GDM exposure and the rate of increase in BMI over
time. Because measurement error in BMI was small, the impact of measurement error on
the bias in the naï ve association estimates was also small, thus there was little need to
correct the bias. However, if the measurement error became larger, our bias-correction
formulae would calibrate the bias in the naï ve estimates more significantly.
139
Table 6.1 Comparison of characteristics of offspring exposed and unexposed to GDM in
utero for children entering the study at five years of age*.
GDM Exposed
GDM
Unexposed
P
Maternal age (yrs) 32.7 (0.08) 29.5 (0.08) <0.0001
Gestational week at delivery
38.7 (0.03) 39.4 (0.02) <0.0001
Birth weight (kg) 3.36 (0.009) 3.38 (0.008) 0.038
Baseline BMI (kg/m
2
)
17.1 (0.09) 16.7 (0.07) 0.0001
Parity
<0.0001
0 32.8 39
1 32 35.8
≥2 35.3 25.2
Mother Education
0.0002
≤ high school graduate 47.5 43.6
some college 26.1 27.8
≥ college graduate 25.7 27.4
Uknown 0.7 1.2
*Means and standard errors are presented for continuous variables, and a pairwise t-test is used to compare the means
of GDM exposed and unexposed offspring. The proportion of subjects in each level is presented for categorical
variables, and chi-squre test is used to test the difference in frequency distribution among GDM exposed and
unexposed offspring.
140
Table 6.2, The association between baseline BMI and the longitudinal change in BMI
between the ages of 5 to 11 years*.
Model I(a)
0.080 0.00004 <0.0001 -0.3 0.0066 0.080 <0.0001
0.05 0.084 <0.0001
0.1 0.088 <0.0001
0.2 0.100 <0.0001
0.3 0.114 <0.0001
0.4 0.133 <0.0001
Model II(a)
0.081 0.00003 <0.0001 -0.3 0.0066 0.084 <0.0001
0.05 0.103 <0.0001
0.1 0.127 <0.0001
0.2 0.183 <0.0001
0.3 0.256 <0.0001
0.4 0.353 <0.0001
is the naï ve estimate of the baseline BMI association with the longitudinal change in BMI using model I(a)
and model II(a).
is the variance of the naive estimate
.
is the p-value for significance test of
the naive estimate
. is the time parameter defined in formula (2.16). is the estimated proportion of
measurement error variance in the entire variance of BMI.
is the bias-corrected estimate of the baseline BMI
association with the longitudinal change in BMI after using bias-correction formulae (2.12) and (2.23) for model I(a)
and model II(a), respectively.
is the significance test of bias-corrected estimate
using the delta method
described in formula (4.1).
Table 6.3, The association between GDM exposure in untero and the rate of increase in BMI between the ages of 5 to 11 after
adjusting for baseline BMI, race, sex, maternal age, parity, education, gestational week at delivery, and birth weight.*
P
GDM× t
P
b_BMI× t
Model I(b)
0.055 0.0013 0.1281 0.078 0.00004 <0.0001 -0.00002 0.44 -0.3 0.0066 0.055 0.130 0.091 <0.0001
0.05 0.053 0.143 0.095 <0.0001
0.1 0.051 0.162 0.100 <0.0001
0.2 0.045 0.212 0.113 <0.0001
0.3 0.038 0.293 0.129 <0.0001
0.4 0.029 0.429 0.151 <0.0001
Model II(b)
0.083 0.0010 0.007 0.082 0.00003 <0.0001 -0.00001 0.44 -0.3 0.0066 0.082 0.008 0.085 <0.0001
0.05 0.074 0.017 0.103 <0.0001
0.1 0.063 0.041 0.127 <0.0001
0.2 0.038 0.214 0.184 <0.0001
0.3 0.007 0.833 0.257 <0.0001
0.4 -0.036 0.246 0.354 <0.0001
is the naï ve estimate of the association between intrauterine GDM exposure and offspring rates of increase in BMI after adjusting for baseline BMI using
model I(a) and model II(a).
is the variance of the naï ve estimate
. P
GDM× t
is the p-value of significance test of the naï ve estimate
.
is the
naï ve estimate of the baseline BMI association with the longitudinal change in BMI after adjusting for intrauterine GDM exposure using model I(a) and model II(a).
is the variance of the naive estimate
.
is the p-value for significance test of the naive estimate
.
is the covariance
between naï ve estimates
and
. is the linear regression coefficient of the association between intrauterine GDM exposure and baseline BMI. is the
time parameter defined in formula (2.16). is the estimated proportion of measurement error variance in the entire variance of BMI.
is the bias-corrected
estimate using correction formulae (3.13) and (3.22) for model I(a) and model II(a), respectively.
is the p-value for significance test of bias-corrected estimate
using the delta method as formula (4.3).
is the bias-corrected estimate of the baseline BMI association with the rate of increase of BMI after using bias-
correction formulae (3.12) and (3.21) for model I(a) and model II(a), respectively.
is the significance test of bias-corrected estimate
using the delta
method described in formula (4.2).
141
142
Figure 6.1 Offspring BMI growth between the ages of 5 to 11. Means and standard errors
of BMI at each visit are presented. n is the number of sample size at each visit.
BMI (kg/m
2
)
Offspring Age (years)
143
Figure 6.2 Offspring BMI growth curves between the ages of 5 to 11 stratified by median
baseline BMI (16.2kg/m
2
). Means and standard errors of BMI at each visit are presented.
n is the number of sample size at each visit.
144
Figure 6.3 Offspring BMI growth curves between the ages of 5 to 11 stratified by
children exposed and unexposed to GMD in utero. Means and standard errors of BMI at
each visit are presented. n is the number of sample size at each visit.
BMI (kg/m
2
)
Offspring Age (years)
GDM unexposed GDM exposed
145
Chapter 7
Summary and future directions
7.1 Overview
We reviewed various mixed effects models for analyzing the associations between
baseline and exposure variables with the longitudinal trajectory when responses are
subject to measurement error. Results from the simulation study showed that naï ve
association estimates were biased if the outcome variable was prone to random
measurement error. We proposed new methods to correct such biases and validated our
method using the simulation study. Two real studies were analyzed as examples for
applying our bias-correction formulae.
When the measurement error in the outcome variable cannot be ignored,
correcting the bias in the naï ve association estimate is necessary to obtain a correct
estimate of the true baseline association with the longitudinal change during the follow-
up. For the baseline-adjusted association between the exposure variable and the
longitudinal change, measurement error in the outcome variable can cause the bias in the
naive association estimate when the exposure variable is correlated with the baseline
value of the outcome.
Although some researchers have performed analysis using Model II(a) and Model
II(b) with the baseline value of the outcome included as both the dependent and the
independent variables (Gardiner et al., 2012; Josephs et al., 2011), we show that these are
146
incorrect models for data analysis which can cause biased association estimates and
inflated type I error. Although Byth’s and Harrison’s methods do not require additional
data to estimate the variance of measurement error, these methods fail to specify the
residual variability due to measurement error or other reasons such as true biological
variations and model misspecifications. In real studies, variations other than measurement
error always exist in the data, thus Byth’s and Harrison’s methods still provide biased
association estimates.
The advantage of our method is that we show the close forms to correct biases in
naï ve estimates due to measurement error in the outcome for Model I(a) and Model I(b).
These bias-correction formulae are easy to apply and can successfully correct the biases
in the naï ve estimates. Also, we proposed two methods (bootstrapping method and delta
method) for testing the significance of bias-corrected estimates. The bootstrapping
method keeps the type I error rate better than the delta method, especially for a small
sample size when the outcome variable is prone to measurement error. However, it is
computationally intense for the large data. The delta method is easier to apply. However,
some assumptions need to be made for the test statistics and the type I error rate is
inflated when the sample size is small. However, the delta method controls the type I
error rate well when the sample size is relatively large (more than 100) if using the mixed
effects model excluding the baseline value from the outcome. It should be noted that an
inflated type I error exists in a small sample size even when no measurement error exists
in the outcome variable. Some previous studies also suggest that the mixed effects model
may result in an inflated type I error rate when the sample size is small or covariance
147
matrices are unequal (Keselman, Algina and Kowalchuk, 2001; Keselman et al., 1999a,
b).
Along with our findings, we also acknowledge some limitations of our methods,
which will be discussed in the following sections. These limitations need to be studied in
the future research.
7.2 Correcting measurement error with no additional data
One limitation of our method is that the additional data is required to estimate
measurement error variation in the real data analysis. However, the additional data is not
always available and the estimate of measurement error variation from the additional data
might not be accurate for the analysis data. Methods not requiring the additional data to
estimate measurement error variation should be studied. Byth’s method and Harrison’s
method have the advantage to address this problem. Their methods use the residual
variation to estimate the measurement error variation in the outcome. However, both
methods are strictly conditioned on the linear growth trajectory. Any violation of the
linear growth trajectory of the underlying true data can cause biases in the association
estimate.
A possible way to improve Byth’s method and Harrison’s method is to relax
the assumption of the linear growth trajectory by using the non-parametric mixed effects
model to fit data more accurately. Because the main problem of Byth’s method and
Harrison’s method is that measurement error variations are overestimated by residual
variations when the longitudinal trajectory of the data is not perfectly linear. If a non-
148
linear model can be used to fit the data, the residual variations will be more likely to only
arise from measurement error. Rice et al. proposed incorporating a B-spline function (Lee,
1982) into the mixed effects model and to model the individual curve by spline functions
with random coefficients (Rice and Wu, 2001). Rice’s model has been used in the
measurement error adjustment of independent variables when no additional data is
available to estimate measurement error variations. In the study of treatment effect on
CD4+ cell counts over time, Liang et al. used Rice’s method to model time-varied
covariates measured with error, and then used regression a calibration approach to correct
biases in the association estimate (Liang, Wu and Carroll, 2003).
Under a linear regression framework, a two-stage least square (TSLS) estimation
has been used to obtain an unbiased association estimate for independent variables
measured with error (Lewbel, 1997). TSLS is generally used to address the problem
when errors in the dependent variable are correlated with the independent variable, which
can cause biased estimates using OLS estimation (Hoffman, 1987). Biased estimates
appear because when measurement error exists in either dependent or independent
variable, the model residual could be correlated with the independent variable. The
procedure of TSLS analysis is to predict the error-prone variable using some instrumental
variables that are uncorrelated with the measurement error and the residual error, but
correlated with the underlying true value of the variable (the first stage), and then uses
those predicted values to estimate a linear regression model of the dependent variable
(the second stage). Several instrumental variables are proposed in Lewbel’s paper to
address the measurement error in the independent variable. Because in our derivations of
149
bias-correction formulae, we separate the mixed effects model into a two-level model,
and then incorporate an estimated measurement error variance into the OLS estimation to
correct the bias in the naï ve estimate, TSLS estimation might be an approach to optimize
the association estimation without addition data to estimate measurement error variation.
7.3 Analysis of non-linear growth trajectories
Although our method is developed based on models using a linear time effect,
non-linear time effects such as quadratic, log, or another power-transformed time effect
can also be applied using our method. Bias-correction formulae for Model I(a) and
Model I(b) are the same. But for Model II(a) and Model II(b), time-related parameter
should be calculated based on power-transformed time variables. However, our bias-
correction formulae cannot be directly applied to the situation when both linear and
higher order polynomial time effects are treated as random effects in the model. Although
bias-correction formulae can be developed using our proposed method including two-
level modeling and least square estimates for the association parameter, the derivation
becomes more complicated when more polynomial time effects are included in the model.
Fortunately, quadratic and cubic time effects are usually enough to describe most of the
non-linear growth trajectory in the real study.
For more complex growth trajectories, a parametric model might not be enough to
fit the data accurately, thus non-parametric mixed effects models including smoothing
functions in both fixed effects and random effects are proposed to model both population-
average and individual-specific curves with the same smoothness property. Such models
150
are known as functional mixed effects models (Chen and Wang, 2011; Guo, 2002; Liu
and Guo, 2012; Rice and Wu, 2001). It will be a challenge to extend our method into
functional mixed effects models because smoothing functions usually includes multiple
time knots, which causes the derivation of correction formulae to be very complex using
our method. However, as we mentioned in 7.1, Byth’s method and Harrison’s method
may be able to extend in the field of non-parametric mixed effects model.
7.4 Adjustment for additional confounders
In real applications, there might be other confounders beyond baseline
measurement to be adjusted for in the association between the exposure of interest and
the longitudinal change of the outcome during the follow-up. In this paper, residual
analysis is performed for additional covariates adjustment, which includes two-stage
analysis. First, marginal residuals (146) are taken out from the mixed effects model
regressing repeated measurements on additional covariates. Second, residuals are further
regressed on baseline values and the exposure variable as Model I(b). However, residual
analysis is not optimal for our problem because correlations among the exposure variable,
baseline values, and additional confounders can influence the bias in naï ve estimates for
the association between the exposure variable and the growth trajectories of the outcome.
Thus, bias-correction formulae need to be developed under a multiple regression
framework with a baseline value and other confounders all in one mixed effects model.
It is not difficult to extend our bias-correction method to address this
aforementioned problem as long as these additional confounders are all treated as fixed
151
effects. We can still use the idea of two-level models and model the first level as (3.6),
and build the second level model as
, where
, and total number of additional
confounders. Let
,
as naï ve association
estimates, and
are the corresponding true association parameters.
Based on biases in naï ve estimates using multiple linear regressions (John, 2010),
(7.1), where
is the covariance
matrix of ;
is the covariance matrix of measurement error in independent variables.
For our problem of adjustment for baseline value,
=
2
(k)
(k)
(k)
(k) (k) (k) (k)
0 0 0
0 0 0 0
0 0 0 0
0000
u
and
is
the measurement error variance existing in the observed baseline value
.
is the
covariance matrix of the measurement error in independent variables and the
measurement error in the dependent variable, which is equal to
for
Model I(b), and is equal to
for Model II(b) ( is described in (2.16)
and is a parameter of follow-up time). Although the derivation idea is simple, the
computation can be very complex when more confounders are included in the model.
The above extension of our method is especially valuable that it also addresses the
issue of analyzing categorical variables. Although in our simulation study, we apply our
bias-correction formulae to the naï ve estimate of categorical exposures, we only analyze
152
the categorical exposure as an ordinal variable. It is well known that categorical variables
are not always meaningful to interpret as ordinal variables. When there is no trend across
different levels of a categorical variable, dummy coding is usually used to represent
specific level of the variable. The number of dummy variables is equal to the number of
levels in the categorical variable minus one. In order to include all the dummy variables
in one model, we need to use the above extension to derive bias-correction formulae for
the association estimate of each dummy variable.
7.5 Exposure variables also contain measurement error
In our work, we only correct the bias in naï ve estimates due to the measurement
error in the outcome variable. It is often that exposure variables are also prone to
measurement error. As in the GDM cohort study, caloric intake estimated from self-
reported questionnaires is well-known to carry measurement error (Stram et al., 2000).
Thus, correcting bias in the naï ve estimate due to measurement error in both the outcome
and the exposure variable is important. Based on the discussion in 7.4, we can address the
problem easily for measurement error in the covariates. The relationship between naï ve
estimates (
) and true parameters (
) can still be expressed as described in 7.4.
Additional measurement error in the exposure variable can be added into covariance
matrices of measurement error
and
. Bias-correction formulae can be developed
by solving the equation (7.1).
153
7.6 Conclusion
In the analysis of baseline prediction of the longitudinal trajectory, or the
association between the exposure variable and the longitudinal trajectory adjusting for the
baseline, random measurement error and biological fluctuations in the outcome variable
can cause biases in naï ve association estimates using traditional mixed effects models.
We developed novel methods to correct such biases under the framework of linear mixed
effects models. Our bias-correction formulae were validated in the simulation study and
were easily applied to real studies due to their simple close forms. Our results highlighted
the importance of correcting biases due to measurement error in the outcome to obtain
correct estimates for this type of association analysis. Compared to Byth’s and Harrison’s
methods, our method is more robust to additional variations in data besides the random
measurement error. Additionally, our results indicate that the baseline value should not be
included in the dependent variable when using our proposed mixed effects models to
analyze the data with more than two visits.
We acknowledge that our methods endure some limitations, and require
additional data to estimate measurement error variations and they cannot be easily
applied to non-linear growth trajectories. Some future directions have been proposed to
extend our method and methods developed by Byth et al. (Byth and Cox, 2005) and
Harrison et al. (Harrison et al., 2009) to satisfy future requirements of measurement error
adjustment without additional data, analyzing non-linear growth trajectories, correcting
biases when multiple confounders are adjusted for, and when exposure variables are also
154
subject to measurement errors. These extensions will broaden the application of bias-
correction methods for more general longitudinal data analyses.
155
Bibliography
Mixed Model Influence Diagnostics. Oliver Schabenberger, SAS Institute Inc., Cary, NC,
29.
United States Department of Health and Human Services. 2008 Physical Activity
Guidelines for Americans. Available from
http://www.health.gov/PAGuidelines/guidelines/default.aspx. Accessed May 1,
2011. .
Altman, D. (1982). Statistics in medical journals. Statistics in Medicine 1.
Altman, D. (1991a). Practical Stattistics for Medical Research: Chapman & Hall/CRC:
London, Boca Raton, FL.
Altman, D. G. (1991b). Statistics in medical journals: Developments in the 1980s.
Statistics in Medicine 10, 1897-1913.
Anderson, B. (1990). Methodological errors in medical research.: London: Blackwell.
Archie, J. (1981). Mathematic Coupling of Data. Annals of Surgery 193, 296-303.
Association, A. D. (2010). Diagnosis and classification of diabetes mellitus. Diabetes
care 33 Suppl 1, S62-69.
Association, A. D. (2011). Standards of Medical Care in Diabetes—2011. Diabetes Care
34, S11-S61.
Baird, D. C. (1995). Experimentation: An Introduction to Measurement Theory and
Experiment Design, 3rd edn.: Prentice Hall: Englewood Cliffs, NJ.
Baptiste-Roberts, K., Nicholson, W. K., Wang, N. Y., and Brancati, F. L. (2012).
Gestational diabetes and subsequent growth patterns of offspring: the National
Collaborative Perinatal Project. Maternal and child health journal 16, 125-132.
Barnett, A. G., van der Pols, J. C., and Dobson, A. J. (2005). Regression to the mean:
what it is and how to deal with it. International Journal of Epidemiology 34, 215-
220.
Bergman, R. N. (1989). Lilly lecture 1989. Toward physiological understanding of
glucose tolerance. Minimal-model approach. Diabetes 38, 1512-1527.
Berkowitz, G. S., Lapinski, R. H., Wein, R., and Lee, D. (1992). Race/ethnicity and other
risk factors for gestational diabetes. Am J Epidemiol 135, 965-973.
Bevington, P. R., and Robinson, D. K. (1992). Data Reduction and Error Analysis for the
Physical Sciences, 2nd edn.: McGraw-Hill: New York.
Black, M. H., Fingerlin, T. E., Allayee, H., et al. (2008). Evidence of Interaction Between
PPARG2 and HNF4A Contributing to Variation in Insulin Sensitivity in Mexican
Americans. Diabetes 57, 1048-1056.
Bland, J. M., and Altman, D. G. (1996). Measurement error. BMJ 312.
Blomqvist, N. (1977). On the Relation Between Change and Initial Value. Journal of the
American Statistical Association 72, 746-749.
Blomqvist, N. (1987). On the bias caused by regression toward the mean in studying the
relation between change and initial value. Journal of Clinical Periodontology 14,
34-37.
156
Blomqvist, N., Cederblad, G., Korsan-Bengtsen, K., and Wallerstedt, S. (1977).
Application of a method for correcting an observed regression between change
and initial value for the bias caused by random errors in the initial value. Clin
Chem 23, 1845-1848.
Boerschmann, H., Pfluger, M., Henneberger, L., Ziegler, A. G., and Hummel, S. (2010).
Prevalence and predictors of overweight and insulin resistance in offspring of
mothers with gestational diabetes mellitus. Diabetes care 33, 1845-1849.
Box, G. E. P. (1950). Problems in the Analysis of Growth and Wear Curves. Biometrics 6,
362-389.
Buchanan, T. A., Xiang, A., Kjos, S. L., et al. (1998). Gestational diabetes: antepartum
characteristics that predict postpartum glucose intolerance and type 2 diabetes in
Latino women. Diabetes 47, 1302-1310.
Buchanan, T. A., Xiang, A., Kjos, S. L., and Watanabe, R. (2007). What Is Gestational
Diabetes? Diabetes Care 30, S105-S111.
Buchanan, T. A., Xiang, A. H., Peters, R. K., et al. (2002). Preservation of Pancreatic β-
Cell Function and Prevention of Type 2 Diabetes by Pharmacological Treatment
of Insulin Resistance in High-Risk Hispanic Women. Diabetes 51, 2796-2803.
Burton, P., Gurrin, L., and Sly, P. (1998). Extending the simple linear regression model
to account for correlated responses: An introduction to generalized estimating
equations and multi-level mixed modelling. Statistics in Medicine 17, 1261-1291.
Byth, K., and Cox, D. R. (2005). On the relation between initial value and slope.
Biostatistics 6, 395-403.
Carroll, R., Ruppert, D., Stefanski, L., and Crainiceanu, C. (2006). Measurement Error in
Nonlinear Models: Chapman & Hall, Boca Raton.
Carroll, R. J., Gallo, P., and Leon Jay, G. (1985). Comparison of Least Squares and
Errors-in-Variables Regression, With Special Reference to Randomized Analysis
of Covariance. Journal of the American Statistical Association 80, 929-932.
Carroll, R. J., Senn, S., and Carroll, R. J. (1990). Covariance analysis in generalized
linear measurement error models. Statistics in Medicine 9, 583-586.
Casella, G., and Berger, R. L. (2002). Statistical Inference, 2nd Edition: Thomson
Learning.
Chambless, L. E., and Davis, V. (2003). Analysis of associations with change in a
multivariate outcome variable when baseline is subject to measurement error.
Statistics in Medicine 22, 1041-1067.
Chan, L., and Mak, T. (1979). Maximum likelihood estimation of a linear structural
relationship with replication. J.R. Statist. Soc. Ser. B 41, 263-268.
Chan, S. F., Macaskill, P., Irwig, L., and Walter, S. D. (2004). Adjustment for baseline
measurement error in randomized controlled trials induces bias. Controlled
Clinical Trials 25, 408-416.
Chen, H., and Wang, Y. (2011). A penalized spline approach to functional mixed effects
model analysis. Biometrics 67, 861-870.
Cnop, M., Vidal, J., Hull, R. L., et al. (2007). Progressive Loss of β-Cell Function Leads
to Worsening Glucose Tolerance in First-Degree Relatives of Subjects With Type
2 Diabetes. Diabetes Care 30, 677-682.
157
Cole, J. W. L., and Grizzle, J. E. (1966). Applications of Multivariate Analysis of
Variance to Repeated Measurements Experiments. Biometrics 22, 810-828.
Cook, J. R., and Stefanski, L. A. (1994). Simulation-Extrapolation Estimation in
Parametric Measurement Error Models. Journal of the American Statistical
Association 89, 1314-1328.
Corbeil, R. R., and Searle, S. R. (1976). Restricted Maximum Likelihood (REML)
Estimation of Variance Components in the Mixed Model. Technometrics 18, 31-
38.
Crume, T. L., Ogden, L., Daniels, S., Hamman, R. F., Norris, J. M., and Dabelea, D.
(2011a). The impact of in utero exposure to diabetes on childhood body mass
index growth trajectories: the EPOCH study. The Journal of pediatrics 158, 941-
946.
Crume, T. L., Ogden, L., West, N. A., et al. (2011b). Association of exposure to diabetes
in utero with adiposity and fat distribution in a multiethnic population of youth:
the Exploring Perinatal Outcomes among Children (EPOCH) Study. Diabetologia
54, 87-92.
Dabelea, D. (2007). The Predisposition to Obesity and Diabetes in Offspring of Diabetic
Mothers. Diabetes Care 30, S169-S174.
Despres, J. P., Pouliot, M. C., Moorjani, S., et al. (1991). Loss of abdominal fat and
metabolic response to exercise training in obese women. American Journal of
Physiology - Endocrinology And Metabolism 261, E159-E167.
Diagnosis, T. E. C. o. t., and Mellitus*, C. o. D. (2003). Follow-up Report on the
Diagnosis of Diabetes Mellitus. Diabetes Care 26, 3160-3167.
Dominici, F., Zeger, S. L., and Samet, J. M. (2000). A measurement error model for time-
series studies of air pollution and mortality. Biostatistics 1, 157-175.
Efron, B. (1982). The Jackknife, the Bootstrap and Other Resampling Plans: SIAM,
Philadelphia.
Efron, B., and Tibshirani, R. (1993). An Introduction to the Bootstrap. Boca Raton, FL:
Chapman & Hall/CRC.
Espino-Hernandez, G., Gustafson, P., and Burstyn, I. (2011). Bayesian adjustment for
measurement error in continuous exposures in an individually matched case-
control study. BMC Medical Research Methodology 11, 67.
Everitt, B. S. (1995). The Analysis of Repeated Measures: A Practical Review with
Examples. Journal of the Royal Statistical Society. Series D (The Statistician) 44,
113-135.
Festa, A., Williams, K., D’Agostino, R., Wagenknecht, L. E., and Haffner, S. M. (2006).
The Natural Course of β-Cell Function in Nondiabetic and Diabetic Individuals.
Diabetes 55, 1114-1120.
Fisher, R. A. (1918). The correlation between relatives on the supposition of Mendelian
inheritance. Transactions of the Royal Society of Edinburgh 52, 399-433.
Fisher, R. A. (1925). Statistical Methods for Research Workers: Oliver and Boyd,
Edinburgh.
Fluke. (1994). Calibration: Philosophy and Practice, 2nd. edn.: Fluke Corporation:
Everett, WA.
158
Fox, J. P., and Glas, C. A. (2003). Bayesian Modeling of Measurement Error in Predictor
Variables Using Item Response Theory. Psychometrika 68, 169-191.
Franks, P. W., Looker, H. C., Kobes, S., et al. (2006). Gestational glucose tolerance and
risk of type 2 diabetes in young Pima Indian offspring. Diabetes 55, 460-465.
Fraser, G. E., and Stram, D. O. (2001). Regression calibration in studies with correlated
variables measured with error. Am J Epidemiol 154, 836-844.
Fuller, W. A. (1987). Measurement Error Models: John Wiley & Sons.
Gardiner, J. C., Luo, Z., and Roman, L. A. (2009). Fixed effects, random effects and GEE:
What are the differences? Statistics in Medicine 28, 221-239.
Gardiner, S. K., Johnson, C. A., and Demirel, S. (2012). Factors Predicting the Rate of
Functional Progression in Early and Suspected Glaucoma. Investigative
Ophthalmology & Visual Science 53, 3598-3604.
Garrett, M. F., Nan, M. L., and James, H. W. (2004). Applied Longitudinal Analysis: John
Wiley & Sons, Inc., Hoboken, New Jersey.
Geisser, S. (1963). Multivariate Analysis of Variance for a Special Covariance Case.
Journal of the American Statistical Association 58, 660-669.
Gill, J. S., Beevers, D. G., Zezulka, A. V., and Davies, P. (1985). RELATION
BETWEEN INITIAL BLOOD PRESSURE AND ITS FALL WITH
TREATMENT. The Lancet 325, 567-569.
Gillman, M. W., Rifas-Shiman, S., Berkey, C. S., Field, A. E., and Colditz, G. A. (2003).
Maternal gestational diabetes, birth weight, and adolescent obesity. Pediatrics 111,
e221-226.
Green, R., and van de Vijver, F. (1993). A simple test of the law of initial value.
Psychophysiology 30, 525-530.
Greenhouse, S. W., and Geisser, S. (1959). On methods in the analysis of profile data.
Psychometrika 24, 95-112.
Guilford, J., and Fruchter, B. (1973). Fundamental statistics in psychology and education:
New York: McGraw-Hill.
Guo, W. (2002). Functional Mixed Effects Models. Biometrics 58, 121-128.
Gustafson, P. (2004). Measurement Error and Misclassification in Statistics and
Epidemiology: Chapman & Hall, Boca Raton.
Halverson, J., and RE., K. (1981). Gastric bypass: analysis of weight loss and factors
determining success. Surgery 90, 446-455.
Hand, D. J., and Taylor, C. C. (1987). Multivariate Analysis of Variance and Repeated
Measures: London: Chapman and Hall.
Harrison, L., Dunn, D. T., Green, H., and Copas, A. J. (2009). Modelling the association
between patient characteristics and the change over time in a disease measure
using observational cohort data. Statistics in Medicine 28, 3260-3275.
Harville, D. A. (1988). Mixed-Model Methodology: Theoretical Justifications and Future
Directions. In Proceedings of the Statistical Computing Section, 41-49: American
Statistical Association, New Orleans.
Harville, D. A. (1990). BLUP (Best Linear Unbiased Prediction), and Beyond. In
Advances in Statistical Methods for Genetic Improvement of Livestock, 239-276:
Springer-Verlag.
159
Hayes, R. J. (1988). Methods for assessing whether change depends on initial value.
Statistics in Medicine 7, 915-927.
Henderson, C. R. (1975). Best Linear Unbiased Estimation and Prediction under a
Selection Model. Biometrics 31, 423-447.
Henderson, C. R. (1984). Application of Linear Models in Animal Breeding: University
of Guelph.
Henderson, C. R., Kempthorne, O., Searle, S. R., and Krosigk, C. M. v. (1959). The
Estimation of Environmental and Genetic Trends from Records Subject to Culling.
Biometrics 15, 192-218.
HENRY, R. R., SCHEAFFER, L., and OLEFSKY, J. M. (1985). Glycemic Effects of
Intensive Caloric Restriction and Isocaloric Refeeding in Noninsulin-Dependent
Diabetes Mellitus. Journal of Clinical Endocrinology & Metabolism 61, 917-925.
Hoffman, D. L. (1987). Two-Step Generalized Least Squares Estimators in Multi-
Equation Generated Regressor Models. The Review of Economics and Statistics
69, 336-346.
Ironson, G., O’Cleirigh, C., Fletcher, M. A., et al. (2005). Psychosocial Factors Predict
CD4 and Viral Load Change in Men and Women With Human Immunodeficiency
Virus in the Era of Highly Active Antiretroviral Treatment. Psychosomatic
Medicine 67, 1013-1021.
ISO (1993). Guide to the Expression of Uncertainty in Measurement. In International
Organization for Standardization (ISO) and the International Committee on
Weights and Measures (CIPM): Switzerland.
Jin, P. (1992). Toward a Reconceptualization of the Law of Initial Value. Psychological
Bulletin 111, 176-184.
John, P. (2010). Measurement Error Models, Methods, and Applications: Chapman &
Hall / CRC.
Josephs, K. A., Whitwell, J. L., Weigand, S. D., et al. (2011). Predicting functional
decline in behavioural variant frontotemporal dementia. Brain 134, 432-448.
Kelley, D. E., Wing, R., Buonocore, C., Sturis, J., Polonsky, K., and Fitzsimmons, M.
(1993). Relative effects of calorie restriction and weight loss in noninsulin-
dependent diabetes mellitus. Journal of Clinical Endocrinology & Metabolism 77,
1287-1293.
Kenward, M. G., and Roger, J. H. (2010). The use of baseline covariates in crossover
studies. Biostatistics 11, 1-17.
Keselman, H. J., Algina, J., and Kowalchuk, R. K. (2001). The analysis of repeated
measures designs: a review. Br J Math Stat Psychol 54, 1-20.
Keselman, H. J., Algina, J., Kowalchuk, R. K., and Wolfinger, R. D. (1999a). The
analysis of repeated measurements: a comparison of mixed-model satterthwaite f
tests and a nonpooled adjusted degrees of freedom multivariate test.
Communications in Statistics - Theory and Methods 28, 2967-2999.
Keselman, H. J., Algina, J., Kowalchuk, R. K., and Wolfinger, R. D. (1999b). A
comparison of recent approaches to the analysis of repeated measurements.
British Journal of Mathematical and Statistical Psychology 52, 63-78.
Kirkwood, B., and JAC., S. (2003). Medical Statistics (2nd edn): Blackwell: Oxford.
160
Knowler WC, Barrett-Connor E, Fowler SE, et al. (2002). Reduction in the Incidence of
Type 2 Diabetes with Lifestyle Intervention or Metformin. New England Journal
of Medicine 346, 393-403.
Kolonel, L. N., Henderson, B. E., Hankin, J. H., et al. (2000). A Multiethnic Cohort in
Hawaii and Los Angeles: Baseline Characteristics. American Journal of
Epidemiology 151, 346-357.
Konstan, M. W., Schluchter, M. D., Xue, W., and Davis, P. B. (2007). Clinical Use of
Ibuprofen Is Associated with Slower FEV1 Decline in Children with Cystic
Fibrosis. American Journal of Respiratory and Critical Care Medicine 176, 1084-
1089.
Kotler, D., Burastero, S., Wang, J., and Pierson, R. (1996). Prediction of body cell mass,
fat-free mass, and total body water with bioelectrical impedance analysis: effects
of race, sex, and disease. The American Journal of Clinical Nutrition 64, 489S-
497S.
Laird, N. M., and Fong, W. (1990). Estimating rates of change in randomized clinical
trials. Controlled Clinical Trials 11, 405-419.
Laird, N. M., and Ware, J. H. (1982). Random-Effects Models for Longitudinal Data.
Biometrics 38, 963-974.
Larson-Meyer, D. E., Heilbronn, L. K., Redman, L. M., et al. (2006). Effect of Calorie
Restriction With or Without Exercise on Insulin Sensitivity, β-Cell Function, Fat
Cell Size, and Ectopic Lipid in Overweight Subjects. Diabetes Care 29, 1337-
1344.
Lawlor, D. A., Fraser, A., Lindsay, R. S., et al. (2010). Association of existing diabetes,
gestational diabetes and glycosuria in pregnancy with macrosomia and offspring
body mass index, waist and fat mass in later childhood: findings from a
prospective pregnancy cohort. Diabetologia 53, 89-97.
Lawlor, D. A., Lichtenstein, P., and Langstrom, N. (2011). Association of maternal
diabetes mellitus in pregnancy with offspring adiposity into early adulthood:
sibling study in a prospective cohort of 280,866 men from 248,293 families.
Circulation 123, 258-265.
Lawrence, J. M., Contreras, R., Chen, W., and Sacks, D. A. (2008). Trends in the
prevalence of preexisting diabetes and gestational diabetes mellitus among a
racially/ethnically diverse population of pregnant women, 1999-2005. Diabetes
Care 31, 899-904.
Lee, E. T. Y. (1982). A simplified B-spline computation routine. Computing 29, 365-371.
Lewbel, A. (1997). Constructing Instruments for Regressions With Measurement Error
When no Additional Data are Available, with An Application to Patents and R&D.
Econometrica 65, 1201-1213.
Liang, H., Wu, H., and Carroll, R. J. (2003). The relationship between virologic and
immunologic responses in AIDS clinical research using mixed-effects varying-
coefficient models with measurement error. Biostatistics 4, 297-312.
Lim, E., Hollingsworth, K., Aribisala, B., Chen, M., Mathers, J., and Taylor, R. (2011).
Reversal of type 2 diabetes: normalisation of beta cell function in association with
decreased pancreas and liver triacylglycerol. Diabetologia 54, 2506-2514.
161
Lindströ m, J., Louheranta, A., Mannelin, M., et al. (2003). The Finnish Diabetes
Prevention Study (DPS). Diabetes Care 26, 3230-3236.
Liu, Z., and Guo, W. (2012). Functional mixed effects models. Wiley Interdisciplinary
Reviews: Computational Statistics 4, 527-534.
MacGregor, G. A., Sagnella, G. A., and K.D., M. (1985). (Letter). Lance i, 926-927.
Matthews, J. N., Altman, D. G., Campbell, M. J., and Royston, P. (1990). Analysis of
serial measurements in medical research. British Medical Journal 300, 230-235.
McLean, R. A., Sanders, W. L., and Stroup, W. W. (1991). A Unified Approach to Mixed
Linear Models. The American Statistician 45, 54-64.
Moreno, L. F., Stratton, H. H., Newell, J. C., and Feustel, P. J. (1986). Mathematical
coupling of data: correction of a common error for linear calculations. Journal of
Applied Physiology 60, 335-343.
Nakamura, T. (1990). Corrected score function for errors-in-variables models:
Methodology and application to generalized linear models. Biometrika 77, 127-
137.
Nö thlings, U., Wilkens, L., Murphy, S., Hankin, J., Henderson, B., and Kolonel, L.
(2007). Body mass index and physical activity as risk factors for pancreatic
cancer: the Multiethnic Cohort Study. Cancer Causes and Control 18, 165-175.
Oldham, P. D. (1962). A note on the analysis of repeated measurements of the same
subjects. Journal of Chronic Diseases 15, 969-977.
Paterno, E. M., and Amemiya, Y. (1996). Random effect and random coefficient analysis
with errors-in-variables. ASA Proceeding of Business and Economic Statistics, 76-
79.
Pettitt, D. J., Baird, H. R., Aleck, K. A., Bennett, P. H., and Knowler, W. C. (1983).
Excessive Obesity in Offspring of Pima Indian Women with Diabetes during
Pregnancy. New England Journal of Medicine 308, 242-245.
Pettitt, D. J., Knowler, W. C., Bennett, P. H., Aleck, K. A., and Baird, H. R. (1987).
Obesity in Offspring of Diabetic Pima Indian Women Despite Normal Birth
Weight. Diabetes Care 10, 76-80.
Pettitt, D. J., Nelson, R. G., Saad, M. F., Bennett, P. H., and Knowler, W. C. (1993).
Diabetes and Obesity in the Offspring of Pima Indian Women With Diabetes
During Pregnancy. Diabetes Care 16, 310-314.
Rice, J. A., and Wu, C. O. (2001). Nonparametric mixed effects models for unequally
sampled noisy curves. Biometrics 57, 253-259.
Robinson, G. K. (1991). That BLUP is a Good Thing: The Estimation of Random Effects.
Statistical Science 6, 15-32.
Rowell, J. G., and Walters, D. E. (1976). Analysing data with repeated observations on
each experimental unit The Journal of Agricultural Science 87, 423-432.
Searle, S. R. (1971). Linear Models: New York: John Wiley & Sons.
Sen, A., Jen, K.-L. C., and Djuric, Z. (2007). Baseline Leptin Levels Predict Change in
Leptin Levels During Weight Loss in Obese Breast Cancer Survivors. The Breast
Journal 13, 180-186.
Senn, S. (1990). Covariance analysis in generalized linear measurement error models.
Stat Med 9, 583-585.
162
Stefanski, L. A., and Cook, J. R. (1995). Simulation-Extrapolation: The Measurement
Error Jackknife. Journal of the American Statistical Association 90, 1247-1256.
Stram, D. O., Hankin, J. H., Wilkens, L. R., et al. (2000). Calibration of the Dietary
Questionnaire for a Multiethnic Cohort in Hawaii and Los Angeles. American
Journal of Epidemiology 151, 358-370.
Stratton, H. H., Feustel, P. J., and Newell, J. C. (1987). Regression of calculated variables
in the presence of shared measurement error. Journal of Applied Physiology 62,
2083-2093.
Su, Z. (2011). Balancing multiple baseline characteristics in randomized clinical trials.
Contemporary Clinical Trials 32, 547-550.
Taylor, J. (1997). An introduction to Error Analysis, 2nd. edn.: University Science Books:
Sausalito, CA.
Tosteson, T. D., Buonaccorsi, J. P., and Demidenko, E. (1998). Covariate measurement
error and the estimation of random effect parameters in a mixed model for
longitudinal data. Statistics in Medicine 17, 1959-1971.
Tsadok, M. A., Friedlander, Y., Paltiel, O., et al. (2011). Obesity and blood pressure in
17-year-old offspring of mothers with gestational diabetes: insights from the
Jerusalem Perinatal Study. Experimental diabetes research 2011, 906154.
Tu, Y.-K., and Gilthorpe, M. S. (2007). Revisiting the relation between change and initial
value: a review and evaluation. Statistics in Medicine 26, 443-457.
Tu, Y.-K., Maddick, I. H., Griffiths, G. S., and Gilthorpe, M. S. (2004). Mathematical
coupling can undermine the statistical assessment of clinical research: illustration
from the treatment of guided tissue regeneration. Journal of Dentistry 32, 133-142.
Twisk, J. W. R. (2004). Longitudinal Data Analysis. A Comparison Between Generalized
Estimating Equations and Random Coefficient Analysis. European Journal of
Epidemiology 19, 769-776.
Utzschneider, K. M., Carr, D. B., Barsness, S. M., Kahn, S. E., and Schwartz, R. S.
(2004). Diet-Induced Weight Loss Is Associated with an Improvement in β-Cell
Function in Older Men. Journal of Clinical Endocrinology & Metabolism 89,
2704-2710.
Vaarasmaki, M., Pouta, A., Elliot, P., et al. (2009). Adolescent manifestations of
metabolic syndrome among children born to women with gestational diabetes in a
general-population birth cohort. American journal of epidemiology 169, 1209-
1215.
Wang, C. Y., Huang, Y., Chao, E. C., and Jeffcoat, M. K. (2008). Expected Estimating
Equations for Missing Data, Measurement Error, and Misclassification, with
Application to Longitudinal Nonignorable Missing Data. Biometrics 64, 85-95.
Wang, N., and Davidian, M. (1996). A Note on Covariate Measurement Error in
Nonlinear Mixed Effects Models. Biometrika 83, 801-812.
Wang, N., Lin, X., Gutierrez, R. G., and Carroll, R. J. (1998). Bias Analysis and SIMEX
Approach in Generalized Linear Mixed Measurement Error Models. Journal of
the American Statistical Association 93, 249-261.
163
Wang, N., Lin, X., and Guttierrez, R. B. (1999). A bias correction regression calibration
approach in generalized linear mixed measurement error models.
Communications in Statistics - Theory and Methods 28, 217-233.
Wattmo, C., Wallin, A., Londos, E., and Minthon, L. (2011). Predictors of long-term
cognitive outcome in Alzheimer's disease. Alzheimer's Research & Therapy 3, 23.
Weyer, C., Tataranni, P. A., Bogardus, C., and Pratley, R. E. (2001). Insulin Resistance
and Insulin Secretory Dysfunction Are Independent Predictors of Worsening of
Glucose Tolerance During Each Stage of Type 2 Diabetes Development. Diabetes
Care 24, 89-94.
WHO (2000). Obesity: preventing and managing the global epidemic. Report of a WHO
consultation. World Health Organ Tech Rep Ser 894, i-xii, 1-253.
WHO (2003). Diet, Nutrition and Prevention of Chronic Diseases. In Report of the joint
WHO/FAO expert consultation. WHO Technical Report Series No 916 (TRS 916).
Winer, B. J. (1971). Statistical principles in experimental design, second edition: New
York: McGraw-Hill.
Wishart, J. (1938). Growth-Rate Determinations in Nutrition Studies with the Bacon Pig,
and Their Analysis. Biometrika 30, 16-28.
Xiang, A. H., Kawakubo, M., Trigo, E., Kjos, S. L., and Buchanan, T. A. (2010).
Declining β-Cell Compensation for Insulin Resistance in Hispanic Women With
Recent Gestational Diabetes Mellitus. Diabetes Care 33, 396-401.
Xiang, A. H., Peters, R. K., Trigo, E., Kjos, S. L., Lee, W. P., and Buchanan, T. A.
(1999). Multiple metabolic defects during late pregnancy in women at high risk
for type 2 diabetes. Diabetes 48, 848-854.
Xiang, A. H., Wang, C., Peters, R. K., Trigo, E., Kjos, S. L., and Buchanan, T. A. (2006).
Coordinate Changes in Plasma Glucose and Pancreatic β-Cell Function in Latino
Women at High Risk for Type 2 Diabetes. Diabetes 55, 1074-1079.
Yanez, N. D., Kronmal, R. A., and Shemanski, L. R. (1998). The effects of measurement
error in response variables and tests of association of explanatory variables in
change models. Statistics in Medicine 17, 2597-2606.
Yates, F. (1935). Complex Experiments. Supplement to the Journal of the Royal
Statistical Society 2, 181-247.
Yoshioka (1998). Use of randomisation in the Medical Research Council's clinical trial of
streptomycin in pulmonary tuberculosis in the 1940s. BMJ 317, 1220-1223.
Abstract (if available)
Abstract
Evaluating the associations between the baseline value and other exposure variables with the longitudinal trajectory constantly occurs in the medical research. Traditional analyses using the observed measurements can produce biased estimates of the true relationships when responses are measured with error. Bias-correction formulas are developed using the linear regression for studies with a baseline and one follow-up visit. When the data has multiple follow-up visits, random intercept and random slope mixed effects models are always used. However, current bias-correction methods developed from aforementioned mixed effects model have their limitations. Motivated by a real data example, we developed a new method to correct the bias in association estimates using the mixed effects model. The validity of formulas was evaluated by the simulation study. Two methods were proposed to test the bias-corrected estimates. Finally, two real study examples were provided to demonstrate the application of our methods.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Effect of measurement error on the association between baseline and longitudinal change
PDF
Relationship between progression of atherosclerosis and coagulation measures in a randomized-controlled trial
PDF
Inference correction in measurement error models with a complex dosimetry system
PDF
Association between endogenous sex hormone levels and cognition: the Women’s Isoflavone Soy Health (WISH) trial
PDF
Bayesian multilevel quantile regression for longitudinal data
PDF
An assessment of necrosis grading in childhood osteosarcoma: the effect of initial treatment on prognostic significance
PDF
Association of carotid artery stiffness with measures of cognition in older adults
PDF
ROC surface in the presence of verification bias
PDF
Comparing robustness to outliers and model misspecification between robust Poisson and log-binomial models
PDF
The association of asthma and asthma-related medications on subclinical atherosclerosis: a cross-sectional analysis of four randomized clinical trials
PDF
Associations between physical activities with bone mineral density in postmenopausal women
PDF
Correcting for shared measurement error in complex dosimetry systems
PDF
Applications of multiple imputations in survival analysis
PDF
The association between self-reported physical activity and cognition in elderly clinical trial participants
PDF
The association of apolipoprotein E genotype with cholesterol levels and atherosclerosis: three randomized clinical trials
PDF
Relationship of blood pressure and antihypertensive medications to cognitive change in the BVAIT, WISH, and ELITE clinical trials
PDF
Shortcomings of the genetic risk score in the analysis of disease-related quantitative traits
PDF
Estimation of nonlinear mixed effects mixture models with individually varying measurement occasions
PDF
Associations between isoflavone soy protein (ISP) supplementation and breast cancer in postmenopausal women in the Women’s Isoflavone Soy Health (WISH) clinical trial
PDF
The evaluation of the long-term effectiveness of zero/low fluoroscopy workflow in ablation procedures for the treatment of paroxysmal and persistent atrial fibrillation
Asset Metadata
Creator
Chen, Zhanghua
(author)
Core Title
Evaluating the associations between the baseline and other exposure variables with the longitudinal trajectory when responses are measured with error
School
Keck School of Medicine
Degree
Doctor of Philosophy
Degree Program
Biostatistics
Publication Date
11/26/2013
Defense Date
10/15/2013
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
baseline adjustment,baseline prediction,longitudinal analysis,measurement error,mixed effects model,OAI-PMH Harvest,regression to the mean
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Watanabe, Richard M. (
committee chair
), Buchanan, Thomas A. (
committee member
), Mack, Wendy Jean (
committee member
), Stram, Daniel O. (
committee member
), Xiang, Anny H. (
committee member
)
Creator Email
zhanghuc@gmail.com,zhanghuc@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c3-351704
Unique identifier
UC11296370
Identifier
etd-ChenZhangh-2187.pdf (filename),usctheses-c3-351704 (legacy record id)
Legacy Identifier
etd-ChenZhangh-2187.pdf
Dmrecord
351704
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Chen, Zhanghua
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
baseline adjustment
baseline prediction
longitudinal analysis
measurement error
mixed effects model
regression to the mean