# 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