-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathScript02
67 lines (49 loc) · 3.15 KB
/
Script02
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
########################################################################################################################################################################
# Script 02 - Meta-analysis
# Haubrock et al. The invasion curve and dynamics of a European-wide introduced species
#
# This script includes the code for the following steps:
# (1) synthesize the results using meta-analytical models,
#
########################################################################################################################################################################
(1) synthesize the results using meta-analytical models
# test run for the overall model
res.rma.Abund <- rma(yi = S_statistic, vi = old.variance, method="REML", data= master_file, control=list(maxiter=1000))
res.rma.Abund # inspect model
#Consideration of explanatory variables
master_file$Latitude<-as.numeric(master_file$Latitude) # ensure data is numeric
master_file$Longitude<-as.numeric(master_file$Longitude) # ensure data is numeric
master_file$StudyLength.s <- decostand(master_file$Study_length, "standardize") # standardize accordingly
# the standardization is to be repeated with every environmental variable included.
# in this works case, this includes: the average temperature over the period, average precipiation over the period, the change in temperature (trends S-statistic), change in precipitation (trends S-statistic),
# years after monitoring started, elevation, elevation drop, slope, ppt, runoff, NH4, NO3, percentage of urban areas in the site's proximity, percentage of open water bodies in the site's proximity
res.rma.Abund <- rma(yi = S_statistic, vi = old.variance, mods=~Lat.s+Lon.s+S_Temp.s+S_Prec.s+Mean_temp.s+Mean_prec.s+Years_after_monitoring_start.s+Elevation.s+Elevation_drop.s+Slope.s+Runoff.s+Urban.s+Water.s-1, method="REML", data= total, control=list(maxiter=1000)) # model under consideration of all covariates
res.rma.Abund
# precautious test for variance inflation
vif(res.rma.Abund)
# this step will identify the best fit covariates for the model.
res <- glmulti(S_statistic ~ Lat.s+Lon.s+S_Temp.s+S_Prec.s+Mean_temp.s+Mean_prec.s+Years_after_monitoring_start.s+Elevation.s+Elevation_drop.s+Slope.s+Runoff.s+Urban.s+Water.s,
data=total,
level=1, method="ML",crit="aicc", confsetsize=128)
print(res)
plot(res)
top <- weightable(res)
top <- top[top$aicc <= min(top$aicc) + 2,]
top
#incluing the best fitting covariates into the meta regression model
res.rma.Abund <- rma(yi = S_statistic, vi = old.variance, mods=~Lat.s + Runoff.s + Water.s-1, method="REML", data= total, control=list(maxiter=1000))
res.rma.Abund # inspect the model result
#export or raw plot for inspection
svg("all_trends.svg")
plot_model(res.rma.Abund, vline.color = "black", transform = NULL)
dev.off()
#export a better version for later display
# forest plot with extra annotations
svg("trends.svg",height=30,width=15)
op <- par(cex=.8, font=2)
forest(res.rma.Abund, xlim=c(-400,400),
ilab=cbind(total$site_id,
round(total$P_value, digits = 3)),
ilab.xpos=c(-250,-170), cex=.8, header="Study & Site",
mlab="")
op <- par(cex=0.6, font=2)