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@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 ______________________________________________ R-help@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.