# S-Plus code for problem 2 on Assignment 9 
# in 2002.  The data are posted as pnonlin.dat
# This file is posted as pnonlin.ssc

# Create a data frame

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

# Get starting values for fitting the model
# Note that Intercept in the following model is
# approximately log(B0)  So a starting value for
# B0 is B0 = exp(Intercept)

lm.out <- lm(log(Y) ~ log(X1)+log(X2), data=pnonlin)
start.val <- coef(lm.out)             
start.val[1] <- exp(start.val[1])  

# Specify the name of the starting values
# This must match the model formula used
# by the nls() function
   
names(start.val) <- c("b0","b1","b2") 
start.val


# Obtain least squares estimates of model parameters
# This procedure uses numerical approximations for
# first partial derivatives

hw9.2.nls <- nls(formula = Y ~ b0*(X1^b1)*(X2^b2),
                    data = pnonlin,
                   start = start.val,
                   trace = T)
summary(hw9.2.nls)

#  Create a function to compute points
# of the estimated surface

coeff <- coef(hw9.2.nls)
fit.model  <- function(x1, x2, coeff) {
    coeff[1] * (x1^coeff[2]) * (x2^coeff[3])}



# Both x1 and x2 have only 3 levels, so 
# you can slice the 3D plot and present 6 xy-plots.

par(mkh=.1,mex=1.5)

X1 <- seq(0,110,length=100)
plot(c(0,110),c(0,500),main="Slices of the 3D plot",xlab="X1",
            ylab="Fitted Value",type="n")

X2 <- 1                   # Fix X2 = 1
Z  <- fit.model(X1,X2,coeff)
lines(X1,Z,lty=1)
points(pnonlin$X1[pnonlin$X2==1 ],pnonlin$Y[pnonlin$X2==1],pch=16)

X2 <- 10                  # Fix X2 = 10
Z  <- fit.model(X1,X2,coeff)
lines(X1,Z,lty=3)
points(pnonlin$X1[pnonlin$X2==10 ],pnonlin$Y[pnonlin$X2==10],pch=17)

X2 <- 100                 # Fix X2 = 100
Z  <- fit.model(X1,X2,coeff)
lines(X1,Z,lty=5)
points(pnonlin$X1[pnonlin$X2==100 ],pnonlin$Y[pnonlin$X2==100],pch=18)

legend(x=5.3374,y=446.2853,c("X2 = 1","X2 = 10","X2 = 100"), 
       lty=c(1,3,5),marks = c(16,17,18), pch = "  ",bty = "n")


par(mkh=.1,mex=1.5)

X1 <- seq(0,110,length=100)
plot(c(0,110),c(0,500),main="Slices of the 3D plot",xlab="X1",
            ylab="Fitted Value",type="n")

X1 <- 1                   # Fix X1 = 1
Z  <- fit.model(X1,X2,coeff)
lines(X2,Z,lty=1)
points(pnonlin$X2[pnonlin$X1==1 ],pnonlin$Y[pnonlin$X1==1],pch=16)



X1 <- 10                  # Fix X1 = 10
Z  <- fit.model(X1,X2,coeff)
lines(X2,Z,lty=3)
points(pnonlin$X2[pnonlin$X1==10 ],pnonlin$Y[pnonlin$X1==10],pch=17)

X1 <- 100                 # Fix X1 = 100
Z  <- fit.model(X1,X2,coeff)
lines(X2,Z,lty=5)
points(pnonlin$X2[pnonlin$X1==100 ],pnonlin$Y[pnonlin$X1==100],pch=18)

legend(x=5.3374,y=446.2853,c("X1 = 1","X1 = 10","X1 = 100"), 
       lty=c(1,3,5),marks = c(16,17,18), pch = "  ",bty = "n")


# Construct residual plots for part (d)

plot(pnonlin$X1, hw9.2.nls$res, xlab="X1", ylab="Residuals", 
     main="Residuals against X1"); abline(h=0,lty=2)
plot(pnonlin$X2, hw9.2.nls$res, xlab="X2", ylab="Residuals", 
     main="Residuals against X2"); abline(h=0,lty=2)

# Construct a normal probability plot for the residuals

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


# Use the deriv3() function available in the
# MASS library to compute first and second partial
# derivatives of the mean function, evaluated at 
# final solution for the parameter estimates.

library(MASS)

yield.2pds <- deriv3(~ b0*x1^b1*x2^b2,c("b0","b1","b2"), 
                       function(x1,x2, b0, b1, b2) NULL)

hw9.2.nls2 <- nls(formula = Y ~ yield.2pds(X1,X2,b0,b1,b2),
                    data = pnonlin,
                   start = start.val,
                   trace = T)


# Construct 95% confidence intervals for b0, b1, and b2.
# Use the built in profile method

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

# Use the standard large sample approach
est.b
est.b(hw9.2.nls) 


# Construct a 95% confidence interval for the mean 
# yield when x1=50 and x2=60.

# Evaluate first partial derviatives of the mean function 
# needed to apply the delta method
#  Compute the estimates of the yields and the large sample
#  covariance matrix and standard errors

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


est.yield(x1=50,x2=60, hw9.2.nls)


# t-test for difference in parameters

      C <- matrix(c(0,1,-1),nrow=1)
      b <- coef(hw9.2.nls) 
   V.Cb <-  C %*% vcov(hw9.2.nls) %*% t(C)
 t.stat <- (C %*% b) / sqrt(V.Cb)
     df <- length(hw9.2.nls$fitted) - length(hw9.2.nls$par)
  p.val <- 2*(1- pt(abs(t.stat),df))
 result <- c("b1-b2"= C %*% b, std.err=sqrt(V.Cb), t.stat=t.stat,df=df,p.val=p.val)
 result