在数据探索的过程中,我们往往会对数据与数据的相关性进行分析,例如我们常用的corrplot
包,或者psych
包中的corr.test
函数,对两两变量间的相关性进行分析。我们常常会看到这样的相关性热图:
但有时变量间的关系并非线性,例如:我们熟知的微生物物种多样性与土壤pH之间的关系,因此在数据探索过程中,仅使用线性相关会导致我们遗漏一些重要信息。二次关系也需要纳入我们考虑的范围,因此本篇博客将对两两变量的二次关系的R代码及可视化进行介绍,助力数据探索。
我们尝试在对数据进行线性相关分析的同时,加入二次关系的分析,数据的可视化显示的结果大致如下:其中菱形表示线性相关显著,圆圈表示线性相关不显著,而二次关系显著。
可视化代码遵循如下原则:
- 若线性相关性显著,无论二次关系是否显著,均使用线性相关的R值和P值;(使用最简单的模型即可)
- 若线性相关性不显著,则查看二次关系是否显著,若显著则显示二次关系的R值和P值;若不显著则均不显示。
1.计算变量的两两二次关系函数
qcorr <- function(x){
pm <- matrix(rep(0,ncol(x)^2),nrow = ncol(x))
rm <- matrix(rep(0,ncol(x)^2),nrow = ncol(x))
rownames(pm) <- colnames(x)
colnames(pm) <- colnames(x)
rownames(rm) <- colnames(x)
colnames(rm) <- colnames(x)
for(i in 1:length(x)){
for(j in 1:length(x)){
if (i != j){
x1 <- x[,i]
x2 <- x[,j]
s <- summary(lm(x1~x2+I(x2^2)))
pm[i,j] <- pf(s$fstatistic[1],s$fstatistic[2],s$fstatistic[3],lower.tail=F)
rm[i,j] <- sqrt(s$r.squared)
}
}
}
return(list(p=pm,r=rm))
}
2.与线性相关性进行整合
env_corr2 <- qcorr(env_div)$r
env_corp2 <- qcorr(env_div)$p
env_corr1 <- corr.test(env_div,adjust = "none")$r
env_corp1 <- corr.test(env_div,adjust = "none")$p
corr2 <- melt(env_corr2, na.rm = TRUE)
corp2 <- melt(env_corp2, na.rm = TRUE)
corr1 <- melt(env_corr1, na.rm = TRUE)
corp1 <- melt(env_corp1, na.rm = TRUE)
type <- ifelse(corp1$value>0.05 & corp2$value < 0.05,"Q","L")
corp1[corp1$value>0.05 & corp2$value < 0.05,] <- corp2[corp1$value >0.05 & corp2$value < 0.05,]
corr1[corp1$value>0.05 & corp2$value < 0.05,] <- corr2[corp1$value >0.05 & corp2$value < 0.05,]
corr <- corr1
corr$p_value <- corp1$value
corr$type <- type
corr$abs_corr <- abs(corr$value)*10
corr <- corr[corr$p_value<0.05 & corr$value < 1,]
3.对整合的结果进行可视化
outline.color = "gray"
colors = c("blue", "white", "red")
shape <- cols <- c("Q" = 21, "L" = 23)
library(ggplot2)
p <- ggplot(data = corr, mapping = aes_string(x = "Var1", y = "Var2"))
p <- p + geom_point(color = outline.color,aes_string(size = "abs_corr", fill = "value",shape = "type")) +
scale_size(range = c(4,10)) +
scale_shape_manual(values = shape) +
guides(size = "none")
p <- p + scale_fill_gradient2(low = colors[1], high = colors[3],
mid = colors[2], midpoint = 0, limit = c(-1, 1), space = "Lab")
p <- p + ggplot2::theme_minimal()
p
这里使用的是ggplot2的包,方便对可视化的结果进一步自定义。
输出结果如下:
测试数据及代码可参考如下连接:https://mbd.pub/o/bread/ZZ6WmZhr