Modelos de distribución de especies

Pseudo-ausencias (i.e. no-avistamientos)

Necesitamos generar datos de pseudo-ausencias para evaluar sitios donde no se ha encontrado tricahue. Para eso generaremos puntos aleatorios en la misma cantidad de ocurrencias descargadas de GBIF.

Luego, combinamos con los datos de variables climáticas obtenidas para generar un objeto con todas las variables climáticaspara los sitios de presencias y de preudo-ausencias.

#NORTEEEE------------------

p_load(dismo)

# setting random seed to always create the same
# random set of points for this example
set.seed(69) 

backgrN <- dismo::randomPoints(raster::stack(Bioclim), nrow(Var_ambientalesN)/2)

## Warning in dismo::randomPoints(raster::stack(Bioclim),
## nrow(Var_ambientalesN)/2): generated random points = 0.388824214202561 times
## requested number

pseudo_ausenciasN <- extract(Bioclim[[names(Var_ambientalesN)]], backgrN) # restringimos extraccion a capas bioclimaticas de interes

#presenciasN <- c(rep(1, nrow(pseudo_ausenciasN)), rep(0, nrow(Var_ambientalesN)))
#mde_dataN <- data.frame(cbind(presenciasN, rbind(Var_ambientalesN,pseudo_ausenciasN)))

presenciasN <- c(rep(1, nrow(Var_ambientalesN)),rep(0, nrow(pseudo_ausenciasN)))
mde_dataN <- data.frame(cbind(presenciasN, rbind(Var_ambientalesN,pseudo_ausenciasN)))

#SUUUURR-----------------

p_load(dismo)

# setting random seed to always create the same
# random set of points for this example
set.seed(69) 

backgrS <- dismo::randomPoints(raster::stack(Bioclim), nrow(Var_ambientalesS)/2)
pseudo_ausenciasS <- extract(Bioclim[[names(Var_ambientalesS)]], backgrS) # restringimos extraccion a capas bioclimaticas de interes

#presenciasS <- c(rep(1, nrow(pseudo_ausenciasS)), rep(0, nrow(Var_ambientalesS)))
#mde_dataS <- data.frame(cbind(presenciasS, rbind(Var_ambientalesS,pseudo_ausenciasS)))

presenciasS <- c(rep(1, nrow(Var_ambientalesS)),rep(0, nrow(pseudo_ausenciasS)))

vc <- rbind(Var_ambientalesS,pseudo_ausenciasS)

# mde_dataS <- data.frame(cbind(presenciasS, rbind(Var_ambientalesS,pseudo_ausenciasS)))
mde_dataS <- data.frame(cbind(presenciasS, vc))

Ahora que tenemos datos climáticos para nuestros puntos de presencia y pseudoausencia en el objeto mde_data, vamos a construir nuestro modelo utilizando sólo una parte de nuestros datos, y utilizaremos los datos de testeo para evaluar después el rendimiento del modelo con el set de validación.

Separamos entonces nuestros datos en un conjunto de entrenamiento, o testeo (i.e. datos utilizados para construir el modelo) y un conjunto de prueba, o validación (i.e. los datos utilizados para evaluar el modelo).

Vamos a reservar el 20% de los datos para las pruebas, por lo que utilizaremos la función folds() del paquete predicts para asignar uniformemente cada punto a un grupo aleatorio. Para asegurarnos de que tenemos una muestra más o menos representativa tanto de puntos de presencia como de pseudoausencia, utilizamos la columna presencias para indicar a R que nuestros datos tienen estos dos subgrupos.

#NORTEEEEE-----
p_load(predicts)
foldN <- folds(x = mde_dataN,
              k = 5,
              by = mde_dataN$presenciasN)
# tabla de frecuencia en cada grupo
table(foldN)

## foldN
##   1   2   3   4   5 
## 205 206 204 206 205

