使用R图形化显示不同组数据之间的相互关系 Graphically analyzing variable interactions in R

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

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

对于一和三两种情况,我们可以考虑柱状图并连线的方式显示相关性,对于第二种情况,考虑采用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)

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

# 我们假设其分布符合高斯分布。而后计算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
Barplot with significant differences and interactions

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

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

Scatterplot Matrices
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)

psych
再例:

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

PerformanceAnalytics

这些类型的图,其实还有很多包可以画,比如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 )

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

2 thoughts on “使用R图形化显示不同组数据之间的相互关系 Graphically analyzing variable interactions in R

  1. 博主老师好!请问用R如何作2个差异基因(2种组织,3个时间点)相对表达的线性回归图?谢谢

    1. 这是个统计学问题。我想你的问题应该是如何做线性回归。你可以试试http://data.princeton.edu/R/linearModels.html

boya888进行回复 取消回复

电子邮件地址不会被公开。 必填项已用*标注