Re: [R-sig-eco] Ecological datasets for teaching statistics

[Eduard Szöcs]

> but the p-values generated were
all the same for each pair of comparison. Is this output normal?

With 199 bootstrap samples the smallest p you can get is 1/199=0,00502512...

Hope this helps,

Eduard

Am 19. Juni 2020, 11:33, um 11:33, Thierry Onkelinx  
schrieb:
>Try to send plain text mail without attachments.
>
>ir. Thierry Onkelinx
>Statisticus / Statistician
>
>Vlaamse Overheid / Government of Flanders
>INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE
>AND
>FOREST
>Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
>thierry.onkel...@inbo.be
>Havenlaan 88 bus 73, 1000 Brussel
>www.inbo.be
>
>///
>To call in the statistician after the experiment is done may be no more
>than asking him to perform a post-mortem examination: he may be able to
>say
>what the experiment died of. ~ Sir Ronald Aylmer Fisher
>The plural of anecdote is not data. ~ Roger Brinner
>The combination of some data and an aching desire for an answer does
>not
>ensure that a reasonable answer can be extracted from a given body of
>data.
>~ John Tukey
>///
>
>
>
>
>Op vr 19 jun. 2020 om 11:23 schreef Kanfra, Xorla <
>xorla.kan...@julius-kuehn.de>:
>
>> Dear Roman,
>>
>> Thank you for the suggestion. I tried a few minutes ago but go this
>> response  ‘’The message's content type was not explicitly allowed’’.
>>
>> Please can you advise?
>>
>> Kind regards
>>
>> Xorla
>>
>>
>>
>>
>>
>>
>>
>> *From:* Roman Luštrik [mailto:roman.lust...@gmail.com]
>> *Sent:* Friday, June 19, 2020 8:40 AM
>> *To:* Kanfra, Xorla
>> *Cc:* Thierry Onkelinx; Manuel Spínola; r-sig-ecology@r-project.org
>> *Subject:* Re: [R-sig-eco] Ecological datasets for teaching
>statistics
>>
>>
>>
>> Your question is likely to get lost in the debris. Try posting it as
>a new
>> topic.
>>
>>
>>
>> Cheers,
>>
>> Roman
>>
>>
>>
>> On Fri, Jun 19, 2020 at 7:40 AM Kanfra, Xorla <
>> xorla.kan...@julius-kuehn.de> wrote:
>>
>> Dear All,
>> Sorry to use this thread for communication. I subscribed to this
>forum
>> just a few days ago but couldn’t manage to communicate across the
>board.
>> Please pardon me. If any of you can help me get around a statistical
>> challenge
>> I have just started using manyglm to analyze nematode with associated
>> microbial community data and I am really happy with how robust the
>> algorithm works. I have a question about the adjusted p- values. I
>run a
>> pairwise comparison for one of my factors but the p-values generated
>were
>> all the same for each pair of comparison. Is this output normal?
>besides,
>> are the p-values adjusted when running the pairwise comparison? Which
>> method does it use? ex Tukey, FDR, Holm, BH, etc. I am a bit confused
>as I
>> do not understand really what is going on behind the scripts.
>> Below is the script I used.
>> Thanks and hoping to hear from you soon
>> Xorla Kanfra
>>
>> mod_pairwise <-anova.manyglm(otutable.glm, nBoot=199, cor.type = "I",
>test
>> = "LR", p.uni="adjusted",pairwise.comp = ~envdata$soil)
>>
>> Analysis of Deviance Table
>>
>> Model: manyglm(formula = otutable.fv, family = "negative.binomial")
>>
>> Multivariate test:
>>   Res.Df Df.diff   Dev Pr(>Dev)
>> (Intercept)   88
>> envdata$soil  81   7 429680.005 **
>> envdata$response  80   1  48240.005 **
>> envdata$soil:envdata$response 73   7 160660.005 **
>> ---
>> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>>
>> Pairwise comparison results:
>>Observed statistic Free Stepdown
>> Adjusted P-Value
>> envdata$soil:H vs envdata$soil:HH6945
>> 0.005 **
>> envdata$soil:E vs envdata$soil:HH6881
>> 0.005 **
>> envdata$soil:HH vs envdata$soil:R6471
>> 0.005 **
>> envdata$soil:HH vs envdata$soil:M6458
>> 0.005 **
>> envdata$soil:E vs envdata$soil:EE6301
>> 0.005 **
>> envdata$soil:HH vs envdata$soil:K6032
>> 0.005 **
>> envdata$soil:R vs envdata$soil:RR5993
>> 0.005 **
>> envdata$soil:EE vs envdata$soil:RR   5855
>> 0.005 **
>> envdata$soil:EE vs envdata$soil:M5835
>> 0.005 **
>>
>>
>>
>> -Original Message-
>> From: R-sig-ecology [mailto:r-sig-ecology-boun...@r-project.org] On
>> Behalf Of Thierry Onkelinx
>> Sent: Thursday, June 18, 2020 9:07 PM
>> To: Manuel Spínola
>> Cc: r-sig-ecology@r-project.org
>> Subject: Re: [R-sig-eco] Ecological datasets for teaching statistics
>>
>> Dear Manuel,
>>
>> Our institute has published 54 datasets under an open data licence at
>GBIF
>> (
>>
>>
>https://www.gbif.org/dataset/search?publishing_org=1cd669d0-80ea-11de-a9d0-f1765f95f18b
>> )
>

Re: [R-sig-eco] How to obtain P value from Monte Carlo sampling for adonis (permanova)?

[Eduard Szöcs]

Dear Ellen,

You can create a Permutation matrix holding ALL possible permutations
for a given permutation design using permute::allPerms() [Checkout the
help therein].

This matrix can be passed to adonis/adonis2 via the permutations= argument.

Hope this helps,

Eduard

On 18.11.2016 09:39, Ellen Pape wrote:
> Dear all,
> 
> I have recently decided to switch from Permanova/Primer to R, because the
> latter is freeware (and I don't know for how long I will still have a
> license). However, if I cannot seem to resolve my problem (see below), I
> might have to go back to using Primer/Permanova.
> 
> If I run pairwise permanova/adonis tests on my data, the number of unique
> permutations is small (I have two groups, each group has 3 observations)
> and the minimum P value I can get is larger than the alpha value I (and
> most people) that I use to determine statistical significance (i.e. 0.05).
> In the manual of the PERMANOVA+ add-on in Primer by Anderson et al. (2008)
>  it is mentioned (page 28 and onwards) that when the number of unique
> permutations is small (<100) than one should preferably interpret the
> Monte-Carlo p value.
> 
> Does anyone know how to do this in R?
> 
> 
> In my internet search I stumbled upon this page:
> http://r.789695.n4.nabble.com/monte-carlo-simulations-in-permanova-in-vegan-package-td4714034.html.
> "JAri Oksanen answered here: 2. If you want to do so, you can generate your
> resampling matrices by hand and use that matrix as the argument of
> permutations=. See the documentations (?adonis) which tells how to do so.
> ", but I don't understand how to this..
> 
> Thank you very much!
> 
> Ellen
> 
>   [[alternative HTML version deleted]]
> 
> ___
> R-sig-ecology mailing list
> R-sig-ecology@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
> 

-- 
Eduard Szöcs

Quantitative Landscape Ecology
Institute for Environmental Sciences
University Koblenz-Landau
Tel: +49 6341 280 31552
Office: Building I, Room 2.27
http://www.uni-koblenz-landau.de/en/campus-landau/faculty7/environmental-sciences/landscape-ecology/Staff/eduardszoecs

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] CCA INTERPRATE

[Eduard Szöcs]

Hai Mahnaz,

12% of the total Variance can be explained by the 3 constraining
environmental variables and 88% residual variance is left.



Cheers,

Eduard

On 18/10/16 09:51, Mahnaz Rabbaniha wrote:
> Dear all
> 
> I have done CCA in vegan on 41 response and 3 environmental variables , i
> have taken this result:
> 
>   Inertia Proportion
> 
> Total  5.8796 1.
> 
> Constrained   0.7079 0.1204
> 
> Unconstrained   5.1717 0.8796
> 
> 
> how do interpate, is this analyses  acceptable or no?
> 
> 
> 
> all the best
> 

-- 
Eduard Szöcs

Quantitative Landscape Ecology
Institute for Environmental Sciences
University Koblenz-Landau
Tel: +49 6341 280 31552
Office: Building I, Room 2.27
http://www.uni-koblenz-landau.de/en/campus-landau/faculty7/environmental-sciences/landscape-ecology/Staff/eduardszoecs

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] get species within sites ordihull polys

[Eduard Szöcs]

Dear Tim,

Sound's like your looking for indicator species ?

Maybe the indicspecies packages suites your needs.

Dufrêne, M and P. Legendre. 1997. Species assemblages and indicator
species: the need for a flexible assymetrical approach. Ecological
Monographs 67: 345 - 366.

De Cáceres, M., P. Legendre, and M. Moretti, 2010. Improving indicator
species analysis by combining groups of sites. Oikos 119: 1674 -1684.

Cheers,

Eduard

On 25/09/15 20:58, Howard, Tim G (DEC) wrote:
> This isn't the best example because species points don't fall in more than 
> one of the black polygons. But, my question: Can I easily identify which blue 
> points (species) fall within the polygon?   Or can I easily identify which 
> species are 'most important' (most abundant?) for defining each of the three 
> groups? 
> 
> Thanks for any pointers
> 
> Tim Howard

-- 
Eduard Szöcs
Quantitative Landscape Ecology
Institute for Environmental Sciences
University Koblenz-Landau
Tel. +49 6341 280 31552
http://www.uni-koblenz-landau.de/en/campus-landau/faculty7/environmental-sciences/landscape-ecology/Staff/eduardszoecs

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] calculate "dispersion index" (betadisper, vegan)

[Eduard Szöcs]

Dear Natalie,

it should be possible to use the betadisper-result to do this.
However, a possibly faster way would be to calculate the area of the
hull directly.

Here's an example:

require(vegan)
data(dune)
data(dune.env)
MDS <- cmdscale(vegdist(dune, method = 'bray'))
plot(MDS, type = 'n')
points(MDS, col = dune.env$Management, pch = 16)


# area only for hull of NM
take <- MDS[dune.env$Management == 'NM', ]
hull <- chull(take)
lines(take[c(hull, hull[1]), ])
require(splancs)
areapl(take[c(hull, hull[1]), ])


HTH,

Eduard

On 29/07/15 19:26, Natalie wrote:
> hi all
> 
> I would like to  calculate a dispersion index or similar from the
> results of the vegan function betadisper.
> i would like to know if there is a possibility to e.g. calculate the
> volume of the community dispersion(area of the hull in the PCoA) from
> the betadisper results.
> common is to show the distance to the centroid.but I have the impression
> the dispersion within communities are not that clear,especially if you
> have similar variance.
> The program primer v5 calculates a disperion index based on  the
> similarity percentage.
> similarity percentage is not an option since it is mainly based on most
> counted species.
> 
> ideas and help are very much appreciated
> cheers natalie
> 
> ___
> R-sig-ecology mailing list
> R-sig-ecology@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
> 

