Come estrarre gli effetti fissi lm dall'osservazione?

Question

Ho un oggetto LME, costruito da alcune misure ripetute dati di assunzione dei nutrienti (due periodi di aspirazione di 24 ore per RespondentID):

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

e posso recuperare con successo gli effetti casuali di RespondentID utilizzando ranef(Male.lme1). Vorrei anche raccogliere il risultato degli effetti fissi per RespondentID. coef(Male.lme1) non fornisce esattamente ciò di cui ho bisogno, come mostrato di seguito.

> 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 

ho aggiunto i risultati montati ai miei dati, head(Male.Data) mostra

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 

La prima coppia di linee da coef(Male.lme1) sono:

$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 

Per dimostrare come i coef risultati riguardano le stime fornite in Male.Data (che sono state afferrate usando Male.Data$lmefits <- fitted(Male.lme1), per il primo RespondentID, che ha l'età Livello di fattore 9-13: - il valore stimato è 15.22633, che equivale - dalle coeffs - (Intercept) + (AgeFactor9-13) = 14.28304 + 0.9432941

C'è un comando intelligente per me di utilizzare che farà voglio voglio automaticamente, che è quello di estrarre l'effetto fisso stima per ogni argomento o mi trovo di fronte a una serie di affermazioni if che tentano di applicare il livello AgeFactor corretto a ciascun soggetto per ottenere la stima dell'effetto fisso corretta, dopo aver dedotto il contributo dell'effetto casuale dall'Intercept?

Aggiornamento, scuse, stava cercando di ridurre l'output che stavo fornendo e ho dimenticato di str(). Output è:

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

Il peso corporeo e sesso non vengono utilizzati (questo è il dato maschi, quindi tutti i valori genere sia lo stesso) e la NutrientID analogamente viene fissato per i dati.

Ho fatto orribili dichiarazioni ifelse sinceramente che ho postato, quindi proverò immediatamente il tuo suggerimento. :)

Update2: questo funziona perfettamente con i miei dati attuali e dovrebbe essere a prova di futuro per i nuovi dati, grazie a DWin per l'ulteriore aiuto nel commento per questo. :)

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

Si sta per essere qualcosa di simile (anche se in realtà dovrebbe ci hai visti i risultati di str (Male.Data) perché il risultato del modello non ci dice i livelli dei fattori per i valori basali :)

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

Tu sei fondamentalmente esegue i dati originali attraverso una funzione match di scegliere il coefficiente corretto (s) da aggiungere alla intercettazione ... che sarà 0 se i dati è il livello base del fattore (la cui ortografia Sto indovinando a.)

EDIT: Ho appena notato che hai inserito un "-1" nella formula, quindi forse tutti i termini di AgeFactor sono elencati nell'output e puoi scoprire lo 0 nel vettore del coefficiente e il livello di AgeFactor inventato nella tabella delle corrispondenze vettore.

Answer 2

Di seguito è come ho sempre trovato più semplice estrapolare gli effetti fissi e componenti a effetti casuali degli individui nel pacchetto lme4. In realtà estrae l'adattamento corrispondente ad ogni osservazione. Supponendo che abbiamo un modello misto effetti della forma:

y = Xb + Zu + e 

dove Xb sono gli effetti fissi e Zu sono gli effetti casuali, è possibile estrarre i componenti (usando lme4 s' sleepstudy come esempio):

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 

So che questo funziona come un approccio generalizzato all'estrazione di componenti da modelli ad effetti misti lineari. Per i modelli non lineari, la matrice del modello X contiene ripetizioni e potrebbe essere necessario adattare un po 'il codice precedente. Ecco qualche uscita di convalida e una visualizzazione con reticolo:

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