class: center, top, inverse background-image: url("img/coati.jpg") background-position: 60% 60% background-size: cover # Modelos de ocupación multi-especie # 🦨🦥 Modelos de comunidad 🦝🦌 --- name: hola class: center, middle ## Gabriel Andrade Ponce 🇨🇴 <img style="border-radius: 50%;" src="img/jc_me.png" width="300px height= 300px" /> ### Estudiante de Doctorado .fade[Instituto de Ecología A.C, Xalapa, Ver., Mex] [<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg"> <path d="M459.37 151.716c.325 4.548.325 9.097.325 13.645 0 138.72-105.583 298.558-298.558 298.558-59.452 0-114.68-17.219-161.137-47.106 8.447.974 16.568 1.299 25.34 1.299 49.055 0 94.213-16.568 130.274-44.832-46.132-.975-84.792-31.188-98.112-72.772 6.498.974 12.995 1.624 19.818 1.624 9.421 0 18.843-1.3 27.614-3.573-48.081-9.747-84.143-51.98-84.143-102.985v-1.299c13.969 7.797 30.214 12.67 47.431 13.319-28.264-18.843-46.781-51.005-46.781-87.391 0-19.492 5.197-37.36 14.294-52.954 51.655 63.675 129.3 105.258 216.365 109.807-1.624-7.797-2.599-15.918-2.599-24.04 0-57.828 46.782-104.934 104.934-104.934 30.213 0 57.502 12.67 76.67 33.137 23.715-4.548 46.456-13.32 66.599-25.34-7.798 24.366-24.366 44.833-46.132 57.827 21.117-2.273 41.584-8.122 60.426-16.243-14.292 20.791-32.161 39.308-52.628 54.253z"></path></svg>@Gatorco_AP](https://twitter.com/Gatorco_AP) [<svg viewBox="0 0 496 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg"> <path d="M165.9 397.4c0 2-2.3 3.6-5.2 3.6-3.3.3-5.6-1.3-5.6-3.6 0-2 2.3-3.6 5.2-3.6 3-.3 5.6 1.3 5.6 3.6zm-31.1-4.5c-.7 2 1.3 4.3 4.3 4.9 2.6 1 5.6 0 6.2-2s-1.3-4.3-4.3-5.2c-2.6-.7-5.5.3-6.2 2.3zm44.2-1.7c-2.9.7-4.9 2.6-4.6 4.9.3 2 2.9 3.3 5.9 2.6 2.9-.7 4.9-2.6 4.6-4.6-.3-1.9-3-3.2-5.9-2.9zM244.8 8C106.1 8 0 113.3 0 252c0 110.9 69.8 205.8 169.5 239.2 12.8 2.3 17.3-5.6 17.3-12.1 0-6.2-.3-40.4-.3-61.4 0 0-70 15-84.7-29.8 0 0-11.4-29.1-27.8-36.6 0 0-22.9-15.7 1.6-15.4 0 0 24.9 2 38.6 25.8 21.9 38.6 58.6 27.5 72.9 20.9 2.3-16 8.8-27.1 16-33.7-55.9-6.2-112.3-14.3-112.3-110.5 0-27.5 7.6-41.3 23.6-58.9-2.6-6.5-11.1-33.3 2.6-67.9 20.9-6.5 69 27 69 27 20-5.6 41.5-8.5 62.8-8.5s42.8 2.9 62.8 8.5c0 0 48.1-33.6 69-27 13.7 34.7 5.2 61.4 2.6 67.9 16 17.7 25.8 31.5 25.8 58.9 0 96.5-58.9 104.2-114.8 110.5 9.2 7.9 17 22.9 17 46.4 0 33.7-.3 75.4-.3 83.6 0 6.5 4.6 14.4 17.3 12.1C428.2 457.8 496 362.9 496 252 496 113.3 383.5 8 244.8 8zM97.2 352.9c-1.3 1-1 3.3.7 5.2 1.6 1.6 3.9 2.3 5.2 1 1.3-1 1-3.3-.7-5.2-1.6-1.6-3.9-2.3-5.2-1zm-10.8-8.1c-.7 1.3.3 2.9 2.3 3.9 1.6 1 3.6.7 4.3-.7.7-1.3-.3-2.9-2.3-3.9-2-.6-3.6-.3-4.3.7zm32.4 35.6c-1.6 1.3-1 4.3 1.3 6.2 2.3 2.3 5.2 2.6 6.5 1 1.3-1.3.7-4.3-1.3-6.2-2.2-2.3-5.2-2.6-6.5-1zm-11.4-14.7c-1.6 1-1.6 3.6 0 5.9 1.6 2.3 4.3 3.3 5.6 2.3 1.6-1.3 1.6-3.9 0-6.2-1.4-2.3-4-3.3-5.6-2z"></path></svg>@gpandradep](https://github.com/gpandradep) <svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg"> <path d="M502.3 190.8c3.9-3.1 9.7-.2 9.7 4.7V400c0 26.5-21.5 48-48 48H48c-26.5 0-48-21.5-48-48V195.6c0-5 5.7-7.8 9.7-4.7 22.4 17.4 52.1 39.5 154.1 113.6 21.1 15.4 56.7 47.8 92.2 47.6 35.7.3 72-32.8 92.3-47.6 102-74.1 131.6-96.3 154-113.7zM256 320c23.2.4 56.6-29.2 73.4-41.4 132.7-96.3 142.8-104.7 173.4-128.7 5.8-4.5 9.2-11.5 9.2-18.9v-19c0-26.5-21.5-48-48-48H48C21.5 64 0 85.5 0 112v19c0 7.4 3.4 14.3 9.2 18.9 30.6 23.9 40.7 32.4 173.4 128.7 16.8 12.2 50.2 41.8 73.4 41.4z"></path></svg> gpandradep@gmail.com --- class: center # En este taller vamos a ver dos aproximaciones ## 1. Modelar la "abundancia" con modelos jerárquicos Predecimos las "abundancias" o intensidad de las especies para usarlas como insumo en los cálculos de diversidad ## 2. Usamos modelos jerárquicos para estimar riqueza Usamos directamente una proximación jerárquica para estimar la riqueza de especies de una localidad --- name: det1 # La detección imperfecta Es cuando la probabilidad de detectar una especie dado que esté presente en nuestro sitio sea menor a 1 .pull-left[ ![](img/det1.png) ] .pull-right[ ![](img/det2.png) ] --- class: midle # La detección imperfecta .pull-left[ <img src= https://zslpublications.onlinelibrary.wiley.com/cms/asset/1875594b-72d8-4b0d-bad4-b0ab9d5c7805/jzo12828-fig-0003-m.jpg width= "70%"> .footnote[Tourani, M., Brøste, E.N., Bakken, S., Odden, J. and Bischof, R. (2020), Sooner, closer, or longer: detectability of mesocarnivores at camera traps. J Zool, 312: 259-270. https://doi.org/10.1111/jzo.12828] ] .pull-right[ Cada especie tiene su propia probabilidad de ser detectada y puede depender de: - El tamaño - La velocidad de paso - La distancia a la cámara - El tiempo de muestreo - Abundancia - Entre otros.... <br>] --- class: center, middle ## Modelos de jerárquicos al rescate Son dos o más modelos probabilísticos, condicionados por las variables de probabilidad que describen. <br> ![](img/guillera2017.jpg) --- class: center ### La función que vamos a utilizar usa como base el modelo Royle-Nichols [(Royle & Nichols 2003)](https://esajournals.onlinelibrary.wiley.com/doi/10.1890/0012-9658%282003%29084%5B0777%3AEAFRPA%5D2.0.CO%3B2) .pull-left[ # Ecológico Describe la variable de estado ecológica de interés <br> En este caso describe la abundancia $$N\sim Poisson (\lambda) $$ ] .pull-right[ # Observacional Describe el proceso mediante el cual se toman los datos $$ p= 1- (1-r)^N $$ $$ y \sim Bernoulli (p)$$ ] --- class: inverse, middle, center background-image: url("img/biodiv.jpg") background-position: center background-size: cover # .big-text[ Diversidad con modelos jerárquicos] --- ## Script `DivOC_script.R` ### Comenzaremos usando la paquetería `DiversityOccupancy` Este paquete estima la diversidad alfa por medio de modelos jerárquicos. ```r # 1. Instalar y cargar el paquete ---- # install.packages("DiversityOccupancy") library(DiversityOccupancy) library(camtrapR) library(tidyverse) library(hillR) # Estimar diversidad library(ggeffects) #gráficas de pred para glm library(beepr) # Opcional para avisar R termine library(tictoc) # Opcional para tomar el tiempo de la función ``` Este paquete te va a pedir instalar también `MuMIn`, `unmarked`, `reshape`, `lattice`, `Rcpp` --- ## Formato de los datos Son necesarios eventos de muestreo repetidos. <br> Vamos a trabajar con una base de datos de 17 especies, 67 sitios (cámaras) y 19 eventos de muestreo. > **Nota**: Todas las especies deben tener una matriz de historias de detección de las mismas dimensiones ```r # Cargamos la tabla de registros de las especies registers <- read.csv("Data/Survey/recordTable_OC.csv") table(registers$Species) ``` ``` ## ## Bassariscus astutus Canis latrans Canis lupus familiaris ## 3 6 4 ## Capra hircus Conepatus leuconotus Dasypus novemcinctus ## 13 79 1 ## Lynx rufus Mephitis macroura Nasua narica ## 50 18 4 ## Odocoileus virginianus Pecari tajacu Procyon lotor ## 227 71 40 ## Puma yagouaroundi Roedores Spilogale angustigrons ## 6 66 13 ## Sylvilagus floridanus Urocyon cinereoargenteus ## 492 173 ``` --- ```r # Cargamos la tabla de operación de cámaras CToperation <- read.csv("Data/Survey/CTtable_OC.csv") # Generamos la matríz de operación de las cámaras camop <- cameraOperation(CTtable= CToperation, # Tabla de operación stationCol= "Station", # Columna que define la estación setupCol= "Setup_date", #Columna fecha de colocación retrievalCol= "Retrieval_date", #Columna fecha de retiro hasProblems= T, # Hubo fallos de cámaras dateFormat= "%Y-%m-%d") # Formato de las fechas ``` --- ```r # Función para generar las historias de detección para todas las especies seleccionadas DetHist_list <- lapply(unique(registers$Species), FUN = function(x) { detectionHistory( recordTable = registers, # Tabla de registros camOp = camop, # Matriz de operación de cámaras stationCol = "Station", speciesCol = "Species", recordDateTimeCol = "DateTimeOriginal", recordDateTimeFormat = "%d/%m/%Y", species = x, # la función reemplaza x por cada una de las especies occasionLength = 10, # Colapso de las historias a 10 días day1 = "station", #inicie en la fecha de cada estación datesAsOccasionNames = FALSE, includeEffort = TRUE, scaleEffort = TRUE, timeZone = "America/Mexico_City" )} ) # Se genera una lista con cada historia de detección y el esfuerzo de muestreo, ahora le colocaremos los nombres para saber a cual especie corresponde names(DetHist_list) <- unique(registers$Species) # Finalmente creamos una lista nueva donde estén solo las historias de detección ylist <- lapply(DetHist_list, FUN = function(x) x$detection_history) # Todas las historias deben tener el mismo número de sitios y de ocasiones de muestreo data <- ylist %>% # Lista con los datos reduce(cbind) # Unimos las historias de captura en un solo dataframe ``` --- .center[## Formato de los datos] Obtenemos una matriz con el mismo número de sitios **67** y *19* eventos X *17* especies = 323 columnas
--- # Cargar las covariables Vamos a leer el archivo .csv de la ruta `data/Covs/` que contiene las covariables de sitio (no usaremos de observación). Todas estas estandarizadas ```r # Ahora cargamos las covariables (previamente estandarizadas) covars <- read.csv("Data/Covs/stdcovs_OC.csv") %>% select(-"...1", -X, -Station, -Cam, -Cluster) # Las gráficas no funcionan con categóricas ``` <table> <thead> <tr> <th style="text-align:right;"> Vertcover_50 </th> <th style="text-align:right;"> Dcrops </th> <th style="text-align:right;"> MSAVI </th> <th style="text-align:right;"> Slope </th> <th style="text-align:right;"> Effort </th> <th style="text-align:right;"> Dpop_G </th> </tr> </thead> <tbody> <tr> <td style="text-align:right;"> -0.4744602 </td> <td style="text-align:right;"> 0.6765368 </td> <td style="text-align:right;"> 0.7316819 </td> <td style="text-align:right;"> -0.2195550 </td> <td style="text-align:right;"> -0.2080680 </td> <td style="text-align:right;"> -0.4260974 </td> </tr> <tr> <td style="text-align:right;"> 0.6519157 </td> <td style="text-align:right;"> 1.4620448 </td> <td style="text-align:right;"> -0.8666152 </td> <td style="text-align:right;"> -0.9927301 </td> <td style="text-align:right;"> -1.4605398 </td> <td style="text-align:right;"> -0.3167450 </td> </tr> <tr> <td style="text-align:right;"> -0.9243340 </td> <td style="text-align:right;"> 2.2193059 </td> <td style="text-align:right;"> -0.5838870 </td> <td style="text-align:right;"> -1.0258694 </td> <td style="text-align:right;"> -0.0719298 </td> <td style="text-align:right;"> -0.1805434 </td> </tr> <tr> <td style="text-align:right;"> -0.9243340 </td> <td style="text-align:right;"> 0.4603495 </td> <td style="text-align:right;"> -0.1430110 </td> <td style="text-align:right;"> 1.4183658 </td> <td style="text-align:right;"> -0.8887592 </td> <td style="text-align:right;"> -0.1922714 </td> </tr> </tbody> </table> > **Importante**: Cada proceso es afectado por diferentes covariables. Para más información mira este **[enlace](https://mmeredith.net/blog/2021/Polar_bear.htm?fbclid=IwAR2XKX0X9Lu-Bqu2skniXNJV54dsMqy0lb7AdFZpoCR4ZN1E5Yt71rLvgF4)** --- # Manos a la obra ### Vamos a utilizar la función `diversityoccu()` Se genera un objeto lista con los modelos para cada una de las especies, cálculo de diversidad y otras cosas.... ```r cam_diver <- diversityoccu(pres = data, # La matriz de datos sitecov = covars, # Covariables de sitio obscov = NULL, # no tenemos covariables de observación, spp = 17, # Número de especies form = ~ Effort + Slope ~ Dcrops, # Formula del modelo p(), lambda() dredge = FALSE # En este primer ejemplo no usaremos AIC ) ``` ``` ## [1] "Species 1 ready!" ## [1] "Species 2 ready!" ## [1] "Species 3 ready!" ## [1] "Species 4 ready!" ## [1] "Species 5 ready!" ## [1] "Species 6 ready!" ## [1] "Species 7 ready!" ## [1] "Species 8 ready!" ## [1] "Species 9 ready!" ## [1] "Species 10 ready!" ## [1] "Species 11 ready!" ## [1] "Species 12 ready!" ## [1] "Species 13 ready!" ## [1] "Species 14 ready!" ## [1] "Species 15 ready!" ## [1] "Species 16 ready!" ## [1] "Species 17 ready!" ``` Se va a generar un objeto lista con los modelos para cada una de las especies, cálculo de diversidad y otras cosas.... --- ### Veamos uno de los modelos 🐺 ```r cam_diver$models[[9]] # Modelo para la especie 9 ``` ``` ## ## Call: ## occuRN(formula = form, data = data2) ## ## Abundance: ## Estimate SE z P(>|z|) ## (Intercept) -0.1543 0.278 -0.555 0.579 ## Dcrops -0.0306 0.190 -0.161 0.872 ## ## Detection: ## Estimate SE z P(>|z|) ## (Intercept) -2.835 0.339 -8.37 5.54e-17 ## Effort 0.909 0.324 2.81 4.99e-03 ## Slope -0.596 0.216 -2.76 5.78e-03 ## ## AIC: 317.6726 ``` --- class: inverse, center # El problema..... No sabemos si todas las especies responden de la misma manera a las covariables que usamos. Debemos escoger de todas las variables cual se ajusta mejor a cada especie .pull-left[ ### ¿Cómo vamos a generar todas las posibles combinaciones de modelos para cada especie? ] -- .pull-right[ ![](https://blogs.unitec.mx/content/dam/blogs/imagenes/corp_samara/por-que-es-tan-dificil-quedarte-en-tu-primera-opcion-7.gif) Tranquilo esta misma función lo hace por ti :D ] --- ## Solamente tenemos que activar `dredge` Dependiendo de tu computador la función puede tardar más o menos. En la mía duró ~57 segs. Pero puede tardar mucho más dependiendo del número de especies, sitios, eventos de muestreo y cantidad de covariables. ```r cam_diver_AIC <- diversityoccu(pres = data, # La matriz de datos sitecov = covars, # Covariables de sitio obscov = NULL, # no tenemos covariables de observación, spp = 17, # Número de especies form = ~ Effort + Slope ~ Dcrops, # Formula del modelo p(), lambda() dredge = TRUE # escoge los mejores modelos con AIC ) ``` --- class: center ## Veamos de nuevo el modelo que seleccionó para la sp9 🐺 .pull-left[ ```r cam_diver$models[[9]] ``` ``` ## ## Call: ## occuRN(formula = form, data = data2) ## ## Abundance: ## Estimate SE z P(>|z|) ## (Intercept) -0.1543 0.278 -0.555 0.579 ## Dcrops -0.0306 0.190 -0.161 0.872 ## ## Detection: ## Estimate SE z P(>|z|) ## (Intercept) -2.835 0.339 -8.37 5.54e-17 ## Effort 0.909 0.324 2.81 4.99e-03 ## Slope -0.596 0.216 -2.76 5.78e-03 ## ## AIC: 317.6726 ``` ] .pull-right[ ```r cam_diver_AIC$models[[9]] ``` ``` ## ## Call: ## occuRN(formula = ~Effort + Slope + 1 ~ 1, data = data2) ## ## Abundance: ## Estimate SE z P(>|z|) ## -0.145 0.272 -0.532 0.595 ## ## Detection: ## Estimate SE z P(>|z|) ## (Intercept) -2.846 0.333 -8.54 1.30e-17 ## Effort 0.925 0.308 3.01 2.62e-03 ## Slope -0.596 0.215 -2.77 5.58e-03 ## ## AIC: 315.6987 ``` ] --- class: center ### Veamos el resultado gráfico para la especie 3 🐗 ```r (responseplot.abund( batch = cam_diver_AIC, # objeto creado con diversityoccu spp = 3, # número o nombre de la sp variable= Dcrops # variable )) ``` <img src="multi_sp_GAP_files/figure-html/unnamed-chunk-14-1.png" width="100%" /> --- class: center # Momento de la diversidad Tenemos un modelo donde se estima la abundancia para cada especie, es hora de modelar la diversidad ```r glm.div <- model.diversity(DivOcc = cam_diver_AIC,# modelos method = "h", # método delta = 2, squared = T # términos cuadráticos ) ``` --- class: center # Momento de la diversidad ```r AICtab <- glm.div$Table ``` <table> <thead> <tr> <th style="text-align:left;"> model </th> <th style="text-align:right;"> aicc </th> <th style="text-align:right;"> weights </th> <th style="text-align:right;"> Delta.AICc </th> </tr> </thead> <tbody> <tr> <td style="text-align:left;"> Diversity ~ 1 + Dcrops + I(Dcrops^2) </td> <td style="text-align:right;"> -331.8899 </td> <td style="text-align:right;"> 0.3229081 </td> <td style="text-align:right;"> 0.0000000 </td> </tr> <tr> <td style="text-align:left;"> Diversity ~ 1 + Vertcover_50 + Dcrops + I(Dcrops^2) </td> <td style="text-align:right;"> -331.3285 </td> <td style="text-align:right;"> 0.2438852 </td> <td style="text-align:right;"> 0.5613404 </td> </tr> <tr> <td style="text-align:left;"> Diversity ~ 1 + Dcrops + MSAVI + I(Dcrops^2) </td> <td style="text-align:right;"> -330.5617 </td> <td style="text-align:right;"> 0.1662141 </td> <td style="text-align:right;"> 1.3281821 </td> </tr> <tr> <td style="text-align:left;"> Diversity ~ 1 + Dcrops + MSAVI + I(Dcrops^2) + I(MSAVI^2) </td> <td style="text-align:right;"> -330.3166 </td> <td style="text-align:right;"> 0.1470432 </td> <td style="text-align:right;"> 1.5732824 </td> </tr> <tr> <td style="text-align:left;"> Diversity ~ 1 + Dcrops + Dpop_G + I(Dcrops^2) </td> <td style="text-align:right;"> -329.9093 </td> <td style="text-align:right;"> 0.1199495 </td> <td style="text-align:right;"> 1.9805941 </td> </tr> </tbody> </table> --- ### Ahora veamos la respuesta gráfica de la diversidad a una variable ```r responseplot.diver(glm.div, Dcrops) ``` <img src="multi_sp_GAP_files/figure-html/unnamed-chunk-19-1.png" width="100%" /> A medida que aumenta el valor de distancia a cultivos (escalado) hay mayor diversidad ????? --- class: inverse, center, middle # Relativamente fácil ![](https://ugc.kn3.net/i/760x/https://media1.tenor.com/images/048cc35b45c1545328f41a03d0177ffa/tenor.gif?itemid=5513707) # Para ser verdad ..... --- # Diversidad.... Si seguimos la viñeta del paquete nunca nos dice que elemento de la diversidad mide o calcula.. - Riqueza? - Shannon, simpson?.... - Número efectivo de especies?? (Hill) -- # Glm? de que tipo? Un Glm puede ser de varias familias (distribuciones) y depende de la naturaleza de los datos: conteos, proporciones, unos y ceros --- class: center # Es importante leer el manual .pull-left[ ![](https://i.gifer.com/7hRo.gif) ] .pull-right[ Hay otro argumento de la función `diversityoccu` y es *"index"*. Este argumento permite escoger que índice utilizar **"shannon"**, **"simpson"** o **"invsimpson"**. ```r cam_diver_sh <- diversityoccu(pres = data, sitecov = covars, obscov = NULL, spp = 17, form = ~ Effort + Slope ~ Dcrops, dredge = TRUE, * index = "shannon" ) ``` ] --- ## Podemos aplicar la función para cada índice...? > También se puede hacer una función que lo haga en automático, pero por simplicidad (~~no se hacerlo bien~~) corremos tres veces la función ```r cam_diver_sim <- diversityoccu(pres = data, sitecov = covars, obscov = NULL, spp = 17, form = ~ Effort + Slope ~ Dcrops, dredge = TRUE, * index = "simpson" ) cam_diver_inv <- diversityoccu(pres = data, sitecov = covars, obscov = NULL, spp = 17, form = ~ Effort + Slope ~ Dcrops, dredge = TRUE, * index = "invsimpson" ) ``` --- # Índices de entropía ..... El problema de estos índices (~~para los ecólogos~~) es que - son adimensionales - no siguen una relación lineal con la riqueza (doble de riqueza `\(\neq\)` doble de diversidad) ![](img/jost2006.png) --- ### Calculemos el número efectivo de especies con las abundancias estimadas ```r # Extraer los datos de abundancia hill_data <- cam_diver_inv[[4]] %>% select(-h) ``` Sospechoso... <table> <thead> <tr> <th style="text-align:right;"> species.1 </th> <th style="text-align:right;"> species.2 </th> <th style="text-align:right;"> species.3 </th> <th style="text-align:right;"> species.4 </th> <th style="text-align:right;"> species.5 </th> <th style="text-align:right;"> species.6 </th> <th style="text-align:right;"> species.7 </th> <th style="text-align:right;"> species.8 </th> <th style="text-align:right;"> species.9 </th> <th style="text-align:right;"> species.10 </th> <th style="text-align:right;"> species.11 </th> <th style="text-align:right;"> species.12 </th> <th style="text-align:right;"> species.13 </th> <th style="text-align:right;"> species.14 </th> <th style="text-align:right;"> species.15 </th> <th style="text-align:right;"> species.16 </th> <th style="text-align:right;"> species.17 </th> </tr> </thead> <tbody> <tr> <td style="text-align:right;"> 0.6082829 </td> <td style="text-align:right;"> 1.179013 </td> <td style="text-align:right;"> 0.9865411 </td> <td style="text-align:right;"> 0.7986634 </td> <td style="text-align:right;"> 0.6344764 </td> <td style="text-align:right;"> 0.470707 </td> <td style="text-align:right;"> 1.73449 </td> <td style="text-align:right;"> 2.942110 </td> <td style="text-align:right;"> 0.8654214 </td> <td style="text-align:right;"> 0.1201388 </td> <td style="text-align:right;"> 8.197757 </td> <td style="text-align:right;"> 0.1485776 </td> <td style="text-align:right;"> 8.906907 </td> <td style="text-align:right;"> 0.0098931 </td> <td style="text-align:right;"> 0.119418 </td> <td style="text-align:right;"> 11.42479 </td> <td style="text-align:right;"> 0.237951 </td> </tr> <tr> <td style="text-align:right;"> 0.6082829 </td> <td style="text-align:right;"> 1.179013 </td> <td style="text-align:right;"> 1.2351996 </td> <td style="text-align:right;"> 0.7986634 </td> <td style="text-align:right;"> 0.6344764 </td> <td style="text-align:right;"> 0.470707 </td> <td style="text-align:right;"> 1.73449 </td> <td style="text-align:right;"> 6.441039 </td> <td style="text-align:right;"> 0.8654214 </td> <td style="text-align:right;"> 0.1201388 </td> <td style="text-align:right;"> 8.197757 </td> <td style="text-align:right;"> 0.1485776 </td> <td style="text-align:right;"> 8.906907 </td> <td style="text-align:right;"> 0.0014365 </td> <td style="text-align:right;"> 0.119418 </td> <td style="text-align:right;"> 11.42479 </td> <td style="text-align:right;"> 0.237951 </td> </tr> <tr> <td style="text-align:right;"> 0.6082829 </td> <td style="text-align:right;"> 1.179013 </td> <td style="text-align:right;"> 1.5340820 </td> <td style="text-align:right;"> 0.7986634 </td> <td style="text-align:right;"> 0.6344764 </td> <td style="text-align:right;"> 0.470707 </td> <td style="text-align:right;"> 1.73449 </td> <td style="text-align:right;"> 13.709316 </td> <td style="text-align:right;"> 0.8654214 </td> <td style="text-align:right;"> 0.1201388 </td> <td style="text-align:right;"> 8.197757 </td> <td style="text-align:right;"> 0.1485776 </td> <td style="text-align:right;"> 8.906907 </td> <td style="text-align:right;"> 0.0002236 </td> <td style="text-align:right;"> 0.119418 </td> <td style="text-align:right;"> 11.42479 </td> <td style="text-align:right;"> 0.237951 </td> </tr> <tr> <td style="text-align:right;"> 0.6082829 </td> <td style="text-align:right;"> 1.179013 </td> <td style="text-align:right;"> 0.9273586 </td> <td style="text-align:right;"> 0.7986634 </td> <td style="text-align:right;"> 0.6344764 </td> <td style="text-align:right;"> 0.470707 </td> <td style="text-align:right;"> 1.73449 </td> <td style="text-align:right;"> 2.371387 </td> <td style="text-align:right;"> 0.8654214 </td> <td style="text-align:right;"> 0.1201388 </td> <td style="text-align:right;"> 8.197757 </td> <td style="text-align:right;"> 0.1485776 </td> <td style="text-align:right;"> 8.906907 </td> <td style="text-align:right;"> 0.0168258 </td> <td style="text-align:right;"> 0.119418 </td> <td style="text-align:right;"> 11.42479 </td> <td style="text-align:right;"> 0.237951 </td> </tr> <tr> <td style="text-align:right;"> 0.6082829 </td> <td style="text-align:right;"> 1.179013 </td> <td style="text-align:right;"> 1.1369200 </td> <td style="text-align:right;"> 0.7986634 </td> <td style="text-align:right;"> 0.6344764 </td> <td style="text-align:right;"> 0.470707 </td> <td style="text-align:right;"> 1.73449 </td> <td style="text-align:right;"> 4.824359 </td> <td style="text-align:right;"> 0.8654214 </td> <td style="text-align:right;"> 0.1201388 </td> <td style="text-align:right;"> 8.197757 </td> <td style="text-align:right;"> 0.1485776 </td> <td style="text-align:right;"> 8.906907 </td> <td style="text-align:right;"> 0.0029269 </td> <td style="text-align:right;"> 0.119418 </td> <td style="text-align:right;"> 11.42479 </td> <td style="text-align:right;"> 0.237951 </td> </tr> <tr> <td style="text-align:right;"> 0.6082829 </td> <td style="text-align:right;"> 1.179013 </td> <td style="text-align:right;"> 1.4022314 </td> <td style="text-align:right;"> 0.7986634 </td> <td style="text-align:right;"> 0.6344764 </td> <td style="text-align:right;"> 0.470707 </td> <td style="text-align:right;"> 1.73449 </td> <td style="text-align:right;"> 10.022276 </td> <td style="text-align:right;"> 0.8654214 </td> <td style="text-align:right;"> 0.1201388 </td> <td style="text-align:right;"> 8.197757 </td> <td style="text-align:right;"> 0.1485776 </td> <td style="text-align:right;"> 8.906907 </td> <td style="text-align:right;"> 0.0004836 </td> <td style="text-align:right;"> 0.119418 </td> <td style="text-align:right;"> 11.42479 </td> <td style="text-align:right;"> 0.237951 </td> </tr> </tbody> </table> --- # Calculemos diversidad con `hillR` ```r # calcular los perfiles de diversidad q0 <- hill_taxa(hill_data, q=0) q1 <- hill_taxa(hill_data, q=1) q2 <- hill_taxa(hill_data, q=2) ``` Ahora unimos las bases de datos y las covariables para modelar ```r # Unir las bases de perfiles de diversidad hill_div <- data.frame(q0=q0, q1=q1, q2=q2) # Unir con las covariables glm_hill <- cbind(hill_div, covars) ``` --- # Obtenemos esta base <table> <thead> <tr> <th style="text-align:right;"> q0 </th> <th style="text-align:right;"> q1 </th> <th style="text-align:right;"> q2 </th> <th style="text-align:right;"> Vertcover_50 </th> <th style="text-align:right;"> Dcrops </th> <th style="text-align:right;"> MSAVI </th> <th style="text-align:right;"> Slope </th> <th style="text-align:right;"> Effort </th> <th style="text-align:right;"> Dpop_G </th> </tr> </thead> <tbody> <tr> <td style="text-align:right;"> 17 </td> <td style="text-align:right;"> 7.290641 </td> <td style="text-align:right;"> 5.283733 </td> <td style="text-align:right;"> -0.4744602 </td> <td style="text-align:right;"> 0.6765368 </td> <td style="text-align:right;"> 0.7316819 </td> <td style="text-align:right;"> -0.2195550 </td> <td style="text-align:right;"> -0.2080680 </td> <td style="text-align:right;"> -0.4260974 </td> </tr> <tr> <td style="text-align:right;"> 17 </td> <td style="text-align:right;"> 7.476420 </td> <td style="text-align:right;"> 5.687832 </td> <td style="text-align:right;"> 0.6519157 </td> <td style="text-align:right;"> 1.4620448 </td> <td style="text-align:right;"> -0.8666152 </td> <td style="text-align:right;"> -0.9927301 </td> <td style="text-align:right;"> -1.4605398 </td> <td style="text-align:right;"> -0.3167450 </td> </tr> <tr> <td style="text-align:right;"> 17 </td> <td style="text-align:right;"> 7.068294 </td> <td style="text-align:right;"> 5.418051 </td> <td style="text-align:right;"> -0.9243340 </td> <td style="text-align:right;"> 2.2193059 </td> <td style="text-align:right;"> -0.5838870 </td> <td style="text-align:right;"> -1.0258694 </td> <td style="text-align:right;"> -0.0719298 </td> <td style="text-align:right;"> -0.1805434 </td> </tr> <tr> <td style="text-align:right;"> 17 </td> <td style="text-align:right;"> 7.204272 </td> <td style="text-align:right;"> 5.173335 </td> <td style="text-align:right;"> -0.9243340 </td> <td style="text-align:right;"> 0.4603495 </td> <td style="text-align:right;"> -0.1430110 </td> <td style="text-align:right;"> 1.4183658 </td> <td style="text-align:right;"> -0.8887592 </td> <td style="text-align:right;"> -0.1922714 </td> </tr> <tr> <td style="text-align:right;"> 17 </td> <td style="text-align:right;"> 7.448817 </td> <td style="text-align:right;"> 5.558398 </td> <td style="text-align:right;"> -1.1509622 </td> <td style="text-align:right;"> 1.1723154 </td> <td style="text-align:right;"> 0.2129457 </td> <td style="text-align:right;"> -0.0597006 </td> <td style="text-align:right;"> -0.0719298 </td> <td style="text-align:right;"> -0.0948151 </td> </tr> <tr> <td style="text-align:right;"> 17 </td> <td style="text-align:right;"> 7.342298 </td> <td style="text-align:right;"> 5.686248 </td> <td style="text-align:right;"> 0.6519157 </td> <td style="text-align:right;"> 1.9052628 </td> <td style="text-align:right;"> 0.5284530 </td> <td style="text-align:right;"> -1.0258694 </td> <td style="text-align:right;"> -0.0719298 </td> <td style="text-align:right;"> 0.0268169 </td> </tr> </tbody> </table> --- # Ahora ajustemos un glm sencillo ```r glm_q1 <- glm(q1~ Dcrops, family = gaussian, data = glm_hill) ``` y usemos ggeffects para graficar .pull-left[ ```r plot_q1 <- ggpredict(glm_q1, terms = "Dpop_G") plot(plot_q1)+ labs(y= "Diversidad q1", x= "Distancia a poblados (estandarizado)")+ theme_classic() ``` ] .pull-right[ <img src="multi_sp_GAP_files/figure-html/unnamed-chunk-30-1.png" width="100%" /> ] --- class: center # Tarea - Modelar los otros perfiles de diversidad - Crear modelos candidatos con otras covariables y escoger con AIC ![](https://www.recreoviral.com/wp-content/uploads/2016/04/Las-tareas-est%C3%A1n-destruyendo-a-nuestros-hijos-y-los-cient%C3%ADficos-piden-prohibirla-4.gif) --- .pull-left[ # Consideraciones finales 1. El paquete es bueno y agiliza muchos pasos del modelado 2. Cuando no hay ninguna covariable que explique la abundancia, se asume que es constante para todos los sitios. Eso puede subestimar o sobrestimar la abundancia... 3. Abundancia constante: diversidad de la localidad y no de la cámara? 4. Si se usan los índices de entropía se debe saber cómo interpretarlos ([Jost 2006](https://onlinelibrary.wiley.com/doi/abs/10.1111/j.2006.0030-1299.14714.x) ; [Jost et al 2010](https://onlinelibrary.wiley.com/doi/full/10.1111/j.1472-4642.2009.00626.x))] .pull-right[ <img src="https://compote.slate.com/images/18065fdf-67ca-448a-9ebd-6236a7dac8e8.jpg" /> ] --- ## Un ejemplo de cómo se pueden usar los modelos Royle-Nichols para calcular diversidad .pull-left[ > Tesis de pregrado premiada por la AMMAC https://www.youtube.com/watch?v=qaD9NRAg3SQ ] .pull-right[ <iframe width="560" height="315" src="https://www.youtube.com/embed/qaD9NRAg3SQ" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>] --- class: inverse, center, middle <img src="https://pa1.narvii.com/6159/c7f0dd0c6a7a81293a754044edc52cd4ee80e310_hq.gif" width="50%" align="right" /> .pull-left[ # Descanso para un cafecito
10
:
00
] --- class: inverse, center # Otra opción elegante .pull-left[ <br> ### Es posible estimar la riqueza de especies directamente de un modelo jerárquico Es mejor modelar directamente la variación espacial de la riqueza mientras se considera la detección, eliminando así varias etapas de análisis ] .pull-right[ ![](img/memediv.jpg) ] --- # Modelo de comunidad Dorazio & Royle (2005) Es una extensión del modelo de una especie y una temporada en donde se combinan las historias de detección de todas las especies encontradas en los sitios de muestreo ![](img/multi.png) --- class: center # ¿Cómo funciona? .pull-left[ ![](img/M-N.png) ] .pull-right[ <br> <br> Para entender el modelo tenemos que imaginar una comunidad con **N** especies. Al realizar un muestreo encontramos **n** especies. Cómo no sabemos que porcentaje de **n** representa **N**, debemos asumir una comunidad **M** lo suficientemente grande para contener a **N**.] --- # El aumentado de datos 🪄 .pull-left[ Es un concepto utilizado bajo el enfoque bayesiano y consiste en agregar un número arbitrario de especies potenciales a la muestra de modo que `\(M>N\)`. <br> <br> Este número se define como `\(nz= M-n\)`, donde **n** es el número de especies observado y **nz** es el número de especies potenciales adicionadas. ] .pull-right[ <br> ![](https://i.pinimg.com/236x/09/1b/50/091b50565f23ae2be2c4c2b533d1bd33--ugly-americans-nerd.jpg) ] --- # Este modelo tiene tres niveles o procesos ### 1. Primern nivel o super parámetro El parámetro aumentado describe la probabilidad de que la especie "real" **k** pertenezca a **N**, dado el total de especies potenciales de **M** $$ \omega{_k} \sim Bernoulli(\Omega) $$ -- Donde `\(\Omega\)` es la probabilidad de que la especie **k** pertenezca a la meta-comunidad de tamaño (desconocido) **N**. En ese sentido `\(N=M\Omega\)`, por lo cual al estimar `\(\Omega\)` conoceremos **N**. Esto es posible ya que conocemos **M**. --- ### 2. Segundo nivel o proceso ecológico Describe la ocupación de las especies en cada sitio *i* siguiendo una distribución *Bernoulli* `$$z_{ik}|w_k\sim Bernoulli(\psi_{k}\omega_k)$$` -- Donde `\(z_{k}\)` son elementos de la matriz de ocupación real y `\(\psi_{k}\)` es la probabilidad de ocupación de cada especie en cada sitio de muestreo dado que este presente en la comunidad. -- ### 3. Tercer nivel o proceso observacional El proceso de detección o logístico es similar al del modelo de ocupación básico. Puede ser *Binomial* o *Bernoulli* dependiendo del formato de los datos `$$ysum_{ik}|z_{ik}\sim Binomial(J_i,p_{k}z_{ik})$$` -- Donde `\(ysum_{ik}\)` indica la detección o no detección de una especie **k** en el sitio **i** en **J** número de eventos para cada sitio, con una probabilidad de `\(p_{k}\)` --- # El modelo asume heterogeneidad en psi y p El modelo asume a cada especie como un **efecto aleatorio**. Es decir, que todas provienen de una misma distribución y son similares pero no iguales. Este supuesto hace que especies con menos datos no difumine la información de la comunidad. Lo que proporciona mayor precisión y habilidad de predicción al parámetro aumentado. En términos de probabilidad los efectos aleatorios se define así para la ocupación `$$logit(\psi_k)\sim Normal(\mu_{lpsi},\sigma^2_{lpsi})$$` Y así para la detección `$$logit(p_k)\sim Normal(\mu_{lp},\sigma^2_{lp})$$` Donde `\(\mu\)` y `\(\sigma\)` son los parámetros de una distribución normal --- # ¿Cómo se estiman las especies faltantes? El modelo usa la información de las especies presentes para estimar las especies sin datos. Así la probabilidad de no encontrar una especie **k** es: `$$m_k=[1-\psi_{ik}[1-(1-p_k)^J]]^S$$` --- class: center ### ¿Cómo se estiman las especies faltantes? A medida que tenemos menor `\(\psi_k\)` o `\(p_k\)` pues `\(N>n\)` <img src="img/mk.png" width="1200" height= "400" /> --- # Manos a la obra ### Script `multioccu_script.R` Recuerden instalar la versión de desarrollo de camtrapR ```r remotes::install_github("jniedballa/camtrapR") ``` Lo primero es cargar todas las librerias necesarias. > Recuerden que para instalar rjags necesitan instalar en su maquina el programa [JAGS](https://sourceforge.net/projects/mcmc-jags/files/) ```r library(camtrapR) # Datos de cámaras y modelos library(rjags) # Para correr el modelo library(SpadeR) # Riqueza Chao2 library(tidyverse) # Manipular datos library(nimble) # LEnguaje BUGS library(nimbleEcology) # Nimble enfocado en jerárquicos library(bayesplot) # gráficos estimaciones bayesianas library(SpadeR) # Riqueza Chao2 library(beepr) # Opcional para avisar R termine library(tictoc) # Opcional para tomar el tiempo de la función library(extrafont) #opcional para cambiar la fuente library(snowfall) ``` --- # Formato de los datos Vamos a usar los mismos datos que del ejemplo anterior. Seguimos un procedimiento similar para cargar los datos. ```r # Cargamos la tabla de registros de las especies registers <- read.csv("Data/Survey/recordTable_OC.csv") # Cargamos la tabla de operación de cámaras CToperation <- read.csv("Data/Survey/CTtable_OC.csv") # Generamos la matríz de operación de las cámaras camop <- cameraOperation(CTtable= CToperation, # Tabla de operación stationCol= "Station", # Columna que define la estación setupCol= "Setup_date", #Columna fecha de colocación retrievalCol= "Retrieval_date", #Columna fecha de retiro hasProblems= T, # Hubo fallos de cámaras dateFormat= "%Y-%m-%d") # Formato de las fechas ``` --- ```r # Función para generar las historias de detección para todas las especies seleccionadas DetHist_list <- lapply(unique(registers$Species), FUN = function(x) { detectionHistory( recordTable = registers, # Tabla de registros camOp = camop, # Matriz de operación de cámaras stationCol = "Station", speciesCol = "Species", recordDateTimeCol = "DateTimeOriginal", recordDateTimeFormat = "%d/%m/%Y", species = x, # la función reemplaza x por cada una de las especies occasionLength = 10, # Colapso de las historias a 10 ías day1 = "station", #inicie en la fecha de cada estación datesAsOccasionNames = FALSE, includeEffort = TRUE, scaleEffort = TRUE, timeZone = "America/Mexico_City" )} ) # Se genera una lista con cada historia de detección y el esfuerzo de muestreo, ahora le colocaremos los nombres para saber a cual especie corresponde names(DetHist_list) <- unique(registers$Species) # Finalmente creamos una lista nueva donde estén solo las historias de detección ylist <- lapply(DetHist_list, FUN = function(x) x$detection_history) ``` Terminaremos con el objeto ylist que cotiente todas las historias de detección --- ### Finalmente las covariables ```r ## Covariables #Cargamos la base de covariables covars <- read.csv("Data/Covs/stdcovs_OC.csv") identical(nrow(ylist[[1]]), nrow(covars)) # Verificar que tengan el mismo número de filas ``` ``` ## [1] TRUE ``` --- ## Finalmente unimos todo en un objeto lista con los datos que requiere el modelo ```r # Generamos la base de datos para el modelo data_list <- list(ylist = ylist, # Historias de detección siteCovs = covars, # Covariables de sitio obsCovs = list(effort = DetHist_list[[1]]$effort)) # agregamos el esfuerzo de muestreo como covariable de observación str(data_list) ``` ``` ## List of 3 ## $ ylist :List of 17 ## ..$ Conepatus leuconotus : num [1:67, 1:19] 0 NA 0 0 1 0 0 0 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Lynx rufus : num [1:67, 1:19] 0 NA 0 0 0 0 0 0 1 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Sylvilagus floridanus : num [1:67, 1:19] 1 NA 0 0 1 0 0 0 1 1 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Urocyon cinereoargenteus: num [1:67, 1:19] 1 NA 0 0 1 0 0 0 1 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Canis latrans : num [1:67, 1:19] 0 NA 0 0 0 0 0 0 1 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Mephitis macroura : num [1:67, 1:19] 0 NA 0 0 0 0 0 0 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Odocoileus virginianus : num [1:67, 1:19] 0 NA 1 0 1 0 0 0 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Canis lupus familiaris : num [1:67, 1:19] 0 NA 0 0 0 0 0 0 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Pecari tajacu : num [1:67, 1:19] 0 NA 0 0 0 0 0 0 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Roedores : num [1:67, 1:19] 0 NA 0 0 1 0 0 0 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Dasypus novemcinctus : num [1:67, 1:19] 0 NA 0 0 0 0 0 0 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Spilogale angustigrons : num [1:67, 1:19] 0 NA 0 0 0 0 0 0 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Bassariscus astutus : num [1:67, 1:19] 0 NA 0 0 0 0 0 0 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Procyon lotor : num [1:67, 1:19] 0 NA 0 0 0 0 0 0 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Capra hircus : num [1:67, 1:19] 0 NA 0 0 0 0 0 0 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Puma yagouaroundi : num [1:67, 1:19] 0 NA 0 0 0 0 0 0 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## ..$ Nasua narica : num [1:67, 1:19] 0 NA 0 0 0 0 0 0 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ## $ siteCovs:'data.frame': 67 obs. of 11 variables: ## ..$ X : int [1:67] 1 2 3 4 5 6 7 8 9 10 ... ## ..$ Station : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## ..$ Cam : chr [1:67] "MoultrieA30" "Moultrie" "MoultrieA30" "Primos" ... ## ..$ ...1 : num [1:67] -1.69 -1.64 -1.59 -1.54 -1.49 ... ## ..$ Vertcover_50: num [1:67] -0.474 0.652 -0.924 -0.924 -1.151 ... ## ..$ Dcrops : num [1:67] 0.677 1.462 2.219 0.46 1.172 ... ## ..$ MSAVI : num [1:67] 0.732 -0.867 -0.584 -0.143 0.213 ... ## ..$ Slope : num [1:67] -0.2196 -0.9927 -1.0259 1.4184 -0.0597 ... ## ..$ Cluster : num [1:67] -1.06 -1.06 -1.06 -1.06 -1.06 ... ## ..$ Effort : num [1:67] -0.2081 -1.4605 -0.0719 -0.8888 -0.0719 ... ## ..$ Dpop_G : num [1:67] -0.4261 -0.3167 -0.1805 -0.1923 -0.0948 ... ## $ obsCovs :List of 1 ## ..$ effort: num [1:67, 1:19] 0.0626 NA 0.0626 0.0626 0.0626 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:67] "C1T2P1" "C1T2P11" "C1T2P21" "C1T4P1" ... ## .. .. ..$ : chr [1:19] "o1" "o2" "o3" "o4" ... ``` --- ## Creando el modelo CamtrapR permite ajustar modelos multi-especie en JAGS y Nimble (es decir en lenguaje BUGS), nosotros vamos a usar JAGS ya que la versión de Nimble aun no permite estimar parámetro N de riqueza de especies .pull-left[ ```r # Usaremos la función ` communityModel` # Generemos el modelo comu_model <- communityModel(data_list, # la lista de datos occuCovs = list(ranef = "Dcrops"), # La covariables de sitio detCovsObservation = list(fixed = "effort"), #Covariables de observación intercepts = list(det = "ranef", occu = "ranef"), augmentation = c(full = 30),# Número aumentado de especies modelFile = "multmod")# Guardamos la especificación en un archivo ``` ] .pull-right[ ```r summary(comu_model) ``` ``` ## commOccu object (for JAGS) ## ## 30 species, 67 stations, 19 occasions ## 758 occasions with effort ## Number of detections (by species): 0 - 146 ## ## Available site covariates: ## X, Station, Cam, ...1, Vertcover_50, Dcrops, MSAVI, Slope, Cluster, Effort, Dpop_G ## ## Used site covariates: ## Dcrops ## ## Available site-occasion covariates: ## effort ``` ] --- ## Hora de correr el modelo En este caso no corran la función porque dura alrededor de 56 min ```r fit.commu <- fit(comu_model, n.iter = 22000, n.burnin = 2000, thin = 2, chains = 3, cores = 3, quiet = T );beep(sound = 4) # Duración 56 min aprox ``` --- ## Nuestro resultado ```r modresult <- summary(fit.commu)[["statistics"]] ```
--- ### Veamos el resultado gráfico .pull-left[ Otra gran ventaja de CamtrapR es que permite gráficar de manera muy sencilla la predicción posterior del modelo. Veamos que pasa con la ocupación de cada especie ```r plot_effects(comu_model, # El modelo fit.commu, # El objeto ajustado submodel = "state") # el parámetro de interés ``` ] .pull-right[ ``` ## $Dcrops ``` <img src="multi_sp_GAP_files/figure-html/unnamed-chunk-45-1.png" width="100%" /> ] --- ### Ahora los coeficientes ```r plot_coef(comu_model, fit.commu, submodel = "state") ``` ``` ## $Dcrops ``` <img src="multi_sp_GAP_files/figure-html/unnamed-chunk-46-1.png" width="100%" /> --- ### Realizamos el mismo procedimiento para el submodelo de detección ```r plot_effects(comu_model, fit.commu, submodel = "det") ``` ``` ## $effort ``` <img src="multi_sp_GAP_files/figure-html/unnamed-chunk-47-1.png" width="100%" /> --- ### y sus respectivos coeficientes ```r plot_coef(comu_model, fit.commu, submodel = "det") ``` ``` ## $effort ``` <img src="multi_sp_GAP_files/figure-html/unnamed-chunk-48-1.png" width="100%" /> --- ### Lo que vinimos a buscar fue la riqueza de especies La riqueza estimada de especies es ~20 sps .pull-left[ ```r (riqueza_est <- modresult["Ntotal",]) ``` ``` ## Mean SD Naive SE Time-series SE ## 19.96107056 2.95338837 0.01705054 0.08235277 ``` ```r # Veamos el gráfico de la distribución posterior mcmc_areas(fit.commu, # objeto jags pars= "Ntotal", # parámetro de interés point_est = "mean", prob = 0.95) # intervalos de credibilidad ``` ] .pull-right[ <img src="multi_sp_GAP_files/figure-html/unnamed-chunk-51-1.png" width="100%" /> ] --- ### La estimación no se ve muy bien, hay que verificar los trace plots Debería verse como un cesped, muy probablemente necesitamos muchas mas iteraciones para este modelo ```r mcmc_trace(fit.commu, pars = "Ntotal") ``` <img src="multi_sp_GAP_files/figure-html/unnamed-chunk-52-1.png" width="100%" /> ```r gd <- as.data.frame(gelman.diag(fit.commu, multivariate = FALSE)[[1]]) gd["Ntotal",] ``` ``` ## Point est. Upper C.I. ## Ntotal 1.001335 1.002693 ``` La prueba de Gelman-Rubin debe ser ~1 para considerar que hay buena convergencia. Aunque tenemos un valor bueno para Ntotal, hay varios valores de omega con NA, eso puede estar causando los problemas. --- class:inverse # Ajustamos nuestro primer modelo de multi-especie. Recordemos que **N** es la riqueza estimada a un **área mayor** de nuestro muestreo (área que no conocemos) .pull-left[ ### **N** depende de si cumplimos los supuestos del área - El muestreo es aleatorio. El área de muestreo debe **representar la región** - En caso contrario N representa el número de especies un área hipotética con **las mismas condiciones** - Si la región es pequeña, **N** puede ser **sobre-estimada** ] .pull-right[ ### y de las especies - Datos de **insectos** no sirven para predecir aves - Predicción a especies que sean **detectados de manera similar** con la metodología usada. .footnote[Guillera-Arroita, G, Kéry, M, Lahoz-Monfort, JJ. Inferring species richness using multispecies occupancy modeling: Estimation performance and interpretation. Ecol Evol. 2019; 9: 780– 792. https://doi.org/10.1002/ece3.4821] ] --- ## ¿Será mejor que un estimador no-paramétrico? ```r ## 5.1 Riqueza con Chao2---- # Formatear los datos a un vector de frecuencia inci_Chao <- ylist %>% # historias de captura map(~rowSums(.,na.rm = T)) %>% # sumo las detecciones en cada sitio reduce(cbind) %>% # unimos las listas t() %>% # trasponer la tabla as_tibble() %>% #formato tibble mutate_if(is.numeric,~(.>=1)*1) %>% #como es incidencia, formateo a 1 y 0 rowSums() %>% # ahora si la suma de las incidencias en cada sitio as_tibble() %>% add_row(value= 67, .before = 1) %>% # el formato requiere que el primer valor sea el número de sitios as.matrix() # Requiere formato de matriz # Calcular la riqueza con estimadores no paramétricos chao_sp <- ChaoSpecies(inci_Chao, datatype = "incidence_freq") NIChao <- chao_sp$Species_table[4,c(1,3,4)] # Extraer valores de IChao Nocu<- mcmc_intervals(fit.commu, pars = "Ntotal", prob = 0.95,prob_outer = 0.99, point_est = "mean")[[1]] %>% # Extraer valores del bayes plot select(m,l,h) %>% # Seleccionar columnas rename("Estimate"= m, # Renombrarlas "95%Lower"= l, "95%Upper"= h) ``` --- # Comparación Veamos de manera gráfica que tanto difieren las estimaciones de la riqueza .pull-left[ ```r # Unir en un solo dataframe Nplotdata <- rbind(IChao=NIChao, DR.mod=Nocu) %>% as.data.frame() %>% rownames_to_column(.) # Gráfico para comparar la riqueza estimada plotN <- ggplot(Nplotdata, aes(x=rowname, y= Estimate, col=rowname))+ geom_point(aes(shape=rowname),size=3)+ geom_errorbar(aes(ymin= `95%Lower`, ymax= `95%Upper`), width=.3, size=1)+ labs(x="Estimador de riqueza", y="Número de especies estimado", title = "Diferencia de los estimadores de riqueza")+ theme_classic()+ theme(text=element_text(size = 13), plot.title = element_text(hjust= 0.5), legend.position = "none") ``` ] .pull-right[ <img src="multi_sp_GAP_files/figure-html/unnamed-chunk-56-1.png" width="100%" /> ] --- class: inverse # Ocupación multi-specie vs IChao .pull-left[ ![](img/pildoras.jpg) .footnote[ Tingley, MW, Nadeau, CP, Sandor, ME. Multi-species occupancy models as robust estimators of community richness. Methods Ecol Evol. 2020; 11: 633– 642. https://doi.org/10.1111/2041-210X.13378] ] .pull-right[ Ambos lo hacen muy mal cuando `\(\psi_k\)` es muy **bajo** - **Chao** siempre estará sesgado al valor **más bajo de N** por no considerar la detección (De hecho Chao lo reconoce) - Los modelos de ocupación son muy **hambrientos** de datos - Tienes que saber "**programar**" para ajustar modelos multi-specie - Todavía hay algunos bmoles con la heterogeneidad y los prior (Guillera-Arroita et al. 2019) ] --- # Las posibilidades son infinitas 1. Estructura de la diversidad- Número efectivo de especies ![](img/dethill.png) --- ## 1. Estructura de la diversidad- Número efectivo de especies .pull-left[ Cómo la presencia o no de leones afecta la riqueza de meso-carnívoros >Curveira-Santos Gonçalo, Sutherland Chris, Tenan Simone, Fernández-Chacón lbert, Mann Gareth K. H., Pitman Ross T.and Swanepoel Lourens H. 2021. Mesocarnivore community structuring in the presence of Africa's apex predator. Proc. R. Soc. B.2882020237920202379. http://doi.org/10.1098/rspb.2020.2379 ] .pull-right[ ![](img/hilllion.png) ] --- # 2. Diversidad funcional y filogenética .pull-left[ Las consecuencias de no considerar a especies no detectadas en los análisis. > Jarzyna, M. A., & Jetz, W. (2016). Detecting the multiple facets of biodiversity. Trends in ecology & evolution, 31(7), 527-538. https://doi.org/10.1016/j.tree.2016.04.002 Les dejo ejemplos divertidos en la carpeta de bibliografía ] .pull-right[ ![](img/funcional.jpg) ] --- class: inverse, center, top background-image: url("img/peca.JPG") background-position: 60% 60% background-size: cover # .big-text[Gracias]