-- 
Eduard Szöcs
Quantitative Landscape Ecology
Institute for Environmental Sciences
University Koblenz-Landau
Tel. +49 6341 280 31552
http://www.uni-koblenz-landau.de/en/campus-landau/faculty7/environmental-sciences/landscape-ecology/Staff/eduardszoecs

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] Using multiple species data for gam

[Eduard Szöcs]

Dear Rajendra,

your post reminds me on constrained additive ordination [1].

[1] Yee, T. W. (2006) Constrained additive ordination. Ecology, 87, 203–213.

All the best,

Eduard Szöcs

On 10/02/15 10:41, Rajendra Mohan panda wrote:
> Dear All
> 
> I have >1000 species with presence and absence (0 or 1) values and with
> seven corresponding predictor variables. If I can run gam/glm for the data
> using all species data simultaneously vs predictors. Data are arranged in
> columns against their GPS locations (see below). I know it is possible to
> do separately for each species.
> 
> Your kind response is highly appreciated.
> 
> Sites  Sp1  Sp2 Sp3 Alt Temp Pptn   Ft
> 1A 0  11 20   30 1000 Evergreen
> 
> With Best Regards
> Rajendra M Panda
> School of Water Resources
> Indian Institute of Technology Kharagpur, India
> 
>   [[alternative HTML version deleted]]
> 
> ___
> R-sig-ecology mailing list
> R-sig-ecology@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
> 

-- 
Eduard Szöcs
Quantitative Landscape Ecology
Institute for Environmental Sciences
University Koblenz-Landau
Tel. +49 6341 280 31552
http://www.uni-koblenz-landau.de/campus-landau/faculty7/environmental-sciences/landscape-ecology/Staff/eduardszoecs

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] Doubt about mle2, optim and convergence

[Eduard Szöcs]

-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Hai Mario,

take a look at Ben Bolkers 'Ecological Models and Data in R'. The
chapter 'Optimization and all that' therein should be a good reference.

Cheers,

Eduard Szöcs



On 12/11/14 13:36, Mario José Marques wrote:
> Hi all!
> 
> I am working with species abundance distribution using biomass as
> abundance measure. For this, I am trying to fit models gamma, weibull
> and lognormal to my data.
> 
> I have problems specially with gamma and weibull. When I try fit this
> models, the mle2 fit without error, but the fit not converge. When I try
> plot profile with profile function, this function return that have found
> a better fit to data. Then I use this mles found by profile and try fit
> with mle2 again. I do this three or four times until no have problems to
> use profile to plot mles found.
> 
> I read Bolker text (2014 bbmle; Maximum likelihood estimation and
> analysis with the bbmle package) and found the parameters control and
> parscale (parameters of optim) that solve the problem to fit. I did this:
> 
> # Search for mle. Without convergence
> g1_gamma <- fitgamma(g1, trueLL = FALSE, start.value = mlegamma(g1))
> 
> # Search for mle. Converged
> g1_gamma <- fitgamma(g1, trueLL = FALSE,
> control=list(parscale=mlegamma(g1)))
> 
> I am using sads package that use mle2 to fit model to data. Function
> mlegamma calculate the start values to mle2 and fitgamma use mle2 with
> dgamma to search mles. (the problem is the same if I use only mle2
> without function of sads package).
> 
> I can not understand what control parameters do. I read something, but
> not understood at all. Could someone explain or recommend me some text
> to understand that?
> 
> Thank you and best regards!
> 
> 
> 
> Mario José
> ...
> Doctoral student in Ecology
> Institute of Biology, Dept. Plant Biology
> State University of Campinas - UNICAMP
> Campinas, São Paulo, Brazil
> 
> _______
> R-sig-ecology mailing list
> R-sig-ecology@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
> 

- -- 
Eduard Szöcs
Quantitative Landscape Ecology
Institute for Environmental Sciences
University Koblenz-Landau
Fortstrasse 7
76829 Landau
Germany
http://www.uni-koblenz-landau.de/campus-landau/faculty7/environmental-sciences/landscape-ecology/Staff/eduardszoecs
-BEGIN PGP SIGNATURE-
Version: GnuPG v1

iQEcBAEBAgAGBQJUY2AvAAoJEK3OBW85la/OMw8H/0Y0YjThHCfgJzFzjEJ1TCaL
gIQxQJNQctobbjeAGbP8hep3e8vGalEo1sYHba21pCCJkNbZ1YhyMAt8AxwoY+Vy
v/eLVZPLLciPjwn5R5kLx3sE9vzdWJltiDibm0PJs2GM0cBppPilpPQDl+A6Z8xN
TQuXuLB6KCMPPXDNOOPAT6eBGt7EnsVUyhyJzKT0yVq7cc7To1O8KRwpD/IpwCyU
QVbT+FbH+f1RR+EgmjeFN4u3ircnNp0n6qaQgbJjTML+2czbPDjq9sLsbYVqrS4k
IYIGiyTj6okZP1LGYLtHO8pgqQZO7somjY0nX/sMm5scgsWeJNZ8etpA3OwBvfg=
=/z/u
-END PGP SIGNATURE-

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [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-boun...@r-project.org 
> [r-sig-ecology-boun...@r-project.org] on behalf of Maas, Kendra 
> [kendr...@mail.ubc.ca]
> Sent: Tuesday, October 28, 2014 4:31 PM
> To: r-sig-ecology@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@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
> 
> ___
> R-sig-ecology mailing list
> R-sig-ecology@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

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

[Eduard Szöcs]

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)
}

# try it
require(mvabund)
# some data
env <- runif(100)
zone <- gl(2, 50)
response <- data.frame(A = rpois(100, 5), B = rpois(100, 1))
# run it
mymv(response, env, zone)
###---------

Hope this helps,

Eduard


-- 
Eduard Szöcs
Quantitative Landscape Ecology
Institute for Environmental Sciences
University Koblenz-Landau
Fortstrasse 7
76829 Landau
Germany
http://www.uni-koblenz-landau.de/campus-landau/faculty7/environmental-sciences/landscape-ecology/Staff/eduardszoecs


On 29/10/14 00:31, Maas, Kendra 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@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
> .
>

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] Errors in running some commands in "dismo" in r

[Eduard Szöcs]

Dear Rajendra,

could you be more explicit about your problem?


For me the example presented in the documentation [here:

http://www.rdocumentation.org/packages/dismo/functions/geoDist]

runs without any problems.

Cheers,

Eduard

On 16/08/14 07:41, Rajendra Mohan panda wrote:
> Dear All/ Dear Prof Hijmans
> 
> I find problems regarding some commands written at page 34 of "dismo"
> package pdf published on dec 2013,
> While I use geoDist fun error comes but using geodist fun it works. I think
> there may be some typological error. However using the following fun cited
> in examples (page 34):
> gd <- geoDist(train, lonlat=FALSE)
> 
> I get following error
> Error in geodist(train, lonlat = FALSE) :
>   unused argument (lonlat = FALSE)
> 
> Kindly advise.
> 
> With best Regards
> Rajendra M Panda
> SWR, Indian Institute of Technology Kharagpur
> 
>   [[alternative HTML version deleted]]
> 
> ___
> R-sig-ecology mailing list
> R-sig-ecology@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] Extract residuals from adonis function in vegan package

[Eduard Szöcs]

Dear Alicia and Jari,

just a thought:
Couldn't be capscale or betadisper be used for this?
 - To obtain the distances to the group centroid?

But than: How to convert this from distances to abundances?

Eduard Szöcs


On 03/18/2014 01:21 PM, Jari Oksanen wrote:
> Dear Alicia Valdés, 
> 
> On 18/03/2014, at 13:53 PM, Alicia Valdés wrote:
>>
>> My problem is that I cannot figure out how to get residual values from the
>> adonis model.
>>
> You cannot get residuals from the output of adonis(). 
> 
> We could change the function so that this is possible, but the current 
> function does not return information for getting residuals. Neither would 
> they be residuals in the traditional meaning of the word as we are dealing 
> with dissimilarities or distances, and these cannot be negative. We got to 
> discuss this with vegan developers.
> 
> Cheers, Jari Oksanen
> ___
> R-sig-ecology mailing list
> R-sig-ecology@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
> 


-- 
Eduard Szöcs
Quantitative Landscape Ecology
Institute for Environmental Sciences
University Koblenz-Landau
Fortstrasse 7
76829 Landau
Germany
http://www.uni-koblenz-landau.de/campus-landau/faculty7/environmental-sciences/landscape-ecology/Staff/eduardszoecs

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] NA error in envfit

[Eduard Szöcs]

Hai,

building and installing from github is quite easy with the devtools
package (thanks to Hadley!):

# install devtools
install.packages('devtools')
require(devtools)

# install vegan from github
install_github('vegan', 'jarioksa')


Cheers,

Eduard Szöcs




On 12/05/2013 02:07 PM, Jari Oksanen wrote:
> Hello,
> 
> I think I saw something like this in autumn (northern hemisphere) when a 
> variable had constant values with no variation, and envfit did not know how 
> to scale the arrow. 
> 
> We fixed this in the development version of vegan in R-Forge on September 29. 
> Unfortunately R-Forge is again dysfunctional and cannot build the package, 
> but if you are able to do that yourself you can try to see if the problem is 
> fixed there. The same files are also in github, but you need to build the 
> package yourself there too. I'm working with a minor release of vegan 
> (2.0-10) which may be published on Monday 9 Dec, but there are no guarantees 
> that it will have this envfit fix or be released like planned (you know, the 
> best laid plans of mice and men...) 
> 
> It may be easiest to see if a constant "variable" is the culprit, and remove 
> that if needed. If this is not the case, we need more info and a 
> *reproducible* example. We haven't got any now, and I cannot reproduce your 
> problem.
> 
> Cheers, Jari Oksanen
> 
> From: r-sig-ecology-boun...@r-project.org 
> [r-sig-ecology-boun...@r-project.org] on behalf of Stephen Sefick 
> [sas0...@auburn.edu]
> Sent: 04 December 2013 22:01
> To: Mitchell, Kendra
> Cc: r-sig-ecology@r-project.org
> Subject: Re: [R-sig-eco] NA error in envfit
> 
> Kendra,
> 
> Something is wrong in X or P; find out what the foreign function call is
>   and then you may be able to track down the offending data problem.
> Maybe a logarithm somewhere? This is probably not much help; I don't
> have much experience with envfit.
> 
> Stephen
> 
> On 12/03/2013 07:06 PM, Mitchell, Kendra wrote:
>> I'm running a bunch of NMS with vectors fitted (slicing and dicing a large 
>> dataset in different ways).  I'm suddenly getting an error  from envfit
>>
>> f.bSBS.org.fit<-envfit(f.bSBS.org.nms, f.bSBS.org.env, permutations=999, 
>> na.rm=TRUE)
>>
>> Error in vectorfit(X, P, permutations, strata, choices, w = w, ...) :
>>NA/NaN/Inf in foreign function call (arg 1)
>> In addition: Warning message:
>> In vectorfit(X, P, permutations, strata, choices, w = w, ...) :
>>NAs introduced by coercion
>>
>> I can plot the NMS and even run ordifit on individual env variables, so 
>> can't figure out what the problem is.   There aren't any NA/NaN/Inf in 
>> either of those data that I can find.  I've tried running it without 
>> na.rm=TRUE and still get the error.  Guidance on how to fix this would be 
>> appreciated.
>>
>> Here's the whole slicing process and str for the data
>>
>>
>> f.bSBS.org<-f.env$zone.hor=="bSBS.1"
>> f.bSBS.org.tyc<-f.tyc[f.bSBS.org,f.bSBS.org]
>> f.bSBS.org.env<-subset(f.env, f.env$zone.hor=="bSBS.1")
>> f.bSBS.org.nms<-metaMDS(as.dist(f.bSBS.org.tyc), k=3, trymin=50, trymax=250, 
>> wascores=FALSE)
>> f.bSBS.org.fit<-envfit(f.bSBS.org.nms, f.bSBS.org.env, permutations=999, 
>> na.rm=TRUE)
>>
>>
>> str(f.bSBS.org.env)
>> 'data.frame':63 obs. of  14 variables:
>>   $ zone : Factor w/ 6 levels "bIDF","bSBS",..: 2 2 2 2 2 2 2 2 2 2 
>> ...
>>   $ site : Factor w/ 18 levels "A7","A8","A9",..: 12 12 12 12 12 12 
>> 12 12 12 12 ...
>>   $ om   : Factor w/ 4 levels "0","1","2","3": 2 2 2 3 3 3 2 2 2 3 
>> ...
>>   $ compaction   : num  1 1 1 1 1 1 1 1 1 1 ...
>>   $ herbicide: num  0 0 0 0 0 0 0 0 0 0 ...
>>   $ horizon  : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
>>   $ Water_content: num  50.3 50.3 50.3 50.1 50.1 ...
>>   $ DNA_ug_g : num  71.2 71.2 71.2 68.6 68.6 ...
>>   $ C: num  30.5 30.5 30.5 28.4 28.4 ...
>>   $ N: num  0.863 0.863 0.863 0.81 0.81 ...
>>   $ pH_H2O   : num  4.63 4.63 4.63 4.49 4.49 ...
>>   $ CN   : num  35.3 35.3 35.3 35.1 35.1 ...
>>   $ f.env$zone   : Factor w/ 6 levels "bIDF","bSBS",..: 2 2 2 2 2 2 2 2 2 2 
>> ...
>>   $ zone.hor : chr  "bSBS.1" "bSBS.1" "bSBS.1" "bSBS.1" ...
>>
>&

Re: [R-sig-eco] modify plot of metaMDS (package vegan)

[Eduard Szöcs]

Hai Elaine,

You can add the sites and species separately to your plot.


So first create an empty plot:

require(vegan)
data(dune)
mds <- metaMDS(dune)
plot(mds, type = 'n')



And then add the information you want (points, text, etc...):

# species as symbols
points(mds, display = 'species', pch = '+', cex = 0.6)

# sites as text
text(mds, display = 'sites')


See also Gavin Simpsons series on ordination plots:
http://www.fromthebottomoftheheap.net/2013/01/12/decluttering-ordination-plots-in-vegan-part-1-ordilabel/


Cheers,

Eduard


On 08/02/2013 08:46 AM, Elaine Kuo wrote:

Hello,

I am running NMDS using metaMDS, based on a matrix.
The matrix has rows as islands and columns as species ID.

I generated a plot of the result of metaMDS using type="t".
However, the island names are mostly covered by species ID and thus become
unclear.

Please kindly advise any method to show island names (row names) only
for the plot of metaMDS.

Thank you

Elaine

[[alternative HTML version deleted]]

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology




--
Eduard Szöcs
Quantitative Landscape Ecology
Institute for Environmental Sciences
University Koblenz-Landau
Fortstrasse 7
76829 Landau
Germany
http://www.uni-koblenz-landau.de/landau/fb7/umweltwissenschaften/landscape-ecology/Staff/eduardszoecs

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] Simper

[Eduard Szöcs]

Hi Maarten,

no there is no such limit in simper().
Computational time increases with groups (for 6 groups there are 15 
comparisons).


HTH,

Eduard Szöcs

On 05/29/2013 01:52 PM, Jong, Maarten de wrote:

Hi all,

I'm wondering how many classes can be used in the simper analysis. I think the 
max is 4 but I have 6 cluster, the analysis works when I merge 2 similar 
clusters.

Kind regards

Maarten de Jong
maarten.dej...@wur.nl<mailto:maarten.dej...@wur.nl>



[[alternative HTML version deleted]]

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology



___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] spatial gls

[Eduard Szöcs]

Hai Viviane,

your link to the data seems to be broken (No such file).

However I tried to reproduce your problem using the data from [1] and 
could not reproduce this:

library(nlme)
snouter.df <- 
read.table('http://www.oikos.ekol.lu.se/appendixdown/snouterdata.txt', 
header=T, sep="\t")
summary(gls(snouter1.1 ~ rain + djungle, data=snouter.df,
correlation=corLin(form=~X+Y)))
dist <- corLin(form=~X+Y)
summary(gls(snouter1.1 ~ rain + djungle, data=snouter.df,
correlation=dist))
summary(lm(snouter1.1 ~ rain + djungle, data=snouter.df))

If you haven't already take a look at [1] and the supplements - very 
helpful.


HTH,

Eduard



[1]
F. Dormann C, M. McPherson J, B. Araújo M, Bivand R, Bolliger J, Carl G, 
et al. Methods to account for spatial autocorrelation in the analysis of 
species distributional data: a review. Ecography 2007;30:609--28.

On 04/10/2013 06:55 PM, Viviane Deslandes wrote:
> Hey all,
>
> I need to perform a spatial gls and I'm using nlme package. I constructed a
>   distance matrix using the corLin command Despite the  R to make the matrix
> it seems that is not using to perform the gls, since that the results of
> GLS with this distance matrix are exactly the same that an OLS. Someone
> could help with this or suggesting another package, tutorial or script?
> Here is the script:
> library(nlme)
>
> dat http://zoo.bio.ufpr.br/pie/data/birds/tyrannidae.dadostyrannidaefinal_semoutlier.txt
> ")
>
> # spatial GLS->peakfreq=NDVI*masse
>
> distgeo=corLin(form=~long+lat)
> peakfreq<-gls(peak_freq~NDVI*masselog10, correlation=distgeo, data=dat)
> summary(peakfreq)
>
> #OLS - peakfreq=NDVI*masse
>
> peakfreq<-lm(peak_freq~NDVI*masselog10,data=dat)
> summary(peakfreq)
>
>   [[alternative HTML version deleted]]
>
> ___
> R-sig-ecology mailing list
> R-sig-ecology@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>


[[alternative HTML version deleted]]

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] Problem with function iucn_summary() of taxize package, returning error in strsplit()

[Eduard Szöcs]

Hai Augusto,

thanks for reporting this bug!

I just pushed a fixed version, where NA is returned when the IUCN-API 
does not return anything.

You can install the newest version from github via:

install.packages("devtools")
require(devtools)
install_github("taxize_",  "ropensci")
require(taxize)


Happy holidays,

Eduard


On 12/19/2012 08:29 PM, Sarah Goslee wrote:

Hi Augusto,

I took a quick look at the code for iucn_summary, and it looks like a bug to me.

I copied the package maintainer on this reply, as is appropriate for
problems with a particular package.

It looks like:
pop <- xpathSApply(h, "//div[@id=\"population\"]/text()[preceding-sibling::br]",
 xmlValue)

is returning an empty list, which then throws an error in the next line:
pop <- do.call(rbind, lapply(strsplit(pop, split = ":"),
 rbind))


It looks like the function works well for complete cases, but has
problems with certain missing/different data (easy mistake I've
committed myself a time or two). You've hit on one of those fringe
cases by accident, and it should be an easy update to the package.

Sarah

On Wed, Dec 19, 2012 at 2:15 PM, Augusto Ribas  wrote:

Hello

I'm using the package taxize and was trying to use the function
iucn_summary() but it work perfectly  with the examples, but when i
try another name i receive the error msg:

#

iucn_summary("Molossus currentium")

Erro em strsplit(pop, split = ":") : argumento modo não caractere
##

My r is in Portuguese, but it say argument mode not a character.

I checked this specie on iucn site and its there, dunno what is the problem.


the examples i cite are:

#
library(taxize)

iucn_summary("Lynx lynx")

iucn_summary("")

iucn_summary("Molossus currentium")
#

first two work properly, but the third give an error.. i have no idea
of what i'm doing wrong.

I'm using ubunto here, with R 2.15.2 from emacs.

Thanks for the attention.

--
Grato
Augusto C. A. Ribas


--
Sarah Goslee
http://www.functionaldiversity.org

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology



___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] Package bio.infer: get.taxonomic() hangs R after editing taxa list

[Eduard Szöcs]

Hai Rich,

Have a look at our recently released taxize-package.
However we some some problems with the ITIS-API, too. (Maybe some server 
issues?)

See https://github.com/ropensci/taxize_/issues/72

However with taxize you can also query the NCBI taxonomy browser as 
alternative:


df <- read.table(header = TRUE, as.is = TRUE, text = 'SVN  
Taxon CountValue

1 WP220110711  Zaitzevia.parvula484
2 WP220110711   Tvetenia109
3 WP220110711Tubificidae   1054
4 WP220110711Sweltsa 11
5 WP220110711 Suwallia.pallidula 32
6 WP220110711 Stempellinella 11 ')
df$Taxon <- gsub("\\.", " ", df$Taxon)

require(taxize)

#query itis
classification(get_tsn(df$Taxon))
### This may fail sometimes. 'Connection reset by peer'
### we are working at it:
### see https://github.com/ropensci/taxize_/issues/72

# query ncbi taxonomy-browser instead
classification(get_uid(df$Taxon))


Hope this helps,

Eduard



On 12/11/2012 11:08 PM, Rich Shepard wrote:

On Mon, 10 Dec 2012, Sarah Goslee wrote:


... get.taxonomic() no longer has an outputFile argument.


  I've read the get.taxonomic() description in bio.infer.pdf and it works
... to a point. There are 3 taxa currently not in the ITIS database (a
freshwater worm, a midge, and one chironomid for which I can't track down
the genus). The first two are going through the validation process and 
will

be added to ITIS Real Soon Now. In the meantime, I don't know how to
proceed. So I've copied Lester Yuan on this message.

  I tell get.taxonomic() that I'm finished editing and it flashes a 
Tcl/Tk

window then ists the taxa not in the data base, but does not return the R
prompt:


jerall <- get.taxonomic(bioinfer)

The following taxa are not in ITIS:
EISENIELLA
RADOTANYPUS
ZIATZEVIA

  I've not found a key combination that lets me gracefully exit R (or 
emacs
for that matter); I can only kill the process. Nothing I've seen in 
the help

files appears relevant.

  Any suggestions?

Rich

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology



___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] Applied Environmental Data Analyses: Water Chemistry and Biota

[Eduard Szöcs]

Hai Richard,


have a look at:

Borcard, Daniel, Francois Gillet, and Pierre Legendre. /Numerical 
Ecology with R (Use R)/. 1st Edition. Springer, 2011.

which is the little brother to Legendre's Numerical Ecology, applying 
these methods with R.


A very nice freely available book using R is:

