######################################################################### # # # Program to implement the 3-step GLS algorithm with # # cmax total iterations, theta known, using the R # # function nls() # # # # Variance proportional to a power 2*theta of the mean # # theta known chosen by the user. # # # # 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,d,b10,b11,b20,b21,b3){ t0 <- 20 b1 <- b10*(1-d) + b11*d b2 <- b20*(1-d) + b21*d piece <- (1-exp(-b2*x))/(1-exp(-b2*t0)) pieceb3 <- piece^b3 f1 <- b1*pieceb3 # compute analytical dervivatives -- create the gradient matrix X fb2 <- f1*b3*(x*exp(-b2*x)/(1-exp(-b2*x)) - t0*exp(-b2*t0)/(1-exp(-b2*t0))) meangrad <- array(0,c(length(x),5), list(NULL,c("b10", "b11","b20","b21","b3"))) meangrad[,"b10"] <- pieceb3*(1-d) meangrad[,"b11"] <- pieceb3*d meangrad[,"b20"] <- (1-d)*fb2 meangrad[,"b21"] <- d*fb2 meangrad[,"b3"] <- f1*log(piece) 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. # # 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, # # we 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 attributes, if calculated using the function # # "meanfunc" 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,d,b10,b11,b20,b21,b3){ t0 <- 20 b1 <- b10*(1-d) + b11*d b2 <- b20*(1-d) + b21*d piece <- (1-exp(-b2*x))/(1-exp(-b2*t0)) pieceb3 <- piece^b3 f1 <- b1*pieceb3 f1 } # The transformed mean function (multiplied by the square root # of the current estimated weights, which are considered fixed) weightfunc <- function(x,d,b10,b11,b20,b21,b3,wt){ t0 <- 20 b1 <- b10*(1-d) + b11*d b2 <- b20*(1-d) + b21*d piece <- (1-exp(-b2*x))/(1-exp(-b2*t0)) pieceb3 <- piece^b3 f1 <- b1*pieceb3 w12 <- sqrt(wt) weightf <- f1*w12 # compute analytical dervivatives -- create the gradient matrix X fb2 <- f1*b3*(x*exp(-b2*x)/(1-exp(-b2*x)) - t0*exp(-b2*t0)/(1-exp(-b2*t0))) weightgrad <- array(0,c(length(x),5), list(NULL,c("b10", "b11","b20","b21","b3"))) weightgrad[,"b10"] <- w12*pieceb3*(1-d) weightgrad[,"b11"] <- w12*pieceb3*d weightgrad[,"b20"] <- w12*(1-d)*fb2 weightgrad[,"b21"] <- w12*d*fb2 weightgrad[,"b3"] <- w12*f1*log(piece) attr(weightf,"gradient") <- weightgrad weightf } ######################################################################### # # # The data # # # ######################################################################### thedata <- scan("trees.dat") thedata <- matrix(thedata,ncol=3,byrow=T) delta <- thedata[,1] x <- thedata[,2] y <- thedata[,3] n <- length(x) p <- 5 theta <- 1.0 bstart <- list(b10=20,b11=20,b20=0.1,b21=0.1,b3=2) outfile <- "trees_known.Rout" ######################################################################### # # # Create the data frame for nls() # # # ######################################################################### thedat <- data.frame(x,y,delta) ######################################################################### # # # Specify the max number of iterations of GLS (C) # # # # Alternatively, we could check for convergence after step (iii) # # # ######################################################################### cmax <- 10 ######################################################################### # # # Step (i) -- initial fit by OLS # # # # A call to nls() is pretty self-explanatory. The first argument # # specifies the model: on the LHS of the ~ is the response variable, # # on the RHS is the mean function (may also be just an expression; # # it need not be a function call). The second argument is the name # # of the dataframe where the data reside, and the third is a list # # containing the starting values for the Gauss-Newton algorithm. # # Additional options are described in the Chambers/Hastie book. # # The algorithm employs a form of step-halving; this also described # # in the book. The call creates an object containing the # # parameter estimates and other summary information from the fit. # # # ######################################################################### olsfit <- nls(y ~ meanfunc(x,delta,b10,b11,b20,b21,b3),thedat,bstart) ######################################################################### # # # Extract the estimate from the object olsfit # # # ######################################################################### bols <- coef(olsfit) ######################################################################### # # # 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 TREE DATA BY GLS",file=outfile,"\n","\n","\n", append=F) cat("OLS estimate ",round(bols,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){ ######################################################################### # # # Step (ii) -- calculate the weights and the transformed response # # to use in the WLS calculation in (iii) # # # ######################################################################### mu <- unweightfunc(x,delta,bgls[1],bgls[2],bgls[3],bgls[4],bgls[5]) 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() # # # ######################################################################### thedat2 <- data.frame(x,ywt,wt,delta) glsfit <- nls(ywt ~ weightfunc(x,delta,b10,b11,b20,b21,b3,wt),thedat2,bstart) ######################################################################### # # # Get the updated GLS estimate to use for constructing weights # # on the next iteration # # # ######################################################################### bgls <- coef(glsfit) ######################################################################### # # # Print results of this iteration to the output file # # # ######################################################################### cat("Iteration ",k,"\n",file=outfile,append=T) 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,delta,bgls[1],bgls[2],bgls[3],bgls[4],bgls[5]) 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(glsfit)) sink() ######################################################################### # # # Additional code to make the plots # # # ######################################################################### d <- "0" postscript("treesdat.ps",width=9) plot(x,y,type="p",xlab="age (years)",ylab="dominant height (m)",cex=0.8, pch=as.character(delta)) postscript("treesfit.ps",width=9) xgrid <- seq(1,80,0.05) fols0 <- unweightfunc(xgrid,0,bols[1],bols[2],bols[3],bols[4],bols[5]) fols1 <- unweightfunc(xgrid,1,bols[1],bols[2],bols[3],bols[4],bols[5]) fgls0 <- unweightfunc(xgrid,0,bgls[1],bgls[2],bgls[3],bgls[4],bgls[5]) fgls1 <- unweightfunc(xgrid,1,bgls[1],bgls[2],bgls[3],bgls[4],bgls[5]) plot(x,y,type="p",xlab="age (years)",ylab="dominant height (m)",cex=0.8, pch=as.character(delta)) lines(xgrid,fols0,lty=2) lines(xgrid,fols1,lty=2) lines(xgrid,fgls0,lty=1) lines(xgrid,fgls1,lty=1) graphics.off()