在R语言中如何将参与者作为随机效应纳入对数线性/泊松模型？

我希望将参与者作为随机效应纳入对数线性模型（泊松模型），以此解释参与者间的方差。我的数据来自列联表，是基于观测值预测期望值的场景，因此无法使用R中的常规混合效应模型（如lme4包的函数）。

简化版数据（仅2名参与者） #

我需要测试awareness与narrative的交互效应，同时将参与者作为随机效应处理。当前的代码未纳入随机效应，请问如何实现？

原代码：

observed <- c(7, 2, 8, 22, 0, 0, 14, 2, 8, 22, 0, 0)
awareness <- factor(c("aware","aware","aware","unaware","unaware","unaware","aware","aware","aware","unaware","unaware","unaware"))
narrative <- factor(c("M1","MU2","MU3","M1","MU2","MU3","M1","MU2","MU3","M1","MU2","MU3"))
# participant <- factor(c("1", "1", "1","1", "1", "1","2", "2", "2","2", "2", "2"))
model1 <- glm(observed~awareness*narrative, poisson)
model2 <- glm(observed~awareness+narrative, poisson)
anova(model1, model2, test="Chi")
summary(model2)

四维列联表数据场景 #

同样的问题，如何将参与者作为随机效应纳入以下模型？

原代码：

y <- structure(c(7.33, 13.5, 7.5, 8.3, 3.16, 2, 3.3, 5.3, 9.5, 4.2, 9.2, 6.5, 8.33, 14.5, 8.5, 9.3, 4.16, 3, 4.3, 6.3, 10.5, 5.2, 10.2, 7.5), 
               dim = c(2L, 2L, 2L, 3L), 
               dimnames = structure(list(c("ParticipantOne", "ParticipantTwo"), 
                                         c("above threshold", "below threshold"), 
                                         c("Healthy controls", "Psychotic observers"), 
                                         c("M1", "MU2", "MU3")), 
                                    names = c("Awareness", "Group", "Narrative")), 
               class = "table")
xy <- aperm(y, c(2, 2, 1, 3))
names(dimnames(xy)) <- c("Participant" ,"Group", "Awareness", "Narrative")
ftable(xy)
fourfoldplot(xy, margin = 2)

narrative <- gl(3,4)
group <- gl(2,1,12)
awareness <- gl(2,2,12)
participant <- gl(???)  # 这里不知道怎么构造
model1 <- glm(as.vector(xy) ~narrative*group*awareness, poisson)
model2 <- update(model1, ~. -narrative:group:awareness)
anova(model1, model2, test="Chi")

解决方案 #

核心思路：使用glmmTMB包拟合带随机效应的泊松对数线性模型

你提到的“不存在待预测变量y”其实是误解：列联表中的观测频数就是模型的响应变量，对数线性模型本质就是以频数为因变量的泊松回归，因此完全可以用混合效应泊松模型拟合。glmmTMB在处理计数数据和复杂随机结构时比lme4更灵活，适合这类场景。

1. 简化版数据的实现

补全participant变量后，直接用glmmTMB加入随机截距：

# 加载包
library(glmmTMB)

# 补全参与者变量
observed <- c(7, 2, 8, 22, 0, 0, 14, 2, 8, 22, 0, 0)
awareness <- factor(c("aware","aware","aware","unaware","unaware","unaware","aware","aware","aware","unaware","unaware","unaware"))
narrative <- factor(c("M1","MU2","MU3","M1","MU2","MU3","M1","MU2","MU3","M1","MU2","MU3"))
participant <- factor(c("1", "1", "1","1", "1", "1","2", "2", "2","2", "2", "2"))

# 带参与者随机截距的交互模型
model1_re <- glmmTMB(observed ~ awareness*narrative + (1|participant), 
                     family = poisson)
# 无交互的模型
model2_re <- glmmTMB(observed ~ awareness+narrative + (1|participant), 
                     family = poisson)

# 模型比较
anova(model1_re, model2_re, test="Chisq")
summary(model2_re)

2. 四维列联表数据的实现

将列联表转换为长格式数据框，避免手动构造变量出错，再拟合模型：

library(glmmTMB)
library(tidyr)

# 处理四维表数据
y <- structure(c(7.33, 13.5, 7.5, 8.3, 3.16, 2, 3.3, 5.3, 9.5, 4.2, 9.2, 6.5, 8.33, 14.5, 8.5, 9.3, 4.16, 3, 4.3, 6.3, 10.5, 5.2, 10.2, 7.5), 
               dim = c(2L, 2L, 2L, 3L), 
               dimnames = structure(list(c("ParticipantOne", "ParticipantTwo"), 
                                         c("above threshold", "below threshold"), 
                                         c("Healthy controls", "Psychotic observers"), 
                                         c("M1", "MU2", "MU3")), 
                                    names = c("Awareness", "Group", "Narrative")), 
               class = "table")

# 转换为长格式数据框，修正变量命名
long_data <- as.data.frame(y) %>%
  rename(Participant = Awareness, 
         Group = Group, 
         Awareness = Narrative, 
         Freq = Freq) %>%
  mutate(across(c(Participant, Group, Awareness, Narrative), factor))

# 带参与者随机截距的三维交互模型
model1_re_4d <- glmmTMB(Freq ~ Narrative*Group*Awareness + (1|Participant), 
                        family = poisson, data = long_data)
# 去掉三维交互的模型
model2_re_4d <- update(model1_re_4d, ~. -Narrative:Group:Awareness)

# 模型比较
anova(model1_re_4d, model2_re_4d, test="Chisq")
summary(model2_re_4d)

关键说明

• (1|participant)表示给每个参与者添加随机截距，用来解释参与者间的个体差异。如果样本量足够（比如参与者数量多于5个），可以尝试更复杂的随机结构，比如(Awareness|Participant)表示随参与者变化的Awareness效应斜率，但简化版数据只有2个参与者，仅能支持随机截距。
• 若数据存在过度离散（泊松模型残差方差远大于均值），可以将family参数改为nbinom2或nbinom1（负二项分布），glmmTMB直接支持这类调整。

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