03_GBIF_MDE.Rmd

---
title: "Nicho y MDE 3"
author: "Horacio Samaniego (horacio.samaniego@gmail.com)"
date: "`r Sys.Date()`"
output: 
  html_document: 
    toc: false
    fig_caption: true
    number_sections: false
    df_print: tibble
editor_options: 
  markdown: 
    wrap: 72
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
	echo = TRUE,
#	message = FALSE,
#	warning = FALSE,
	cache = FALSE,
	tidy = TRUE,
	tidy.opts = list(blank = FALSE, width.cutoff = 80)
)

require(pacman)
pacman::p_load(rgbif,
               rworldxtra,
               sf , 
               terra, 
               ggplot2,
               tidyverse, 
               kableExtra, 
               mapview, 
               geodata,
               ggcorrplot,
               predictions,
               dismo,
               curl)

# # For dev version
# # install.packages("devtools")
# devtools::install_github("haozhu233/kableExtra")

options("kableExtra.html.bsTable" = T)
```

Del ejercicio pasado, tenemos que importar los objetos nuevamente

```{r}
# datos 
datos = st_read("tricahue.gpkg")
# 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.

```{r}
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.

```{r}
p_load(predicts)
fold <- folds(x = mde_data,
              k = 5,
              by = mde_data$presencias)

# tabla de frecuencia en cada grupo
table(fold)
```

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

```{r}
testing <- mde_data[fold == 1, ]
training <- mde_data[fold != 1, ]
```

## GLM

### Ajustar modelo

```{r glm, echo=TRUE}

modelo_glm <- step(glm(presencias ~ ., data=training, family = binomial()),trace=0)
summary(modelo_glm)
```

### Predicción

Con el modelo podemos generar

```{r glm-prediccion, echo=TRUE}
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)](https://www.dropbox.com/scl/fi/s9a0e56osliczbl3csuhg/Fielding_Bell_1997_A-review-of-methods-for-the-assessment-of-prediction-errors-in-conservation.pdf?rlkey=n1iib2m1dcgpc7ad81ykdy3np&dl=0) -
[Liu et al. (2011)](https://doi.org/10.1111/j.1600-0587.2010.06354.x)

### 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](https://onlinelibrary.wiley.com/doi/full/10.1111/j.1365-2486.2005.01000.x)).

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.

```{r}
# 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).

```{r}


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*

```{r}

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`.

```{r}
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]

<!-- -->

a.  ¿En cuántas Regiones encontramos a esta especie? [2pts]
b.  Remueve los "outliers". ¿En cuales comunas de Chile está ahora? [3pts]

### Análisis de variables independientes

3.  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)
    a.  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]
4.  Describe estadisticamente el espacio bioclimático en que ocurre tu
    *especie de preferencia* [10pts]

<!-- -->

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

<!-- -->

a.  Vuelve a a hacer 3 y 4 [3pts]

### Análisis de distribución

6.  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 

<!-- -->

a.  ¿Cuál(es) es(son) la(s) variable(s) independiente(s) que mejor se
    asocian con la presencia ? [5pts]
b.  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]
c.  Describe las diferencias. ¿En qué medida son diferentes los nichos? [5pts]
d.  ¿En qué se diferencian las evaluaciones generadas por el ROC y las
    métricas umbral dependientes, como Kapa o prevalencias? [10pts]
e.  (bonus) Reconstruye la matriz de confusión para los 2 grupos (las
    ocurrencias del norte y las del sur) [10pts]

### Predicción

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

8.  Discute tus resultados desde las siguientes perspectivas:

<!-- -->

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