EXTENDING LONGITUDINAL ANALYSIS

bagimpertinentUrban and Civil

Nov 16, 2013 (3 years and 11 months ago)

250 views

EXTENDINGLONGITUDINAL ANALYSIS
TOESTIMATION OF COEFFICIENTS IN HIV-1 DYNAMICS
A Thesis
Presented to the
Faculty of
San Diego State University
In Partial Fulfillment
of the Requirements for the Degree
Master of Science
in
Statistics
by
Suzanne Louise Papp
Fall 2010
iii
Copyright 2010
by
Suzanne Louise Papp
iv
DEDICATION
To my husband,my mother,and the memory of my father.
v
I amonly one,
But still I amone.
I cannot do everything,
But still I can do something;
And because I cannot do everything
I will not refuse to do the something that I can do.
–Edward Everett Hale
vi
ABSTRACT OF THE THESIS
Extending Longitudinal Analysis
to Estimation of Coefficients in HIV-1 Dynamics
by
Suzanne Louise Papp
Master of Science in Statistics
San Diego State University,2010
Deterministic dynamic equations have proved very important in describing the
response of HIV-1 to treatment in vivo.This thesis develops estimated local linear
mixed-effects methods to estimate the coefficients of a linear differential equation that models
the rate of HIV-1 RNA concentration in the patients’ plasma.Two methods are developed.In
the first,the coefficients are assumed to be constant.In the second,the coefficients are
allowed to vary with time.For both methods the first step is to estimate the viral load and the
first time derivative of the viral load fromthe measured value of viral load,where the
measured value may have error.A local polynomial kernel estimation with mixed-effects is
used for these estimations.The estimates of viral load and the first time derivative of viral
load are then used in estimating the coefficients.Where the the coefficients are assumed
constant,the second step estimates these coefficients with a linear mixed-effects equation.
Where the coefficients are assumed to be time-varying,a local polynomial kernel estimation
with mixed-effects is used to estimate the coefficients.We use a real data application and also
numerical simulations to illustrate the two estimation methods.
vii
TABLE OF CONTENTS
PAGE
ABSTRACT....................................................................................vi
LIST OF TABLES..............................................................................ix
LIST OF FIGURES............................................................................x
ACKNOWLEDGMENTS.....................................................................xi
CHAPTER
1 INTRODUCTION.....................................................................1
1.1 History of HIV-1................................................................1
1.2 An HIV-1 Differential Equation Model........................................2
1.3 Clinical Study Example.........................................................3
1.4 Goals and Organization of the Thesis..........................................4
2 LITERATURE REVIEW..............................................................6
2.1 Virus Dynamics.................................................................6
2.2 HIV-1 and Mixed Models for Longitudinal Analysis.........................9
3 ESTIMATING CONSTANT COEFFICIENTS......................................11
3.1 Proposed Methods...............................................................11
3.1.1 Estimating X
0
(t) and X(t) with Mixed-Effects Curve Modeling.......11
3.1.2 Estimating a
i1
and a
i2
with Linear Mixed-Effects Model...............12
3.2 Application to Real Data........................................................15
3.2.1 Estimating X
0
(t) and X(t) in a Real Data Example.....................15
3.2.2 Estimating a
1
and a
2
in a Real Data Example............................19
3.3 Simulation for a Simplified Model with Constant a
1
..........................21
4 ESTIMATING TIME-VARYING COEFFICIENTS.................................27
4.1 Proposed Methods...............................................................27
4.2 Application to a Real Data Example...........................................30
4.3 Simulation for a Simplified Model with a
1
(t).................................34
5 BANDWIDTH SELECTION..........................................................41
5.1 Choosing a Bandwidth for
^
X
i
(t
ij
) and
^
X
0
i
(t
ij
)................................41
5.2 Choosing a Bandwidth for a
1
(t
ij
) and a
2
(t
ij
).................................43
viii
6 DISCUSSION AND CONCLUSION.................................................48
BIBLIOGRAPHY..............................................................................50
APPENDICES
A R-CODE FOR STEP 1.................................................................53
B R-CODE FOR TIME-INVARIANT STEP 2.........................................59
C R-CODE FOR CONSTANT a
1
SIMULATION MODEL...........................62
D R-CODE FOR TIME-VARYING STEP 2............................................77
E R-CODE FOR a
1
(t) SIMULATION MODEL.......................................83
F R-CODE FOR GCV STEP 1..........................................................102
G R-CODE FOR GCV STEP 2..........................................................109
ix
LIST OF TABLES
PAGE
Table 3.1 Estimates of Time-Invariant Coefficients...........................................20
Table 3.2 Time-Invariant ^a
1
and ^a
2
Estimated With and Without Outlier...................21
Table 3.3 The Estimates for a
1
fromSimulation Models.....................................23
Table 4.1 RASE for a
1
(t) Estimation..........................................................40
Table 5.1 GCV Values for
^
X
i
(t
ij
) and
^
X
0
i
(t
ij
)................................................44
Table 5.2 GCV Values for a
1
(t
ij
) and a
2
(t
ij
).................................................47
x
LIST OF FIGURES
PAGE
Figure 1.1 The days on which data were collected for all patients...........................3
Figure 1.2 The number of samples per patient................................................4
Figure 1.3 A scatterplot matrix fromthe real data example..................................5
Figure 3.1 Estimation of log
10
(RNA) and the first derivative................................16
Figure 3.2 Estimation of log
10
(RNA) and the first derivative without patient 262874......17
Figure 3.3 Plots of log
10
(RNA) vs.day and the first derivative of log
10
(RNA)
vs.day for representative patients.....................................................18
Figure 3.4 Plots of log
10
(RNA) vs.day and the first derivative of log
10
(RNA)
vs.day for patient outlier...............................................................19
Figure 3.5 Simulation points of log
10
(RNA) for simulation with constant a
1
...............24
Figure 3.6 Estimation of log
10
(RNA) for simulation with constant a
1
.......................25
Figure 3.7 Estimation of the first derivative of log
10
(RNA) for simulation with
constant a
1
..............................................................................26
Figure 4.1 Fixed plot and patient plots for log
10
(RNA),a
1
,and a
2
..........................31
Figure 4.2 Fixed plot and patient plots for log
10
(RNA),a
1
,and a
2
estimated
without patient 262874.................................................................32
Figure 4.3 Plots of a
1
(t) vs.day and a
2
(t) vs.day for representative patients..............33
Figure 4.4 Plots of a
1
(t) vs.day and a
2
(t) vs.day for patient outlier........................34
Figure 4.5 Simulation points of log
10
(RNA) for simulation with a
1
(t)......................36
Figure 4.6 Estimation of log
10
(RNA) for simulation with a
1
(t).............................37
Figure 4.7 Estimation of the first derivative of log
10
(RNA) for simulation with a
1
(t)......38
Figure 4.8 Estimation of a
1
(t) for the simulation.............................................39
Figure 5.1 The GCV vs.bandwidth for step 1................................................45
Figure 5.2 The GCV vs.bandwidth for step 2................................................47
xi
ACKNOWLEDGMENTS
I would like to thank Dr.Jianwei Chen for his help and encouragement.I appreciate
the time and effort of Dr.John Alcaraz and Dr.Chii-Dean Lin in serving on my committee.
The advice and suggestions of Kevin M.Johnson and Dr.Jeanette Papp were invaluable.I
would like to thank Cynthia Royce for editing my thesis.Finally,I especially appreciate the
patience and support of my husband and mother.
1
CHAPTER 1
INTRODUCTION
This thesis applies local polynomial mixed effects methods to estimate the coefficients
of two dynamic models for the rate of the viral load of Human Immunodeficiency Virus Type
1 (HIV-1).Mixed effects will allow us to extend these equations frommodeling a single
patient to modeling all patients in a clinical study.The remaining sections of the introduction
are as follows.Section 1.1 gives a brief history of the HIV-1 virus and virus dynamics.
Section 1.2 introduces the dynamic models this thesis addresses.Section 1.3 describes the
clinical study that provides data for our real application estimates.Section 1.4 states the goals
of this thesis and outlines the organization of the remaining chapters.
1.1 HISTORY OF HIV-1
In 1980 physicians in Los Angeles,Atlanta,New York,and other U.S.cities were
puzzled to find relatively young and previously healthy men coming down with opportunistic
infections that generally indicate a compromised immune system[11,22].In 1983 the HIV-1
virus was identified by [1] and [10].In 1985,a related virus,HIV-2 was discovered [14,22].
HIV-2 is much less infective and less virulent than HIV-1 and is currently largely confined to
West Africa [14].When HIV-1 was first identified there was great optimismthat an effective
vaccine would be discovered within two years;unfortunately,development of a vaccine has
proved remarkably difficult [20].As of 2010 work on an effective HIV-1 vaccine is ongoing.
The development of antiretroviral drugs has been successful in extending the lives of
HIV infected patients and delaying the onset of Acquired Immune Deficiency Syndrome
(AIDS) [8].The first drug successfully used for AIDS treatment was azidothymidine (AZT)
which was licensed to treat AIDS in 1987 [11].Treatment with AZT is initially beneficial.
However,patients commonly developed AZT-resistant HIV-1 strains [11,22] and thus,acting
alone,it lost its effectiveness against the virus.The HIV-1 virus is particularly adept at
developing drug resistance because of its high rate of mutation and short life span [22].A
combination of different antiretroviral drugs,termed Highly Active Antiretroviral Therapy
(HAART) was developed in the mid 1990s [8,22].HAART has been effective in slowing
drug resistance and increasing life expectancy for HIV infected individuals [8].The
development of drug resistance still remains a problem,which is why the development of new
antiretrovirals is essential in reducing illness and extending lives [8].The AIDS epidemic
2
remains a major world problemwith an estimated 33.4 million HIV infected individuals and
2.0 million AIDS related deaths in 2008 [29].
Virus dynamics have had a significant effect on both the understanding and the
treatment of HIV-1 infection.The average time for AIDS to develop fromtime of infection is
10 years,but virus dynamics helped uncover the fact that many important processes in the
course of the infection were occurring within hours,weeks,or months and assign quantitative
values to these processes [22,23].Models of the reaction between the virus and the immune
systemhave influenced the treatment of HIV infection with antiretrovirals [22,23].
1.2 AN HIV-1 DIFFERENTIAL EQUATION MODEL
The rate of change of viral load for the patients in a longitudinal HIV-1 study for m
patients with n
i
time points each can be modeled by the equations,
Y
i
(t
ij
) = X
i
(t
ij
) +
i
(t
ij
);
d
dt
X
i
(t
ij
) = a
i1
+a
i2
CD4
i
(t
ij
) cX
i
(t
ij
);
i = 1;:::;m;j = 1;:::;n
i
;
(1.1)
where X
i
(t
ij
) is log
10
of the viral load in plasma,Y
i
(t
ij
) is the measured value of X
i
(t
ij
),
CD4
i
(t
ij
) is the concentration of CD4+ T cells,and c is the clearance rate of free virus.If we
allow the coefficients,a
i1
and a
i2
to vary with time,this becomes,
Y
i
(t
ij
) = X
i
(t
ij
) +
i
(t
ij
);
d
dt
X
i
(t
ij
) = a
i1
(t
ij
) +a
i2
(t
ij
)CD4
i
(t
ij
) cX
i
(t
ij
);
i = 1;:::;m;j = 1;:::;n
i
:
(1.2)
This thesis advances two estimated linear mixed-effects models.The first estimates the
coefficients a
i1
and a
i2
for Equation (1.1).The second estimates a
i1
(t
ij
) and a
i2
(t
ij
) for
Equation (1.2).For both equations,the first step is to use a local polynomial mixed-effects
method developed in [33] to derive the estimates,
^
X
i
(t
ij
) and
^
X
0
i
(t
ij
),of X
i
(t
ij
) and
dX
i
(t
ij
)=dt,respectively.The second step for Equation (1.1) is to use these results in fitting a
linear mixed model of dX
i
(t
ij
)=dt.Fromthis model we get an estimate of a
i1
and a
i2
The
second step for Equation (1.2) extends mixed-effects methods to the local estimation method
for estimating a
1
(t) and a
2
(t) that was developed in [3].We use the estimated values of
^
X
i
(t
ij
) and
^
X
0
i
(t
ij
) to fit a local polynomial mixed-effects model of dX
i
(t
ij
)=dt to get a
time-varying estimate of a
i1
(t
ij
),and a
2i
(t
ij
).Determining an appropriate bandwidth is
essential in local polynomial methods.We will use a generalized cross-validation (GCV)
method to choose an appropriate bandwidth for both steps.
3
1.3 CLINICAL STUDY EXAMPLE
We apply the methods developed to a real data example in Section 3.2 and Section 4.2.
The data is fromthe ACTG 398 study,a Phase II clinical study of the protease inhibitor,
Amprenavir,as part of a treatment for HIV-1 patients.The patients enrolled in this study had
stopped responding to their current treatment,as measured by a plasma HIV-1 RNA
concentration of at least 1,000 copies/ml.In addition,the patients’ history of treatment was
required to conformto specified protease inhibitors.Clinical visits to gather data were
scheduled at baseline,at weeks 2,4,8,16,24,32,40,48,56,64,and 72.Some variation from
this schedule took place in practice.Figure 1.1 is a histogramof the days on which data was
collected.Many more data collections were made at the beginning of the study.Data
collection tails off towards the end of the study.Figure 1.2 is a histogramof the number of
sample collections per patient.Most patients had between five and twelve sample collections.
The largest number of collections for a patient was fifteen.The treatment included
Amprenavir as part of a cocktail of drugs.The study had four different treatment arms,
including the control treatment.The 35 patients used in this thesis were those in the first arm,
treated with Amprenavir,Saquinavir,Abacavir,Efavirenz,and Adefovir.For this thesis we
use the variables day,CD4,and log
10
(RNA).A scatterplot matrix of these variables can be
seen in Figure 1.3.There is a negative relationship between CD4+ T cells and log
10
(RNA).If
there is any relationship between day and log
10
(RNA) or day and CD4+ T cells,it is not seen
in a simple scatterplot.We hope to expose a relationship between themin the following
chapters.
Day
Frequency
0
100
200
300
400
500
0
5
15
25
35
Figure 1.1.The days on which data were
collected for all patients.
4
Number of samples
Frequency
0
5
10
15
0
1
2
3
4
5
6
Figure 1.2.The number of samples per
patient.
1.4 GOALS AND ORGANIZATION OF THE THESIS
The goal of this thesis is to apply mixed effect longitudinal analysis to Equation 1.2
and Equation 1.1 and assess the performance of the two models.Mixed effects model the
dependent variables as a function of the independent variables where the residuals are
normally distributed but are not necessarily independent or constant.In our model,the
residuals for each individual patient are correlated.The function consists of a fixed part,
applicable to all the patients,and a randompart,uniquely applicable to each patient.By using
mixed modeling we can use all the data froma clinical trial while considering the differences
between individual patients.This method gives more information than simply modeling a
single patient.
The remainder of this thesis is organized as follows.Chapter 2 is a literature review of
virus dynamics and of mixed models applied to the HIV-1 viral load.Chapter 3 extends
longitudinal methods to Equation (1.1).Chapter 4 extends longitudinal methods to
Equation (1.2).Chapter 5 develops the generalized cross-validation (GCV) method for
choosing the two bandwidths.The first bandwidth is require to estimate the viral load and the
rate of change of viral load.The second bandwidth is required to estimate the tim-varying
coefficients of Equation (1.2).Chapter 6 discusses the results of the thesis,draws conclusions,
and suggests further areas of study.
5
Figure 1.3.A scatterplot matrix fromthe real data
example.
6
CHAPTER 2
LITERATURE REVIEW
2.1 VIRUS DYNAMICS
In 1995 two papers appeared in the same issue of Nature that may be said to have
started the study of HIV-1 dynamics.Both papers used the effect of antiretroviral drugs in
vivo and a mathematical model of the virus to quantify the interaction between the immune
systemand HIV-1.Wei and Shaw [30] presented equations modeling infected cells and viral
load in plasma,
T

