# Theophylline Self-Start Example
#
# INTRODUCTION
# This script includes the examples below.
# 1) naive pooled analysis using nls()
# 2) standard two-stage analysis using nlsList()
# 3) population mixed effects analysis using nlme()
# 4) modifying a population model with the update function
# 5) comparing population models
# 6) modeling heteroscedastic error
#
# A self-starting nonlinear regression model is used,
# eliminating the need to provide initial estimates.
# 
# 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 THEOPHYLLINE DATA
# The object Theoph 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(Theoph)
plot(Theoph,outer=~1,inner=~Subject,aspect="fill")
#
# NAIVE POOLED ANALYSIS
#
# The function SSfol is a self-starting function for
# a one compartment, first order absorption, first order
# elimination model. Initial estimates are obtained by
# curve peeling. The model is parametrized in terms of the
# logarithms of Cl and 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 SSfol in the Search dialogue.
# You can print the function by typing SSfol at the command line.
# 
# SSfol
# selfStart(~Dose * exp(lKe + lKa - lCl) * (exp(-exp(lKe) * input) -
# exp(-exp(lKa) * input))/(exp(lKa) - exp(lKe)), 
# function(Dose, input, lKe, lKa, lCl)
# Create an nls object named 'Theonaive.nls', the output
# object includes the model and results.  
Theonaive.nls <- nls(conc~SSfol(Dose, Time, lKe, lKa, lCl),
data=Theoph)

# Apply nls methods to plot and list the results. 
plot(Theonaive.nls)
# note the biased results using naive pooled analysis
plot(Theonaive.nls,Subject~resid(.),abline=0)
coef(Theonaive.nls)
summary(Theonaive.nls)

#
# STANDARD TWO-STAGE ANALYSIS
#
# Create an nlsList object named 'TheoSTS.nls', the output
# object includes the model and results.  
TheoSTS.nls <- nlsList(conc~SSfol(Dose, Time, lKe, lKa, lCl),
data=Theoph)
# Apply nls methods to plot and list the results. 
# plot the residuals
plot(TheoSTS.nls)
# note less bias as compared to the results using naive pooled analysis
plot(TheoSTS.nls,Subject~resid(.),abline=0)
qqnorm(TheoSTS.nls)
# predictions augmented with observed values
plot(augPred(TheoSTS.nls))
# deviations of coefficients from average
plot(ranef(TheoSTS.nls))
# 95% CI on the coefficients
plot(intervals(TheoSTS.nls))
# scatter-plot matrix of coefficients or random effects
pairs(TheoSTS.nls,id=.15)
pairs(coef(TheoSTS.nls))

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

#
# POPULATION MIXED EFFECTS ANALYSIS
#
# Create an nlme object named 'Theo.nlme', the output
# object includes the model and results.
# Three fixed effects and 3 random effects are specified,
# with no intercept (the ~ 1)  
Theo.nlme <- nlme(conc ~ SSfol(Dose, Time, lKe, lKa, lCl),
	fixed = lKe + lKa + lCl ~ 1,
	random = lKe + lKa + lCl ~ 1,
	data=Theoph)
# Alternatively, use the nlsList standard-two stage as the
# model
Theo.nlme <- nlme(TheoSTS.nls,
	fixed = lKe + lKa + lCl ~ 1,
	random = lKe + lKa + lCl ~ 1,
	data=Theoph)

# Apply nlme methods to plot and list the results. 
# plot the residuals
plot(Theo.nlme)
plot(Theo.nlme,Subject~resid(.),abline=0)
qqnorm(Theo.nlme, abline=c(0,1))
# predictions augmented with observed values
plot(augPred(Theo.nlme))
# random effects estimates
plot(ranef(Theo.nlme))
qqnorm(Theo.nlme, ~ranef(.))
# random effects covariance
pairs(Theo.nlme,id=.1)

# 95% CI on the coefficients
intervals(Theo.nlme)
# deviations of coefficients from average
ranef(Theo.nlme)
# coefficients from individual nls fits
coef(Theo.nlme)
# fitted values from individual nls fits
fitted(Theo.nlme)
# print details 
summary(Theo.nlme)

#COMPARE POPULATION AND STS ANALYSIS
plot(compareFits(coef(Theo.nlme), coef(TheoSTS.nls)))
pairs(compareFits(coef(Theo.nlme), coef(TheoSTS.nls)))

# UPDATE THE POPULATION MODEL
#update the model fit and specify random effects for only Ka and Cl
#name the model Theo2.nlme
Theo2.nlme <- update(Theo.nlme, random = lKa + lCl ~ 1)
	
plot(Theo2.nlme)
plot(augPred(Theo2.nlme))
plot(ranef(Theo2.nlme))

# COMPARE THE MODELS
anova(Theo.nlme,Theo2.nlme)
plot(compareFits(coef(Theo.nlme), coef(Theo2.nlme)))
pairs(compareFits(coef(Theo.nlme), coef(Theo2.nlme)))

# HETEROSCEDASTIC ERROR PATTERN
# The within-group variance increases with the concentration of theophylline
# A power variance function can't be used if at time = 0 (where conc = 0)
# are retained in the data because the weights are undefined at t=0.
# The constant plus power function, varConstPower, can be used. A constant
# is added to the power of the fitted value.   
Theo3.nlme <- update(Theo2.nlme, weights=varConstPower(power=0.1))
plot(Theo3.nlme)