# 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)