如何通过观察提取 lmer 固定效应？答案

【问题标题】：How do I extract lmer fixed effects by observation?如何通过观察提取 lmer 固定效应？
【发布时间】：2011-12-30 18:50:28
【问题描述】：

我有一个 lme 对象，由一些重复测量的营养摄入数据构成（每个 RespondentID 有两个 24 小时摄入期）：

Male.lme2 <- lmer(BoxCoxXY ~ -1 + AgeFactor + IntakeDay + (1|RespondentID),
    data = Male.Data, 
    weights = SampleWeight)

我可以使用ranef(Male.lme1) 成功检索RespondentID 的随机效果。我也想通过RespondentID收集固定效果的结果。 coef(Male.lme1) 没有提供我所需要的，如下所示。

> summary(Male.lme1)
Linear mixed model fit by REML 
Formula: BoxCoxXY ~ AgeFactor + IntakeDay + (1 | RespondentID) 
   Data: Male.Data 
  AIC   BIC logLik deviance REMLdev
  9994 10039  -4990     9952    9980
Random effects:
 Groups       Name        Variance Std.Dev.
 RespondentID (Intercept) 0.19408  0.44055 
 Residual                 0.37491  0.61230 
Number of obs: 4498, groups: RespondentID, 2249

Fixed effects:
                    Estimate Std. Error t value
(Intercept)         13.98016    0.03405   410.6
AgeFactor4to8        0.50572    0.04084    12.4
AgeFactor9to13       0.94329    0.04159    22.7
AgeFactor14to18      1.30654    0.04312    30.3
IntakeDayDay2Intake -0.13871    0.01809    -7.7

Correlation of Fixed Effects:
            (Intr) AgFc48 AgF913 AF1418
AgeFactr4t8 -0.775                     
AgeFctr9t13 -0.761  0.634              
AgFctr14t18 -0.734  0.612  0.601       
IntkDyDy2In -0.266  0.000  0.000  0.000

我已将拟合结果附加到我的数据中，head(Male.Data) 显示

   NutrientID RespondentID Gender Age SampleWeight  IntakeDay IntakeAmt AgeFactor BoxCoxXY  lmefits
2         267       100020      1  12    0.4952835 Day1Intake 12145.852     9to13 15.61196 15.22633
7         267       100419      1  14    0.3632839 Day1Intake  9591.953    14to18 15.01444 15.31373
8         267       100459      1  11    0.4952835 Day1Intake  7838.713     9to13 14.51458 15.00062
12        267       101138      1  15    1.3258785 Day1Intake 11113.266    14to18 15.38541 15.75337
14        267       101214      1   6    2.1198688 Day1Intake  7150.133      4to8 14.29022 14.32658
18        267       101389      1   5    2.1198688 Day1Intake  5091.528      4to8 13.47928 14.58117

coef(Male.lme1) 的前几行是：

$RespondentID
       (Intercept) AgeFactor4to8 AgeFactor9to13 AgeFactor14to18 IntakeDayDay2Intake
100020    14.28304     0.5057221      0.9432941        1.306542          -0.1387098
100419    14.00719     0.5057221      0.9432941        1.306542          -0.1387098
100459    14.05732     0.5057221      0.9432941        1.306542          -0.1387098
101138    14.44682     0.5057221      0.9432941        1.306542          -0.1387098
101214    13.82086     0.5057221      0.9432941        1.306542          -0.1387098
101389    14.07545     0.5057221      0.9432941        1.306542          -0.1387098

为了演示 coef 结果与 Male.Data 中的拟合估计值之间的关系（使用 Male.Data$lmefits <- fitted(Male.lme1) 抓取，对于 AgeFactor 级别为 9-13 的第一个 RespondentID： - 拟合值为15.22633，等于 - 从系数 - (Intercept) + (AgeFactor9-13) = 14.28304 + 0.9432941

我是否有一个聪明的命令可以自动使用，即提取每个主题的固定效应估计值，或者我是否面临一系列 if 语句试图应用正确的 AgeFactor在从 Intercept 中扣除随机效应贡献后，对每个受试者进行水平以获得正确的固定效应估计？

更新，抱歉，试图减少我提供的输出并忘记了 str()。输出是：

>str(Male.Data)
'data.frame':   4498 obs. of  11 variables:
 $ NutrientID  : int  267 267 267 267 267 267 267 267 267 267 ...
 $ RespondentID: Factor w/ 2249 levels "100020","100419",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ Gender      : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Age         : int  12 14 11 15 6 5 10 2 2 9 ...
 $ BodyWeight  : num  51.6 46.3 46.1 63.2 28.4 18 38.2 14.4 14.6 32.1 ...
 $ SampleWeight: num  0.495 0.363 0.495 1.326 2.12 ...
 $ IntakeDay   : Factor w/ 2 levels "Day1Intake","Day2Intake": 1 1 1 1 1 1 1 1 1 1 ...
 $ IntakeAmt   : num  12146 9592 7839 11113 7150 ...
 $ AgeFactor   : Factor w/ 4 levels "1to3","4to8",..: 3 4 3 4 2 2 3 1 1 3 ...
 $ BoxCoxXY    : num  15.6 15 14.5 15.4 14.3 ...
 $ lmefits     : num  15.2 15.3 15 15.8 14.3 ...

未使用 BodyWeight 和 Gender（这是男性数据，因此所有 Gender 值都相同），并且 NutrientID 对数据同样是固定的。

