Apéndice

En este apéndice incluimos una rutina de R que ilustra los tres tipos de escenarios (factores que promueven ensambles neutrales y ensambles convergentes o divergentes basados en el nicho ecológico), así como la forma de evaluarlos (Fig. 4). Los tres primeros ejemplos representan cada uno de los tres escenarios y para eso utilizan ensambles inventados simulando datos de presencia y ausencia de especies. El modelo neutral utilizado en este caso es el propuesto por Chase et al. (2011), que a su vez es una ligera modificación del índice de dissimilitud de Raup Crick. El cuarto ejemplo utiliza ensambles con datos de abundancia para cada especie, y el modelo neutral utilizado es el mismo que propusieron Stegen et al. (2013b), que es básicamente el modelo de Chase et al. (2011), pero incluyendo las abundancias de las especies.

Escenario 1: factor que promueve ensambles divergentes basados en el nicho

Para el escenario 1, imaginemos que un grupo de ecólogos quiere una excusa para conocer las montañas, y entonces deciden analizar si los procesos que estructuran los ensambles de pastizales de llanura (“Llanura”) difieren de los de montaña (“Montaña”). Para eso, toman muestras de vegetación en la llanura, y otras 10 en la montaña, a lo largo de un gradiente altitudinal. Supongamos que en los dos tipos de pastizales se encuentran las mismas especies, pero que en la llanura cualquier especie puede estar en cualquier parcela, mientras que en la montaña hay una clara zonación de especies con la altura. Con el siguiente código generamos los ensambles de los dos tipos de pastizales.

#crea una matrix de 20 especies para la llanura, asignando la presencia de especies al azar
com.llan <- matrix(sample(c(1,0),200, replace=T),ncol=20, nrow=10)
#en este tratamiento hay sólo 4 especies por parcela y una mínima superposición entre parcelas
com.mont <- matrix(rep(c(rep(1,3),rep(0,19)),9),ncol=20, nrow=10, byrow = T)
## Warning in matrix(rep(c(rep(1, 3), rep(0, 19)), 9), ncol = 20, nrow = 10, :
## la longitud de los datos [198] no es un submúltiplo o múltiplo del número
## de filas [10] en la matriz
#crea un vector que define el tipo de pastizal
past<-c(rep("Llanura",10),rep("Montaña",10))
#une en una misma matriz los ensambles de los dos tipos de pastizales
com.past<-rbind(com.llan,com.mont)

Luego de procesar las muestras y pasar los datos, se sientan a analizarlos. * Antes de empezar, cargan los paquetes que van a usar para los análisis y los gráficos

library(vegan)
library(ggplot2)
library(cowplot)
#calcula la dissimilitud de Bray Curtis (que se usará para estimar la diversidad beta)
dist.bray<-vegdist(com.past, method="bray")
#calcula la dissimilitud de con el índice de Raup Crick (compara lo observado con un mondelo neutral)
dist.raup<-vegdist(com.past, method="raup")
#calcula la distancia media a los centroides (diversidad beta) en los dos tipos de pastizales
beta.pre<-betadisper(dist.bray, past)
#calcula la distancia post comparación con el modelo nulo
beta.post<-betadisper(dist.raup, past)
#evalúa las diferencias en diversidad beta
permutest(beta.pre, pairwise = T)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq Mean Sq     F N.Perm Pr(>F)    
## Groups     1 0.37670 0.37670 84.36    999  0.001 ***
## Residuals 18 0.08038 0.00447                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##            Llanura Montaña
## Llanura              0.001
## Montaña 3.2512e-08
#evalúa las diferencias post modelo nulo
permutest(beta.post, pairwise = T)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)    
## Groups     1 0.108756 0.108756 30.874    999  0.001 ***
## Residuals 18 0.063406 0.003523                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##            Llanura Montaña
## Llanura              0.001
## Montaña 2.8303e-05
#primero la figura de diversidad beta
#crea una matriz con la distancia a los centroides (para usar en la figura)
df.pre<-data.frame(beta.pre$distances, past)
#crea una figura con las distancias (la media en rojo) al centroide (diversidad beta) en los dos tipos de pastizales
fig.bray<-ggplot(data=df.pre, aes(past,beta.pre.distances))+
  geom_jitter(width=.2, height=0,size=4, alpha=0.75)+
  stat_summary(fun.data = "mean_cl_boot", color="red", size = 1, alpha=0.75)+
  scale_x_discrete(name="")+
  labs(y="")+
  ggtitle("Bray-Curtis\ndissimilarity")
#luego la figura post modelo neutral
df.post<-data.frame(beta.post$distances, past)
#crea una figura con las distancias de Raup Crick (la media en rojo). A las distancias
#les resta 0,5 y las multiplica por dos para llevarlas a una escala con rango de -1 a 1
fig.raup<-ggplot(data=df.post, aes(past,(beta.post.distances-0.5)*2))+
  geom_hline(yintercept = 0, color="blue", linetype=2, size=1.5, alpha=.75)+
  geom_jitter(width=.2, height=0,size=4, alpha=0.75)+
  stat_summary(fun.data = "mean_cl_boot", color="red", size = 1, alpha=0.75)+
  scale_x_discrete(name="")+
  ggtitle("Raup-Crick\ndissimilarity")+
  labs(y="")
#la tercer figura es el NMDS post modelo neutral
#calcula los valores para un NMDS de dos dimensiones
nmds.post<-metaMDS(dist.raup,trymax = 250, autotransform =FALSE, wascores = FALSE, expand = FALSE)
#crea una matriz de datos para usar en los gráficos
NMDS.post<-data.frame(MDS1= nmds.post$points[,1], MDS2 = nmds.post$points[,2], past=past)
#crea un NMDS con los valores post modelo nulo. Cuanto más grandes son las elipses, 
#mayor es la divergencia entre parcelas
fig.post<-ggplot(data = NMDS.post, aes(MDS1, MDS2, shape=past))+
  stat_ellipse(aes(linetype = past), size=1)+
  scale_shape_manual(values=c(15,0), name = "")+
  geom_point(size=4)+
  guides(linetype=F)
#y por último la figura que integra las tres anteriores
# une las tres figuras anteriores en una sola para poder compararlas
plot_grid(fig.bray+theme(legend.position = "none"),
          fig.raup+theme(legend.position = "none"),
          fig.post+theme(legend.position = "top"),
          ncol=3, rel_widths = c(1,1,2))

En la figura del medio, valores cercanos a 0 (línea azul) reflejan ensambles neutrales. Cuanto más negativos son los valores, mayor es la influencia del nicho ecológico, llevando a ensambles convergentes. Cuanto más positivos son los valores, la influencia del nicho también es mayor, pero los ensambles en este caso difieren más que lo esperable por azar. Es decir, reflejan condiciones heterogéneas entre parcelas de un mismo tratamiento. Como era esperable, los resultados muestran que la diversidad beta fue mayor en la montaña, y que luego de comparar los datos con el modelo neutral, el pastizal de llanura muestra valores cercanos a ensambles puramente neutrales, mientras que el de montaña muestra valores mayores de dispersión, sugiriendo la influencia de factores relacionados al nicho ecológico que llevan a ensambles divergentes.

Escenario 2: ensambles neutrales

Para el escenario 2 imaginemos que los ecólogos, con tal de seguir viajando, deciden evaluar si los tipos de ensambles de insectos difieren entre dos especies de bromelias, una de hojas chicas (“Chica”) y la otra de hojas grandes (“Grande”). Supongamos que, para los insectos que viven en esas bromelias, el principal factor limitante es el tamaño de la hoja, y que no hay especies dominantes. Con el siguiente código generamos los ensambles de las dos especies de bromelias.

#crea una matrix de 20 especies para la especie de hojas chicas, asignando la presencia de especies al azar
com.ch <- matrix(sample(c(1,0),200, replace=T),ncol=20, nrow=10)
#lo mismo para la especie de hojas grandes, pero con mayor riqueza
com.gr <- matrix(sample(c(1,1,1,1,0),200, replace=T),ncol=20, nrow=10)
#crea un vector con el tamaño de hoja correspondiente a cada fila
hoja<-c(rep("Chica",10),rep("Grande",10))
#une en una misma matriz los ensambles de las dos especies
com.bro<-rbind(com.ch,com.gr)

Al regresar del viaje, los ecólogos procesan las muestras y pasan los datos. * Antes de analizar los datos cargan los paquetes necesarios para los análisis y los gráficos.

