######################################################################### # # # Program to implement the 3-step GLS algorithm with # # cmax total iterations, theta known, using the R # # function nls(), estimating theta by PL # # # # Details on the nls() function may be found in Chapter 10 # # of the book "Statistical Models in S" edited by J.M. Chambers # # and T. Hastie (1993, Chapman and Hall) # # Variance proportional to a power 2*theta of the mean # # # # The program may be used for any problem by changing the # # code defining the mean function and "weights" # # # ######################################################################### ######################################################################### # # # Define the mean function f and the gradient matrix # # of its partial derivatives with respect to beta (n x p) as an # # attribute. The nls() function will know to use analytic derivatives # # when it spots the presence of the attribute "gradient" defined # # along with the function # # # ######################################################################### meanfunc <- function(x,b1,b2,b30,b31,d){ b3 <- b30*(1-d)+b31*d piece <- 1+(b2/x)^b3 f1 <- b1/piece # compute analytical dervivatives -- create the gradient matrix X meangrad <- array(0,c(length(x),4),list(NULL,c("b1","b2", "b30","b31"))) meangrad[,"b1"] <- 1/piece meangrad[,"b2"] <- (-b1/piece^2)*((b2/x)^b3)*(b3/b2) meangrad[,"b30"] <- (-b1/piece^2)*((b2/x)^b3)*log(b2/x)*(1-d) meangrad[,"b31"] <- (-b1/piece^2)*((b2/x)^b3)*log(b2/x)*d attr(f1,"gradient") <- meangrad f1 } ######################################################################### # # # To implement step (iii), we wish to do weighted least squares # # with known weights. This is accomplished by transforming the # # the response and mean function to a problem with constant variance, # # as shown in the notes. nls() is then called to # # do OLS on the transformed problem, thereby doing WLS # # # # Because the weights in this case are a function of the mean, # # I also define the mean function with no gradient attribute. # # This is because the transformed mean function will depend # # on the current estimated weights, which are found by evaluating # # the mean function at the current estimate. Because of the nature # # of R attributes, if calculated using the function # # "indofunc" above, the weights will carry along the gradient # # attribute, which will confuse things when we wish to calculate # # the gradient of the transformed mean function assuming the # # weights are constant! There are more elegant ways around this, # # but doing it this way here is meant to highlight the issue so # # you will be aware of it. # # # ######################################################################### unweightfunc <- function(x,b1,b2,b30,b31,d){ b3 <- b30*(1-d)+b31*d piece <- 1+(b2/x)^b3 f1 <- b1/piece f1 } # The transformed mean function (multiplied by the square root # of the current estimated weights, which are considered fixed) weightfunc <- function(x,b1,b2,b30,b31,d,wt){ b3 <- b30*(1-d)+b31*d piece <- 1+(b2/x)^b3 f1 <- b1/piece w12 <- sqrt(wt) weightf <- f1*w12 # compute analytical dervivatives -- create the gradient matrix X weightgrad <- array(0,c(length(x),4),list(NULL,c("b1","b2","b30","b31"))) weightgrad[,"b1"] <- w12*(1/piece) weightgrad[,"b2"] <- w12*(-b1/piece^2)*((b2/x)^b3)*(b3/b2) weightgrad[,"b30"] <- w12*(-b1/piece^2)*((b2/x)^b3)*log(b2/x)*(1-d) weightgrad[,"b31"] <- w12*(-b1/piece^2)*((b2/x)^b3)*log(b2/x)*d attr(weightf,"gradient") <- weightgrad weightf } ######################################################################### # # # The data # # # ######################################################################### thedata <- scan("snake.dat") thedata <- matrix(thedata,ncol=3,byrow=T) theta <- 0.8 outfile <- "snake_known.Rout" x <- thedata[,1] y <- thedata[,2] delta <- thedata[,3] n <- length(x) p <- 4 ######################################################################### # # # Define list of starting values for beta # # # ######################################################################### bstart <- list(b1=650,b2=300,b30=1,b31=1) ######################################################################### # # # Create the data frame for nls() # # # ######################################################################### indat <- data.frame(x,y) ######################################################################### # # # Specify the max number of iterations of GLS (C) # # # ######################################################################### cmax <- 10 ######################################################################### # # # Step (i) -- initial fit by OLS # # # ######################################################################### ols.fit <- nls(y ~ meanfunc(x,b1,b2,b30,b31,delta),indat,bstart) ######################################################################### # # # Extract the estimate # # # ######################################################################### bols <- coef(ols.fit) mu <- unweightfunc(x,bols[1],bols[2],bols[3],bols[4],delta) resid <- y-mu sigols <- sqrt(sum(resid*resid)/(n-p)) ######################################################################### # # # Print out the results to a file -- we round the results to 6 # # decimal places to the right of the decimal # # # ######################################################################### cat("FIT OF THE SNAKE DATA BY GLS",file=outfile,"\n", append=F) cat("WITH POWER VARIANCE, THETA KNOWN AND EQUAL TO 0.8", file=outfile,"\n","\n","\n",append=T) cat("OLS estimate ",round(bols,6),file=outfile,"\n","\n",append=T) cat("Estimate of assumed constant SD of y ",round(sigols,6),file=outfile, "\n","\n","\n",append=T) ######################################################################### # # # Use the OLS estimator as the preliminary estimator # # # ######################################################################### bgls <- bols ######################################################################### # # # Begin iterating between steps (ii) and (iii) # # # ######################################################################### for (k in 1:cmax){ mu <- unweightfunc(x,bgls[1],bgls[2],bgls[3],bgls[4],delta) wt <- 1/(mu^(2*theta)) ywt <- y*sqrt(wt) ######################################################################### # # # Step (iii) -- update estimation of beta by WLS with the weights # # held fixed. First, create the updated data frame of transformed # # responses and weights for use by nls() # # # ######################################################################### indat2 <- data.frame(x,ywt,wt) gls.fit <- nls(ywt ~ weightfunc(x,b1,b2,b30,b31,delta,wt),indat2,bstart) ######################################################################### # # # Get the updated GLS estimate to use for constructing weights # # on the next iteration # # # ######################################################################### bgls <- coef(gls.fit) ######################################################################### # # # Print results of this iteration to the output file # # # ######################################################################### cat("GLS estimate of beta ",round(bgls,6),"\n","\n",file=outfile, append=T) } ######################################################################### # # # Finished iteration loop -- now compute the estimate of sigma^2 # # based on the final GLS estimate. Use the "adjusted" version # # # ######################################################################### mu <- unweightfunc(x,bgls[1],bgls[2],bgls[3],bgls[4],delta) resid <- y-mu g <- mu^theta sigma2 <- sum((resid/g)**2)/(n-p) sigma <- sqrt(sigma2) ######################################################################### # # # Print out the final estimate of sigma and the summary provided by # # the nls() function # # # ######################################################################### cat("Final estimate of sigma ",round(sigma,6),"\n","\n", file=outfile,append=T) sink(outfile,append=T) print(summary(gls.fit)) sink() ######################################################################### # # # Additional code to make the plots # # # ######################################################################### # separate side-by-side plots for each type of snake postscript("snake.ps",width=6) par(mfrow=c(2,1),cex=0.75,oma=c(0,0,2,0),mar=c(4,4,3,2)) # grid for plotting fits ymax <- max(y) xgrid <- seq(31,2000,1) y0 <- y[delta==0] x0 <- x[delta==0] y1 <- y[delta==1] x1 <- x[delta==1] # Northern (delta=0) f0ols <- unweightfunc(xgrid,bols[1],bols[2],bols[3],bols[4],0) f0gls <- unweightfunc(xgrid,bgls[1],bgls[2],bgls[3],bgls[4],0) plot(x0,y0,,lty=1,xlab="substrate concentration",ylab="reaction rate",ylim=c(0,ymax),pch="0",cex=0.8) title("Northern") lines(xgrid,f0ols,lty=1) lines(xgrid,f0gls,lty=2) # Southern (delta=0) f1ols <- unweightfunc(xgrid,bols[1],bols[2],bols[3],bols[4],1) f1gls <- unweightfunc(xgrid,bgls[1],bgls[2],bgls[3],bgls[4],1) plot(x1,y1,,lty=1,xlab="substrate concentration",ylab="reaction rate",ylim=c(0,ymax),pch="0",cex=0.8) title("Southern") lines(xgrid,f1ols,lty=1) lines(xgrid,f1gls,lty=2) graphics.off()