How to extract fixed lmer effects when observed?

Question

I have an lme object built from several repetitive data on nutrient intake (two 24-hour intake periods on the RespondentID):

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

and I can successfully get random effects using the RespondentID using ranef(Male.lme1) . I would also like to get the result of fixed effects on the RespondentID . coef(Male.lme1) does not give exactly what I need, as shown below.

 > 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 

I attached the results to my data, head(Male.Data) shows

  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 

The first two lines from 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 

To demonstrate how the coef results relate to the established grades in Male.Data (which were grabbed using Male.Data$lmefits <- fitted(Male.lme1) , for the first RespondentID, whose level is 9-13 AgeFactor: - the set value is 15.22633 , which is equal to - from the coefficients - (Intercept) + (AgeFactor9-13) = 14.28304 + 0.9432941

Is there a smart command to use that will be needed, I want to automatically extract a fixed effect estimate for each object, or do I have to deal with a series of if trying to apply the correct AgeFactor level for each subject to get the correct fixed effect estimate, after subtracting the contribution of a random effect with Intercept?

The update, apologies, tried to shorten the result that I provided and forgot about str (). Exit:

 >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 and Gender are not used (this is data about men, therefore all Gender values ​​are the same), and NutrientID is similarly fixed for the data.

I make terrible statements made by me that I published, so I'll try your suggestion right away. :)

Update2: this works great with my current data and should be promising for new data, thanks to DWin for additional help with comments for this. :)

 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

It will be something like this (although you really should have given us the results of str (Male.Data), because the model output does not tell us about the coefficients for the base 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") )] 

Basically you run the source data using the match function to select the correct coefficient (s) to add to the interception ... which will be 0 if the data is the base level of the factor (whose spelling I guess at.)

EDIT: I just noticed that you put “-1” in the formula, so maybe all of your AgeFactor terms are listed in the output, and you can scan 0 in the coefficient vector and the invented level of AgeFactor in the match table vector.

Answer 2

The following shows how it has always been easier for me to extract the fixed effect components and the random effects components in the lme4 package. It actually extracts the corresponding correspondence to each observation. Assuming we have a mixed effects model:

 y = Xb + Zu + e 

where Xb are fixed effects and Zu are random effects, we can extract the components (using the lme4 sleepstudy example as an example):

 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 

I know that this works as a generalized approach to extracting components from linear models of mixed effects. For nonlinear models, the model matrix X contains repetitions, and you may need to tweak the code above a bit. Here's the output of validation, as well as visualization using the grid:

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