Getting Generalized Least Squares Means for fixed effects in nlme or lme4

Question

Least Squares Means with their standard errors for aov object can be obtained with model.tables function:

npk.aov <- aov(yield ~ block + N*P*K, npk)
model.tables(npk.aov, "means", se = TRUE)

I wonder how to get the generalized least squares means with their standard errors from nlme or lme4 objects:

library(nlme)
data(Machines)
fm1Machine <- lme(score ~ Machine, data = Machines, random = ~ 1 | Worker )

Any comment and hint will be highly appreciated. Thanks

Answer 1

lme and nlme fit through maximum likelihood or restricted maximum likelihood (the latter is the default), so your results will be based on either of those methods

summary(fm1Machine) will provide you with the output that includes the means and standard errors:

....irrelevant output deleted
Fixed effects: score ~ Machine 
               Value Std.Error DF  t-value p-value
(Intercept) 52.35556  2.229312 46 23.48507       0
MachineB     7.96667  1.053883 46  7.55935       0
MachineC    13.91667  1.053883 46 13.20514       0
 Correlation: 
....irrelevant output deleted

Because you have fitted the fixed effects with an intercept, you get an intercept term in the fixed effects result instead of a result for MachineA. The results for MachineB and MachineC are contrasts with the intercept, so to get the means for MachineB and MachineC, add the value of each to the intercept mean. But the standard errors are not the ones you would like.

To get the information you are after, fit the model so it doesn't have an intercept term in the fixed effects (see the -1 at the end of the fixed effects:

fm1Machine <- lme(score ~ Machine-1, data = Machines, random = ~ 1 | Worker )

This will then give you the means and standard error output you want:

....irrelevant output deleted
Fixed effects: score ~ Machine - 1 
            Value Std.Error DF  t-value p-value
MachineA 52.35556  2.229312 46 23.48507       0
MachineB 60.32222  2.229312 46 27.05867       0
MachineC 66.27222  2.229312 46 29.72765       0
....irrelevant output deleted

Answer 2

To quote Douglas Bates from

http://markmail.org/message/dqpk6ftztpbzgekm

"I have a strong suspicion that, for most users, the definition of lsmeans is "the numbers that I get from SAS when I use an lsmeans statement". My suggestion for obtaining such numbers is to buy a SAS license and use SAS to fit your models."