(t) = T

(0)e
at
X(t) = X(0)
ue
at
ae
ut
u a
:
(2.1)
where T

are the infected cells,X is the viral load in plasma,a is the rate constant for the
decline in the infected cells,and u is the rate constant of the decline in free virus.Meanwhile,
Ho and Markowitz [12] used the following differential equation to model the rate of change
virus load:
d
dt
X(t) = P cX(t);(2.2)
where X is the viral load in plasma,P is the virus production rate,and c is the clearance rate
of free virus.Before administering antiretrovirals,we assume dX=dt to be approximately
zero,so P = cX(t).If antiretroviral treatment perfectly suppresses virus production then
dX=dt = cX.
A model of short termvirus dynamics on the scale of a week was introduced by
Perelson,et al.[24],
X(t) = X
0
e
ct
+
cX
0
c 

c
c 

e
t
e
ct

te
ct

;(2.3)
where c is the clearance rate of free virus and  is the death rate of infected cells.If short term
measurements of viral load are available,Equation (2.3) can be used to estimate  and c.
A set of deterministic dynamic equations in the absence of treatment was proposed by
Nowak and Banghamin [21],
7
d
dt
T(t) =  T(t) kT(t)X(t)
d
dt
T

(t) = kT(t)X(t) T

(t)
d
dt
X(t) = BT

(t) cX(t);
(2.4)
where T(t) is the concentration of uninfected CD4+ T cells,T

(t) is the concentration of
infected CD4+ T cells,and X(t) is the concentration of free virus.When modeling HIV-1
infection the cells are CD4+ T cells.The  parameter is generation rate of new CD4+ T cells
cells; is the death rate of CD4+ T cells,k is the infection rate of CD4+ T cells;r(t) is
anti-viral drug efficacy; is the death rate of infected CD4+ T cells;B is the burst size (the
number of new virions produced froman infected CD4+ T cell);c is the clearance rate of free
virions.Adding the effect of anti-viral treatment,r(t),we may write Equation (2.5):
d
dt
T(t) =  T(t) k [1 r(t)] T(t)X(t)
d
dt
T

(t) = k [1 r(t)] T(t)X(t) T

(t)
d
dt
X(t) = BT

(t) cX(t):
(2.5)
Wu and Ding [32] and Huang,Liu,and Wu [13] apply nonlinear mixed-effects methods to
these equations to estimate time-varying coefficients.Li,et al.[16] proposed a method using
cubic splines to estimate time-varying coefficients of a set of deterministic dynamic equations.
Liang and Wu [18] used local smoothing and a pseudo-least squares approach to estimate
constant coefficients in (2.5).
Following [3] we combine the second and third equation of 2.5 to get
d
dt
X(t) = B
d
dt
T

(t) +Bk [1 r(t)] X(t)T(t) cX(t):(2.6)
Let
a
1
(t) = B
d
dt
T

