如何基于rmCorr输出构建相关矩阵并使用p.adjust完成p值校正

你提出的需求完全可以实现，以下是具体实现思路和可直接运行的代码：

实现思路 #

• 先提取需要计算重复测量相关的12个变量的列名，生成所有两两不重复的变量组合
• 遍历所有组合调用rmcorr()计算，自动提取每对变量的相关系数r和原始p值
• 汇总所有p值后调用p.adjust()完成FDR（错误发现率）校正，再分别填充到相关系数矩阵、原始p值矩阵、校正后p值矩阵的下三角，和你要求的结构完全一致

完整实现代码 #

第一步：安装加载所需包 #

# 首次运行先安装包，已安装可跳过这行
install.packages("rmcorr")
# 加载rmcorr包
library(rmcorr)

第二步：导入你的数据 #

注意导入数据时删除原示例结构里的.internal.selfref = <pointer: 0x7f8ba48128e0>部分，否则会报错：

mydata <- structure(list(subjectID = c(1, 1, 1, 1, 1, 1), DSTSpeed = c(5.4225, 
6.8532, 5.6649, 5.6137, 6.5338, 6.9774), DSTError = c(0.060606, 
0.11111, 0.032258, 0.0625, 0.068966, 0.11538), CRTSpeed = c(0.46195, 
0.5066, 0.53191, 0.48758, 0.50286, 0.47727), CRTError = c(0.017241, 
0.034483, 0, 0, 0.033898, 0.016949), KSS = c(4L, 4L, 8L, 8L, 
8L, 6L), SIQPhys = c(1.4, 2, 2.8, 3.6, 3.4, 2.2), SIQCog = c(4, 
3.8, 4.2, 4.2, 5.2, 3.8), TotalSleep = c(7.66416666666667, 7.49611111111111, 
7.28944444444444, 7.78611111111111, 7.46916666666667, 12.8872222222222
), SleepEfficiency = c(0.85775, 0.75881, 0.69097, 0.80629, 0.84559, 
0.73939), ProportionSWS = c(0.063709, 0.31109, 0.2135, 0.2107, 
0.46937, 0.2988), EDA = c(0.77086, 1.4112, 1.5735, 2.168, 1.0156, 
1.7074), WakingEDA = c(0.031424, 0.020836, 0.022987, 0.022799, 
0.020879, 0.28959), temp = c(34.904, 35.414, 35.056, 35.248, 
35.39, 35.105), WakingTemp = c(35.999, 35.636, 35.749, 35.336, 
35.66, NA)), row.names = c(NA, -6L), class = c("data.table", 
"data.frame"))

第三步：核心计算逻辑 #

# 1. 提取要计算相关的变量列名，这里默认取除了subjectID之外的所有列，你有12个变量也可以手动指定：
# var_names <- c("变量1列名", "变量2列名", ..., "变量12列名")
var_names <- setdiff(colnames(mydata), "subjectID")
n_var <- length(var_names)

# 2. 生成所有两两不重复的变量组合
var_pairs <- combn(var_names, 2, simplify = FALSE)

# 3. 遍历所有组合计算rmcorr，提取r和p值
res_list <- lapply(var_pairs, function(pair) {
  fit <- rmcorr(participant = subjectID, var1 = .data[[pair[1]]], var2 = .data[[pair[2]]], dataset = mydata)
  return(data.frame(var1 = pair[1], var2 = pair[2], r = fit$r, p = fit$p))
})
res_df <- do.call(rbind, res_list)

# 4. 对所有p值做FDR校正
res_df$p_adjust <- p.adjust(res_df$p, method = "fdr")

# 5. 构建三个目标矩阵，下三角存储结果
# 初始化矩阵
r_mat <- matrix(NA, nrow = n_var, ncol = n_var, dimnames = list(var_names, var_names))
p_mat <- matrix(NA, nrow = n_var, ncol = n_var, dimnames = list(var_names, var_names))
p_adj_mat <- matrix(NA, nrow = n_var, ncol = n_var, dimnames = list(var_names, var_names))
diag(r_mat) <- 1 # 变量与自身的相关系数为1

# 填充下三角
for (i in 1:nrow(res_df)) {
  row_idx <- which(rownames(r_mat) == res_df$var1[i])
  col_idx <- which(colnames(r_mat) == res_df$var2[i])
  r_mat[col_idx, row_idx] <- res_df$r[i]
  p_mat[col_idx, row_idx] <- res_df$p[i]
  p_adj_mat[col_idx, row_idx] <- res_df$p_adjust[i]
}

结果说明 #

• r_mat：相关系数矩阵，下三角为两两重复测量相关的r值，对角线为1
• p_mat：原始p值矩阵，结构和你要求的完全一致，下三角为每对变量的原始p值
• p_adj_mat：FDR校正后的p值矩阵，下三角为校正后的p值

内容的提问来源于stack exchange，提问作者drbroom