plyr and mvabund, conceptual issue

I'm trying to run mvabund (generate glms for each species and do univariate anova to determine "indicator species" that respond to my treatments) on a lot of subsets of my data.  I'm having theoretical difficulty with how to use plyr on multiple dataframes or lists and outputting lists.  Previously I've run this series of commands using text editor to change the selected zone, I know that this is what plyr is designed for but I'm getting stuck

"b.red.otu" is a sample by species dataframe, "b.env" is a sample by factor dataframe containing the variables "zone" and "om"

> b.oJP.mva<-mvabund(subset(b.red.otu,zone=="oJP"))
> b.oJP.nb<-manyglm(b.oJP.mva~subset(b.env, b.env$zone=="oJP")$om, family="negative.binomial")
> b.oJP.nb.anova<-anova(b.oJP.nb,p.uni="adjusted", resamp="perm.resid")


This code works, it's just really ugly and requires a lot of copy and paste/find and replace for every possible zone (I have tens of subsets that I want to look at)


Here is how I'm attempting to work out my code with the Tasmania data packaged with mvabund.  I convert it to dataframes because I much more comfortable with them.


> tas.env <- data.frame(Tasmania$treatment, Tasmania$block)

> tas.abund <- data.frame(Tasmania$abund)

> mva.out <- dlply(tas.abund, ~tas.env$block, function(x) {
      mvabund(x~tas.env$treatment)
       })

Which returns an empty list[0]

I tried creating vectors for treatment and block

> block <- Tasmania$block
> treatment <- Tasmania$treatment

> mva.out <- dlply(tas.abund, ~block, function(x) {
     mvabund(x~treatment)
     })

Error in as.data.frame.default(x[[i]], optional = TRUE) :
  cannot coerce class ""formula"" to a data.frame

I tried llply since mvabund puts out a list and Tasmania is already a list

> mva.out <- llply(Tasmania$abund, ~Tasmania$block, function(x) {
      mvabund(x~Tasmania$treatment)
      })

Error in llply(Tasmania$abund, ~Tasmania$block, function(x) { :
  .fun is not a function.


I'm sure this is possible to do with plyr, I just can't figure out how.  Suggestions please.

thanks

Kendra
_______________________________________________
R-sig-ecology mailing list
R-sig-ecology at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology