[R-sig-eco] R-sig-ecology Digest, Vol 79, Issue 21
Ludovico Frate
ludovicofrate at hotmail.it
Thu Oct 30 12:08:45 CET 2014
Thank you Tim for your suggestion! I'll try with mixed model!
Ludovico
Frate
PhD student (University of Molise - Italy)
Environmetrics Lab
http://www.distat.unimol.it/STAT/environmetrica/organico/collaboratori/ludovico-frate-1
Department of Biosciences and Territory - DiBT
Università del Molise.
Contrada Fonte
Lappone,
86090 - Pesche (IS)
ITALIA.
Cel: ++39
3333767557
Fax: ++39 (0874) 404123
E-mail ludovico.frate at unimol.it
ludovicofrate at hotmail.it
> From: r-sig-ecology-request at r-project.org
> Subject: R-sig-ecology Digest, Vol 79, Issue 21
> To: r-sig-ecology at r-project.org
> Date: Thu, 30 Oct 2014 12:00:01 +0100
>
> Send R-sig-ecology mailing list submissions to
> r-sig-ecology at r-project.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
> or, via email, send a message with subject or body 'help' to
> r-sig-ecology-request at r-project.org
>
> You can reach the person managing the list at
> r-sig-ecology-owner at r-project.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of R-sig-ecology digest..."
>
>
> Today's Topics:
>
> 1. Re: plyr and mvabund, conceptual issue (Hadley Wickham)
> 2. Re: plyr and mvabund, conceptual issue (Hadley Wickham)
> 3. draw 95% confidence limits for non significant regression
> lines (Paolo Piras)
> 4. species richness, GLM and negative values (Ludovico Frate)
> 5. Re: plyr and mvabund, conceptual issue (Maas, Kendra)
> 6. Re: plyr and mvabund, conceptual issue (Krzysztof Sakrejda)
> 7. Re: species richness, GLM and negative values (Martin Weiser)
> 8. Re: species richness, GLM and negative values (Bob O'Hara)
> 9. Re: plyr and mvabund, conceptual issue (Maas, Kendra)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Wed, 29 Oct 2014 08:38:46 -0500
> From: Hadley Wickham <h.wickham at gmail.com>
> To: Eduard Sz?cs <eduardszoecs at gmail.com>
> Cc: "r-sig-ecology at r-project.org" <r-sig-ecology at r-project.org>
> Subject: Re: [R-sig-eco] plyr and mvabund, conceptual issue
> Message-ID:
> <CABdHhvHcFL_mVCL-Smi5_2mUK-a-s=o5sk1rvdk_B6+mqLqPHg at mail.gmail.com>
> Content-Type: text/plain; charset=UTF-8
>
> On Wed, Oct 29, 2014 at 4:23 AM, Eduard Sz?cs <eduardszoecs at gmail.com> wrote:
> > Hai Kendra,
> >
> > i've used a simple for-loop to do this in the past.
> >
> > Something along these lines:
> >
> >
> > ###-----------------------------------------------------------------
> > mymv <- function(response, env, zone) {
> > df <- data.frame(env = env, zone = zone)
> > out <- NULL
> > for (i in levels(zone)) {
> > rsp <- mvabund(response[zone == i, ])
> > out[[i]]$mod <- manyglm(rsp ~ env, data = df[zone == i, ])
> > out[[i]]$anova <- anova(out[[i]]$mod, p.uni = "adjusted", resamp =
> > "perm.resid", nBoot = 10)
> > }
> > return(out)
> > }
>
> If you're going to use a for-loop, you need to preallocate the output:
>
> out <- vector("list", length(levels(zone))
>
> Otherwise each iteration of the loop needs to copy all previous output
> to a new location, making things rather slow.
>
> Hadley
>
> --
> http://had.co.nz/
>
>
>
> ------------------------------
>
> Message: 2
> Date: Wed, 29 Oct 2014 08:41:54 -0500
> From: Hadley Wickham <h.wickham at gmail.com>
> To: "Maas, Kendra" <kendrami at mail.ubc.ca>
> Cc: "r-sig-ecology at r-project.org" <r-sig-ecology at r-project.org>
> Subject: Re: [R-sig-eco] plyr and mvabund, conceptual issue
> Message-ID:
> <CABdHhvGDzaPMHeCU3pjXq4XSW4Htt2OZb0nq-sCr-7ZZfoBccg at mail.gmail.com>
> Content-Type: text/plain; charset=UTF-8
>
> 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/
>
>
>
> ------------------------------
>
> Message: 3
> Date: Wed, 29 Oct 2014 13:49:14 +0000
> From: Paolo Piras <paolo.piras at uniroma3.it>
> To: "r-sig-ecology at r-project.org" <r-sig-ecology at r-project.org>
> Subject: [R-sig-eco] draw 95% confidence limits for non significant
> regression lines
> Message-ID: <1414590683967.64008 at uniroma3.it>
> Content-Type: text/plain; charset="iso-8859-1"
>
> Dear list,
> I write for a surely trivial doubt about plotting 95% CI of regression lines.
> I saw that someone draw them even for non significant regression lines. I think that this is not correct
> as , in the majority of cases, Ho is that beta is not significantly different from 0 and consequently the observed beta is just due to chance.
> On the other hand one could argue that 95% CI limits are still informative about the distribution of the data.
>
> Thanks in advance for an opinion.
> regards
> paolo
>
>
> ------------------------------
>
> Message: 4
> Date: Wed, 29 Oct 2014 16:27:25 +0100
> From: Ludovico Frate <ludovicofrate at hotmail.it>
> To: "r-sig-ecology at r-project.org" <r-sig-ecology at r-project.org>
> Subject: [R-sig-eco] species richness, GLM and negative values
> Message-ID: <DUB126-W62B34E6E01FDB7291648FDD69C0 at phx.gbl>
> Content-Type: text/plain; charset="UTF-8"
>
> Dear all,I'am trying to fit a very simple linear model. I am analyzing the differences in the number of species (DS) found in several permanent plots in two year of observations.
> Firstly, I have calculated the differences per plot (i.e. number of species in Plot 1 in Time A - number of species in Plot 1 in Time B and so on for all the plots).Secondly, those differences were tested for deviation from zero by means of a linear model
> M2<-lm(DC~1, data = gransasso)summary(M2)E2<-residuals(M2)qqnorm(E2, pch = 19, col = "blue"); qqline(E2, col = "red")
> The qqnorm has shown that residuals were not normally distributed, thus I need to use a GLM. However GLM (poisson family) does not work with negative values (DS has negative values).I've tried to add a constant value to these differences (i.e. +100) but the result is misleading since I am testing for deviation from zero.
> Do you have any suggestions?
> Regards,Ludovico
>
>
> Ludovico
> Frate
>
> PhD student (University of Molise - Italy)
> Environmetrics Lab
> http://www.distat.unimol.it/STAT/environmetrica/organico/collaboratori/ludovico-frate-1
> Department of Biosciences and Territory - DiBT
> Universit?el Molise.
> Contrada Fonte
> Lappone,
> 86090 - Pesche (IS)
> ITALIA.
> Cel: ++39
> 3333767557
> Fax: ++39 (0874) 404123
> E-mail ludovico.frate at unimol.it
> ludovicofrate at hotmail.it
>
> [[alternative HTML version deleted]]
>
>
>
> ------------------------------
>
> Message: 5
> Date: Wed, 29 Oct 2014 17:09:22 +0000
> From: "Maas, Kendra" <kendrami at mail.ubc.ca>
> To: Hadley Wickham <h.wickham at gmail.com>
> Cc: "r-sig-ecology at r-project.org" <r-sig-ecology at r-project.org>
> Subject: Re: [R-sig-eco] plyr and mvabund, conceptual issue
> Message-ID:
> <CC055ABE91B74741AB0ADA913A0061AF019205BB90 at s-itsv-mbx03p.ead.ubc.ca>
> Content-Type: text/plain; charset="us-ascii"
>
> 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/
>
>
>
> ------------------------------
>
> Message: 6
> Date: Wed, 29 Oct 2014 13:27:12 -0400
> From: Krzysztof Sakrejda <krzysztof.sakrejda at gmail.com>
> To: "Maas, Kendra" <kendrami at mail.ubc.ca>
> Cc: "r-sig-ecology at r-project.org" <r-sig-ecology at r-project.org>
> Subject: Re: [R-sig-eco] plyr and mvabund, conceptual issue
> Message-ID:
> <CAJCSVaDxa+JySR_ER8kXJc=4ssY9tFUAUV59RuTFi+w3haiD_g at mail.gmail.com>
> Content-Type: text/plain; charset=UTF-8
>
> 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
>
>
>
> ------------------------------
>
> Message: 7
> Date: Wed, 29 Oct 2014 20:32:18 +0100
> From: Martin Weiser <weiser2 at natur.cuni.cz>
> To: Ludovico Frate <ludovicofrate at hotmail.it>
> Cc: "r-sig-ecology at r-project.org" <r-sig-ecology at r-project.org>
> Subject: Re: [R-sig-eco] species richness, GLM and negative values
> Message-ID: <1414611138.3848.39.camel at zouzel.natur.cuni.cz>
> Content-Type: text/plain; charset=UTF-8
>
> Ludovico Frate p??e v St 29. 10. 2014 v 16:27 +0100:
> > Dear all,I'am trying to fit a very simple linear model. I am analyzing the differences in the number of species (DS) found in several permanent plots in two year of observations.
> > Firstly, I have calculated the differences per plot (i.e. number of species in Plot 1 in Time A - number of species in Plot 1 in Time B and so on for all the plots).Secondly, those differences were tested for deviation from zero by means of a linear model
> > M2<-lm(DC~1, data = gransasso)summary(M2)E2<-residuals(M2)qqnorm(E2, pch = 19, col = "blue"); qqline(E2, col = "red")
> > The qqnorm has shown that residuals were not normally distributed, thus I need to use a GLM. However GLM (poisson family) does not work with negative values (DS has negative values).I've tried to add a constant value to these differences (i.e. +100) but the result is misleading since I am testing for deviation from zero.
> > Do you have any suggestions?
> > Regards,Ludovico
> >
> >
> > Ludovico
> > Frate
> >
> > PhD student (University of Molise - Italy)
> > Environmetrics Lab
> > http://www.distat.unimol.it/STAT/environmetrica/organico/collaboratori/ludovico-frate-1
> > Department of Biosciences and Territory - DiBT
> > Universit del Molise.
> > Contrada Fonte
> > Lappone,
> > 86090 - Pesche (IS)
> > ITALIA.
> > Cel: ++39
> > 3333767557
> > Fax: ++39 (0874) 404123
> > E-mail ludovico.frate at unimol.it
> > ludovicofrate at hotmail.it
> >
> > [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-ecology mailing list
> > R-sig-ecology at r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>
> Hi Ludovico,
>
> correct me if I am wrong, but I think you want to say something like:
> "In X out of Y plots, number of species did substantially (statistically
> significantly?) changed".
> For this, I would separate the plots where number of species is bigger
> and those where opposite is true in the second census and perform
> binomial test:
> binom.test(x=abs(difference), n=more.species.i.e.bigger.of.the.numbers,
> p=0.5) per each plot. Maybe you need better p (something like mean
> probability of species to disappear,...)
>
> HTH.
> Martin
>
>
>
>
>
>
> --
>
> ------------------------------
> Pokud je tento e-mail sou??st? obchodn?ho jedn?n?, P??rodov?deck? fakulta
> Univerzity Karlovy v Praze:
> a) si vyhrazuje pr?vo jedn?n? kdykoliv ukon?it a to i bez uveden? d?vodu,
> b) stanovuje, ?e smlouva mus? m?t p?semnou formu,
> c) vylu?uje p?ijet? nab?dky s dodatkem ?i odchylkou,
> d) stanovuje, ?e smlouva je uzav?ena teprve v?slovn?m dosa?en?m shody na
> v?ech n?le?itostech smlouvy.
>
>
>
> ------------------------------
>
> Message: 8
> Date: Wed, 29 Oct 2014 23:01:55 +0100
> From: "Bob O'Hara" <bohara at senckenberg.de>
> To: r-sig-ecology at r-project.org
> Subject: Re: [R-sig-eco] species richness, GLM and negative values
> Message-ID: <545163D3.2020809 at senckenberg.de>
> Content-Type: text/plain; charset="UTF-8"
>
> On 10/29/2014 04:27 PM, Ludovico Frate wrote:
> > Dear all,I'am trying to fit a very simple linear model. I am analyzing the differences in the number of species (DS) found in several permanent plots in two year of observations.
> > Firstly, I have calculated the differences per plot (i.e. number of species in Plot 1 in Time A - number of species in Plot 1 in Time B and so on for all the plots).Secondly, those differences were tested for deviation from zero by means of a linear model
> > M2<-lm(DC~1, data = gransasso)summary(M2)E2<-residuals(M2)qqnorm(E2, pch = 19, col = "blue"); qqline(E2, col = "red")
> > The qqnorm has shown that residuals were not normally distributed, thus I need to use a GLM. However GLM (poisson family) does not work with negative values (DS has negative values).I've tried to add a constant value to these differences (i.e. +100) but the result is misleading since I am testing for deviation from zero.
> > Do you have any suggestions?
> > Regards,Ludovico
> Use the Poisson to model the number of species in each sample, so use
> data like this:
>
> Plot Time DS
> 1 A 4
> 1 B 7
> 2 A 0
> 2 B 1
> 3 A 23
> 3 B 7
> ...
>
> Then you fit the model
>
> Mgood <- glm(DS ~ Plot + Time, family=poisson())
>
> where Plot and Time are factors (use Plot <- factor(Plot), for example,
> if you need to).
>
> You're interested in the Time effect, which is the average difference
> between the numbers of species in the plots in the different times. The
> Plot effect controls for different plots having different numbers of
> species overall. If you look at summary(Mgood), the Time effect is the
> log of the ratio of the species richnesses in times A and B. It will be
> written as TimeB, which means it's log(E.B/E.A) (where E.A and E.B are
> the expected species richnesses at times A and B). So, for example, an
> estimate of 0.4 would mean that at Time B there are exp(0.4)=1.49 times
> more species at time B than time A.
>
> Bob
>
>
> >
> > Ludovico
> > Frate
> >
> > PhD student (University of Molise - Italy)
> > Environmetrics Lab
> > http://www.distat.unimol.it/STAT/environmetrica/organico/collaboratori/ludovico-frate-1
> > Department of Biosciences and Territory - DiBT
> > Universit? del Molise.
> > Contrada Fonte
> > Lappone,
> > 86090 - Pesche (IS)
> > ITALIA.
> > Cel: ++39
> > 3333767557
> > Fax: ++39 (0874) 404123
> > E-mail ludovico.frate at unimol.it
> > ludovicofrate at hotmail.it
> >
> > [[alternative HTML version deleted]]
> >
> >
> >
> > _______________________________________________
> > R-sig-ecology mailing list
> > R-sig-ecology at r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>
>
> --
> Bob O'Hara
>
> Biodiversity and Climate Research Centre
> Senckenberganlage 25
> D-60325 Frankfurt am Main,
> Germany
>
> Tel: +49 69 7542 1863
> Mobile: +49 1515 888 5440
> WWW: http://www.bik-f.de/root/index.php?page_id=219
> Blog: http://blogs.nature.com/boboh
> Journal of Negative Results - EEB: www.jnr-eeb.org
>
>
> [[alternative HTML version deleted]]
>
>
>
> ------------------------------
>
> Message: 9
> Date: Thu, 30 Oct 2014 00:15:22 +0000
> From: "Maas, Kendra" <kendrami at mail.ubc.ca>
> To: Krzysztof Sakrejda <krzysztof.sakrejda at gmail.com>
> Cc: "r-sig-ecology at r-project.org" <r-sig-ecology at r-project.org>
> Subject: Re: [R-sig-eco] plyr and mvabund, conceptual issue
> Message-ID:
> <CC055ABE91B74741AB0ADA913A0061AF019205BC7F at s-itsv-mbx03p.ead.ubc.ca>
> Content-Type: text/plain; charset="us-ascii"
>
> 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
>
>
> ------------------------------
>
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>
>
> End of R-sig-ecology Digest, Vol 79, Issue 21
> *********************************************
[[alternative HTML version deleted]]
More information about the R-sig-ecology
mailing list