[R] how to get "lsmeans"?

[Frank E Harrell Jr]

Chuck Cleland wrote:
> 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))

And with the Design package you can do:

f <- ols(y ~ ...)
ds <- gendata(A=c('1','2'),B=c('1','2'),C=c('1','2'))
predict(f, ds)

But this sets the covariate to the median (nothing unreasonable about 
that, just will not agree with SAS).  To set to mean add d=mean(dat$d) 
in gendata.

Frank

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


-- 
Frank E Harrell Jr   Professor and Chair           School of Medicine
                      Department of Biostatistics   Vanderbilt University