Del ejercicio pasado, tenemos que importar los objetos nuevamente

# datos
datos = st_read("tricahue.gpkg")
## Reading layer `tricahue' from data source 
##   `/mnt/data/bin/CBIT241-SDM/tricahue.gpkg' using driver `GPKG'
## Simple feature collection with 2258 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -73.83926 ymin: -40.77264 xmax: -70.0613 ymax: -23.64955
## Geodetic CRS:  WGS 84
# raster con (todo) Bioclim para el área
Bioclim = rast("Bioclim_raster.tif")
plot(Bioclim[["bio6"]])

Var_ambientales = read.csv("Var_ambientales.csv", row.names = 1)

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.

p_load(dismo)
# setting random seed to always create the same random set of points for this
# example
set.seed(69)
backgr <- dismo::randomPoints(raster::stack(Bioclim), nrow(Var_ambientales)/2)
pseudo_ausencias <- extract(Bioclim[[names(Var_ambientales)]], backgr)  # restringimos extraccion a capas bioclimaticas de interes
presencias <- c(rep(1, nrow(pseudo_ausencias)), rep(0, nrow(Var_ambientales)))
mde_data <- data.frame(cbind(presencias, rbind(Var_ambientales, pseudo_ausencias)))

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.

p_load(predicts)
fold <- folds(x = mde_data, k = 5, by = mde_data$presencias)
# tabla de frecuencia en cada grupo
table(fold)
## fold
##   1   2   3   4   5 
## 678 677 677 677 678

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

testing <- mde_data[fold == 1, ]
training <- mde_data[fold != 1, ]

GLM

Ajustar modelo

modelo_glm <- step(glm(presencias ~ ., data = training, family = binomial()), trace = 0)
summary(modelo_glm)
## 
## Call:
## glm(formula = presencias ~ bio1 + bio5 + bio9 + bio10, family = binomial(), 
##     data = training)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.69229    0.36072 -10.236  < 2e-16 ***
## bio1         0.23648    0.10132   2.334   0.0196 *  
## bio5         0.10645    0.07421   1.435   0.1514    
## bio9         0.82919    0.10940   7.580 3.46e-14 ***
## bio10       -0.95430    0.21588  -4.421 9.85e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3437.4  on 2701  degrees of freedom
## Residual deviance: 2770.2  on 2697  degrees of freedom
##   (7 observations deleted due to missingness)
## AIC: 2780.2
## 
## Number of Fisher Scoring iterations: 8

Predicción

Con el modelo podemos generar

predictions <- predict(Bioclim[[names(Var_ambientales)]], modelo_glm, type = "response")
plot(predictions, colNA = "black")

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_eval, "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
glm_threshold <- glm_eval@thresholds$max_spec_sens
plot(predictions > glm_threshold)

Algebra de mapas en R

para comparar mapas raster usando la libreía terra.

st = predictions - predictions
plot(st)

Reporte

Con estos ejercicios iremos haciendo un reporte, que será entregado el 26 de junio 2024. [75 pts totales]

Descripción de datos y análisis preliminar

Mapeo y representación gráfica

  1. Haz un mapa de la distribución de tu C. patagus para Chile [5pts]
  1. ¿En cuántas Regiones encontramos a esta especie? [2pts]
  2. Remueve los “outliers”. ¿En cuales comunas de Chile está ahora? [3pts]

Análisis de variables independientes

  1. Construye una base de datos (tabla), con los valores de variables bioclimáticas donde ocurre tu C.patagus en Chile. (opcional: puedes incluir temperatura, pp y radiación, para una descripción mas precisa del nicho)
    1. Elige y justifica 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. [5pt]
  2. Describe estadisticamente el espacio bioclimático en que ocurre tu especie de preferencia [10pts]
  1. Separa los puntos de ocurrencia en 2 sets que representen las poblaciones disjuntas de la especie. [2pts]
  1. Vuelve a a hacer 3 y 4 [3pts]

Análisis de distribución

  1. General un modelo de distribución para cada grupo de ocurrencias. 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
  1. ¿Cuál(es) es(son) la(s) variable(s) independiente(s) que mejor se asocian con la presencia ? [5pts]
  2. Genera una predicción de distribución, una con calibración, la del norte y la del sur para predecir sobre toda el área [3pts]
  3. Describe las diferencias. ¿En qué medida son diferentes los nichos? [5pts]
  4. ¿En qué se diferencian las evaluaciones generadas por el ROC y las métricas umbral dependientes, como Kapa o prevalencias? [10pts]
  5. (bonus) Reconstruye la matriz de confusión para los 2 grupos (las ocurrencias del norte y las del sur) [10pts]

Predicción

  1. Indica cual es la predicción para la distribución de C. patagus [2pts]

  2. Discute tus resultados desde las siguientes perspectivas:

  1. Técnicas la construcción del modelo elegido [10pts]
  2. Biológica y de conservación del loro Tricahue. [10pts]