Kindt, Roeland, and Richard Coe. /Tree Diversity Analysis: [a Manual and 
Software for Common Statistical Methods for Ecological and Biodiversity 
Studies]/. Nairobi and Kenya: World Agrofirestry Centre, 2005. 
Available here: 
http://www.worldagroforestry.org/downloads/publications/PDFs/MN08242.PDF


HTH,

Eduard




On 12/06/2012 02:26 AM, Rich Shepard wrote:
>   I left academia and basic ecological research decades ago and now work
> with environmental data collected by companies in compliance with 
> regulatory
> permit requirements. I've bought and read (mostly) all the books I could
> find on ecological analyses using R (all but one of Alain Zuur's books,
> Legendre & Legendre's 2nd English edition, Mike McCarthy's Bayesian 
> methods
> book, and Ben Bolker's book) but cannot find any references to 
> 'communities'
> in the indices. I'd greatly appreciate pointers to sources appropriate 
> for
> environmental data (which is much sloppier than ecological research 
> data).
>
>   The last time I addressed community analyses was my post-doc research
> which I published in Freshwater Biology in 1984. Only within the past 
> year
> have my clients needed to address issues using benthic macroinvertebrate
> assemblages (and fish) in streams. And, since I work by myself, I've 
> no one
> with whom to share ideas and discuss approaches; perhaps there's a better
> forum than this mail list for this.
>
>   The available benthic data has little taxonomic consistency below the
> family level. I want to use functional feeding groups rather than taxa as
> the basis of comparison because those better reflect conditions in each
> stream (collections of biota are made only once per year), and I want to
> examine correlations and cause-and-effect relations between biotic
> assemblages and water chemistry. There are only a few fish collections,
> tool, in the available data.
>
>   All ideas are certainly welcome!
>
> TIA,
>
> Rich
>


[[alternative HTML version deleted]]

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] Simper analysis with Morisita-Horn

[Eduard Szöcs]

Hai Vesna,


your code looks reasonable.

However it will fail, when you try running the permutations:

'Error: object 'gpa' not found'

This comes from your second inserted code block (should be 'ga' ect).

However in the current version it is not intended that the 'normal' user 
knows about the permutations. And one must modify the function (change 
the second line) in order to get the permutations running.


So you can skip your second code block when you don't want to run the 
permutations (and leave the function as is).


Also have a look at the development-version of vegan at github:
https://github.com/jarioksa/vegan/blob/master/R/simper.R
for a slightly changed version.


Hope this helps,

Eduard




On 12/03/2012 10:18 AM, Vesna Gagic wrote:

Dear ecology fellows,

I tried to implement Morisita-Horn distance (instead of Bray that is in the 
current version)  in the code for the Simper analysis in vegan. I would be very 
grateful if someone can check if the code  is right.

function (comm, group, ...)
{
 if (any(rowSums(comm, na.rm = TRUE) == 0))
 warning("you have empty rows: results may be meaningless")
 permutations <- 0
 trace <- FALSE
 comm <- as.matrix(comm)
 comp <- t(combn(unique(as.character(group)), 2))
 outlist <- NULL
 P <- ncol(comm)
 nobs <- nrow(comm)
 if (length(permutations) == 1) {
 perm <- shuffleSet(nobs, permutations, ...)
 }
 else {
 perm <- permutations
 }
 if (ncol(perm) != nobs)
 stop(gettextf("'permutations' have %d columns, but data have %d rows",
 ncol(perm), nobs))
 nperm <- nrow(perm)
 if (nperm > 0)
 perm.contr <- matrix(nrow = P, ncol = nperm)
 for (i in 1:nrow(comp)) {
 group.a <- comm[group == comp[i, 1], ]
 group.b <- comm[group == comp[i, 2], ]
 n.a <- nrow(group.a)
 n.b <- nrow(group.b)
 contr <- matrix(ncol = P, nrow = n.a * n.b)
 for (j in 1:n.b) {
 for (k in 1:n.a) {

  ### Morisita-Horn

 aN <- sum(group.a[k, ])

 bN <- sum(group.b[j, ])
 da <- sum(group.a[k, ]^2)/aN^2
 db <- sum(group.b[j, ]^2)/bN^2
 top <- (group.a[k, ] * group.b[j, ])
 contr[(j - 1) * n.a + k, ] <- 2 * top/((da + db) * aN * bN)
#
 }
 }
 average <- colMeans(contr)
 if (nperm > 0) {
 if (trace)
 cat("Permuting", paste(comp[i, 1], comp[i, 2],
   sep = "_"), "\n")
 contrp <- matrix(ncol = P, nrow = n.a * n.b)
 for (p in 1:nperm) {
 groupp <- group[perm[p, ]]
 ga <- comm[groupp == comp[i, 1], ]
 gb <- comm[groupp == comp[i, 2], ]
 for (j in 1:n.b) {
   for (k in 1:n.a) {
  
  ### Morisita-Horn


 aNp <- sum(gpa[k, ])
 bNp <- sum(gpb[j, ])
 dap <- sum(gpa[k, ]^2)/aNp^2
 dbp <- sum(gpb[j, ]^2)/bNp^2
 topp <- (gpa[k, ] * gpb[j, ])
 contrp[(j - 1) * n.a + k, ] <- 2 * topp/((dap + dbp) * aNp * bNp)
#
   }
 }
 perm.contr[, p] <- colMeans(contrp)
 }
 p <- (apply(apply(perm.contr, 2, function(x) x >=
 average), 1, sum) + 1)/(nperm + 1)
 }
 else {
 p <- NULL
 }
 overall <- sum(average)
 sdi <- apply(contr, 2, sd)
 ratio <- average/sdi
 ava <- colMeans(group.a)
 avb <- colMeans(group.b)
 ord <- order(average, decreasing = TRUE)
 cusum <- cumsum(average[ord]/overall)
 out <- list(species = colnames(comm), average = average,
 overall = overall, sd = sdi, ratio = ratio, ava = ava,
 avb = avb, ord = ord, cusum = cusum, p = p)
 outlist[[paste(comp[i, 1], "_", comp[i, 2], sep = "")]] <- out
 }
 attr(outlist, "permutations") <- nperm
 class(outlist) <- "simper"
 outlist
}

Regards,
Vesna
[[alternative HTML version deleted]]

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology



___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] summing F stats and permutation

[Eduard Szöcs]

On 11/29/2012 05:41 PM, Steve Brewer wrote:
> [...]
>
> Incidentally, I was initially very skeptical of the sumF approach until I
> started comparing sumF results to perMANOVA results (using PC-Ord) on
> several different datasets of mine and got strikingly similar results
> between the two techniques. It seems to go against what a lot of
> statisticians are saying about the unique value of multivariate stats. I
> don't have a satisfactory explanation for why it seems to work with my
> data.
I think this is similar to the approach of David Warton:

" The problem of testing for an effect of a difference in mean
abundance between a set of multivariate samples (Fig. 4c) can
be addressed using the simple approach of fitting separate
models to each taxon, then summing the test statistics, and
using resampling to assess the significance of this multivariate
test statistic (Warton 2011). This 'sum-of-LR' statistic is analogous
to the 'sum-of-F' approach based on analysis of variance
(Edgington & Onghena 2007), which was previously shown to
have better power properties than distance-based analyses
(Warton & Hudson 2004). The key advantage however of the
sum-of-LR approach is that we can correctly model the mean--
variance relationship, leading to better power properties than
sum-of-F statistics when sampling is unbalanced (Warton
2008b). The effectiveness of the approach is demonstrated in
Fig. 4c, where compared with PERMANOVA, power of the
sum-of-LR statistic was usually much higher and tended to
vary much less as a function of the variance of the effect taxon.
See Warton (2011) for further simulations demonstrating
power advantages. "

- cited from:
Warton, David I., Stephen T. Wright, and Yi Wang. "Distance-based 
Multivariate Analyses Confound Location and Dispersion Effects." 
/Methods in Ecology and Evolution/ (2011).

Have a look at the mvabund-package, implementing their methods:
Wang, Yi, Ulrike Naumann, Stephen T. Wright, and David I. Warton. 
"Mvabund- an R Package for Model-based Analysis of Multivariate 
Abundance Data." /Methods in Ecology and Evolution/ (February 2012). 
doi:10./j.2041-210X.2012.00190.x.

Eduard

>
> Steve
>
>
> J. Stephen Brewer
> Professor
> Department of Biology
> PO Box 1848
>   University of Mississippi
> University, Mississippi 38677-1848
>   Brewer web page - http://home.olemiss.edu/~jbrewer/
> FAX - 662-915-5144
> Phone - 662-915-1077
>
>
>
>
> On 11/29/12 9:43 AM, "Jari Oksanen"  wrote:
>
>> Steve,
>>
>> This is R, so it is not about whether this can be done, but how this can
>> be done. Unfortunately, doing exactly this requires some fluency in R.
>> Doing something similar is very simple.
>>
>> The description of your problem sounds very much like the description of
>> permutation test in redundancy analysis (RDA). The difference is that in
>> RDA you sum up nominators and denominators before getting the ratio, but
>> in your model you sum up the ratios. So in RDA test you have (num_1 +
>> num_2 + ... + num_p)/(den_1 + den_2 + ... + den_p), and in your
>> description you have num_1/den_1 + num_2/den_2 + ... + num_p/den_p. The
>> former test in canned for instance in the vegan package, but the latter
>> you must develop yourself (and the former method of summing variances
>> instead of their ratios feels sounder). It would not be too many lines of
>> code to fit your code, though. Please note that RDA works by fitting
>> linear models for each species independently so that you can get all
>> needed statistics from a fitted RDA in the vegan package (function rda).
>> The following function extracts F-values by species from a fitted
>> vegan:::rda() result object:
>>
>> spF <-
>> function (object)
>> {
>>inherits(object, "rda") || stop("needs an rda() result object")
>>df1 <- object$CCA$qrank
>>df2 <- nrow(object$CA$Xbar) - df1 - 1
>>num <- colSums(predict(object, type="working", model="CCA")^2)
>>den <- colSums(predict(object, type="working", model="CA")^2)
>>(num/df1)/(den/df2)
>> }
>>
>> HTH, Jari Oksanen
>> 
>> From: r-sig-ecology-boun...@r-project.org
>> [r-sig-ecology-boun...@r-project.org] on behalf of Steve Brewer
>> [jbre...@olemiss.edu]
>> Sent: 29 November 2012 16:42
>> To: r-sig-ecology@r-project.org
>> Subject: [R-sig-eco] summing F stats and permutation
>>
>> Dear Colleagues,
>>
>> I'm wondering if anyone in this group has developed code for doing a sumF
>> test for examining community responses in an experiment. For those not
>> familiar, sumF is a simple univariate alternative to MANOVA and perMANOVA,
>> wherein univariate ANOVAs and their associated F statistics are calculated
>> for each species' abundance and then the F statistic for each effect is
>> summed over all species. The significance of the resulting summed F
>> statistic is then evaluated using random permutation. The summed F
>> statistic
>> is interpreted as an overall community response to the treatments, whereas
>> the F sta

Re: [R-sig-eco] Use cor.test for multiple parameters?

[Eduard Szöcs]