(t)
a
2
(t) = Bk [1 r(t)] X(t)T(t);
and assume the uninfected CD4+ T cells are proportional to the measured CD4+ T cells,then
we have the equation:
d
dt
X(t) = a
1
(t) +a
2
(t)CD4(t) cX(t):(2.7)
8
An efficient method for estimating a
1
(t) and a
2
(t) using a two-step,local polynomial
estimation was developed by Chen and Wu in [3,4].A brief overview of those portions of [3]
as it relates to this thesis is given in the remainder of this section.
If the actual log of the viral load is X(t) then the measured value,Y (t),includes some
error:
Y (t) = X(t) +
1
(t):(2.8)
In the first step of this method,X(t) and dX(t)=dt are estimated by making use of a Taylor
expansion and local polynomial kernel estimators.The following locally weighted least
squares function is minimized:
n
X
i=1
fY (t
i
) [
0
(t) +
1
(t) (t
i
t)]g
2
K
h
1
(t
i
t);
(2.9)
where K
h
1
is a kernel function,K
h
1
(t) = K(t=h
1
)=h
1
,and h
1
is the bandwidth.The most
usual Kernel functions are symmetric probability functions,such as the Gaussian,
Epanechnikov,and Uniform[7].Local polynomial modeling is generally not sensitive to the
Kernel function selected.Let X
0
(t) = dX(t)=dt.Then for a given t,X(t) and X
0
(t),the
parameters
^

0
and
^

1
are estimated as follows:
^
(t) =

^

0
(t);
^

1
(t)

T
=

Z
1
T
K
1
Z
1

1
Z
1
T
K
1
Y;(2.10)
where the weighting function,K
1
= diag(K
h
1
(t
1
t);:::;K
h
1
(t
n
t)) and the design
matrix,Z
1
,is:
Z
1
=
0
B
B
@
1 t
1
t
.
.
.
.
.
.
1 t
n
t
1
C
C
A
:
In the second step,we estimate a
1
and a
2
.A substitution of
^
X(t) for X(t) and
^
X
0
(t) for X
0
(t)
is made in Equation (2.7).As these are estimates,not the actual values,an error term,
2
(t)
must be added to give,
^
X
0
(t) = a
1
(t) +a
2
(t)CD4 c
^
X(t) +
2
(t):(2.11)
Again,we make use of the Taylor expansion and local polynomial kernel estimators.
A local estimate for a
1
and a
2
is made at t by minimizing the locally weighted least squares
function,
n
X
i=1
"
^
X
0
(t
i
) +c
^
X(t
i
) 
1
X
r=0
(a
1
(t) +CD4(t
ij
)a
2
(t)) (t
i
t)
r
#
2
K
h
2
(t
i
t);(2.12)
9
where K is a kernel function,K
h
2
(t) = K(t=h
2
)=h
2
,and h
2
is the bandwidth.The local
linear estimator of a
1
and a
2
is,

a
1
(t)
a
2
(t)
!
=

1 0 0 0
0 0 1 0
!

Z
2
T
K
2
Z
2

1
Z
2
T
K
2
(U
1
+cU
0
);(2.13)
where the weighting function is K
2
= diag(K
h
2
(t
1
t):::K
h
2
(t
n
t)),and U
0
,U
1
,and Z
2
are defined as:
U
0
=
0
B
B
@
^
X(t
1
)
.
.
.
^
X(t
n
)
1
C
C
A
U
1
=
0
B
B
@
^
X
0
(t
1
)
.
.
.
^
X
0
(t
n
)
1
C
C
A
Z
2
=
0
B
B
@
1 (t
1
t) CD4(t
1
) CD4(t
1
)(t
1
t)
.
.
.
.
.
.
.
.
.
.
.
.
1 (t
n
t) CD4(t
n
) CD4(t
n
)(t
n
t)
1
C
C
A
:
This two-step method will be extended in this thesis to longitudinal data using mixed-effects.
2.2 HIV-1 AND MIXED MODELS FOR
LONGITUDINAL ANALYSIS
In the mixed-effects model for longitudinal data described by Laird and Ware [15] the
i
th
individual is modeled as,
Y
i
= X
i
 +Z
i
b
i
+
i
  N(0;R
i
)
b
i
 N(0;D);
i = 1;:::;m;
(2.14)
where Y
i
is the response vector, is the vector of fixed effects,b
i
is the vector of random
effects,X
i
is the fixed effects design matrix,Z
i
is the randomeffects design matrix,and 
i
is
the error.
Shi and Taylor [27] develop a method of fitting linear combinations of B-splines to
both the fixed and randompart of unbalanced longitudinal data in which the randompart is
fitted to each patient.Rice and Wu [26] also presented a spline fitting method of applying
mixed-effects to unbalanced longitudinal data
In this thesis we use local polynomial kernel regression [7,9].Lin and Carroll [19]
introduced a local polynomial kernel regression method to longitudinal data.Their method
assumes no within-patient error correlation.
10
Wu and Zhang [33] make no requirements on the within-patient correlation and so
develop a more general method that can more accurately model longitudinal data.They model
the measured value,Y
i
,as a mixed effect model:
Y
i
(t
ij
) = X(t
ij
) +
i
(t
ij
) +
i
(t
ij
);i = 1;:::;m;j = 1;:::;n
i
;(2.15)
for mpatients who each have n
i
time points.The fixed and randomeffects for any t are
estimated using a first degree Taylor expansion,
X(t
ij
)  X(t) +X
0
(t)(t
ij
t) = Z
T
ij
;

i
(t
ij
)  
i
(t) +
i
0
(t) (t
ij
t) = Z
T
ij
b
i
;
i = 1;  ;m;j = 1;  ;n
i
;
(2.16)
where X
0
(t) = dX(t)=dt and 
0
i
(t) = d
i
(t)=dt.They can then estimate Y
i
as:
Y
i
(t
ij
) = Z
i
(t
ij
) ((t
ij
) +b
i
(t
ij
)) +
i
(t
ij
);i = 1;:::;m;j = 1;:::;n
i
:(2.17)
To estimate  and b
i
they minimize the objective function,(see [6]):
n
X
i=1
[Y
i
Z
i
( +b
i
)]
T
K
1=2
ih
R
i
1
K
1=2
ih
[Y
i
Z
i
( +b
i
)] +
b
T
i
D
1
b
i
+logjDj +logjR
i
j;
(2.18)
where the kernel weight,K
ih
= diag (K
h
(t
i1
t):::K
h
(t
in
i
t)),K
h
(t) = K(t=h)=h,and
Z
i
=
0
B
B
@
1 (t
i1
t)
.
.
.
.
.
.
1 (t
in
i
t)
1
C
C
A
;i = 1;  ;m:
Solving for  and b they get the equations (2.19) and (2.20):
^
 = f
m
X
i=1
Z
i
T


ih
Z
i
g
1
f
m
X
i=1
Z
i
T


ih
Y
i
g (2.19)
^
b
i
= fZ
i
T

ih
Z
i
+D
1
g
1
Z
i
T

ih
(Y
i
Z
i
^
);(2.20)
where


ih
= K
1=2
ih
V
1
i
K
1=2
ih
V
i
= K
1=2
ih
Z
T
i
K
1=2
ih
+R
i

