[R] Passing data to aov

David Winsemius dwinsemius at comcast.net
Thu Aug 26 16:14:01 CEST 2010


On Aug 26, 2010, at 8:47 AM, Bryan Hanson wrote:

> Hello Again Gurus and Lurkers:
>
> I’m trying to build a very user-friendly function which does aov  
> without
> having the user type in a formula (which would be tedious in this  
> case).
> The idea is to take the response from a PCA score matrix, and the  
> factors
> from a list.  A simple example is the function given below, along  
> with test
> data and a sample call to the function.
>
> I'm certainly having trouble understanding the proper ways to work  
> with
> formulas and related items, but I think what I do in the function  
> should
> work (it's built on ideas dug out of the archives).  However, when  
> the data
> is passed to aov (directly or via manova), something in the bowels  
> of aov
> complains with the following error:
>
> Error in model.frame.default(formula = form, drop.unused.levels =  
> TRUE) :
>  object is not a matrix

Actually I got a different error (also on a Mac. but with a version  
that is two month later than yours and with quite a smaller number of  
packages loaded):
 > hypTestScores(mylist = td2, score.matrix = td1,
+    fac = c("f1", "f2"))
Error in eval(expr, envir, enclos) : object 'scores' not found

Which went away if I removed the "env=" argument. My guess is that you  
were telling manova to look further up the lexical tree than it should  
have (or that attach() operation messed up the environment somehow):

 > hypTestScores(mylist = td2, score.matrix = td1,
+    fac = c("f1", "f2"))
           Df  Pillai approx F num Df den Df Pr(>F)
f1         1 0.31663  0.61777      3      4 0.6392
f2         1 0.36652  0.77145      3      4 0.5673
f1:f2      1 0.15687  0.24808      3      4 0.8593
Residuals  6

 > detach("mylist") # needed if there is an error
Error in detach("mylist") : invalid 'name' argument

You might want to do inside the function something like

res =try(  your function )
if (res="try-error") {stop() } else{
     <process>}

>
> To me, the formula looks legitimate, and the variables in the  
> formula are
> all in the environment (I think: The way I am doing this is  
> basically that
> described in ?summary.manova where only a formula is passed, no data
> argument).  Based upon reading the archives, the problem might arise  
> with
> one of the deparse statements in aov, but I can't resolve it.  It  
> might also
> be one of scoping/environment, but again, this is only an idea.
>
> TIA for any assistance.  Bryan
> *************
> Bryan Hanson
> Professor of Chemistry & Biochemistry
> DePauw University, Greencastle IN USA
>
>
> hypTestScores <-
>    function(mylist, score.matrix, pcs = 1:3, fac = NULL, ...) {
>
>    scores <- score.matrix[,pcs]
> #    str(scores)    # looks correct to me
>    form <- as.formula(paste("scores", "~", paste(fac, collapse =  
> "*")),
>        env = parent.frame())
> #    str(form)       # looks correct to me
>    attach(mylist)
>    if (length(pcs) > 1) out <- manova(formula = form, ...)
>    if (length(pcs) == 1) out <- aov(formula = form, ...)
>    print(summary(out))
>    detach(mylist)
>    invisible(out)
>    }
>
> # test data
> td1 <- matrix(rnorm(50), ncol = 5) # like PCA scores
> td2 <- list(
>    f1 = as.factor(sample(c("A", "B"), 10, replace = TRUE)),
>    f2 = as.factor(sample(c("C", "D"), 10, replace = TRUE)))
>
> # test call
> hypTestScores(mylist = td2, score.matrix = td1,
>    fac = c("f1", "f2"))
> detach("mylist") # needed if there is an error
>
>> sessionInfo()
> R version 2.11.0 (2010-04-22)
> x86_64-apple-darwin9.8.0
>
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8
>
> attached base packages:
> [1] datasets  tools     grid      graphics  grDevices utils     stats
> [8] methods   base
>
> other attached packages:
> [1] faraway_1.0.4      GGally_0.2         xtable_1.5-6
> [4] mvbutils_2.5.1     ggplot2_0.8.8      proto_0.3-8
> [7] reshape_0.8.3      ChemoSpec_1.45     R.utils_1.4.0
> [10] R.oo_1.7.2         R.methodsS3_1.2.0  rgl_0.91
> [13] lattice_0.18-5     mvoutlier_1.4      plyr_1.0.3
> [16] RColorBrewer_1.0-2 chemometrics_0.8   som_0.3-5
> [19] robustbase_0.5-0-1 rpart_3.1-46       pls_2.1-0
> [22] pcaPP_1.8-1        mvtnorm_0.9-9      nnet_7.3-1
> [25] mclust_3.4.4       MASS_7.3-5         lars_0.9-7
> [28] e1071_1.5-23       class_7.3-2
>
> ______________________________________________
> 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.

David Winsemius, MD
West Hartford, CT



More information about the R-help mailing list