R语言rstanarm包中分类变量的唯一截距法应用咨询

当然可以在rstanarm里复刻McElearth在《统计反思》158-159页提到的索引变量替代虚拟编码的线性回归方法！这个思路核心是用整数索引直接映射分类变量的每个水平到独立的模型参数，替代传统的虚拟编码（K-1个哑变量），下面我一步步给你演示：

1. 准备数据集并生成索引变量 #

首先我们不需要依赖rethinking包的coerce_index函数，用base R就能生成等价的索引变量——本质就是把分类变量转成从1开始的整数标识：

library(rstanarm)
# 加载milk数据集（来自rethinking包）
data(milk, package = "rethinking")
d <- milk

# 生成clade的整数索引，和coerce_index效果完全一致
d$clade_id <- as.integer(factor(d$clade))

factor(d$clade)会把原始的分类字符串转成有序因子，as.integer则提取它的水平序号，得到和coerce_index一样的索引结果。

2. 拟合对应模型 #

原书里的方法是给每个clade水平分配独立的固定效应参数（即每个组有自己的截距），对应到rstanarm里有两种实现方式，看你需要哪种：

方式一：固定效应版本（完全对应原书的索引变量思路） #

我们用无截距模型，直接让每个clade索引对应的水平作为独立参数：

# 无截距模型，每个clade组对应一个独立参数
model_fixed <- stan_lm(
  kcal.per.g ~ 0 + factor(clade_id),
  data = d,
  refresh = 0 # 关闭迭代过程输出，让结果更整洁
)

# 查看模型结果
print(model_fixed, digits = 3)

这个模型的系数就是每个clade组对应的kcal.per.g的预测均值，和原书里map模型输出的a[clade_id]参数完全对应。

方式二：分层（随机效应）版本（带正则化的索引变量思路） #

如果你想给组间参数加上收缩正则化（这也是Bayesian方法常用的技巧），可以用分层模型，语法上用(1 | clade_id)表示每个clade组有独立的随机截距：

# 分层线性模型，给clade组的截距加正则化
model_hierarchical <- stan_lmer(
  kcal.per.g ~ 1 + (1 | clade_id),
  data = d,
  refresh = 0
)

# 查看分层模型的结果
print(model_hierarchical, digits = 3)

这个版本的结果会输出全局截距和组间的随机效应偏差，和原书里如果给a[clade_id]加弱先验的效果类似。

补充说明 #

• 索引变量方法和虚拟编码的核心区别：虚拟编码是用K-1个哑变量以某个组为基准，而索引变量是给K个组各分配一个独立参数，两者可以通过线性变换互相转换，但索引变量的参数解读更直接（每个组的绝对均值）。
• 在rstanarm里，两种方式都完全支持，你可以根据是否需要正则化来选择固定或分层模型。

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