/
EW_nullmodel.Rmd
135 lines (101 loc) · 5.17 KB
/
EW_nullmodel.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
---
title: "E-W structure using null models"
output:
html_document:
df_print: paged
---
First, read the csv file with the list of number of individuals per OTU in each mountain. | Primero, se lee el documento csv con la lista de individuos de cada OTU por montaña
```{r}
OTUs <- read.csv("OTUs.csv")
row.names(OTUs)<-OTUs$OTU
OTUs<-OTUs[1:44,2:8] #delete column OTU, not needed.
OTUs
```
Convert the csv file to matrix and tranpose the matrix so that otus are in columns and sites in rows. Then, calculate a presence-absence matrix, where each cell will be first converted to T if there are one or more individuals in a site, and F if there are 0 individuals. Then, T are transformed to 1, and F to 0. | Se traspone la matriz para que las especies sean las columnas (variables) y los sitios los renglones. Después, se calcula una matriz de presencia-ausencia, en donde cada celda es transformada en T si hay al menos un individuo o más en el sitio, y F si hay 0 individuos. Después se transforman las T en 1 y las F en 0.
```{r}
OTUs<-as.matrix(OTUs) #La convertimos en matriz
OTUs<-t(OTUs) #Trasponemos la matriz para que las especies sean las columnas (variables) y los sitios los renglones.
image(t(log(OTUs+1)), axes=FALSE, ylab="sites", xlab="OTU") #Se ve en la imagen los valores de abundancia, los colores rojos indican abundancias bajas, amarillo a blancos indica abundancias altas.
OTUs.PA<-OTUs>0 #T and F instead of number of individuals
mode(OTUs.PA)<-"integer" # T=1 y F=0
head(OTUs)
OTUs.PA
```
In a null model, empirical data are compared to simulated data using the function oecosimu. The presence-absence matrix (explained in diversity.Rmd) is used to calculate beta diversity. I made two groups representing eastern and western mountains. Since there are 7 mountains, 34 combinations of eastern-western mountains different than the real one can exist. After making two groups of mountains in 34 possible combinations, data are randomized between the two groups. This way, 34 randomized groups are compared with the actual E-W structure we are testing. | Se comparan los datos empíricos con datos simulados con la función oecosimu. Se usa la base de datos de presencia-ausencia construida para calcular la diversidad beta. Hice dos grupos que representan montañas del Este y montañas del Oeste. Debido a que hay 7 montañas, existen 34 posibles combinaciones de agrupamientos E-O. Después de hacer dos grupos de montañas con las 34 posibles combinaciones, los datos fueron aleatorizados entre los dos grupos. De esta manera, los 34 grupos aleatorios son comárados con el agrupamiento que representa la estructura E-O que se está poniendo a prueba.
```{r}
library(vegan)
#Group1: E=An,Bl,To W=Aj,Iz,Ma,Pe
E1<-OTUs.PA[1:3,] # vector E
E1<-colSums(E1) #Add all species from these 3 mountains
W1<-OTUs.PA[4:7,] #vector W
W1<-colSums(W1)
EW1<-matrix(data=(c(E1,W1)),nrow=2,ncol=44,byrow=T,dimnames=list(c("E1","W1"))) #Build matrix E-W
EW1
sim1st<-oecosimu(EW1,nestedbetasor,method="r0",nsimul=999) #randomize sites
sim1st
c(sim1st$statistic[1:3], sim1st$oecosimu$pval[1:3]) # get diversity statistics and their pvalue
#plot
dplot1st<-densityplot(permustats(sim1st))
dplot1st
```
Repeat with all 35 combinations | Se repite el mismo procedimiento con los 35 agrupamientos de montañas.
```{r, warning=FALSE}
## Read file with the 35 groupings
EW_groups<-read.delim("EW_groups.txt", stringsAsFactors = FALSE)
## Perform analysis in a loop
div_results<-as.numeric() # to save div results
for(i in c(1:35)){
# get character vector of group W and E
W<-strsplit(EW_groups[i,2], ",")[[1]]
print(paste("results for grouping", i , "W:"))
print(W)
E<-strsplit(EW_groups[i,3], ",")[[1]]
print("and E:")
print(E)
# Get matrix for group i
W<-OTUs.PA[W,]
E<-OTUs.PA[E,]
W<-colSums(W)
E<-colSums(E)
EW<-matrix(data=(c(E,W)),nrow=2,ncol=44,byrow=T,dimnames=list(c("E","W")))
# Run analysis
sim1st<-oecosimu(EW,nestedbetasor,method="r0",nsimul=999) #randomize sites
print(sim1st)
x<-c(sim1st$statistic[1:3], sim1st$oecosimu$pval[1:3]) # get diversity statistics and their pvalue
div_results<-rbind(div_results,x) # save in table outside loop
# #plot
dplot1st<-densityplot(permustats(sim1st))
print(dplot1st)
}
colnames(div_results)[4:6]<-c("p_turnoever", "p_nestedness", "p_sorensen")
rownames(div_results)<-NULL
div_results<-round(div_results[,],4)
EW_results_null_model<-cbind(EW_groups, div_results)
EW_results_null_model
```
Perform an analysis of variance with permultations | Análisis de varianza con permutaciones
```{r}
permu_results<-numeric()
for(i in 1:35){
# get character vector of group W and E
W<-strsplit(EW_groups[i,2], ",")[[1]]
print(paste("results for grouping", i , "W:"))
print(W)
E<-strsplit(EW_groups[i,3], ",")[[1]]
print("and E:")
print(E)
## Generate EW vector for each grouping
vec<-c(W, E)
vec<-sub("An|Bl|To", "W", vec)
vec<-sub("Aj|Iz|Ma|Pe", "E", vec)
## run permutation analyses
set.seed(100) #set seed to get same results every run
OTUsdist<-betadiver(OTUs.PA,method=1)
AV<-adonis(OTUsdist~as.factor(vec),permutation=999)
AV
permu_results<-rbind(permu_results, AV$aov.tab$Pr[1])
}
EW_results_permu<-cbind(EW_groups, permu_results)
colnames(EW_results_permu)[4]<-"p_value"
EW_results_permu
```