R语言 类别值与连续值可以做相关性分析吗 r语言求相关性

R语言的偏相关分析过程

• 偏相关分析介绍

• 阶偏相关分析：
• 偏相关系数：

• Pearson相关系数
• Spearman相关系数
• Kendall等级相关系数

• R语言实现偏相关分析

• pcor()
• pcor.test()
• 示例数据简介

• R包加载和导入数据
• pcor、pcor.test函数的使用及其结果解读
• 编写函数完成批量计算偏相关系数
• 绘制热图来展示解释变量与响应变量间的偏相关系数和显著性

   当我们想要定量描述两个变量之间的联系时，一般可以使用相关分析来判断二者之间是正相关，负相关，还是不存在相关？再计算二者的相关系数。但在现实中，两变量的相关性往往受其他变量影响，导致两者间的相关程度不能被真实体现。当存在可能会影响两变量之间相关性的因素时，就需要采用偏相关分析，控制其他因素对解释变量和响应变量之间的影响，使得到的结论更加真实可靠。本期主要分享R语言进行偏相关分析时的方法：

1) 使用ppcor包中的pcor.test() 与pcor()函数 计算一个矩阵中等长度两个变量的偏相关性；

2) 编写函数实现：在确定响应变量与需要控制的变量后，计算剩余所有变量与响应变量变量之间的偏相关性，并进行绘图。

偏相关分析介绍

   在多个变量中，排除了一个或若干个控制变量的影响后，直接计算两个变量之间的相关程度。这更有助于准确理解变量间的相互影响、相互作用的关系。通俗点说，如果我们的数据集中有A,B,C,D,E五个变量，其中E为响应变量。当我们分析A,E间的相关性时，可能B的变化会同时引起A,E的变化。消除了B影响后的相关，称为偏相关。
   偏相关分析的主要作用在于，在所有的变量中，判断哪些解释变量对响应变量的影响较大，从而选择作为必需的解释变量，至于那些对响应变量影响较小的自变量就可以舍去而不予考虑。最直接的做法便是，轮流控制某个（一阶）/某类（高阶）变量，看剩余的解释变量与响应变量的相关程度。如果出现：在控制某些变量后，其他变量与响应变量的相关系数骤减，由显著到不显著。而控制其他变量，他们与响应变量的相关性基本不变，就说明这些变量更值得关注。

阶偏相关分析：

   偏相关分析的阶数表示，在对两个变量进行偏相关分析时，被控制的变量数目。当数目为0时，称为零阶偏相关或者（简单）相关，本期主要讲解一阶偏相关分析。

偏相关系数：

   相关系数通常用来表示相关性大小中最常用的指标。pcor系数介于-1~1之间，越接近0相关性越弱，越接近-1 or 1则相关性越强，并存在对应的负正相关。根据分析时使用的方法，可以分为：

Pearson相关系数

适用于两个符合正态分布的连续变量；

Spearman相关系数

适用于不满足Pearson相关系数正态分布要求的连续变量或有序变量间的相关性比较；

Kendall等级相关系数

适用于两个有序变量间的相关性比较

R语言实现偏相关分析

本文中所有的分析均为“Spearman”相关系数。主要利用 pcor包中的pcor()、 pcor.test() 函数实现偏相关分析。以环境因子调控土壤微生物残体积累”的案例进行讲解。例子讲解前，我们先用?pcor()、?pcor.test()了解一下函数的用法和命令。

函数的基本命令：

pcor(x= data, method= )

pcor.test(x=data$"解释变量", y=data$"响应变量"， z=data$“控制变量”, method= )

对应参数的描述：

pcor()

值得注意的是：pcor()会将两个变量视作一组，然后自动控制其他变量计算二者的偏相关系数，所以结果的格式与rcorr()类似，有兴趣的同学可以对比一下。相较之下，pcor.test()函数在控制变量的选择上更加自由，因此在进行偏相关分析时，更倾向使用pcor.test()函数。

• 
  
• 
  

pcor.test()

• 
  
• 
  
• 
  

示例数据简介

