extract fixed effect for first level in lme function nlme package

Question

I would like estimate Best Linear Unbiased Estimates. I have data like this.

set.seed(1234)
data1 <- data.frame (entry=rep(1:20, 3), repl = rep(1:3, each=20), 
yld = rnorm(60)+50)
require(nlme)
data1$entry <- as.factor (data1$entry)
data1$repl <- as.factor (data1$repl)

# with intercept 
fm1 <- lme(yld ~  entry, random = ~ 1|repl, data=data1 ) 
fixed.effects(fm1)

(Intercept)      entry2      entry3      entry4      entry5      entry6      entry7      entry8      entry9     entry10 
50.12550625 -0.55280604 -0.19599669 -0.84775000 -0.54515146 -0.76239402 -0.49460706 -1.06626407 -0.49331238 -0.89978504 
    entry11     entry12     entry13     entry14     entry15     entry16     entry17     entry18     entry19     entry20 
-0.51914838 -0.81085806 -0.99036743 -0.60942666 -0.40280922 -0.36378922 -0.47325014 -1.13402026  0.03264178  0.13853678 

You can see I did not get fixed effect of first level of the variable. If I remove the intercept term it gives me fixed effect for first level of the entry variable.

# without intercept 
fm2 <- lme(yld ~  -1 + entry, random = ~ 1|repl, data=data1 ) 
fixed.effects(fm2)

entry1   entry2   entry3   entry4   entry5   entry6   entry7   entry8   entry9  entry10  entry11  entry12  entry13 
50.12551 49.57270 49.92951 49.27776 49.58035 49.36311 49.63090 49.05924 49.63219 49.22572 49.60636 49.31465 49.13514 
 entry14  entry15  entry16  entry17  entry18  entry19  entry20 
49.51608 49.72270 49.76172 49.65226 48.99149 50.15815 50.26404 

But I do not think my case is intercept is 0 case. So I want to fit the first model with intercept still get the fixed effects. Am I missing something little?

Answer 1

Generally, I wouldn't drop the intercept term if there isn't a reason for it (see also these answers. If your problem is just getting the same values for the fixed effects as seen in fm2, simply subtract or add the entry2, entry3, $...$ , entry20 from the intercept value (which in this case is the estimate for entry1), e.g.:

entry2: $50.12550625-0.55280604 = 49.5727$

entry20: $50.12550625+0.13853678 = 50.26404$

In a regression output with categorical predictors, the intercept is the expected mean value of $Y$ when all $X=0$.

But as @Roland pointed out, better use the lsmeans() function in the lsmeans package:

> library(lsmeans)
> lsmeans(fm1,"entry")
 entry   lsmean        SE df lower.CL upper.CL
 1     50.12551 0.6053462  2 47.52091 52.73010
 2     49.57270 0.6053462  2 46.96811 52.17729
 3     49.92951 0.6053462  2 47.32491 52.53410
 4     49.27776 0.6053462  2 46.67316 51.88235
 5     49.58035 0.6053462  2 46.97576 52.18495
 6     49.36311 0.6053462  2 46.75852 51.96771
 7     49.63090 0.6053462  2 47.02630 52.23549
 8     49.05924 0.6053462  2 46.45465 51.66384
 9     49.63219 0.6053462  2 47.02760 52.23679
 10    49.22572 0.6053462  2 46.62113 51.83032
 11    49.60636 0.6053462  2 47.00176 52.21095
 12    49.31465 0.6053462  2 46.71005 51.91924
 13    49.13514 0.6053462  2 46.53054 51.73973
 14    49.51608 0.6053462  2 46.91148 52.12067
 15    49.72270 0.6053462  2 47.11810 52.32729
 16    49.76172 0.6053462  2 47.15712 52.36631
 17    49.65226 0.6053462  2 47.04766 52.25685
 18    48.99149 0.6053462  2 46.38689 51.59608
 19    50.15815 0.6053462  2 47.55355 52.76274
 20    50.26404 0.6053462  2 47.65945 52.86864

Confidence level used: 0.95