非配对数据关联分析方法探讨——含土壤样本实例

问题背景 #

需要结合采样点（location）开展不同土壤类型间的变量关联分析，但各土壤类型在同一采样点的样本数量不一致，属于非配对数据。尝试两种方法均存在问题：通过inner_join强制配对会导致观测值重复，疑似引入偏差；按采样点取均值实现配对则会丢失样本变异信息，且两种方法的关联结果差异显著。

示例数据 #

df <- structure(list(location = c("A", "A", "A", "A", "B", "B", "B","B", "C", "C", "C", "C", "C", "C", "D", "D"), 
                     type_of_soil = c("type1", "type1", "type2", "type3", "type1", "type2", "type2", "type3", "type1", "type2", "type2", "type2", "type3", "type3", "type1", "type2"), 
                     var1 = c(29, 7, 11, 20, 12, 15, 22, 27, 29, 32, 31, 23, 22, 20, 20, 15), 
                     var2 = c(56, 92, 111, 15, 230, 23, 47, 50, 190, 10, 6, 2, 15, 73, 223, 50), 
                     var3 = c(1000, 1643, 1242, 6352, 111, 3552, 837, 5612, 999, 572, 127, 593, 3232, 524, 41, 2)), 
                class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, -16L))

两种尝试方法的问题分析 #

方法1：宽格式转换+inner_join强制配对 #

用户先将数据转为宽格式，再编写函数通过inner_join按采样点配对：

df$ID <- 1:nrow(df) # 为每个原始样本分配ID
df_longer <- pivot_longer(data = df, cols = -c(1,2,6) , names_to = "var", values_to = "vals")
df_cor <- unite(df_longer, var, type_of_soil, sep = '_', col = 'var_type')
df_wider <- pivot_wider(df_cor, names_from = var_type, values_from = vals)

bicorplot <- function(df_forcor, soil1, soil2){
  coord1 <- grep(soil1, colnames(df_forcor))
  coord2 <- grep(soil2, colnames(df_forcor))
  
  df1 <- na.omit(df_forcor[,c(1, coord1)])
  df2 <- na.omit(df_forcor[,c(1, coord2)])
  
  coord_eq1 <- df1[df1$location %in% df2$location,]
  coord_eq2 <- df2[df2$location %in% df1$location,]
  
  temp_df <- inner_join(coord_eq1, coord_eq2, by = "location")
  temp_df1 <- temp_df[1:4]
  temp_df2 <- temp_df[c(1,5:7)]
  
  cor_mat <- cor(x = temp_df1[2:4], y = temp_df2[2:4], method = "spearman")
  
  corrplot(cor_mat)
}

核心问题：当同一采样点下某土壤类型存在多个样本时，inner_join会生成笛卡尔积式的重复观测。例如某采样点有2个type1样本和2个type3样本，配对后会生成4行数据，该采样点的权重被放大4倍，直接扭曲了关联分析的结果，引入严重偏差。

方法2：按采样点取均值后分析 #

用户先按采样点合并均值，再进行关联分析：

df_mean <- df_wider %>%
  group_by(location) %>%
  summarise_all(mean, na.rm = T)

bicorplot(df_wider, 'type1', 'type3')
bicorplot(df_mean, 'type1', 'type3')

核心问题：取均值会完全丢失同一采样点内的样本变异信息，而这些变异往往包含变量关联的关键细节。同时，若不同采样点的样本数量差异较大，均值会被样本多的采样点过度主导，导致关联结果过于平滑，无法反映真实的变量关系。

正确分析思路：考虑采样点聚类的混合效应/稳健关联分析 #

由于数据存在嵌套结构（采样点下包含多个不同土壤类型的样本），需控制采样点带来的组内非独立性，避免偏差。推荐两种严谨的方法：

方法1：混合效应模型（推荐） #

使用混合效应模型将采样点作为随机效应，控制组内变异，直接在原始数据上分析变量关联：

library(lme4)
library(dplyr)
library(tidyr)

# 筛选同时包含两种目标土壤类型的采样点
valid_locs <- df %>%
  group_by(location) %>%
  filter(all(c("type1", "type3") %in% type_of_soil)) %>%
  pull(location) %>%
  unique()

# 提取目标数据并转为配对宽格式
df_paired <- df %>%
  filter(location %in% valid_locs, type_of_soil %in% c("type1", "type3")) %>%
  pivot_wider(
    names_from = type_of_soil,
    values_from = c(var1, var2, var3),
    names_sep = "_"
  )

# 构建混合效应模型，以location为随机截距
model_var1 <- lmer(var1_type1 ~ var1_type3 + (1|location), data = df_paired)
summary(model_var1)

# 批量分析所有变量对并可视化
cor_mat <- sapply(c("var1", "var2", "var3"), function(v1) {
  sapply(c("var1", "var2", "var3"), function(v2) {
    formula <- as.formula(paste0(v1, "_type1 ~ ", v2, "_type3 + (1|location)"))
    coef(summary(lmer(formula, data = df_paired)))[2,1]
  })
})

corrplot::corrplot(cor_mat, method = "color", title = "混合效应模型下的变量关联系数")

方法2：稳健聚类关联分析 #

使用psych包的corr.test函数，指定采样点为聚类变量，控制非独立性：

library(psych)
library(corrplot)

# 生成所有土壤类型配对的观测数据
df_pair <- inner_join(
  df %>% filter(type_of_soil == "type1") %>% rename_with(~paste0(., "_type1"), starts_with("var")),
  df %>% filter(type_of_soil == "type3") %>% rename_with(~paste0(., "_type3"), starts_with("var")),
  by = "location"
)

# 执行聚类稳健的关联分析
corr_result <- corr.test(
  df_pair %>% select(starts_with("var")),
  cluster = df_pair$location,
  method = "spearman"
)

# 可视化结果
corrplot(corr_result$r, p.mat = corr_result$p, sig.level = 0.05, insig = "blank")

内容的提问来源于stack exchange，提问作者Igor Salerno Filgueiras