library(vegan)
library(ggplot2)
library(cowplot)
#calcula la dissimilitud de Bray Curtis (que se usará para estimar la diversidad beta)
dist.bray<-vegdist(com.bro, method="bray")
#calcula la dissimilitud de con el índice de Raup Crick (compara lo observado con un mondelo neutral)
dist.raup<-vegdist(com.bro, method="raup")
#calcula la distancia media a los centroides (diversidad beta) en las dos especies
beta.pre<-betadisper(dist.bray, hoja)
#calcula la distancia post comparación con el modelo nulo
beta.post<-betadisper(dist.raup, hoja)
#evalúa las diferencias en diversidad beta
permutest(beta.pre, pairwise = T)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)    
## Groups     1 0.213058 0.213058 64.602    999  0.001 ***
## Residuals 18 0.059365 0.003298                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##             Chica Grande
## Chica              0.001
## Grande 2.2926e-07
#evalúa las diferencias post modelo nulo
permutest(beta.post, pairwise = T)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.000316 0.0003164 0.0185    999  0.899
## Residuals 18 0.307112 0.0170618                     
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##          Chica Grande
## Chica           0.898
## Grande 0.89319
#primero la figura de diversidad beta
#crea una matriz con la distancia a los centroides (para usar en la figura)
df.pre<-data.frame(beta.pre$distances, hoja)
#crea una figura con las distancias (la media en rojo) al centroide (diversidad beta) para las dos especies
fig.bray<-ggplot(data=df.pre, aes(hoja,beta.pre.distances))+
  geom_jitter(width=.2, height=0,size=4, alpha=0.75)+
  stat_summary(fun.data = "mean_cl_boot", color="red", size = 1, alpha=0.75)+
  scale_x_discrete(name="")+
  labs(y="")+
  ggtitle("Bray-Curtis\ndissimilarity")
#luego la figura post modelo neutral
df.post<-data.frame(beta.post$distances, hoja)
#crea una figura con las distancias de Raup Crick (la media en rojo). A las distancias
#les resta 0,5 y las multiplica por dos para llevarlas a una escala con rango de -1 a 1
fig.raup<-ggplot(data=df.post, aes(hoja,(beta.post.distances-0.5)*2))+
  geom_hline(yintercept = 0, color="blue", linetype=2, size=1.5, alpha=.75)+
  geom_jitter(width=.2, height=0,size=4, alpha=0.75)+
  stat_summary(fun.data = "mean_cl_boot", color="red", size = 1, alpha=0.75)+
  scale_x_discrete(name="")+
  ggtitle("Raup-Crick\ndissimilarity")+
  labs(y="")
#la tercer figura es el NMDS post modelo neutral
#calcula los valores para un NMDS de dos dimensiones
nmds.post<-metaMDS(dist.raup,trymax = 250, autotransform =FALSE, wascores = FALSE, expand = FALSE)
#crea una matriz de datos para usar en los gráficos
NMDS.post<-data.frame(MDS1= nmds.post$points[,1], MDS2 = nmds.post$points[,2], hoja=hoja)
#crea un NMDS con los valores post modelo nulo. Cuanto más grandes son las elipses, 
#mayor es la divergencia entre parcelas
fig.post<-ggplot(data = NMDS.post, aes(MDS1, MDS2, shape=hoja))+
  stat_ellipse(aes(linetype = hoja), size=1)+
  scale_shape_manual(values=c(15,0), name = "")+
  geom_point(size=4)+
  guides(linetype=F)
#y por último la figura que integra las tres anteriores
# une las tres figuras anteriores en una sola para poder compararlas
plot_grid(fig.bray+theme(legend.position = "none"),
          fig.raup+theme(legend.position = "none"),
          fig.post+theme(legend.position = "top"),
          ncol=3, rel_widths = c(1,1,2))

Como era esperable, los resultados mostraron que las bromelias de hojas chicas presentaron una mayor diversidad beta (figura de la izquierda; hubo un menor número de especies de insecto en las plantas de hoja chica, y como la diversidad regional es la misma para las dos especies, eso implica menores probabilidades de solapamiento de especies entre plantas). Sin embargo, los procesos de estructuración de los ensambles en ambas especies fueron los mismos (valores de dispersión similares entre sí; figura de la derecha), y próximos a los esperados para ensambles neutrales (figura central).

Escenario 3: factor que promueve ensambles convergentes basados en el nicho ecológico

