# Fit the model in problem 3 on
# assignment 9 in Spring 2002

kinetics <- read.table("c:/courses/st511/snew/kinetics.dat",
            col.names=c("Y","X"))

# Compute starting values

YY <- 1/kinetics$Y
XX <- 1/kinetics$X
lm.out <- lm(YY ~ XX)

start.val <- coef(lm.out)                 # Note that Intercept = 1/b0
start.val[1] <- 1/start.val[1]            # Therefore, b0 = 1/Intercept
start.val[2] <- start.val[1]*start.val[2] # b1 = b0 * (b1/b0)
names(start.val) <- c("b0","b1")          # specify the name of each starting
                                          # value ( to use this in nls())
start.val

# Fit a nonlinear model

hw9.3.nls <- nls(formula = Y ~ b0*X/(b1+X),
                    data = kinetics,
                   start = start.val,
                   trace = T)
summary(hw9.3.nls)

# Plot the points and the curve

coeff <- coef(hw9.3.nls)  # estimated coefficient
fit.model  <- function(x) (coeff[1] * x) / (coeff[2] + x)
x <- seq(1,40, length=100)
y <- fit.model(x)

par(mkh=.1,mex=1.5)
plot(c(0,40),c(0,25),main="Estimated Curve",xlab="X",
              ylab="Fitted Value & Obs.",type="n")
lines(x,y,lty=1)
points(kinetics$X,kinetics$Y,pch=16)

# Make residual plots

plot(kinetics$X, hw9.3.nls$res, xlab="X", ylab="Residuals", 
      main="Residuals against X"); abline(h=0,lty=2)

# Normal probability plot for residuals

qqnorm(hw9.3.nls$res,main="Normal Probability Plot of Residuals")
qqline(hw9.3.nls$res)



# Make 95% confidence intervals for b0 and b1.

# Profile method  You must add the MASS library

library(MASS)

rbind(b0 = confint.nls(hw9.3.nls,parm="b0"),
      b1 = confint.nls(hw9.3.nls,parm="b1"))

# Large sample normal approximation
est.b
est.b(hw9.3.nls)  

# Estimate the concentration (x) at which the expeceted 
# velocity of the reaction is 25.  Report a 95% C.I.

# First evaluate partial derviatives of the mean function  
# needed to apply the delta method

x.pds <- deriv(~ b1/(b0/y -1),
                c("b0","b1"), function(y, b0, b1) NULL)

#  list the function

x.pds

#  Compute the estimates of the times and the large sample
#  covariance matrix and standard errors

est.x <- function(y, obj,level=.95) 
{      b <- coef(obj)
   expr2 <- (b[1]/y) - 1
   x.hat <-  b[2]/expr2
       G <- matrix(0,nrow=1,ncol=length(b))
     G[1,1] <-  - ((b[2] * (1/y))/(expr2^2))
     G[1,2] <-  - 1/(expr2)
       v <- (G %*% vcov(obj) %*% t(G))
       n <- length(obj$fitted.values) - length(obj$parameters)
     low <- x.hat - qt(1 - (1 - level)/2, n)*sqrt(v)
    high <- x.hat + qt(1 - (1 - level)/2, n)*sqrt(v)
  result <- c(x.hat = x.hat, std.err = sqrt(v), lower.bound=low, upper.bound=high)
  return(G.hat=G, V.hat= vcov(obj), S2 = v, result=result)
}
est.x(y=25, hw9.3.nls)