Làm thế nào để trích xuất các hiệu ứng cố định lmer bằng cách quan sát?

Question

Tôi có một đối tượng LME, được xây dựng từ một số biện pháp lặp đi lặp lại dữ liệu hấp thụ chất dinh dưỡng (hai giai đoạn tiêu thụ 24 giờ mỗi RespondentID):

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

và tôi thành công có thể lấy lại những tác động ngẫu nhiên bằng cách sử dụng RespondentIDranef(Male.lme1). Tôi cũng muốn thu thập kết quả của các hiệu ứng cố định bằng RespondentID. coef(Male.lme1) không cung cấp chính xác những gì tôi cần, như tôi hiển thị bên dưới.

> 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 

Tôi đã gắn kết quả được trang bị để dữ liệu của tôi, head(Male.Data) lãm

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 

Cặp đôi đầu tiên của dòng từ coef(Male.lme1) là:

$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 

Để chứng minh như thế nào coef kết quả liên quan đến các ước tính được trang bị trong Male.Data (được lấy bằng cách sử dụng Male.Data$lmefits <- fitted(Male.lme1), cho RespondentID đầu tiên, người có Độ tuổi mức yếu tố 9-13: - giá trị được trang bị là 15.22633, bằng - từ coeffs - (Intercept) + (AgeFactor9-13) = 14.28304 + 0.9432941

Có một câu lệnh thông minh cho tôi để sử dụng mà sẽ làm muốn tôi muốn tự động, mà là để trích xuất ảnh hưởng ổn định ước tính cho mỗi chủ đề, hoặc tôi phải đối mặt với một loạt các câu lệnh if cố gắng áp dụng mức độ AgeFactor chính xác cho từng đối tượng để có được ước tính hiệu quả cố định chính xác, sau khi trừ phần đóng góp ngẫu nhiên khỏi Intercept?

Cập nhật, xin lỗi, đã cố gắng cắt giảm sản lượng mà tôi đã cung cấp và quên về str(). Đầu ra là:

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

Các trọng lượng cơ thể và giới tính không được sử dụng (đây là dữ liệu nam giới, vì vậy tất cả các giá trị giới đều giống nhau) và NutrientID được tương tự cố định cho dữ liệu.

Tôi đã thực hiện các câu lệnh ifelse khủng khiếp mà tôi đã đăng, vì vậy hãy thử đề xuất của bạn ngay lập tức. :)

Update2: hoạt động này hoàn hảo với dữ liệu hiện tại của tôi và phải tương lai cho dữ liệu mới, nhờ DWin để được trợ giúp thêm trong nhận xét về điều này. :)

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

Nó sẽ là một cái gì đó như thế này (mặc dù bạn thực sự nên đã ban cho chúng ta kết quả của str (Male.Data) vì mô hình đầu ra không cho chúng tôi biết mức độ yếu tố cho các giá trị cơ bản :)

#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"))] 

Bạn đang về cơ bản chạy các dữ liệu ban đầu thông qua một hàm match để chọn hệ số đúng (s) để thêm vào đánh chặn ... mà sẽ là 0 nếu dữ liệu là mức độ cơ bản của yếu tố (mà chính tả tôi đang đoán.)

EDIT: Tôi chỉ nhận thấy rằng bạn đặt một "-1" trong công thức vì vậy có lẽ tất cả các điều khoản AgeFactor của bạn được liệt kê trong đầu ra và bạn có thể tale ra 0 trong vector hệ số và mức độ AgeFactor được phát minh trong bảng đối sánh vector.

Answer 2

Dưới đây là cách tôi luôn tìm thấy cách dễ nhất để trích xuất các hiệu ứng cố định của cá nhân và các thành phần hiệu ứng ngẫu nhiên trong gói lme4. Nó thực sự chiết xuất sự phù hợp tương ứng với mỗi quan sát. Giả sử chúng ta có một mô hình hỗn hợp tác dụng của hình thức:

y = Xb + Zu + e 

nơi XB là những hiệu ứng cố định và Zu là những hiệu ứng ngẫu nhiên, chúng ta có thể trích xuất các thành phần (sử dụng lme4 's sleepstudy là một ví dụ):

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 

Tôi biết rằng điều này hoạt động như một cách tiếp cận tổng quát để trích xuất các thành phần từ mô hình hiệu ứng hỗn hợp tuyến tính. Đối với các mô hình phi tuyến tính, ma trận mô hình X chứa lặp lại và bạn có thể phải điều chỉnh mã trên một chút. Dưới đây là một số lượng xác nhận cũng như một hình dung sử dụng mạng:

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