ih
= K
1=2
ih
R
1
i
K
1=2
ih
i = 1;  ;m:
11
CHAPTER 3
ESTIMATINGCONSTANT COEFFICIENTS
In this chapter we extend longitudinal methods to Equation (1.1).The equation that
models the measured value of viral load is,
Y
i
(t
ij
) = X
i
(t
ij
) +
i1
(t
ij
);
i = 1;:::;m;j = 1;:::;n
i
:
(3.1)
The equation that models the rate of change of viral load is,
d
dt
X
i
(t
ij
) = a
i1
+a
i2
CD4
i
(t
ij
) cX
i
(t
ij
);
i = 1;:::;m;j = 1;:::;n
i
;
(3.2)
where X
i
(t
ij
) is log
10
of the viral load in plasma,Y
i
(t
ij
) is the measured value of X
i
(t
ij
),
CD4
i
(t
ij
) is the concentration of CD4+ T cells,and c is the clearance rate of free virus.
Section 3.1 develops the methods,Section 3.2 applies the methods to the Amprenavir clinical
study,and Section 3.3 further illustrates the methods with numerical simulations.
3.1 PROPOSED METHODS
Section 3.1.1 develops the equations for Step 1,estimating dX
i
(t
ij
)=dt and X
i
(t
ij
).
The method follows [33] for the special case where within-patient correlation is assumed to
be negligible and within-patient errors are homoscedastic.Section 3.1.2 describes the
equation for Step 2,estimating a
i1
and a
i2
.
3.1.1 Estimating X
0
(t) and X(t) with Mixed-Effects
Curve Modeling
We model the measured value,Y
i
,fromEquation (3.1) using the method of [33]
outlined in Section 2.2:
Y
i
(t
ij
) = X(t
ij
) +
i
(t
ij
) +
i1
(t
ij
);
X
i
(t
ij
) = X(t
ij
) +
i
(t
ij
);
i = 1;:::;m;j = 1;:::;n
i
:
(3.3)
We will assume the special case where the error,
i1
,is independent and identically
distributed,that is,
i1
 N(0;
1
I
n
i
),where I
n
i
is the n
i
n
i
identity matrix.The fixed and
randomeffects for any t are estimated using a first degree Taylor expansion as in
12
Equation (2.16).The equations are repeated here:
X(t
ij
)  X(t) +X
0
(t)(t
ij
t) = Z
T
ij1

1
;

i
(t
ij
)  
i
(t) +
i
0
(t)(t
ij
t) = Z
T
ij1
b
i1
;
i = 1;:::;m;j = 1;:::;n
i
;
where X
0
(t) = dX(t)=dt,
0
i
(t) = d
i
(t)=dt,Z
ij1
= (1;t
ij
t)
T
,the vector of fixed effects is

1
= (X(t);X
0
(t))
T
,and the vector of randomeffects is b
i1
= (
i
(t);
i
0
(t))
T
.Furthermore,
we assume the distribution of b
i1
,i = 1;  ;m,to be normal,where b
i1
 N(0;D
1
) and,
D
1
=


11

12

12

22
!
The response matrix is,
Y
i
=
0
B
B
@
Y
i
(t
i1
)
.
.
.
Y
i
(t
in
i
)
1
C
C
A
;i = 1;  ;m:
The design matrix,Z
i1
,for any given t is,
Z
i1
(t) =
0
B
B
@
1 (t
i1
t)
.
.
.
.
.
.
1 (t
in
i
t)
1
C
C
A
;i = 1;  ;m:(3.4)
To estimate 
1
and b
i1
we minimize the objective function,Equation (2.18),repeated here:
m
X
i=1
[Y
i
Z
i1
(
1
+b
i1
)]
T
K
1=2
ih
1
R
1
i2
K
1=2
ih
1
[Y
i
Z
i1
(
1
+b
i1
)] +
b
T
i1
D
1
1
b
i1
+logjD
1
j +logjR
i1
j;
where the kernel weight K
ih
1
= diag(K
h
1
(t
i1
t);:::;K
h
1
(t
in
i
t)),K
h
1
= K(t=h
1
)=h
1
,and
K is a kernel function.Then we can use R software [28] to estimate
^

1
(t) and
^
b
i1
(t) for any
given t.We use the R function,lme,fromthe nlme package [25] to estimate the viral load and
the first derivative of the viral load as a linear mixed-effects equation.
3.1.2 Estimating a
i1
and a
i2
with Linear
Mixed-Effects Model
In the second step we estimate the constant coefficients,a
i1
and a
i2
,from
Equations (3.1) and (3.2),repeated here:
Y
i
(t
ij
) = X
i
(t
ij
) +
i1
(t
ij
)
d
dt
X
i
(t
ij
) = a
i1
+a
i2
CD4
i
(t
ij
) cX
i
(t
ij
)
i = 1;:::;m;j = 1;:::;n
i
:
13
Equation (3.2) can be rearranged as:
d
dt
X
i
(t
ij
) +cX
i
(t
ij
) = a
i1
+a
i2
CD4
i
(t
ij
);
i = 1;:::;m;j = 1;:::;n
i
:
(3.5)
We will make the fixed and randompart of the parameters a
i1
and a
i2
explicit as follows:
d
dt
X
i
(t
ij
) +cX
i
(t
ij
) = a
1
+
i1
+CD4
i
(t
ij
)(a
2
+
i2
);
i = 1;:::;m;j = 1;:::;n
i
:
(3.6)
However dX
i
(t
ij
)=dt and X
i
(t
ij
) are unknown,so we will substitute the estimates derived in
Section 3.1.1 to get an estimated linear mixed-effects model:
^
X
i
0
(t
ij
) +c
^
X
i
(t
ij
) = a
1
+
i1
+CD4
i
(t
ij
)(a
2
+
i2
) +
i2
(t
ij
);
i = 1;:::;m;j = 1;:::;n
i
;
(3.7)
where 
i2
(t) is an error termintroduced as we are using the estimates for dX
i
(t
ij
)=dt and
X
i
(t
ij
) rather than the (unknown) actual values.We will assume the error to be independent
and identically distributed.The distributions of 
i1
and 
i2
are thus:


i1

i2
!
 N(0;D
2
);D
2
=


11

12

12

22
!
;
i = 1;:::;m;j = 1;:::;n
i
;
and 
i2
(t
ij
)  N(0;R
i2
),where R
i2
= 
2
I
n
i
.The response vector for Equation (3.7) is:
^
Y
i
=
0
B
B
@
^
X
0
i
(t
i1
) +c
^
X
i
(t
i1
)
.
.
.
^
X
0
i
(t
in
i
) +c
^
X
i
(t
in
i
)
1
C
C
A
;i = 1;:::;m:(3.8)
We will assume the response variable is distributed as,
^
Y
i
 N(Z
i2

2
;Z
i2
D
2
Z
T
i2
+R
i2
);i = 1;:::;m:(3.9)
The variance,V
i
= Z
i2
D
2
Z
T
i2
+R
i2
is a function of  = (
2
;
11
;
12
;
22
)
T
.The fixed effects
design matrix and the randomeffects design matrix for Equation (3.7) are identical.We will
designate this matrix as Z
i2
,
Z
i2
=
0
B
B
@
1 CD4
i
(t
i1
)
.
.
.
.
.
.
1 CD4
i
(t
in
i
)
1
C
C
A
;i = 1;:::;m:(3.10)
14
The vector of fixed effects is 
2
= (a
1
;a
2
)
T
.The vector of randomeffects is b
i
= (
i1
;
i2
)
T
.
Our aimis to estimate 
2
,b
i2
,D
2
,and 
2
.The best linear unbiased estimator (BLUE) of 
2
is
(see [31]),
^

2;BLUE
=

m
X
i=1
Z
T
i2
V
i
1
Z
i2
!
1
m
X
i=1
Z
T
i2
V
i
1
^y
i
;(3.11)
where
^
Y
i
=
^
X
0
(t
ij
) +c
^
X(t
ij
).V
i
is unknown so we will replace it with an estimate,
^
V
i
,giving
us an estimated best linear unbiased estimator (EBLUE),
^

2;EBLUE
=

m
X
i=1
Z
T
i2
^
V
1
i
Z
i2
!
1
m
X
i=1
Z
T
i2
^
V
1
i
^y
i
:(3.12)
The best linear unbiased predictor (BLUP) of b
i2
is (see [31]),
^
b
i2;BLUP
= D
2
Z
T
i2
V
i
1
(^y
i
Z
i2
^

2
);i = 1;:::;m:(3.13)
As D
2
,V
i
,and 
2
are unknown we replace themwith their estimates to get an estimated best
linear unbiased predictor (EBLUP):
^
b
i2;EBLUP
=
^
D
2
Z
T
i2
^
V
1
i
(^y
i
Z
i2
^

2
);i = 1;:::;m:(3.14)
We can estimate 
2
and b
i2
through an iterative process.To obtain an initial estimate of V
i
,we
set
^
R
i2
= I
n
i
and
^
D
2
= I
2
.Thus
^
V
i
= Z
i2
Z
T
i2
+I
n
i
.Then we use these point estimates to
estimate 
2
fromEquation (3.12) and b
i2
fromEquation (3.14).Using these estimates we can
estimate 
2
2
and D
2
as follows:
^
2
2
=
1
N 2
m
X
i=1
h
^
Y
i
Z
i2

^

2
+
^
b
i2
i
T
h
^
Y
i
Z
i2

^

2
+
^
b
i2
i
(3.15)
^
D
2
=
1
m
m
X
i=1
^
b
i2
^
b
T
i2
;(3.16)
where N =
P
m
i=1
n
i
.Now we repeat,estimating first 
2
and b
i2
and then 
2
2
and D
2
until
these variables converge,as follows:
1.Estimate
^

2;EBLUE
as a function of ^
2
2
and
^
D
2
fromEquation (3.12).Then using
^

2;EBLUE
,estimate
^
b
i2;EBLUP
as a function of
^

2;EBLUE
,^
2
2
,and
^
D
2
using
Equation (3.14).
2.Estimate ^
2
2
as a function of
^

2;EBLUE
and
^
b
i2;EBLUP
fromEquation (3.15) and
^
D
2
as a
function of
^
b
i2;EBLUP
fromEquation (3.16).
This can be implemented using software such as the lme function fromthe R package,
nlme [25],or the MIXED procedure fromSAS software.
15
3.2 APPLICATION TO REAL DATA
This section describes the results of applying the methods of this chapter to the data
fromthe first armof the ACTG 398 study described in Section 1.3.
3.2.1 Estimating X
0
(t) and X(t) in a Real Data
Example
The data we used are fromthe first armof the ACTG 398 study.The measure of time
is day = (t
11
;  ;t
1n
1
;  ;t
m1
;  ;t
mn
m
)
T
.The response vector is measured
log
10
(RNA) = (Y
11
;  ;Y
1n
1
;  ;Y
m1
;  ;Y
mn
m
)
T
.We chose a grid width of 1 and
estimate values of X over the timespan of the ACTG 398 data.We estimated log
10
(RNA) for
t = 0;1;2;:::;524 using the methods developed in Section 3.1.1,as follows:
^
X
i
(t
ij
) =
^
X(t
ij
) + ^
i
(t
ij
);i = 1;  ;m;j = 1;  ;n
i
;(3.17)
^
X
i
0
(t
ij
) =
^
X
0
(t
ij
) + ^
i
0
(t
ij
);i = 1;  ;m;j = 1;  ;n
i
:(3.18)
The bandwidth chosen for the estimates was 70 days.It will be seen in Section 5.1 that the
bandwidth chosen by generalized cross-validation (GCV) was 20 days.However,70 days
resulted in a sufficiently smooth curve.Graphs of the estimates for log
10
(RNA) can be seen in
Figure 3.1.
In the plots of the estimated log
10
(RNA) one patient is a noticeable outlier,262874.
The estimate for this patient increases fromless than 2 at day 300 to 6 at day 400.The patient
at day 300 is one of the patients with the lowest viral load and then by day 400 is one of the
patients with the highest viral load.There is a corresponding increase in the plot of the first
derivative of log
10
(RNA).The outlier can easily be seen in Figure 3.1 in both the patient plots
of log
10
(RNA) and the patient plots of the first derivative of log
10
(RNA).A second set of
plots in Figure 3.2 estimates log
10
(RNA) and the first derivative of log
10
(RNA) without
patient 262874.We see that the fixed effect estimation for log
10
(RNA) in Figure 3.1 increases
fromabout 400 days to the conclusion of the study in comparison to the estimation without
the patient 262874,in Figure 3.2.However,overall the outlier does not have a strong influence
on the estimates of the fixed effect of log
10
(RNA) or the first derivative of log
10
(RNA).
Figure 3.3 shows plots of log
10
(RNA) and first derivative of log
10
(RNA) for three
selected patients.Patient 231041 is an example of a patient who responds well to the
treatment.Patient 12894 is an example of a patient who initially responds to treatment and
then experiences a rebound of viral load.Patient 630681 is an example of a patient who does
not show a strong response to the treatment.The plots of log
10
(RNA) and first derivative of
log
10
(RNA) for the outlier patient are shown in Figure 3.4.The code for estimating X
i
(t
ij
)
and dX(t
ij
)=dt is in Appendix A.
16
0
100
200
300
400
500
1
2
3
4
5
6
(a)
Day
X(t)
0
100
200
300
400
500
1
2
3
4
5
6
(b)
Day
X(t)
0
100
200
300
400
500
-0.02
0.00
0.02
0.04
(c)
Day
X'(t)
0
100
200
300
400
500
-0.02
0.00
0.02
0.04
(d)
Day
X'(t)
Figure 3.1.Estimation of log
10
(RNA) and the first derivative.(a) is the
graph of the fixed effect of
^
X
i
(t).The dotted lines indicate the 95%
confidence interval.(b) is the graph of
^
X
i
(t) for each patient.(c) is the
graph of the fixed effect of
^
X
0
i
(t).The dotted lines indicate the 95%
confidence interval.(d) is the graph of
^
X
0
i
(t) for each patient.
17
0
100
200
300
400
500
1
2
3
4
5
6
Fixed X(t)
Day
X(t)
0
100
200
300
400
500
1
2
3
4
5
6
Patient X(t)
Day
X(t)
0
100
200
300
400
500
-0.02
0.00
0.02
0.04
Fixed X'(t)
Day
X'(t)
0
100
200
300
400
500
-0.02
0.00
0.02
0.04
Patient X'(t)
Day
X'(t)
Figure 3.2.Estimation of log
10
(RNA) and the first derivative without
patient 262874.(a) is the graph of the fixed effect of
^
X
i
(t).The dotted
lines indicate the 95%confidence interval.(b) is the graph of
^
X
i
(t) for
the 34 other patients.(c) is the graph of the fixed effect of
^
X
0
i
(t).The
dotted lines indicate the 95%confidence interval.(d) is the graph of
^
X
0
i
(t) for the 34 other patients.
18
0
100
200
300
400
500
1
2
3
4
5
6
Patient 12894
Day
X(t)