自从我发布以来，我一直在做可怕的 ifelse 语句，所以会立即尝试你的建议。 :)

Update2：这与我当前的数据完美配合，并且对于新数据应该是面向未来的，感谢 DWin 在评论中提供的额外帮助。 :)

AgeLevels <- length(unique(Male.Data$AgeFactor))
Temp <- as.data.frame(fixef(Male.lme1)['(Intercept)'] + 
c(0,fixef(Male.lme1)[2:AgeLevels])[
      match(Male.Data$AgeFactor, c("1to3", "4to8", "9to13","14to18",  "19to30","31to50","51to70","71Plus") )] + 
c(0,fixef(Male.lme1)[(AgeLevels+1)])[
      match(Male.Data$IntakeDay, c("Day1Intake","Day2Intake") )])
names(Temp) <- c("FxdEffct")

【问题讨论】：

标签： r glm

【解决方案1】：

以下是我一直认为最容易在 lme4 包中提取个人的固定效果和随机效果组件的方法。它实际上为每个观察提取了相应的拟合。假设我们有一个形式的混合效应模型：

y = Xb + Zu + e

其中Xb是固定效应，Zu是随机效应，我们可以提取成分（以lme4的sleepstudy为例）：

library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)

# Xb 
fix <- getME(fm1,'X') %*% fixef(fm1)
# Zu
ran <- t(as.matrix(getME(fm1,'Zt'))) %*% unlist(ranef(fm1))
# Xb + Zu
fixran <- fix + ran

我知道这是一种从线性混合效应模型中提取组件的通用方法。对于非线性模型，模型矩阵 X 包含重复，您可能需要稍微修改上述代码。这是一些验证输出以及使用 lattice 的可视化：

> head(cbind(fix, ran, fixran, fitted(fm1)))
         [,1]      [,2]     [,3]     [,4]
[1,] 251.4051  2.257187 253.6623 253.6623
[2,] 261.8724 11.456439 273.3288 273.3288
[3,] 272.3397 20.655691 292.9954 292.9954
[4,] 282.8070 29.854944 312.6619 312.6619
[5,] 293.2742 39.054196 332.3284 332.3284
[6,] 303.7415 48.253449 351.9950 351.9950

# Xb + Zu
> all(round((fixran),6) == round(fitted(fm1),6))
[1] TRUE

# e = y - (Xb + Zu)
> all(round(resid(fm1),6) == round(sleepstudy[,"Reaction"]-(fixran),6))
[1] TRUE

nobs <- 10 # 10 observations per subject
legend = list(text=list(c("y", "Xb + Zu", "Xb")), lines = list(col=c("blue", "red", "black"), pch=c(1,1,1), lwd=c(1,1,1), type=c("b","b","b")))
require(lattice)
xyplot(
    Reaction ~ Days | Subject, data = sleepstudy,
    panel = function(x, y, ...){
        panel.points(x, y, type='b', col='blue')
        panel.points(x, fix[(1+nobs*(panel.number()-1)):(nobs*(panel.number()))], type='b', col='black')
        panel.points(x, fixran[(1+nobs*(panel.number()-1)):(nobs*(panel.number()))], type='b', col='red')
    },
    key = legend
)

【讨论】：

这太棒了，除了 fixran 似乎并不总是与 lme4 1.1-12 很好的近似。你可以尝试复制吗？
y和Xb+Zu有什么区别？不应该一样吗？
and：随机效果不应该只有Zu而不是Xb+Zu吗？
抱歉，距离这个答案已经四年多了，我会尝试回到这个@smci 看看有什么变化。 @本； 1) 没错，它们是相同的，或者更准确地说，y = Xb+Zu+e。我用两种方法提取了这两个向量；原始数据，以及模型项的总和，以说明这两者确实是相等的。 2）没错，严格来说，Zu 是唯一一个随机效应的比例。但是，如果您忽略 Xb，您只剩下零均值正态分布随机效应的总和，由于总体平均值被排除在外，因此更难解释。因此，我在这里补充它们。

【解决方案2】：

这将是这样的（虽然你真的应该给了我们 str(Male.Data) 的结果，因为模型输出没有告诉我们基线值的因子水平： )

#First look at the coefficients
fixef(Male.lme2)

#Then do the calculations
fixef(Male.lme2)[`(Intercept)`] + 
c(0,fixef(Male.lme2)[2:4])[
          match(Male.Data$AgeFactor, c("1to3", "4to8", "9to13","14to18") )] + 
c(0,fixef(Male.lme2)[5])[
          match(Male.Data$IntakeDay, c("Day1Intake","Day2Intake") )]

您基本上是通过match 函数运行原始数据以选择正确的系数以添加到截距...如果数据是因子的基本水平（我猜是其拼写），则该系数将为 0在。）

编辑：我刚刚注意到您在公式中添加了“-1”，因此您的所有 AgeFactor 项可能都列在输出中，您可以找出系数向量中的 0 和匹配中发明的 AgeFactor 级别表向量。

【讨论】：

感谢您的帮助，我刚刚修改了 (Intercept) 名称周围的引号。我正在创建一个适用于所有年龄组的通用 R 分析，当前数据框只有孩子，当我不一定知道模型中有多少年龄因素水平时，如何调整搜索的列索引？我正在尝试尽可能地自动化分析
length(unique(Male.Data$AgeFactor)) 将为您提供级别数，您可以使用该数字加 1 而不是 4 来获取 AgeFactor 的索引，您显然需要添加使用适当更高的值IntakeDay 效应的指数也是如此。