EXTENDINGLONGITUDINAL ANALYSIS

TOESTIMATION OF COEFFICIENTS IN HIV-1 DYNAMICS

A Thesis

Presented to the

Faculty of

San Diego State University

In Partial Fulﬁllment

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 Coefﬁcients 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 coefﬁcients of a linear differential equation that models

the rate of HIV-1 RNA concentration in the patients’ plasma.Two methods are developed.In

the ﬁrst,the coefﬁcients are assumed to be constant.In the second,the coefﬁcients are

allowed to vary with time.For both methods the ﬁrst step is to estimate the viral load and the

ﬁrst 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 ﬁrst time derivative of viral

load are then used in estimating the coefﬁcients.Where the the coefﬁcients are assumed

constant,the second step estimates these coefﬁcients with a linear mixed-effects equation.

Where the coefﬁcients are assumed to be time-varying,a local polynomial kernel estimation

with mixed-effects is used to estimate the coefﬁcients.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 Simpliﬁed 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 Simpliﬁed 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 Coefﬁcients...........................................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 ﬁrst derivative................................16

Figure 3.2 Estimation of log

10

(RNA) and the ﬁrst derivative without patient 262874......17

Figure 3.3 Plots of log

10

(RNA) vs.day and the ﬁrst derivative of log

10

(RNA)

vs.day for representative patients.....................................................18

Figure 3.4 Plots of log

10

(RNA) vs.day and the ﬁrst 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 ﬁrst 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 ﬁrst 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 coefﬁcients

of two dynamic models for the rate of the viral load of Human Immunodeﬁciency 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 ﬁnd 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 identiﬁed 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 conﬁned to

West Africa [14].When HIV-1 was ﬁrst identiﬁed there was great optimismthat an effective

vaccine would be discovered within two years;unfortunately,development of a vaccine has

proved remarkably difﬁcult [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 Deﬁciency Syndrome

(AIDS) [8].The ﬁrst drug successfully used for AIDS treatment was azidothymidine (AZT)

which was licensed to treat AIDS in 1987 [11].Treatment with AZT is initially beneﬁcial.

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 signiﬁcant 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 inﬂuenced 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 coefﬁcients,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 ﬁrst estimates the

coefﬁcients 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 ﬁrst 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 ﬁtting 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 ﬁt 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 speciﬁed 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 ﬁve and twelve sample collections.

The largest number of collections for a patient was ﬁfteen.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 ﬁrst 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 ﬁxed 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 ﬁrst 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

coefﬁcients 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) = BT

(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 efﬁcacy; 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) = BT

(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 coefﬁcients.Li,et al.[16] proposed a method using

cubic splines to estimate time-varying coefﬁcients of a set of deterministic dynamic equations.

Liang and Wu [18] used local smoothing and a pseudo-least squares approach to estimate

constant coefﬁcients 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 efﬁcient 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 ﬁrst 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 deﬁned 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 ﬁxed effects,b

i

is the vector of random

effects,X

i

is the ﬁxed effects design matrix,Z

i

is the randomeffects design matrix,and

i

is

the error.

Shi and Taylor [27] develop a method of ﬁtting linear combinations of B-splines to

both the ﬁxed and randompart of unbalanced longitudinal data in which the randompart is

ﬁtted to each patient.Rice and Wu [26] also presented a spline ﬁtting 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 ﬁxed and randomeffects for any t are

estimated using a ﬁrst 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 ﬁxed and

randomeffects for any t are estimated using a ﬁrst 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 ﬁxed 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 ﬁrst 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 coefﬁcients,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 ﬁxed 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 ﬁxed 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 ﬁxed 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 ﬁrst

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 ﬁrst 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 ﬁrst 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 sufﬁciently 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 ﬁrst

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 ﬁrst derivative of log

10

(RNA).A second set of

plots in Figure 3.2 estimates log

10

(RNA) and the ﬁrst derivative of log

10

(RNA) without

patient 262874.We see that the ﬁxed 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 inﬂuence

on the estimates of the ﬁxed effect of log

10

(RNA) or the ﬁrst derivative of log

10

(RNA).

Figure 3.3 shows plots of log

10

(RNA) and ﬁrst 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 ﬁrst 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 ﬁrst derivative.(a) is the

graph of the ﬁxed effect of

^

X

i

(t).The dotted lines indicate the 95%

conﬁdence interval.(b) is the graph of

^

X

i

(t) for each patient.(c) is the

graph of the ﬁxed effect of

^

X

0

i

(t).The dotted lines indicate the 95%

conﬁdence 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 ﬁrst derivative without

patient 262874.(a) is the graph of the ﬁxed effect of

^

X

i

(t).The dotted

lines indicate the 95%conﬁdence interval.(b) is the graph of

^

X

i

(t) for

the 34 other patients.(c) is the graph of the ﬁxed effect of

^

X

0

i

(t).The

dotted lines indicate the 95%conﬁdence 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 ﬁrst 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 ﬁxed effect.The plots to the right

are the ﬁrst 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 ﬁrst 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 ﬁxed effect.The plot to the right

is the ﬁrst 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 ﬁxed effects and randomeffects

can be seen in Table 3.1.Larger values of c resulted in larger absolute value of the estimates of

the coefﬁcients.The standard error for the ﬁxed effect estimate of a

1

is small in proportion to

the size of the estimate,while the standard error of the ﬁxed effect estimate of a

2

is relatively

large in proportion to the size of the estimate.Table 3.2 lists the ﬁxed effects estimated with

all patients and also the ﬁxed 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 inﬂuenced by the outlier.For a

2

,we see that there is a proportionally larger

inﬂuence on the value,the standard error of the ﬁxed 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 Coefﬁcients

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 ﬁrst-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

n1

1

n2

.

.

.

.

.

.

.

.

.

.

.

.

n1

n2

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%conﬁdence 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 conﬁdence 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 ﬁxed X(t) and the estimated patient X(t) for the same three models.The

conﬁdence intervals for ﬁxed X(t) can be seen to be very narrow,even with noisy data or

small n and m.Figure 3.7 shows the estimated ﬁxed ﬁrst derivative of X(t) and the estimated

patient ﬁrst derivative of X(t) for the same three models.Again,we see the conﬁdence

intervals for ﬁxed 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)

## Comments 0

Log in to post a comment