############################################################################### # NAME: Chris Bilder # # DATE: 9-5-03, Updated last on 10-29-04 # # PURPOSE: Calculate coverage, RB, and RE for Bilder and Tebbs (2004) paper # # # # NOTES: 1) See MAIN PROGRAM for settings of p, s, n, and alpha # # 2) Like RB_RE.R, this program also calculates the RB and RE for the # # empirical Bayes estimators! The reason for two different programs# # is the HPD calculations add a lot of time here. You can use the # # RB_RE.R when you are not interested in coverage or if you want # # RB and RE for the Burrows estimator. # # # ############################################################################### #Find start time start.time<-proc.time() ########################################################################### #Calculate the coverage and other summary measures # ########################################################################### check.in2<-function(ci, prob, p, prob.out){ save<-ifelse(p>ci$low, ifelse(pci$up])/(1-prob.out) #check if p is ABOVE upper limit low<-sum(prob[p mat.sort[mat.sort[, 2] == max(mat.sort[, 2]), ][[1]], ] matlowq <- cbind(matlow.sort[, 1], abs(matlow.sort[, 2] - q)) matupq <- cbind(matup.sort[, 1], abs(matup.sort[, 2] - q)) matlowq.sort <- matlowq[sort.list(matlowq[, 2]), ] matupq.sort <- matupq[sort.list(matupq[, 2]), ] low <- matlowq.sort[, 1][[1]] up <- matupq.sort[, 1][[1]] c(low, up) } ########################################################################### # Posterior distribution # ########################################################################### f.p<-function(p, n, s, t, beta) { s*exp(lgamma(n+beta/s+1) - lgamma(n-t+beta/s) - lgamma(t+1)) * (1-p)^(s*(n-t)+beta-1) * (1-(1-p)^s)^t } ########################################################################### # EB Interval # ########################################################################### eb<-function(info, beta) { p<-info[1] s<-info[2] n.group<-info[3] alpha<-info[4] theta<-1-(1-p)^s prob<-dbinom(1:(n.group-1), n.group, theta) prob.out<-dbinom(0, n.group, theta)+dbinom(n.group, n.group, theta) t<-1:(n.group-1) #Find the EB estimates of p and find bias and mse p.eb1<-1-exp(lgamma(n.group-t+beta/s+1/s) + lgamma(n.group+beta/s+1) - lgamma(n.group+beta/s+1+1/s) - lgamma(n.group-t+beta/s)) p.eb2<-1-(1-(t+1)/(n.group+beta/s+1))^(1/s) bias.mse.eb1<-bias.mse(p.eb1, p, prob) bias.mse.eb2<-bias.mse(p.eb2, p, prob) ############################################################ #MOM - could do MOM outside of function beta.mom<-s*(n.group-t)/t p.eb3<-1-exp(lgamma(n.group-t+beta.mom/s+1/s) + lgamma(n.group+beta.mom/s+1) - lgamma(n.group+beta.mom/s+1+1/s) - lgamma(n.group-t+beta.mom/s)) bias.mse.eb3<-bias.mse(p.eb3, p, prob) ############################################################ #Get HPD hpd.all<-matrix(NA, length(t), 5) #t, p<-low, p<-up, f(p)<-low, f(p)<-up hpd.all.mom<-matrix(NA, length(t), 5) for (a in 1:length(t)) { U<-rbeta(10000, t[a]+1, n.group-t[a]+beta[a]/s) p.obs<-1-(1-U)^(1/s) f<-f.p(p.obs, n.group, s, t[a], beta[a]) hpd<-HPD.func(p.obs,f,alpha) f.low<-f.p(hpd[1], n.group, s, t[a], beta[a]) #Use as check - should be about same as f.up f.up<-f.p(hpd[2], n.group, s, t[a], beta[a]) hpd.all[a,]<-c(t[a], hpd, f.low, f.up) #Using MOM #U<-rbeta(10000, t[a]+1, n.group-t[a]+beta.mom[a]/s) #p.obs<-1-(1-U)^(1/s) #f<-f.p(p.obs, n.group, s, t[a], beta.mom[a]) #hpd<-HPD.func(p.obs,f,alpha) #f.low<-f.p(hpd[1], n.group, s, t[a], beta.mom[a]) #Use as check - should be about same as f.up #f.up<-f.p(hpd[2], n.group, s, t[a], beta.mom[a]) #hpd.all.mom[a,]<-c(t[a], hpd, f.low, f.up) } ci<-data.frame(low=hpd.all[,2], up=hpd.all[,3]) save.hpd<-check.in2(ci, prob, p, prob.out) ci<-data.frame(low=hpd.all.mom[,2], up=hpd.all.mom[,3]) save.hpd.mom<-check.in2(ci, prob, p, prob.out) #Get Credible interval low<-1-(1-qbeta(alpha/2, t+1, n.group-t+beta/s))^(1/s) up<-1-(1-qbeta(1-alpha/2, t+1, n.group-t+beta/s))^(1/s) ci<-data.frame(low=low, up=up) save.it<-check.in2(ci, prob, p, prob.out) save.it2<-c(save.it[1:3], bias.mse.eb1, bias.mse.eb2, bias.mse.eb3, save.hpd[1:3], save.hpd.mom[1:3], save.it[4:5], save.hpd[4:5]) save.it2 } ########################################################################### #Individual coverage plot with Wald # ########################################################################### plot.it2<-function(title, name) { wald.x<-switch(name, "p"=save.all$p, "n"=save.all$n.group, "s"=save.all$s) plot(wald.x, save.all$coverage, ylim=c(0.85,1), type="n", main=paste(title ), xlab=name, ylab="Coverage", lab=c(10,10)) lines(wald.x, save.all$coverage, type="l", lty=1, col=1, lwd=1) abline(h = 1-alpha, lty=2, col=1) invisible() } ########################################################################### #Control the vary p plots # ########################################################################### plot.vary<-function(name, n.use, s.use, alpha=alpha) { #Range for P(T=0) and P(T=n) > 0.01 low<-1-(0.01)^(1/(s.use*n.use)) up<-1-(1-0.01^(1/n.use))^(1/s.use) #Coverage plots - these will always be "p" now - originally set up program to vary n and s as well wald.x<-switch(name, "p"=save.all$p, "n"=save.all$n.group, "s"=save.all$s) eb.x<-switch(name, "p"=save.all.eb$p, "n"=save.all.eb$n.group, "s"=save.all.eb$s) win.graph(width = 7, height = 7, pointsize = 12) par(mfrow=c(2,1)) plot.it2("Equal-tail", name) lines(eb.x, save.all.eb$coverage, type="l", lty=1, col=4, lwd=3, panel.first = grid(col = "black", lty = "dotted")) abline(v=low, lty=5) abline(v=up, lty=5) plot.it2("HPD", name) lines(eb.x, save.all.eb$coverage.hpd, type="l", lty=1, col=2, lwd=3, panel.first = grid(col = "black", lty = "dotted")) abline(v=low, lty=5) abline(v=up, lty=5) #Mean length plot win.graph(width = 7, height = 7, pointsize = 12) par(mfrow=c(1,1)) plot(wald.x, save.all$mean.length, type="n", main=paste("Mean length plot with n=", n.use, "and s=", s.use), xlab="p", ylab="Mean length", lab=c(10,10)) lines(wald.x, save.all$mean.length, type="l", lty=1, col=1) lines(eb.x, save.all.eb$mean.length, type="l", lty=1, col=4) lines(eb.x, save.all.eb$mean.length.hpd, type="l", lty=1, col=2) lines(eb.x, save.all.eb$mean.length.hpd.mom, type="l", lty=1, col=5) names<-c("Wald", "Equal-tail", "HPD") legend(x=0.05, y=0.03, names, lty=c(1, 1, 1), col=c(1, 3, 6)) #RB and RE plots win.graph(width = 7, height = 7, pointsize = 10) par(mfrow=c(2,1)) min.y<-min(save.all$bias.mle/save.all.eb$bias.eb1,save.all$bias.mle/save.all.eb$bias.eb2,save.all$bias.mle/save.all.eb$bias.eb3) max.y<-max(save.all$bias.mle/save.all.eb$bias.eb1,save.all$bias.mle/save.all.eb$bias.eb2,save.all$bias.mle/save.all.eb$bias.eb3) plot(save.all$p, save.all$bias.mle/save.all.eb$bias.eb2, type="n", main=paste("n=", n.use, "and s=", s.use), xlab="p", ylab="Relative bias", lab=c(10,10), ylim=c(min.y, max.y), panel.first=grid(col = 1)) lines(save.all$p, save.all$bias.mle/save.all.eb$bias.eb1, type="l", lty=1, col=2) lines(save.all$p, save.all$bias.mle/save.all.eb$bias.eb2, type="l", lty=1, col=4) lines(save.all$p, save.all$bias.mle/save.all.eb$bias.eb3, type="l", lty=1, col="green") abline(h = 1, lty=1, col=1) min.y<-min(save.all$mse.mle/save.all.eb$mse.eb1,save.all$mse.mle/save.all.eb$mse.eb2,save.all$mse.mle/save.all.eb$mse.eb3) max.y<-max(save.all$mse.mle/save.all.eb$mse.eb1,save.all$mse.mle/save.all.eb$mse.eb2,save.all$mse.mle/save.all.eb$mse.eb3) plot(save.all$p, save.all$mse.mle/save.all.eb$mse.eb1, type="n", main=paste("n=", n.use, "and s=", s.use), xlab="p", ylab="Relative efficiency", lab=c(10,10), ylim=c(min.y, max.y), panel.first=grid(col = 1)) lines(save.all$p, save.all$mse.mle/save.all.eb$mse.eb1, type="l", lty=1, col=2) lines(save.all$p, save.all$mse.mle/save.all.eb$mse.eb2, type="l", lty=1, col=4) lines(save.all$p, save.all$mse.mle/save.all.eb$mse.eb3, type="l", lty=1, col="green") abline(h = 1, lty=1, col=1) names<-c("EB1", "EB2", "EB3") legend(x=0.005, y=min.y+(max.y-min.y)/2, names, lty=c(1, 1), col=c(2, 4, "green"), bg="white") #Where interval miss plot win.graph(width = 7, height = 7, pointsize = 10) par(mfrow=c(2,2)) plot(save.all$p, save.all$low, type="l", main=paste("Wald C.I. miss, n=", n.use, "and s=", s.use), xlab="p", ylab="Probability", lab=c(10,10), ylim=c(0.0, 0.15), col=3, panel.first = grid(col = "black", lty = "dotted")) lines(save.all$p, save.all$high, type="l", lty=1, col=2) legend(x=0.07, y=0.15, c("low", "high"), lty=c(1, 1), col=c(3, 2), bg="white") #equal tail plot(save.all.eb$p, save.all.eb$low.eq, type="l", main=paste("Equal-tail miss, n=", n.use, "and s=", s.use), xlab="p", ylab="Probability", lab=c(10,10), ylim=c(0.0, 0.15), col=3, panel.first = grid(col = "black", lty = "dotted")) lines(save.all.eb$p, save.all.eb$high.eq, type="l", lty=1, col=2) legend(x=0.07, y=0.15, c("low", "high"), lty=c(1, 1), col=c(3, 2), bg="white") #HPD plot(save.all.eb$p, save.all.eb$low.hpd, type="l", main=paste("HPD miss, n=", n.use, "and s=", s.use), xlab="p", ylab="Probability", lab=c(10,10), ylim=c(0.0, 0.15), col=3, panel.first = grid(col = "black", lty = "dotted")) lines(save.all.eb$p, save.all.eb$high.hpd, type="l", lty=1, col=2) legend(x=0.07, y=0.15, c("low", "high"), lty=c(1, 1), col=c(3, 2)) invisible() } ########################################################################### #MAIN PROGRAM # ########################################################################### #Settings p<-seq(0.001, 0.10, 0.0005) #Settings for paper p<-seq(0.01, 0.10, 0.02) #Try these settings first to see how long it will take to run s<-10 n.group<-40 alpha<-c(0.05) #Put all settings in one matrix sim.set<-expand.grid(p=p, s=s, n.group=n.group, alpha=alpha) sim.set.mat<-as.matrix.data.frame(sim.set) #Wald res<-apply(sim.set.mat, 1, wald) save.all<-data.frame(sim.set.mat, coverage=round(res[1,],4), mean.length=round(res[2,],4), ck.zero=round(res[3,],4), bias.mle=res[4,], mse.mle=res[5,], low=res[6,], high=res[7,]) #EB #Derivative of marginal distribution root.func<-function(beta, t, s, n.group) { 1/beta + digamma(n.group - t + beta/s)/s - digamma(n.group + beta/s + 1)/s } #Find beta vector outside of loop since do not need p. #WARNING: If changed program for vary over different s and n values, will need to do this inside of EB function #If beta > upper bound, a warning message will be produced! beta<-numeric(length(1:(n.group-1))) for (t in 1:(n.group-1)) { beta[t]<-uniroot(root.func, c(0.00001, 10000), t=t, s=s, n.group=n.group, tol=10^(-8))$root } res.eb<-apply(sim.set.mat, 1, eb, beta) save.all.eb<-data.frame(sim.set.mat, coverage=round(res.eb[1,],4), mean.length=round(res.eb[2,],4), ck.zero=round(res.eb[3,],4), bias.eb1=res.eb[4,], mse.eb1=res.eb[5,], bias.eb2=res.eb[6,], mse.eb2=res.eb[7,], bias.eb3=res.eb[8,], mse.eb3=res.eb[9,], coverage.hpd=round(res.eb[10,],4), mean.length.hpd=round(res.eb[11,],4), coverage.hpd.mom=round(res.eb[13,],4), mean.length.hpd.mom=round(res.eb[14,],4), low.eq=res.eb[16,], high.eq=res.eb[17,],low.hpd=res.eb[18,], high.hpd=res.eb[19,] ) #Plots for varying p - adjustments to the program could vary n and s as well vary<-"p" plot.vary(vary, n.group, s) #Save results for bias and mse bias<-data.frame(p=save.all$p, n=rep(factor(n.group), length(save.all$p)), s=rep(factor(s), length(save.all$p)), mle=save.all$bias.mle, eb1=save.all.eb$bias.eb1, eb2=save.all.eb$bias.eb2, eb3=save.all.eb$bias.eb3, rb.eb1=save.all$bias.mle/save.all.eb$bias.eb1, rb.eb2=save.all$bias.mle/save.all.eb$bias.eb2, rb.eb3=save.all$bias.mle/save.all.eb$bias.eb3) mse<-data.frame(p=save.all$p, n=rep(factor(n.group), length(save.all$p)), s=rep(factor(s), length(save.all$p)), mle=save.all$mse.mle, eb1=save.all.eb$mse.eb1, eb2=save.all.eb$mse.eb2, eb3=save.all.eb$mse.eb3, re.eb1=save.all$mse.mle/save.all.eb$mse.eb1, re.eb2=save.all$mse.mle/save.all.eb$mse.eb2, re.eb3=save.all$mse.mle/save.all.eb$mse.eb3) epsilon<-0.00001 print.bias<-bias[(bias$p>0.01-epsilon & bias$p<0.01+epsilon) | (bias$p>0.02-epsilon & bias$p<0.02+epsilon) | (bias$p>0.05-epsilon & bias$p<0.05+epsilon) | (bias$p>0.10-epsilon & bias$p<0.10+epsilon),] print.bias print.mse<-mse[(mse$p>0.01-epsilon & mse$p<0.01+epsilon) | (mse$p>0.02-epsilon & mse$p<0.02+epsilon) | (mse$p>0.05-epsilon & mse$p<0.05+epsilon) | (mse$p>0.10-epsilon & mse$p<0.10+epsilon),] print.mse #Summary measure - what method is closest closeness3<-function(coverage) { round(mean(abs(coverage-(1-alpha))),6) } wald.close<-closeness3(save.all$coverage) eb.close<-closeness3(save.all.eb$coverage) eb.hpd.close<-closeness3(save.all.eb$coverage.hpd) names<-c("Wald", "Empirical Bayes equal tail", "HPD") close.measures<-c(wald.close, eb.close, eb.hpd.close) save.close<-data.frame(names, close.measures) save.close[order(close.measures),] #Percent >= 1-alpha wald.above<-mean(save.all$coverage>(1-alpha)) eb.above<-mean(save.all.eb$coverage>(1-alpha)) eb.hpd.above<-mean(save.all.eb$coverage.hpd>(1-alpha)) names<-c("Wald", "Empirical Bayes equal tail", "HPD") above.measures<-c(wald.above, eb.above, eb.hpd.above) save.above<-data.frame(names, above.measures) save.above[rev(order(above.measures)),] #Find end time and total time elapsed end.time<-proc.time() save.time<-end.time-start.time cat("\n Number of minutes running:", save.time[3]/60, "\n \n")