#================================================================================================== # The following are certain R packages and file locations/destinations needed in the code below. # Note: We have left examples in the code. USER should make changes where # ######## CHANGE LOCATION ######## appears on the right-hand side of the code #================================================================================================== library(MASS) #-- Libraries needed when calling WinBUGS from R library(lattice) library(boot) library(coda) library(R2WinBUGS) #-- USER must specify where WinBUGS executable file is located BUGdir <- "C:/Program Files (x86)/WinBUGS14/" # <-- Example; comment out #BUGdir <- ".../WinBUGS14/" ######## CHANGE LOCATION ######## #-- USER must specify where file containing WinBUGS model is located #-- Leave "RP_Int_Est_Probit_Model_WinBUGS.odc" unchanged in the code below. BUGmod<-"I:/LD50 estimation/Simulation Study Paper/J Biopharm Stat submission/Submission 1 files/RP-Intervals-WinBUGS.odc" # <-- Example; comment out #BUGmod<-".../RP_Int_Est_Probit_Model_WinBUGS.odc" ######## CHANGE LOCATION ######## #-- USER must specify where MCMC results are to be stored WORKdir <- "I:/LD50 estimation/Simulation Study Paper/J Biopharm Stat submission/Submission 1 files/" # <-- Example; comment out #WORKdir <- ".../" ######## CHANGE LOCATION ######## #================================================================================================== # Creating the Example data, and creating objecsts that will be used later on. #================================================================================================== #-- The example data ExmplDat <- as.data.frame( matrix( c( 1, 0, 10, 6.30, 0.7992, 0, 6, 1, 0, 30, 6.87, 0.8371, 1, 6, 1, 0, 50, 7.30, 0.8633, 0, 6, 1, 0, 70, 7.75, 0.8895, 5, 6, 1, 0, 90, 8.46, 0.9274, 6, 6, 0, 1, 10, 7.56, 0.8784, 0, 6, 0, 1, 30, 8.25, 0.9163, 1, 6, 0, 1, 50, 8.76, 0.9425, 3, 6, 0, 1, 70, 9.31, 0.9687, 4, 6, 0, 1, 90, 10.15, 1.0066, 5, 6), 10,7,byrow=TRUE)) colnames(ExmplDat) <- c('CNTL','TRT','TargetP','X', 'logX', 'Y','N') #-- Printing the data ExmplDat #-- First, obtaining the data objects needed for performing the probit regressions CNTL <- ExmplDat[,1] #-- Identifier for 'Control' group TRT <- ExmplDat[,2] #-- Identifier for 'Treatment' group logX <- ExmplDat[,5] #-- log_10 of dose X Y <- ExmplDat[,6] #-- Number of dead animals N <- ExmplDat[,7] #-- Number of animals tested #================================================================================================== # Fitting the relative potency Model (Display 1): output --> alpha0, alpha1*, beta, and V #================================================================================================== #-- Fitting the DRF Model, and picking off the desired results model.rp <- glm( cbind(Y,N-Y) ~ TRT + logX , family=binomial(link='probit')) #-- Making objects out of the estimated coefficients alpha0 <- model.rp$coefficients[1] alpha1 <- model.rp$coefficients[2] beta <- model.rp$coefficients[3] #-- Obtaining the variance-covariance matrix of estimated model coefficients V <- vcov(model.rp) #-- Printing alpha0, alpha1*, beta, and V alpha0; alpha1; beta; V #-- MLE of theta = log(relative potency) theta <- -alpha1 / beta #-- Significance level & standard normal critical value signif.level <- 0.05 crit.z <- qnorm(1-signif.level/2) #================================================================================================== # Wald's method for "Margin of Error" (MoE) for log(RP) = theta #================================================================================================== MoE <- crit.z*sqrt( (theta^2)*( (V[2,2]/alpha1^2) + (V[3,3]/beta^2) - (2*V[2,3]/(alpha1*beta)))) WaldEst <- c( theta, theta-MoE, theta+MoE) #================================================================================================== # Fieller's method for "Margin of Error" for log(RP) = theta #================================================================================================== g <- -V[2,3]*crit.z^2 - alpha1*beta addend <- sqrt( g^2 - (beta^2 - V[3,3]*crit.z^2) * (alpha1^2 - V[2,2]*crit.z^2) ) d <- beta^2 - V[3,3]*crit.z^2 FiellerEst <- c(theta, (g - addend)/d, (g + addend)/d ) #================================================================================================== # LRT's method for CI of log(RP) = theta #================================================================================================== #---inverted loglikelihood function inv.lrt <- function(THETA0, n, x, r, y, maxloglik, conf.coef){ new.x <- x-THETA0*r new.m1 <- glm( cbind(y,n-y) ~ new.x, family=binomial(link='probit')) logLik(new.m1) - maxloglik+qchisq(conf.coef,1)/2 } maxLL <- logLik(model.rp) lrt.lo <- uniroot(inv.lrt, c(-999, theta), n=N, x=logX, r=TRT, y=Y, maxloglik=maxLL, 1-signif.level)$root lrt.up <- uniroot(inv.lrt, c(theta, 999), n=N, x=logX, r=TRT, y=Y, maxloglik=maxLL, 1-signif.level)$root LRTEst <- c(theta, lrt.lo, lrt.up ) #================================================================================================== # Bootstrap method for CI of log(RP) = theta # -- Note: with 5000 bootstrap samples, this takes no more than 30 seconds on a 64-bit desktop # computer having a 3.4 GHz processor # -- Note: due mainly to the small size of the example data set, there will likely be bootstrap # samples that fail to converge. In that instance, you will see a warning message from R. #================================================================================================== #-- Number of datalines in ExamplDat n.obs <- nrow(ExmplDat) #-- X-matrix for original data X.mat <- matrix( c(CNTL+TRT, TRT, logX), nrow(ExmplDat), 3, byrow=FALSE) #-- Estimated probits from original data probits <- X.mat %*% c(alpha0, alpha1, beta) #-- Desired number of bootstrapped data sets bootM <- 5000 #-- Matrix to keep desired results from analysis of bootstrapped data sets bootkeep <- matrix( rep(c(-999,-999,-999,-999,-999,-999),bootM), bootM, 6) #-- Performing the parametric bootstrap for(iter in 1:bootM){ Y.boot <- rbinom( n.obs, N, pnorm( probits ) ) #-- Parametrically generating Y data based on original data design boot.m1 <- glm( cbind(Y.boot,N-Y.boot) ~ TRT + logX , family=binomial(link='probit')) #-- Analyzing bootstrapped data bootkeep[iter,] <- c( boot.m1$coef[1:3], #-- Bootstrapped model estimates -boot.m1$coef[1]/boot.m1$coef[3], #-- mu0 -(boot.m1$coef[1] + boot.m1$coef[2]) /boot.m1$coef[3], #-- mu1 -boot.m1$coef[2]/boot.m1$coef[3]) #-- theta rm(Y.boot) } BootEst <- quantile(bootkeep[,6], probs = c(.5, signif.level/2, 1-signif.level/2)) #================================================================================================== # Bayes method log(RP) = theta # # NOTE: Due to the high thinning rate & large number of iterations, this section takes about # 2 minutes to run 64-bit desktop computer having a 3.4 GHz processor.# #================================================================================================== #-- Parameters (a.k.a. "nodes" in WinBUGS-speak) you want to estimate. nodes <- c('theta','mu0','mu1','alpha','beta') #-- Setting WinBUGS parameters N.Iterations <- 20000 #-- This is the number of MCMC iterations you want to #-- estimate your posterior distributions *after* allowing #-- for thinning (to reduce autocorrelation), "burn-in" #-- (to give the posterior distributions time to converge), #-- and the use of N.Chains of independent MCMCs. N.Chains <- 3 #-- Number of independent MCMC to estimate the posteriors. 1 chain is too few. N.Thin <- 100 #-- This "thinning rate" is very high. For this example data set, there is a # high degree of autocorrelation among posterior draws for the probit # regression parameters - alpha0, alpha1, & beta. N.Burnin <- 1000*N.Thin N.Iter <- ceiling(N.Iterations *N.Thin /N.Chains) + N.Burnin N.Save <- floor( (N.Iter - N.Burnin)*N.Chains / N.Thin ) #-- Initial values for *3* chains. This will need modifying if N.Chains <> 3. #-- Based off MLEs from classical probit regression BUGinits<-function(){c( list(beta=beta*0.98, alpha=c(alpha0,alpha1)/0.98), #-- Low-end of initial values list(beta=beta, alpha=c(alpha0,alpha1)), #-- Middle initial values list(beta=beta/0.98, alpha=c(alpha0,alpha1)*0.98) #-- High-end of initial values )} #-- Data needed for WinBUGS to perform the analysis BUGdata <- list( TRT=TRT, logX=logX, N=N, Y=Y, N.obs=NROW(Y), mu.beta=0, #-- Prior for beta (Gaussian mean) p.beta=0.001, #-- Prior for beta (Gaussian precision = 1/var) mu.alpha=0, #-- Prior for alpha (Gaussian mean) p.alpha=0.001) #-- Prior for alpha (Gaussian precision = 1/var) #-- Fitting the RP model with a Bayesian probit regression bayes.m1 <- bugs(data=BUGdata, inits=BUGinits, parameters.to.save=nodes, model.file=BUGmod, n.chains=N.Chains, n.iter=N.Iter, n.burnin=N.Burnin, n.thin=N.Thin, debug=TRUE, #-- debug=TRUE means you will have to manually quit # the WinBUGS program to obtain the results. Setting # this to FALSE will automate WinBUGS closure. bugs.directory=BUGdir, working.directory=WORKdir, codaPkg=TRUE) #-- Reads the posterior distributions from WinBUGS into R #-- This will increase/decrease if N.Chains increases/decreases from 3. chain123 <- rbind(read.coda(bayes.m1[1], paste(WORKdir,'codaIndex.txt', sep="/")), read.coda(bayes.m1[2], paste(WORKdir,'codaIndex.txt', sep="/")), read.coda(bayes.m1[3], paste(WORKdir,'codaIndex.txt', sep="/"))) #-- Computes the median and (signif.level/2) and (1-signif.level/2) quantiles from the posterior distributions BayesEst1 <- apply( chain123, 2, function(x) quantile( x ,probs=c(.5, signif.level/2, 1-signif.level/2) ) ) BayesEst <- BayesEst1[,7] #-- USER will need to look at the columns of "BayesEst1" to determine which column contains "theta". #================================================================================================== # Here we plot the different confidence intervals, centered on the maximum likelihood point estimate # of theta. #================================================================================================== #-- Plot preliminaries .x <- 0 ymin <- min(WaldEst[2],FiellerEst[2],LRTEst[2],BootEst[2],BayesEst[2]) ymax <- max(WaldEst[3],FiellerEst[3],LRTEst[3],BootEst[3],BayesEst[3]) del <- (ymax-ymin)/10 #-- Begin plot plot(1:5, seq(ymin,ymax,length.out=5), type='n', xlab='Method', ylab='theta estimate', xlim=c(0,5.5), ylim=c(ymin-del,ymax+del), axes=FALSE ) axis(1, at=1:5, label=c('Wald','Fieller', 'LRT', 'Bootstrap', 'Bayes')) axis(2, las=2) #-- Wald est0 <- WaldEst; .x <- .x+1 ; e <- est0[1]; l <- est0[2]; u <- est0[3] abline(h=e, col='red',lty=3) points(.x, e, pch=19); arrows(.x,l,.x,u, lwd=2, angle=90, length=.1, code=3); text(.x, l, pos=1, round(u-l,4)) #-- Fiellers est0 <- FiellerEst; .x <- .x+1 ; e <- est0[1]; l <- est0[2]; u <- est0[3] points(.x, e, pch=19); arrows(.x,l,.x,u, lwd=2, angle=90, length=.1, code=3); text(.x, l, pos=1, round(u-l,4)) #-- LRT est0 <- LRTEst; .x <- .x+1 ; e <- est0[1]; l <- est0[2]; u <- est0[3] points(.x, e, pch=19); arrows(.x,l,.x,u, lwd=2, angle=90, length=.1, code=3); text(.x, l, pos=1, round(u-l,4)) #-- Bootstrap est0 <- BootEst; .x <- .x+1 ; e <- est0[1]; l <- est0[2]; u <- est0[3] points(.x, e, pch=19); arrows(.x,l,.x,u, lwd=2, angle=90, length=.1, code=3); text(.x, l, pos=1, round(u-l,4)) #-- Bayes est0 <- BayesEst; .x <- .x+1 ; e <- est0[1]; l <- est0[2]; u <- est0[3] points(.x, e, pch=19); arrows(.x,l,.x,u, lwd=2, angle=90, length=.1, code=3); text(.x, l, pos=1, round(u-l,4)) #-- Figure keys text(0.5, e, pos=3,'ML Est', col='red') text(4, min(WaldEst[2],FiellerEst[2],LRTEst[2],BootEst[2],BayesEst[2]), 'Values are CI length')