Para el escenario 3, imaginemos que los ecólogos se gastaron todos los viáticos, y por lo tanto deciden estudiar los tipos de ensamble que se generan en etapas tempranas de sucesiones secundarias en un intermareal rocoso cercano. Para ello, deciden realizar un experimento en el que limpian de animales y plantas a 10 unidades experimentales de 40 x 40 cm en rocas ubicadas en la zona baja de un intermareal rocoso (“Abajo”) y otras 10 en la zona alta (“Arriba”), más lejana al agua, y por lo tanto, con mayor estrés por desecación. Supongamos además que en el ensamble de la zona baja la colonización es mayormente al azar y, al ser un ambiente favorable, prácticamente todas las especies son capaces de sobrevivir. En la zona alta la colonización tambien es al azar pero el estrés es tan grande que la mitad de las especies del intermareal se ven excluidas porque no sobreviven a esas condiciones ambientales. Luego de 1 mes, registran las especies presentes en cada unidad experimental, obteniendo datos similares a los generados a continuación.

#crea una matrix de 20 especies para la zona baja, asignando la presencia de especies al azar
com.ab <- matrix(sample(c(1,0),200, replace=T),ncol=20, nrow=10)
#lo mismo para la zona alta, pero sólo están presentes las primeras 10 especies
com.ar <- matrix(c(sample(c(1,0),100, replace=T),rep(0,100)),ncol=20, nrow=10)
#crea un vector que define la zona del intermareal
zona<-c(rep("Abajo",10),rep("Arriba",10))
#une en una misma matriz los ensambles de las dos zonas
com.zona<-rbind(com.ab,com.ar)
library(vegan)
library(ggplot2)
library(cowplot)
#calcula la dissimilitud de Bray Curtis (que se usará para estimar la diversidad beta)
dist.bray<-vegdist(com.zona, method="bray")
#calcula la dissimilitud de con el índice de Raup Crick (compara lo observado con un mondelo neutral)
dist.raup<-vegdist(com.zona, method="raup")
#calcula la distancia media a los centroides (diversidad beta) para las dos zonas
beta.pre<-betadisper(dist.bray, zona)
#calcula la distancia post comparación con el modelo nulo
beta.post<-betadisper(dist.raup, zona)
#evalúa las diferencias en diversidad beta
permutest(beta.pre, pairwise = T)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.010818 0.010818 0.7382    999  0.397
## Residuals 18 0.263782 0.014655                     
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##          Abajo Arriba
## Abajo           0.395
## Arriba 0.40153
#evalúa las diferencias post modelo nulo
permutest(beta.post, pairwise = T)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
## Groups     1 0.25244 0.252445 13.815    999  0.003 **
## Residuals 18 0.32893 0.018274                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##            Abajo Arriba
## Abajo             0.004
## Arriba 0.0015784
#primero la figura de diversidad beta
#crea una matriz con la distancia a los centroides (para usar en la figura)
df.pre<-data.frame(beta.pre$distances, zona)
#crea una figura con las distancias (la media en rojo) al centroide (diversidad beta) para las zonas
fig.bray<-ggplot(data=df.pre, aes(zona,beta.pre.distances))+
  geom_jitter(width=.2, height=0,size=4, alpha=0.75)+
  stat_summary(fun.data = "mean_cl_boot", color="red", size = 1, alpha=0.75)+
  scale_x_discrete(name="")+
  labs(y="")+
  ggtitle("Bray-Curtis\ndissimilarity")
#luego la figura post modelo neutral
df.post<-data.frame(beta.post$distances, zona)
#crea una figura con las distancias de Raup Crick (la media en rojo). A las distancias
#les resta 0,5 y las multiplica por dos para llevarlas a una escala con rango de -1 a 1
fig.raup<-ggplot(data=df.post, aes(zona,(beta.post.distances-0.5)*2))+
  geom_hline(yintercept = 0, color="blue", linetype=2, size=1.5, alpha=.75)+
  geom_jitter(width=.2, height=0,size=4, alpha=0.75)+
  stat_summary(fun.data = "mean_cl_boot", color="red", size = 1, alpha=0.75)+
  scale_x_discrete(name="")+
  ggtitle("Raup-Crick\ndissimilarity")+
  labs(y="")
