# Fit Michaelis-Menten data from Bates & Watts

conc.treated = c(rep(.02,2),rep(.06,2),rep(.11,2),
                        rep(.22,2),rep(.56,2),rep(1.1,2))
conc.untreated = conc.treated[-12]                       
vel.treated = c(76,47,97,107,123,139,159,152,191,201,207,200)
vel.untreated = c(67,51,84,86,98,115,131,124,144,158,160)
Micmen.treated = data.frame(conc=conc.treated, vel=vel.treated)
Micmen.untreated = data.frame(conc=conc.untreated, vel=vel.untreated) 

# Get starting values from transformed model with linear response

Micmen = Micmen.treated
attach(Micmen)

lin.params = lsfit(x=1/Micmen$conc, y=1/Micmen$vel)$coef
names(lin.params) = NULL #Without this, the names 'X' and 'intercept' would
                          #be associated with 'theta1' and 'theta2'
theta1.start = 1/lin.params[1]; theta2.start = lin.params[2]/lin.params[1]
starting.values = list(theta1.start,theta2.start)
names(starting.values) = c("theta1","theta2")


cat("Starting values are","\n"); print(starting.values) 
cat("\n","Trace of sum of squares function and parameter values:","\n","\n")
fit1 = nls(vel~theta1*conc/(theta2+conc),Micmen,
                  start = starting.values, trace=T)
cat("\n","'fit1' is the call to nls( ... ) ","\n","\n") 
print(fit1)
cat("\n","summary(fit1) gives the following output:","\n")
print(summary(fit1))

# Define the response function AND its gradient:
velocity = function(conc,theta1,theta2)  {
        velocity = theta1*conc/(theta2+conc)
       # Compute the derivatives w.r.t. theta1 and theta2: 
        dvel.theta1 = conc/(theta2+conc)
        dvel.theta2 = -theta1*conc/((theta2+conc)^2)
       # Bind them together into an nxp matrix, as an 'attribute' of vel:         
        attr(velocity, "gradient") = cbind(dvel.theta1, dvel.theta2)        
        velocity
   }          

cat("\n","\n","\n",
       "'fit2' is the call to nls( ... ) using the gradient","\n","\n") 

fit2 = nls(vel~velocity(conc,theta1,theta2),Micmen,
            start=starting.values, trace=T)
cat("\n","\n")
print(fit2)
cat("\n","\n"); print(summary(fit2))


#Get the derivative matrix, as a function and 
# evaluated at the final estimates:

V = function(theta1,theta2) {
        attr(velocity(conc,theta1,theta2),"gradient")
  }

coefs = coef(fit2)
resids = residuals(fit2)

V.hat = V(coefs[1], coefs[2])

cat("\n","\n",
"The derivative matrix evaluated at the final parameter estimates is","\n","\n")
print(V.hat)

cat("\n","A linear regression of residuals(fit2) on the columns of the 
derivative matrix V.hat gives the following output:","\n","\n")
fin.lin.fit = lsfit(x=V.hat,y=resids,intercept=F)
ls.print(fin.lin.fit)

#Prepare plot of estimated vel vs. conc on a grid of values, with data points
#on the same axes; include simultaneous confidence bands.

conc.grid = seq(from=0,to=1.2,by=.02)
vel.plot = velocity(conc.grid,coefs[1],coefs[2])
 
vecnorm = function(vec) sqrt(crossprod(vec))

  # Get half width of 100*(1 - alpha)% confidence band:
half.width = function(x, alpha)  {
    quant = qf(1-alpha, 2, 10)
    sigmahat = vecnorm(resids)/sqrt(10) 
    R1 = qr.R(fin.lin.fit$qr)
    v0 = attr(velocity(x,coefs[1],coefs[2]),"gradient")
    norm = vecnorm(solve(t(R1),t(v0)))
    half.width = sigmahat*norm*sqrt(2*quant)
    half.width
  }

alpha = .05
half = NULL; for (x in conc.grid) half = c(half, half.width(x, alpha))
upper.limits = vel.plot + half
lower.limits = vel.plot - half

par(oma=c(8,0,0,0))
plot(x=conc.grid, y=vel.plot, type = "l", ylim = c(0,250), 
        xlab = "concentration", ylab = "velocity")
points(conc,vel,pch="*")
lines(conc.grid, upper.limits, lty=2)
lines(conc.grid, lower.limits, lty=2)
mtext("Puromycin data, estimated Michaelis-Menten response, 95% inference band",
  side=1, outer=T)