[R] interaction plot with error bars based on TukeyHSD intervals

[Manuel López-Ibáñez]

Hello,

maybe the previous message was not very clear. Thus, I will try to 
explain the problem in a better way.

I would like to have an interaction plot with error bars based on 
TukeyHSD intervals.

In an experiment with two factors "A" and "B", each of them with two 
levels "A={a1,a2}" and "B={b1,b2}", I would like to do a plot like the 
following:


<attached image prueba.png>



I calculate the corresponding interval with:

tk <- TukeyHSD(aov(X$response ~ (factor(X$A) +  factor(X$B))2, data=X) , 
conf.level=0.95)

HSDfactor <- 
max(abs(tk$"factor(X$A):factor(X$B)"[,2]-tk$"factor(X$A):factor(X$B)"[,3]))


Finally, I modified interaction.plot() to add error bars with:

......................................................................
++ ylim <- c(min(cells)-(HSDfactor*0.5), max(cells)+(HSDfactor*0.5))
                                                                                                                                              


  matplot(xvals, cells, ..., type = type,  xlim = xlim, ylim = ylim,
          xlab = xlab, ylab = ylab, axes = axes, xaxt = "n",
          col = col, lty = lty, pch = pch)
                                                                                                                                              


++  ly <- cells[,1]+(HSDfactor*0.5)
++   uy <- cells[,1]-(HSDfactor*0.5)
                                                                                                                                              


++  errbar(xvals,cells[,1],ly,uy,add=TRUE, lty=3, cap=0, lwd=2)
                                                                                                                                             

++  ly <- cells[,2]+(HSDfactor*0.5)
++   uy <- cells[,2]-(HSDfactor*0.5)
                                                                                                                                             

++  errbar(xvals,cells[,2],ly,uy,add=TRUE, lty=3, cap=0, lwd=2)

 if(legend) {
        yrng <- diff(ylim)
        yleg <- ylim[2] - 0.1 * yrng

.............................................................................


However, the resulting intervals are much bigger than they should be, 
thus I think I did something wrong.

I have searched on Google, CRAN and the R mail archive but I only found 
related problems [1, 2], but not any answer...

Andrew Robinson [1] shows that there is a problem when using TukeyHSD 
for interaction terms. However, nobody suggested any work-around.

In another different thread, Martin Henry H. Stevens [2] suggests a 
work-around for the problem, but he thinks that the results obtained are 
"slightly inaccurate". As well, there is not feedback to solve the problem.


Any idea?

Thank you very much.

Manuel.

[1] http://finzi.psych.upenn.edu/R/Rhelp02/archive/12926.html
[2] http://finzi.psych.upenn.edu/R/Rhelp02/archive/32849.html



-------------- next part --------------
A non-text attachment was scrubbed...
Name: prueba.png
Type: image/png
Size: 18151 bytes
Desc: not available
Url : https://stat.ethz.ch/pipermail/r-help/attachments/20040611/5a021a43/prueba.png