[R-sig-eco] repeated measures NMDS?
Kay Cecil Cichini
Kay.Cichini at uibk.ac.at
Thu Nov 11 13:03:54 CET 2010
thanks a lot for the illustrative example.
..referring to your quote:
"...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."
..so for the time beeing we assume the former case -
and for the latter there is no way out.
yours,
kay
ps: in germany/austria there are two alternative spellings for the name
kai/kay - beeing a male name opposed to the english kay.
Gavin Simpson schrieb:
> [Apologies - I replied with this only to Kay. Hopefully she won't mind
> receiving it twice!]
>
> On Thu, 2010-11-11 at 10:32 +0100, Kay Cecil Cichini wrote:
>> thanks a lot for your elaborations.
>>
>> of course, envfit(..,strata=rep.mes) does it.
>>
>> then, at least, i consider it exercise for other cases, were you really
>> might need a handmade permutation
>>
>> so, as to round this off, i actually can't analyse this very design in
>> such a way, with the right NULL concerned - but were to go from here?
>
> You could hook up the code in 'permute'. It contains a new
> permuted.index() function (and currently no NAMESPACE, so will overwrite
> (mask) the vegan version if loaded after vegan) which will break all the
> permutation code in vegan).
>
> Here is your example, modified to use the code in permute. I post this
> to illustrate how you'd use the new code. There are lots of examples
> in ?permuted.index (for the permute package, not the vegan package
> version), but *don't* touch the permutation t-test example code as it
> uses permCheck() and it might call allPerms() and allPerms() *IS*
> *WRONG* for some designs --- this is the last bit I need to fix/get
> working before we can make our first release of this code.
>
> HTH
>
> G
>
> Here's the example script:
>
> ## Load packages
> require(vegan)
> require(permute)
>
> ## Data
> set.seed(123)
> sp <- matrix(runif(24*3*5, 0, 100), nrow = 24 * 3, ncol = 5)
> env <- rnorm(24*3, 10, 2)
> rep.mes <- gl(24, 3)
>
> ### NMDS:
> sol <- metaMDS(sp, trymax = 5)
> fit <- envfit(sol~env, permutations = 0) ## perms now won't work!
>
> B <- 999 ## number of perms
>
> ### setting up frame for population of r2 values:
> pop <- rep(NA, B + 1)
> pop[1] <- fit$vectors$r
>
> ## set-up a Control object:
> ctrl <- permControl(strata = rep.mes,
> within = Within(type = "series", mirror = FALSE))
> ## we turn off mirroring as time should only flow in one direction
>
> ## Number of observations
> nobs <- nrow(sp)
>
> ## check it works
> matrix(permuted.index(nobs, control = ctrl), ncol = 3, byrow = TRUE)
> ## Yep - Phew!!!
>
> ### loop:
> set.seed(1)
> for(i in 2:(B+1)){
> idx <- permuted.index(nobs, control = ctrl)
> fit.rand <- envfit(sol ~ env[idx], permutations = 0)
> pop[i] <- fit.rand$vectors$r
> }
>
> ### p-value:
> pval <- sum(pop >= pop[1]) / (B + 1)
> pval
>
> I get:
>
>> pval
> [1] 0.286
>
> Now to compare with the actual permutation you'd have gotten from
> env.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 <szoe8822 at uni-landau.de>:
>>>>>>
>>>>>>> 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 at r-project.org
>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>>>>>>>
>>>>>> _______________________________________________
>>>>>> R-sig-ecology mailing list
>>>>>> R-sig-ecology at r-project.org
>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>
More information about the R-sig-ecology
mailing list