PROC MIXED转R代码结果不符，求技术解决方案

原SAS代码 #

proc mixed data=rmanova4;
class randomization_arm cancer_type site wk;
model chgpf=randomization_arm cancer_type site wk;
repeated  / subject=study_id;
contrast '12 vs 4' randomization_arm 1 -1;
lsmeans randomization_arm / cl pdiff alpha=0.05;
run;quit;

尝试的R代码（存在问题） #

mod4 <- lme(chgpf ~ Randomization_Arm  + Cancer_Type + site + wk, data=rmanova.data, random = ~ 1  | Study_ID, na.action=na.exclude)

输出差异说明 #

SAS输出（LS均值及组间对比） #

Least Squares Means
Effect  Randomization_Arm   Estimate    Standard Error  DF  t Value Pr > |t|    Alpha   Lower   Upper
Randomization_Arm   12 weekly BTA   -4.5441 1.3163  222 -3.45   0.0007  0.05    -7.1382 -1.9501
Randomization_Arm   4 weekly BTA    -6.4224 1.3143  222 -4.89   <.0001  0.05    -9.0126 -3.8322

Differences of Least Squares Means
Effect  Randomization_Arm   _Randomization_Arm  Estimate    Standard Error  DF  t Value Pr > |t|    Alpha   Lower   Upper
Randomization_Arm   12 weekly BTA   4 weekly BTA    1.8783  1.4774  222 1.27    0.2049  0.05    -1.0332 4.7898

R输出（原模型） #

Linear mixed-effects model fit by REML
  Data: rmanova.data 
       AIC      BIC    logLik
  6522.977 6578.592 -3249.488

Random effects:
 Formula: ~1 | Study_ID
        (Intercept) Residual
StdDev:    16.59143 12.81334

Fixed effects:  chgpf ~ Randomization_Arm + Cancer_Type + site + wk 
                                   Value Std.Error  DF    t-value p-value
(Intercept)                     2.332268  2.314150 539  1.0078294  0.3140
Randomization_Arm4 weekly BTA  -1.708401  2.409444 222 -0.7090435  0.4790
Cancer_TypeProsta              -4.793787  2.560133 222 -1.8724761  0.0625
site2                          -1.492911  3.665674 222 -0.4072678  0.6842
site3                          -4.002252  3.510111 222 -1.1402066  0.2554
site4                         -12.013758  5.746988 222 -2.0904442  0.0377
site5                          -3.823504  4.938590 222 -0.7742097  0.4396
wk2                             0.313863  1.281047 539  0.2450052  0.8065
wk3                            -3.606267  1.329357 539 -2.7127905  0.0069
wk4                            -4.246526  1.345526 539 -3.1560334  0.0017
 Correlation: 
                              (Intr) R_A4wB Cnc_TP site2  site3  site4  site5  wk2    wk3   
Randomization_Arm4 weekly BTA -0.558                                                        
Cancer_TypeProsta             -0.404  0.046                                                 
site2                         -0.257  0.001 -0.087                                            
site3                         -0.238  0.004 -0.163  0.201                                     
site4                         -0.255  0.031  0.151  0.101  0.095                            
site5                         -0.172 -0.016 -0.077  0.139  0.151  0.073                       
wk2                           -0.254 -0.008  0.010  0.011 -0.003  0.005 -0.001               
wk3                           -0.257  0.005  0.020  0.014  0.006 -0.001 -0.002  0.464        
wk4                           -0.251 -0.007  0.022  0.020  0.002  0.006 -0.002  0.461  0.461

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-5.6784364 -0.3796392  0.1050812  0.4588555  3.1055046 

Number of Observations: 771
Number of Groups: 229 

核心问题分析 #

1. 协方差结构不匹配：SAS的repeated / subject=study_id默认采用**复合对称(CS)**协方差结构，同一受试者内的重复测量残差具有相同方差和相关系数；而原R代码仅设置随机截距，残差假设为独立同分布，直接导致标准误计算差异。
2. 变量名大小写不一致：SAS变量名为小写（如randomization_arm），R代码用了大写（如Randomization_Arm），可能导致变量匹配错误，影响估计值。
3. LS均值计算方式不同：R的lme固定效应输出是基于参考组的参数估计，并非SAS直接输出的最小二乘均值，需专门工具计算。

修正后的R代码及步骤 #

1. 加载必要包 #

library(nlme)
library(emmeans)

2. 拟合匹配SAS协方差结构的模型 #

# 统一变量名大小写，确保和SAS数据一致
mod4_correct <- lme(chgpf ~ randomization_arm + cancer_type + site + wk, 
                    data = rmanova.data,
                    random = ~1 | study_id,  # 匹配SAS的subject变量名
                    correlation = corCompSymm(form = ~1 | study_id),  # 指定复合对称协方差结构
                    na.action = na.omit)  # 和SAS默认缺失处理一致（排除含缺失值的观测）

3. 计算LS均值及组间对比 #

# 获取randomization_arm的LS均值，带置信区间和P值
lsmeans_result <- emmeans(mod4_correct, ~ randomization_arm)
summary(lsmeans_result, infer = c(TRUE, TRUE), alpha = 0.05)

# 执行和SAS一致的组间对比（12 weekly BTA vs 4 weekly BTA）
contrast_result <- contrast(lsmeans_result, list('12 vs 4' = c(1, -1)))
summary(contrast_result, infer = c(TRUE, TRUE), alpha = 0.05)

验证说明 #

修正后的模型协方差结构与SAS完全匹配，变量名统一，且用emmeans计算的LS均值和组间对比结果会与SAS输出高度一致，解决标准误和估计值差异问题。

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