#la tercer figura es el NMDS post modelo neutral
#calcula los valores para un NMDS de dos dimensiones
nmds.post<-metaMDS(dist.raup,trymax = 250, autotransform =FALSE, wascores = FALSE, expand = FALSE)
#crea una matriz de datos para usar en los gráficos
NMDS.post<-data.frame(MDS1= nmds.post$points[,1], MDS2 = nmds.post$points[,2], zona=zona)
#crea un NMDS con los valores post modelo nulo. Cuanto más grandes son las elipses, 
#mayor es la divergencia entre parcelas
fig.post<-ggplot(data = NMDS.post, aes(MDS1, MDS2, shape=zona))+
  stat_ellipse(aes(linetype = zona), size=1)+
  scale_shape_manual(values=c(15,0), name = "")+
  geom_point(size=4)+
  guides(linetype=F)
#y por último la figura que integra las tres anteriores
# une las tres figuras anteriores en una sola para poder compararlas
plot_grid(fig.bray+theme(legend.position = "none"),
          fig.raup+theme(legend.position = "none"),
          fig.post+theme(legend.position = "top"),
          ncol=3, rel_widths = c(1,1,2))

Como era esperable, los resultados muestran que, a pesar de que ambas zonas presentan una diversidad beta similar (figura de la izquierda), al comparar los datos con el modelo neutral observamos que la zona alta parece tener mayor influencia de factores relacionados al nicho ecológico que llevan a ensambles convergentes (figuras del centro y de la derecha). En otras palabras, el estrés por desecación actúa como un filtro que elimina una parte de las especies del intermareal y lleva a ensambles más convergentes.

Ejemplo con datos de abundancias por especie

El último ejemplo sirve para mostrar el escenario 3 (podría haber sido cualquiera de los tres) utilizando datos de abundancias por especie. Dado que no hay un paquete de R que haga un modelo neutral con abundancias, lo primero es cargar la función utilizada por Stegen et al. (2013b) para el cálculo de los modelos nulos.

