Hi,

The problem is that I would like to do an interaction plot with 
intervals based on Tukey's honestly significant difference (HSD) 
procedure, but I do not know how to do it in R.

I have 3 factors "A", "B" and "C" and a response variable "response".
I would like to study a model where there are main effects and second 
order interaction effects.

For instance, for factors B and C, I would like to plot a line for each 
level of C connecting the least square means for every level of B. 
Additionally, I would like to include the 95,0% HSD intervals for
the means in such a way that any two intervals which do not overlap 
correspond to a
pair of means which have a statistically significant difference.

For comparison, in STATGRAPHICS Plus, I choose Analysis of Variance -> 
Multifactor ANOVA...
Then I plot the second order interactions and I have the option to plot 
intervals. I choose to plot Tukey HSD intervals at 95% and I obtain an 
interaction plot like using interaction.plot on R, but addittionally an 
interval is plotted on every point.

How can I do this on R?

I have already looked in Google, in the mail archive and in some books, 
but I didn´t find the answer.

I tried to calculate the intervals using:

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

HSDfactor1 <- 
max(abs(tk$"factor(X$A):factor(X$B)"[,2]-tk$"factor(X$A):factor(X$B)"[,3]))
HSDfactor2 <- 
max(abs(tk$"factor(X$A):factor(X$C)"[,2]-tk$"factor(X$A):factor(X$C)"[,3]))
HSDfactor3 <- 
max(abs(tk$"factor(X$B):factor(X$C)"[,2]-tk$"factor(X$B):factor(X$C)"[,3]))


And I modified the function interaction.plot() adding the following 
lines of code:

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

*++ 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

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

|Finally, I call this modified function as:

|interaction.plot2( factor(X$A), factor(X$B), response,las=3, type="b", 
HSDfactor = HSDfactor1, lwd=3)


|However, the resulting intervals are much bigger than the ones 
calculated by STATGRAPHICS, thus I think I did something wrong.

I am not an expert in Statistics nor in R, thus if anyone has any 
suggestion...

Thank you very much,

Manuel.




	[[alternative HTML version deleted]]