On 11/15/2012 01:49 PM, Kathleen Regan wrote:

In the version of R
that I am using, there is no such thing as "correlation.r".

Hi Kathleen,

I suppose 'correlation.r' is a file with R-Code written by your 
colleague (probably containing the definition for the function 'cor_test()'.
If you don't have this file anymore, it wouldn't be hard to write some 
R-Code doing what you want (perhabs ~15Lines of Code).



Eduard

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] IndVal groups & Plotting an RDA with labels

[Eduard Szöcs]

Hai Caitlin,

there is a missing " in your  *#plot.cca (Vegan package) * - Code, hence 
the error.

plot.cca() takes the rownames of the species table as labels in the triplot.


So something like this might work for you:

require(vegan)
data(dune)
data(dune.env)
# plot.cca takes the rownames of the species table as labels in a triplot
rownames(dune) <- LETTERS[1:nrow(dune.env)]

RDA <- rda(dune ~ ., dune.env)
plot(RDA, type = "text")


Also have a look at Gavin Simpsons blog-post about customizing 
ordination plots.
http://ucfagls.wordpress.com/2012/04/11/customising-vegans-ordination-plots/

HTH,

Eduard





On 03/07/12 15:12, Caitlin Porter wrote:
> Dear all,
>
>
> I am working with a large data set to examine plant community structure on
> barrens habitats and have a couple of questions regarding Indval (labdsv
> package) on Wards Clusters (hclust function in stats package) and plotting
> axes of an RDA result with plot, plot.cca or esqplot functions (vegan,
> MASS, graphics packages)
>
> 1.  I have run IndVal (package labdsv) on a Ward's Cluster Analysis to
> determine indicator species for groups. I selected k=5 groups as my most
> meaningful and conservative (sample size, etc) classification. However,
> according to my average silhouette width, k=12 is the most optimal number
> of groups.  I ran an IndVal on k=12 out of curiosity and notice some
> interesting patterns I would like to explore further but I am having
> difficulty understanding the output. The groups are labeled in the
> indicator species output as "group 1,2,3,412". I do not understand how to
> determine which group is associated with which cluster (e.g.. in my Ward's
> dendrogram -- which is group 11?). It is somewhat obvious when k is
> relatively small because of the order the groups are clustered in, branch
> height in the dendrogram, and species frequency and abundances, however I'd
> like to know for larger groups. It would be ideal if I could even label the
> dendrogram with these groups. I've seen examples of these in a couple of
> papers with color coded boxes, but I can't  seem to figure out how to code
> it myself.
>
>   2.  My second question relates to plotting an RDA. I have been able to
> run an RDA in vegan package successfully but unable to plot it in a way
> that I can interpret. I need to label sites, species (response matrix) and
> environmental variables/PCA axes (explanatory matrix).
>
>   So far, I've only been able to label either the response matrix or the
> explanatory matrix in my graphs, but not all 3 sets of points. I've tried
> modifying plot function and code from Borcard et al 2011, (Numerical
> ecology with R), esqplot code for MASS package and plot.cca in Vegan
> package.  I would prefer to use esqplot since I understand  already how to
> better customize it, but I'm just looking to get any graph I can read at
> this point.  When I use the plot function from Borcard et al. I see PC axes
> names only. When I use esqplot, I see species names only.  I also tried
> plot.cca in vegan package but wasn't able to call up a graph. This code
> looks like a great way to do it, but I'm not sure what I'm doing, *e.g*.
> what to put in for const. or what the 'unexpected symbol' error means.
>
>   This old thread asks a similar question (
> https://stat.ethz.ch/pipermail/r-help/2009-February/188282.html), but I'm
> not sure I understand its solution and have approximately 300 species so
> providing a separate name for each individually might not be feasible. This
> other thread asks another similar question (
> http://r.789695.n4.nabble.com/RDA-Triplot-td3055474.html) but the author
> finds an error is generated specifying that biplot is not an appropriate
> method.
>
> I have included my code for question 2 below.Any help would be very much
> appreciated!
>
>
> Sincere thanks,
>
>
> Caitlin
>
>
> *#esqplot (MASS package) *
>
> library(MASS)
>
> #subset species and sites scores from the rda for first 10 RDA axes
>
> sr <- scores(c1.rda, display = c("sites", "species"), choices =
> c(1,2,3,4,5,6), scaling = 2)
>
> sites.only<- as.data.frame(sr$sites)
>
> srsp <- as.data.frame(sr$species) # data frame with just the species in it
>
> c1.site <- c1[,1] # object with just the site names from the original data
> set
>
> cp.m <- merge(c1.site, sites.only, by=0, sort=FALSE) # merged site object
> with site names
>
> eqscplot(cp.m$RDA1, cp.m$RDA2, xlim=c(-1, 1), ylim=c(-1, 1), col="blue",
> xlab="RDA Axis 1", ylab="RDA Axis 2", cex=0.3) # defining the variables &
> limits of plot and what the symbols look like
>
> text(cp.m$RDA1, cp.m$RDA2, labels=cp.m$x, col="black", cex=0.3) #adding
> names on plots
>
> text(srsp$RDA1, srsp$RDA2, labels=rownames(srsp),col= "red", cex=0.3) #adding
> names/species
>
> # Error in text.default(cp.m$RDA1, cp.m$RDA2, labels = cp.m$x, col =
> "black",  :   zero length 'labels'
>
>
>
> *#Borcard et al. 2011 plot function*
>
> plot(c1.rda, main= "Triplot RDA - scaling 2 - wa scores")
>

Re: [R-sig-eco] Help with function to webscrap

[Eduard Szöcs]

Hai Augusto,

regarding question #3:
You could use the red list API with RCurl and XML packages.
Here is an example:

> require(RCurl)
> require(XML)
> get_IUCN_status <- function(x) {
+   spec <- tolower(x)
+   spec <- gsub(" ", "-", spec)
+   url <- paste("http://api.iucnredlist.org/go/";, spec, sep="")
+   get <- getURL(url, followlocation = TRUE)
+   h <- htmlParse(get)
+   status <- xpathSApply(h, '//div[@id ="red_list_category_code"]', 
xmlValue)

+   return(status)
+ }
>
> get_IUCN_status("Panthera uncia")
[1] "EN"

For more resources just type 'webscraping R' in your favourite search 
engine.


HTH,

Eduard

On 26/06/12 20:57, Augusto Ribas wrote:

Hello.
I'm haveing problems with a function to do webscrap.
I have a long list like this example:

data<-data.frame(especie=c("Rana pipiens","Rana vaillanti","Ctenosaura
similis","Bos taurus"),group=c("sapo","sapo","reptil","mamifero"))

And, as some species names are out of data, i trying to make a
function to check catalogue of life (http://www.catalogueoflife.org/)
and get the current names.
This have some problems, like species name that split, but help as a
first check.

So i made this function to web scrap the data.
Its simple, it search the site, makeing a link with the keywords, then
enter the first link of the list of results produced and get the
accepted name and author, giveing the results as a list.
for example:


sp.check("Rana pipiens")

$sp.aceito
[1] "Lithobates pipiens"

$autor
[1] "Schreber, 1782"

But sometimes the function cannot acess the internet, and give a error.

I'm made this function trying to copy some examples on foruns, but i
have some doubts:

01) How do i supress the readlines() warnings?

02) How can i make the function try again when cannot acess internet,
or just print something like "Cant acess internet", or when i try
something like:

data$check<-NA
for(i in 1:nrow(data)) {
  data$check[i]<-sp.check(data$especie[i])
  }

the loop dont stop.
I made a short list, but when with 500 or more lines it usually stop
in the middle.

03) Anyone have an example how to scrap http://www.iucnredlist.org/
the status of species, as it does not use the keyword in the link? Is
there any tutorial simple for someone without any background on
programing or computer science?


Well thanks for the attention.

#função sp.check

sp.check<-function(especie) {
#split species name
especie<-as.character(especie)

gen<-strsplit(especie,"\\ ")[[1]][1]
esp<-strsplit(especie,"\\ ")[[1]][2]

#makeing first link
link<-paste("http://www.catalogueoflife.org/col/search/all/key/",gen,"+",esp,"/match/1",sep="";)
link <- iconv(link, 'latin1', 'UTF-8')
Encoding(link) <- 'bytes'

#reading table of results
pagina <- readLines(url(link))

n.linhas<-which(pagina%in%"")

#is there any results?
if(length(n.linhas)>0) {

pag.sp<-strsplit(pagina[n.linhas[1]+1],'\\"')[[1]][2]

#second link
link2 <- paste( "http://www.catalogueoflife.org",pag.sp,sep="";)
link2 <- iconv(link2, 'latin1', 'UTF-8')
Encoding(link2) <- 'bytes'
link2

#read
pagina2 <- readLines(url(link2))

#get line of interest
linha2<-grep('(accepted name)',pagina2)
sp.final<-pagina2[linha2]

#get species name
corte1<-strsplit(sp.final,'')[[1]][2]
sp.aceito<-strsplit(corte1,'')[[1]][1]

#get author
corte2<-strsplit(sp.final,'\\(')[[1]][2]
autor<-strsplit(corte2,')')[[1]][1]
}else {
sp.aceito<-c("Não encontrado")
autor<-c("Não encontrado")
}
return(list(sp.aceito=sp.aceito,autor=autor))
}

--
Grato
Augusto C. A. Ribas

Site Pessoal: http://augustoribas.heliohost.org
Lattes: http://lattes.cnpq.br/7355685961127056

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology



___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] CCA and related ordination techniques

[Eduard Szöcs]

Hai Marc,

have a look at this Minireview:

@article{ramette2007multivariate,
  title={Multivariate analyses in microbial ecology},
  author={Ramette, A.},
  journal={FEMS microbiology ecology},
  volume={62},
  number={2},
  pages={142--160},
  year={2007},
  publisher={Wiley Online Library}
}


I also like Tree Diversity Analysis:
http://www.worldagroforestry.org/resources/databases/tree-diversity-analysis



HTH,

Eduard



On 21/03/12 13:22, Marc Taylor wrote:

Dear R-sig-ecology group,

Can anyone recommend a good reference on ordination techniques? In
particular, I could use some clear explanations as to the
differences/similarities between Canonical Correlation Analysis, Redundancy
Analysis, and Canonical Correspondence Analysis. I'm hoping to explore some
of the methods included in the vegan package but am lacking some basic
knowledge beyond those methods offered in the program PRIMER.

Cheers,
Marc

[[alternative HTML version deleted]]

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology



___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] test for for differences in floristic composition between factors with UNBALANCED sample sizes

[Eduard Szöcs]

I have no experience with it, but have a look at the mvabund-package:

http://onlinelibrary.wiley.com/doi/10./j.2041-210X.2012.00190.x/abstract


Cheers,

Edi





On 09/03/12 11:05, cristabel.duran wrote:
> Dear R-ecology-list,
>
> I have been looking already for a while in Internet for an alternative
> to Permanova, when the data to analyze has unbalanced sample size
> between factors, but I still do not find good information about it.
> It seems that an alternative should be anosim, but Jari does not
> recommend the use of this function.
>
> I really would appreciate the help you can give concerning this topic.
> Thanks in advanced,
>
> Cristabel.
>
>
>
> ___
> R-sig-ecology mailing list
> R-sig-ecology@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology


