在R中计算Spearman系数时如何处理个体内重复测量聚类问题？

你的数据存在个体内重复测量的聚类结构，直接用cor.test会忽略同一受试者多次测量的相关性，导致统计推断偏倚。下面提供两种实用的解决方法：

方法1：秩变换+混合效应模型 #

Spearman相关本质是变量秩次的Pearson相关，因此可以先对测量值做秩变换，再用混合效应模型控制个体聚类效应，以此得到校正后的关联结果。

代码实现 #

# 安装并加载混合效应模型包（首次运行需执行安装）
install.packages("lme4")
library(lme4)

# 生成你的模拟数据
id <- rep(c(1:100), times=3)
df1 <- data.frame(id)
df1$var1 <- sample(500:1000, length(df1$id))
df1$var2 <- sample(500:1000, length(df1$id))

# 对两个测量变量做秩变换
df1$rank_var1 <- rank(df1$var1)
df1$rank_var2 <- rank(df1$var2)

# 拟合带个体随机截距的混合效应模型
model <- lmer(rank_var2 ~ rank_var1 + (1 | id), data = df1)

# 查看结果
summary(model)

结果解释 #

模型中rank_var1的系数就是考虑个体聚类后的Spearman相关系数估计值，模型输出的标准误和p值已经修正了个体内重复测量的相关性。

方法2：聚类Bootstrap法 #

通过按个体而非单个观测值进行抽样，重复计算Spearman相关系数，以此获得稳健的置信区间和p值，避免聚类带来的偏倚。

代码实现 #

# 安装并加载bootstrap工具包（首次运行需执行安装）
install.packages("boot")
library(boot)

# 生成你的模拟数据
id <- rep(c(1:100), times=3)
df1 <- data.frame(id)
df1$var1 <- sample(500:1000, length(df1$id))
df1$var2 <- sample(500:1000, length(df1$id))

# 定义bootstrap计算函数：按个体抽样
spearman_clustered <- function(data, indices) {
  # 抽取指定个体的所有观测
  sampled_ids <- unique(data$id)[indices]
  sampled_data <- data[data$id %in% sampled_ids, ]
  # 计算Spearman相关系数
  cor(sampled_data$var1, sampled_data$var2, method = "spearman")
}

# 执行1000次聚类bootstrap抽样
boot_result <- boot(data = df1, statistic = spearman_clustered, R = 1000, strata = df1$id)

# 查看结果及稳健置信区间
print(boot_result)
boot.ci(boot_result, type = "bca")

结果解释 #

这种方法保留了个体内的聚类结构，得到的置信区间和p值更贴合数据的真实相关性特征，适合验证结果的稳健性。

额外提示 #

如果数据中存在时间趋势，可在模型中加入时间变量作为协变量，进一步控制混杂因素。

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