Manuel López-Ibáñez
2004-Jun-08 15:41 UTC
[R] interaction plot with intervals based on TukeyHSD
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]]