#########################################################################
#                                                                       #
#  Program to implement normal theory ML using profiling and            #
#  the "trick" to estimate beta and theta (and sigma)                   #
#                                                                       #
#  Details on the nls() function may be found in Chapter 10             #
#  of the book "Statistical Models in S" edited by J.M. Chambers        #
#  and T. Hastie (1993, Chapman and Hall)                               #
#                                                                       #
#  Variance proportional to a power 2*theta of the mean                 #
#                                                                       #
#########################################################################

#########################################################################
#                                                                       #
#  Define the mean function f                                           #
#                                                                       #
#########################################################################

# 4 parameter logistic model

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

#  note that we must be careful at x=0 -- the function approaches
#  b2 as x goes to 0

   xx <-(x>0)*x+(x==0)*0.0001
   denom <- 1+exp(b4*(log(xx)-b3))
   f1 <- b1 + (b2-b1)/denom
   fpl <- (x>0)*f1+(x==0)*b2

#  I have removed the derivatives, because R will get confused.  The 
#  mean function is part of the "trick function."  When the "fake 
#  regression problem is run, the nls function needs the derivative of
#  the trick function.  "meanfunc" will carry around the derivatives
#  below as an attribute, and this will cause nls to screw up computation
#  of derivatives for "trkfunc," as it will attempt to make use of these
#  derivatives in the wrong way.  
        
#   meangrad <- array(0,c(length(x),4),list(NULL,c("b1","b2","b3","b4")))
#   meangrad[,"b1"] <- (x>0)*(1-1/denom)+(x==0)*0
#   meangrad[,"b2"] <- (x>0)*(1/denom)+(x==0)*1
#   meangrad[,"b3"] <- (x>0)*(b4*(b2-b1)*exp(b4*(log(xx)-b3)))/(denom**2)
#        +(x==0)*0
#   meangrad[,"b4"] <- -(x>0)*(b2-b1)*
#          exp(b4*(log(xx)-b3))*(log(xx)-b3)/(denom**2)+(x==0)*0 
#   attr(fpl,"gradient") <- meangrad
   fpl
}

#########################################################################
#                                                                       #
#  Define the variance function g                                       #
#                                                                       #
#########################################################################
  
varfunc <- function(x,b1,b2,b3,b4,theta){
     g1 <- meanfunc(x,b1,b2,b3,b4)^theta
     g1
}

#########################################################################
#                                                                       #
#  Define the mean function for the "trick"                             #
#                                                                       #
#########################################################################
  
trkfunc <- function(y,x,b1,b2,b3,b4,theta){

#  define the geometric mean of the variance function

   gdot <- exp( mean( log(varfunc(x,b1,b2,b3,b4,theta)) ) )
   trk <- (y-meanfunc(x,b1,b2,b3,b4))*gdot/varfunc(x,b1,b2,b3,b4,theta)	

#  too lazy to find analytic derivatives -- let nls compute them numerically!

   trk
}

#########################################################################
#                                                                       #
#  The data                                                             #
#                                                                       #
#########################################################################

thedat <- scan("assay.dat")
thedat <- matrix(thedat,ncol=3,byrow=T)
y <- thedat[,3]
x <- thedat[,1]

n <- length(x)

p <- 4

#########################################################################
#                                                                       #
#  Define list of starting values for beta and theta                    #
#                                                                       #
#########################################################################

allstart <- list(b1= 600, b2=6000, b3=2, b4=1, theta=1)

outfile <- "assay_ml.Rout"

#########################################################################
#                                                                       #
#  Create the data frame for nls() with "dummy" responses               #
#                                                                       #
#########################################################################

dummy <- rep(0,n)
thedat <- data.frame(dummy,x,y)

#########################################################################
#                                                                       #
#  Call nls() with the trick function and dummy data
#                                                                       #
#########################################################################
	
mlefit <- nls(dummy ~ trkfunc(y,x,b1,b2,b3,b4,theta),thedat,allstart,trace=T)
	
#########################################################################
#                                                                       #
#  Extract the estimates from the object mlefit and compute             #
#  sigma                                                                #
#                                                                       #
#########################################################################

bt <- coef(mlefit)
bmle <- bt[1:4]
theta <- bt[5]
mu <- meanfunc(x,bmle[1],bmle[2],bmle[3],bmle[4])
resid <- y-mu
sigmle <- sqrt(sum((resid/varfunc(x,bmle[1],bmle[2],bmle[3],bmle[4],theta))^2)/n)
	
#########################################################################
#                                                                       #
#  Print out the results to a file -- we round the results to 6         #
#  decimal places to the right of the decimal                           #
#                                                                       #
#########################################################################

cat("FIT OF THE ASSAY DATA BY NORMAL ML USING THE TRICK",
	file=outfile,"\n","\n","\n",append=F)
	
cat("beta ",round(bmle,6),file=outfile,"\n","\n","\n",append=T)
cat("theta ",round(theta,6),file=outfile,"\n","\n","\n",append=T)
cat("sigma ",round(sigmle,6),file=outfile,"\n","\n","\n",append=T)