[[alternative HTML version deleted]]

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] SIMPER analysis

[Eduard Szöcs]

Hai Gian,

This is implemented in the newest version (2.0-3) of vegan.
Also note this about simper:
https://github.com/jarioksa/vegan/issues/3

Let me know if your facing any problems with this function.


Cheers,

Edi



On 09/03/12 09:51, Gian Maria Niccolò Benucci wrote:

Hi members,

Do anyone know if SIMPER analysis (Clarke 1993) is implemented in R?
It is used to address which species are responsible for the differences
between different groups.
Thank you so much,



___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] split-split-plot design and adonis

[Eduard Szöcs]

Also have a look at the new permute package:
http://cran.r-project.org/web/packages/permute/index.html

AFAIK permute is not fully hooked up in vegan at the moment (at least 
not in adonis),
but you could hack the adonis function to use shuffle() instead of 
permuted.index().


The permutation design, as proposed by Kay, would then look like this:

ctrl <- permControl(strata = Site,
within = Within(type = "none"),
blocks = Blocks(type = "free"))
table(Site, perm.Sediment = Sediment[shuffle(nobs, control = ctrl)])


HTH,

Edi

 





On 05/03/12 11:45, Kay Cichini wrote:

Hi Arnaud and list-members,

what about the following:

library(vegan)

Sediment<-as.factor(rep(c("Sed1","Sed2","Sed3"),each=36))
Site<-as.factor(rep(c("Site1","Site2","Site3","Site4","Site5","Site6","Site7","Site8","Site9"),each=12))
Hydrology<-as.factor(rep(rep(c("Dry","Wet"),each=6),9))
Depth<-as.factor(rep(rep(c("D1","D2"),each=3),18))

nobs = length(Site)

# depth and hydrology are nested within sites:
table(Site, Depth)
table(Site, Hydrology)

# so the following call is ok as it restricts permutation on sites -
# that is, levels of depth and hydrology are shuffled only within and not
# across sites:
set.seed(123)
community<- matrix(rnorm(nobs * 3, 10, 1), nobs, 3)
adonis(community ~ Hydrology * Depth, perm = 999, strata = Site)

# sediment is different:
table(Site, Sediment)

# to allow for this structure you would need to shuffle sites as
# a whole. the below restriction will randomly assign levels of sediment
# to the sites - like this:
ctrl = permControl(strata = Site, permute.strata = TRUE)
table(Site, perm.Sediment = Sediment[permuted.index2(nobs, control = ctrl)])

### to test sediment you'd need a customized code:
(fit<- adonis(community ~ Sediment, permutations = 1))

### number of perms
B<- 999

### setting up frame which will be populated by
### random r2 values:
pop<- rep(NA, B + 1)

### the first entry will be the true r2:
pop[1]<- fit$aov.tab[1, 5]

for(i in 2:(B+1)){
  idx<- permuted.index2(nobs, control = ctrl)
  fit.rand<- adonis(community ~ Sediment[idx], permutations = 1)
  pop[i]<- fit.rand$aov.tab[1,5]
}

### get the p-value (H0: Sediment has no effect on location of communities):
print(pval<- sum(pop>= pop[1]) / (B + 1))

### make a histogram to see random R2-values and the true one:
hist(pop, xlab = "Population R2")
abline(v = pop[1], col = 2, lty = 3)

What do you think?
Kay


2012/3/4 Arnaud Foulquier


I have run the two analysis. The only differences are the p-values that
appear in the anova tables. Other values (Df, SumsOfSqs, MeanSqs, F.Model)
are not modified because they are calculated on the "original" (not
permuted) data.

The p-value is calculated by comparing the value of F obtained with the
original data to the distribution obtained by permutations. If i understand
well, the strata argument is used to constrain the permutations of the
original data within the groups specified in the strata argument. The
distribution of the F-statistic that is generated by permutations is not
the
same according to whether the strata argument is specified or not. This is
why the p-values are modified when the strata argument is specified.

However, I'm not sure that strata=Site_Hydro correctly specify the nested
structure of my data.





--
View this message in context:
http://r-sig-ecology.471788.n2.nabble.com/split-split-plot-design-and-adonis-tp7332029p7342690.html
Sent from the r-sig-ecology mailing list archive at Nabble.com.

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology


[[alternative HTML version deleted]]

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology



___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] betadisper() - tests for multivariate homogenity

[Eduard Szöcs]

Hai Gian,

On 02/03/12 17:30, Gian Maria Niccolò Benucci wrote:
> Hi Members,
>
> On 24 samples from two different fungal ectomycorrhizal communtites
> (natural vs. cultivated), I run betadisper() to test multivariate
> homogeneity of group dispersions.
> To test if one group is more variable than the other, I run anova() of the
> distances to group centroids. Moreover I tried o run a permutest() to
> generate a permutation distribution of F under the Null hypothesis of no
> difference in dispersion between groups.
> Maybe is a stupid question but I can't explain why with both F and
> non-parametric statistics I got same results?
What do you mean, with the same results? The p-values are slightly 
different. In the parametric test, this is assed by the F-Distribution 
[see pf(18.485, 1, 22, lower.tail = FALSE)].

In the Permutationtest, you generate a distribution of F-values by 
shuffling the group vector and compare your original F with this 
generated distribution. Of course the F is the same in both methods, 
only how the p-values is assessed differs.

HTH,

Edi

>
>> betadisper2_cent<- betadisper(dist, groups, type = c("centroid"))
>> anova(betadisper2_cent)
> Analysis of Variance Table
>
> Response: Distances
>Df  Sum Sq Mean Sq F valuePr(>F)
> Groups 1 0.20906 0.20906  18.485 0.0002903 ***
> Residuals 22 0.24881 0.01131
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
>> permutest(betadisper2_cent, control=permControl(nperm=), pairwise=T)
> Permutation test for homogeneity of multivariate dispersions
>
> No. of permutations: 
>
>  STRATA 
> Permutations are unstratified
>
>  SAMPLES 
> Permutation type: free
> Mirrored permutations for Samples?: No
>
> Response: Distances
>Df  Sum Sq Mean Sq  F N.Perm Pr(>F)
> Groups 1 0.20906 0.20906 18.485     3e-04 ***
> Residuals 22 0.24881 0.01131
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Thank you so much,
>
>
>
> ___
> R-sig-ecology mailing list
> R-sig-ecology@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology


[[alternative HTML version deleted]]

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] adjusting for empty rows in metaMDS

[Eduard Szöcs]

Hai Tracy,

You could also add a dummy-species , see [1].

Eduard

[1] K. Robert Clarke, Paul J. Somerfield, M. Gee Chapman (2006) .
On resemblance measures for ecological studies, including taxonomic 
dissimilarities and a zero-adjusted Bray–Curtis coefficient for denuded 
assemblages. Journal of Experimental Marine Biology and Ecology.

Volume 330. Issue 1


On 25/01/12 22:07, Pinney, Tracy A wrote:

Hello,

I am trying to use metaMDS (or isoMDS) for a community data set with 10 
species.  However, I have several sites with no species present (empty rows), 
and this seems to present a problem in the dissimilarity matrix in the form of 
missing values.  Does anyone have suggestions on how to deal with this?  Is it 
appropriate to create a dummy variable in the original community matrix?  If 
so, what is the best way to do this?  Thanks for any help.



Tracy Pinney



PhD candidate

Department of Biology

Baylor University

One Bear Place 97388-7388

Waco, TX  76798-7388


[[alternative HTML version deleted]]

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology



___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] The problems on BIO-ENV procedure by [R]

[Eduard Szöcs]

Regarding your second question:

Did you run metaMDS on the dissimilarity matrix? - Than metaMDS and 
ordiplot have no information about the species.
Try running metaMDS on your community matrix and specify the distance 
argument.



Hope this helps,

Szöcs Eduard


On 12/05/2011 06:57 PM, Dave Roberts wrote:

Shun Tsuobi,

   isoMDS (and routines based on it) will not allow zero 
dissimilarity, which implies perfect replicates.  One alternative is 
to remove one of the replicates, but that may have unsatisfactory 
effect on further analyses.  Alternatively, if you are sure you want 
to keep the duplicate plots you can change the zero to a small value


d[d==0] <- 0.0001

and run the isoMDS on the resulting revised dissimilarity matrix.

Undoubtedly someone from the vegan group will respond to your second 
question.


Dave Roberts

On 12/05/2011 01:28 AM, 坪井　隼 wrote:

Dear Madam / Sir,

I have two questions for the use of the “R” program for the ecological
research. I am studying the relationship between the community
structures of environmental microbes and some environmental conditions.
For this objective, I have known that the BIOENV procedure, which was
developed by Clarke&  Ainsworth (1993), is available on the “R” 
software.



Fist question;
I attempted the use of the procedure to analyze the relationship between
the variation of the microbial community structures and the
environmental factors. However, I can not analyze the relationship based
on isoMDS function. The isoMDS was inacceptable for my dataset. The
command for the BIOENV procedure, which I programmed, and the error
massage I gained was as follow;


