BD in amphibians, Norte de Santander. A GLM aproach
Acevedo, et.al.,
24 de Octubre de 2018
- hacer mapa por especie sacar coeficientes del modelo 4 (mod_lmer4)
BD
Read packages
knitr::opts_chunk$set(echo = TRUE)
library(sp)
library(rgdal)## rgdal: version: 1.3-4, (SVN revision 766)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/diego.lizcano/Documents/R/win-library/3.4/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/diego.lizcano/Documents/R/win-library/3.4/rgdal/proj
## Linking to sp version: 1.3-1library(raster)
library(rgdal)
library(spatial.tools)
library(dismo)
library(rasterVis)## Loading required package: lattice## Loading required package: latticeExtra## Loading required package: RColorBrewerlibrary(xtable)
library(sf)## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3library(mapview)
library(corrgram)##
## Attaching package: 'corrgram'## The following object is masked from 'package:latticeExtra':
##
## panel.ellipse## The following object is masked from 'package:lattice':
##
## panel.filllibrary (corrplot) # correlations## corrplot 0.84 loadedlibrary(lme4)## Loading required package: Matrixlibrary(nlme)##
## Attaching package: 'nlme'## The following object is masked from 'package:lme4':
##
## lmList## The following object is masked from 'package:raster':
##
## getDatalibrary(ggplot2)##
## Attaching package: 'ggplot2'## The following object is masked from 'package:latticeExtra':
##
## layerlibrary(sjPlot)
library(mgcv) # GAM## This is mgcv 1.8-17. For overview type 'help("mgcv-package")'.library(MuMIn) # model average
library(mgcv) # visual GAM
library(ggeffects) # to predict modelsRead Data
# Path to working directory
path <- "C:/Users/diego.lizcano/Documents/GitHub/Bd_N_Santander/data"
shp <- "C:/Users/diego.lizcano/Documents/GitHub/Bd_N_Santander/Data2"
bio <- "C:/Users/diego.lizcano/Documents/GitHub/Bd_N_Santander/Capas2"
BD4<-readOGR(dsn=shp, layer="Bd_generos1")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\diego.lizcano\Documents\GitHub\Bd_N_Santander\Data2", layer: "Bd_generos1"
## with 150 features
## It has 5 fields
## Integer64 fields read as strings: Diagnostic# define geo coords and put to SpatialPointDataFrames
crs.geo <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") # geographical, datum WGS84
proj4string(BD4) <- crs.geo # define projection system of our data
BD <- st_as_sf(BD4) # Convert foreign object to an sf object
pos <- grep(" ", BD$Species)
keep <- substr(pos, 1, pos-1)
########################
# BD$genus <-
###cargando variables ambientales
files <- list.files(path=bio, pattern='tif', full.names=TRUE )View files
# View(files)
# Agrupar las variables ambientales en un objeto stack
predictors <- stack(files)
### clac slope aspect
slope_aspect <- terrain(predictors$DEM30, opt = c("slope", "aspect"), unit = "degrees")
predictors <- addLayer(predictors, slope_aspect)
deptos <-raster::getData ('GADM', country='COL', level=1)
# names(predictors)
# plot(predictors$Bio1)
# plot(BD1, pch=18, add=TRUE, col="blue", main="Batrachochytrium dendrobatidis")
mapview(BD, zcol = "Diagnostic", legend = TRUE, map.types = c("Esri.WorldShadedRelief", "Esri.WorldImagery"), color = "grey40")levelplot(predictors, main="all predictors")
# plot(BD["Diagnostic"], add=T)
# levelplot(predictors)
### correlation
jnk=layerStats(predictors, 'pearson', na.rm=T)
corr_matrix=jnk$'pearson correlation coefficient'
# corrgram(corr_matrix, order=NULL, lower.panel=panel.shade,
# upper.panel=NULL, text.panel=panel.txt,
# main="Car")
##### correlation plot
# corrplot(corr_matrix, order = "FPC",type = "lower", tl.cex = 0.8) # tl.pos = "lower",
####### function to add significance
# mat : is a matrix of data
# ... : further arguments to pass to the native R cor.test function
cor.mtest <- function(mat, ...) {
mat <- as.matrix(mat)
n <- ncol(mat)
p.mat<- matrix(NA, n, n)
diag(p.mat) <- 0
for (i in 1:(n - 1)) {
for (j in (i + 1):n) {
tmp <- cor.test(mat[, i], mat[, j], ...)
p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
}
}
colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
p.mat
}
# matrix of the p-value of the correlation
# use the function
p.mat <- cor.mtest(predictors)
# head(p.mat[, 1:5])
##### correlation plot
# Leave blank on no significant coefficient
corrplot(corr_matrix, p.mat = p.mat, sig.level = 0.05, order="FPC", type = "lower", tl.cex = 0.8) # tl.pos = "lower",
Remove higly correlated layers
## remove correlated layers
predictors2 <- dropLayer(predictors, c(1,2,3,4,9,15,16,18,19))
### correlation
jnk=layerStats(predictors2, 'pearson', na.rm=T)
corr_matrix=jnk$'pearson correlation coefficient'
# matrix of the p-value of the correlation
p.mat <- cor.mtest(predictors2)
corrplot(corr_matrix, p.mat = p.mat, sig.level = 0.05, order="FPC", type = "lower", tl.cex = 0.8) # tl.pos = "lower",
print("predictores no correlacionados")## [1] "predictores no correlacionados"names(predictors2)## [1] "Bio13" "Bio14" "Bio15" "Bio16" "Bio18" "Bio19" "Bio2"
## [8] "Bio3" "Bio4" "Bio7" "DEM30" "slope" "aspect"Building GLM models
Now only significant explanatory covars (layers) from Final Model glm1
##########################################
## now only covars from Final Model glm1
##########################################
glm2 <- glm(pa ~ Bio13 + Bio16 + Bio19 + Bio2 + Bio4 + DEM30,
family = binomial(link = "logit"), data=training)
glm3 <- glm(pa ~ Bio13 + Bio19 + Bio2 + DEM30,
family = binomial(link = "logit"), data=training)
glm4 <- glm(pa ~ Bio13 + Bio2 + DEM30,
family = binomial(link = "logit"), data=training)
glm5 <- glm(pa ~ Bio13 + Bio19 + DEM30, # only precipitation
family = binomial(link = "logit"), data=training)
glm6 <- glm(pa ~ Bio13 + Bio19 + genus, # only precipitation
family = binomial(link = "logit"), data=training)
AICc(glm2, glm3, glm4, glm5)# AICc should be used instead AIC when sample size is small in comparison to the number of estimated parameters (Burnham & Anderson 2002 recommend its use when n / K < 40). 105/6 = 17.5
summary(glm2)##
## Call:
## glm(formula = pa ~ Bio13 + Bio16 + Bio19 + Bio2 + Bio4 + DEM30,
## family = binomial(link = "logit"), data = training)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3859 -0.9372 0.4960 0.8673 1.9627
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.8295 0.3453 2.402 0.016284 *
## Bio13 -2.1483 1.0090 -2.129 0.033245 *
## Bio16 2.0441 1.0584 1.931 0.053437 .
## Bio19 0.5698 0.2745 2.076 0.037939 *
## Bio2 -1.6674 0.4871 -3.423 0.000619 ***
## Bio4 1.0658 0.5854 1.821 0.068646 .
## DEM30 -1.6701 0.5330 -3.133 0.001730 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 144.41 on 104 degrees of freedom
## Residual deviance: 115.59 on 98 degrees of freedom
## AIC: 129.59
##
## Number of Fisher Scoring iterations: 4#### multiple group effects with multiple group effect terms factor species.
# #### Random-effects
# mod_lmer1 <- glmer(pa ~ Bio13 + Bio16 + Bio19 + Bio2 + Bio4 + DEM30 + (1|genus), data=training, #fulldat,
# family = binomial(link = "logit"))
# #### fixed-effects
# mod_lmer2 <- glmer(pa ~ I(Bio13^2) : I(Bio16^2) + I(Bio19^2) + I(Bio2^2) + (1|genus), data=training, #fulldat,
# family = binomial(link = "logit"))
#
# mod_lmer3 <- glmer(pa ~ Bio16 + DEM30 + (Bio2|genus), data=training, #fulldat,
# family = binomial(link = "logit"))
#
# mod_lmer4 <- glmer(pa ~ DEM30 + Bio2 + (Bio16|genus), data=training, #fulldat,
# family = binomial(link = "logit"))
#
# mod_lmer5 <- glmer(pa ~ Bio16 : Bio2 * (1|genus), data=training, #fulldat,
# family = binomial(link = "logit"))
mod_lmer1 <- glmer(pa ~ Bio13 + Bio2 + DEM30 + (1|genus), data=training, #fulldat,
family = binomial(link = "logit"))
mod_lmer2 <- glmer(pa ~ Bio14 + Bio4 + DEM30 + (1|genus), data=training, #fulldat,
family = binomial(link = "logit"))
mod_lmer3 <- glmer(pa ~ Bio19 + Bio2 + Bio4 + (1|genus), data=training, #fulldat,
family = binomial(link = "logit"))
mod_lmer4 <- glmer(pa ~ Bio2 + Bio14 + DEM30 + (1|genus), data=training, #fulldat,
family = binomial(link = "logit"))
mod_lmer5 <- glmer(pa ~ Bio4 + Bio13 + DEM30 + (1|genus), data=training, #fulldat,
family = binomial(link = "logit"))
mod_lmer6 <- glmer(pa ~ Bio2 + DEM30 + (1|genus), data=training, #fulldat,
family = binomial(link = "logit"))
summary(mod_lmer1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: pa ~ Bio13 + Bio2 + DEM30 + (1 | genus)
## Data: training
##
## AIC BIC logLik deviance df.resid
## 120.1 133.4 -55.1 110.1 100
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1374 -0.4869 0.1277 0.3650 2.6837
##
## Random effects:
## Groups Name Variance Std.Dev.
## genus (Intercept) 10.46 3.235
## Number of obs: 105, groups: genus, 21
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.64059 1.17641 1.395 0.16314
## Bio13 0.04847 0.30354 0.160 0.87314
## Bio2 -1.68854 0.63172 -2.673 0.00752 **
## DEM30 -1.37954 0.67739 -2.037 0.04170 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Bio13 Bio2
## Bio13 -0.078
## Bio2 -0.176 0.313
## DEM30 -0.092 0.411 0.822summary(mod_lmer2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: pa ~ Bio14 + Bio4 + DEM30 + (1 | genus)
## Data: training
##
## AIC BIC logLik deviance df.resid
## 128.1 141.4 -59.0 118.1 100
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1500 -0.6440 0.1535 0.3124 1.8985
##
## Random effects:
## Groups Name Variance Std.Dev.
## genus (Intercept) 9.054 3.009
## Number of obs: 105, groups: genus, 21
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.9911 1.1041 1.803 0.0713 .
## Bio14 0.6339 0.3790 1.673 0.0944 .
## Bio4 0.7121 0.6407 1.111 0.2664
## DEM30 0.3713 0.3928 0.945 0.3445
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Bio14 Bio4
## Bio14 0.228
## Bio4 0.277 0.099
## DEM30 0.229 0.556 0.345summary(mod_lmer3)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: pa ~ Bio19 + Bio2 + Bio4 + (1 | genus)
## Data: training
##
## AIC BIC logLik deviance df.resid
## 123.2 136.5 -56.6 113.2 100
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0433 -0.5947 0.1349 0.3591 1.7744
##
## Random effects:
## Groups Name Variance Std.Dev.
## genus (Intercept) 7.819 2.796
## Number of obs: 105, groups: genus, 21
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.0220 1.0708 1.888 0.0590 .
## Bio19 0.5174 0.3534 1.464 0.1431
## Bio2 -1.0017 0.4243 -2.361 0.0182 *
## Bio4 1.4044 0.7472 1.879 0.0602 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Bio19 Bio2
## Bio19 0.060
## Bio2 -0.266 -0.363
## Bio4 0.331 0.496 -0.370summary(mod_lmer4) ## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: pa ~ Bio2 + Bio14 + DEM30 + (1 | genus)
## Data: training
##
## AIC BIC logLik deviance df.resid
## 118.4 131.6 -54.2 108.4 100
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.98110 -0.41960 0.09507 0.30828 2.82645
##
## Random effects:
## Groups Name Variance Std.Dev.
## genus (Intercept) 12.83 3.582
## Number of obs: 105, groups: genus, 21
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.8503 1.2986 1.425 0.15420
## Bio2 -1.6903 0.6150 -2.749 0.00599 **
## Bio14 0.5055 0.3936 1.284 0.19911
## DEM30 -1.1184 0.6691 -1.672 0.09461 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Bio2 Bio14
## Bio2 -0.138
## Bio14 0.155 -0.013
## DEM30 0.006 0.747 0.327AICc(mod_lmer1, mod_lmer2, mod_lmer3, mod_lmer4, mod_lmer5)#, mod_lmer6 )G A M Models
#### G A M Models
# ####
# gam1 <- gam (pa ~ s(Bio13) + s(Bio16) + s(Bio19) +
# s(Bio2) + s(Bio4) + DEM30,
# family = binomial(link = "logit"), data=training)
#
# gam2 <- gam (pa ~ te(Bio13, Bio16) +
# s(Bio19) + te(Bio2, Bio4) ,
# family = binomial(link = "logit"), method="REML",
# select=TRUE, data=training)
#
# gam3 <- gam (pa ~ Bio13 + Bio16 +
# Bio19 + Bio2 + Bio4 + s(DEM30),
# family = binomial, data=training)
#
# gam4 <- gam (pa ~ s(DEM30) + te(Bio16, Bio2),
# family = binomial, data=training)
#
#
#
# AICc(gam1, gam2, gam3, gam4)
#
# #vis.gam(gam1, view=c("Bio6", "Bio13"), type="response")
# gam.ME <- gamm(pa ~ s(Bio13) + s(Bio16) + s(Bio19) +
# s(Bio2) + s(Bio4) + s(DEM30) , method="REML",
# family = binomial(link = "logit"), data=training )#, random=list(genus=~1))
#
################ plot using sjPlot
# plot_model(mod_lmer1, vline.color = "gray",
# sort.est = TRUE)
plot_model(mod_lmer4, type = "re", sort.est = TRUE) #plots conditionals of random effects## Sorting each group of random effects ('sort.all') is not possible when 'facets = TRUE'.
plot_model(mod_lmer4, type = "est", #fixed effects
show.ci = TRUE, sort.est = TRUE)
dat <-ggpredict(mod_lmer4, c( "DEM30", "Bio2","genus") )
ggplot(dat, aes(x = x, y = predicted, colour = group)) +
stat_smooth( se = FALSE, fullrange = TRUE) +
facet_wrap(~facet)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
# dat2 <- ggpredict(mod_lmer6, terms = c("Bio2", "DEM30", "genus"))
# ggplot(dat2, aes(x = x, y = predicted, colour = group)) +
# stat_smooth( fullrange = TRUE) +
# facet_wrap(~facet)
# Probabilities of fixed effects depending on grouping level (random intercept)
# plot_model(mod_lmer1, type = "ri.pc",
# show.se = TRUE)
# sjp.lmer(mod_lmer1, type = "fe",
# show.ci = TRUE, sort.est = TRUE)
# sjp.lmer(mod_lmer1, type = "slope")
## evaluacion
ge2 <- evaluate(testing[testing$pa==1,], testing[testing$pa==0,], glm2)
ge2## class : ModelEvaluation
## n presences : 30
## n absences : 15
## AUC : 0.6911111
## cor : 0.2996497
## max TPR+TNR at : 0.7853403ge_lmer <- evaluate(testing[testing$pa==1,], testing[testing$pa==0,], mod_lmer4, allow.new.levels = TRUE) ####### change by genus
ge_lmer## class : ModelEvaluation
## n presences : 30
## n absences : 15
## AUC : 0.7266667
## cor : 0.3379722
## max TPR+TNR at : 0.3646983ge3 <- evaluate(testing[testing$pa==1,], testing[testing$pa==0,], glm3, allow.new.levels = TRUE) ####### change by genus
ge3## class : ModelEvaluation
## n presences : 30
## n absences : 15
## AUC : 0.7133333
## cor : 0.3529656
## max TPR+TNR at : -0.1701741par(mfrow=c(1,3))
plot(ge2, 'ROC', main="GLM")
plot(ge_lmer, 'ROC', main="GLMM")
plot(ge3, 'ROC', main="GLMM")
# dev.off()
##mapeo del modelo GLMM 4
var_orig <- scale(predictors2, center=TRUE, scale=TRUE)
pg4 <- raster::predict(scale(predictors2, center=TRUE, scale=TRUE), mod_lmer4, re.form=~0, type = "response")
plot(pg4)
list_ras <- list()
for(i in 1:length(unique(BD$Genero))){
list_ras[[i]] <- raster::predict(scale(predictors2, center=TRUE, scale=TRUE), mod_lmer4, const=(data.frame(genus=unique(BD$Genero)[i])),
type = "response", allow.new.levels=T)
}
genus_maps <- stack(list_ras)
names(genus_maps) <- unique(BD$Genero)
# library(ggiraphExtra)
for(i in 1:length(unique(BD$Genero))){
par(mfrow=c(1,2))
plot(genus_maps[[i]],
main=paste('Pred. prob. BD for', names(genus_maps[[i]])))
plot(BD4, add=TRUE, col='grey50')
plot(deptos, add=TRUE, border='white')
tr <- threshold(ge_lmer, 'spec_sens')
plot(genus_maps[[i]]>tr,
main=paste('BD pres/abs for', names(genus_maps[[i]])))
plot(BD4, add=TRUE, col='grey50')
plot(deptos, add=TRUE, border='white')
# points(testing, pch='+')
# dev.off()
}
##mapeo del modelo lmer1
###########################
###### no es posible... falta genero
############################
# pg_lmer <- predict(scale(predictors2, center=TRUE, scale=TRUE), mod_lmer1, ext=predictors2)
# par(mfrow=c(1,2))
# plot(pg_lmer, main='GLM/binomial, BD')
# plot(BD4, add=TRUE, border='dark grey')
# tr <- threshold(ge_lmer, 'spec_sens')
# plot(pg>tr, main='presence/backgence')
# plot(BD4, add=TRUE, border='dark grey')
# points(testing, pch='+')
# dev.off()Comentarios
BIO1 = Annual Mean Temperature BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp)) BIO3 = Isothermality (BIO2/BIO7) (* 100) BIO4 = Temperature Seasonality (standard deviation *100) BIO5 = Max Temperature of Warmest Month BIO6 = Min Temperature of Coldest Month BIO7 = Temperature Annual Range (BIO5-BIO6) BIO8 = Mean Temperature of Wettest Quarter BIO9 = Mean Temperature of Driest Quarter BIO10 = Mean Temperature of Warmest Quarter BIO11 = Mean Temperature of Coldest Quarter BIO12 = Annual Precipitation BIO13 = Precipitation of Wettest Month BIO14 = Precipitation of Driest Month BIO15 = Precipitation Seasonality (Coefficient of Variation) BIO16 = Precipitation of Wettest Quarter BIO17 = Precipitation of Driest Quarter BIO18 = Precipitation of Warmest Quarter BIO19 = Precipitation of Coldest Quarter