######################################################################### # # # Program to implement the 3-step GLS algorithm with # # cmax total iterations, theta unknown, using the Splus # # function nls() # # # # 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) # # # # Applied here to the indomethacin data assuming variance proprotional # # to a power 2*theta of the mean response ("Power-of-mean") # # model. Theta is estimated from the data using the pseudo- # # likelihood method, implemented by "tricking" the nonlinear # # regression routine. The mean model f is the biexponential # # model parameterized to enforce positivity # # # # The program may be used for any problem by changing the # # code defining the mean function and "weights" # # # ######################################################################### ######################################################################### # # # Define the bioexponential 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 # # # ######################################################################### indofunc <- function(time,b1,b2,b3,b4){ eb1 <- exp(b1) eb2 <- exp(b2) eb3 <- exp(b3) eb4 <- exp(b4) indof <- eb1*exp(-eb2*time)+eb3*exp(-eb4*time) # compute analytical dervivatives -- create the gradient matrix X indograd <- array(0,c(length(time),4),list(,NULL,c("b1","b2","b3","b4"))) indograd[,"b1"] <- eb1*exp(-eb2*time) indograd[,"b2"] <- -eb1*eb2*time*exp(-eb2*time) indograd[,"b3"] <- eb3*exp(-eb4*time) indograd[,"b4"] <- -eb3*eb4*time*exp(-eb4*time) attr(indof,"gradient") <- indograd indof } ######################################################################### # # # 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, # # 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 Splus 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(time,b1,b2,b3,b4){ eb1 <- exp(b1) eb2 <- exp(b2) eb3 <- exp(b3) eb4 <- exp(b4) uwtf <- eb1*exp(-eb2*time)+eb3*exp(-eb4*time) uwtf } # The transformed mean function (multiplied by the square root # of the current estimated weights, which are considered fixed) weightfunc <- function(time,b1,b2,b3,b4,wt){ eb1 <- exp(b1) eb2 <- exp(b2) eb3 <- exp(b3) eb4 <- exp(b4) pred <- unweightfunc(time,b1,b2,b3,b4) w12 <- sqrt(wt) weightf <- pred*w12 # compute analytical dervivatives -- create the gradient matrix X weightgrad <- array(0,c(length(time),4),list(,NULL,c("b1","b2","b3","b4"))) weightgrad[,"b1"] <- eb1*exp(-eb2*time)*w12 weightgrad[,"b2"] <- -eb1*eb2*time*exp(-eb2*time)*w12 weightgrad[,"b3"] <- eb3*exp(-eb4*time)*w12 weightgrad[,"b4"] <- -eb3*eb4*time*exp(-eb4*time)*w12 attr(weightf,"gradient") <- weightgrad weightf } ######################################################################### # # # To estimate theta, we use the "trick" of turning the PL estimation # # problem into a "nonlinear regression" problem. Here, the function # # trkfunc() is the mean function for the "dummy" regression problem # # that is solved to estimated theta. # # # ######################################################################### trkfunc <- function(resid,mudot,mu,theta){ trk <- resid*((mudot/mu)**theta) trkgrad <- array(0,c(length(mu),1),list(,NULL,c("theta"))) trkgrad[,"theta"] <- trk*log(mudot/mu) attr(trk,"gradient") <- trkgrad # analytic derivative trk } ######################################################################### # # # The data -- alternatively, we could read them from a file, of course # # Note time has already been put into hours # # # ######################################################################### # subject 5 time <- c(0.25,0.50,0.75,1.00,1.25,2.00,3.00,4.00,5.00,6.00,8.00) conc <- c(2.05,1.04,0.81,0.39,0.30,0.23,0.13,0.11,0.08,0.10,0.06) n <- length(conc) p <- 4 ######################################################################### # # # Create the data frame for nls() # # # ######################################################################### indodat <- data.frame(time,conc) ######################################################################### # # # Specify the max number of iterations of GLS (C) # # # ######################################################################### 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 Splus object containing the # # parameter estimates and other summary information from the fit. # # # ######################################################################### indo.olsfit <- nls(conc ~ indofunc(time,b1,b2,b3,b4),indodat, list(b1=0.69, b2=0.69, b3=-1.6, b4=-1.6)) ######################################################################### # # # Extract the estimate from the Splus object indo.olsfit # # # ######################################################################### bols <- indo.olsfit$parameters ######################################################################### # # # 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 INDOMETHACIN DATA BY GLS",file="indogls2.sout","\n","\n","\n", append=F) cat("OLS estimate ",round(bols,6),file="indogls2.sout","\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) -- estimate theta and calculate the weights and the # # transformed response to use in the WLS calculation in (iii) # # # ######################################################################### # compute the geometric mean and predicted values mu <- unweightfunc(time,bgls[1],bgls[2],bgls[3],bgls[4]) mudot <- prod(mu)**(1/n) mudot <- rep(mudot,n) resid <- conc-mu dummy <- rep(0,length(time)) pldat <- data.frame(resid,mu,mudot) # Now estimate theta using the PL trick indo.plfit <- nls(dummy ~ trkfunc(resid,mudot,mu,theta),pldat, list(theta=1)) theta <- indo.plfit$parameters cat("Iteration ",k,"\n",file="indogls2.sout",append=T) cat("Estimate of theta ",round(theta,4),"\n", file="indogls2.sout",append=T) mu <- unweightfunc(time,bgls[1],bgls[2],bgls[3],bgls[4]) wt <- 1/(mu^(2*theta)) concwt <- conc*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() # # # ######################################################################### indodat2 <- data.frame(time,concwt,wt) indo.glsfit <- nls(concwt ~ weightfunc(time,b1,b2,b3,b4,wt),indodat2, list(b1=0.69, b2=0.69, b3=-1.6, b4=-1.6)) ######################################################################### # # # Get the updated GLS estimate to use for constructing weights # # on the next iteration # # # ######################################################################### bgls <- indo.glsfit$parameters ######################################################################### # # # Print results of this iteration to the output file # # # ######################################################################### cat("GLS estimate of beta ",round(bgls,6),"\n","\n",file="indogls2.sout", append=T) } ######################################################################### # # # Finished iteration loop -- now compute the estimate of sigma^2 # # based on the final GLS estimate. Use the "adjusted" version # # # ######################################################################### mu <- unweightfunc(time,bgls[1],bgls[2],bgls[3],bgls[4]) resid <- conc-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="indogls2.sout",append=T) sink("indogls2.sout",append=T) print(summary(indo.glsfit)) sink()