如何通過觀察提取lmer固定效果？

Question

我有一個LME對象，從一些重複測量構造營養攝入數據（每RespondentID兩個24小時攝取期間）：

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

，我可以成功地檢索通過使用RespondentIDranef(Male.lme1)隨機效應。我也想通過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)爲第一個RespondentID抓取，該年齡爲Age因子水平9-13： - 擬合值是15.22633，它等於 - 從coeffs - (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 ... 

的體重和性別不被使用（這是男性數據，因此，所有的性別值相同）與NutrientID同樣固定的數據。

我一直在做可怕的ifelse陳述sinced我張貼，所以會立即嘗試你的建議。 :)

更新2：這與我目前的數據完美合作，並應該是未來的新數據，這要歸功於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") 

Answer 1

這是怎麼回事，儘管你真的應該給我們的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條款都列在輸出中，你可以在係數向量和發明的AgeFactor級別匹配表中輸出0向量。

Answer 2

下面是我總是發現在包裝中提取個人固定效果和隨機效果組件最容易的方法。它實際上爲每個觀察提取相應的擬合。假設我們有形式的混合效應模型：

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包含重複項，您可能需要修改上述代碼。下面是一些驗證輸出，以及一個可視化的格子：

> 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 
)