数据来自一项关于环境因子调控土壤微生物残体（Amino sugars）积累的研究数据，记录了38个观测点的各项指标：生态系统（Ecosystem）、干旱指数（Aridity index）、植物地上生物量（AGB）、植物地下生物量（BGB）、土壤有机碳含量（SOC）、土壤氮含量（N）、土壤pH、土壤粘粒组分（Clay）、磷脂脂肪酸含量（PLFA）、氨基糖含量（Amino sugars）信息。接下来，通过分析，找到对残体积累影响最明显的因子。

注：本数据在原始数据的基础上经过调整与删改，分析的结果无任何参考价值。

R包加载和导入数据

加载R包：

getwd() #若当前路径不是文件所在位置，则运行第二条代码
# ctrl + shift + h 快速设定工作路径
rm(list = ls()) # 清除当前工作环境
library(tidyverse)
library(ggplot2)
library(Hmisc)
#install.packages("ppcor")
library(ppcor)

数据导入：

test_data<- openxlsx::read.xlsx("test_data.xlsx", 1)[1:38, -c(1:2)]
head(test_data)

数据格式：

pcor、pcor.test函数的使用及其结果解读

首先，我们先使用Hmisc中的rcorr()函数计算干旱指数与残体浓度的简单相关系数。结果显示：二者之间的相关系数为0.31，p值为0.056。

cor_uncontrolled<- rcorr(as.matrix(test_data), type= "spearman")$r %>%
  subset(, select= Amino.sugars) %>% data.frame() # 相关系数
p_uncontrolled<- rcorr(as.matrix(test_data), type= "spearman")$P %>%
  subset(, select= Amino.sugars) %>% data.frame() # p值
uncontrolled<- cbind(cor_uncontrolled, p_uncontrolled) # 合并数据
uncontrolled

然后我们使用ppcor包中的pcor.test()函数，在控制其他变量的情况下，计算干旱指数与残体浓度的一阶与二阶偏相关系数。

由于函数的特殊性，在使用该函数进行高阶偏相关分析时，需要注意对控制变量的调用方式。

一般建议使用： z = data[, c(“var1-name”, “var2-name”,…)] 或data[, c(1:2,5)] 这两种方法调用需控制的变量。

#控制AGB时的一阶偏相关系数
pcor.test(x= test_data$Aridity.index, #相关变量
       y= test_data$Amino.sugars, #目的变量
       z= test_data$AGB, method= "spearman") #控制的变量

#控制BGB时的二阶偏相关系数
pcor.test(x= test_data$Aridity.index,
       y= test_data$Amino.sugars,
       z= test_data[,c("AGB", "BGB")],
       method= "spearman")

结果如下：

从结果可知，estimate正是在控制AGB（BGB）后二者的一阶（二阶）偏相关系数，与简单相关相比，相关系数下降了0.22，表明植物地上部生物量对二者的相关性存在很强的影响。同理，还可以使用pcor()函数完成二阶或者高阶的偏相关分析，但由于pcor()的特殊性，为了保证结果一致，需要新建一个只有干旱指数、AGB、BGB、Amino sugars的数据框。

newdata<- test_data[,c("Aridity.index", "AGB", "BGB", "Amino.sugars")]
pcor(newdata, method = "spearman")

以上结果显示，在统计学范畴内，部分变量（如AGB等）的存在对干旱指数与残体浓度的相关性存在很大的影响。这表明，如果不考虑其他变量的影响，仅通过干旱指数与残体浓度数据计算得到的相关系数是不可靠的。在数据分析时，我们要结合变量在生态学中的联系，合理控制可能对结果有影响的其他变量。

编写函数完成批量计算偏相关系数

在使用pcor与pcor.test计算偏相关系数时，我们会发现：在确定解释变量、响应变量与需要控制的变量后，每次计算只能得到一组变量间的偏相关系数，希望批量完成；而pcor虽然能得到数据集中所有“变量对”的偏相关系数，但不能自由的选择需要控制的解释变量。因此，我们编写了一个函数用于批量计算偏相关系数。

