[R] Cumulative Points and Confidence Interval Manipulation in barplot2

[Marc Schwartz]

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,
> Bret

Bret,

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