可视化 R 中事后 Tukey 的临界值/成对比较答案

【问题标题】：Visualize critical values / pairwise comparisons from posthoc Tukey in R可视化 R 中事后 Tukey 的临界值/成对比较
【发布时间】：2011-12-05 20:57:44
【问题描述】：

我正在尝试对我从 posthoc Tukey 获得的关键值进行细粒度的可视化。有一些good guidelines out there 用于可视化成对比较，但我需要更精致的东西。我的想法是，我会有一个图，其中每个小方块都代表下面矩阵中的一个临界值，编码方式如下：

如果该值大于或等于 5.45 - 它是一个黑色方块；
如果该值小于或等于-5.45 - 它是一个灰色方块；
如果值在 -5.65 和 5.65 之间 - 它是一个白色方块。

数据矩阵为here。

或者您可能有更好的建议如何可视化这些关键值？

编辑：继 @Aaron 和 @DWin 的 cmets 之后，我想为上述数据提供更多背景信息，并为我的问题提供理由。我正在查看七个虚拟角色的可接受性平均评分，每个角色都有 5 个不同级别的动画。所以，我有两个因素——角色（7 个级别）和动作（5 个级别）。因为我发现了这两个因素之间的相互作用，所以我决定查看所有角色对于所有运动级别的平均值之间的差异，这导致了这个巨大的矩阵，作为 posthoc Tukey 的输出。现在可能细节太多了，但请不要把我扔到交叉验证，他们会活活吃掉我......

【问题讨论】：

更好的建议取决于您想用这些数据讲述什么故事。通常按照平均值对级别进行排序，但看起来您可能在级别中有一些具有特殊含义的模式；特别是，似乎每五个可能以某种特定方式相关（或不相关）。我也对可视化临界值持谨慎态度，因为这是衡量统计显着性的指标，因此会隐藏结果的实际意义（或缺乏实际意义）。
来吧，你可以做到的！发布您的第一个交叉验证问题...我敢于挑战您！（真的，我们在那里并没有那么可怕。）只要解释你的数据，你想知道什么，最重要的是，你为什么想知道它。你不必有答案；这就是你问这个问题的原因。

标签： r matrix data-visualization

【解决方案1】：

image 相当简单：

d <- as.matrix(read.table("http://dl.dropbox.com/u/2505196/postH.dat"))    
image(x=1:35, y=1:35, as.matrix(d), breaks=c(min(d), -5.45, 5.45, max(d)), 
      col=c("grey", "white", "black"))

对于只有一半，使用d[upper.tri(d)] <- NA 将一半设置为缺失并将na.rm=TRUE 添加到 min 和 max 函数。

【讨论】：

【解决方案2】：

这是一个 ggplot2 解决方案。我敢肯定有更简单的方法可以做到这一点——我想我被冲昏了头脑！

library(ggplot2)

# Load data.
postH = read.table("~/Downloads/postH.dat")
names(postH) = paste("item", 1:35, sep="") # add column names.
postH$item_id_x = paste("item", 1:35, sep="") # add id column.

# Convert data.frame to long form.
data_long = melt(postH, id.var="item_id_x", variable_name="item_id_y")

# Convert to factor, controlling the order of the factor levels.
data_long$item_id_y = factor(as.character(data_long$item_id_y), 
                        levels=paste("item", 1:35, sep=""))
data_long$item_id_x = factor(as.character(data_long$item_id_x), 
                        levels=paste("item", 1:35, sep=""))

# Create critical value labels in a new column.
data_long$critical_level = ifelse(data_long$value >= 5.45, "high",
                             ifelse(data_long$value <= -5.65, "low", "middle"))

# Convert to labels to factor, controlling the order of the factor levels.
data_long$critical_level = factor(data_long$critical_level,
                                  levels=c("high", "middle", "low"))

# Named vector for ggplot's scale_fill_manual
critical_level_colors = c(high="black", middle="grey80", low="white")

# Calculate grid line positions manually.
x_grid_lines = seq(0.5, length(levels(data_long$item_id_x)), 1)
y_grid_lines = seq(0.5, length(levels(data_long$item_id_y)), 1)

# Create plot.
plot_1 = ggplot(data_long, aes(xmin=as.integer(item_id_x) - 0.5,
                               xmax=as.integer(item_id_x) + 0.5,
                               ymin=as.integer(item_id_y) - 0.5,
                               ymax=as.integer(item_id_y) + 0.5,
                               fill=critical_level)) +
     theme_bw() +
     opts(panel.grid.minor=theme_blank(), panel.grid.major=theme_blank()) +
     coord_cartesian(xlim=c(min(x_grid_lines), max(x_grid_lines)),
                     ylim=c(min(y_grid_lines), max(y_grid_lines))) +
     scale_x_continuous(breaks=seq(1, length(levels(data_long$item_id_x))),
                        labels=levels(data_long$item_id_x)) +
     scale_y_continuous(breaks=seq(1, length(levels(data_long$item_id_x))),
                        labels=levels(data_long$item_id_y)) +
     scale_fill_manual(name="Critical Values", values=critical_level_colors) +
     geom_rect() +
     geom_hline(yintercept=y_grid_lines, colour="grey40", size=0.15) +
     geom_vline(xintercept=x_grid_lines, colour="grey40", size=0.15) +
     opts(axis.text.y=theme_text(size=9)) +
     opts(axis.text.x=theme_text(size=9, angle=90)) +
     opts(title="Critical Values Matrix")

# Save to pdf file.
pdf("plot_1.pdf", height=8.5, width=8.5)
print(plot_1)
dev.off()

【讨论】：

哇，这真的是一个巨大的解决方案@bdemarest，但非常感谢你的努力。代码描述也很棒。一个问题 - 你如何只显示这个矩阵的一半，对角线切割它？
@Geek On Acid：我已经考虑了一段时间，但我没有看到一个简单的方法来做到这一点。看来我的解决方案不是很通用/可定制！也许其他人对此有想法......
你不能改用ggfluctuation()吗？

【解决方案3】：

如果您将 findInterval 设置为bg、col 和/或pch 参数的索引（尽管目前它们都是正方形），您应该会发现代码相当紧凑且易于理解.

您需要先获取长格式数据；这是一种方法：

d <- as.matrix(read.table("http://dl.dropbox.com/u/2505196/postH.dat"))
dat <- within(as.data.frame(as.table(d)), 
              { Var1 <- as.numeric(Var1)  
                Var2 <- as.numeric(Var2) })

那么代码如下； pch=22 使用填充正方形，bg 设置正方形的填充颜色，col 设置边框颜色，cex=1.5 只是让它们比默认值大一点。

plot(dat$Var1, dat$Var2, 
     bg = c("grey", "white", "black")[1+findInterval(dat$Freq, c(-5.45,5.45))],
     col="white", cex=1.5, pch = 22)

你需要1+，因为值是 0,1,2，你的索引需要从 1 开始。

【讨论】：

介意我是否编辑添加结果图片（代码稍作修改？）

【解决方案4】：

为了结束这里，我使用了@DWin 和@Aaron 的大部分建议来创建下面的图。最浅的灰度代表不重要的值。我还使用rect 在轴名称上方创建线条以更好地区分条件：

d <- as.matrix(read.table("http://dl.dropbox.com/u/2505196/postH.dat"))
#remove upper half of the values (as they are mirrored values)
d[upper.tri(d)] <- NA
dat <- within(as.data.frame(as.table(d)),{
Var1 <- as.numeric(Var1)
Var2 <- as.numeric(Var2)})
par(mar=c(6,3,3,6))
colPh=c("gray50","gray90","black")
plot(dat$Var1,dat$Var2,bg = colPh[1+findInterval(dat$Freq, c(-5.45,5.45))],
    col="white",cex=1.2,pch = 21,axes=F,xlab="",ylab="")
labDis <- rep(c("A","B","C","D","E"),times=7)
labChar <- c(1:7)
axis(1,at=1:35,labels=labDis,cex.axis=0.5,tick=F,line=-1.4)
axis(1,at=seq(3,33,5),labels=labChar, tick=F)
#drawing lines above axis for better identification
rect(1,0,5,0,angle=90);rect(6,0,10,0,angle=90);rect(11,0,15,0,angle=90);
rect(16,0,20,0,angle=90);rect(21,0,25,0,angle=90);rect(26,0,30,0,angle=90);
rect(31,0,35,0,angle=90)
axis(4,at=1:35,labels=labDis,cex.axis=0.5,tick=F,line=-1.4)
axis(4,at=seq(3,33,5),labels=labChar,tick=F)
#drawing lines above axis for better identification
rect(36,1,36,5,angle=90);rect(36,6,36,10,angle=90);rect(36,11,36,15,angle=90);
rect(36,16,36,20,angle=90);rect(36,21,36,25,angle=90);rect(36,26,36,30,angle=90);
rect(36,31,36,35,angle=90)
legend("topleft",legend=c("not significant","p<0.01","p<0.05"),pch=16,
col=c("gray90","gray50","black"),cex=0.7,bty="n")

【讨论】：