[R] how to get "lsmeans"?

[Chuck Cleland]

Liaw, Andy wrote:
> I verified the result from the following with output from JMP 6 on the
> same data (don't have SAS: don't need it):
> 
> set.seed(631)
> n <- 100
> dat <- data.frame(y=rnorm(n), A=factor(sample(1:2, n, replace=TRUE)),
>                   B=factor(sample(1:2, n, replace=TRUE)),
>                   C=factor(sample(1:2, n, replace=TRUE)),
>                   d=rnorm(n))
> fm <- lm(y ~ A + B + C + d, dat)
> ## Form a data frame of points to predict: all combinations of the
> ## three factors and the mean of the covariate.
> p <- data.frame(expand.grid(A=1:2, B=1:2, C=1:2))
> p[] <- lapply(p, factor)
> p <- cbind(p, d=mean(dat$d))
> p <- cbind(yhat=predict(fm, p), p)
> ## lsmeans for the three factors:
> with(p, tapply(yhat, A, mean))
> with(p, tapply(yhat, B, mean))
> with(p, tapply(yhat, C, mean))

  Using Andy's example data, these are the LSMEANS and intervals I get
from SAS:

A        y LSMEAN      95% Confidence Limits
1       -0.071847       -0.387507     0.243813
2       -0.029621       -0.342358     0.283117

B        y LSMEAN      95% Confidence Limits
1       -0.104859       -0.397935     0.188216
2        0.003391       -0.333476     0.340258

C        y LSMEAN      95% Confidence Limits
1       -0.084679       -0.392343     0.222986
2       -0.016789       -0.336374     0.302795

  One way of reproducing the LSMEANS and intervals from SAS using
predict() seems to be the following:

> dat.lm <- lm(y ~ A + as.numeric(B) + as.numeric(C) + d, data = dat)
> newdat <- expand.grid(A=factor(c(1,2)),B=1.5,C=1.5,d=mean(dat$d))
> cbind(newdat, predict(dat.lm, newdat, interval="confidence"))
  A   B   C          d         fit        lwr       upr
1 1 1.5 1.5 0.09838595 -0.07184709 -0.3875070 0.2438128
2 2 1.5 1.5 0.09838595 -0.02962086 -0.3423582 0.2831165

  However, another possibility seems to be:

> dat.lm <- lm(y ~ A + as.numeric(B) + as.numeric(C) + d, data = dat)
> newdat <-
expand.grid(A=factor(c(1,2)),B=mean(as.numeric(dat$B)),C=mean(as.numeric(dat$C)),d=mean(dat$d))
> cbind(newdat, predict(dat.lm, newdat, interval="confidence"))
  A    B    C          d         fit        lwr       upr
1 1 1.43 1.48 0.09838595 -0.08078243 -0.3964661 0.2349012
2 2 1.43 1.48 0.09838595 -0.03855619 -0.3479589 0.2708465

  The predictions directly above match what effect() in the effects
package by John Fox returns:

library(effects)

> effect("A", fm, xlevels=list(d = mean(dat$D)))

 A effect
A
          1           2
-0.08078243 -0.03855619

  But for some reason the predict() and effect() intervals are a little
different:

> effect("A", fm, xlevels=list(d = mean(dat$D)))$lower
          [,1]
101 -0.3924451
102 -0.3440179

> effect("A", fm, xlevels=list(d = mean(dat$D)))$upper
         [,1]
101 0.2308802
102 0.2669055

  I would be interested in any comments on these different approaches
and on the difference in intervals returned by predict() and effect().

hope this helps,

Chuck

> Andy 
> 
> From: Xingwang Ye
>> Dear all, 
>>       
>>     I search the mail list about this topic and learn that no 
>> simple way is available to get "lsmeans" in R as in SAS.
>>     Dr.John Fox and Dr.Frank E Harrell have given very useful 
>> information about "lsmeans" topic.    
>>     Dr. Frank E Harrell suggests not to think about lsmeans, 
>> just to think about what predicted values wanted
>>     and to use the predict function. However, after reading 
>> the R help file for a whole day, I am still unclear how to do it.
>>     Could some one give me a hand? 
>>  
>> for example:
>>   
>> A,B and C are binomial variables(factors); d is a continuous 
>> variable ; The response variable Y is  a continuous variable too.  
>>
>> To get lsmeans of Y according to A,B and C, respectively, in 
>> SAS, I tried proc glm data=a;  class A B C;  model Y=A B C d; 
>>  lsmeans A B C/cl; run;  
>>
>> In R, I tried this:  
>>  library(Design)
>>  ddist<-datadist(a)
>>  options(datadist="ddist")
>>  f<-ols(Y~A+B+C+D,data=a,x=TRUE,y=TRUE,se.fit=TRUE)  
>>
>> then how to get the "lsmeans" for A, B, and C, respectively 
>> with predict function?
>>
>>  
>>
>> Best wishes
>> yours, sincerely 
>> Xingwang Ye    
>> PhD candidate     
>> Research Group of Nutrition Related Cancers and Other Chronic 
>> Diseases      
>> Institute for Nutritional Sciences,  
>> Shanghai Institutes of Biological Sciences,     
>> Chinese Academy of Sciences     
>> P.O.Box 32     
>> 294 Taiyuan Road     
>> Shanghai 200031     
>> P.R.CHINA
>>
>>
> 
> 
> ------------------------------------------------------------------------------
> Notice:  This e-mail message, together with any attachments,...{{dropped}}
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Chuck Cleland, Ph.D.
NDRI, Inc.
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894