-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathR code.R
127 lines (119 loc) · 5.81 KB
/
R code.R
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
###################################
#R code to reproduce models tested in Table 2
#Bastazini, V.A.G., Ferreira, P.M., Azambuja, B.O., Casas, G., Debastiani, V.J., Guimarães, P.R. & Pillar, V.D. 2017. Untangling the Tangled Bank: A Novel Method for Partitioning the Effects of Phylogenies and Traits on Ecological Networks.
#Evolutionary Biology 44(3): 312–324. DOI: 10.1007/s11692-017-9409-8
#Last updated: 2020-03-19
#Contacts:
#bastazini.vinicius@gmail.com
#vanderleidebastiani@yahoo.com.br
####################################
###Packages
require(vegan)
require(SYNCSA)
####################################
### Function mantel.residuals
# it extracts residuals from a mantel matrix correlation (see Fig. 2)
####################################
mantel.residuals <- function(A, B){
a <- as.vector(as.dist(A))
b <- as.vector(as.dist(B))
mod <- lm(a~b)
RESIDUALS <- residuals(mod)
N <- dim(as.matrix(A))[1]
RES <- matrix(NA, N, N)
OR <- 0
for(j in 1:N){
for(i in 1:N){
if(i==j){
RES[i,j] <- 0
}
if(i>j){
OR <- OR+1
RES[i,j] <- RESIDUALS[OR]
}
}
}
R <- sqrt(summary(mod)$r.squared)
dis.residuals <- as.dist(RES, diag = TRUE)
results <- list(R = R, Residuals = dis.residuals)
return(results)
}
####################################
### Function MATRIX.P1
#this is a modification of the function matrix.P from SYNCSA
#which allows the computation of matrices with some rows that sum to 0
####################################
matrix.p1 <- function (comm, dist.spp) {
matrix.w <- sweep(comm, 1, rowSums(comm), "/")
for(i in 1:dim(comm)[1]){
for(j in 1:dim(comm)[2]){
if(is.nan(matrix.w[i,j])){
matrix.w[i,j] <- 0
}
}
}
similar.phy <- 1 - (dist.spp/max(dist.spp))
matrix.phy <- 1/colSums(similar.phy)
matrix.q <- sweep(similar.phy, 1, matrix.phy, "*")
matrix.P <- matrix.w %*% matrix.q
return(list(matrix.w = matrix.w, matrix.q = matrix.q, matrix.P = matrix.P))
}
####################################
######### Data entry #########
####################################
plants_traits <- as.matrix(read.csv("./Data/plants_traits.csv", header = TRUE, row.names = 1))
birds_traits <- as.matrix(read.csv("./Data/birds_traits.csv", header = TRUE, row.names = 1))
interactions <- as.matrix(read.csv("./Data/interactions.csv", header = TRUE, row.names = 1))
birds_phylogeny <- as.matrix(read.csv("./Data/birds_phylogeny.csv", header = TRUE, row.names = 1))
plants_phylogeny <- as.matrix(read.csv("./Data/plants_phylogeny.csv", header = TRUE, row.names = 1))
temporal_occurrence <- as.matrix(read.csv("./Data/temporal_occurrence.csv", header = TRUE, row.names = 1))
birds_traits #Bird trait data with traits as rows and species as colums
plants_traits #Plant trait data with traits as rows and species as colums
interactions #Interaction matrix, with plants as rows and birds as colums
birds_phylogeny #A species by species phylogenetic distance matrix
plants_phylogeny #A species by species phylogenetic distance matrix
temporal_occurrence # A co-occurrence matrix, with plants as rows and birds as colums
############################################################################
######### REMOVING PHYLOGENETIC EFFECT FROM TRAITS #########
############################################################################
sa1 <- vegdist(birds_traits, method = "gower", na.rm = TRUE); sa1
pa1 <- vegdist(birds_phylogeny, method = "gower", na.rm = TRUE); pa1
sp1 <- vegdist(plants_traits, method = "gower", na.rm = TRUE); sp1
pp1 <- vegdist(plants_phylogeny, method = "gower", na.rm = TRUE); pp1
mant1 <- mantel.residuals(sa1, pa1); mant1
mant2 <- mantel.residuals(sp1, pp1); mant2
sa <- as.matrix(mant1$Residuals); sa # distance matriz between species described by their traits, removing the effect of the phylogeny
sb <- as.matrix(mant2$Residuals); sb
##################################################################################################################################
#### PROBABILISTIC INTERACTION MATRICES WEIGHTED BY TRAITS AFTER REMOVING PHYLOGENETIC SIGNAL and accounting for temporal variation
##################################################################################################################################
ua <- matrix.p1(t(interactions), sa)
ya <- ua$matrix.P; ya
up <- matrix.p1(interactions, sb)
yp <- up$matrix.P; yp
ya <- t(ya)*temporal_occurrence
yp <- yp*temporal_occurrence
ya <- t(ya)
yp <- t(yp)
##################################################################################################################################
#### PROBABILISTIC INTERACTION MATRICES WEIGHTED BY PHYLOGENETIC RESEMBLANCE and accounting for temporal variation
##################################################################################################################################
dpa <- matrix.p1(t(interactions), birds_phylogeny)
pa <- dpa$matrix.P
dpp <- matrix.p1(interactions, plants_phylogeny)
pp <- dpp$matrix.P
pat <- t(pa)*temporal_occurrence
pat <- t(pat)
ppt <- (pp)*temporal_occurrence
ppt <- t(ppt)
##############################
#### MODELS TESTED IN TABLE 2
##############################
#Functional component, removing phylogenetic signal
ya.yp <- cor.matrix(mx1 = ua$matrix.w, mx2 = ua$matrix.q, x = ya, y = yp, method="pearson", dist = "euclidean", permutations = 9999, norm = FALSE); ya.yp
#Phylogenetic signal in bird species
ya.pa <- cor.matrix(mx1 = ua$matrix.w, mx2 = ua$matrix.q, x = ya, y = pa, method="pearson", dist = "euclidean", permutations = 9999, norm = FALSE); ya.pa
#Phylogenetic signal in plant species
yp.pp <- cor.matrix(mx1 = up$matrix.w, mx2 = up$matrix.q, x = t(yp), y = pp, method = "pearson", dist = "euclidean", permutations = 9999, norm = FALSE); yp.pp
#Phylogenetic component
pa.pp <- cor.matrix(mx1 = ua$matrix.w, mx2 = ua$matrix.q, x = pat, y = ppt, method = "pearson", dist = "euclidean", permutations = 9999, norm = FALSE); pa.pp