多组数据两两比较相关性

假设我们有多组数据，需要两两比较它们之间的相关性，并使用图形将它们显示出来。

这个问题可以分为几种情况，第一，每组数据数据量都很小，比如只有三个或者四个。第二，每组数据都有一定的数据量，比如十几二十甚至更多，并且成对出现。第三，每组数据都有一定数据量，但没有一一对应的关系，等等。

对于一和三两种情况，我们可以考虑柱状图并连线的方式显示相关性，对于第二种情况，考虑采用pair plot的形式。

首先大约示例一下第一种情况：

> library(gplots) #载入包，因为要使用barplot2函数 > summary(iris) #使用iris数据做为示例数据 Sepal.Length Sepal.Width Petal.Length Petal.Width Species Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100 setosa :50 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300 versicolor:50 Median :5.800 Median :3.000 Median :4.350 Median :1.300 virginica :50 Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800 Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500 > hist.panel = function (x, ...) { # 展示数据的分布，之后我们使用高斯分布，还是其它分布，需要画图确认 + hist(x, + col = "light gray", + probability = TRUE, + ...) + lines(density(x, na.rm=TRUE), + col = "red", + lwd = 1) + rug(x) + } > op <- par(mfrow=c(2,2)) > for(i in 1:4) hist.panel(iris[,i], main=paste("Histogram of", colnames(iris)[i]), xlab=NA) > par(op)

从上图可以看出，有一些是附合高斯分布的，但并不是全部。在这个示例中，我们假设数据都应该符合高斯分布。

# 我们假设其分布符合高斯分布。而后计算95%的置信区间 > avg <- apply(iris[,1:4], 2, mean) # 计算平均值 > s <- apply(iris[,1:4], 2, sd) # 计算标准差 > n <- nrow(iris) # 计数数据量 > avg Sepal.Length Sepal.Width Petal.Length Petal.Width 5.843333 3.057333 3.758000 1.199333 > s Sepal.Length Sepal.Width Petal.Length Petal.Width 0.8280661 0.4358663 1.7652982 0.7622377 > n [1] 150 > error <- qnorm(.95)*s/sqrt(n) # 计算95%置信区间的大小 > error Sepal.Length Sepal.Width Petal.Length Petal.Width 0.1112107 0.0585376 0.2370826 0.1023698 > avg <- matrix(avg, nrow=2, ncol=2, byrow=TRUE) # 对数据进行格式化，以方便使用barplot2画图 > colnames(avg) <- c("Length", "Width") > rownames(avg) <- c("Sepal", "Petal") > avg Length Width Sepal 5.843333 3.057333 Petal 3.758000 1.199333 > error <- matrix(error, nrow=2, ncol=2, byrow=TRUE, dimnames=list(c("Sepal","Petal"), c("Length", "Width"))) > error Length Width Sepal 0.1112107 0.0585376 Petal 0.2370826 0.1023698 > ci.l <- avg - error # 计算置信区间下限 > ci.u <- avg + error # 计算置信区间上限 > ci.l Length Width Sepal 5.732123 2.998796 Petal 3.520917 1.096963 > ci.u Length Width Sepal 5.954544 3.115871 Petal 3.995083 1.301703 > mp <- barplot2(avg, beside = TRUE, # 绘制barplot + col = c("grey12", "grey82"), + legend = rownames(avg), ylim = c(0, 10), + cex.names = 1.2, plot.ci = TRUE, ci.l = ci.l, ci.u = ci.u) > mp # 查看mp，这个barplot2的返回值是每个样品绘制时的x轴位置 [,1] [,2] [1,] 1.5 4.5 [2,] 2.5 5.5 > y.cord<-rbind(c(ci.u[1,]*1.05),c(apply(ci.u,2,max)*1.10), + c(apply(ci.u,2,max)*1.10),c(ci.u[2,]*1.05)) > x.cord<-apply(mp,2,function(x) rep(x,each=2)) # 计算坐标，以方便画线连接要比较的样品 > sapply(1:ncol(x.cord),function(x) lines(x.cord[,x],y.cord[,x])) # 画线连接样品 > x.text<-colMeans(mp) # 计算画*号的x坐标 > y.text<-apply(ci.u,2,max)*1.2 # 计算画*号的y坐标 > panel.cor <- function(x, y, x.cord, y.cord, col=2) # 在相应的位置使用cor.text来计算它的可信度，并用*号表示 + { + test <- cor.test(x,y) + # borrowed from printCoefmat + Signif <- symnum(test$p.value, corr = FALSE, na = FALSE, + cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), + symbols = c("***", "**", "*", ".", " ")) + + text(x.cord, y.cord, Signif, col=col) + } > colnames(iris) [1] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" "Species" > for(i in 1:ncol(avg)) panel.cor(iris[, i], iris[, i+nrow(avg)], x.text[i], y.text[i]) # 打星号

Barplot with significant differences and interactions

对于第二类情况，我们可以使用的包有很多，但是比较常用的基本上只有两个，分别是psych以及PerformanceAnalytics。我们首先看一眼不使用统计分析的pairs绘图的效果。

> pairs(iris[-5], bg=iris$Species, pch=21)

 

Scatterplot Matrices

从图中我们可以看到，pairs把图分成了多个部分，lower.panel, upper.panel, diag.panel以及text.panel。其实最后两个部分都是针对对角线上的小图的，而后就是前面的下三角区域和上三角区域了。两psych以及PerformanceAnalytics就是使用了pairs来绘制相应的图。

 

> library(psych) > pairs.panels(iris[,-5], bg=iris$Species, pch=21)

再例：

> library(PerformanceAnalytics) > chart.Correlation(iris[-5], bg=iris$Species, pch=21)

这些类型的图，其实还有很多包可以画，比如ggplot2等。但这里就不一一介绍了。

除了使用pairs这种类型的图以外，还有一种比较常用的表示相互关系的图–热图（heatmap）。

> image( + z = cor( x = iris[-5], method = "spearman" ), + axes = FALSE, + zlim = c( -1.0, 1.0 ), + col=colorRampPalette(c("green","black","red"))(100)) > axis( + side = 1, + labels = names( iris[-5] ), + at = seq( 0, 1, length = length( names( iris[-5] ) ) ), + cex.axis = 0.8 ) > axis( + side = 2, + labels = names( iris[-5] ), + at = seq( 0, 1, length = length( names( iris[-5] ) ) ), + cex.axis = 0.8 )

但是这种热图表示方式有些失于简单，不过可以一目了然。