如何在 merMod 对象中找到具有交互项的原始变量（配备 lme4 包）？

Question

我想计算交互效应，如this blog posting 所示，所以我可以计算y = (b0 + (b1 * xa) + (b3 * xa * xb))，其中b0 是截距的估计值，b1 是预测变量 A 的估计值，@987654325 @ 是预测变量 A 和 B 之间交互作用的估计值。

为此，我需要

1. 拟合模型的估计值
2. 预测变量的原始值（即所用变量的“data.frame”，其中每个值都插入到上面显示的公式中）

对于“简单”的广义线性模型，我可以使用数据矩阵参数将原始数据值存储在拟合模型对象中（glm(... x=TRUE)。示例：

fit <- glm(f1care ~ c12hour + neg_c_7 + c172 + 
           sex1 + sex2 + c172:neg_c_7 + sex1:c12hour,
           data = mydf,
           x = TRUE,
           family = binomial("logit"))

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)    1.4674300  0.1759883   8.338  < 2e-16 ***
c12hour        0.0016157  0.0016315   0.990   0.3220    
neg_c_7       -0.0549614  0.0108868  -5.048 4.45e-07 ***
c1721         -0.1575547  0.2752344  -0.572   0.5670    
sex12         -0.0320619  0.1286984  -0.249   0.8033    
sex22          0.2091354  0.0874629   2.391   0.0168 *  
neg_c_7:c1721  0.0090043  0.0202720   0.444   0.6569    
c12hour:sex12 -0.0003121  0.0018246  -0.171   0.8642  

> head(fit$x)
  (Intercept) c12hour neg_c_7 c1721 sex12 sex22 neg_c_7:c1721 c12hour:sex12
1           1      10       9     0     1     1             0            10
2           1       4      13     1     1     1            13             4
3           1      12      21     0     1     1             0            12
4           1      60      14     0     1     1             0            60
5           1      40      17     0     1     1             0            40
6           1      50      17     0     1     1             0            50

如您所见，系数的名称与模型数据框中的列名称相同 (fit$x)。

如果我找到一个交互项（例如c12hour:sex12），我可以在冒号处拆分名称，拥有c12hour 和sex12，并通过模型矩阵列名称查找值（名称相同）。

现在我的问题是，如何使用 merMod 对象执行此操作？ 用于将模型 (fit@frame) 与原始值拟合的数据框的列名称如下所示：

library(lme4)
fit <- glmer(f1care ~ c12hour + neg_c_7 + c172 + 
             sex1 + sex2 + 
             c172:neg_c_7 + sex1:c12hour + (1|g2ctry),
             data = mydf,
             family = binomial("logit"))

> colnames(fit@frame)[-1]
[1] "c12hour" "neg_c_7" "c172"    "sex1"    "sex2"    "g2ctry" 

通过summary 或fixef 访问的系数如下所示：

Fixed effects:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)    1.2388127  0.2382648   5.199    2e-07 ***
c12hour        0.0018070  0.0016450   1.098 0.272018    
neg_c_7       -0.0387817  0.0115511  -3.357 0.000787 ***
c1721          0.1187357  0.2842743   0.418 0.676181    
sex12         -0.0305578  0.1306499  -0.234 0.815069    
sex22          0.1580400  0.0897806   1.760 0.078358 .  
neg_c_7:c1721 -0.0106958  0.0209961  -0.509 0.610458    
c12hour:sex12 -0.0009486  0.0018350  -0.517 0.605206   

> cbind(fixef(fit))
                       [,1]
(Intercept)    1.2388126568
c12hour        0.0018069626
neg_c_7       -0.0387817065
c1721          0.1187357405
sex12         -0.0305578499
sex22          0.1580400407
neg_c_7:c1721 -0.0106958049
c12hour:sex12 -0.0009485618

如您所见，列名 sex1 和 sex2 与系数“名称”sex12 和 sex22 不同（与第一个简单的 glm 示例不同，两者相同)。

我目前的做法是

• 取列名（例如sex1）
• 检查是否为factor，如果是，则检索levels（例如1和2）
• concat 列名和因子级别（所以我们有sex11 和sex12）
• 找到包含连接“预测器名称”的交互项（因此sex11 将不匹配，但是sex12 将匹配第二个交互项c12hour:sex12）

我想知道这是否是最好的或唯一的解决方案，或者在某些情况下我的“检测”交互术语的方式不起作用？

Answer 1

好的，model.matrix 成功了，并返回了一个数据框，其中列名对应于模型摘要中的术语名称。