Capítulo 7 Análisis de ocupación, metodo ML
Ya que hemos entendido como funcionan e interactúan los dos procesos; el ecológico y el observacional para producir los datos de ocupación. Luego de generar varios sets de datos, ahora solo nos resta analizarlos. La forma más directa e intuitiva es usar la función occu del paquete unmarked (Fiske & Chandler, 2011). Posteriormente podremos usar un modelo de tipo bayesiano en el lenguaje BUGS para analizar los mismos datos y al final comparar cuál de los dos estimadores, Máxima Verosimilitud o Bayesiano, se acerca más a los parámetros verdaderos.
7.0.1 Limpiando la memoria de R
Antes de continuar, y como ya hemos generado una gran cantidad de datos y modelos en los pasos previos, vamos a borrar la memoria de R antes de comenzar. Esto lo hacemos con el comando:
Una vez hayamos hecho esto debemos volver a correr el código de la función que genera los datos que hemos creado el Capítulo 6.
7.1 Generando los datos
Esta vez recurriremos a un diseño tipo TEAM (http://www.teamnetwork.org) con 60 sitios de muestreo y 30 visitas repetidas, que equivalen a los 30 días en que las cámaras permanecen activas en campo. De nuevo nuestra especie es el venado de cola blanca. Para este ejemplo asumiremos que la detección es 0.6, la ocupación 0.8 y las interacciones son mucho más sencillas con la altitud como la única covariable que explica la ocupación. Sin embargo para la detección hay una relación más compleja, asumiendo que hay una leve interacción entre las covariables de la observación. Para la observación la altitud y temperatura interactúan entre sí. También observe como la altitud influye en direcciones opuestas con un signo positivo en la altitud para la detección y negativo para la ocupación.
# Data generation
# Lets build a model were elevation explain occupancy and p has interactions
datos2<-data.fn(M = 60, J = 30, show.plot = FALSE,
mean.occupancy = 0.8, beta1 = -1.5, beta2 = 0, beta3 = 0,
mean.detection = 0.6, alpha1 = 2, alpha2 = 1, alpha3 = 1.5
)
# Function to simulate occupancy measurements replicated at M sites during J occasions.
# Population closure is assumed for each site.
# Expected occurrence may be affected by elevation (elev),
# forest cover (forest) and their interaction.
# Expected detection probability may be affected by elevation,
# temperature (temp) and their interaction.
# Function arguments:
# M: Number of spatial replicates (sites)
# J: Number of temporal replicates (occasions)
# mean.occupancy: Mean occurrence at value 0 of occurrence covariates
# beta1: Main effect of elevation on occurrence
# beta2: Main effect of forest cover on occurrence
# beta3: Interaction effect on occurrence of elevation and forest cover
# mean.detection: Mean detection prob. at value 0 of detection covariates
# alpha1: Main effect of elevation on detection probability
# alpha2: Main effect of temperature on detection probability
# alpha3: Interaction effect on detection of elevation and temperature
# show.plot: if TRUE, plots of the data will be displayed;
# set to FALSE if you are running simulations.
#To make the objects inside the list directly accessible to R, without having to address
#them as data$C for instance, you can attach datos2 to the search path.
attach(datos2) # Make objects inside of 'datos2' accessible directly
#Remember to detach the data after use, and in particular before attaching a new data
#object, because more than one data set attached in the search path will cause confusion.
# detach(datos2) # Make clean up7.2 Poniendo los datos en unmarked
Unmarked es el paquete de R que usamos para analizar los datos de ocupación (Fiske & Chandler, 2011). Para lograr esto debemos primero preparar los datos y juntarlos en un objeto de tipo unmarkedFrame. En este caso usamos la función unmarkedFrameOccu que es específica para análisis de ocupación de una sola época o estación. Más sobre unkarked en: https://sites.google.com/site/unmarkedinfo/home
library(unmarked)
siteCovs <- as.data.frame(cbind(forest,elev)) # occupancy covariates
obselev<-matrix(rep(elev,J),ncol = J) # make elevetion per observation
obsCovs <- list(temp= temp,elev=obselev) # detection covariates
# umf is the object joining observations (y), occupancy covariates (siteCovs)
# and detection covariates (obsCovs). Note that obsCovs should be a list of
# matrices or dataframes.
umf <- unmarkedFrameOccu(y = y, siteCovs = siteCovs, obsCovs = obsCovs)El paquete unmarked nos permite ver gráficamente como están dispuestos los datos en los sitios de muestreo con la función plot.
Figure 7.1: Inspección grafica del objeto umf.
7.3 Ajustando los modelos
El siguiente paso es ajustar los modelos que se requerían variando las co-variables. Esto se logra con la función occu().
Tenga en cuenta que en el proceso de construcción de modelos su modelo debe tener un significado biológico.
# detection first, occupancy next
fm0 <- occu(~1 ~1, umf) # Null model
fm1 <- occu(~ elev ~ 1, umf) # elev explaining detection
fm2 <- occu(~ elev ~ elev, umf) # elev explaining detection and occupancy
fm3 <- occu(~ temp ~ elev, umf)
fm4 <- occu(~ temp ~ forest, umf)
fm5 <- occu(~ elev + temp ~ 1, umf)
fm6 <- occu(~ elev + temp + elev:temp ~ 1, umf)
fm7 <- occu(~ elev + temp + elev:temp ~ elev, umf)
fm8 <- occu(~ elev + temp + elev:temp ~ forest, umf)7.4 Selección de modelos
Unmarked permite hacer selección de modelos basándose en el AIC de cada uno. De forma tal que el menor AIC es el modelo más parsimonioso de acuerdo a nuestros datos (Burnham & Anderson, 2004).
models <- fitList( # here e put names to the models
'p(.)psi(.)' = fm0,
'p(elev)psi(.)' = fm1,
'p(elev)psi(elev)' = fm2,
'p(temp)psi(elev)' = fm3,
'p(temp)psi(forest)' = fm4,
'p(temp+elev)psi(.)' = fm5,
'p(temp+elev+elev:temp)psi(.)' = fm6,
'p(temp+elev+elev:temp)psi(elev)' = fm7,
'p(temp+elev+elev:temp)psi(forest)' = fm8)
modSel(models) # model selection procedure## nPars AIC delta
## p(temp+elev+elev:temp)psi(elev) 6 1630.77 0.00
## p(temp+elev+elev:temp)psi(.) 5 1636.27 5.49
## p(temp+elev+elev:temp)psi(forest) 6 1638.24 7.47
## p(temp+elev)psi(.) 4 1672.55 41.78
## p(elev)psi(elev) 4 1741.34 110.57
## p(elev)psi(.) 3 1746.83 116.06
## p(temp)psi(elev) 4 1898.56 267.79
## p(temp)psi(forest) 4 1906.03 275.26
## p(.)psi(.) 2 1969.25 338.48
## AICwt cumltvWt
## p(temp+elev+elev:temp)psi(elev) 9.2e-01 0.92
## p(temp+elev+elev:temp)psi(.) 5.9e-02 0.98
## p(temp+elev+elev:temp)psi(forest) 2.2e-02 1.00
## p(temp+elev)psi(.) 7.8e-10 1.00
## p(elev)psi(elev) 9.0e-25 1.00
## p(elev)psi(.) 5.8e-26 1.00
## p(temp)psi(elev) 6.5e-59 1.00
## p(temp)psi(forest) 1.6e-60 1.00
## p(.)psi(.) 2.9e-74 1.007.5 Predicción en graficas
El modelo con menor AIC puede ser usado para predecir resultados esperados de acuerdo a un nuevo set de datos. Por ejemplo, uno podría preguntar la abundancia de venados que se espera encontrar en un sitio con mayor altitud. Las predicciones también son otra forma de presentar los resultados de un análisis. Aquí ilustraremos como se ve la predicción de \(\psi\) y p sobre el rango de las covariables estudiadas. Note que estamos usando covariables estandarizadas. Si estuviéramos usando covariables en su escala real, tendríamos que tener en cuenta que hay que transformarlas usando la media y la desviación estándar.
Antes de usar el modelo para predecir es buena idea verificar los parámetros del modelo y sus errores, posteriormente verificar gráficamente que el modelo ajusta bien con la función parboot, la cual hace un remuestreo del modelo. Esta grafica se interpreta como que el modelo tiene buen ajuste, cuando la media (línea punteada) está entre los intervalos del histograma. Si la línea se encuentra muy lejos del histograma el modelo no podría ser bueno para predecir.
##
## Call:
## occu(formula = ~elev + temp + elev:temp ~ elev, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 1.41 0.366 3.85 0.000117
## elev -1.63 0.656 -2.48 0.013082
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.53 0.0698 7.60 3.04e-14
## elev 1.80 0.1277 14.10 3.56e-45
## temp 1.16 0.1209 9.57 1.08e-21
## elev:temp 1.33 0.2211 6.00 1.96e-09
##
## AIC: 1630.771
## Number of sites: 60
## optim convergence code: 0
## optim iterations: 31
## Bootstrap iterations: 0## t0 = 409.3793
Figure 7.2: Evaluación grafica del ajuste del modelo fm7.
Ahora que sabemos que nuestro mejor modelo tiene buen ajuste, podemos usarlo para predecir la ocupación en el rango de la altitud para ver su comportamiento en una gráfica.
elevrange<-data.frame(elev=seq(min(datos2$elev),max(datos2$elev),length=100)) # newdata
pred_psi <-predict(fm7,type="state",newdata=elevrange,appendData=TRUE)
plot(Predicted~elev, pred_psi,type="l",col="blue",
xlab="elev",
ylab="psi")
lines(lower~elev, pred_psi,type="l",col=gray(0.5))
lines(upper~elev, pred_psi,type="l",col=gray(0.5))
Figure 7.3: Grafica de la ocupación con respecto a la elevación.
7.6 Predicción espacialmente explicita
Podemos también usar el mejor modelo para predecir de forma espacialmente explicita sí tenemos los mapas. Como ilustración vamos a construir mapas simulados para cada una de nuestras covariables. Los mapas surgen de un patrón aleatorio de puntos con distribución Poisson. Luego estos puntos los convertimos en una superficie interpolada.
# lets make random maps for the three covariates
library(raster)
library(spatstat)
set.seed(24) # Remove for random simulations
# CONSTRUCT ANALYSIS WINDOW USING THE FOLLOWING:
xrange=c(-2.5, 1002.5)
yrange=c(-2.5, 502.5)
window<-owin(xrange, yrange)
# Build maps from random points and interpole in same line
elev <- density(rpoispp(lambda=0.6, win=window)) #
forest <- density(rpoispp(lambda=0.2, win=window)) #
temp <- density(rpoispp(lambda=0.5, win=window)) #
# Convert covs to raster and Put in the same stack
mapdata.m<-stack(raster(elev),raster(forest), raster(temp))
names(mapdata.m)<- c("elev", "forest", "temp") # put names to raster
# lets plot the covs maps
plot(mapdata.m)
Figure 7.4: Mapa simulado de elevación, bosque y temperatura.
Una vez tenemos nuestros mapas de covariables, los usamos para predecir con el mejor modelo. De esta forma podemos tener un mapa con predicciones de la ocupación y la probabilidad de detección.
# make the predictions
predictions_psi <- predict(fm7, type="state", newdata=mapdata.m) # predict psi## doing row 1000 of 16384
## doing row 2000 of 16384
## doing row 3000 of 16384
## doing row 4000 of 16384
## doing row 5000 of 16384
## doing row 6000 of 16384
## doing row 7000 of 16384
## doing row 8000 of 16384
## doing row 9000 of 16384
## doing row 10000 of 16384
## doing row 11000 of 16384
## doing row 12000 of 16384
## doing row 13000 of 16384
## doing row 14000 of 16384
## doing row 15000 of 16384
## doing row 16000 of 16384## doing row 1000 of 16384
## doing row 2000 of 16384
## doing row 3000 of 16384
## doing row 4000 of 16384
## doing row 5000 of 16384
## doing row 6000 of 16384
## doing row 7000 of 16384
## doing row 8000 of 16384
## doing row 9000 of 16384
## doing row 10000 of 16384
## doing row 11000 of 16384
## doing row 12000 of 16384
## doing row 13000 of 16384
## doing row 14000 of 16384
## doing row 15000 of 16384
## doing row 16000 of 16384# put in the same stack
predmaps<-stack(predictions_psi$Predicted,predictions_p$Predicted)
names(predmaps)<-c("psi_predicted", "p_predicted") # put names
plot(predmaps)
Figure 7.5: Modelos de ocupación y detección espacialmente explícitos.
References
Burnham, K. P., & Anderson, D. R. (2004). Information and likelihood theory: A basis for model selection and inference. In K. P. Burnham & D. R. Anderson (Eds.), Model selection and multimodel inference: A practical information-theoretic approach (pp. 49–97). Retrieved from http://link.springer.com/chapter/10.1007%2F978-0-387-22456-5_2
Fiske, I., & Chandler, R. (2011). unmarked : An R Package for fitting hierarchical models of wildlife occurrence and abundance. Journal of Statistical Software, 43(10), 1–23. doi: 10.18637/jss.v043.i10