partialCor<- function(data = test_data, x = NULL,
                      y, method= "spearman"){
  require(ppcor)
  require(tidyverse)
  
  dat<- as.data.frame(data[, -which(colnames(data) %in% c(x , y))])
  if (ncol(dat) <= 1) {
    warning("Small number of remained variable may causes an invalid var1 name")
    colnames(dat) = colnames(data)[-match(c(x,y),colnames(data))]
  }
  dat1<- data.frame()
  
  if(is.null(x) == T){
    for(i in 1:ncol(dat)){
      for(j in 1:ncol(dat)){
        if(i == j) next
        dat2 <- pcor.test(dat[[j]], data[,y], dat[[i]], method = method) %>%
          as_tibble() %>%
          mutate(var2 = y,
                 var1 = names(dat[j]),
                 control_var = names(dat[i]),
                 signif = ifelse(p.value< 0.001, "***",
                                 ifelse(p.value< 0.01, "**",
                                        ifelse(p.value< 0.05, "*", "")))) %>%
          dplyr::select(var1,var2,
                        control_var,estimate,
                        p.value,signif,everything())
        
        dat1<- rbind(dat1,dat2)
      }
    }
  }
  else{
    for(i in 1:ncol(dat)){
      dat2 <- pcor.test(dat[[i]], data[,y], data[,x], method = method) %>%
        as_tibble() %>%
        mutate(var2 = y,
               var1 = names(dat[i]),
               control_var = str_c(x, collapse = "-"),
               signif = ifelse(p.value< 0.001, "***",
                               ifelse(p.value< 0.01, "**",
                                      ifelse(p.value< 0.05, "*", "")))) %>%
        dplyr::select(var1,var2,
                      control_var,estimate,
                      p.value,signif,everything())
      
      dat1<- rbind(dat1,dat2)
    }
  }
  dat1
}

注：以上代码粘贴复制到R或者Rstudio后，直接在最后一个大括号“}”后点击运行即可

参数介绍：

• 
  
• 
  

1）不提供x参数时，输出结果如下：

result<- partialCor(data = test_data, y="Amino.sugars", method = "spearman")
result

2）提供x参数时，输出的结果如下：

partialCor(data = test_data, x= c("SOC", "N", "pH"), y="Amino.sugars", method = "spearman")

其中，var1是解释变量，var2是响应变量，control_var是控制的变量，estimate代表偏相关系数。在使用该函数进行分析时，可以根据这些统计结果，结合实际数据，解释可能存在的生态学意义即可。

绘制热图来展示解释变量与响应变量间的偏相关系数和显著性

使用ggplot2包中的geom_tile函数对环境变量与残体浓度偏相关分析的结果进行绘图：

result<- partialCor(data = test_data,
                    y = "Amino.sugars",
                    method= "spearman")
result1<- subset(result, select = c(var1,var2, control_var, estimate, signif))
result1$estimate<- round(result1$estimate, digits = 2)
result1$var1<- factor(result1$var1, levels = c("PLFA", "pH", "Clay", "N", "SOC",
                                               "BGB", "AGB", "Aridity.index"))
result1$control_var<- factor(result1$control_var,
                             levels = c("Aridity.index", "AGB", "BGB", "SOC",
                                        "N", "Clay", "pH", "PLFA"))

ggplot(data=result1) +
  geom_tile(aes(control_var, var1, fill= estimate))+
  geom_text(aes(control_var, var1, label= estimate),vjust=0, size=2)+
  geom_text(aes(control_var, var1, label= signif),vjust=1, size=3)+
  scale_fill_gradientn(limit=c(-0.5,0.5),
                       colors=c('#4987D3', "#FBF7C7",'#DB562D'))+
  facet_grid(~var2)+
  theme(panel.grid = element_blank(),
        strip.text = element_text(colour = 'black', size = rel(1)),
        strip.background = element_rect(fill = "grey", colour = 'grey', size = rel(1)),
        panel.background = element_rect(fill = "#f2f2f2", color = "#f2f2f2"),
        axis.text = element_text(vjust=1, size = 6),
        axis.title = element_text(vjust = 1, size = 8),
        legend.text = element_text(size = 4),
        legend.title = element_text(size = 6),
        legend.key.size = unit(10, "pt"))+
  labs(x = 'Controlled variables', y = '', fill = "Correlations\nCoefficient r", size = 3)+
  scale_y_discrete(expand = c(0, 0))+
  scale_x_discrete(expand = c(0, 0))

export::graph2ppt(file="partialheatmap.ppt", width=6 ,heigh=3, append= T)