raup_crick_abundance = function(spXsite, plot_names_in_col1=FALSE, classic_metric=TRUE, split_ties=TRUE, reps=9999, set_all_species_equal=FALSE, as.distance.matrix=TRUE, report_similarity=FALSE){
  ##expects a species by site matrix for spXsite, with row names for plots, or optionally plots named in column 1.  By default calculates a modification of the Raup-Crick metric (standardizing the metric to range from -1 to 1 instead of 0 to 1). Specifying classic_metric=TRUE instead calculates the original Raup-Crick metric that ranges from 0 to 1. The option split_ties (defaults to TRUE) adds half of the number of null observations that are equal to the observed number of shared species to the calculation- this is highly recommended.  The argument report_similarity defaults to FALSE so the function reports a dissimilarity (which is appropriate as a measure of beta diversity).  Setting report_similarity=TRUE returns a measure of similarity, as Raup and Crick originally specified.  If ties are split (as we recommend) the dissimilarity (default) and similarity (set report_similarity=TRUE) calculations can be flipped by multiplying by -1 (for our modification, which ranges from -1 to 1) or by subtracting the metric from 1 (for the classic metric which ranges from 0 to 1). If ties are not split (and there are ties between the observed and expected shared number of species) this conversion will not work. The argument reps specifies the number of randomizations (a minimum of 999 is recommended- default is 9999).  set_all_species_equal weights all species equally in the null model instead of weighting species by frequency of occupancy.  
  ##Note that the choice of how many plots (rows) to include has a real impact on the metric, as species and their occurrence frequencies across the set of plots is used to determine gamma and the frequency with which each species is drawn from the null model 
  ##this section moves plot names in column 1 (if specified as being present) into the row names of the matrix and drops the column of names
  if(plot_names_in_col1){
    row.names(spXsite)<-spXsite[,1]
    spXsite<-spXsite[,-1]
  }
  ## count number of sites and total species richness across all plots (gamma)
  n_sites<-nrow(spXsite)
  gamma<-ncol(spXsite)
  ##build a site by site matrix for the results, with the names of the sites in the row and col names:
  results<-matrix(data=NA, nrow=n_sites, ncol=n_sites, dimnames=list(row.names(spXsite), row.names(spXsite)))
  ##make the spXsite matrix into a new, pres/abs. matrix:
  ceiling(spXsite/max(spXsite))->spXsite.inc
  ##create an occurrence vector- used to give more weight to widely distributed species in the null model:
  occur<-apply(spXsite.inc, MARGIN=2, FUN=sum)
  ##create an abundance vector- used to give more weight to abundant species in the second step of the null model:
  abundance<-apply(spXsite, MARGIN=2, FUN=sum)
  ##make_null:
  ##looping over each pairwise community combination:
  for(null.one in 1:(nrow(spXsite)-1)){
    for(null.two in (null.one+1):nrow(spXsite)){
      null_bray_curtis<-NULL
      for(i in 1:reps){
        ##two empty null communities of size gamma:
        com1<-rep(0,gamma)
        com2<-rep(0,gamma)
        ##add observed number of species to com1, weighting by species occurrence frequencies:
        com1[sample(1:gamma, sum(spXsite.inc[null.one,]), replace=FALSE, prob=occur)]<-1
        com1.samp.sp = sample(which(com1>0),(sum(spXsite[null.one,])-sum(com1)),replace=TRUE,prob=abundance[which(com1>0)]);
        com1.samp.sp = cbind(com1.samp.sp,1); # head(com1.samp.sp);
        com1.sp.counts = as.data.frame(tapply(com1.samp.sp[,2],com1.samp.sp[,1],FUN=sum)); colnames(com1.sp.counts) = 'counts'; # head(com1.sp.counts);
        com1.sp.counts$sp = as.numeric(rownames(com1.sp.counts)); # head(com1.sp.counts);
        com1[com1.sp.counts$sp] = com1[com1.sp.counts$sp] + com1.sp.counts$counts; # com1;
        #sum(com1) - sum(spXsite[null.one,]); ## this should be zero if everything work properly
        rm('com1.samp.sp','com1.sp.counts');            
        ##same for com2:
        com2[sample(1:gamma, sum(spXsite.inc[null.two,]), replace=FALSE, prob=occur)]<-1
        com2.samp.sp = sample(which(com2>0),(sum(spXsite[null.two,])-sum(com2)),replace=TRUE,prob=abundance[which(com2>0)]);
        com2.samp.sp = cbind(com2.samp.sp,1); # head(com2.samp.sp);
        com2.sp.counts = as.data.frame(tapply(com2.samp.sp[,2],com2.samp.sp[,1],FUN=sum)); colnames(com2.sp.counts) = 'counts'; # head(com2.sp.counts);
        com2.sp.counts$sp = as.numeric(rownames(com2.sp.counts)); # head(com2.sp.counts);
        com2[com2.sp.counts$sp] = com2[com2.sp.counts$sp] + com2.sp.counts$counts; # com2;
        # sum(com2) - sum(spXsite[null.two,]); ## this should be zero if everything work properly
        rm('com2.samp.sp','com2.sp.counts');
        
        null.spXsite = rbind(com1,com2); # null.spXsite;
        ##calculate null bray curtis
        null_bray_curtis[i] = vegdist(null.spXsite,method='bray');
      }; # end reps loop
      ## empirically observed bray curtis
      obs.bray = vegdist(spXsite[c(null.one,null.two),],method='bray');
      ##how many null observations is the observed value tied with?
      num_exact_matching_in_null = sum(null_bray_curtis==obs.bray);
      ##how many null values are smaller than the observed *dissimilarity*?
      num_less_than_in_null = sum(null_bray_curtis<obs.bray);
      rc = (num_less_than_in_null )/reps; # rc;
      if(split_ties){
        rc = ((num_less_than_in_null +(num_exact_matching_in_null)/2)/reps)
      };
      if(!classic_metric){
        ##our modification of raup crick standardizes the metric to range from -1 to 1 instead of 0 to 1
        rc = (rc-.5)*2
      };
      results[null.two,null.one] = round(rc,digits=2); ##store the metric in the results matrix
      print(c(null.one,null.two,date()));
    }; ## end null.two loop
  }; ## end null.one loop
  if(as.distance.matrix){ ## return as distance matrix if so desired
    results<-as.dist(results)
  } 
  return(results)
}; ## end function

Una vez hecho esto, se analizan los datos al igual que en los ejemplos anteriores. Se cargan los paquetes y los datos necesarios,

library(vegan)
library(ggplot2)
library(cowplot)
#cargamos datos ya disponibles en R 
data("mite")
data("mite.env")

se crean las matrices de dissimilitud,

se calcula la dispersión de los ensambles,

