求助：SAS多向列联表分析转R后无法生成一致输出，求正确脚本

我需要将以下SAS代码转换为R脚本，要求输出与SAS完全一致（包括Type III Wald检验结果、预测值）。SAS结果已验证正确，但之前尝试的R实现无法匹配输出。

原SAS代码 #

data new; 
input prey $ stage $ habitat $ number;
cards ; 
chr a l 69
chr j p 8
chr j l 43
chr a p 21
amp a l 29
amp j p 1
amp j l 19
amp a p 6
zoo a l 6
zoo j p 28
zoo j l 3
zoo a p 72
ost a l 10
ost j p 4
ost j l 6
ost a p 4
;

proc genmod data=new; 
class prey stage habitat/param=effect; 
model number=prey|stage|habitat / dist=poisson type3 wald;
output out=out p=p; 
run;

proc sort data=out; 
by prey stage habitat; 
run;

proc print data=out; 
run;

R解决方案 #

要匹配SAS的输出，核心是对齐编码规则、检验类型和预测逻辑，以下是完整实现：

# 1. 导入数据集
new <- data.frame(
  prey = c("chr", "chr", "chr", "chr", "amp", "amp", "amp", "amp",
           "zoo", "zoo", "zoo", "zoo", "ost", "ost", "ost", "ost"),
  stage = c("a", "j", "j", "a", "a", "j", "j", "a",
            "a", "j", "j", "a", "a", "j", "j", "a"),
  habitat = c("l", "p", "l", "p", "l", "p", "l", "p",
              "l", "p", "l", "p", "l", "p", "l", "p"),
  number = c(69, 8, 43, 21, 29, 1, 19, 6, 6, 28, 3, 72, 10, 4, 6, 4)
)

# 2. 设置效应编码（匹配SAS的param=effect）
options(contrasts = c("contr.sum", "contr.poly"))

# 3. 拟合包含三阶交互项的泊松模型
model <- glm(number ~ prey * stage * habitat, family = poisson(link = "log"), data = new)

# 4. 计算Type III Wald检验（匹配SAS的type3 wald）
library(car)
Anova(model, type = "III", test.statistic = "Wald")

# 5. 生成预测值并排序打印
out <- cbind(new, p = predict(model, type = "response"))
out_sorted <- out[order(out$prey, out$stage, out$habitat), ]
print(out_sorted)

关键匹配点说明 #

• 效应编码：SAS的param=effect对应R的contr.sum偏差编码，必须设置此规则才能让模型参数估计与SAS一致。
• Type III检验：R默认anova()输出Type II检验结果，需用car::Anova()指定type="III"和test.statistic="Wald"来对齐SAS的检验逻辑。
• 预测值：predict(model, type="response")生成响应尺度的预测值，对应SAS的output p=p。
• 排序打印：用order()按指定变量排序，匹配SAS的proc sort和proc print流程。

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