######################################################################### # # # Program to implement the 3-step GLS algorithm with # # cmax total iterations, theta unknown, using the R # # function nls(), estimating theta by ID # # # # 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(t,W0,W1,beta0,beta1,tau0,tau1,v){ log2 <- log(2) t0 <- 20 f10 <- W0*(1-exp(-log2*(t/tau0)^beta0)) joinpt <- W0*(1-exp(-log2*(t0/tau0)^beta0)) f11 <- (t>t0)*(joinpt + W1*(1-exp(-log2*((t-t0)*(t>t0)/tau1)^beta1))) f1 <- (1-v)*f10 + v*( (t<=t0)*f10 + (t>t0)* f11 ) stuff <- (t+0.1*(t==t0)-t0)/tau1 piece <- log(abs(stuff)) # compute analytical dervivatives -- create the gradient matrix X meangrad <- array(0,c(length(t),6),list(NULL,c("W0","W1","beta0","beta1", "tau0","tau1" ))) meangrad[,"W0"] <- ( (1-v)*(1-exp(-log2*(t/tau0)^beta0)) + v*( (t<=t0)*(1-exp(-log2*(t/tau0)^beta0)) + (t>t0)*(1-exp(-log2*(t0/tau0)^beta0)) )) meangrad[,"W1"] <- v*(t>t0)*(1-exp(-log2*((t-t0)*(t>t0)/tau1)^beta1)) meangrad[,"beta0"] <- ( (1-v)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(log(t/tau0))*log2 + v*( (t<=t0)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(log(t/tau0))*log2 + (t>t0)*W0*exp(-log2*(t0/tau0)^beta0)*((t0/tau0)^beta0)*(log(t0/tau0))*log2 ) ) meangrad[,"beta1"] <- v*(t>t0)*W1*exp(-log2*((t-t0)*(t>t0)/tau1)^beta1)*(((t-t0)*(t>t0)/tau1)^beta1)*piece*log2 meangrad[,"tau0"] <- -( (1-v)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(beta0/tau0)*log2 + v*( (t<=t0)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(beta0/tau0)*log2 + (t>t0)*W0*exp(-log2*(t0/tau0)^beta0)*((t0/tau0)^beta0)*(beta0/tau0)*log2 ) ) meangrad[,"tau1"] <- -v*(t>t0)*W1*exp(-log2*((t-t0)*(t>t0)/tau1)^beta1)*(((t-t0)*(t>t0)/tau1)^beta1)*(beta1/tau1)*log2 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 on p. 1.6. 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 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(t,W0,W1,beta0,beta1,tau0,tau1,v){ log2 <- log(2) t0 <- 20 f10 <- W0*(1-exp(-log2*(t/tau0)^beta0)) joinpt <- W0*(1-exp(-log2*(t0/tau0)^beta0)) f11 <- (t>t0)*(joinpt + W1*(1-exp(-log2*((t-t0)*(t>t0)/tau1)^beta1))) f1 <- (1-v)*f10 + v*( (t<=t0)*f10 + (t>t0)* f11 ) f1 } # The transformed mean function (multiplied by the square root # of the current estimated weights, which are considered fixed) weightfunc <- function(t,W0,W1,beta0,beta1,tau0,tau1,v,wt){ w12 <- sqrt(wt) log2 <- log(2) t0 <- 20 # no vibration, v=0 f10 <- W0*(1-exp(-log2*(t/tau0)^beta0)) joinpt <- W0*(1-exp(-log2*(t0/tau0)^beta0)) f11 <- (t>t0)*(joinpt + W1*(1-exp(-log2*((t-t0)*(t>t0)/tau1)^beta1))) f1 <- ( (1-v)*f10 + v*( (t<=t0)*f10 + (t>t0)* f11 ) ) weightf <- f1*w12 # compute analytical dervivatives -- create the gradient matrix X stuff <- (t+0.1*(t==20)-t0)/tau1 piece <- log(abs(stuff)) # compute analytical dervivatives -- create the gradient matrix X weightgrad <- array(0,c(length(t),6),list(NULL,c("W0","W1","beta0","beta1", "tau0","tau1" ))) weightgrad[,"W0"] <- w12*( (1-v)*(1-exp(-log2*(t/tau0)^beta0)) + v*( (t<=t0)*(1-exp(-log2*(t/tau0)^beta0)) + (t>t0)*(1-exp(-log2*(t0/tau0)^beta0)) )) weightgrad[,"W1"] <- w12*v*(t>t0)*(1-exp(-log2*((t-t0)*(t>t0)/tau1)^beta1)) weightgrad[,"beta0"] <- w12*( (1-v)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(log(t/tau0))*log2 + v*( (t<=t0)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(log(t/tau0))*log2 + (t>t0)*W0*exp(-log2*(t0/tau0)^beta0)*((t0/tau0)^beta0)*(log(t0/tau0))*log2 ) ) weightgrad[,"beta1"] <- w12*v*(t>t0)*W1*exp(-log2*((t-t0)*(t>t0)/tau1)^beta1)*(((t-t0)*(t>t0)/tau1)^beta1)*piece*log2 weightgrad[,"tau0"] <- -w12*( (1-v)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(beta0/tau0)*log2 + v*( (t<=t0)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(beta0/tau0)*log2 + (t>t0)*W0*exp(-log2*(t0/tau0)^beta0)*((t0/tau0)^beta0)*(beta0/tau0)*log2 ) ) weightgrad[,"tau1"] <- -w12*v*(t>t0)*W1*exp(-log2*((t-t0)*(t>t0)/tau1)^beta1)*(((t-t0)*(t>t0)/tau1)^beta1)*(beta1/tau1)*log2 attr(weightf,"gradient") <- weightgrad weightf } ######################################################################### # # # To estimate theta, we use the "trick" of turning the ID 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 <- sqrt(abs(resid)*((mudot/mu)**theta)) trkgrad <- array(0,c(length(mu),1),list(NULL,c("theta"))) trkgrad[,"theta"] <- 0.5*trk*log(mudot/mu) attr(trk,"gradient") <- trkgrad # analytic derivative trk } ######################################################################### # # # Read in the data # # # ######################################################################### thedata <- scan("dissolution.dat") thedata <- matrix(thedata,ncol=3,byrow=T) outfile <- "diss_id.Rout" v <- thedata[,1] x <- thedata[,2] y <- thedata[,3] n <- length(x) p <- 6 ######################################################################### # # # Define list of starting values for beta and theta # # # ######################################################################### bstart <- list(W0=60,W1=30,beta0=1,beta1=1, tau0=10,tau1=10) tstart <- 1 ######################################################################### # # # Create the data frame for nls() # # # ######################################################################### indat <- data.frame(x,y) ######################################################################### # # # Specify the max number of iterations of GLS (C) # # # ######################################################################### cmax <- 20 ######################################################################### # # # 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. # # 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. # # # ######################################################################### ols.fit <- nls(y ~ meanfunc(x,W0,W1,beta0,beta1,tau0,tau1,v),indat,bstart) ######################################################################### # # # Extract the estimate from the object olsfit # # # ######################################################################### bols <- coef(ols.fit) mu <- unweightfunc(x,bols[1],bols[2],bols[3],bols[4],bols[5],bols[6],v) 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 DISSOLUTION DATA BY GLS",file=outfile,"\n", append=F) cat("WITH POWER VARIANCE, THETA UNKNOWN AND ESTIMATED BY IDENTITY",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){ ######################################################################### # # # 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(x,bgls[1],bgls[2],bgls[3],bgls[4],bgls[5],bgls[6],v) mudot <- prod(mu)**(1/n) mudot <- rep(mudot,n) resid <- y-mu dummy <- rep(0,length(x)) pldat <- data.frame(resid,mu,mudot) # Now estimate theta using ID "trick" thetafit <- nls(dummy ~ trkfunc(resid,mudot,mu,theta),pldat, list(theta=tstart)) theta <- coef(thetafit) cat("Iteration ",k,"\n",file=outfile,append=T) cat("Estimate of theta ",round(theta,4),"\n", file=outfile,append=T,"\n") mu <- unweightfunc(x,bgls[1],bgls[2],bgls[3],bgls[4],bgls[5],bgls[6],v) 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,W0,W1,beta0,beta1,tau0,tau1,v,wt),indat2,bstart, control=nls.control(maxiter=200)) ######################################################################### # # # 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("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],v) resid <- y-mu g <- mu^theta sigma2 <- sum((resid/g)**2)/(n-p) sigma <- sqrt(sigma2) ######################################################################### # # # Print out the final estimatea of sigma and theta # # # ######################################################################### cat("Final estimate of sigma ",round(sigma,6),"\n","\n", file=outfile,append=T) cat("Final estimate of theta ",round(theta,6),"\n","\n", file=outfile,append=T) predols <- unweightfunc(x,bols[1],bols[2],bols[3],bols[4],bols[5],bols[6],v) residols <- y-predols sink(outfile,append=T) print(summary(gls.fit)) sink()