# Phenobarbital recursive form Example
#
# INTRODUCTION
# This script includes the examples below.
# 1) PK modeling of multiple dose data using superposition
# 2) specification of covariates
# 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 PHENOBARBITAL DATA
# The object Phenobarb 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
#
# 1-compartment open model with IV administration and 1st order elimination
trellis.device()
plot(Phenobarb)
options(width=68, digits=5)
#
# THE MODEL
# The C function 'nlme_one_comp_first' (provided with S-Plus 2000, nlme3.3)
# implements a 1-compartment open model with IV drug administration and 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. V and Cl are reparameterized as the log(V) and the 
# log(Cl) to obtain positive estimates without constraints.
# The C function is called by the S-Plus phenoModel 
# function. The model is not self-starting, so initial estimates are required.
# The phenoModel function can be printed to the screen by typing phenoModel at the command line
phenoModel
#
# The C source code for 'nlme_one_comp_first' is available on the website. 
#
# Print one subject's data, note the format, the data consists of dose records or conc records
Phenobarb[Phenobarb$Subject==3,]
#
# 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 records.
# Create an nlme object called Pheno1.nlme with the model and results
# Assume a diagonal to avoid convergence problems due to the sparse data
Pheno1.nlme <- nlme(conc~phenoModel(Subject,time,dose,lCl,lV),data=Phenobarb,
fixed=lCl+lV~1,random=pdDiag(lCl+lV~1),start=c(-5,0),
na.action=na.include,naPattern=~!is.na(conc))

# Apply nlme methods to plot and list the results. 
# plot random effects vs Wt and ApgarInd and add a loess smoother for the continuous covars
# note both Cl and V increase with birth weight
Pheno1.ranef <- ranef(Pheno1.nlme,augFrame=T)
plot(Pheno1.ranef, form=lCl~Wt+ApgarInd)
plot(Pheno1.ranef, form=lV~Wt+ApgarInd)

# ADD COVARIATES TO THE MODEL
# use a linear model for increase in V and Cl with Wt
# 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"))
Pheno2.nlme <- update(Pheno1.nlme,
fixed=list(lCl~Wt,lV~Wt+ApgarInd),
start=c(-5.0935,0,0.34259,0,0),
control=list(pnlsTol=1e-6))
#reduce the pnlsTol tolerance to prevent a convergence problem in PNLS step
summary(Pheno1.nlme)
fixef(Pheno2.nlme)
fixef(Pheno1.nlme)

# drop ApgarInd from the model due to lack of significance 
Pheno3.nlme <- update(Pheno2.nlme,
fixed=lCl+lV~Wt,start=fixef(Pheno2.nlme)[-5])
summary(Pheno3.nlme)
plot(fm3Pheno.nlme,conc~fitted(.),abline=c(0,1))