[R] eval-parse and lme in a loop

[Greg Snow]

> fortune(106)

If the answer is parse() you should usually rethink the question.
   -- Thomas Lumley
      R-help (February 2005)

Here is one approach that seems to be working:

for (meanCol in paste("Mean_", 1:3, sep="")) {

	cat('\n######',meanCol,'\n\n')

	f <- formula( paste( meanCol, '~ Group + c1 + c2 + c3' ) )
	
	means.lmeWithCovariate <- lme(f, data=df, random = ~ 1 | Subject)
  	print(summary(means.lmeWithCovariate))
  	print(anova(means.lmeWithCovariate))
}



-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at imail.org
801.408.8111


> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
> project.org] On Behalf Of Connolly, Colm
> Sent: Thursday, August 05, 2010 3:18 PM
> To: r-help at r-project.org
> Subject: [R] eval-parse and lme in a loop
> 
> Hi everybody,
> 
> I'm having trouble getting an eval-parse combination to work with lme
> in a for loop.
> 
> Consider the code below. The call to lme outside the loop works. The
> call to aov inside the loop works, but the call to to lme inside the
> loop does not.
> The latter fails with
> 
> Error in model.frame.default(formula = ~meanCol + Group + c1 + c3 +
> Subject,  :
>   variable lengths differ (found for 'Group')
> 
> and I don't know why the eval-parse portion of the code is node being
> replaced as it is in the call to aov. Anybody got any idea what's
> wrong?
> 
> rm(list=ls())
> 
> library('nlme')
> 
> df=data.frame(
>   Group=rep(c("ctrl", "short", "long"), each=9, times=1),
>   Subject=paste("S", 1:27, sep=""),
>   Mean_1=c(-2.335482,  8.269832,  6.704734, -2.355664,  0.724432,
> 11.212987,  5.966652, 10.741065,  3.567522,  5.359832,  4.899951,
> 4.084452,  5.976326, 8.634508,  5.750942,  7.455622,  6.823843,
> 0.903801, 13.198157,  5.238985,  5.171453, 10.273715, 5.184785,
> 6.225874,  8.889997, 12.096968, 9.843028),
>   Mean_2=c(16.244444, -1.486571,  1.817303,  3.814328, -0.343775,
> 19.323227, -0.390764,  7.149855,  4.517766, -0.967120, 11.219844,
> 0.561103,  5.188478, -4.819444,  6.253271,  7.062638,  4.502291, -
> 0.469567, 12.739646,  5.832679,  9.193511,  7.664636, 6.083930,
> 9.493136,  7.478424,  7.642287, 7.205822),
>   Mean_3=c(7.213662,  2.147219, -1.345795, -0.495567,  3.733844,
> 3.544084,  0.949573,  4.318238,  2.914248,  3.913391,  9.240128,
> 8.315122, 12.323612, 14.622557, 11.169847, 10.857631, 14.637243,
> 7.203096,  7.777472, -1.457673,  2.766380,  9.675148, -5.047292,
> 1.982447,  3.865502,  5.045913, 10.660579),
>   c1=rnorm(9, mean=0.1, sd=0.4),
>   c2=rnorm(9, mean=0.5, sd=0.7),
>   c3=rnorm(9, mean=0.9, sd=1)
>   );
> 
> 
> ## this works
> means.lmeNoCovariate=lme(Mean_1 ~ Group, data=df, random = ~ 1 |
> Subject)
> print(summary(means.lmeNoCovariate))
> print(anova(means.lmeNoCovariate))
> 
> 
> for (meanCol in paste("Mean_", 1:3, sep="")) {
>   ## LMEs go here
> 
> 
> cat("##################################################################
> ##################################\n")
>   cat("### AOV\n")
> 
> cat("##################################################################
> ##################################\n")
>   ## this works
>   means.aov<-aov(eval(parse(text=meanCol)) ~ Group, data=df);
>   means.aov.summary=summary(means.aov)
>   print(means.aov.summary)
> 
> 
> cat("##################################################################
> ##################################\n")
>   cat("### LME\n")
> 
> cat("##################################################################
> ##################################\n")
> 
>   ## this doe not work
>   means.lmeWithCovariate=lme(eval(parse(text=meanCol)) ~ Group + c1 +
> c1 + c3, data=df, random = ~ 1 | Subject)
>   print(summary(means.lmeWithCovariate))
>   print(anova(means.lmeWithCovariate))
> 
> }
> 
> Regards,
> --
> Colm G. Connolly, Ph. D.
> Dept of Psychiatry
> University of California, San Diego
> 9500 Gilman Dr. #0738
> La Jolla, CA 92093-0738
> Tel: +1-858-246-0684
> 
> ______________________________________________
> R-help at r-project.org 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.