当前位置：首页 > news >正文

相关性及P值计算过程

news 2025/9/4 16:42:45

# 载入所需包-------------------------
library(piecewiseSEM) # 用于结构方程模型分析
library(nlme) # 用于线性混合效应模型
library(semPlot) # 绘制结构方程模型的图形
library(lavaan) # 主要用于验证性因子分析(cfa)
library(ggplot2) # 绘图包
library(readxl) # 读取Excel文件
library(devEMF) # 用于绘图的EMF文件输出
library(dplyr) # 数据操作
library(Hmisc) # 相关性计算

# 设置工作目录-------------------------
setwd("E:/study/LW/YEBCD/zuotu2/3权衡协同/MGWR对") # 设置工作目录
getwd() # 输出当前工作目录

# 1. 数据准备及处理-------------------------
theme_set(theme_bw()) # 设置主题样式为白底黑字

# 1.1 数据分类及读取-------------------------
YEBQDYS = read_excel("SEM.xlsx", sheet = "QDL") # 从Excel文件读取数据

# 1.2 将数据按LX分为上、中、下游区域---------------
data_up <- YEBQDYS %>% filter(LX == "1L") # 上游
data_mid <- YEBQDYS %>% filter(LX == "2L") # 中游
data_down <- YEBQDYS %>% filter(LX == "3L") # 下游

# 2. 选择需要计算相关性的列-----------------
selected_columns <- c("Tem", "Pre", "ET", "NDVI", "DEM", "Slope", "Tsand", "Tsilt", "Tclay",
"Toc", "PCL", "PWL", "PGL", "PWS", "PCS", "RD", "GPR", "FIAV", "SIAV",
"TIAV", "NL", "POP") # 选择的列名

# 3. 计算相关性、显著性p值-----------------
# 3.1 计算相关性矩阵-------------------------
correlation_matrix_up <- cor(data_up[, selected_columns], data_up[, c("WY", "SC", "CS")], use = "complete.obs") # 上游相关性
correlation_matrix_mid <- cor(data_mid[, selected_columns], data_mid[, c("WY", "SC", "CS")], use = "complete.obs") # 中游相关性
correlation_matrix_down <- cor(data_down[, selected_columns], data_down[, c("WY", "SC", "CS")], use = "complete.obs") # 下游相关性

# 3.2 计算P值-------------------------
# 将数据转换为数值型
data_up_numeric <- data_up[, selected_columns] %>% mutate(across(everything(), as.numeric))
data_mid_numeric <- data_mid[, selected_columns] %>% mutate(across(everything(), as.numeric))
data_down_numeric <- data_down[, selected_columns] %>% mutate(across(everything(), as.numeric))

# 确保服务列（WY, SC, CS）为数值型
data_up_services <- data_up[, c("WY", "SC", "CS")] %>% mutate(across(everything(), as.numeric))
data_mid_services <- data_mid[, c("WY", "SC", "CS")] %>% mutate(across(everything(), as.numeric))
data_down_services <- data_down[, c("WY", "SC", "CS")] %>% mutate(across(everything(), as.numeric))

# 3.3 计算相关性并获取P值-------------------------
cor_test_up <- rcorr(as.matrix(data_up_numeric), as.matrix(data_up_services)) # 上游
cor_test_mid <- rcorr(as.matrix(data_mid_numeric), as.matrix(data_mid_services)) # 中游
cor_test_down <- rcorr(as.matrix(data_down_numeric), as.matrix(data_down_services)) # 下游

# 输出相关系数和P值-------------------------
cat("上游相关系数:\n")
print(cor_test_up$r)
cat("上游P值:\n")
print(cor_test_up$P)

cat("中游相关系数:\n")
print(cor_test_mid$r)
cat("中游P值:\n")
print(cor_test_mid$P)

cat("下游相关系数:\n")
print(cor_test_down$r)
cat("下游P值:\n")
print(cor_test_down$P)

# 4. 计算描述性统计（均值、标准差、最小值、最大值）-------------------------
summary_stats_up <- data_up[, selected_columns] %>% summarise_all(list(mean = ~mean(.), sd = ~sd(.), min = ~min(.), max = ~max(.)))
summary_stats_mid <- data_mid[, selected_columns] %>% summarise_all(list(mean = ~mean(.), sd = ~sd(.), min = ~min(.), max = ~max(.)))
summary_stats_down <- data_down[, selected_columns] %>% summarise_all(list(mean = ~mean(.), sd = ~sd(.), min = ~min(.), max = ~max(.)))

# 输出描述性统计-------------------------
cat("上游描述性统计:\n")
print(summary_stats_up)

cat("中游描述性统计:\n")
print(summary_stats_mid)

cat("下游描述性统计:\n")
print(summary_stats_down)

# 5. 保存相关性和P值结果-------------------------
# 5.1 保存相关性结果为CSV-------------------------
write.csv(cor_test_up$r, "correlation_up.csv")
write.csv(cor_test_mid$r, "correlation_mid.csv")
write.csv(cor_test_down$r, "correlation_down.csv")

# 5.2 保存P值结果为CSV-------------------------
write.csv(cor_test_up$P, "p_values_up.csv")
write.csv(cor_test_mid$P, "p_values_mid.csv")
write.csv(cor_test_down$P, "p_values_down.csv")

# 6. 提取相关矩阵并重新命名列-------------------------
# 提取相关性矩阵并重命名列
cor_up <- cor_test_up$r[, c("WY", "SC", "CS")]
colnames(cor_up) <- c("CorWYup", "CorSCup", "CorCSup")

cor_mid <- cor_test_mid$r[, c("WY", "SC", "CS")]
colnames(cor_mid) <- c("CorWYmid", "CorSCmid", "CorCSmid")

cor_down <- cor_test_down$r[, c("WY", "SC", "CS")]
colnames(cor_down) <- c("CorWYdown", "CorSCdown", "CorCSdown")

# 提取P值并重命名列
pval_up <- cor_test_up$P[, c("WY", "SC", "CS")]
colnames(pval_up) <- c("PValWYup", "PValSCup", "PValCSup")

pval_mid <- cor_test_mid$P[, c("WY", "SC", "CS")]
colnames(pval_mid) <- c("PValWYmid", "PValSCmid", "PValCSmid")

pval_down <- cor_test_down$P[, c("WY", "SC", "CS")]
colnames(pval_down) <- c("PValWYdown", "PValSCdown", "PValCSdown")

# 6.1 合并相关性矩阵-------------------------
cor_combined <- cbind(cor_up, cor_mid, cor_down) # 合并相关性矩阵
pval_combined <- cbind(pval_up, pval_mid, pval_down) # 合并P值矩阵

# 6.2 获取行名并将它们转换为长列-------------------------
cor_row_names <- rep(rownames(cor_combined), times = ncol(cor_combined)) # 重复行名
pval_row_names <- rep(rownames(pval_combined), times = ncol(pval_combined)) # 重复行名

# 将矩阵转换为长列-------------------------
cor_long <- as.vector(cor_combined) # 相关性长列
pval_long <- as.vector(pval_combined) # P值长列

# 创建最终数据框-------------------------
final_df <- data.frame(RowName = c(cor_row_names, pval_row_names),
Correlation = cor_long,
PValue = pval_long)

# 输出为CSV文件-------------------------
write.csv(final_df, "correlation_and_pvalues_with_row_names.csv", row.names = FALSE)

查看全文

http://www.xdnf.cn/news/37045.html