library(MASS)
library(vegan)
communitydat<-read.table("C:/Documents and Settings/shuntsuboi/desktop

/bray.txt", head- er=T)

environdat<-read.csv("C:/Documents and Settings shuntsuboi/desktop/ev.

csv",header=T)

env<-environdat[,c("variablesA","variablesB.","variablesC","variablesD

","variablesE")]

d<- vegdist(communitydat, "bray")
isoMDS(d)

error   isoMDS(d) : zero or negative distance between objects 1 and 2

As mentioned above, I can not run the program because the error, which
is “isoMDS(d) : zero or negative distance between objects 1 and 2”,
occurred. What are the ways to solve this problem ? On isoMDS function,
what are the ways that the zero distance of “Bray-Curtis distance” is
acceptable in the function ?

Second question;
Based on the command as above, I ran the metaMDS function. However,
although I could automatically describe the two dimensional ordination
plot figure, I could not gain the X and Y value of the respective plots.
Then, the error massage was shown as follow;
“In ordiplot(x, choices = choices, type = type, display = display,  :
Species scores not available”

What are the ways to solve this problem ?

Sincerely yours,
Shun Tsuboi

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology




___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] Re : Problem to use varpart

[Eduard Szöcs]

Hai Momadou,

have a look at the relaimpo-package.


HTH,

Edi


Am 07.07.2011 14:32, schrieb Gavin Simpson:

On Thu, 2011-07-07 at 12:26 +0100, momadou sow wrote:

Hi,

Thank Gavin for your answer but for your request “but why would you
want to?” : to partition the explained variance in my linear
regression. Thus, I want have the contribution of each explanatory
variable.
Momadou

Don't cherry-pick. What I actually said was why would you want to use
RDA to fit a linear model? Not why do you want to do variance
partitioning.

Anyway, as I said, you can't do what you want with varpart() so you need
to look at how you decompose a linear model (lm()).

G

De : Gavin Simpson
À : momadou sow
Cc : "R-sig-ecology@r-project.org"
Envoyé le : Jeudi 7 Juillet 2011 23h00
Objet : Re: [R-sig-eco] Problem to use varpart

On Thu, 2011-07-07 at 06:33 +0100, momadou sow wrote:

Hi,
I did a linear regression with 5 explanatory variables. Now, to see

the contribution of each variable, I use varpart  from vegan library.
But varpart don’t accepts my 5 explanatory variables, it accept only
4.

1- How must I  do to use my 5 explanatory variables?
2- Is it the sum of variance fraction of each variable must be equal

to 1.

Thanks for your help.

varpart is for redundancy analysis - aka reduced rank regression. I
suppose one can do linear regression with RDA, but why would you want
to?

Anyway, you can't do what you want with varpart(). Look at
interpretting
the linear model directly.

G

--
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
Dr. Gavin Simpson[t] +44 (0)20 7679 0522
ECRC, UCL Geography,  [f] +44 (0)20 7679 0565
Pearson Building,[e] gavin.simpsonATNOSPAMucl.ac.uk
Gower Street, London  [w] http://www.ucl.ac.uk/~ucfagls/
UK. WC1E 6BT.[w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%






___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] New to R

[Eduard Szöcs]

Hai Mary,

afaik there is no "Indval" package.

Have a look at indval() in the labdsv package and the indicspecies package.

Edi


Am 28.06.2011 09:45, schrieb mary anto:

Dear Group,



I am new to R and I was wondering if someone could  suggest the appropriate
   R package  for my analysis.

I  have sampled 166  butterfly species present/sighted  along   line
transects  distributed across  5 habitat types in 3 locations during  August
2008 to May 2011 (area located in the Western Ghats biodiversity hotspot),.

  I would like to find out species indicative of  particular habitats. I plan
to use the  Indval  package. Any other suggestions or comments?

Thanks,

Mary Anto


Project Scientist

Forest Entomology

Division of Forest Health

Kerala Forest Research Institute

Peechi, Thrisur, Kerala

[[alternative HTML version deleted]]

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology



___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] Multiple count if style "queries"

[Eduard Szöcs]

Like this?

dfc2 <- cast(dfm, ECO_NAME ~ variable + value, length)

Eduard


Am 21.05.2011 00:39, schrieb Chris Mcowen:

Hi Eduard,

That is a handy package, i hadn't come across.

The ouput from dcf is perfect.

But for dfm:

 ECO_NAMEOrder   Family
134   Central Range montane rain forests ZingiberalesZingiberaceae
269   Central Range montane rain forests ZingiberalesZingiberaceae
364Biak-Numfoor rain forests  Asparagales  Orchidaceae
500   Central Range montane rain forests  Asparagales  Orchidaceae
566   Central Range montane rain forests   Poales Centrolepidaceae
567  Huon Peninsula montane rain forests   Poales Centrolepidaceae
755   Central Range montane rain forests  Asparagales Lomandraceae
906  Huon Peninsula montane rain forests  Asparagales  Orchidaceae
1079  Central Range montane rain forests ZingiberalesZingiberaceae
1175  Central Range montane rain forests  Asparagales  Orchidaceae
1271  Central Range montane rain forests   Poales   Cyperaceae
1421  Central Range montane rain forests   PoalesEriocaulaceae
1581 Huon Peninsula montane rain forests  Asparagales  Orchidaceae
1621   Biak-Numfoor rain forests  Alismatales  Araceae
1622  Central Range montane rain forests  Alismatales  Araceae
1811  Central Range montane rain forests   Poales   Cyperaceae
1928  Central Range montane rain forests CommelinalesCommelinaceae
2141  Central Range montane rain forests   Poales  Poaceae
2142 Huon Peninsula montane rain forests   Poales  Poaceae
2276  Central Range montane rain forests   Poales  Poaceae
2391 Huon Peninsula montane rain forests  Asparagales  Orchidaceae
2525  Central Range montane rain forests  Asparagales  Orchidaceae
2526 Huon Peninsula montane rain forests  Asparagales  Orchidaceae
2660  Central Range montane rain forests  Asparagales  Orchidaceae

I would want to count how many time a order appears in each eco_name and then 
how many times a family appears in each eco_name and how many times a genus 
appears in each eco_name.

For example to ask the question: How many  Zingiberaceae family are there in 
Central Range montane rain forests? But to query this for each eco_name.

I hope this makes sense.

Chris




On 20 May 2011, at 23:23, Eduard Szöcs wrote:

Using the reshape package:

require(reshape)
dfm<- melt(data2, id = "ECO_NAME")
dfc<- cast(dfm, ECO_NAME~variable, function(x) length(unique(x)))


Eduard

Am 20.05.2011 23:19, schrieb Chris Mcowen:

Dear List,

I am looking to calculate two things from my data frame and was after some 
advice. For the example below i want to know.

1. How many unique Orders/Families and Genera there are per eco-name

2. How many incidences are there for each Order/Family and Genus there are per 
eco-region

I have 650 econame.

I.e for Biak-Numfoor rain forests there are 2 orders, 2 families and two 
genera. Also, Alismatales are represented once, Asparagales once etc etc.

Thanks for any advice.

Chris


ECO_NAMEOrder   
Family  Genus
Biak-Numfoor rain forests   Alismatales 
Araceae Homalomena
Biak-Numfoor rain forests   Asparagales 
Orchidaceae Bromheadia
Central Range montane rain forests  Alismatales Araceae 
Homalomena
Central Range montane rain forests  Asparagales Lomandraceae
Cordyline
Central Range montane rain forests  Asparagales Orchidaceae 
Thelymitra
Central Range montane rain forests  Asparagales Orchidaceae 
Dendrobium
Central Range montane rain forests  Asparagales Orchidaceae 
Vanda
Central Range montane rain forests  Asparagales Orchidaceae 
Bulbophyllum
Central Range montane rain forests  Asparagales Orchidaceae 
Bulbophyllum
Central Range montane rain forests  Asparagales Orchidaceae 
Dendrobium
Central Range montane rain forests  Asparagales Orchidaceae 
Dendrobium
Central Range montane rain forests  CommelinalesCommelinaceae   
Murdannia
Central Range montane rain forests  Poales  
CentrolepidaceaeCentrolepis
Central Range montane rain forests  Poales  Cyperaceae  
Machaerina
Central Range montane rain forests  Poales  Cyperaceae  
Eleocharis
Central Range montane rain forests  Poales  Eriocaulaceae   
Eriocaulon
Ce

Re: [R-sig-eco] Multiple count if style "queries"

[Eduard Szöcs]

Using the reshape package:

require(reshape)
dfm <- melt(data2, id = "ECO_NAME")
dfc <- cast(dfm, ECO_NAME~variable, function(x) length(unique(x)))


Eduard

Am 20.05.2011 23:19, schrieb Chris Mcowen:

Dear List,

I am looking to calculate two things from my data frame and was after some 
advice. For the example below i want to know.

1. How many unique Orders/Families and Genera there are per eco-name

2. How many incidences are there for each Order/Family and Genus there are per 
eco-region

I have 650 econame.

I.e for Biak-Numfoor rain forests there are 2 orders, 2 families and two 
genera. Also, Alismatales are represented once, Asparagales once etc etc.

Thanks for any advice.

Chris


ECO_NAMEOrder   
Family  Genus
Biak-Numfoor rain forests   Alismatales 
Araceae Homalomena
Biak-Numfoor rain forests   Asparagales 
Orchidaceae Bromheadia
Central Range montane rain forests  Alismatales Araceae 
Homalomena
Central Range montane rain forests  Asparagales Lomandraceae
Cordyline
Central Range montane rain forests  Asparagales Orchidaceae 
Thelymitra
Central Range montane rain forests  Asparagales Orchidaceae 
Dendrobium
Central Range montane rain forests  Asparagales Orchidaceae 
Vanda
Central Range montane rain forests  Asparagales Orchidaceae 
Bulbophyllum
Central Range montane rain forests  Asparagales Orchidaceae 
Bulbophyllum
Central Range montane rain forests  Asparagales Orchidaceae 
Dendrobium
Central Range montane rain forests  Asparagales Orchidaceae 
Dendrobium
Central Range montane rain forests  CommelinalesCommelinaceae   
Murdannia
Central Range montane rain forests  Poales  
CentrolepidaceaeCentrolepis
Central Range montane rain forests  Poales  Cyperaceae  
Machaerina
Central Range montane rain forests  Poales  Cyperaceae  
Eleocharis
Central Range montane rain forests  Poales  Eriocaulaceae   
Eriocaulon
Central Range montane rain forests  Poales  Poaceae 
Schizostachyum
Central Range montane rain forests  Poales  Poaceae 
Poa
Central Range montane rain forests  ZingiberalesZingiberaceae   
Alpinia
Central Range montane rain forests  ZingiberalesZingiberaceae   
Curcuma
Central Range montane rain forests  ZingiberalesZingiberaceae   
Amomum
Huon Peninsula montane rain forests Asparagales Orchidaceae 
Taeniophyllum
Huon Peninsula montane rain forests Asparagales Orchidaceae 
Corybas
Huon Peninsula montane rain forests Asparagales Orchidaceae 
Thelymitra
Huon Peninsula montane rain forests Asparagales Orchidaceae 
Glomera
Huon Peninsula montane rain forests Poales  
CentrolepidaceaeCentrolepis
Huon Peninsula montane rain forests Poales  Poaceae 
Poa

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology



___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] Is adonis the right test for me?

[Eduard Szöcs]

Hai Kari,

as far as i know SimPer is not available in vegan.

I have written some code for myself a few months ago to do this.
I think its doing the job, but coding is poor (3 for-loops...).


- Eduard



## SIMPER  , Clarke (1993)   
   
## V 0.3   29.04.2011  Szöcs, Eduard 





###
 Code #
###

simper <- function(comm, groups, select){
  group.a <- as.matrix(comm[groups == select[1], ])
  group.b <- as.matrix(comm[groups == select[2], ])
  n.a <- nrow(group.a)
  n.b <- nrow(group.b)
  P <- ncol(comm)
  me <- matrix(ncol = P)
  md <- matrix(ncol = P)
  contr <- matrix(ncol = P, nrow = n.a * n.b)
  for(j in 1:n.b) {
  for(k in 1:n.a) {
  for(s in 1:P) {
  md[s] <- abs(group.a[k, s] - group.b[j, s])
  me[s] <- group.a[k, s] + group.b[j, s]
  a <- rowSums(me)
  c <- md / a
  contr[(j-1)*n.a+k, ] <- md / a
  }}}
  av.contr <- apply(contr, 2, mean)*100
  ov.av.dis <- sum(av.contr)
  sdi <- apply(contr, 2, sd)
  sdi.av <- sdi / av.contr
  av.a <- colMeans(group.a)
  av.b <- colMeans(group.b)
  dat <- data.frame(av.contr, sdi, sdi.av, av.a, av.b)
  dat <- dat[order(dat$av.contr, decreasing = TRUE),]
  cum <-  cumsum(dat$av.contr / ov.av.dis)*100
  out <- data.frame(dat, cum)
  out
}


###
 Usage 
###

#  simper(comm, groups, select)

# Arguments :


#   * comm : commumity data. species in columns, observations in rows
#   * groups : vector containing treatment per observation
#   * select : vector with length 2, indicating the two groups to compare


# Output :
##

#   * av.contr : average contribution to overall similarity
#   * sdi : standard deviation of contribution
#   * sdi.va : ratio mean / sd = sdi / av.contr.
#   * av.a : average abundance in first group
#   * av.b : average abundance in second group
#   * cum : cumulative contribution


#
 Example 
#

require(vegan)
data(dune)
data(dune.env)
adonis(dune ~ Management, data=dune.env, permutations=99)
levels(dune.env$Management)
simper(dune, dune.env$Management, c("BF", "SF"))




Am 29.04.2011 06:34, schrieb Kari Lintulaakso:

Dear list,

I'm not so familiar with ecological statistics (though I should be),
so I'm looking for a confirmation and support to my decision to move
on with the adonis path.

I've been trying to find a method to analyse differences between
mammalian communities and I think adonis would do the work for me.

My data consists of different localities (each of assigned into a
vegetation class) and number of species in different ecomorphological
classes (groups) I have generated.
There are nine vegetation classes, 52 localities, six mammalian groups
having some 200-300 species in them.
My idea is to use these mammalian groups as "species" in the data
matrix which would look like this:

L   V   G1  G2  G3  ... G6
1   1   0   1   2   ... 8
2   2   1   2   1   ... 3
3   2   1   3   2   ... 4
4   3   5   2   2   ... 9
5   3   5   3   3   ... 8
6   3   4   2   3   ... 7
... 
52  9   8   8   1   ... 0

L=locality id
V = vegetation class
G1...G6 mammalian groups, and in the matrix the counts of species in the group.

First I'm trying to test if there are any differences between the
vegetation groups. For this I'm planning to use a pairwise adonis
comparisons between the different vegetation groups, i.e V1 vs V2, V1
vs V3, ... V8 vs V9. From each of these comparisons I'm taking the Df,
F, P>F etc values into a matrix so that I'm able to see the values
later.

When I get the P-values I suppose I can say that there are differences
between the vegetation classes when the P-value is<  0.05, right?
And if so, I need to find a way to explain which mammalian groups are
responsible for these differences? Does anyone have a good suggestion
for this? I've found one approach which is called SimPer (Similarity
Percentage Analysis), but I'm not sure is it in vegan.

And finally, is there any requirement, how many different data rows I
need for one type of vegetation? In some cases I only have one which
would mean that the within group variation is zero?

Any suggestions are welcome!
-Kari

Kari Lintulaakso, M.Sc.(Biosciences)

Doctoral student
Paleontology and Paleoecology
Department of Geosciences and Geography
University of Helsinki

* Email: kari.lintulaa...@helsinki.fi
* Post: Department of Geology, Gustaf Hällströmin katu 2a (P.O Box
64), 00014 University of Helsinki
* Web page: http://blogs.helsinki.fi/lintulaa/

___

Re: [R-sig-eco] Which analysis should I take?

[Eduard Szöcs]

Hai Yong,

have a look at the vegan package: with the function vegdist() you can 
calculate dissimilarity indices.
There are also some good tutorials on vegan available on the net.


HTH,

Eduard





Am 28.12.2010 06:53, schrieb Yong Zhang:
> Hello all,
>
> Thing is like that now I have 14 sampling sites with the species and 
> environment data. Which analysis is the best choice for determining the 
> difference/dissimilarity among those sites. I have installed the R (2.12.2 
> version), which soft package or function should I take to achieve my 
> objective.
>
> Thanks very much for your time and suggestion.
>
> Best Regards
>
> Yong
> 2010-12-28
>
>
>
>
> Yong Zhang, Ph.D.
> Lab of aquatic insects&  stream ecology
> Dept.of Entonology, Nanjing Agricultural University
> Nanjing, 210095,China
> Phone number:  (+86) -25-84395241
> E-mail:2010202...@njau.edu.cn
>
>
> ___
> R-sig-ecology mailing list
> R-sig-ecology@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology


[[alternative HTML version deleted]]

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: [R-sig-eco] repeated measures NMDS?

[Eduard Szöcs]

.fit, you first need:

detach(package:permute)

Then run:

set.seed(1)
fit2 <- envfit(sol~env, permutations = 999, strata = rep.mes)
fit2

***VECTORS

  NMDS1   NMDS2 r2 Pr(>r)
env 0.28727 0.95785 0.0315  0.321
P values based on 999 permutations, stratified within strata.

a simplistic approach could be, averaging sites, yielding n=24 for 
further testing.


yours,
kay


Gavin Simpson schrieb:

On Thu, 2010-11-11 at 09:50 +0100, Kay Cecil Cichini wrote:

gavin,

sorry - of course it should be permute.strata=F, permuting 
within individual sites!

but despite of this the code should work, doesn't it?
Yes, it should - i.e the permutation will be random within the 
blocks.
Whether it does or not is another matter entirely. AFAICR, this 
option

in permuted.index2() did work.

*But*, this is doing exactly what the original permuted.index() 
does if

you set argument 'strata' to be the grouping factor - in your case:

envfit(sol ~ env, strata = rep.mes, perm = 999)

for example. So there is no need to code up the analysis by hand.

This of course doesn't account for any temporal correlation 
within sites
nor, if the observations on the 24 blocks were made at the same 
times,

that you might want to have the same permutation within each block.

In the former there are 3^24 possible permutations (time series 
within
blocks), so 999 random perms is reasonable, *but* some of these 
random
perms (in permuted.index()) will not respect the temporal 
ordering and

thus you aren't really exploring the correct NULL.

With the latter constraint (same temporal perm within blocks), 
there are
3 random permutations, so good luck getting a reasonable p-value 
from

that.

The two restricted permutations /should/ work correctly, /but/ if 
you

are doing this by hand, I'd look at the functions in the 'permute'
package - only on R-Forge, on the Vegan R-Forge area - as I know the
code to generate these permutations in that package *does* work. 
(It is

the helper cruft around it that needs more work.)

https://r-forge.r-project.org/R/?group_id=68

I've had a busy Summer and not made as much progress as I would have
liked, but I hope to finish this soon and make an initial release to
CRAN so we can start grafting it into vegan.

In the meantime, I can help people try to link the two packages if
needed, but I don't have much time till the end of this month.

G


thanks,
kay



Gavin Simpson schrieb:

On Wed, 2010-11-10 at 23:33 +0100, Kay Cecil Cichini wrote:

hi eduard,

i faced similar problems recently and came to the below solution.
i only try to address the pseudoreplication with an 
appropiate  permutation scheme.
when it comes to testing the interactions, things may get more 
complicated.


the code is in no way approven of, but at least it maybe good 
enough  for a starter.


best,
kay

Hi Kay,

I don't think you have this right.
If you have measured repeatedly, say 5 times, on the same 10
individuals, or if you have ten fields and you take 5 quadrats 
from
each, you need to permute *within* the individuals/fields, not 
permute
the individuals/fields which is what permute.strata does. 
permute.strata

would be useful in evaluating factors that vary at the block
(individuals/fields) level, not at the sample levels.

From what Eduard and you describe, the code you show is not the 
correct

permutation. But I may have misunderstood your intention.

Also, be careful with permuted.index2 - there are reasons why 
it hasn't
been integrated (design goals changed and we felt it would work 
best in
a separate package that others could draw upon without loading 
all of
vegan) and the code has festered a bit and may contain bugs - 
buyer

beware!

G


library(vegan)

### species matrix with 5 sp.
### one env.variable
### a factor denoting blocks of repeated measurments

sp<-matrix(runif(24*3*5,0,100),24*3,5)
env<-rnorm(24*3,10,2)
rep.mes<-gl(24,3)

### NMDS:
sol<-metaMDS(sp,trymax=5)

fit<-envfit(sol~env)
plot(sol)
plot(fit)

### testing code for appropiate randomization,
### permuting blocks of 3 as a whole:
permuted.index2(nrow(sp),permControl(strata = 
rep.mes,permute.strata=T))

correctly, this should say:
### testing code for appropiate randomization,
### permuting within sites:
permuted.index2(nrow(sp),permControl(strata = rep.mes))


B=4999

### setting up frame for population of r2 values:
pop<-rep(NA,B+1)
pop[1]<-fit$vectors$r


and:

### loop:
for(i in 2:(B+1)){
fit.rand<-envfit(sol~env[permuted.index2(nrow(sp),
  permControl(strata = rep.mes))])
   pop[i]<-fit.rand$vectors$r
}

### p-value:

>> print(pval<-sum(pop>pop[1])/B+1)

here a bracket was missing:
print(pval<-sum(pop>pop[1])/(B+1))


### compare to anti-conservative p-value from envfit(),
### not restricting permutations:

envfit(sol~env,perm=B)


Zitat von Eduard Szöcs :


Thanks, that helped.

permuted.index2() gen

Re: [R-sig-eco] repeated measures NMDS?

[Eduard Szöcs]

Thanks, that helped.

permuted.index2() generates these types of permutations. But envfit() 
does not use this yet.
What if I modify vectorfit() (used by envfit() ) in such a way that it 
uses permuted.index2() instead of permuted.index()?



Eduard Szöcs





Am 08.11.2010 22:01, schrieb Gavin Simpson:

On Mon, 2010-11-08 at 15:39 +0100, Eduard Szöcs wrote:

Hi listers,

I have species and environmental data for 24 sites that were sampled
thrice. If I want to analyze the data with NMDS I could run metaMDS on
the whole dataset (24 sites x 3 times = 72) and then fit environmental
data, but this would be some kind of pseudoreplication given that the
samplings are not independent and the gradients may be overestimated,
wouldn`t it?

For environmental data a factor could be included for the sampling
dates - but this would not be possible for species data.

Is there an elegant way either to aggregate data before ordination or
to conduct sth. like a repeated measures NMDS?

Thank you in advance,
Eduard Szöcs

Depends on how you want to fit the env data - the pseudo-replication
isn't relevant o the nMDS. If you are doing it via function `envfit()`,
then look at argument `'strata'` which should, in your case, be set to a
factor with 24 levels. This won't be perfect because your data are a
timeseries and, strictly, one should permute them whilst maintaining
their ordering in time, but as yet we don't have these types of
permutations hooked into vegan.

If you are doing the fitting some other way you'll need to include
"site" as a fixed effect factor to account for the within site
correlation.

You don't need to worry about the species data and accounting for
sampling interval. You aren't testing the nMDS "axes" or anything like
that, and all the species info has been reduced to dissimilarities and
thence to a set of nMDS coordinates. You need to account for the pseudo
rep at the environmental modelling level, not the species level.

HTH

G



___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

[R-sig-eco] repeated measures NMDS?

[Eduard Szöcs]

Hi listers,

I have species and environmental data for 24 sites that were sampled thrice. If 
I want to analyze the data with NMDS I could run metaMDS on the whole dataset 
(24 sites x 3 times = 72) and then fit environmental data, but this would be 
some kind of pseudoreplication given that the samplings are not independent and 
the gradients may be overestimated, wouldn`t it?

For environmental data a factor could be included for the sampling dates - but 
this would not be possible for species data.

Is there an elegant way either to aggregate data before ordination or to 
conduct sth. like a repeated measures NMDS?

Thank you in advance,
Eduard Szöcs

___
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology