能否在BradleyTerry2 R包中计算每名被试对每个刺激的能力值？

我开展了一项任务，让被试对12个刺激进行两两比较，所有被试都会完成全部比较，且比较顺序完全随机。后续分析中，我需要提取每名被试针对每个刺激的能力值，请问是否可行？

现有未考虑被试因素的模型 #

model <- BTm(chosen_left, left_stimulus, right_stimulus, data = data) 

查阅文档后的疑问 #

我已查阅BradleyTerry2包的文档，其中有使用CEMS数据的示例，但该示例是纳入被试水平协变量，而非仅设置被试随机截距。且运行该示例模型后，BTabilities()函数仍仅输出每个“玩家”的整体能力值：

data("CEMS", package = "BradleyTerry2")

table8.model <- BTm(outcome = cbind(win1.adj, win2.adj),
player1 = school1, player2 = school2, formula = ~ .. +
WOR[student] * LAT[..] + DEG[student] * St.Gallen[..] +
STUD[student] * Paris[..] + STUD[student] * St.Gallen[..] + ENG[student] * St.Gallen[..] + FRA[student] * London[..] +
FRA[student] * Paris[..] + SPA[student] * Barcelona[..] +
ITA[student] * London[..] + ITA[student] * Milano[..] +
SEX[student] * Milano[..],
refcat = "Stockholm", data = CEMS)

完全可以提取每名被试对每个刺激的能力值，核心是在模型中纳入被试与刺激的关联效应，以下是两种实操思路：

思路1：分层混合效应BT模型（推荐） #

通过BradleyTerry2结合混合效应框架，为每个被试设置针对刺激的随机偏差，既利用群体信息稳定估计，又能得到个体层面的能力值：

library(BradleyTerry2)
library(lme4)

# 拟合包含被试-刺激随机效应的模型
mixed_bt_model <- BTm(
  chosen_left, left_stimulus, right_stimulus,
  # 公式中(1 | subject_id)表示被试层面的随机截距，与刺激关联
  formula = ~ .. + (1 | subject_id),
  data = data  # data需包含subject_id（被试ID变量）
)

# 提取刺激的固定效应能力值
stim_fixed_abilities <- BTabilities(mixed_bt_model)
# 提取每个被试对各刺激的随机偏差
subject_random_deviations <- ranef(mixed_bt_model)$subject_id

# 计算每名被试对每个刺激的最终能力值
subject_stim_abilities <- expand.grid(
  subject_id = rownames(subject_random_deviations),
  stimulus = names(stim_fixed_abilities)
)

subject_stim_abilities$ability <- apply(subject_stim_abilities, 1, function(row) {
  stim_fixed_abilities[row$stimulus] + subject_random_deviations[row$subject_id, row$stimulus]
})

思路2：为每个被试单独拟合BT模型 #

如果被试数量不多，可直接为每个被试单独拟合模型，直接得到该被试视角下的刺激能力值：

library(BradleyTerry2)

# 按被试拆分数据集
subject_data_list <- split(data, data$subject_id)

# 循环拟合模型并提取能力值
subject_stim_results <- lapply(subject_data_list, function(subj_data) {
  subj_model <- BTm(chosen_left, left_stimulus, right_stimulus, data = subj_data)
  data.frame(
    subject_id = unique(subj_data$subject_id),
    stimulus = names(BTabilities(subj_model)),
    ability = as.numeric(BTabilities(subj_model))
  )
})

# 合并所有被试的结果
final_results <- do.call(rbind, subject_stim_results)

关键注意点 #

• 分层模型（思路1）更高效，能通过群体信息收缩极端估计值，适合被试数量较多的场景；
• 独立模型（思路2）更直观，但当被试比较次数不足时（你这里是全比较，无此问题），估计结果可能不稳定；
• BTabilities()默认会将某一刺激设为参考类别（能力值为0），若需要绝对能力值，可自行标准化（例如将所有刺激能力值的均值设为0）。

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