#SUUUUUUR------
p_load(predicts)
foldS <- folds(x = mde_dataS,
              k = 5,
              by = mde_dataS$presenciasS)

# tabla de frecuencia en cada grupo
table(foldS)

## foldS
##   1   2   3   4   5 
## 415 416 415 416 415

Dejaremos los puntos en el grupo 1 para validación y los otros grupos (i.e. 2, 3, 4 y 5) para entrenar el modelo

#NORTEEEE-----
testingN <- mde_dataN[foldN == 1, ]
trainingN <- mde_dataN[foldN != 1, ]

#SUR----
testingS <- mde_dataS[foldS == 1, ]
trainingS <- mde_dataS[foldS != 1, ]

GLM

Ajustar modelo

#NORTEEEE-----
modelo_glmN <- step(glm(presenciasN ~ ., data=trainingN, family = binomial()),trace=0)

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

summary(modelo_glmN)

## 
## Call:
## glm(formula = presenciasN ~ bio14 + bio3 + bio2 + bio15, family = binomial(), 
##     data = trainingN)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.40976   0.00032   0.00198   0.01636   1.58534  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 55.41562   15.33328   3.614 0.000301 ***
## bio14       -7.98066    2.11850  -3.767 0.000165 ***
## bio3        -0.48936    0.20658  -2.369 0.017840 *  
## bio2        -2.28741    0.63968  -3.576 0.000349 ***
## bio15        0.07345    0.02830   2.595 0.009456 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 728.484  on 814  degrees of freedom
## Residual deviance:  18.309  on 810  degrees of freedom
##   (6 observations deleted due to missingness)
## AIC: 28.309
## 
## Number of Fisher Scoring iterations: 17

#SUR----
modelo_glmS <- step(glm(presenciasS ~ ., data=trainingS, family = binomial()),trace=0)

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

summary(modelo_glmS)

## 
## Call:
## glm(formula = presenciasS ~ bio4 + bio6 + bio15 + bio8, family = binomial(), 
##     data = trainingS)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.1920  -0.0003   0.0614   0.1819   1.8106  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -39.773314   4.140240  -9.607  < 2e-16 ***
## bio4          0.063141   0.006507   9.704  < 2e-16 ***
## bio6          1.114972   0.137439   8.113 4.96e-16 ***
## bio15         0.212606   0.019668  10.810  < 2e-16 ***
## bio8         -1.039943   0.149996  -6.933 4.12e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2115.8  on 1661  degrees of freedom
## Residual deviance:  273.1  on 1657  degrees of freedom
## AIC: 283.1
## 
## Number of Fisher Scoring iterations: 10

Validación

Con la función pa_evaluate(), pasamos datos que “sabemos” cuál debería ser la respuesta correcta para estos cálculos de probabilidad. Es decir, el modelo m1 debería predecir valores cercanos a 1 para aquellas filas que pasemos al argumento p (porque sabemos que los tricahues se dan en esos lugares) y debería predecir valores cercanos a 0 para aquellas filas que pasemos al argumento a. Utilizamos esta información sobre el rendimiento del modelo para determinar el valor de probabilidad que se utilizará como límite para decir si una ubicación concreta es adecuada o inadecuada para los tricahues.

#NORTEEEE-----

# Use testing data for model evaluation
glm_evalN <- pa_evaluate(p = testingN[testingN$presenciasN == 1, ],
                        a = testingN[testingN$presenciasN == 0, ],
                        model = modelo_glmN,
                        type = "response")
#SUR----

# Use testing data for model evaluation
glm_evalS <- pa_evaluate(p = testingS[testingS$presenciasS == 1, ],
                        a = testingS[testingS$presenciasS == 0, ],
                        model = modelo_glmS,
                        type = "response")

Predicción

Con el modelo podemos generar

#NORTEEEE-----
predictionsN <- predict(Bioclim[[names(Var_ambientalesN)]], modelo_glmN, type = "response")

