BD N de Sant. GLM
Acevedo, et.al.,
23 de enero de 2018
BD
Read packages
knitr::opts_chunk$set(echo = TRUE)
library(sp)## Warning: package 'sp' was built under R version 3.4.3library(rgdal)## Warning: package 'rgdal' was built under R version 3.4.3## rgdal: version: 1.2-16, (SVN revision 701)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
## 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.2-5library(raster)## Warning: package 'raster' was built under R version 3.4.3library(rgdal)
library(dismo)## Warning: package 'dismo' was built under R version 3.4.1library(rasterVis)## Warning: package 'rasterVis' was built under R version 3.4.3## Loading required package: lattice## Loading required package: latticeExtra## Loading required package: RColorBrewerlibrary(xtable)Read Data
# Path to working directory
path <- "C:/Users/diego.lizcano/Documents/GitHub/Bd_N_Santander/data"
BD1<-readOGR(dsn=path, layer="BD")## OGR data source with driver: ESRI Shapefile
## Source: "C:/Users/diego.lizcano/Documents/GitHub/Bd_N_Santander/data", layer: "BD"
## with 54 features
## It has 4 fields###cargando variables ambientales
files <- list.files(path=path, pattern='tif', full.names=TRUE )
# View(files)
# Agrupar las variables ambientales en un objeto stack
predictors <- stack(files)
# names(predictors)
plot(predictors$Bio1)
plot(BD1, pch=18, add=TRUE, col="blue", main="Batrachochytrium dendrobatidis")
levelplot(predictors)
#Extraer la informacion de variables ambientales por ptos. con presencia de BD
# head(BD)
BD <- BD1[,-1] # quito la primera columna
head(BD)## X Y Diagnostic
## 1 -73.06614 8.380760 1
## 2 -72.62444 7.936156 1
## 3 -72.63474 7.920514 1
## 4 -72.63305 7.542250 1
## 5 -72.54150 7.815944 1
## 6 -72.57844 7.727333 1names(BD)[1]<-paste("x")
names(BD)[2]<-paste("y") # re nombro encabezados de columnas
# head(BD)
presvals <- (BD)
###presvals <- extract(predictors, BD)
set.seed(2000)
backg <- randomPoints(predictors, n=500, extf = 1.25)
plot(predictors$DEM30)
points(backg, pch=18, col="red")
plot(BD1, pch=18, add=TRUE, col="blue", main="Batrachochytrium dendrobatidis")
#generando datos de prueba y entrenamiento con los datos de presencia y pseudo-ausencia con particion 80-20%
indexes_pres = sample(1:nrow(presvals), size=0.2*nrow(presvals))
indexes_backg = sample(1:nrow(backg), size=0.2*nrow(backg))
###Presencias
#identifica el 20% de los datos como datos prueba/test
pres_test = presvals[indexes_pres,]
dim(pres_test)## [1] 10 3#identifica el 80% de los datos como datos entrenamiento/train
pres_train = presvals[-indexes_pres,]
dim(pres_train)## [1] 44 3# View(pres_train)
###pseudo-ausencia
#identifica el 20% de los datos como datos prueba/test
backg_test = backg[indexes_backg,]
dim(backg_test)## [1] 100 2#identifica el 80% de los datos como datos entrenamiento/train
backg_train = backg[-indexes_backg,]
dim(backg_train)## [1] 400 2# mapeando los datos de entrenamiento y prueba de presencias y backg
r = raster(predictors, 1)
plot(!is.na(r), col=c('white', 'light grey'), legend=FALSE, main="datos de entrenamiento y prueba ")
points(backg_train, pch='-', cex=0.5, col='red')
points(backg_test, pch='-', cex=0.5, col='black')
points(pres_train, pch= '+', col='green')
points(pres_test, pch='+', col='blue')
trainpres <- data.frame( extract(predictors, pres_train) )
trainbackg <- data.frame( extract(predictors, backg_train) )
train <- rbind(trainpres, trainbackg)
pb_train <- c(rep(1, nrow(trainpres)), rep(0, nrow(trainbackg)))
envtrain <- data.frame( cbind(pa=pb_train, train) )
#### extract from stack
testpres <- data.frame( extract(predictors, pres_test) )
testbackg <- data.frame( extract(predictors, backg_test) )
# head(envtrain)
#Ajuste del modelo GLM y predicc??n
# GLM binomial
gm1 <- glm(pa ~ Bio1 + Bio2 + Bio3 + Bio4 + Bio5 + Bio6 + Bio7 + Bio8 + Bio9 + Bio10 + Bio11 + Bio12 + Bio13 + Bio14 + Bio15 + Bio16 + Bio17 + Bio18 + Bio19 + DEM30,family = binomial(link = "logit"), data=envtrain)
summary(gm1)##
## Call:
## glm(formula = pa ~ Bio1 + Bio2 + Bio3 + Bio4 + Bio5 + Bio6 +
## Bio7 + Bio8 + Bio9 + Bio10 + Bio11 + Bio12 + Bio13 + Bio14 +
## Bio15 + Bio16 + Bio17 + Bio18 + Bio19 + DEM30, family = binomial(link = "logit"),
## data = envtrain)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8262 -0.3164 -0.1262 -0.0170 3.3506
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.064e+02 6.167e+01 1.725 0.08455 .
## Bio1 1.782e+01 6.596e+00 2.702 0.00689 **
## Bio2 8.517e+00 7.276e+00 1.171 0.24177
## Bio3 -9.305e-01 7.025e-01 -1.325 0.18531
## Bio4 -1.636e-01 1.670e-01 -0.980 0.32727
## Bio5 4.118e+05 2.125e+06 0.194 0.84635
## Bio6 -4.118e+05 2.125e+06 -0.194 0.84635
## Bio7 -4.118e+05 2.125e+06 -0.194 0.84635
## Bio8 4.979e+00 1.891e+00 2.633 0.00847 **
## Bio9 -7.806e+00 1.730e+00 -4.512 6.41e-06 ***
## Bio10 -1.073e+01 6.323e+00 -1.697 0.08966 .
## Bio11 -2.781e+00 6.782e+00 -0.410 0.68176
## Bio12 -1.019e-03 2.510e-03 -0.406 0.68481
## Bio13 -7.774e-03 1.105e-02 -0.704 0.48157
## Bio14 -4.569e-02 2.321e-02 -1.968 0.04903 *
## Bio15 -1.542e-02 7.796e-02 -0.198 0.84322
## Bio16 1.251e-03 6.157e-03 0.203 0.83901
## Bio17 8.391e-03 1.021e-02 0.822 0.41110
## Bio18 -2.301e-03 2.867e-03 -0.802 0.42228
## Bio19 5.943e-03 2.050e-03 2.899 0.00374 **
## DEM30 -4.150e-03 1.390e-03 -2.985 0.00283 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 286.91 on 443 degrees of freedom
## Residual deviance: 173.29 on 423 degrees of freedom
## AIC: 215.29
##
## Number of Fisher Scoring iterations: 8coef(gm1)## (Intercept) Bio1 Bio2 Bio3 Bio4
## 1.063801e+02 1.782204e+01 8.516710e+00 -9.304713e-01 -1.636418e-01
## Bio5 Bio6 Bio7 Bio8 Bio9
## 4.117560e+05 -4.117581e+05 -4.117646e+05 4.978916e+00 -7.806320e+00
## Bio10 Bio11 Bio12 Bio13 Bio14
## -1.073107e+01 -2.781057e+00 -1.018781e-03 -7.773728e-03 -4.568756e-02
## Bio15 Bio16 Bio17 Bio18 Bio19
## -1.541791e-02 1.250792e-03 8.391070e-03 -2.300673e-03 5.943058e-03
## DEM30
## -4.150483e-03# GLM gausiano
gm2 <- glm(pa ~ Bio1 + Bio2 + Bio3 + Bio4 + Bio5 + Bio6 + Bio7 + Bio8 + Bio9 + Bio10 + Bio11 + Bio12 + Bio13 + Bio14 + Bio15 + Bio16 + Bio17 + Bio18 + Bio19 + DEM30,family = gaussian(link = "identity"), data=envtrain)
summary(gm2)##
## Call:
## glm(formula = pa ~ Bio1 + Bio2 + Bio3 + Bio4 + Bio5 + Bio6 +
## Bio7 + Bio8 + Bio9 + Bio10 + Bio11 + Bio12 + Bio13 + Bio14 +
## Bio15 + Bio16 + Bio17 + Bio18 + Bio19 + DEM30, family = gaussian(link = "identity"),
## data = envtrain)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.41261 -0.14746 -0.06648 0.02656 0.94159
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.241e-01 2.908e+00 0.318 0.750775
## Bio1 4.390e-01 2.971e-01 1.478 0.140282
## Bio2 -1.365e-01 3.345e-01 -0.408 0.683459
## Bio3 1.677e-02 3.356e-02 0.500 0.617505
## Bio4 3.833e-03 6.328e-03 0.606 0.545017
## Bio5 8.142e+04 1.555e+05 0.523 0.600910
## Bio6 -8.142e+04 1.555e+05 -0.523 0.600910
## Bio7 -8.142e+04 1.555e+05 -0.523 0.600910
## Bio8 1.920e-01 8.794e-02 2.183 0.029606 *
## Bio9 -2.058e-01 7.771e-02 -2.649 0.008378 **
## Bio10 -5.171e-01 2.810e-01 -1.840 0.066464 .
## Bio11 1.942e-02 2.746e-01 0.071 0.943663
## Bio12 6.289e-05 1.513e-04 0.416 0.677868
## Bio13 -1.569e-04 5.906e-04 -0.266 0.790655
## Bio14 -2.502e-03 1.494e-03 -1.675 0.094654 .
## Bio15 -1.606e-03 4.799e-03 -0.335 0.738074
## Bio16 -1.167e-04 3.445e-04 -0.339 0.735040
## Bio17 4.443e-05 5.975e-04 0.074 0.940757
## Bio18 -1.478e-04 1.683e-04 -0.878 0.380572
## Bio19 3.181e-04 1.001e-04 3.178 0.001592 **
## DEM30 -3.334e-04 9.517e-05 -3.503 0.000508 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.0774165)
##
## Null deviance: 39.640 on 443 degrees of freedom
## Residual deviance: 32.747 on 423 degrees of freedom
## AIC: 146.51
##
## Number of Fisher Scoring iterations: 2coef(gm2)## (Intercept) Bio1 Bio2 Bio3 Bio4
## 9.241218e-01 4.390238e-01 -1.365051e-01 1.677058e-02 3.833013e-03
## Bio5 Bio6 Bio7 Bio8 Bio9
## 8.141698e+04 -8.141697e+04 -8.141690e+04 1.919517e-01 -2.058373e-01
## Bio10 Bio11 Bio12 Bio13 Bio14
## -5.170827e-01 1.941800e-02 6.289456e-05 -1.568926e-04 -2.501821e-03
## Bio15 Bio16 Bio17 Bio18 Bio19
## -1.605737e-03 -1.166529e-04 4.443371e-05 -1.477622e-04 3.180856e-04
## DEM30
## -3.334307e-04## evaluacion
ge1 <- evaluate(testpres, testbackg, gm1)
ge1## class : ModelEvaluation
## n presences : 10
## n absences : 100
## AUC : 0.867
## cor : 0.2735291
## max TPR+TNR at : -1.809294ge2 <- evaluate(testpres, testbackg, gm2)
ge2## class : ModelEvaluation
## n presences : 10
## n absences : 100
## AUC : 0.822
## cor : 0.3157062
## max TPR+TNR at : 0.2317657par(mfrow=c(1,2))
plot(ge1, 'ROC', main="ge1")
plot(ge2, 'ROC', main="ge2")
##mapeo del modelo
pg <- predict(predictors, gm2, ext=predictors)
par(mfrow=c(1,2))
plot(pg, main='GLM/gaussian, BD')
plot(BD, add=TRUE, border='dark grey')
tr <- threshold(ge2, 'spec_sens')
plot(pg>tr, main='presence/backgence')
plot(BD, add=TRUE, border='dark grey')
points(pres_train, pch='+')
dev.off()## null device
## 1Comentarios
En terminos generales el codigo esta bien mis comentarios son mas a los datos y la seleccion del modelo que predice.
Creo que es mejor si usa los datos de presencia ausencia reales en lugar de pseudo ausencias al azar. Si tiene y puede recuperar los datos de las coordenadas de los bichos con deteccion negativa es es mucho mejor.
En mi opinion es mejor que haga una preseleccion de variables para trabajar con las que no estan correlacionadas. En teoria la altitud y la temperatura lo estan, asi que indicarian lo mismo y al predecir esto introduce ruido de sobreprediccion.
Es raro que haga un modelo gausiano (normal) con datos de presencia ausencia. Es mejor uno binomial. En el binomial vale la pena probar con combinaciones cuadraticas y exponenciales de alguna variables ya que no siempre la relacion es lineal.
Si el ciriterio de seleccion del mejor modelo es el AUC hay que seleccionar el mayor AUC para predecir.
En mi opinion la ventana de prediccion es demasiado grande comparado con la extension de los puntos de presencia.
Un modelo interesante podria ser: incluir presencias y ausencias verdadras. Incluir las especies (o familias) como un fix effect en el GLM