[R] plotting results from tapply; using indexing to incorporate error bars

Marc Schwartz marc_schwartz at comcast.net
Mon Jan 29 18:42:26 CET 2007


Mike,

Using interaction.plot():

with(ex.dat, interaction.plot(x1, x2, response = y1,
                               type = "b", pch = 21,
                               col = c("red", "blue"),
                               ylim = range(ex.dat$y1)))

arrows(1:2, xbar + sem, 1:2, xbar - sem, 
       col = rep(c("red", "blue"), each = 2),
       angle = 90, code = 3, length = .05)


Note that you need to draw two pairs of error bars (4). By specifying
only two colors in the call to arrows() as you do below, they will be
recycled as:

  red, blue, red, blue

This is defined in ?arrows:

   The graphical parameters 'col', 'lty' and 'lwd' can be vectors of
   length greater than one and will be recycled if necessary.

What you actually want is:

  red, red, blue, blue

which you get with:

> rep(c("red", "blue"), each = 2)
[1] "red"  "red"  "blue" "blue"


Try the above and it should work.

HTH,

Marc

On Mon, 2007-01-29 at 12:12 -0500, Michael Rennie wrote:
> Hi, Mark
> 
> Thanks for the examples- this is great, and has helped me understand alot 
> more what's going on in the plotting functions.
> 
> Now that I'm trying to work error bars into this, I was curious if someone 
> might give me a hand indexing this properly so that I can format my error 
> bars to match the formatting of the grouping variables that match the lines 
> in the plot.
> 
> Here's a re-worked example from my first submission where I calculate 
> standard errors for plotting on the figure, plot the error bars, and try to 
> index them to match the appropriate lines. The values appear to be in the 
> right place (I've turned on the "legend" option for the interaction plot to 
> help verify this), however, my attempts at matching the colours on the bars 
> to the colours on the lines fails miserably (as you'll see if you execute 
> the code below). Is there any way to assign my colours to match them in a 
> way that makes more sense?
> 
> Thanks for your help so far,
> 
> Mike
> 
> y1<-rnorm(40, 2)
> x1<-rep(1:2, each=20)
> x2<-rep(1:2, each=10, times=2)
> 
> ex.dat<-data.frame(cbind(y1,x1,x2))
> 
> ex.dat$x1<-factor(ex.dat$x1, labels=c("A", "B"))
> ex.dat$x2<-factor(ex.dat$x2, labels=c("C", "D"))
> 
> attach(ex.dat)
> 
> xbar<-tapply(ex.dat$y1, ex.dat[,-1], mean)
> xbar
> s <- tapply(ex.dat$y1, ex.dat[,-1], sd)
> n <- tapply(ex.dat$y1, ex.dat[,-1], length)
> sem <- s/sqrt(n)
> sem
> 
> row.names(xbar)
> xbar[,1]
> 
> 
> 
> #from Marc Schwartz#
> 
> #par(mfcol = c(1, 3))
> 
> 
> options(graphics.record = TRUE)
> 
> #using plot
> 
> with(ex.dat, plot(1:2, xbar[, 1], ylim = range(y1),
>                    type = "b", col = "red",
>                    lty = c("dashed", "solid"),
>                    xaxt = "n", xlab = "x1",
>                    ylab = "mean of y1"))
> 
> 
> with(ex.dat, points(1:2, xbar[, 2], col = "blue",
>                      type = "b"))
> 
> 
> axis(1, at = 1:2, labels = c("A", "B"))
> 
> 
> #using matplot
> 
> matplot(1:2, xbar, col = c("red", "blue"),
>          pch = 21, type = "b", ylim = range(y1),
>          lty = c("dashed", "solid"),
>          xaxt = "n", xlab = "x1",
>          ylab = "mean of y1")
> 
> 
> axis(1, at = 1:2, labels = c("A", "B"))
> arrows(1:2,xbar+sem, 1:2,xbar-sem, col = c("red", "blue"), angle=90, 
> code=3, length=.1)
> 
> #using interaction.plot
> 
> with(ex.dat, interaction.plot(x1, x2, response = y1,
>                                type = "b", pch = 21,
>                                col = c("red", "blue"),
>                                ylim = range(ex.dat$y1)))
> 
> arrows(1:2,xbar+sem, 1:2,xbar-sem, col = c("red", "blue"), angle=90, 
> code=3, length=.05)
> 
> #as you can see, though the values for standard errors match the 
> appropriate means, the assignment of colours does not.
> 
> 
> 
> At 12:46 PM 26/01/2007, Marc Schwartz wrote:
> >On Fri, 2007-01-26 at 11:50 -0500, Michael Rennie wrote:
> > > Hi, there
> > >
> > > I'm trying to plot what is returned from a call to tapply, and can't 
> > figure
> > > out how to do it. My guess is that it has something to do with the
> > > inclusion of row names when you ask for the values you're interested in,
> > > but if anyone has any ideas on how to get it to work, that would be
> > > stellar.  Here's some example code:
> > >
> > > y1<-rnorm(40, 2)
> > > x1<-rep(1:2, each=20)
> > > x2<-rep(1:2, each=10 times=2)
> > >
> > > ex.dat<-data.frame(cbind(y1,x1,x2))
> > >
> > > ex.dat$x1<-factor(ex.dat$x1, labels=c("A", "B"))
> > > ex.dat$x2<-factor(ex.dat$x2, labels=c("C", "D"))
> > >
> > > attach(ex.dat)
> > >
> > > xbar<-tapply(ex.dat$y1, ex.dat[,-1], mean)
> > > xbar
> > >
> > > #values I'd like to plot:
> > > row.names(xbar) #levels of x1
> > > xbar[,1] #mean response of y1 for group C (x2) across x1
> > >
> > > #plot mean response y1 for group C against x1 (i.e., using x2 as a 
> > grouping
> > > factor).
> > > plot(row.names(xbar), xbar[,1], ylim=range[y1])
> > >
> > > #This is where things break down. The error message states that I need
> > > "finite xlim values" but I haven't assigned anything to xlim. If I just
> > > plot the data:
> > >
> > > plot(x1, y1)
> > >
> > > #This works fine.
> > >
> > > #And, I can get this to work:
> > >
> > > stripchart(xbar[1,]~row.names(xbar), vert=T)
> > >
> > > #However, I'd like to then add the values for the second set of means
> > > (i.e., mean values for group D against x1, or (xbar[,2])) to the plot.
> > > #I tried following up the previous command with:
> > >
> > > points(row.names(xbar), xbar[,2])
> > >
> > > #But that returns an error as well (NAs introduced by coercion).
> > >
> > >
> > >
> > > Any suggestions?
> > >
> > > Cheers,
> > >
> > > Mike
> > >
> > > PS- some of you might suggest for me to use interaction.plot, since 
> > this is
> > > essentially what I'm building here. True, but I can't get error bars using
> > > interaction.plot. I'm therefore trying to build my own version where I can
> > > specify the inclusion of error bars. Presumably the interaction.plot has
> > > figured out how to do what I'm attempting, so I have some faith that I am
> > > on the right track....
> >
> >Michael,
> >
> >The problem is that when you are using the rownames for 'xbar', they are
> >a character vector:
> >
> > > str(rownames(xbar))
> >  chr [1:2] "A" "B
> >
> >When you attempt to use the same values from 'ex.dat$x1', they are
> >factors, which are being coerced to their numeric equivalents, so that
> >they can work as x axis values.
> >
> > > str(ex.dat$x1)
> >  Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
> >
> > > as.numeric(ex.dat$x1)
> >  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2
> >[35] 2 2 2 2 2 2
> >
> >
> >
> >I might note using the following:
> >
> >par(mfcol = c(1, 3))
> >
> >with(ex.dat, plot(1:2, xbar[, 1], ylim = range(y1),
> >                   type = "b", col = "red",
> >                   lty = c("dashed", "solid"),
> >                   xaxt = "n", xlab = "x1",
> >                   ylab = "mean of y1"))
> >
> >with(ex.dat, points(1:2, xbar[, 2], col = "blue",
> >                     type = "b"))
> >
> >axis(1, at = 1:2, labels = c("A", "B"))
> >
> >
> >matplot(1:2, xbar, col = c("red", "blue"),
> >         pch = 21, type = "b", ylim = range(y1),
> >         lty = c("dashed", "solid"),
> >         xaxt = "n", xlab = "x1",
> >         ylab = "mean of y1")
> >
> >axis(1, at = 1:2, labels = c("A", "B"))
> >
> >
> >with(ex.dat, interaction.plot(x1, x2, response = y1,
> >                               type = "b", pch = 21,
> >                               col = c("red", "blue"),
> >                               ylim = range(ex.dat$y1),
> >                               legend = FALSE))
> >
> >
> >
> >You get basically the same 3 plots, with the principal difference in
> >interaction.plot() being a different x axis range.
> >
> >Using interaction.plot(), you can get the proper basic plot and then add
> >the CI's if you wish, since you know the x axis values of the mean
> >estimates, which is 1:NumberOfGroups (as in a boxplot).
> >
> >interaction.plot() actually uses matplot() internally.
> >
> >HTH,
> >
> >Marc Schwartz
> 
> Michael Rennie
> Ph.D. Candidate, University of Toronto at Mississauga
> 3359 Mississauga Rd. N.
> Mississauga, ON  L5L 1C6
> Ph: 905-828-5452  Fax: 905-828-3792
> www.utm.utoronto.ca/~w3rennie



More information about the R-help mailing list