days = c(1,2,3,4,5,7)
oxy = c(8.3,10.3,19,16,15.6,19.8)
bod = data.frame(days,oxy)

# Compute starting values.  The value '20' was chosen as an approximate
# asymptote, upon examining the data
theta1 = 20
theta2 = -1*lsfit(x=days,y=log(1-oxy/theta1),intercept=F)$coef
starting.values = c(theta1,theta2)

# Define the response function and its gradient

V1 = function(theta2) 1-exp(-theta2*days) 
V2 = function(theta1,theta2) theta1*days*exp(-theta2*days)
V = function(theta1,theta2)  cbind(V1(theta2),V2(theta1,theta2))

oxygen = function(days,theta1,theta2)  {
 value = theta1*(1-exp(-theta2*days))
 attr(value,"gradient") = V(theta1,theta2)
 value
     }

# Do the nonlinear fit

names(starting.values) = c("theta1","theta2")

fit.bod = nls(oxy~oxygen(days,theta1,theta2),bod,start = starting.values,trace=T)
print(summary(fit.bod))

coefs = coef(fit.bod)
resids = residuals(fit.bod)

theta1=coefs[1];  theta2=coefs[2]

opar = par(mfrow = c(1,2), oma = c(3, 0, 0, 0), las = 1)

plot(days, oxy)
newdata = seq(from = 1, to = 7, length = 50)
lines(newdata, predict(fit.bod, list(days=newdata)))

# Plot S(theta)

grid1 = seq(0,60,length=50); grid2 = seq(0,6,length=50)
S.theta = function(theta.one,theta.two)  {
     sum((oxy-theta.one*(1-exp(-theta.two*days)))^2)
    }    # S.theta is the sum of squares function
S.mat = matrix(0,length(grid1),length(grid2))
for (i in 1:length(grid1))  {
     for (j in 1:length(grid2))
       S.mat[i,j] = S.theta(grid1[i],grid2[j])
     }
p = length(coefs); n = length(oxy)
S.max = function(conf.level)  {
     S.theta(theta1,theta2)*(1+(p/(n-p))*qf(conf.level,p,n-p))
    }

contour(grid1,grid2,S.mat,levels=c(S.max(.95)),xlab="",ylab="theta2")
points(theta1,theta2,pch="+")
title(sub="\n theta1 \n 80% and 95% l'hood regions;\n LSE shown as '+'", line=4)
contour(grid1,grid2,S.mat,levels=c(S.max(.80)),add=T, lty=4)
par(opar)

win.graph()
pr.bod = profile(fit.bod)
opar = par(mfrow = c(2,2), oma = c(1.1, 0, 1.1, 0), las = 1)
plot(pr.bod, conf = c(95, 90, 80, 50)/100)
plot(pr.bod, conf = c(95, 90, 80, 50)/100, absVal = FALSE)
mtext("Confidence intervals based on the profile sum of squares",
      side = 3, outer = TRUE)
mtext("BOD data - confidence levels of 50%, 80%, 90% and 95%",
      side = 1, outer = TRUE)
par(opar)