######################################################################### # # # Program to implement the 3-step GLS algorithm with # # cmax total iterations, theta known, using the Splus # # function nls() # # # # Use it here to fit logistic regression # # # # 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) # # # # 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,b3,b4,b5,b6,b7,b8){ temp <- exp(b1 + b2*x[,1] + b3*x[,2] + b4*x[,3] + b5*x[,4] + b6*x[,5] + b7*x[,6] + b8*x[,7]) fpl <- temp/(1+temp) # compute analytical dervivatives -- create the gradient matrix X meangrad <- array(0,c(nrow(x),8),list(NULL,c("b1","b2","b3","b4","b5", "b6","b7","b8"))) meangrad[,"b1"] <- fpl*(1-fpl) meangrad[,"b2"] <- fpl*(1-fpl)*x[,1] meangrad[,"b3"] <- fpl*(1-fpl)*x[,2] meangrad[,"b4"] <- fpl*(1-fpl)*x[,3] meangrad[,"b5"] <- fpl*(1-fpl)*x[,4] meangrad[,"b6"] <- fpl*(1-fpl)*x[,5] meangrad[,"b7"] <- fpl*(1-fpl)*x[,6] meangrad[,"b8"] <- fpl*(1-fpl)*x[,7] attr(fpl,"gradient") <- meangrad fpl } ######################################################################### # # # 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,b1,b2,b3,b4,b5,b6,b7,b8){ temp <- exp(b1 + b2*x[,1] + b3*x[,2] + b4*x[,3] + b5*x[,4] + b6*x[,5] + b7*x[,6] + b8*x[,7]) fpl <- temp/(1+temp) fpl } # The transformed mean function (multiplied by the square root # of the current estimated weights, which are considered fixed) weightfunc <- function(x,b1,b2,b3,b4,b5,b6,b7,b8,wt){ temp <- exp(b1 + b2*x[,1] + b3*x[,2] + b4*x[,3] + b5*x[,4] + b6*x[,5] + b7*x[,6] + b8*x[,7]) pred1 <- temp/(1+temp) w12 <- sqrt(wt) fpl <- pred1*w12 # compute analytical dervivatives -- create the gradient matrix X weightgrad <- array(0,c(nrow(x),8),list(NULL,c("b1","b2","b3","b4","b5", "b6","b7","b8"))) weightgrad[,"b1"] <- pred1*(1-pred1)*w12 weightgrad[,"b2"] <- pred1*(1-pred1)*x[,1]*w12 weightgrad[,"b3"] <- pred1*(1-pred1)*x[,2]*w12 weightgrad[,"b4"] <- pred1*(1-pred1)*x[,3]*w12 weightgrad[,"b5"] <- pred1*(1-pred1)*x[,4]*w12 weightgrad[,"b6"] <- pred1*(1-pred1)*x[,5]*w12 weightgrad[,"b7"] <- pred1*(1-pred1)*x[,6]*w12 weightgrad[,"b8"] <- pred1*(1-pred1)*x[,7]*w12 attr(fpl,"gradient") <- weightgrad fpl } ######################################################################### # # # The data # # # ######################################################################### # data thedata <- scan("trial.dat") thedata <- matrix(thedata,ncol=8,byrow=T) # need to create the 2 dummy variables for smoking y <- thedata[,2] covs <- thedata[,3:8] x3 <- covs[,3]==1 x4 <- covs[,3]==2 n <- length(y) p <- 8 # form the appropriate vector x as in problem x <- cbind(covs[,1],covs[,2],x3,x4,covs[,4:6]) # starting value for beta xstar <- cbind(rep(1,n),x) ystar <- (y+1/2)/2 wstar <- diag(ystar*(1-ystar)) eta <- log(ystar/(1-ystar)) zstar <- eta + (y-ystar)/(ystar*(1-ystar)) bs <- solve(t(xstar)%*%wstar%*%xstar)%*%t(xstar)%*%wstar%*%zstar bstart = list(b1=bs[1], b2=bs[2], b3=bs[3], b4=bs[4], b5=bs[5], b6=bs[6], b7=bs[7], b8=bs[8]) outfile <- "trial.gls" ######################################################################### # # # Create the data frame for nls() # # # ######################################################################### thedat <- data.frame(x,y) ######################################################################### # # # Specify the max number of iterations of GLS (C) # # # # Alternatively, we could check for convergence after step (iii) # # # ######################################################################### cmax <- 50 ######################################################################### # # # 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. # # # ######################################################################### # Note that we have used the nls.control option to change the # convergence criterion to something more stringent than the # default of 1e-5 olsfit <- nls(y ~ meanfunc(x,b1,b2,b3,b4,b5,b6,b7,b8),thedat,bstart, nls.control(tol=1e-8)) ######################################################################### # # # 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 TRIAL 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){ # alternatively can use the most recent iterate as the starting value # bstart <- list(b1=bgls[1], b2=bgls[2], b3=bgls[3], b4=bgls[4], # b5=bgls[5], b6=bgls[6], b7=bgls[7], b8=bgls[8]) ######################################################################### # # # Step (ii) -- calculate the weights and the transformed response # # to use in the WLS calculation in (iii) # # # ######################################################################### mu <- unweightfunc(x,bgls[1],bgls[2],bgls[3],bgls[4],bgls[5],bgls[6], bgls[7],bgls[8]) wt <- 1/(mu*(1-mu)) 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) glsfit <- nls(ywt ~ weightfunc(x,b1,b2,b3,b4,b5,b6,b7,b8,wt),thedat2, bstart,nls.control(tol=1e-8)) ######################################################################### # # # 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,bgls[1],bgls[2],bgls[3],bgls[4],bgls[5],bgls[6], bgls[7],bgls[8]) resid <- y-mu ######################################################################### # # # Print out summary provided by the nls() function for the final # # iteration # # # ######################################################################### sink(outfile,append=T) print(summary(glsfit)) sink()