[R-sig-eco] Testing difference between diversity indices with vegan::oecosimu

[Ivailo]

On Thu, Apr 26, 2012 at 12:17 AM, Kay Cichini <kay.cichini at gmail.com> wrote:
> Hello all,
>
> I'd like to test if total diversity differs between two communities. For
> each community several samples were taken and abundances collapsed over
> groups to compute total diversity for each group. I tried to use
> vegan::oecosimu to test non-randomness of my statisitc (difference in
> Simpson-Diversity indices of collapsed abundances) - however, I am not
> quite sure if I oversee posssible pitfalls:
>
> library(vegan)
> data(dune)
>
> # a grouping variable:
> gr <- gl(2, nrow(dune)/2)
>
> divdiff <- function(x) abs(diversity(colSums(x[gr == "1", ]), "simp") -
>                           diversity(colSums(x[gr == "2", ]), "simp"))
> # testing function:
> divdiff(dune)
>
> oecosimu(dune, divdiff, "r2dtable", nsimul = 1999)
> # oecosimu with 1999 simulations
> # simulation method r2dtable
> # alternative hypothesis: true mean is not equal to the statistic
> #           statistic        z     2.5%      50% 97.5% Pr(sim.)
> # statistic   0.00275 -0.20996  0.00013  0.00280  0.01     0.98

Dear Kay,

I am not sure about any possible pitfalls with your approach, but I
have tested the same data using the randomisation functions of the
"rich" library, and found that neither the Simpson diversity nor the
simple species richness differ significantly among the defined groups.

Here are the results following your example:

library(rich)

# prepare data
one <- as.data.frame(dune[gr == "1", ])
two <- as.data.frame(dune[gr == "2", ])

data <- list(one, two)

# compare cumulative species richness
c2cv(com1=data[[1]],com2=data[[2]],nrandom=1999)
#$res
#
#cv1                  27.0000
#cv2                  28.0000
#cv1-cv2              -1.0000
#p                     0.4220 # N.S.
#quantile 0.025       -4.0000
#quantile 0.975        4.0000
#randomized cv1-cv2    0.0225
#nrandom            1999.0000

# compare the Simpson diversity
simp.one <- diversity(dune[gr == "1", ], "simp")
simp.two <- diversity(dune[gr == "2", ], "simp")
c2m(pop1=simp.one,pop2=simp.two,nrandom=1999,verbose=FALSE)
#done.
#$res
#
#mv1                 8.630e-01
#mv2                 8.773e-01
#mv1-mv2            -1.439e-02
#p                   2.440e-01 # N.S.
#quantile 0.025     -3.456e-02
#quantile 0.975      3.351e-02
#randomized mv1-mv2  3.899e-04
#nrandom             1.999e+03
#########################

The possible pitfalls might be hidden under the different results ;-)

Cheers,
Ivailo
-- 
UBUNTU: a person is a person through other persons.