#calcula la distancia media a los centroides (diversidad beta) en tres densidades de arbustos
beta.pre<-betadisper(dist.bray, mite.env$Shrub)
#calcula la distancia post comparación con el modelo nulo
beta.post<-betadisper(dist.raup, mite.env$Shrub)

se las evalúa estadísticamente,

#evalúa las diferencias en diversidad beta
permutest(beta.pre, pairwise = T)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     2 0.03153 0.015767 1.0123    999  0.355
## Residuals 67 1.04355 0.015575                     
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##          None      Few  Many
## None          0.620000 0.509
## Few  0.619711          0.056
## Many 0.511939 0.053032
#evalúa las diferencias post modelo nulo
permutest(beta.post, pairwise = T)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df Sum Sq Mean Sq      F N.Perm Pr(>F)    
## Groups     2 1.4484 0.72418 28.686    999  0.001 ***
## Residuals 67 1.6914 0.02525                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##            None        Few  Many
## None            1.0000e-03 0.001
## Few  1.7584e-08            0.933
## Many 5.3278e-09 9.2926e-01

y se hacen los gráficos,

#primero la figura de diversidad beta
#crea una matriz con la distancia a los centroides (para usar en las figuras)
df.pre<-data.frame(beta.pre$distances, mite.env$Shrub)
#crea una figura con las distancias (la media en rojo) al centroide (diversidad beta) para las tres densidades de arbustos
fig.bray<-ggplot(data=df.pre, aes(mite.env$Shrub,beta.pre.distances))+
  geom_jitter(width=.2, height=0,size=4, alpha=0.75)+
  stat_summary(fun.data = "mean_cl_boot", color="red", size = 1, alpha=0.75)+
  scale_x_discrete(name="")+
  labs(y="")+
  ggtitle("Bray-Curtis\ndissimilarity")
#luego la figura post modelo neutral
#crea una matriz con la distancia a los centroides (para usar en las figuras)
df.post<-data.frame(beta.post$distances, mite.env$Shrub)
#crea una figura con las distancias de Raup Crick (la media en rojo). A las distancias
#les resta 0,5 y las multiplica por dos para llevarlas a una escala con rango de -1 a 1
fig.raup<-ggplot(data=df.post, aes(mite.env$Shrub,(beta.post.distances-0.5)*2))+
  geom_hline(yintercept = 0, color="blue", linetype=2, size=1.5, alpha=.75)+
  geom_jitter(width=.2, height=0,size=4, alpha=0.75)+
  stat_summary(fun.data = "mean_cl_boot", color="red", size = 1, alpha=0.75)+
  scale_x_discrete(name="")+
  ggtitle("Raup-Crick\ndissimilarity")+
  labs(y="")
#la tercer figura es el NMDS post modelo neutral
#calcula los valores para un NMDS de dos dimensiones
nmds.post<-metaMDS(dist.raup,trymax = 250, autotransform =FALSE, wascores = FALSE, expand = FALSE)
#crea una matriz de datos para usar en los gráficos
NMDS.post<-data.frame(MDS1= nmds.post$points[,1], MDS2 = nmds.post$points[,2], trt=mite.env$Shrub)
#crea un NMDS con los valores post modelo nulo. Cuanto más grandes son las elipses, 
#mayor es la divergencia entre parcelas
fig.post<-ggplot(data = NMDS.post, aes(MDS1, MDS2, shape=trt))+
  stat_ellipse(aes(linetype = trt), size=1)+
  scale_shape_manual(values=c(15,0,8), name = "")+
  geom_point(size=4)+
  guides(linetype=F)
#y por último la figura que integra las tres anteriores
# une las tres figuras anteriores en una sola para poder compararlas
plot_grid(fig.bray+theme(legend.position = "none"),
          fig.raup+theme(legend.position = "none"),
          fig.post+theme(legend.position = "top"),
          ncol=3, rel_widths = c(1,1,2))

En el ejemplo, la falta de arbustos lleva una menor riqueza de oribátidos, sin diferencias en diversidad beta entre densidades de arbustos (figura de la izquierda). Sin embargo, al tener en cuenta la mayor dispersión esperable por azar (por el hecho de tener menor riqueza), se hace evidente que los ensambles de oribátidos cuando no hay arbustos son marcadamente convergentes (figuras del centro y de la derecha).