#  fit the nonlinear mixed model to the HCV data using nlme( )

library(nlme)

meanfunc <- function(x,b1,b2,b3){

    v0 <- exp(b1)
    cc <- exp(b2)
    e <- exp(b3)
    t0 <- 0.20      #  t0 is fixed
    tt <- (x>t0)*(x-t0)
    fpl <- v0*(1-e+e*exp(-cc*tt))

#  compute analytical dervivatives

   meangrad <- array(0,c(length(x),3),list(NULL,c("b1","b2","b3")))
   meangrad[,"b1"] <- fpl
   meangrad[,"b2"] <- -v0*e*exp(-cc*tt)*cc*tt
   meangrad[,"b3"] <- v0*e*(exp(-cc*tt)-1)
   attr(fpl,"gradient") <- meangrad
   fpl
}

thedat <- matrix(scan("hcvmix.dat"),ncol=3,byrow=T)

id <- thedat[,1]
days <- thedat[,2]
vl <- thedat[,3]

vdat <- data.frame(id,days,vl)

outfile <- "hcvnlme.Rout"

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

sink(outfile,append=T)

hcv.mlfit <- nlme(vl ~ meanfunc(days,b1,b2,b3),
   fixed=list(b1 ~ 1, b2 ~ 1, b3 ~ 1),
   random=list(b1 ~ 1, b2 ~ 1, b3 ~ 1),
   groups = ~id,
   data=vdat,
   start=list(fixed=c(1.5,1.5,-0.15)),
   method="ML",verbose=T,weights=varPower(1.0))

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

sds <- c(0.19534581,0.11546963,0.01854371 )
sdmat <- diag(sds)
corrmat <- matrix(c(1,-0.418,-0.117,
                    -0.418,1,0.352,
                    -0.117,0.352,1),3,3)

cat("\n","\n",file=outfile,append=T)
cat("D matrix","\n",file=outfile,append=T)


print(sdmat%*%corrmat%*%sdmat)

sink()