## 
|---------|---------|---------|---------|
=========================================

plot(predictionsN, colNA='black', main ="Norte")

writeRaster(predictionsN, "PrediccionConDatosDelNorte.tif", overwrite=TRUE)

#SUR----
predictionsS <- predict(Bioclim[[names(Var_ambientalesS)]], modelo_glmS, type = "response")

## 
|---------|---------|---------|---------|
=========================================

plot(predictionsS, colNA='black',main ="Sur")

writeRaster(predictionsS, "PrediccionConDatosDelSur.tif", overwrite=TRUE)

Evaluación

Podemos mirar el Area Bajo la Curva

#NORTEEEE-----
presence_dataN = filter(mde_dataN, presenciasN == 1)
absence_dataN = filter(mde_dataN, presenciasN == 0)

evaluationN <- evaluate(presence_dataN, absence_dataN, modelo_glmN)
tr <- dismo::threshold(evaluationN, 'spec_sens')

plot(evaluationN, 'ROC')

#SUR----
presence_dataS = filter(mde_dataS, presenciasS == 1)
absence_dataS = filter(mde_dataS, presenciasS == 0)

evaluationS <- evaluate(presence_dataS, absence_dataS, modelo_glmS)
tr <- dismo::threshold(evaluationS, 'spec_sens')

plot(evaluationS, 'ROC')

### Umbrales de distribución

Evaluación del modelo

Artículos relevantes: - Fieldings et al (1997) - Liu et al. (2011)

Validación

El estadístico de evaluación más común es el AUC: el área bajo la curva receiver-operating characteristic (ROC). Las curvas ROC se generan calculando la sensibilidad (tasa de verdaderos positivos) y la especificidad (tasa de verdaderos negativos) para muchos umbrales a lo largo de toda la gama de probabilidades previstas. Como veremos en el gráfico mas abajo, (1-especificidad) se traza en el eje de abscisas frente a la sensibilidad en el eje de ordenadas. El área bajo esta curva se denomina AUC. Cuanto más se desvíe la curva generada de la línea 1:1 hacia la esquina superior izquierda, mejor predice el modelo la presencia/ausencia de una especie. Si tomamos una presencia aleatoria y una ausencia aleatoria de nuestras observaciones y hacemos predicciones, el AUC puede interpretarse como la probabilidad de asignar una probabilidad de ocurrencia predicha más alta al punto de presencia que al de ausencia. Normalmente, consideramos que un AUC > 0,7 indica que las predicciones son correctas (Araujo et al. 2005).

El objeto glm_eval, nos muestra la calidad de la evaluación. Sin embargo, el modelo que generamos nos entrega valores continuos como predictores (i.e. la probabilidad de ocurrencia). Por eso necesitamos encontrar un umbral para decidir donde ocurre realmente la especie.

# Use testing data for model evaluation
#glm_eval <- pa_evaluate(p = testing[testing$presencias == 1, ],
#                        a = testing[testing$presencias == 0, ],
 #                       model = modelo_glm,
 #                       type = "response")

Entonces, con la función pa_evaluate(), le pasamos datos que “sabemos” cuál debería ser la respuesta correcta para estos cálculos de probabilidad. Es decir, el modelo modelo_glm debería predecir valores cercanos a 1 para aquellas filas que pasemos al argumento p (porque sabemos que los tricahues se dan en esos lugares) y debería predecir valores cercanos a 0 para aquellas filas que pasemos al argumento a. Utilizamos esta información sobre el rendimiento del modelo para determinar el valor de probabilidad que se utilizará como límite para decir si una ubicación concreta es adecuada o inadecuada para los tricahues.

El elemento thresholds de glm_eval guarda la información que nos permite determinar el umbral de corte. Aquí elegimos max_spec_sens, que establece “el umbral en el que la suma de la sensibilidad (tasa de verdaderos positivos) y la especificidad (tasa de verdaderos negativos) es más alta”. Para más información, consulte la documentación de la función pa_evaluate().

