################################################################################## ############################################################### ############################################################### #### INSTALL AND LOAD PACKAGES ############################################################### ############################################################### ################################################################################## install.packages("rareNMtests", dep=T) library(rareNMtests) ### Any of the following analysis will take several weeks to complete ### Note that the original analyses were run at the Bioportal server, a web-based portal ### for data analysis that allows for parallel processing. The Bioportal is now closed. ################################################################################## ################################################################################## ################################################################################## ###### SCENARIOS 1 TO 8 FOR INDIVIDUAL-BASED RAREFACTION CURVES ################################################################################## ################################################################################## ################################################################################## ################################################################################## ################################################################################## ###### SCENARIO 1: DRAWING TWO RANDOM SAMPLES FROM THE SAME METACOMMUNITY ################################################################################## ################################################################################## n <- 750 niter <- 200 simres1 <- data.frame(spn = NULL, Nind = NULL, sdlog = NULL, Ss = NULL, p.ecoS0 = NULL, p.ecoS1 = NULL, p.ecoS2 = NULL, p.ecoC0 = NULL, p.ecoC1 = NULL, p.ecoC2 = NULL, p.bioS0L = NULL, p.bioS1L = NULL, p.bioS2L = NULL, p.bioC0L = NULL, p.bioC1L = NULL, p.bioC2L = NULL, p.bioS0G = NULL, p.bioS1G = NULL, p.bioS2G = NULL, p.bioC0G = NULL, p.bioC1G = NULL, p.bioC2G = NULL, p.bioS0S = NULL, p.bioS1S = NULL, p.bioS2S = NULL, p.bioC0S = NULL, p.bioC1S = NULL, p.bioC2S = NULL) for (i in 1:n) { ### Simulate metacommunity spn <- round(runif(1, 10, 200)) sdlog <- runif(1, log(1.1), log(33)) com <- round(rlnorm(spn, 2, sdlog)) ### Simulate random sampling Ss <- round(runif(1, 50, 1000)) sample1 <- sample(paste("sp",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) sample1 <- data.frame(table(sample1)) colnames(sample1) <- c("species", "sample1") sample2 <- sample(paste("sp",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", "sample2") df <- merge(sample1, sample2, by="species", all=TRUE) rownames(df) <- df$species df[is.na(df)] <- 0 df <- df[,2:3] ### Apply ecological null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers try(ecoS0 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS1 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS2 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC0 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC1 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC2 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) ### Apply biogeographical null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers ### Using the log-normal for the distribution of random assemblages try(bioS0L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS1L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS2L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC0L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC1L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC2L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) ### Using the geometric series for the distribution of random assemblages try(bioS0G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS1G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS2G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC0G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC1G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC2G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) ### Using the brocken-stick for the distribution of random assemblages try(bioS0S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS1S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS2S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC0S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC1S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC2S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(simres <- data.frame(spn, Nind = sum(com), sdlog, Ss, p.ecoS0 =ecoS0$pval[2], p.ecoS1 =ecoS1$pval[2], p.ecoS2 =ecoS2$pval[2], p.ecoC0 =ecoC0$pval[2], p.ecoC1 =ecoC1$pval[2], p.ecoC2 =ecoC2$pval[2], p.bioS0L =bioS0L$pval[2], p.bioS1L =bioS1L$pval[2], p.bioS2L =bioS2L$pval[2], p.bioC0L =bioC0L$pval[2], p.bioC1L =bioC1L$pval[2], p.bioC2L =bioC2L$pval[2], p.bioS0G =bioS0G$pval[2], p.bioS1G =bioS1G$pval[2], p.bioS2G =bioS2G$pval[2], p.bioC0G =bioC0G$pval[2], p.bioC1G =bioC1G$pval[2], p.bioC2G =bioC2G$pval[2], p.bioS0S =bioS0S$pval[2], p.bioS1S =bioS1S$pval[2], p.bioS2S =bioS2S$pval[2], p.bioC0S =bioC0S$pval[2], p.bioC1S =bioC1S$pval[2], p.bioC2S =bioC2S$pval[2])) try(simres1 <- rbind(simres1, simres)) } ################################################################################## ################################################################################## ###### SCENARIO 2: DRAWING A RANDOM NUMBER OF SAMPLES FROM THE SAME METACOMMUNITY ################################################################################## ################################################################################## n <- 250 niter <- 200 simres2 <- data.frame(spn = NULL, Nind = NULL, sdlog = NULL, nplots = NULL, Ss = NULL, p.ecoS0 = NULL, p.ecoS1 = NULL, p.ecoS2 = NULL, p.ecoC0 = NULL, p.ecoC1 = NULL, p.ecoC2 = NULL, p.bioS0L = NULL, p.bioS1L = NULL, p.bioS2L = NULL, p.bioC0L = NULL, p.bioC1L = NULL, p.bioC2L = NULL, p.bioS0G = NULL, p.bioS1G = NULL, p.bioS2G = NULL, p.bioC0G = NULL, p.bioC1G = NULL, p.bioC2G = NULL, p.bioS0S = NULL, p.bioS1S = NULL, p.bioS2S = NULL, p.bioC0S = NULL, p.bioC1S = NULL, p.bioC2S = NULL) for (i in 1:n) { ### Simulate metacommunity spn <- 100 nplots <- round(runif(1, 2, 15)) sdlog <- log(11.6) com <- round(rlnorm(spn, 2, sdlog)) ### Simulate random sampling Ss <- 500 samples <- sample(paste("sp",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) samples <- data.frame(table(samples)) colnames(samples) <- c("species", "sample1") for(t in 2:(nplots)) { sample2 <- sample(paste("sp",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", paste("sample", t, sep="")) samples <- merge(samples, sample2, by="species", all=TRUE) } rownames(samples) <- samples$species samples[is.na(samples)] <- 0 df <- samples[,-1] ### Apply ecological null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers try(ecoS0 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS1 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS2 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC0 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC1 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC2 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) ### Apply biogeographical null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers ### Using the log-normal for the distribution of random assemblages try(bioS0L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS1L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS2L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC0L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC1L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC2L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) ### Using the geometric series for the distribution of random assemblages try(bioS0G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS1G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS2G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC0G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC1G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC2G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) ### Using the brocken-stick for the distribution of random assemblages try(bioS0S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS1S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS2S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC0S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC1S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC2S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(simres <- data.frame(spn, Nind = sum(com), sdlog, nplots, Ss, p.ecoS0 =ecoS0$pval[2], p.ecoS1 =ecoS1$pval[2], p.ecoS2 =ecoS2$pval[2], p.ecoC0 =ecoC0$pval[2], p.ecoC1 =ecoC1$pval[2], p.ecoC2 =ecoC2$pval[2], p.bioS0L =bioS0L$pval[2], p.bioS1L =bioS1L$pval[2], p.bioS2L =bioS2L$pval[2], p.bioC0L =bioC0L$pval[2], p.bioC1L =bioC1L$pval[2], p.bioC2L =bioC2L$pval[2], p.bioS0G =bioS0G$pval[2], p.bioS1G =bioS1G$pval[2], p.bioS2G =bioS2G$pval[2], p.bioC0G =bioC0G$pval[2], p.bioC1G =bioC1G$pval[2], p.bioC2G =bioC2G$pval[2], p.bioS0S =bioS0S$pval[2], p.bioS1S =bioS1S$pval[2], p.bioS2S =bioS2S$pval[2], p.bioC0S =bioC0S$pval[2], p.bioC1S =bioC1S$pval[2], p.bioC2S =bioC2S$pval[2])) try(simres2 <- rbind(simres2, simres)) } ################################################################################## ################################################################################## ###### SCENARIO 3: DRAWING TWO SAMPLES FROM COMMUNITIES WITH DIFFERENT COMPOSITION ###### BUT SIMILAR DISTRIBUTION OF RELATIVE ABUNDANCES AND SPECIES RICHNESS ################################################################################## ################################################################################## n <- 750 niter <- 200 simres3 <- data.frame(spn = NULL, Nind = NULL, sdlog = NULL, Ss = NULL, p.ecoS0 = NULL, p.ecoS1 = NULL, p.ecoS2 = NULL, p.ecoC0 = NULL, p.ecoC1 = NULL, p.ecoC2 = NULL, p.bioS0L = NULL, p.bioS1L = NULL, p.bioS2L = NULL, p.bioC0L = NULL, p.bioC1L = NULL, p.bioC2L = NULL, p.bioS0G = NULL, p.bioS1G = NULL, p.bioS2G = NULL, p.bioC0G = NULL, p.bioC1G = NULL, p.bioC2G = NULL, p.bioS0S = NULL, p.bioS1S = NULL, p.bioS2S = NULL, p.bioC0S = NULL, p.bioC1S = NULL, p.bioC2S = NULL) for (i in 1:n) { ### Simulate metacommunity spn <- round(runif(1, 10, 200)) sdlog <- runif(1, log(1.1), log(33)) com <- round(rlnorm(spn, 3, sdlog)) ### Simulate random sampling Ss <- round(runif(1, 50, 1000)) sample1 <- sample(paste("sp",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) sample1 <- data.frame(table(sample1)) colnames(sample1) <- c("species", "sample1") sample2 <- sample(paste("sp",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", "sample2") sample2$species <- paste("spx", c(1:length(rownames(sample2)))) df <- merge(sample1, sample2, by="species", all=TRUE) rownames(df) <- df$species df[is.na(df)] <- 0 df <- df[,2:3] ### Apply ecological null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers try(ecoS0 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS1 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS2 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC0 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC1 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC2 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) ### Apply biogeographical null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers ### Using the log-normal for the distribution of random assemblages try(bioS0L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS1L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS2L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC0L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC1L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC2L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) ### Using the geometric series for the distribution of random assemblages try(bioS0G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS1G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS2G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC0G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC1G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC2G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) ### Using the brocken-stick for the distribution of random assemblages try(bioS0S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS1S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS2S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC0S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC1S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC2S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(simres <- data.frame(spn, Nind = sum(com), sdlog, Ss, p.ecoS0 =ecoS0$pval[2], p.ecoS1 =ecoS1$pval[2], p.ecoS2 =ecoS2$pval[2], p.ecoC0 =ecoC0$pval[2], p.ecoC1 =ecoC1$pval[2], p.ecoC2 =ecoC2$pval[2], p.bioS0L =bioS0L$pval[2], p.bioS1L =bioS1L$pval[2], p.bioS2L =bioS2L$pval[2], p.bioC0L =bioC0L$pval[2], p.bioC1L =bioC1L$pval[2], p.bioC2L =bioC2L$pval[2], p.bioS0G =bioS0G$pval[2], p.bioS1G =bioS1G$pval[2], p.bioS2G =bioS2G$pval[2], p.bioC0G =bioC0G$pval[2], p.bioC1G =bioC1G$pval[2], p.bioC2G =bioC2G$pval[2], p.bioS0S =bioS0S$pval[2], p.bioS1S =bioS1S$pval[2], p.bioS2S =bioS2S$pval[2], p.bioC0S =bioC0S$pval[2], p.bioC1S =bioC1S$pval[2], p.bioC2S =bioC2S$pval[2])) try(simres3 <- rbind(simres3, simres)) } ################################################################################## ################################################################################## ###### SCENARIO 4: DRAWING A RANDOM NUMBER OF SAMPLES FROM TWO METACOMMUNITIES ###### WITH DIFFERENT SPECIES COMPOSITION BUT SIMILAR RICHNESS AND DISTRIBUTION ###### OF RELATIVE ABUNDANCES ################################################################################## ################################################################################## n <- 250 niter <- 200 simres4 <- data.frame(spn = NULL, Nind = NULL, sdlog = NULL, nplots = NULL, ngrs1 = NULL, ngrs2 = NULL, Ss = NULL, p.ecoS0 = NULL, p.ecoS1 = NULL, p.ecoS2 = NULL, p.ecoC0 = NULL, p.ecoC1 = NULL, p.ecoC2 = NULL, p.bioS0L = NULL, p.bioS1L = NULL, p.bioS2L = NULL, p.bioC0L = NULL, p.bioC1L = NULL, p.bioC2L = NULL, p.bioS0G = NULL, p.bioS1G = NULL, p.bioS2G = NULL, p.bioC0G = NULL, p.bioC1G = NULL, p.bioC2G = NULL, p.bioS0S = NULL, p.bioS1S = NULL, p.bioS2S = NULL, p.bioC0S = NULL, p.bioC1S = NULL, p.bioC2S = NULL) for (i in 1:n) { ### Simulate metacommunity spn <- 100 nplots <- round(runif(1, 2, 15)) ngrs1 <- round(runif(1, 1, max(nplots)-1)) ngrs2 <- nplots-ngrs1 sdlog <- log(11.6) com <- round(rlnorm(spn, 6, sdlog)) ### Simulate random sampling Ss <- 500 samples.a <- sample(paste("sp",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) samples.a <- data.frame(table(samples.a)) colnames(samples.a) <- c("species", "sample1") for(t in 2:(ngrs1)) { sample2 <- sample(paste("sp",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", paste("sample", t, sep="")) samples.a <- merge(samples.a, sample2, by="species", all=TRUE) } samples.b <- sample(paste("spx",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) samples.b <- data.frame(table(samples.b)) colnames(samples.b) <- c("species", paste("sample", ngrs1+1, sep="")) for(t in 2:(ngrs2)) { sample2 <- sample(paste("spx",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", paste("sample", ngrs1+t, sep="")) samples.b <- merge(samples.b, sample2, by="species", all=TRUE) } df <- merge(samples.a, samples.b, by="species", all=TRUE) rownames(df) <- df$species df[is.na(df)] <- 0 df <- df[,-1] ### Apply ecological null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers try(ecoS0 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS1 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS2 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC0 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC1 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC2 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) ### Apply biogeographical null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers ### Using the log-normal for the distribution of random assemblages try(bioS0L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS1L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS2L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC0L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC1L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC2L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) ### Using the geometric series for the distribution of random assemblages try(bioS0G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS1G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS2G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC0G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC1G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC2G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) ### Using the brocken-stick for the distribution of random assemblages try(bioS0S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS1S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS2S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC0S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC1S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC2S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(simres <- data.frame(spn, Nind = sum(com), sdlog, nplots, ngrs1, ngrs2, Ss, p.ecoS0 =ecoS0$pval[2], p.ecoS1 =ecoS1$pval[2], p.ecoS2 =ecoS2$pval[2], p.ecoC0 =ecoC0$pval[2], p.ecoC1 =ecoC1$pval[2], p.ecoC2 =ecoC2$pval[2], p.bioS0L =bioS0L$pval[2], p.bioS1L =bioS1L$pval[2], p.bioS2L =bioS2L$pval[2], p.bioC0L =bioC0L$pval[2], p.bioC1L =bioC1L$pval[2], p.bioC2L =bioC2L$pval[2], p.bioS0G =bioS0G$pval[2], p.bioS1G =bioS1G$pval[2], p.bioS2G =bioS2G$pval[2], p.bioC0G =bioC0G$pval[2], p.bioC1G =bioC1G$pval[2], p.bioC2G =bioC2G$pval[2], p.bioS0S =bioS0S$pval[2], p.bioS1S =bioS1S$pval[2], p.bioS2S =bioS2S$pval[2], p.bioC0S =bioC0S$pval[2], p.bioC1S =bioC1S$pval[2], p.bioC2S =bioC2S$pval[2])) try(simres4 <- rbind(simres4, simres)) } ################################################################################## ################################################################################## ###### SCENARIO 5: DRAWING TWO SAMPLES FROM TWO METACOMMUNITIES WITH DIFFERENT ###### SPECIES RICHNESS BUT SIMILAR DISTRIBUTION OF RELATIVE ABUNDANCES ################################################################################## ################################################################################## n <- 750 niter <- 200 simres5 <- data.frame(Nind1 = NULL, Nind2 = NULL, sdlog = NULL, Ss = NULL, p.ecoS0 = NULL, p.ecoS1 = NULL, p.ecoS2 = NULL, p.ecoC0 = NULL, p.ecoC1 = NULL, p.ecoC2 = NULL, p.bioS0L = NULL, p.bioS1L = NULL, p.bioS2L = NULL, p.bioC0L = NULL, p.bioC1L = NULL, p.bioC2L = NULL, p.bioS0G = NULL, p.bioS1G = NULL, p.bioS2G = NULL, p.bioC0G = NULL, p.bioC1G = NULL, p.bioC2G = NULL, p.bioS0S = NULL, p.bioS1S = NULL, p.bioS2S = NULL, p.bioC0S = NULL, p.bioC1S = NULL, p.bioC2S = NULL) for (i in 1:n) { ### Simulate metacommunities spn1 <- 50 spn2 <- 150 sdlog <- runif(1, log(1.1), log(33)) com1 <- round(rlnorm(spn1, 6, sdlog)) com2 <- round(rlnorm(spn2, 6, sdlog)) ### Simulate random sampling Ss <- round(runif(1, 50, 1000)) sample1 <- sample(paste("sp",1:length(com1)), rpois(1, Ss), replace=TRUE, prob=com1) sample1 <- data.frame(table(sample1)) colnames(sample1) <- c("species", "sample1") sample2 <- sample(paste("sp",1:length(com2)), rpois(1, Ss), replace=TRUE, prob=com2) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", "sample2") df <- merge(sample1, sample2, by="species", all=TRUE) rownames(df) <- df$species df[is.na(df)] <- 0 df <- df[,2:3] ### Apply ecological null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers try(ecoS0 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS1 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS2 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC0 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC1 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC2 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) ### Apply biogeographical null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers ### Using the log-normal for the distribution of random assemblages try(bioS0L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS1L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS2L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC0L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC1L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC2L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) ### Using the geometric series for the distribution of random assemblages try(bioS0G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS1G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS2G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC0G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC1G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC2G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) ### Using the brocken-stick for the distribution of random assemblages try(bioS0S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS1S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS2S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC0S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC1S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC2S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(simres <- data.frame(Nind1 = sum(com1), Nind2 = sum(com2), sdlog, Ss, p.ecoS0 =ecoS0$pval[2], p.ecoS1 =ecoS1$pval[2], p.ecoS2 =ecoS2$pval[2], p.ecoC0 =ecoC0$pval[2], p.ecoC1 =ecoC1$pval[2], p.ecoC2 =ecoC2$pval[2], p.bioS0L =bioS0L$pval[2], p.bioS1L =bioS1L$pval[2], p.bioS2L =bioS2L$pval[2], p.bioC0L =bioC0L$pval[2], p.bioC1L =bioC1L$pval[2], p.bioC2L =bioC2L$pval[2], p.bioS0G =bioS0G$pval[2], p.bioS1G =bioS1G$pval[2], p.bioS2G =bioS2G$pval[2], p.bioC0G =bioC0G$pval[2], p.bioC1G =bioC1G$pval[2], p.bioC2G =bioC2G$pval[2], p.bioS0S =bioS0S$pval[2], p.bioS1S =bioS1S$pval[2], p.bioS2S =bioS2S$pval[2], p.bioC0S =bioC0S$pval[2], p.bioC1S =bioC1S$pval[2], p.bioC2S =bioC2S$pval[2])) try(simres5 <- rbind(simres5, simres)) } ################################################################################## ################################################################################## ###### SCENARIO 6: DRAWING A RANDOM NUMBER OF SAMPLES FROM TWO METACOMMUNITIES ###### WITH DIFFERENT SPECIES RICHNESS BUT SIMILAR DISTRIBUTION OF RELATIVE ABUNDANCES ################################################################################## ################################################################################## n <- 250 niter <- 200 simres6 <- data.frame(nplots = NULL, ngrs1 = NULL, ngrs2 = NULL, Ss = NULL, p.ecoS0 = NULL, p.ecoS1 = NULL, p.ecoS2 = NULL, p.ecoC0 = NULL, p.ecoC1 = NULL, p.ecoC2 = NULL, p.bioS0L = NULL, p.bioS1L = NULL, p.bioS2L = NULL, p.bioC0L = NULL, p.bioC1L = NULL, p.bioC2L = NULL, p.bioS0G = NULL, p.bioS1G = NULL, p.bioS2G = NULL, p.bioC0G = NULL, p.bioC1G = NULL, p.bioC2G = NULL, p.bioS0S = NULL, p.bioS1S = NULL, p.bioS2S = NULL, p.bioC0S = NULL, p.bioC1S = NULL, p.bioC2S = NULL) for (i in 1:n) { ### Simulate metacommunities spn1 <- 50 spn2 <- 150 nplots <- round(runif(1, 2, 15)) ngrs1 <- round(runif(1, 1, max(nplots)-1)) ngrs2 <- nplots-ngrs1 sdlog <- log(11.6) com1 <- round(rlnorm(spn1, 2, sdlog)) com2 <- round(rlnorm(spn2, 2, sdlog)) ### Simulate random sampling Ss <- 500 samples.a <- sample(paste("sp",1:length(com1)), rpois(1, Ss), replace=TRUE, prob=com1) samples.a <- data.frame(table(samples.a)) colnames(samples.a) <- c("species", "sample1") for(t in 2:(ngrs1)) { sample2 <- sample(paste("sp",1:length(com1)), rpois(1, Ss), replace=TRUE, prob=com1) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", paste("sample", t, sep="")) samples.a <- merge(samples.a, sample2, by="species", all=TRUE) } samples.b <- sample(paste("sp",1:length(com2)), rpois(1, Ss), replace=TRUE, prob=com2) samples.b <- data.frame(table(samples.b)) colnames(samples.b) <- c("species", paste("sample", ngrs1+1, sep="")) for(t in 2:(ngrs2)) { sample2 <- sample(paste("sp",1:length(com2)), rpois(1, Ss), replace=TRUE, prob=com2) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", paste("sample", ngrs1+t, sep="")) samples.b <- merge(samples.b, sample2, by="species", all=TRUE) } df <- merge(samples.a, samples.b, by="species", all=TRUE) rownames(df) <- df$species df[is.na(df)] <- 0 df <- df[,-1] ### Apply ecological null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers try(ecoS0 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS1 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS2 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC0 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC1 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC2 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) ### Apply biogeographical null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers ### Using the log-normal for the distribution of random assemblages try(bioS0L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS1L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS2L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC0L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC1L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC2L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) ### Using the geometric series for the distribution of random assemblages try(bioS0G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS1G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS2G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC0G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC1G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC2G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) ### Using the brocken-stick for the distribution of random assemblages try(bioS0S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS1S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS2S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC0S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC1S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC2S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(simres <- data.frame(nplots, ngrs1, ngrs2, Ss, p.ecoS0 =ecoS0$pval[2], p.ecoS1 =ecoS1$pval[2], p.ecoS2 =ecoS2$pval[2], p.ecoC0 =ecoC0$pval[2], p.ecoC1 =ecoC1$pval[2], p.ecoC2 =ecoC2$pval[2], p.bioS0L =bioS0L$pval[2], p.bioS1L =bioS1L$pval[2], p.bioS2L =bioS2L$pval[2], p.bioC0L =bioC0L$pval[2], p.bioC1L =bioC1L$pval[2], p.bioC2L =bioC2L$pval[2], p.bioS0G =bioS0G$pval[2], p.bioS1G =bioS1G$pval[2], p.bioS2G =bioS2G$pval[2], p.bioC0G =bioC0G$pval[2], p.bioC1G =bioC1G$pval[2], p.bioC2G =bioC2G$pval[2], p.bioS0S =bioS0S$pval[2], p.bioS1S =bioS1S$pval[2], p.bioS2S =bioS2S$pval[2], p.bioC0S =bioC0S$pval[2], p.bioC1S =bioC1S$pval[2], p.bioC2S =bioC2S$pval[2])) try(simres6 <- rbind(simres6, simres)) } ################################################################################## ################################################################################## ###### SCENARIO 7: DRAWING TWO SAMPLES FROM TWO METACOMMUNITIES WITH DIFFERENT ###### DISTRIBUTION OF RELATIVE ABUNDANCES BUT SIMILAR RICHNESS AND COMPOSITION ################################################################################## ################################################################################## n <- 750 niter <- 200 simres7 <- data.frame(spn = NULL, sdlog1 = NULL, sdlog2 = NULL,Ss = NULL, p.ecoS0 = NULL, p.ecoS1 = NULL, p.ecoS2 = NULL, p.ecoC0 = NULL, p.ecoC1 = NULL, p.ecoC2 = NULL, p.bioS0L = NULL, p.bioS1L = NULL, p.bioS2L = NULL, p.bioC0L = NULL, p.bioC1L = NULL, p.bioC2L = NULL, p.bioS0G = NULL, p.bioS1G = NULL, p.bioS2G = NULL, p.bioC0G = NULL, p.bioC1G = NULL, p.bioC2G = NULL, p.bioS0S = NULL, p.bioS1S = NULL, p.bioS2S = NULL, p.bioC0S = NULL, p.bioC1S = NULL, p.bioC2S = NULL) for (i in 1:n) { ### Simulate metacommunities spn <- round(runif(1, 10, 200)) sdlog1 <- runif(1, log(1.1), log(33)) sdlog2 <- ifelse(sdlog1 < 1.5, sdlog1 + 1, sdlog1 - 1) com1 <- round(rlnorm(spn, 3, sdlog1)) com2 <- round(rlnorm(spn, 3, sdlog2)) ### Simulate random sampling Ss <- round(runif(1, 50, 1000)) sample1 <- sample(paste("sp",1:length(com1)), rpois(1, Ss), replace=TRUE, prob=com1) sample1 <- data.frame(table(sample1)) colnames(sample1) <- c("species", "sample1") sample2 <- sample(paste("sp",1:length(com2)), rpois(1, Ss), replace=TRUE, prob=com2) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", "sample2") df <- merge(sample1, sample2, by="species", all=TRUE) rownames(df) <- df$species df[is.na(df)] <- 0 df <- df[,2:3] ### Apply ecological null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers try(ecoS0 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS1 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS2 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC0 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC1 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC2 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) ### Apply biogeographical null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers ### Using the log-normal for the distribution of random assemblages try(bioS0L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS1L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS2L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC0L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC1L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC2L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) ### Using the geometric series for the distribution of random assemblages try(bioS0G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS1G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS2G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC0G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC1G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC2G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) ### Using the brocken-stick for the distribution of random assemblages try(bioS0S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS1S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS2S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC0S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC1S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC2S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(simres <- data.frame(spn, sdlog1, sdlog2, Ss, p.ecoS0 =ecoS0$pval[2], p.ecoS1 =ecoS1$pval[2], p.ecoS2 =ecoS2$pval[2], p.ecoC0 =ecoC0$pval[2], p.ecoC1 =ecoC1$pval[2], p.ecoC2 =ecoC2$pval[2], p.bioS0L =bioS0L$pval[2], p.bioS1L =bioS1L$pval[2], p.bioS2L =bioS2L$pval[2], p.bioC0L =bioC0L$pval[2], p.bioC1L =bioC1L$pval[2], p.bioC2L =bioC2L$pval[2], p.bioS0G =bioS0G$pval[2], p.bioS1G =bioS1G$pval[2], p.bioS2G =bioS2G$pval[2], p.bioC0G =bioC0G$pval[2], p.bioC1G =bioC1G$pval[2], p.bioC2G =bioC2G$pval[2], p.bioS0S =bioS0S$pval[2], p.bioS1S =bioS1S$pval[2], p.bioS2S =bioS2S$pval[2], p.bioC0S =bioC0S$pval[2], p.bioC1S =bioC1S$pval[2], p.bioC2S =bioC2S$pval[2])) try(simres7 <- rbind(simres7, simres)) } ################################################################################## ################################################################################## ###### SCENARIO 8: DRAWING A RANDOM NUMBER OF SAMPLES FROM TWO METACOMMUNITIES ###### WITH DIFFERENT DISTRIBUTION OF RELATIVE ABUNDANCES BUT SIMILAR RICHNESS ###### AND COMPOSITION ################################################################################## ################################################################################## n <- 250 niter <- 200 simres8 <- data.frame(nplots = NULL, ngrs1 = NULL, ngrs2 = NULL, sdlog1 = NULL, sdlog2 = NULL, p.ecoS0 = NULL, p.ecoS1 = NULL, p.ecoS2 = NULL, p.ecoC0 = NULL, p.ecoC1 = NULL, p.ecoC2 = NULL, p.bioS0L = NULL, p.bioS1L = NULL, p.bioS2L = NULL, p.bioC0L = NULL, p.bioC1L = NULL, p.bioC2L = NULL, p.bioS0G = NULL, p.bioS1G = NULL, p.bioS2G = NULL, p.bioC0G = NULL, p.bioC1G = NULL, p.bioC2G = NULL, p.bioS0S = NULL, p.bioS1S = NULL, p.bioS2S = NULL, p.bioC0S = NULL, p.bioC1S = NULL, p.bioC2S = NULL) for (i in 1:n) { ### Simulate metacommunities spn <- 100 nplots <- round(runif(1, 2, 15)) ngrs1 <- round(runif(1, 1, max(nplots)-1)) ngrs2 <- nplots-ngrs1 sdlog1 <- runif(1, log(1.1), log(33)) sdlog2 <- ifelse(sdlog1 < 1.5, sdlog1 + 1, sdlog1 - 1) com1 <- round(rlnorm(spn, 2, sdlog1)) com2 <- round(rlnorm(spn, 2, sdlog2)) ### Simulate random sampling Ss <- 500 samples.a <- sample(paste("sp",1:length(com1)), rpois(1, Ss), replace=TRUE, prob=com1) samples.a <- data.frame(table(samples.a)) colnames(samples.a) <- c("species", "sample1") for(t in 2:(ngrs1)) { sample2 <- sample(paste("sp",1:length(com1)), rpois(1, Ss), replace=TRUE, prob=com1) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", paste("sample", t, sep="")) samples.a <- merge(samples.a, sample2, by="species", all=TRUE) } samples.b <- sample(paste("sp",1:length(com2)), rpois(1, Ss), replace=TRUE, prob=com2) samples.b <- data.frame(table(samples.b)) colnames(samples.b) <- c("species", paste("sample", ngrs1+1, sep="")) for(t in 2:(ngrs2)) { sample2 <- sample(paste("sp",1:length(com2)), rpois(1, Ss), replace=TRUE, prob=com2) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", paste("sample", ngrs1+t, sep="")) samples.b <- merge(samples.b, sample2, by="species", all=TRUE) } df <- merge(samples.a, samples.b, by="species", all=TRUE) rownames(df) <- df$species df[is.na(df)] <- 0 df <- df[,-1] ### Apply ecological null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers try(ecoS0 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS1 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoS2 <- EcoTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC0 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC1 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) try(ecoC2 <- EcoTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE)) ### Apply biogeographical null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers ### Using the log-normal for the distribution of random assemblages try(bioS0L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS1L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioS2L <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC0L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC1L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) try(bioC2L <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="lnorm")) ### Using the geometric series for the distribution of random assemblages try(bioS0G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS1G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioS2G <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC0G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC1G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) try(bioC2G <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="geom")) ### Using the brocken-stick for the distribution of random assemblages try(bioS0S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS1S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioS2S <- BiogTest.individual(df, niter=niter, method = "sample-size", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC0S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 0, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC1S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 1, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(bioC2S <- BiogTest.individual(df, niter=niter, method = "coverage", q = 2, powerfun = 0.9, log.scale = TRUE, MARGIN=2, trace=FALSE, distr="stick")) try(simres <- data.frame(nplots, ngrs1, ngrs2, sdlog1, sdlog2, p.ecoS0 =ecoS0$pval[2], p.ecoS1 =ecoS1$pval[2], p.ecoS2 =ecoS2$pval[2], p.ecoC0 =ecoC0$pval[2], p.ecoC1 =ecoC1$pval[2], p.ecoC2 =ecoC2$pval[2], p.bioS0L =bioS0L$pval[2], p.bioS1L =bioS1L$pval[2], p.bioS2L =bioS2L$pval[2], p.bioC0L =bioC0L$pval[2], p.bioC1L =bioC1L$pval[2], p.bioC2L =bioC2L$pval[2], p.bioS0G =bioS0G$pval[2], p.bioS1G =bioS1G$pval[2], p.bioS2G =bioS2G$pval[2], p.bioC0G =bioC0G$pval[2], p.bioC1G =bioC1G$pval[2], p.bioC2G =bioC2G$pval[2], p.bioS0S =bioS0S$pval[2], p.bioS1S =bioS1S$pval[2], p.bioS2S =bioS2S$pval[2], p.bioC0S =bioC0S$pval[2], p.bioC1S =bioC1S$pval[2], p.bioC2S =bioC2S$pval[2])) try(simres8 <- rbind(simres8, simres)) } ################################################################################## ################################################################################## ################################################################################## ###### SCENARIOS 9 TO 12 FOR SAMPLE-BASED RAREFACTION CURVES ###### ONLY TWO SAMPLE-BASED RAREFACTION CURVES ARE COMPARED AT A TIME ################################################################################## ################################################################################## ################################################################################## ################################################################################## ################################################################################## ###### SCENARIO 9: DRAWING TWO RANDOM SET OF SAMPLES FROM THE SAME METACOMMUNITY ################################################################################## ################################################################################## n <- 750 niter <- 200 simres9 <- data.frame(nplots = NULL, ngrs1 = NULL, ngrs2 = NULL, sdlog = NULL, p.ecoS0 = NULL, p.ecoS1 = NULL, p.ecoS2 = NULL, p.ecoC0 = NULL, p.ecoC1 = NULL, p.ecoC2 = NULL, p.bioS0L = NULL, p.bioS1L = NULL, p.bioS2L = NULL, p.bioC0L = NULL, p.bioC1L = NULL, p.bioC2L = NULL, p.bioS0G = NULL, p.bioS1G = NULL, p.bioS2G = NULL, p.bioC0G = NULL, p.bioC1G = NULL, p.bioC2G = NULL, p.bioS0S = NULL, p.bioS1S = NULL, p.bioS2S = NULL, p.bioC0S = NULL, p.bioC1S = NULL, p.bioC2S = NULL) for (i in 1:n) { ### Simulate metacommunities spn <- 100 nplots <- round(runif(1, 10, 50)) ngrs1 <- round(runif(1, 5, max(nplots)-5)) ngrs2 <- nplots-ngrs1 sdlog <- runif(1, log(1.1), log(33)) com <- round(rlnorm(spn, 3, sdlog)) ### Simulate random sampling Ss <- 50 samples <- sample(paste("sp",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) samples <- data.frame(table(samples)) colnames(samples) <- c("species", "sample1") for(t in 2:(nplots)) { sample2 <- sample(paste("sp",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", paste("sample", t, sep="")) samples <- merge(samples, sample2, by="species", all=TRUE) } rownames(samples) <- samples$species samples[is.na(samples)] <- 0 df <- samples[,-1] ### Apply ecological null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers try(ecoS0 <- EcoTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoS1 <- EcoTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoS2 <- EcoTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoC0 <- EcoTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoC1 <- EcoTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoC2 <- EcoTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) ### Apply biogeographical null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers ### Using the log-normal for the distribution of random assemblages try(bioS0L <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS1L <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS2L <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC0L <- BiogTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC1L <- BiogTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC2L <- BiogTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) ### Using the geometric series for the distribution of random assemblages try(bioS0G <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS1G <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS2G <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC0G <- BiogTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC1G <- BiogTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC2G <- BiogTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) ### Using the brocken-stick for the distribution of random assemblages try(bioS0S <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS1S <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS2S <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC0S <- BiogTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC1S <- BiogTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC2S <- BiogTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(simres <- data.frame(nplots, ngrs1, ngrs2, sdlog, p.ecoS0 =ecoS0$pval[2], p.ecoS1 =ecoS1$pval[2], p.ecoS2 =ecoS2$pval[2], p.ecoC0 =ecoC0$pval[2], p.ecoC1 =ecoC1$pval[2], p.ecoC2 =ecoC2$pval[2], p.bioS0L =bioS0L$pval[2], p.bioS1L =bioS1L$pval[2], p.bioS2L =bioS2L$pval[2], p.bioC0L =bioC0L$pval[2], p.bioC1L =bioC1L$pval[2], p.bioC2L =bioC2L$pval[2], p.bioS0G =bioS0G$pval[2], p.bioS1G =bioS1G$pval[2], p.bioS2G =bioS2G$pval[2], p.bioC0G =bioC0G$pval[2], p.bioC1G =bioC1G$pval[2], p.bioC2G =bioC2G$pval[2], p.bioS0S =bioS0S$pval[2], p.bioS1S =bioS1S$pval[2], p.bioS2S =bioS2S$pval[2], p.bioC0S =bioC0S$pval[2], p.bioC1S =bioC1S$pval[2], p.bioC2S =bioC2S$pval[2])) try(simres9 <- rbind(simres9, simres)) } ################################################################################## ################################################################################## ###### SCENARIO 10: DRAWING TWO SET OF RANDOM SAMPLES FROM COMMUNITIES WITH ###### DIFFERENT COMPOSITION BUT SIMILAR DISTRIBUTION OF RELATIVE ABUNDANCES AND ###### SPECIES RICHNESS ################################################################################## ################################################################################## n <- 750 niter <- 200 simres10 <- data.frame(nplots = NULL, ngrs1 = NULL, ngrs2 = NULL, sdlog = NULL, p.ecoS0 = NULL, p.ecoS1 = NULL, p.ecoS2 = NULL, p.ecoC0 = NULL, p.ecoC1 = NULL, p.ecoC2 = NULL, p.bioS0L = NULL, p.bioS1L = NULL, p.bioS2L = NULL, p.bioC0L = NULL, p.bioC1L = NULL, p.bioC2L = NULL, p.bioS0G = NULL, p.bioS1G = NULL, p.bioS2G = NULL, p.bioC0G = NULL, p.bioC1G = NULL, p.bioC2G = NULL, p.bioS0S = NULL, p.bioS1S = NULL, p.bioS2S = NULL, p.bioC0S = NULL, p.bioC1S = NULL, p.bioC2S = NULL) for (i in 1:n) { ### Simulate metacommunities spn <- 100 nplots <- round(runif(1, 10, 50)) ngrs1 <- round(runif(1, 5, max(nplots)-5)) ngrs2 <- nplots-ngrs1 sdlog <- runif(1, log(1.1), log(33)) com <- round(rlnorm(spn, 2, sdlog)) ### Simulate random sampling Ss <- 50 samples.a <- sample(paste("sp",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) samples.a <- data.frame(table(samples.a)) colnames(samples.a) <- c("species", "sample1") for(t in 2:(ngrs1)) { sample2 <- sample(paste("sp",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", paste("sample", t, sep="")) samples.a <- merge(samples.a, sample2, by="species", all=TRUE) } samples.b <- sample(paste("spx",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) samples.b <- data.frame(table(samples.b)) colnames(samples.b) <- c("species", paste("sample", ngrs1+1, sep="")) for(t in 2:(ngrs2)) { sample2 <- sample(paste("spx",1:length(com)), rpois(1, Ss), replace=TRUE, prob=com) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", paste("sample", ngrs1+t, sep="")) samples.b <- merge(samples.b, sample2, by="species", all=TRUE) } df <- merge(samples.a, samples.b, by="species", all=TRUE) rownames(df) <- df$species df[is.na(df)] <- 0 df <- df[,-1] ### Apply ecological null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers try(ecoS0 <- EcoTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoS1 <- EcoTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoS2 <- EcoTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoC0 <- EcoTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoC1 <- EcoTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoC2 <- EcoTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) ### Apply biogeographical null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers ### Using the log-normal for the distribution of random assemblages try(bioS0L <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS1L <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS2L <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC0L <- BiogTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC1L <- BiogTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC2L <- BiogTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) ### Using the geometric series for the distribution of random assemblages try(bioS0G <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS1G <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS2G <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC0G <- BiogTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC1G <- BiogTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC2G <- BiogTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) ### Using the brocken-stick for the distribution of random assemblages try(bioS0S <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS1S <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS2S <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC0S <- BiogTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC1S <- BiogTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC2S <- BiogTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(simres <- data.frame(nplots, ngrs1, ngrs2, sdlog, p.ecoS0 =ecoS0$pval[2], p.ecoS1 =ecoS1$pval[2], p.ecoS2 =ecoS2$pval[2], p.ecoC0 =ecoC0$pval[2], p.ecoC1 =ecoC1$pval[2], p.ecoC2 =ecoC2$pval[2], p.bioS0L =bioS0L$pval[2], p.bioS1L =bioS1L$pval[2], p.bioS2L =bioS2L$pval[2], p.bioC0L =bioC0L$pval[2], p.bioC1L =bioC1L$pval[2], p.bioC2L =bioC2L$pval[2], p.bioS0G =bioS0G$pval[2], p.bioS1G =bioS1G$pval[2], p.bioS2G =bioS2G$pval[2], p.bioC0G =bioC0G$pval[2], p.bioC1G =bioC1G$pval[2], p.bioC2G =bioC2G$pval[2], p.bioS0S =bioS0S$pval[2], p.bioS1S =bioS1S$pval[2], p.bioS2S =bioS2S$pval[2], p.bioC0S =bioC0S$pval[2], p.bioC1S =bioC1S$pval[2], p.bioC2S =bioC2S$pval[2])) try(simres10 <- rbind(simres10, simres)) } ################################################################################## ################################################################################## ###### SCENARIO 11: DRAWING TWO RANDOM SETS OF SAMPLES FROM TWO METACOMMUNITIES ###### WITH DIFFERENT SPECIES RICHNESS BUT SIMILAR DISTRIBUTION OF RELATIVE ABUNDANCES ################################################################################## ################################################################################## n <- 750 niter <- 200 simres11 <- data.frame(spn1 = NULL, spn2 = NULL, nplots = NULL, ngrs1 = NULL, ngrs2 = NULL, sdlog = NULL, p.ecoS0 = NULL, p.ecoS1 = NULL, p.ecoS2 = NULL, p.ecoC0 = NULL, p.ecoC1 = NULL, p.ecoC2 = NULL, p.bioS0L = NULL, p.bioS1L = NULL, p.bioS2L = NULL, p.bioC0L = NULL, p.bioC1L = NULL, p.bioC2L = NULL, p.bioS0G = NULL, p.bioS1G = NULL, p.bioS2G = NULL, p.bioC0G = NULL, p.bioC1G = NULL, p.bioC2G = NULL, p.bioS0S = NULL, p.bioS1S = NULL, p.bioS2S = NULL, p.bioC0S = NULL, p.bioC1S = NULL, p.bioC2S = NULL) for (i in 1:n) { ### Simulate metacommunities spn1 <- 50 spn2 <- 150 nplots <- round(runif(1, 10, 50)) ngrs1 <- round(runif(1, 5, max(nplots)-5)) ngrs2 <- nplots-ngrs1 sdlog <- runif(1, log(1.1), log(33)) com1 <- round(rlnorm(spn1, 3, sdlog)) com2 <- round(rlnorm(spn2, 3, sdlog)) ### Simulate random sampling Ss <- 50 samples.a <- sample(paste("sp",1:length(com1)), rpois(1, Ss), replace=TRUE, prob=com1) samples.a <- data.frame(table(samples.a)) colnames(samples.a) <- c("species", "sample1") for(t in 2:(ngrs1)) { sample2 <- sample(paste("sp",1:length(com1)), rpois(1, Ss), replace=TRUE, prob=com1) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", paste("sample", t, sep="")) samples.a <- merge(samples.a, sample2, by="species", all=TRUE) } samples.b <- sample(paste("sp",1:length(com2)), rpois(1, Ss), replace=TRUE, prob=com2) samples.b <- data.frame(table(samples.b)) colnames(samples.b) <- c("species", paste("sample", ngrs1+1, sep="")) for(t in 2:(ngrs2)) { sample2 <- sample(paste("sp",1:length(com2)), rpois(1, Ss), replace=TRUE, prob=com2) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", paste("sample", ngrs1+t, sep="")) samples.b <- merge(samples.b, sample2, by="species", all=TRUE) } df <- merge(samples.a, samples.b, by="species", all=TRUE) rownames(df) <- df$species df[is.na(df)] <- 0 df <- df[,-1] ### Apply ecological null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers try(ecoS0 <- EcoTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoS1 <- EcoTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoS2 <- EcoTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoC0 <- EcoTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoC1 <- EcoTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoC2 <- EcoTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) ### Apply biogeographical null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers ### Using the log-normal for the distribution of random assemblages try(bioS0L <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS1L <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS2L <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC0L <- BiogTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC1L <- BiogTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC2L <- BiogTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) ### Using the geometric series for the distribution of random assemblages try(bioS0G <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS1G <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS2G <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC0G <- BiogTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC1G <- BiogTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC2G <- BiogTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) ### Using the brocken-stick for the distribution of random assemblages try(bioS0S <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS1S <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS2S <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC0S <- BiogTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC1S <- BiogTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC2S <- BiogTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(simres <- data.frame(spn1, spn2, nplots, ngrs1, ngrs2, sdlog, p.ecoS0 =ecoS0$pval[2], p.ecoS1 =ecoS1$pval[2], p.ecoS2 =ecoS2$pval[2], p.ecoC0 =ecoC0$pval[2], p.ecoC1 =ecoC1$pval[2], p.ecoC2 =ecoC2$pval[2], p.bioS0L =bioS0L$pval[2], p.bioS1L =bioS1L$pval[2], p.bioS2L =bioS2L$pval[2], p.bioC0L =bioC0L$pval[2], p.bioC1L =bioC1L$pval[2], p.bioC2L =bioC2L$pval[2], p.bioS0G =bioS0G$pval[2], p.bioS1G =bioS1G$pval[2], p.bioS2G =bioS2G$pval[2], p.bioC0G =bioC0G$pval[2], p.bioC1G =bioC1G$pval[2], p.bioC2G =bioC2G$pval[2], p.bioS0S =bioS0S$pval[2], p.bioS1S =bioS1S$pval[2], p.bioS2S =bioS2S$pval[2], p.bioC0S =bioC0S$pval[2], p.bioC1S =bioC1S$pval[2], p.bioC2S =bioC2S$pval[2])) try(simres11 <- rbind(simres11, simres)) } ################################################################################## ################################################################################## ###### SCENARIO 12: DRAWING TWO SETS OF RANDOM SAMPLES FROM TWO METACOMMUNITIES ###### WITH DIFFERENT DISTRIBUTION OF RELATIVE ABUNDANCES BUT SIMILAR RICHNESS ###### AND COMPOSITION ################################################################################## ################################################################################## n <- 750 niter <- 200 simres12 <- data.frame(nplots = NULL, ngrs1 = NULL, ngrs2 = NULL, sdlog1 = NULL, sdlog2 = NULL, p.ecoS0 = NULL, p.ecoS1 = NULL, p.ecoS2 = NULL, p.ecoC0 = NULL, p.ecoC1 = NULL, p.ecoC2 = NULL, p.bioS0L = NULL, p.bioS1L = NULL, p.bioS2L = NULL, p.bioC0L = NULL, p.bioC1L = NULL, p.bioC2L = NULL, p.bioS0G = NULL, p.bioS1G = NULL, p.bioS2G = NULL, p.bioC0G = NULL, p.bioC1G = NULL, p.bioC2G = NULL, p.bioS0S = NULL, p.bioS1S = NULL, p.bioS2S = NULL, p.bioC0S = NULL, p.bioC1S = NULL, p.bioC2S = NULL) for (i in 1:n) { ### Simulate metacommunities spn <- 100 nplots <- round(runif(1, 10, 50)) ngrs1 <- round(runif(1, 5, max(nplots)-5)) ngrs2 <- nplots-ngrs1 sdlog1 <- runif(1, log(1.1), log(33)) sdlog2 <- ifelse(sdlog1 < 1.5, sdlog1 + 1, sdlog1 - 1) com1 <- round(rlnorm(spn, 3, sdlog1)) com2 <- round(rlnorm(spn, 3, sdlog2)) ### Simulate random sampling Ss <- 50 samples.a <- sample(paste("sp",1:length(com1)), rpois(1, Ss), replace=TRUE, prob=com1) samples.a <- data.frame(table(samples.a)) colnames(samples.a) <- c("species", "sample1") for(t in 2:(ngrs1)) { sample2 <- sample(paste("sp",1:length(com1)), rpois(1, Ss), replace=TRUE, prob=com1) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", paste("sample", t, sep="")) samples.a <- merge(samples.a, sample2, by="species", all=TRUE) } samples.b <- sample(paste("sp",1:length(com2)), rpois(1, Ss), replace=TRUE, prob=com2) samples.b <- data.frame(table(samples.b)) colnames(samples.b) <- c("species", paste("sample", ngrs1+1, sep="")) for(t in 2:(ngrs2)) { sample2 <- sample(paste("sp",1:length(com2)), rpois(1, Ss), replace=TRUE, prob=com2) sample2 <- data.frame(table(sample2)) colnames(sample2) <- c("species", paste("sample", ngrs1+t, sep="")) samples.b <- merge(samples.b, sample2, by="species", all=TRUE) } df <- merge(samples.a, samples.b, by="species", all=TRUE) rownames(df) <- df$species df[is.na(df)] <- 0 df <- df[,-1] ### Apply ecological null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers try(ecoS0 <- EcoTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoS1 <- EcoTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoS2 <- EcoTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoC0 <- EcoTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoC1 <- EcoTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(ecoC2 <- EcoTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, by=c(rep("A", ngrs1), rep("B", ngrs2)))) ### Apply biogeographical null model tests to compare individual-based ### and coverage-based rarefaction curves for different Hill numbers ### Using the log-normal for the distribution of random assemblages try(bioS0L <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS1L <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS2L <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC0L <- BiogTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC1L <- BiogTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC2L <- BiogTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, distr="lnorm", by=c(rep("A", ngrs1), rep("B", ngrs2)))) ### Using the geometric series for the distribution of random assemblages try(bioS0G <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS1G <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS2G <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC0G <- BiogTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC1G <- BiogTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC2G <- BiogTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, distr="geom", by=c(rep("A", ngrs1), rep("B", ngrs2)))) ### Using the brocken-stick for the distribution of random assemblages try(bioS0S <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 0, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS1S <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 1, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioS2S <- BiogTest.sample(df, niter=niter, method = "sample-size", q = 2, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC0S <- BiogTest.sample(df, niter=niter, method = "coverage", q = 0, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC1S <- BiogTest.sample(df, niter=niter, method = "coverage", q = 1, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(bioC2S <- BiogTest.sample(df, niter=niter, method = "coverage", q = 2, MARGIN=2, trace=FALSE, distr="stick", by=c(rep("A", ngrs1), rep("B", ngrs2)))) try(simres <- data.frame(nplots, ngrs1, ngrs2, sdlog1, sdlog2, p.ecoS0 =ecoS0$pval[2], p.ecoS1 =ecoS1$pval[2], p.ecoS2 =ecoS2$pval[2], p.ecoC0 =ecoC0$pval[2], p.ecoC1 =ecoC1$pval[2], p.ecoC2 =ecoC2$pval[2], p.bioS0L =bioS0L$pval[2], p.bioS1L =bioS1L$pval[2], p.bioS2L =bioS2L$pval[2], p.bioC0L =bioC0L$pval[2], p.bioC1L =bioC1L$pval[2], p.bioC2L =bioC2L$pval[2], p.bioS0G =bioS0G$pval[2], p.bioS1G =bioS1G$pval[2], p.bioS2G =bioS2G$pval[2], p.bioC0G =bioC0G$pval[2], p.bioC1G =bioC1G$pval[2], p.bioC2G =bioC2G$pval[2], p.bioS0S =bioS0S$pval[2], p.bioS1S =bioS1S$pval[2], p.bioS2S =bioS2S$pval[2], p.bioC0S =bioC0S$pval[2], p.bioC1S =bioC1S$pval[2], p.bioC2S =bioC2S$pval[2])) try(simres12 <- rbind(simres12, simres)) } ################################################################################## ################################################################################## ###### SAVE ALL RESULTS ################################################################################## ################################################################################## save(simres1, simres2, simres3, simres4, simres5, simres6, simres7, simres8, simres9, simres10, simres11, simres12, file="simres.RData")