0
100
200
300
400
500
-0.02
0.00
0.02
0.04
Patient 12894
Day
X'(t)













0
100
200
300
400
500
1
2
3
4
5
6
Patient 231041
Day
X(t)












0
100
200
300
400
500
-0.02
0.00
0.02
0.04
Patient 231041
Day
X'(t)












0
100
200
300
400
500
1
2
3
4
5
6
Patient 630681
Day
X(t)











0
100
200
300
400
500
-0.02
0.00
0.02
0.04
Patient 630681
Day
X'(t)











Figure 3.3.Plots of log
10
(RNA) vs.day and the first derivative of
log
10
(RNA) vs.day for representative patients.The plots to the left are
log
10
(RNA).The dots are the estimated values.The triangles are the
actual values.The dashed line is the fixed effect.The plots to the right
are the first derivative of log
10
(RNA).
19
0
100
200
300
400
500
1
2
3
4
5
6
Patient 262874
Day
X(t)











0
100
200
300
400
500
-0.02
0.00
0.02
0.04
Patient 262874
Day
X'(t)











Figure 3.4.Plots of log
10
(RNA) vs.day and the first derivative of
log
10
(RNA) vs.day for patient outlier.The plot to the left is
log
10
(RNA).The dots are the estimated values.The triangles are the
actual values.The dashed line is the fixed effect.The plot to the right
is the first derivative of log
10
(RNA).
3.2.2 Estimating a
1
and a
2
in a Real Data Example
We use the R function,lme,fromthe R package,nlme [25].The viral clearance rate,
c,was estimated with a range of values.The estimates of the fixed effects and randomeffects
can be seen in Table 3.1.Larger values of c resulted in larger absolute value of the estimates of
the coefficients.The standard error for the fixed effect estimate of a
1
is small in proportion to
the size of the estimate,while the standard error of the fixed effect estimate of a
2
is relatively
large in proportion to the size of the estimate.Table 3.2 lists the fixed effects estimated with
all patients and also the fixed effects estimated without the outlier,patient 262874.When
comparing the estimates using all patients to the estimates without patient 262874,we see that
a
1
is not strongly influenced by the outlier.For a
2
,we see that there is a proportionally larger
influence on the value,the standard error of the fixed effect,and the standard deviation of the
randomeffect compared to a
1
.The code for Section 3.2.2 is in Appendix B.
20
Table 3.1.Estimates of Time-Invariant Coefficients
Parameters Clearance Rate
c = 2.8 c = 3.0 c = 3.2 c = 3.5
Fixed Effects a
1
11.201 12.001 12.802 14.003
SE(a
1
) 0.589 0.631 0.673 0.736
a
2
-0.0039 -0.0042 -0.0044 -0.0049
SE(a
2
) 0.0010 0.0011 0.0012 0.0013
MS 1.60 1.71 1.82 1.99
Patient
12894 a
1
+ 
i1
8.554 9.165 9.777 10.694
a
2
+ 
i2
0.00041 0.00044 0.00047 0.00051
SE 1.45 1.55 1.65 1.81
60943 a
1
+ 
i1
8.212 8.799 9.387 10.268
a
2
+ 
i2
-0.00024 -0.00026 -0.00028 -0.00031
SE 1.41 1.52 1.62 1.77
82656 a
1
+ 
i1
10.961 11.744 12.527 13.701
a
2
+ 
i2
-0.00300 -0.00322 -0.00343 -0.00376
SE 1.79 1.92 2.04 2.24
150556 a
1
+ 
i1
13.727 14.707 15.687 17.157
a
2
+ 
i2
0.00044 0.00047 0.00050 0.00055
SE 1.68 1.80 1.92 2.10
181754 a
1
+ 
i1
12.939 13.864 14.788 16.174
a
2
+ 
i2
-0.00263 -0.00282 -0.00301 -0.00329
SE 0.24 0.25 0.27 0.30
220802 a
1
+ 
i1
9.303 9.969 10.635 11.633
a
2
+ 
i2
-0.00601 -0.00644 -0.00687 -0.00752
SE 0.77 0.82 0.88 0.96
231041 a
1
+ 
i1
13.999 14.998 15.998 17.497
a
2
+ 
i2
-0.00317 -0.00340 -0.00363 -0.00397
SE 2.23 2.39 2.55 2.79
242673 a
1
+ 
i1
14.422 15.452 16.483 18.028
a
2
+ 
i2
-0.00337 -0.00361 -0.00386 -0.00422
SE 1.53 1.64 1.75 1.92
580899 a
1
+ 
i1
8.717 9.341 9.965 10.901
a
2
+ 
i2
-0.01182 -0.01267 -0.01351 -0.01478
SE 1.96 2.10 2.24 2.45
630681 a
1
+ 
i1
14.000 15.000 16.000 17.500
a
2
+ 
i2
-0.00269 -0.00289 -0.00308 -0.00337
SE 1.03 1.11 1.18 1.29
21
Table 3.2.Time-Invariant ^a
1
and ^a
2
Estimated With and Without
Outlier
Parameters Clearance Rate
c = 2.8 c = 3.0 c = 3.2 c = 3.5
All Samples a
1
11.201 12.001 12.802 14.003
SE(a
1
) 0.589 0.631 0.673 0.736
a
2
-0.0039 -0.0042 -0.0044 -0.0049
SE(a
2
) 0.0010 0.0011 0.0012 0.0013
MS 1.60 1.71 1.82 1.99
Without 262874 a
1
11.0459 11.835 12.625 13.808
SE(a
1
) 0.602 0.645 0.688 0.752
a
2
-0.0027 -0.0029 -0.0031 -0.0034
SE(a
2
) 0.00078 0.00083 0.00089 0.00097
MS 1.43 1.53 1.64 1.79
3.3 SIMULATION FOR A SIMPLIFIED MODEL
WITH CONSTANT a
1
Simulated values for log
10
(RNA) will be generated following the methods in [2].In
these equations,a
1
is modeled.CD4 and a
2
are omitted for simplicity.The equation used is:
d
dt
X(t) = cX(t) +a
1
;(3.19)
where c is set to 3.5 (taken from[3]).This equation has the solution:
X(t) = e
ct

