############################################################
#                                                          #
#  Fit nonlinear  effects model with intra-plot power      #
#  variance function using nlme                            #
#                                                          #
############################################################

library(nlme)

############################################################
#                                                          #
#  Read in the data from the data frame                    #
#                                                          #
############################################################

thedata <- scan("argconc.dat")
thedata <- matrix(thedata,ncol=5,byrow=T);

indiv <- thedata[,2]   # individual indicator
ds <- thedata[,3]      # dose
xs <- thedata[,4]      # time
ys <- thedata[,5]      # concentration

sdat <- data.frame(indiv,xs,ys,ds)
	
############################################################
#                                                          #
#  Define the mean response function.  Here, we program    #
#  both the function and its derivatives with respect to   #
#  each parameter.  If analytic derivatives are not        #
#  provided by the user, nlme will use numeric derivatives #
#                                                          #
############################################################	
	
meanfunc <- function(x,b1,b2,dose){

	tinf <- 240
	cl <- exp(b1)
	v <- exp(b2)
        t1 <- x<=tinf
        t2 <- tinf*(1-t1)+t1*x
        f1 <- (dose/cl)*(1-exp(-cl*t2/v))*exp(-cl*(1-t1)*(x-tinf)/v)

#  compute analytical dervivatives -- create the gradient matrix X

        t3 <- (1-t1)*(x-tinf)
	temp1 <- (dose/cl)*exp(-cl*t3/v)
        temp2 <- (dose/(v^2))*exp(-cl*t3/v)
        
   meangrad <- array(0,c(length(x),2),list(NULL,c("b1","b2")))
   meangrad[,"b1"] <- temp1*(-(1-exp(-cl*t2/v))*(1/cl + t3/v)+(t2/v)*exp(-cl*t2/v))*cl 
   meangrad[,"b2"] <- temp2*((1-exp(-cl*t2/v))*t3-exp(-cl*t2/v)*t2)*v
 
   attr(f1,"gradient") <- meangrad
   f1

}

############################################################
#                                                          #
#  Fit the model with ML estimation variance components    #
#                                                          #
#  -  verbose=T prints out history of each iteration       #
#                                                          #
#  -  Because the model is nonlinear, starting values      #
#     for the fixed effects (the vector beta) are provided #
#                                                          #
#  -  power variance function with power parameter         #
#     estimated from the data; starting value = 0.5        # 
#                                                          #
############################################################

outfile <- "argnlme.Rout"

cat("ARGATROBAN PK DATA",file=outfile,"\n","\n",append=F)
cat("ML fit with estimation of theta",file=outfile,"\n","\n",append=T)
	
sink(outfile,append=T)

arg.mlfit <- nlme(ys ~ meanfunc(xs,b1,b2,ds),
    fixed = list(b1 ~ 1,b2 ~1),
    random = list(b1 ~ 1,b2 ~ 1),
    groups = ~ indiv,
    data = sdat,
    start = list(fixed = c(-6.0,-2.0)),
    method="ML",verbose=T,weights=varPower(0.5))

cat("\n","\n","\n",file=outfile,append=T)
print(summary(arg.mlfit))

cat("Estimate of sigma",arg.mlfit$sigma,file=outfile,"\n",append=T)
cat("\n","\n",file=outfile,append=T)
sink()