# Indomethicin Self Start Example
#
# INTRODUCTION
# This script includes the examples below.
# 1) PK modeling using a self-start nonlinear regression model
# a) standard two-stage analysis using nlsList()
# b) population mixed effects analysis using nlme()
# 2) PK modeling using a custom model function
# a) population mixed effects analysis using nlme()
# The examples also illustrate general functions relevant
# to PK modeling and introduce users to nlme() random effect options
# 
# 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 INDOMETHICIN DATA
# The object Indometh 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(Indometh)
plot(Indometh,outer=~1,inner=~Subject,aspect="fill")

options(width=68, digits=5)

# The function SSbiexp is a self-starting function for
# a biexponential, a linear combination of two negative
# exponential terms. Initial estimates are obtained by
# curve peeling. The model is parametrized in terms of the
# logarithms of the rate constants to ensure
# positive estimates while keeping the optimization
# unconstrained. 
# The function can be used for full profile data.
# More information can be found in the online help
# by entering SSbiexp in the Search dialogue or ?SSbiexp at the command line.
# You can print the function by typing SSbiexp at the command line.
# 
# SSbiexp - linear combo of two negative exponential terms
#  selfStart(~ A1 * exp(-exp(lrc1)*input) + A2 * exp(-exp(lrc2) * input), 
# function(input, A1, lrc1, A2, lrc2)
# 
# STANDARD TWO STAGE
# Create an nlsList object named 'IndoSTS.lis', the output
# object includes the model and results.  
IndoSTS.lis <- nlsList(conc~SSbiexp(time,A1,lrc1,A2,lrc2),data=Indometh)

# Apply nlsList methods to plot and list the results. 
# boxplot of residuals by subject
plot(IndoSTS.lis,Subject~resid(.),abline=0)
# 95% CI on biexp parameters for each individual to chart the variability 
plot(intervals(IndoSTS.lis))
#check normality assumption
qqnorm(IndoSTS.lis,abline=c(0,1))
# predictions augmented with observed values
plot(augPred(IndoSTS.lis))
# deviations of coefficients from average
plot(ranef(IndoSTS.lis))
# 95% CI on the coefficients
plot(intervals(IndoSTS.lis))
# scatter-plot matrix of coefficients or random effects
# use the id statement to identify any extremes 
pairs(IndoSTS.lis,id=.1)
pairs(coef(IndoSTS.lis))

# deviations of coefficients from average
ranef(IndoSTS.lis)
# coefficients from individual nls fits
coef(IndoSTS.lis)
# fitted values from individual nls fits
fitted(IndoSTS.lis)
summary(IndoSTS.lis)

#
# POPULATION MIXED EFFECTS ANALYSIS USING A SELF-START FUNCTION
#
# Create an nlme object named 'Indo.nlme', the output
# object includes the model and results.
# include random effects for all parameters and
# use a diagonal initially to prevent convergence problems from an overparameterized model,
# due to the large number of random parameters relative to number of individuals (24 vs 6)
# use the nlsList object to specify the model to fit, parameters to estimate,
# and starting estimates for fixed effects
Indo.nlme <- nlme(IndoSTS.lis, 
random=pdDiag(A1+lrc1+A2+lrc2~1))
# print the results
summary(Indo.nlme)
# pdDiag is from the class of pdMat - positive definite matrices
# Type ?pdMat at the command line for more information
?pdMat
# use a general positive definite matrix with 3 random effects
Indo2.nlme <- update(Indo.nlme, random=A1+lrc1+A2~1)
# use pairs() to evaluate correlations, A1 and lrc1 random effects are correlated
pairs(Indo2.nlme)

#use a blocked diagonal matrix to represent var-cover structure of random effects
#based on correlation between random effects
Indo3.nlme <- update(Indo2.nlme,
random = pdBlocked(list(A1+lrc1~1,A2~1)))

#compare models 
anova(Indo.nlme, Indo3.nlme)
anova(Indo2.nlme, Indo3.nlme)

# evaluate fit 
# use id statement to identify most outlying obs, 2 subjects labeled
plot(Indo3.nlme, id=0.05, adj=-1)
qqnorm(Indo3.nlme, abline=c(0,1))
plot(augPred(Indo3.nlme,level=0:1))
summary(Indo3.nlme)

#compare nlme and nlsList (STS)
plot(compareFits(coef(Indo3.nlme), coef(IndoSTS.lis)))

#
# POPULATION MIXED EFFECTS ANALYSIS USING A MODEL FUNCTION
#
# write a model function using the same formula as SSbiexp
#  selfStart model (~ A1 * exp(-exp(lrc1)*input) + A2 * exp(-exp(lrc2) * input), 
#
# run the model definition section before using the model in nlme
logtwoCbolus.func <- function(time, A, alpha, B, beta)
A*exp(-exp(alpha)*time) + B*exp(-exp(beta)*time)

# use nlme with the function and initial starting values
initial1Indom.nlme <- nlme(model = conc ~ logtwoCbolus.func(time, A, alpha, B, beta),
fixed = A + alpha + B + beta ~ 1, random = pdBlocked(list(A+alpha~1,B~1)), data=Indometh,
start=c(2.3,0.8,0.6,-1.08))
summary(initial1Indom.nlme)
coef(initial1Indom.nlme)
fixef(initial1Indom.nlme)

# write a model function using the formula below
# 2C-bolus 
# y = A * exp(-Alpha * t) + B * exp(-Beta * t)
#
twoCbolus.func <- function(time, A, alpha, B, beta)
A*exp(-alpha*time) + B*exp(-beta*time)
# nlme
initial2Indom.nlme <- nlme(model = conc ~ twoCbolus.func(time, A, alpha, B, beta),
fixed = A + alpha + B + beta ~ 1, 
random=pdDiag(A+alpha+B~1),
data=Indometh,
start=c(2.3, 2, 0.6, .3))
initial2Indom.nlme
coef(initial2Indom.nlme)
fixef(initial2Indom.nlme)

# compare the fixed effects estimates using the different parametrization 
fixef(initial2Indom.nlme)
fixef(initial1Indom.nlme)
alpha <- exp(0.876)
alpha
beta <- exp(-1.28)
beta