Establecimiento de Alisos en el Noroeste de Argentina
Volvemos al ejemplo de establecimiento de alisos en el NOA. Carguemos los datos y grafiquemos cómo fueron los niveles de lluvia en el pasado.
# Opción 1:
datos <- read.table("https://sites.google.com/site/modelosydatos/alisos.txt",
header = TRUE)
# Opción 2: datos <- read.table('alisos.txt', header = T)
yr <- datos$yr
ne <- datos$ne
lluvia <- datos$lluvia
cuenca <- datos$cuenca
op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(81 - yr[1:81], lluvia[1:81], type = "l", lwd = 2, xlab = "Año (hacia el pasado)",
ylab = "Lluvia (desvíos)")
par(op)Ahora veamos cómo cambia el número de establecimientos de Aliso en las distintas cuencas en función de la variabilidad en la precipitación.
ncuencas <- length(unique(cuenca))
op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(lluvia, ne, xlab = "lluvia (desvíos)", ylab = "Número de Establecimientos")
for (i in 1:ncuencas) {
tmp <- lm(ne[cuenca == i] ~ lluvia[cuenca == i]) #No pooling
points(lluvia[cuenca == i], ne[cuenca == i], pch = 21, bg = i + 1)
abline(tmp, col = i + 1, lwd = 3)
}
tmp <- lm(ne ~ lluvia) #Complete pooling
abline(tmp, lwd = 4)
par(op)Los colores de puntos y líneas corresponden a las distintas cuencas (datos y regresiones lineales simples de no-pooling) y la línea negra gruesa corresponde a una regresión simple para todos los datos (complete pooling). Vamos a ajustar un modelo jerárquico teniendo en cuenta la variabilidad entre cuencas.
A diferencia del ejemplo de las frambuesas, vamos a modelar tanto la ordenada al origen como la pendiente ya que cada cuenca tiene su propio gradiente de precipitaciones.
Análisis jerárquico
Escribimos el modelo:
cat(file = "aliso.bug", "
model{
# modelo de datos
for(i in 1:nobs){
ne[i]~dpois(lambda[i]) # establecimientos por año
lambda[i] <- exp( a[cuenca[i]] + b[cuenca[i]]*lluvia[i] )
}
# previas a nivel de cuencas
for(j in 1:ncuencas){
a[j]~dnorm(mua,taua)
b[j]~dnorm(mub,taub)
}
# previas a nivel poblacional
mua~dnorm(0,0.001) \t
sa~dunif(0,100)\t\t
taua<-1/(sa*sa) \t
mub~dnorm(0,0.001)\t
sb~dunif(0,100)
taub<-1/(sb*sb)
}
")Luego definimos los datos, la función para los valores iniciales y los parámetros que queremos guardar.
nobs = length(yr)
aliso.data <- list("ne", "lluvia", "cuenca", "nobs", "ncuencas")
inits <- function() {
list(a = rnorm(ncuencas, 0, 0.5), b = rnorm(ncuencas, 0, 0.5), mua = rnorm(1,
0, 0.5), mub = rnorm(1, 0, 0.5), sa = runif(1), sb = runif(1))
}
parameters <- c("a", "b", "mua", "mub", "sa", "sb")Ahora llamamos a JAGS
library(jagsUI)
aliso.sim <- jags(data = aliso.data, inits, parameters, model.file = "aliso.bug",
n.thin = 1, n.chains = 3, n.iter = 5000)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 567
## Unobserved stochastic nodes: 18
## Total graph size: 3438
##
## Initializing model
##
## Adaptive phase.....
## Adaptive phase complete
##
## No burn-in specified
##
## Sampling from joint posterior, 5000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.Vemos que hayan convergido y que los tamaños efectivos de las muestras de las posteriores (n.eff) sean razonables.
print(aliso.sim)## JAGS output for model 'aliso.bug', generated by jagsUI.
## Estimates based on 3 chains of 5000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 0 iterations and thin rate = 1,
## yielding 15000 total samples from the joint posterior.
## MCMC ran for 0.278 minutes at time 2019-03-05 22:06:40.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat
## a[1] 1.497 0.047 1.404 1.497 1.589 FALSE 1.000 1.001
## a[2] 1.121 0.063 0.993 1.121 1.242 FALSE 1.000 1.001
## a[3] 0.912 0.072 0.768 0.913 1.049 FALSE 1.000 1.000
## a[4] 1.654 0.048 1.559 1.654 1.745 FALSE 1.000 1.001
## a[5] 1.009 0.067 0.875 1.009 1.136 FALSE 1.000 1.001
## a[6] 1.155 0.060 1.036 1.155 1.272 FALSE 1.000 1.001
## a[7] 2.162 0.036 2.090 2.162 2.233 FALSE 1.000 1.001
## b[1] -0.037 0.046 -0.127 -0.037 0.054 TRUE 0.785 1.001
## b[2] -0.146 0.066 -0.277 -0.146 -0.020 FALSE 0.990 1.002
## b[3] 0.365 0.066 0.238 0.365 0.496 FALSE 1.000 1.001
## b[4] 0.124 0.047 0.031 0.124 0.218 FALSE 0.996 1.000
## b[5] 0.119 0.063 -0.006 0.119 0.244 TRUE 0.970 1.001
## b[6] 0.155 0.060 0.037 0.156 0.274 FALSE 0.996 1.000
## b[7] 0.218 0.036 0.147 0.218 0.287 FALSE 1.000 1.003
## mua 1.358 0.235 0.882 1.359 1.821 FALSE 1.000 1.001
## mub 0.115 0.097 -0.077 0.115 0.304 TRUE 0.907 1.002
## sa 0.577 0.234 0.299 0.522 1.168 FALSE 1.000 1.003
## sb 0.226 0.104 0.107 0.202 0.476 FALSE 1.000 1.008
## deviance 3625.686 9.412 3609.659 3624.973 3645.836 FALSE 1.000 1.000
## n.eff
## a[1] 3385
## a[2] 3411
## a[3] 7305
## a[4] 15000
## a[5] 2549
## a[6] 6806
## a[7] 3302
## b[1] 3533
## b[2] 962
## b[3] 1933
## b[4] 6196
## b[5] 2398
## b[6] 14210
## b[7] 761
## mua 15000
## mub 5149
## sa 1197
## sb 1360
## deviance 7694
##
## Successful convergence based on Rhat values (all < 1.1).
## Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## For each parameter, n.eff is a crude measure of effective sample size.
##
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
##
## DIC info: (pD = var(deviance)/2)
## pD = 44.3 and DIC = 3669.977
## DIC is an estimate of expected predictive error (lower is better).Podemos ver que para algunas cuencas como la \(3\), \(4\), \(6\) y \(7\) los establecimientos aumentan con la precipitación (b \(> 0\)), mientras que en otras cuencas como la \(2\) el efecto es negativo. En general, la tendencia es a un aumento de los establecimientos con las precipitaciones (alrededor del \(90\)% de la posterior de mub es mayor que cero). Podemos ver esto (más o menos) gráficamente:
r <- seq(min(lluvia), max(lluvia), 0.001)
op <- par(mar = c(4, 4, 1, 1) + 0.1, las = 1, cex.lab = 1.5, cex.axis = 1.3)
plot(lluvia, ne, xlab = "lluvia (desvíos)", ylab = "Número de Establecimientos")
for (i in 1:ncuencas) {
points(lluvia[cuenca == i], ne[cuenca == i], pch = 16, col = i + 1)
lines(r, exp(aliso.sim$mean$a[i] + aliso.sim$mean$b[i] * r), col = i + 1,
lwd = 3)
}
abline(aliso.sim$mean$mua, aliso.sim$mean$mub, lwd = 4)
par(op)Pero es posible que no observemos todos los establecimientos. Un grupo de árboles puede desaparecer antes de que los encontremos. Entonces, modelamos como variable oculta al número de establecimientos por año y modelamos la probabilidad de observar un establecimeinto en función del tiempo p[i] <- exp(-d * yr[i]).
Este es el primero de varios modelos que vamos a ver que tienen variables ocultas. A veces es interesante tener en cuenta que no observamos completamente todos los procesos que ocurren en la naturaleza, como en este caso ocurre con los árboles que desaparecen.
Miremos el modelo de abajo y prestemos atención línea por línea. N es una variable oculta que representa el número de arboles que se establecen en la cuenca. Como dijimos arriba esos establecimientos los observamos de manera imperfecta por eso es una variable oculta. Ese N se distribuye con una Poisson como ocurría en el modelo anterior. Ahora bien, de donde sale N? Ahora los establecimiento observados los modelamos con una distribución Binomial asumiendo que nosotros observamos una dada cantidad de arboles establecidos de un N total. Finalmente la probabilidad de observar un establecimiento (p) decae con el tiempo.
cat(file = "alisop.bug", "
model{
for(i in 1:(nobs)){
N[i] ~ dpois(lambda[i]) # variable oculta
lambda[i] <- exp( a[cuenca[i]] + b[cuenca[i]]*lluvia[i])
ne[i] ~ dbin(p[i],N[i]) # establecimientos observados
p[i] <- exp(-d * yr[i]) # probabilidad de observar un establecimiento en función del tiempo
}
#Previas a nivel de cuencas\t
for(j in 1:ncuencas){
a[j]~dnorm(mua,taua)
b[j]~dnorm(mub,taub)
}
# previas para todos los parámetros\t
mua~dnorm(0,0.001)
sa~dunif(0,100)
taua<-1/(sa*sa)
mub~dnorm(0,0.001)
sb~dunif(0,100)
taub<-1/(sb*sb)
d~dunif(0,100)
}
")Luego definimos datos, initis y parámetros para poder ajustar el modelo
alisop.data <- list("yr", "ne", "lluvia", "cuenca", "nobs", "ncuencas")
initsp <- function() {
list(a = rnorm(ncuencas, 6, 0.5), b = rnorm(ncuencas, 1, 0.1), mua = rnorm(1,
0, 0.5), mub = rnorm(1, 0, 0.5), sa = runif(1), sb = runif(1), d = 1e-05)
}
parametersp <- c("a", "b", "mua", "mub", "sa", "sb", "d")
alisop.sim <- jags(data = alisop.data, initsp, parametersp, model.file = "alisop.bug",
n.thin = 5, n.chains = 3, n.iter = 5000)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 567
## Unobserved stochastic nodes: 586
## Total graph size: 4738
##
## Initializing model
##
## Adaptive phase.....
## Adaptive phase complete
##
## No burn-in specified
##
## Sampling from joint posterior, 5000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.print(alisop.sim)## JAGS output for model 'alisop.bug', generated by jagsUI.
## Estimates based on 3 chains of 5000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 0 iterations and thin rate = 5,
## yielding 3000 total samples from the joint posterior.
## MCMC ran for 0.569 minutes at time 2019-03-05 22:06:57.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat
## a[1] 2.216 0.058 2.103 2.216 2.330 FALSE 1.000 1.002
## a[2] 1.979 0.070 1.843 1.978 2.116 FALSE 1.000 1.014
## a[3] 1.766 0.079 1.612 1.764 1.919 FALSE 1.000 1.006
## a[4] 2.507 0.057 2.393 2.506 2.612 FALSE 1.000 1.004
## a[5] 1.858 0.072 1.717 1.857 2.001 FALSE 1.000 1.009
## a[6] 2.002 0.069 1.861 2.004 2.131 FALSE 1.000 1.013
## a[7] 2.934 0.049 2.836 2.935 3.029 FALSE 1.000 1.004
## b[1] -0.185 0.050 -0.281 -0.185 -0.087 FALSE 1.000 1.000
## b[2] -0.362 0.065 -0.488 -0.361 -0.237 FALSE 1.000 1.023
## b[3] 0.119 0.066 -0.008 0.118 0.249 TRUE 0.965 1.000
## b[4] -0.098 0.047 -0.190 -0.100 -0.005 FALSE 0.979 1.003
## b[5] -0.112 0.065 -0.241 -0.112 0.016 TRUE 0.959 1.007
## b[6] -0.069 0.058 -0.185 -0.070 0.041 TRUE 0.877 1.007
## b[7] 0.024 0.037 -0.049 0.025 0.097 TRUE 0.750 1.002
## mua 2.179 0.242 1.717 2.180 2.628 FALSE 1.000 1.003
## mub -0.098 0.088 -0.270 -0.097 0.070 TRUE 0.902 1.001
## sa 0.543 0.235 0.275 0.483 1.142 FALSE 1.000 1.010
## sb 0.207 0.095 0.098 0.187 0.424 FALSE 1.000 1.016
## d 0.024 0.001 0.022 0.024 0.026 FALSE 1.000 1.006
## deviance 2132.543 44.602 2046.146 2132.979 2218.650 FALSE 1.000 1.006
## n.eff
## a[1] 727
## a[2] 147
## a[3] 332
## a[4] 423
## a[5] 224
## a[6] 155
## a[7] 699
## b[1] 2575
## b[2] 106
## b[3] 3000
## b[4] 1065
## b[5] 270
## b[6] 262
## b[7] 1216
## mua 3000
## mub 2906
## sa 1043
## sb 797
## d 396
## deviance 320
##
## Successful convergence based on Rhat values (all < 1.1).
## Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## For each parameter, n.eff is a crude measure of effective sample size.
##
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
##
## DIC info: (pD = var(deviance)/2)
## pD = 989.1 and DIC = 3121.618
## DIC is an estimate of expected predictive error (lower is better).Luego de chequear convergencia y “n.eff” vemos que ahora la relación de los establecimientos con la lluvia es mayormente negativa!
Podemos ver como cambia la probabilidad de encontrar un establecimiento con el tiempo
t = 0:80
plot(t, exp(-alisop.sim$mean$d * t), lwd = 2, type = "l", xlab = "tiempo (años)",
ylab = "Pr Supervivencia al Presente")
Y el establecimiento en función de la lluvia
r = seq(min(lluvia), max(lluvia), 0.001)
plot(lluvia, ne, xlab = "lluvia (desvíos)", ylab = "Número de Establecimientos")
for (i in 1:ncuencas) {
points(lluvia[cuenca == i], ne[cuenca == i], col = i)
lines(r, exp(alisop.sim$mean$a[i] + alisop.sim$mean$b[i] * r), col = i,
lwd = 2)
}
abline(alisop.sim$mean$mua, alisop.sim$mean$mub, lwd = 4)
par(op)Ejercicio:
Ajustar este modelo con variable oculta pero ahora asumiendo que la relación entre el número de árboles que se establecen y la variabilidad en la precipitación tiene un óptimo.
Juan Manuel Morales. 30 de Septiembre del 2015. Última actualización: 2019-03-05