raster-project/fit-species-richness-model.R at master · ricardog/raster-project · GitHub

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
139
140
141
142
143
144
145
146
147
148
149
150
151
#!/usr/bin/env Rscript

##
## "Parse" command-line arguments
##
library(argparser, quietly=TRUE)
p <- arg_parser("Run PREDICTS species richness model")
p <- add_argument(p, "--site-data", help="diversity rds file")
p <- add_argument(p, "--out", help="filename prefix for output files",
                  default='sr-model')
argv <- parse_args(p)

if (is.na(argv$site_data) || is.na(argv$out)) {
    print(p)
    stop()
}

##
## Load packages
##
cat("****\n* Loading packages\n")
library(roquefort, quietly = TRUE)

##
## Read stored site data
##
sites.div <- readRDS(argv$site_data)

##
## Run the model
##
cat("****\n* Running species richness model\n")
model.data <- sites.div[,c('Species_richness', 'UI', 'logHPD.rs', 'logDistRd.rs',
                           'LandUse', 'UseIntensity', 'SS', 'SSB', 'SSBS')]
model.data <- na.omit(model.data)
saveRDS(model.data, file = paste(argv$out, '-data.rds', sep=''))
model.data <- model.data[1:floor(nrow(model.data) / 4), ]
sModel <- glmer(Species_richness ~ UI +
                    poly(logHPD.rs, 2) +
                    poly(logDistRd.rs, 2) +
                    LandUse:poly(logHPD.rs, 2) +
                    LandUse:poly(logDistRd.rs, 2) +
                    (1 + LandUse | SS) +
                    (1 | SSB) +
                    (1 | SSBS),
                family="poisson",data=model.data,
                control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
stopifnot(sModel$converged)
save(sModel, file = paste(argv$out, '.rds', sep=''))

##
## Plot species richness graphs (uses functions in roquefort package)
##
cat("****\n* Plotting species richness results\n")
png("species-richness-response-plot.png",
    width=8.5, height=15, units="cm", res=1200)
par(mfrow=c(3,1))

PlotErrBar(model=sModel, data=sModel@frame, responseVar="Richness",
           logLink="e", catEffects="UI", forPaper=TRUE)
mtext("A", side=3, line=-1, adj=-0.15)

PlotContEffects(model=sModel, data=model.data, effects="logHPD.rs",
                otherContEffects="logDistRd.rs",
                otherFactors=list(UI="Primary Vegetation Minimal use"),
                xlab="Human population density",
                ylab="Species richness",
                byFactor="LandUse",
                seMultiplier=0.5,
                logLink="e")
mtext("B", side=3, line=-0.1, adj=-0.19)

PlotContEffects(model=sModel, data=model.data, effects="logDistRd.rs",
                otherContEffects="logHPD.rs",
                otherFactors=list(LandUse="Primary Vegetation",
                                  UI="Primary Vegetation Minimal use"),
                xlab="Distance to nearest road",
                ylab="Species richness",
                byFactor="LandUse",
                seMultiplier=0.5,
                logLink="e")
mtext("C", side=3, line=-0.1, adj=-0.19)
dev.off()

##
## Extract human population density coefficients
##
hpd.coefs <- numeric()
TempMod <- lm(sModel@frame$`poly(logHPD.rs, 2)`[,1]~model.data$logHPD.rs)
hpd.coefs['Poly1_I'] <- TempMod$coefficients['(Intercept)']
hpd.coefs['Poly1_1'] <- TempMod$coefficients['model.data$logHPD.rs']
TempMod <- lm(sModel@frame$`poly(logHPD.rs, 2)`[,2]~model.data$logHPD.rs+I(model.data$logHPD.rs^2))
hpd.coefs['Poly2_I'] <- TempMod$coefficients['(Intercept)']
hpd.coefs['Poly2_1'] <- TempMod$coefficients['model.data$logHPD.rs']
hpd.coefs['Poly2_2'] <- TempMod$coefficients['I(model.data$logHPD.rs^2)']
write.csv(hpd.coefs, file = paste(argv$out, '-hpd-coefs.csv', sep=''))

##
## Extract road-distance coefficients
##
rddist.coefs <- numeric()
TempMod <- lm(sModel@frame$`poly(logDistRd.rs, 2)`[,1]~model.data$logDistRd.rs)
rddist.coefs['Poly1_I'] <- TempMod$coefficients['(Intercept)']
rddist.coefs['Poly1_1'] <- TempMod$coefficients['model.data$logDistRd.rs']
TempMod <- lm(sModel@frame$`poly(logDistRd.rs, 2)`[,2]~model.data$logDistRd.rs+I(model.data$logDistRd.rs^2))
rddist.coefs['Poly2_I'] <- TempMod$coefficients['(Intercept)']
rddist.coefs['Poly2_1'] <- TempMod$coefficients['model.data$logDistRd.rs']
rddist.coefs['Poly2_2'] <- TempMod$coefficients['I(model.data$logDistRd.rs^2)']
write.csv(rddist.coefs, file = paste(argv$out, '-rddist-coefs.csv', sep=''))

##
## Extract land-use coefficients
##
coefs <- numeric()
coefs['PriMin'] <- 0.0
coefs['PriLt'] <- (fixef(sModel))['UIPrimary Vegetation Light use']
coefs['PriInt'] <- (fixef(sModel))['UIPrimary Vegetation Intense use']
coefs['SecMin'] <- (fixef(sModel))['UISecondary Vegetation Minimal use']
coefs['SecLt'] <- (fixef(sModel))['UISecondary Vegetation Light use']
coefs['CrpMin'] <- (fixef(sModel))['UICropland Minimal use']
coefs['CrpLt'] <- (fixef(sModel))['UICropland Light use']
coefs['CrpInt'] <- (fixef(sModel))['UICropland Intense use']
coefs['PasMin'] <- (fixef(sModel))['UIPasture Minimal use']
coefs['PasLt'] <- (fixef(sModel))['UIPasture Light use']
coefs['PasInt'] <- (fixef(sModel))['UIPasture Intense use']
coefs['UrbMin'] <- (fixef(sModel))['UIUrban Minimal use']
coefs['UrbInt'] <- (fixef(sModel))['UIUrban Intense use']
coefs['hpd_1'] <- (fixef(sModel))['poly(logHPD.rs, 2)1']
coefs['hpd_2'] <- (fixef(sModel))['poly(logHPD.rs, 2)2']
coefs['rddist_1'] <- (fixef(sModel))['poly(logDistRd.rs, 2)1']
coefs['rddist_2'] <- (fixef(sModel))['poly(logDistRd.rs, 2)2']

coefs['hpd_1.Sec'] <- (fixef(sModel))['poly(logHPD.rs, 2)1:LandUseSecondary Vegetation'] + coefs['hpd_1']
coefs['hpd_1.Crp'] <- (fixef(sModel))['poly(logHPD.rs, 2)1:LandUseCropland'] + coefs['hpd_1']
coefs['hpd_1.Pas'] <- (fixef(sModel))['poly(logHPD.rs, 2)1:LandUsePasture'] + coefs['hpd_1']
coefs['hpd_1.Urb'] <- (fixef(sModel))['poly(logHPD.rs, 2)1:LandUseUrban'] + coefs['hpd_1']
coefs['hpd_2.Sec'] <- (fixef(sModel))['poly(logHPD.rs, 2)2:LandUseSecondary Vegetation'] + coefs['hpd_2']
coefs['hpd_2.Crp'] <- (fixef(sModel))['poly(logHPD.rs, 2)2:LandUseCropland'] + coefs['hpd_2']
coefs['hpd_2.Pas'] <- (fixef(sModel))['poly(logHPD.rs, 2)2:LandUsePasture'] + coefs['hpd_2']
coefs['hpd_2.Urb'] <- (fixef(sModel))['poly(logHPD.rs, 2)2:LandUseUrban'] + coefs['hpd_2']

coefs['rddist_1.Sec'] <- (fixef(sModel))['poly(logDistRd.rs, 2)1:LandUseSecondary Vegetation'] + coefs['rddist_1']
coefs['rddist_1.Crp'] <- (fixef(sModel))['poly(logDistRd.rs, 2)1:LandUseCropland'] + coefs['rddist_1']
coefs['rddist_1.Pas'] <- (fixef(sModel))['poly(logDistRd.rs, 2)1:LandUsePasture'] + coefs['rddist_1']
coefs['rddist_1.Urb'] <- (fixef(sModel))['poly(logDistRd.rs, 2)1:LandUseUrban'] + coefs['rddist_1']
coefs['rddist_2.Sec'] <- (fixef(sModel))['poly(logDistRd.rs, 2)2:LandUseSecondary Vegetation'] + coefs['rddist_2']
coefs['rddist_2.Crp'] <- (fixef(sModel))['poly(logDistRd.rs, 2)2:LandUseCropland'] + coefs['rddist_2']
coefs['rddist_2.Pas'] <- (fixef(sModel))['poly(logDistRd.rs, 2)2:LandUsePasture'] + coefs['rddist_2']
coefs['rddist_2.Urb'] <- (fixef(sModel))['poly(logDistRd.rs, 2)2:LandUseUrban'] + coefs['rddist_2']

write.csv(coefs, file = paste(argv$out, '-coefs.csv', sep=''))