X(0) +
Z
0
t
e
cs
a
1
(s)ds

;(3.20)
where X(0) = 1:0.The time-invariant,mixed-effects formof (3.19) is:
dX
i
(t)
dt
= cX
i
(t) +a
1
+
i1
:(3.21)
with the solution:
X
i
(t) = e
ct
X
i
(0) +
a
1
+
i1
c
(1 e
ct
);i = 1;:::;m (3.22)
We set the covariance structure of the randomeffects of X
i
(t) to be the first-order
autoregressive structure (AR(1)).For this simulation we will use a balanced design.The
equations for X
i
(t) are thus,
Y
i
(t) = X(t) +
i
(t) +
i
(t)X
i
(t) = X(t) +
i
(t) (3.23)
22
where 
i1
 N(0;
3
),
i
 N(0;
1
),
i
(t)  N(0;),and
 = 
2
2
0
B
B
B
B
@
1     
n1
 1    
n2
.
.
.
.
.
.
.
.
.
.
.
.

n1

n2
   1
1
C
C
C
C
A
;i = 1;:::;m:
Putting together Equation 3.22 and Equation 3.23 we obtain an equation to simulate points for
X(t),
X
i
(t) = e
ct
X
i
(0) +
a
1
+
i1
c
(1 e
ct
) +
i
(t) +
i
(t);
i = 1;:::;m:
(3.24)
The number of patients (m),number of time points per patient (n),and data quality (
1
2
,
2
2
,
and 
3
2
) may affect the accuracy of the estimates of X(t),dX(t)=dt,and a
1
.We simulated
several versions of the model by varying these variables.The constants were set as X(0) = 1,
c = 3:5,a
1
= 2:0,and  = 0:45.The value for 
1
,
2
,and 
3
where the data quality is good is
0:01.Where data quality is noisy,the value is 0:02.The resulting estimates for the models can
be seen in Table 3.3.The estimate of a
1
is very good.For most of the models it was slightly
biased upward.The table also shows SE(^a
1
) and
^
SE(a
1
) where,
SE(^a
1
) =
s
P
K
k=1
(^a
1k
^a
1
)
2
K 1
^a
1
=
P
K
k=1
^a
1k
K
^
SE(a
1
) =
P
K
k=1
SE(a
1k
)
K
;
where K is the number of iterations.The 95%CI coverage is the ratio of the K simulations
where a
1
falls within the 95%confidence interval derived from^a
1k
and the SE(^a
1
).The
standard error is taken fromthe linear mixed-effects solution to Step 2.This standard error
assumes that the dependent variable,
^
X
0
(t) +c
^
X(t) is exact,when in fact it has error.As a
result of this,the confidence intervals are sometimes anti-conservative.Figure 3.5 shows the
simulated X(t) for representative simulations fromthree separate models.Figure 3.6 shows
the estimated fixed X(t) and the estimated patient X(t) for the same three models.The
confidence intervals for fixed X(t) can be seen to be very narrow,even with noisy data or
small n and m.Figure 3.7 shows the estimated fixed first derivative of X(t) and the estimated
patient first derivative of X(t) for the same three models.Again,we see the confidence
intervals for fixed X
0
(t) are narrow,even with noisy data or small n and m.The code for the
simulation is in Appendix C.
23
Table 3.3.The Estimates for a
1
fromSimulation Models for 200
Simulations
Model:(m;n;)
h
Mean ^a
1
SE(^a
1
)
^
SE(a
1
)
95%CI Coverage
10,15,0.01
0.07
2.02
0.028
0.024
0.92
10,15,0.02
0.10
2.02
0.055
0.049
0.93
40,15,0.01
0.06
2.01
0.013
0.012
0.93
40,15,0.02
0.08
2.02
0.026
0.024
0.91
80,15,0.01
0.06
2.00
0.010
0.009
0.95
80,15,0.02
0.07
2.01
0.021
0.018
0.94
10,50,0.01
0.05
2.01
0.019
0.019
0.93
10,50,0.02
0.08
2.03
0.036
0.038
0.86
40,50,0.01
0.05
2.01
0.010
0.010
0.90
40,50,0.02
0.06
2.01
0.022
0.019
0.94
80,50,0.01
0.05
2.01
0.008
0.007
0.91
80,50,0.02
0.05
2.01
0.015
0.014
0.91
10,100,0.01
0.04
2.01
0.015
0.016
0.93
10,100,0.02
0.05
2.01
0.032
0.031
0.94
40,100,0.01
0.04
2.01
0.008
0.008
0.92
40,100,0.02
0.05
2.01
0.008
0.008
0.78
80,100,0.01
0.04
2.01
0.008
0.006
0.83
80,100,0.02
0.05
2.01
0.012
0.011
0.91
24






















































































































































0.0
0.2
0.4
0.6
0.8
0.4
0.6
0.8
1.0
1.2
(a)
Time
X(t)
























































































































































































































































































































































































































































































































































































































0.0
0.2
0.4
0.6
0.8
0.4
0.6
0.8
1.0
1.2
(b)
Time
X(t)