Una vez defiinido ese umbral, podemos mirar el AUC, o área bajo la curva. Como sabemos, un AUC=0.5, nos dice que el modelo se comporta igual que si hubieramos elegido variables por simple azar, cuando AUC=1.0, la predicción es perfecta (y poco creible, probablemente sobre-ajustada).

plot(glm_evalN, 'ROC')

veamos la evaluación en la proxima sección

Evaluación de modelo

Detección de umbrales

Podemos usar distintos umbrales para definir presencias.

Se define:

kappa: el umbral en el que kappa es mayor
no_omission: el umbral más alto en el que no hay omisión
prevalence : la prevalencia modelizada es la más cercana a la prevalencia observada
equal_sens_spec : equal sensitivity and specificity
max_spec_sens: el umbral en el que la suma de la sensibilidad (tasa de verdaderos positivos) y la especificidad (tasa de verdaderos negativos) es mayor

#NORTEEEE-----
glm_thresholdN <- glm_evalN@thresholds$max_spec_sens

plot(predictionsN > glm_thresholdN)

#SUR----
glm_thresholdS <- glm_evalS@thresholds$max_spec_sens

plot(predictionsS > glm_thresholdS)

Algebra de mapas en R

para comparar mapas raster usando la libreía terra.

st = predictionsN - predictionsS

plot(st,col = terrain.colors(3))

Reporte

Con estos ejercicios iremos haciendo un reporte, que será entregado el 24 de junio 2024.

Descripción de datos y análisis preliminar

Mapeo y representación gráfica

Haz un mapa de la distribución de tu C. patagus para Chile

¿En cuántas Regiones encontramos a esta especie?
Remueve los “outliers”. ¿En cuales comunas de Chile está ahora?

Análisis de variables independientes

Construye una base de datos (tabla), con los valores de temperatura, pp, radiación y variables bioclimáticas donde ocurre tu C.patagus en Chile.
1. Elige las variables mas idoneas para modelar el nicho de C. patagus. Justifica tu elección en términos de la biología de la especie.
Describe estadisticamente el espacio bioclimático en que ocurre tu especie de preferencia

Rangos de T y PP, promedio, moda, desviaciones…

Separa los puntos de ocurrencia en 2 sets que representen las poblaciones disjuntas de la especie.

Vuelve a a hacer 3 y 4

Análisis de distribución

General 2 modelos de distribución para cada grupo de ocurrencias. Esto es, usa el modelos estadísticos glm para explicar la distribución de C. patagus en su conjunto y los 2 grupos de ocurrencia seleccionado en 5

¿Cuál(es) es(son) la(s) variable(s) independiente(s) que mejor se asocian con la presencia ?
Genera 2 predicciones de distribucion, una con calibración, la del norte y la del sur para predecir sobre toda el área
Describe las diferencias. ¿En qué medida son diferentes los nichos?
¿En qué se diferencian las evaluaciones generadas por el ROC y las métricas umbral dependientes, como Kapa o prevalencias?
(bonus) Reconstruye la matriz de confusión para los 2 grupos (las ocurrencias del norte y las del sur)

Predicción

Indica cual es la predicción para la distribución de C. patagus
Discute tus resultados desde las siguientes perspectivas:

Técnicas la construcción del modelo elegido
Biológica y de conservación de tu especie de preferencia.

Nicho y MDE 3

Horacio Samaniego (horacio.samaniego@gmail.com)

2024-06-22

Modelos de distribución de especies

Pseudo-ausencias (i.e. no-avistamientos)

GLM

Ajustar modelo

Validación

Predicción

Evaluación

Evaluación del modelo

Validación

Evaluación de modelo

Detección de umbrales

Algebra de mapas en R

Reporte

Descripción de datos y análisis preliminar

Mapeo y representación gráfica

Análisis de variables independientes

Análisis de distribución

Predicción