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

Thu Oct 30 01:15:22 CET 2014

ok, thanks.  So properly written out, I want to split each block out of the dataset and run mvabund, manyglm, and anova on each subset, pulling coefficients out of the manyglm object and p values out of the anova.  I can split the species matrix and env matrix by hand within the commands, so for block one it would look like this.  My question was how to use plyr and run each block.  In the mean time I'll just use Edward's for loop because slower and functioning is still faster than using find/replace in a text editor.

install.packages('plyr')
install.packages('mvabund')

library(mvabund)
library(plyr)

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

#But here is where I'm stuck writing a reproducible example because not knowing how to combine plyr and mvabund is exactly my problem, so for one block I'd write:

tas.1.mva<- mvabund(tas.abund[tas.env$Tasmania.block==1,])
tas.1.glm <- manyglm(tas.1.mva~tas.env[tas.env$Tasmania.block==1, ]$Tasmania.treatment, family="negative.binomial")
tas.1.anova <- anova(tas.1.glm, p.uni="adjusted", resamp="perm.resid")

write.csv(tas.1.anova$uni.p, file="tas.1.p.csv")
write.csv(tas.1.glm$coef, file="tas.1.coef.csv")

--
Kendra Maas, Ph.D.
Post Doctoral Research Fellow
University of British Columbia
604-822-5646

________________________________________
From: Krzysztof Sakrejda [krzysztof.sakrejda at gmail.com]
Sent: Wednesday, October 29, 2014 10:27 AM
To: Maas, Kendra
Cc: r-sig-ecology at r-project.org
Subject: Re: [R-sig-eco] plyr and mvabund, conceptual issue

Since you asked, this is a minimally reproducible example:

install.packages('plyr')
install.packages('mvabund')

library(mvabund)
library(plyr)
data(Tasmania)
tas.abund <- data.frame(Tasmania$abund)
tas.env <- data.frame(Tasmania$treatment, Tasmania$block)
mva.out <- dlply(tas.abund, ~tas.env$block, function(x) {
    mvabund(x~tas.env$treatment)
})

It looks like Hadley didn't quite hit the nail on the head:

mva.out <- dlply(Tasmania, "block", function(df) {
   mvabund(block ~ treatment, data = df)
})

You would typically also state what output you are expecting, or maybe
show a successful mvabund call on a subset since not everyone will be
familiar with mvabund.

Often people who can spot the solution to your problem in 2 seconds
don't want to be bothered by trying to figure out exactly where the
data/packages came from (it isn't always CRAN) or how to call a
package they don't know.

Hope that helps,

Krzysztof

Krzysztof Sakrejda

Division of Biostatistics and Epidemiology / Organismic and Evolutionary Biology
University of Massachusetts, Amherst
319 Morrill Science Center South
611 N. Pleasant Street
Amherst, MA 01003

work #: 413-325-6555
email: sakrejda at umass.edu
-----------------------------------------------

On Wed, Oct 29, 2014 at 1:09 PM, Maas, Kendra <kendrami at mail.ubc.ca> wrote:
> I tried to replicate my problem with the data that is supplied with the mvabund package.  How is that not a minimally reproducible example?  Because I stop after the first step?
>
>
>
> ________________________________________
> From: Hadley Wickham [h.wickham at gmail.com]
> Sent: Wednesday, October 29, 2014 6:41 AM
> To: Maas, Kendra
> Cc: r-sig-ecology at r-project.org
> Subject: Re: [R-sig-eco] plyr and mvabund, conceptual issue
>
> It's really hard to help with out a minimal reproducible example, but
> normally a dlply call would look more like this:
>
> mva.out <- dlply(Tasmania, "block", function(df) {
>   mvabund(block ~ treatment, data = df)
> })
>
> Hadley
>
> On Tue, Oct 28, 2014 at 6:31 PM, Maas, Kendra <kendrami at mail.ubc.ca> wrote:
>> 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
>
>
>
> --
> http://had.co.nz/
>
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology