6

いくつかの繰り返し測定された栄養摂取データ (RespondentID ごとに 2 つの 24 時間摂取期間) から構築された lme オブジェクトがあります。

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

RespondentIDを使用して、ランダム効果を正常に取得できますranef(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)AgeFactorレベルが9-13の最初のRespondentIDについて、 を使用して取得されました: - 当てはめ値は15.22633で、これは - 係数から -(Intercept) + (AgeFactor9-13) = 14.28304 + 0.9432941

各被験者の固定効果推定値を抽出するために、自動的に使用したい巧妙なコマンドがありますか、ifそれとも正しい AgeFactor レベルを各被験者に適用して取得しようとする一連のステートメントに直面していますか?切片から変量効果の寄与を差し引いた後の、正しい固定効果の推定値?

更新、申し訳ありませんが、私が提供していた出力を削減しようとしていて、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")
4

2 に答える 2

9

以下は、 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には繰り返しが含まれているため、上記のコードを少し調整する必要がある場合があります。いくつかの検証出力と、ラティスを使用した視覚化を次に示します。

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

ここに画像の説明を入力してください

于 2013-02-28T14:12:42.883 に答える
3

It is going to be something like this (although you really should have given us the results of str(Male.Data) because model output does not tell us the factor levels for the baseline values:)

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

You are basically running the original data through a match function to pick the correct coefficient(s) to add to the intercept ... which will be 0 if the data is the factor's base level (whose spelling I am guessing at.)

EDIT: I just noticed that you put a "-1" in the formula so perhaps all of your AgeFactor terms are listed in the output and you can tale out the 0 in the coefficient vector and the invented AgeFactor level in the match table vector.

于 2011-12-30T20:26:25.707 に答える