如何在R中用参考缩放平均法估计部分重复设计的Sequence效应均值

问题对应需求解析 #

你需要执行的计算对应SAS语句：estimate 'average' intercept 3 Sequence 1 1 1/devisor=3;，核心是对Sequence效应的均值做加权平均估计。在部分重复设计的生物等效性研究中，参考缩放平均法需匹配FDA孕酮指南的分析逻辑，以下是R中的落地实现。

R实现步骤 #

1. 数据与包准备 #

确保数据集包含Subject（受试者）、Sequence（序列）、Period（周期）、Treatment（处理组）、Response（响应变量，如药代参数），加载所需包：

library(nlme)
library(emmeans)

2. 拟合匹配SAS逻辑的混合模型 #

采用nlme拟合与SAS PROC MIXED对齐的模型，保留嵌套随机效应和固定效应结构：

# 模型结构：固定效应含Sequence/Period/Treatment，随机效应为嵌套的Subject(Sequence)
fit <- lme(Response ~ Sequence + Period + Treatment,
           random = ~1 | Subject/Sequence,
           data = your_dataset,
           control = lmeControl(opt = "optim"))

3. 手动计算Sequence效应加权均值 #

直接提取模型系数并按需求计算，对应SAS的ESTIMATE语句：

# 提取Sequence对应的固定效应系数
seq_coefs <- coef(fit)$fixed[grepl("Sequence", names(coef(fit)$fixed))]

# 计算加权平均（权重1:1:1，除数3）
avg_seq_effect <- sum(seq_coefs * c(1,1,1)) / 3
names(avg_seq_effect) <- "Average Sequence Effect"

# 输出结果
avg_seq_effect

4. 用emmeans简化实现（更直观） #

借助emmeans包直接生成Sequence的边际均值，再计算加权平均：

# 获取Sequence的边际均值
seq_means <- emmeans(fit, ~ Sequence)

# 计算指定权重的平均均值，对应SAS的ESTIMATE逻辑
avg_seq <- emmeans(seq_means, "grandmean", weights = c(1,1,1)/3)
avg_seq

关键匹配SAS的细节 #

• 随机效应结构严格使用Subject/Sequence嵌套，确保与SAS的RANDOM Subject(Sequence)一致；
• 若需匹配SAS的自由度计算（DDFM=SATTERTHWAITE），可改用lmerTest包拟合模型并指定ddf="Satterthwaite"；
• 若响应变量为对数转换（药代研究常见），拟合模型时需先转换，结果反转换需遵循FDA指南要求。

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