如何扩展代码合并各地点的lsmeans与CLD分析结果

需求：按location分析Factor1、Factor2对响应变量perf的影响，最终将原始数据集与各地点的lsmeans及Tukey分组字母关联。现有代码仅能处理单个地点，需扩展以合并所有地点的lsmeans和CLD输出。

原始代码 #

library(multcomp)
library(lsmeans)
library(dplyr)
library(broom)

output <- as.data.frame(dataset)
output$perf <- as.numeric(output$perf)
output$factor1 <- as.factor(output$factor1)
output$factor2 <- as.factor(output$factor2)
output <- output %>% 
           group_by(location) %>% 
           mutate(row=row_number()) 

model <- lm(perf ~ factor1 + factor2, data = output)
lsmeans <- lsmeans(model, ~ factor1)
CLD = cld(lsmeans,
            alpha   = 0.05,
            Letters = letters,
            adjust  = "tukey")         

output <- left_join(output, CLD, by = "factor1")

扩展解决方案 #

使用dplyr的group_modify实现按地点分组处理，每个地点单独建模、计算LS均值和Tukey分组，最后自动合并所有结果到原始数据集：

library(multcomp)
library(lsmeans)
library(dplyr)
library(broom)

# 数据预处理统一用mutate简化
output <- as.data.frame(dataset) %>%
  mutate(
    perf = as.numeric(perf),
    factor1 = as.factor(factor1),
    factor2 = as.factor(factor2)
  )

# 按location分组执行建模、LS均值计算、CLD生成，并关联原始数据
final_result <- output %>%
  group_by(location) %>%
  group_modify(function(group_data, group_key) {
    # 针对当前地点数据构建模型
    loc_model <- lm(perf ~ factor1 + factor2, data = group_data)
    # 计算factor1的LS均值
    loc_lsmeans <- lsmeans(loc_model, ~ factor1)
    # 生成Tukey分组字母，保留关键列
    loc_cld <- cld(loc_lsmeans,
                   alpha = 0.05,
                   Letters = letters,
                   adjust = "tukey") %>%
      as.data.frame() %>%
      select(factor1, lsmean, SE, df, lower.CL, upper.CL, .group)
    
    # 将当前地点的原始数据与分组结果关联
    group_data %>%
      left_join(loc_cld, by = "factor1")
  }) %>%
  ungroup()

# 查看合并后的结果
head(final_result)

关键说明 #

• group_modify会遍历每个location分组，对每组数据单独执行建模和统计分析，避免全局模型忽略地点差异的问题
• 若需要分析factor1与factor2的交互效应，只需修改lsmeans的公式为~ factor1 * factor2；若单独分析factor2，则改为~ factor2
• 输出结果包含原始数据的所有列，以及对应地点的LS均值、标准误、置信区间和Tukey分组字母（.group列）

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