# Quinidine recursive form Example
#
# INTRODUCTION
# This script includes the examples below.
# 1) PK modeling of multiple dose data using superposition and a steady state equation
# 2) specification of covariates
# 3) modeling heteroscedastic error
# The examples also illustrate general functions relevant
# to PK modeling.
#
# REFERENCE 
# 'Mixed Effects Modeling Using S and S-Plus' by Pinheiro and Bates
# Please find the link to the book the Examples page and the Resources page.
# The example is worked in detail with a more theoretical presentation
# in the book.
#
# RUNNING THE SCRIPT
# Please refer to the 'GENERAL INSTRUCTIONS' link
# on the website. The data and functions are available
# with S-Plus 2000 for Windows. 
#  
#
# PLOT THE QUINIDINE DATA
# The object Quinidine is a groupedData object
# Please refer to the groupedData example for more information
# on how to create and use groupedData objects.
# Plots can be quickly generated using the 'plot' method
#
trellis.device()
plot(Quinidine,layout=c(3,7))
#
# THE MODEL
# The C function 'nlme_one_comp_open' (provided with S-Plus 2000, nlme3.3)
# implements a 1-compartment open model with first order absorption, first
# order elimination that can handle multiple dose data. The model is expressed
# in recursive form and sums the individual contributions of each dose by reading
# the dose records from the data. The data also includes an 'interval' variable for
# conc measurements taken after adminstration of the same dose at regular intervals
# (interval = 8, one dose received every 8 hours). The model converges to a steady state
# model after long trm repeated dosing at regular intervals. The data set includes either a
# dose record (conc and interval missing (NA)), a conc record (dose and interval missing),
# or an interval-dose record (steady-state, conc missing). The C program uses the steady state
# model if interval and dose are not missing.
# V, ka, and Cl are reparameterized as log(ka), log(V) 
# and log(Cl) to obtain positive estimates without constraints.
# The C function is called by the S-Plus quinModel 
# function. The model is not self-starting, so initial estimates are required.
# The quinModel function can be printed to the screen by typing quinModel at the command line
quinModel
#
# The C source code for 'nlme_one_comp_open' is available on the website. 
#
# Print one subject's data, note the format, the data consists of dose records or conc records or
# dose and interval (steady state) records.
Quinidine[Quinidine[["Subject"]] == 3, 1:8]
#
# The nlme option na.action is used to include rows with a dose 
# and naPattern is used to remove these records for calculation
# of the objective function because of the NA (missing) for conc
# Create an nlme object called Quinidine1.nlme with the model and results
# Assume a diagonal to avoid convergence problems due to the sparse data
options(width=68, digits=5)
Quinidine1.nlme <-
nlme(conc ~ quinModel(Subject, time, conc, dose, interval, lV, lKa, lCl),
data=Quinidine, fixed = lV + lKa + lCl ~ 1,
random = pdDiag(lV + lCl ~ 1), groups = ~ Subject,
start = list(fixed = c(5, -0.3, 2)),
na.action = na.include, naPattern = ~ !is.na(conc) )
summary(Quinidine1.nlme)

# EXAMINE COVARIATES
# extract the random effects using ranef() and plot the random effects and
# covariates.
# note Cl decreases with glycoprotien and age and
# Cl increases with CrCl and Wt
Quinidine1.nlmeRE <- ranef(Quinidine1.nlme,aug=T)
Quinidine1.nlmeRE[1:3,]
plot(Quinidine1.nlmeRE, form = lCl ~ Age + Smoke + Ethanol + Weight + Race + 
Height + glyco + Creatinine + Heart,
control = list(cex.axis = 0.7))

# ADD COVARIATES TO THE MODEL USING A STEPWISE MODEL-BUILDING APPROACH
#
# Add glycoprotein to the model, glycoprotein can change with time for a subject, 
# so index by patient and time
# Extract the fixed effects from the base model to use as initial estimates
Quinidine1.fix <- fixef(Quinidine1.nlme)
Quinidine2.nlme <- update(Quinidine1.nlme, fixed = list(lCl ~ glyco, lKa + lV ~ 1),
start = c(Quinidine1.fix[3], 0, Quinidine1.fix[2:1]))
# Summarize the results and note glycoprotein is significant
summary(Quinidine2.nlme)

# Add creatinine to the model
# use treatment contrasts, refer to the online help and S-Plus manual for more
# information on the contrast options
options(contrasts=c("contr.treatment", "contr.poly"))
# Extract the fixed effects from the previous model to use as initial estimates
Quinidine2.fix <- fixef(Quinidine2.nlme)
Quinidine3.nlme <- update(Quinidine2.nlme, fixed = list(lCl ~ glyco + Creatinine, lKa + lV ~ 1),
start = c(Quinidine2.fix[1:2], 0.2, Quinidine2.fix[3:4]))
summary(Quinidine3.nlme)

# add weight to the model
Quinidine3.fix <- fixef(Quinidine3.nlme)
Quinidine4.nlme <- update(Quinidine3.nlme, fixed = list(lCl ~ glyco + Creatinine + Weight, lKa + lV ~ 1),
start = c(Quinidine3.fix[1:3], 0, Quinidine3.fix[4:5]))
summary(Quinidine4.nlme)

# VARIANCE OF THE WITHIN-GROUP ERROR
# plot the residuals and note a heteroscedastic pattern, the within
# group errors increase with quinidine conc
# use the varPower class to calculate the weights and model the within-group
# heteroscedastic error
# The conc in the data are not so close to zero that the power variance function
# can't be used
plot(Quinidine4.nlme, xlim=c(0,6.2))
Quinidine5.nlme <- update(Quinidine4.nlme, weights = varPower() )
summary(Quinidine5.nlme)

# compare the models
anova(Quinidine4.nlme,Quinidine5.nlme)
plot(Quinidine5.nlme, xlim=c(0,8.2))