[R-sig-eco] plyr and mvabund, conceptual issue

[Eduard Szöcs]

Forgotten to mention, that mymv() returns a list with two components
(the model and the anova).

You can then extract the information you need from this list, maybe like
this:

###-----------------------------------------------
per_zone <- mymv(response, env, zone)
# p-values from univariate GLMs
sapply(per_zone, function(y) y$anova$uni.p[2, ])
# coefs for env
sapply(per_zone, function(y) y$mod$coef[2, ])
###-----------------------------------------------


Cheers,

Eduard




On 29/10/14 00:56, Maas, Kendra wrote:
> Forgot the output that I need from these, so using my example for my data
> 
>> 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")
> 
> 
>> write.csv(b.oJP.nb.anova$uni.p, file="b.oJP.nb.anova")
>> write.csv(b.oJP.nb$coef, file="b.oJP.nb.coef")
> 
> 
> 
> 
> 
> ________________________________________
> From: r-sig-ecology-bounces at r-project.org [r-sig-ecology-bounces at r-project.org] on behalf of Maas, Kendra [kendrami at mail.ubc.ca]
> Sent: Tuesday, October 28, 2014 4:31 PM
> To: r-sig-ecology at r-project.org
> Subject: [R-sig-eco] 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
> 
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>