-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy path12_SEM_subset_CP_TSFI_83_isl.R
139 lines (100 loc) · 4.32 KB
/
12_SEM_subset_CP_TSFI_83_isl.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
127
128
129
130
131
132
133
134
135
136
137
138
# Compute SEM Models for 407 world oceanic islands
# with at least 4 alien bird species
rm(list=ls())
library(lavaan)
library(piecewiseSEM)
library(lme4)
library(semEff)
library(tidyverse)
################ Load data ###############
# Clean data were obtained with script 00_create_clean_database
# Load final dataset with all scaled variables used in SEMs
all_ok <- read.csv2("Data/Marino_et_al_DATA_407_ISL.csv")
################ Include CP and TSFI variables ################
# For a subset of islands (n=96), we could obtain data regarding
# the colonization pressure (CP), that is the total number of species
# introduced to an area (Blackburn et al., 2020)
# We thus perform SEM on those 96 islands with and without the CP variable
# to check how this factor can influence our results.
# Filter the islands with CP & TFSI data
isl_CP <- all_ok %>% filter(!is.na(CP)) %>% filter(!is.na(TSFI))
table(isl_CP$Realm)
table(all_ok$Realm)
###### First models: without CP nor TSFI in the subset of 96 islands ######
# Already performed in the script "11_SEM_subset_CP_96_isl.R"
###### Second models: with CP and TSFI in the SEM ######
pop <- lm(pop ~ Area + Lat + SLMP + Dist + TSFI, data = isl_CP)
GDP <- lm(GDP ~ Area + pop + Dist + Elev + Lat, data = isl_CP)
connect <- lm(connect ~ Area + Lat + SLMP + Dist + pop, data = isl_CP)
static <- lm(static_modif ~ GDP + pop + modif_change + Area + Lat + connect, data = isl_CP)
modif <- lm(modif_change ~ GDP + pop + Elev + SLMP + connect, data = isl_CP)
CP <- lm(CP ~ pop + GDP + connect + Dist + Area + static_modif + TSFI, data = isl_CP)
TSFI <- lm(TSFI ~ Dist + Area, data = isl_CP)
# biotic context depends on biogeographic variables
# native SR drives native FD & PD
nat_SR <- lm(SR_nat ~ Area + Dist + Lat + Elev + SLMP, data = isl_CP)
nat_PD <- lm(PD_nat ~ Area + Dist + Lat + Elev + SLMP + SR_nat, data = isl_CP)
nat_FD <- lm(FD_nat ~ Area + Dist + Lat + Elev + SLMP + SR_nat, data = isl_CP)
# alien diversities depend on all contextual variables
alien_SR <- lm(SR_alien ~ Area + Dist + Lat + Elev + SLMP + SR_nat + TSFI+
pop + GDP + modif_change + static_modif + connect + CP, data = isl_CP)
alien_PD <- lm(PD_alien ~ Area + Dist + Lat + Elev + SLMP + SR_nat + PD_nat + TSFI+
pop + GDP + modif_change + static_modif + connect + CP + SR_alien,
data = isl_CP)
alien_FD <- lm(FD_alien ~ Area + Dist + Lat + Elev + SLMP + SR_nat + FD_nat + TSFI+
pop + GDP + modif_change + static_modif + connect + CP + SR_alien,
data = isl_CP)
# build total model
# one for FD and one for pd
cp_pd_mod <- psem(pop, GDP, connect, static, modif, CP, TSFI,
nat_SR, nat_PD, alien_SR, alien_PD,
SR_nat %~~% connect,
SR_nat %~~% TSFI,
SR_nat %~~% CP,
PD_nat %~~% pop,
PD_nat %~~% connect,
PD_nat %~~% GDP,
data = isl_CP)
#basisSet(cp_pd_mod)
d1 <- dSep(cp_pd_mod)
fisherC(cp_pd_mod) # FisherC = 74.307, df = 60, p = 0.101
d1[which.max(abs(d1$Crit.Value)),]
cp_FD_mod <- psem(pop, GDP, connect, static, modif, CP, TSFI,
nat_SR, nat_FD, alien_SR, alien_FD,
SR_nat %~~% connect,
SR_nat %~~% TSFI,
SR_nat %~~% CP,
FD_nat %~~% pop,
FD_nat %~~% static_modif,
FD_nat %~~% connect,
data = isl_CP)
d1 <- dSep(cp_FD_mod)
fisherC(cp_FD_mod) # FisherC = 68.7, df = 60, p = 0.206
d1[which.max(abs(d1$Crit.Value)),]
rsquared(cp_FD_mod)
rsquared(cp_pd_mod)
########## Get direct/ indirect effects of the 2 models with CP
R_sample = 10000
Sys.time()
cp_boot_mod_fd <- bootEff(cp_FD_mod, R = R_sample)
cp_eff_fd <- semEff(cp_boot_mod_fd)
Sys.time()
cp_boot_mod_pd <- bootEff(cp_pd_mod, R = R_sample)
cp_eff_pd <- semEff(cp_boot_mod_pd)
Sys.time()
##### Save outputs
# fd
saveRDS(cp_FD_mod, "Output/12_SEM_CP_TSFI_FD_model.rds")
saveRDS(cp_eff_fd, "Output/12_SEM_CP_TSFI_FD_boot_eff.rds")
# pd
saveRDS(cp_pd_mod, "Output/12_SEM_CP_TSFI_PD_model.rds")
saveRDS(cp_eff_pd, "Output/12_SEM_CP_TSFI_PD_boot_eff.rds")
############## Simple model outputs
# for PD
rsquared(cp_pd_mod)
coefs <- coefs(cp_pd_mod)
coefs[coefs$P.Value<0.05,]
# for FD
rsquared(cp_FD_mod)
coefs <- coefs(cp_FD_mod)
coefs[coefs$P.Value<0.05,]