SPLUS CODE:
######################### Beginning of Splus Code ########################
# Input Data: assume that the data at the end is stored in an ASCII file
# named "data".
workd<-read.table(file="data",header=T,row.names=NULL)

# The measurements before treatment initiation are not used in the analysis.
workd<-workd[workd$Day>0,]

# Impute below detectable viral measurement by 50.
workd$RNA[workd$RNA==100]<-50

# Define functions representing Model (12) and Model (13)
# of Wu and Ding (1999).
exp.model12<-function(p0, p1, d1, X)
{
P0 <- exp(p0)
P1 <- exp(p1)
P0 + P1 * exp( - d1 * X)
}

exp.model13<-function(p1, d1, p2, d2, X)
{
P1 <- exp(p1)
P2 <- exp(p2)
P1 * exp( - d1 * X) + P2 * exp( - d2 * X)
}

# Fit Model (12) and Model (13) by NLME
actg315.model12<-nlme(Y ~ log10(exp.model12(p0, p1, d1, X)),
fixed = list(p0~., p1~., d1~.),
random = list(p0~., p1~., d1~.),
cluster = ~ Z,
data = data.frame(Y=log10(workd$RNA), X=workd$Day, Z=workd$ID),
start = list(fixed= c(5,11,0.5)),verbose=T)

actg315.model13<-nlme(Y ~ log10(exp.model13(p1, d1, p2, d2, X)),
fixed = list(p1~., d1~., p2~., d2~.),
random = list(p1~., d1~., p2~., d2~.),
cluster = ~ Z,
data = data.frame(Y=log10(workd$RNA), X=workd$Day, Z=workd$ID),
start = list(fixed= c(12,0.4,7,0.03)),verbose=T)

anova(actg315.model12,actg315.model13)
###################### End of Splus Code ######################################

RESULTS:
--------

Model Df AIC BIC Loglik Test Lik.Ratio P value
actg315.model12 1 10 433.03 470.69 -206.52
actg315.model13 2 15 257.98 314.46 -113.99 1 vs. 2 185.05 0

Fixed Effects Estimates from Model (12):
LogP0 LogP1 delta_p
5.666933 11.10344 0.2246409

Fixed Effects Estimates from Model (13):
LogP1 delta_p LogP2 lamba_l
12.32537 0.4759455 7.945189 0.04134005

**********************************************************************************