Aves carroñeras y rutas
Sergio Lambertucci y colaboradores pusieron cadáveres de ovejas a distintas distancias de rutas y observaron la actividad de aves carroñeras. Estudian la detección y el uso de carroña en relación a la distancia a las rutas por cóndores y otras aves rapaces. El trabajo completo está aquí
Estos datos presentan exceso de ceros (sitios donde no se registraron rapaces), entonces necesitamos escribir una función para la regresión binomial "inflada" en ceros. Para el total de aves observadas en cada sitio de carroña (volando y en tierra) podemos asumir una Poisson con exeso de ceros. Para la fracción de aves que bajan hasta la carroña asumimos una Binomial.
Además, vamos a usar la función plogis para pasar de valores reales (positivos y negativos) a valores entre cero y uno (este es un "logit link" en la jerga de GLM). Para ver como funciona esta función hacemos:
tmp <- seq(from = -5, to = 5, by = 0.01)
op <- par(cex.lab = 1.3, font.lab = 2, cex.axis = 1.3, bty = "n", las = 1, lwd = 3)
plot(tmp, plogis(tmp), type = "l")
par(op)De esta forma, podemos usar covariables para la probabilidad de éxito de una binomial. Todo junto queda:
library(bbmle)
condor <- read.table("https://sites.google.com/site/modelosydatos/condor.txt",
header = T)
head(condor) # revisamos los datos## obs Dist Condor CoFl Aguila AGFI Carancho CaFI Chimango ChiFI JoteN JNFI
## 1 1 50 1 57 0 11 58 54 14 33 12 28
## 2 2 90 0 45 1 46 7 12 0 3 16 31
## 3 3 100 0 9 0 0 20 4 3 1 0 4
## 4 4 120 3 19 0 0 30 0 14 4 0 0
## 5 5 150 0 170 0 24 4 5 2 2 0 1
## 6 6 200 0 0 0 2 32 4 10 7 0 0
## JoteC JCFI
## 1 Oto\xf1o Oto\xf1o
## 2 1 12
## 3 2 1
## 4 14 2
## 5 Oto\xf1o Oto\xf1o
## 6 0 1
zibinregc = function(p) {
av = p[1] # probabilidad (sin transformar) de detectar la carroña
ad = p[2] # ordenada al origen de la probabilidad de bajar a comer
bd = p[3] # la pendiente de la probabilidad de bajar a comer en función de la distancia a la ruta
an = p[4] # ordenada al origen de abundancia
bn = p[5] # pendiente de la abundancia en función de la distancia a la ruta
v = plogis(av) # probabilidad de detectar una carroña
pr = plogis(ad + bd * condor$Dist/1000) # probabilidad de bajar a comer
lambda = exp(an + bn * condor$Dist/1000) # Poisson en función de la distancia a la ruta
nllp = -sum(log(ifelse(k + f == 0, (1 - v) + v * dpois(0, lambda), v * dpois(k +
f, lambda))))
# Poisson likelihood inflada en ceros
nllb = -sum(dbinom(k, k + f, prob = pr, log = T)) # likelihood Binomial
return(nllp + nllb) # el resultado es la suma de las dos negative log-likelihoods
}
# definir los nombres de los parámetros para que los use mle2
parnames(zibinregc) = c("av", "ad", "bd", "an", "bn")
# Ajustamos la funcion para los cóndores
fit.C <- mle2(minuslogl = zibinregc, start = c(av = 1, ad = 2, bd = 0, an = 2,
bn = 0), data = list(k = condor$Condor, f = condor$CoFl, Dist = condor$Dist))
summary(fit.C) # ver los resultados## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = zibinregc, start = c(av = 1, ad = 2, bd = 0,
## an = 2, bn = 0), data = list(k = condor$Condor, f = condor$CoFl,
## Dist = condor$Dist))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## av 2.833206 1.028988 2.7534 0.005898 **
## ad -2.169019 0.132829 -16.3294 < 2.2e-16 ***
## bd 0.348259 0.065061 5.3528 8.661e-08 ***
## an 3.901983 0.043495 89.7107 < 2.2e-16 ***
## bn -0.138260 0.028856 -4.7914 1.656e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 1028.362
par_cof <- fit.C@coef # extraer los coeficientes (MLE de los parámetros)
d = 0:7000 # variable para graficar distancia
# gráfico de los ajustes
op <- par(mfrow = c(2, 1), cex.lab = 1.3, font.lab = 2, cex.axis = 1.3, bty = "n",
las = 1, lwd = 3)
plot(condor$Dist, (condor$Condor + condor$CoFl), ylab = "# de Cóndores", xlab = "")
lines(d, exp(par_cof[4] + par_cof[5] * d/1000))
plot(condor$Dist, condor$Condor/(condor$Condor + condor$CoFl), ylab = "Prop. Cóndores Comiendo",
xlab = "Distancia a la Ruta")
lines(d, plogis(par_cof[2] + par_cof[3] * d/1000))
par(op)
# Lo mismo para Jote negro
fit.J <- mle2(zibinregc, start = c(av = 2, ad = 2, bd = 0, an = 2, bn = 0),
data = list(k = condor$JoteN, f = condor$JNFI, Dist = condor$Dist))
summary(fit.J)## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = zibinregc, start = c(av = 2, ad = 2, bd = 0,
## an = 2, bn = 0), data = list(k = condor$JoteN, f = condor$JNFI,
## Dist = condor$Dist))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## av 0.013669 0.503786 0.0271 0.978354
## ad -0.643355 0.249231 -2.5814 0.009841 **
## bd -1.544973 0.853916 -1.8093 0.070408 .
## an 3.231078 0.118742 27.2108 < 2.2e-16 ***
## bn -1.628786 0.337199 -4.8303 1.363e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 185.9395
par_cof <- fit.J@coef
op <- par(mfrow = c(2, 1), cex.lab = 1.3, font.lab = 2, cex.axis = 1.3, bty = "n",
las = 1, lwd = 3)
plot(condor$Dist, (condor[, 11] + condor[, 12]), ylab = "# de Jote Negro", xlab = "")
lines(d, exp(par_cof[4] + par_cof[5] * d/1000))
plot(condor$Dist, condor[, 11]/(condor[, 11] + condor[, 12]), ylab = "Prop. Jote Negro Comiendo",
xlab = "Distancia a la Ruta")
lines(d, plogis(par_cof[2] + par_cof[3] * d/1000))
par(op)sessionInfo()## R version 3.2.2 (2015-08-14)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.10.5 (Yosemite)
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] bbmle_1.0.17 knitr_1.11
##
## loaded via a namespace (and not attached):
## [1] lattice_0.20-33 digest_0.6.8 MASS_7.3-44
## [4] mime_0.4 grid_3.2.2 knitrBootstrap_1.0.0
## [7] formatR_1.2 magrittr_1.5 evaluate_0.7.2
## [10] stringi_0.5-5 rmarkdown_0.8 tools_3.2.2
## [13] stringr_1.0.0 markdown_0.7.7 numDeriv_2014.2-1
## [16] yaml_2.1.13 htmltools_0.2.6