“用R进行相关分析”的版本间的差异

2022年11月14日 (一) 12:35的最新版本

脚本与注释

datafilename <- "http://personality-project.org/r/datasets/maps.mixx.epi.bfi.data"   # 指定文件名与路径
mydata <-  read.table(datafilename,header=TRUE)                                      # 读取数据到mydata
attach(mydata)                                                                       # 激活mydata
library(psych)                                                                       # 启用psych包
epi <- data.frame(epiE ,epiS ,epiImp ,epilie ,epiNeur)                               # 从mydata中提取数据子集，只包括要分析的变量
corr.test(epi)                                                                       # 相关分析

结果

> corr.test(epi)
Call:corr.test(x = epi)
Correlation matrix    # 相关分析的结果，可以复制到WORD，然后文字转换成表格
         epiE  epiS epiImp epilie epiNeur
epiE     1.00  0.85   0.80  -0.22   -0.18
epiS     0.85  1.00   0.43  -0.05   -0.22
epiImp   0.80  0.43   1.00  -0.24   -0.07
epilie  -0.22 -0.05  -0.24   1.00   -0.25
epiNeur -0.18 -0.22  -0.07  -0.25    1.00
Sample Size 
[1] 231
Probability values (Entries above the diagonal are adjusted for multiple tests.) 
        epiE epiS epiImp epilie epiNeur
epiE    0.00 0.00   0.00   0.00    0.02
epiS    0.00 0.00   0.00   0.53    0.00
epiImp  0.00 0.00   0.00   0.00    0.53
epilie  0.00 0.43   0.00   0.00    0.00
epiNeur 0.01 0.00   0.26   0.00    0.00

 To see confidence intervals of the correlations, print with the short=FALSE option

结果整理

用R不仅可以完成统计分析，而且可以将统计分析的结果以自己想用的方式输出。在OBHRM研究中，描述性统计一般需要报告：平均值、标准差和相关系数，相关系数还需要标上显著性。下面我们来看，如何汇总得到这些结果。

cor.result<-corr.test(epi)                       # 将相关分析的结果保存为cor.result

# 此处需要根据你的需要修改
digits<-3                                        # 设定所有数据显示为3位小数

# 此处需要修改epi为你自己的数据框
m_sd<-format(round(psych::describe(epi)[3:4],digits),,nsmall=digits)  # 计算平均值和标准差。其中的3:4是只保留平均值和标准差两列


# 后面的代码基本不用修改

m_sd$names<-rownames(m_sd)                      # 新增一列，命名为names，从rownames提取names的值
m_sd<-m_sd[,c(3,1,2)]                           # 调整顺序，将names放在第1列


tempR<-cor.result$r                              # 从cor.result提取相关系数
tempR<-format(round(tempR,digits),nsmall=digits) # 将所有相关系数调整为3位小数
tempR[upper.tri(tempR, diag =T)]<-""             # 将相关系数矩阵上半角和对角线赋值为空白，即""

library(gtools)                                  # 启用gtools
tempP<-cor.result$p                              # 从cor.result提取p值
tempP<-gtools::stars.pval(tempP)                 # 根据p值，转换为相应的星号（*），为简单起见，直接采用stars.pval函数。
tempP[upper.tri(tempP, diag =T)]<-""             # 将p值矩阵上半角和对角线赋值为空白，即""

result<-as.data.frame(matrix(paste(tempR,tempP,sep=""),nrow=nrow(tempR), dimnames=dimnames(tempR)))   # 合并相关系数矩阵和p值矩阵（已根据p值调整为星号）

result<-cbind(m_sd,result)                      # 合并平均值、标准差和相关系数矩阵（已经根据p值标上星号），得到result，可以输出到csv或excel文件

也可以直接用下面的函数，来完成结果的计算和整理。

oCor<-function(data,digits=2,p10=T)
{
	temp<-psych::corr.test(data)

	m_sd.result<-format(round(psych::describe(data)[c(3:4)],digits),nsmall=digits)
	m_sd.result$names<-rownames(m_sd.result)
	m_sd.result<-m_sd.result[,c(3,1,2)]

	tempR<-temp$r 
	tempR<-format(round(tempR,digits),nsmall=digits)
	tempR[upper.tri(tempR, diag =T)]<-""

	tempP<-temp$p
        if(p10==T)
         {
           tempP<-ifelse(tempP < .001, "***", ifelse(tempP < .01, "** ", ifelse(tempP < .05, "*  ", ifelse(tempP < .1, "+   ", "    "))))
          }else
         {
           tempP<-ifelse(tempP < .001, "***", ifelse(tempP < .01, "** ", ifelse(tempP < .05, "*  ", "    ")))
          }      

	tempP[upper.tri(tempP, diag =T)]<-""

	result<-as.data.frame(matrix(paste(tempR,tempP,sep=""),nrow=nrow(tempR), dimnames=dimnames(tempR))) 
	result<-cbind(m_sd.result,result)
	return(result)
}

a<-oCor(epi)