Bret Collier
2005-Apr-12 15:14 UTC
[R] Cumulative Points and Confidence Interval Manipulation in barplot2
R-Users, I am working with gplots (in gregmisc bundle) plotting some posterior probabilities (using barplot2) of harvest bag limits for discrete data (x-axis from 0 to 12, data is counts) and I ran into a couple of questions whose solutions have evaded me. 1) When I create and include the confidence intervals, the lower bound of the confidence intervals for several of the posterior probabilities is below zero, and in those specific cases I only want to show the upper limit for those CI's so they do not extend below the x-axis (as harvest can not be <0). Also, comments on a better technique for CI construction when the data is bounded to be >=0 would be appreciated. 2) I would also like to show the cumulative probability (as say a point or line) across the range of the x-axis on the same figure at the top, but I have been unable to figure out how to overlay a set of cumulative points over the barplot across the same range as the x-axis. Below is some example code showing the test data I am working on (xzero): xzero <- table(factor(WWNEW[HUNTTYPE=="DOVEONLY"], levels=0:12))> xzero0 1 2 3 4 5 6 7 8 9 10 11 12 179 20 9 2 2 0 1 0 0 0 0 0 0> n <- sum(xzero) > k <- sum(table(xzero)) > meantheta1 <-((2*xzero + 1)/(2*n + k)) > vartheta1<-((2*(((2*n)+k)-((2*xzero)+1)))*((2*xzero)+1))/((((2*n)+k)^2)*(((2*n)+k)+2))> stderr <- sqrt(vartheta1) > cl.l <- meantheta1-(stderr*2)#Fake CI: Test > cl.u <- meantheta1+(stderr*2)#Fake CI: Test > barplot2(meantheta1, xlab="WWD HARVEST DOVE ONLY 2001",ylab="Probability", ylim=c(0, 1),xpd=F, col="blue", border="black", axis.lty=1,plot.ci=TRUE, ci.u = cl.u, ci.l = cl.l)> title(main="WHITE WING DOVE HARVEST PROBABILITIES: DOVE HUNT ONLY")I would greatly appreciate any direction or assistance, Thanks, Bret platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major 2 minor 0.1 year 2004 month 11 day 15 language R *Note: I am working in Exmacs
Marc Schwartz
2005-Apr-12 16:26 UTC
[R] Cumulative Points and Confidence Interval Manipulation in barplot2
On Tue, 2005-04-12 at 10:14 -0500, Bret Collier wrote:> R-Users, > I am working with gplots (in gregmisc bundle) plotting some posterior > probabilities (using barplot2) of harvest bag limits for discrete data > (x-axis from 0 to 12, data is counts) and I ran into a couple of > questions whose solutions have evaded me. > > 1) When I create and include the confidence intervals, the lower bound > of the confidence intervals for several of the posterior probabilities > is below zero, and in those specific cases I only want to show the upper > limit for those CI's so they do not extend below the x-axis (as harvest > can not be <0). Also, comments on a better technique for CI > construction when the data is bounded to be >=0 would be appreciated. > > 2) I would also like to show the cumulative probability (as say a > point or line) across the range of the x-axis on the same figure at the > top, but I have been unable to figure out how to overlay a set of > cumulative points over the barplot across the same range as the x-axis. > > Below is some example code showing the test data I am working on > (xzero): > > xzero <- table(factor(WWNEW[HUNTTYPE=="DOVEONLY"], levels=0:12)) > > xzero > > 0 1 2 3 4 5 6 7 8 9 10 11 12 > 179 20 9 2 2 0 1 0 0 0 0 0 0 > > > n <- sum(xzero) > > k <- sum(table(xzero)) > > meantheta1 <-((2*xzero + 1)/(2*n + k)) > > vartheta1 > <-((2*(((2*n)+k)-((2*xzero)+1)))*((2*xzero)+1))/((((2*n)+k)^2)*(((2*n)+k)+2)) > > stderr <- sqrt(vartheta1) > > cl.l <- meantheta1-(stderr*2)#Fake CI: Test > > cl.u <- meantheta1+(stderr*2)#Fake CI: Test > > barplot2(meantheta1, xlab="WWD HARVEST DOVE ONLY 2001", > ylab="Probability", ylim=c(0, 1),xpd=F, col="blue", border="black", > axis.lty=1,plot.ci=TRUE, ci.u = cl.u, ci.l = cl.l) > > title(main="WHITE WING DOVE HARVEST PROBABILITIES: DOVE HUNT ONLY") > > > I would greatly appreciate any direction or assistance, > Thanks, > BretBret, If you replace the lower bound of your confidence intervals as follows, you can get just the upper bound plotted: cl.l.new <- ifelse(cl.l >= 0, cl.l, meantheta1) This will set the lower bound to meantheta1 in those cases, thus plotting the upper portion and you can remove the 'xpd=F' argument. Use 'ci.l = cl.l.new' here: barplot2(meantheta1, xlab="WWD HARVEST DOVE ONLY 2001", ylab="Probability", ylim=c(0, 1), col="blue", border="black", axis.lty=1,plot.ci=TRUE, ci.u = cl.u, ci.l = cl.l.new) I would defer to others with more Bayesian experience on alternatives for calculating bounded CI's for the PP's. With respect to the cumulative probabilities, if I am picturing the same thing you are, you can use the cumsum() function and then add points and/or a line as follows: points(cumsum(meantheta1), pch = 19) lines(cumsum(meantheta1), lty = "solid") See ?cumsum, ?points and ?lines for more information. BTW, some strategically placed spaces would help make your code a bit more readable for folks. HTH, Marc Schwartz