Daño floral en función de frecuencia de visitas

La frambuesa es un cultivo muy importante en la Comarca Andina y depende fuertemente de polinizadores. Para que un fruto de frambuesa quede bien formado es necesario que sea polinizado y en general se supone que la presencia de polinizadores en el cultivo es beneficiosa. Sin embargo, quizás demasiados polinizadores resulten perjudiciales. Un exceso de visitas puede dañar los pistilos de las flores y producir frutos de menor calidad. ¿Pero cuánto es un exceso de visitas? Con esta pregunta en mente, Agustín Saez estimó la tasa de vistas en \(16\) cultivos de la Comarca Andina y además colectó flores de distintas plantas en cada cultivo para ver cuántos pistilos sanos y rotos tenían. Veamos los datos y grafiquemos cómo cambia la proporción de pistilos rotos en función de la frecuencia de visitas.

datos <- read.table("https://sites.google.com/site/modelosydatos/Datosframbuesas.txt", 
    header = TRUE)

Empecemos con un gráfico:

Nsitio <- as.numeric(datos$Sitio)

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")

plot(1, type = "n", ylim = c(0, 1.2), xlim = c(0, 3.5), pch = 21, xlab = "Frecuencia de Visitas", 
    ylab = "Prop. de pistilos rotos")

colores <- rainbow(16, alpha = 0.5)
for (i in 1:max(Nsitio)) {
    idx <- which(Nsitio == i)
    points((datos$Frecuencia[idx]), jitter(datos$Rotos[idx]/datos$Total[idx]), 
        pch = 21, bg = colores[i], cex = 1.5)
}

par(op)

Los colores de los puntos corresponden a los distintos sitios. Como se puede ver, hay una sola estimación de frecuencia de visitas por sitio. Además de variar por la frecuencia de vistas, cada sitio puede tener una proporción de pistilos rotos por otras causas como por ejemplo heladas. Vamos a ajustar un modelo de daño de pistilos en función de la frecuencia de visitas por flor donde permitamos que cada sitio tenga valores diferentes de ordenada al origen. Es decir, el daño esperado cuando tenemos cero visitas. Vamos a asumir que la respuesta en cuanto a la frecuencia de visitas es la misma para todos los sitios.

Ejercicio

• Escriban (en “papel y lápiz”) el modelo que acabamos de describir.

---

Ahora vamos a ajustar usando JAGS el modelo jerárquico que acabamos de describir. Para eso, escribimos el modelo en lenguaje BUGS:

cat(file = "mod1.bug", "
model {
  for(i in 1:n) {
    rotos[i] ~ dbin(p[i], N[i])
    logit(p[i]) <- b0[Nsitio[i]] + b1 * freq[i]
    }
    
  for( j in 1:S){
    b0[j] ~ dnorm(mub0, taub0)
  }

  b1 ~ dnorm(0, 1)
  mub0  ~ dnorm(0, 0.01)
  taub0 <- 1/(sdb0 * sdb0)
  sdb0 ~ dnorm(0, 1)T(0, )
  }
")

Luego definimos los datos, la función para los valores iniciales y los parámetros que queremos guardar.

rotos <- datos$Rotos
freq <- datos$Frecuencia
n <- length(rotos)
N <- datos$Total
Nsitio <- as.numeric(datos$Sitio)
S <- max(Nsitio)

data <- list("rotos", "freq", "N", "n", "S", "Nsitio")
params <- c("b0", "b1", "mub0", "sdb0")
inits <- function() {
    list(b0 = rnorm(S, 0, 1), b1 = rnorm(1, 0, 1), mub0 = rnorm(1, 0, 1), sdb0 = runif(1, 
        0, 10))
}

Ahora definimos cuántas iteraciones queremos correr, el “thining”, el “burn in”, número de cadenas:

ni <- 5000  # número de iteraciones
nt <- 1  # si descartamos 'thin' algunos valores
nb <- 2500  # cuantas iteraciones usamos de 'burn in'
nc <- 3  # y cuantas cadenas corremos

Ahora llamamos a JAGS

library(jagsUI)
m1.sim <- jags(data, inits, params, model.file = "mod1.bug", n.chains = nc, 
    n.iter = ni, n.burnin = nb, n.thin = nt)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 477
##    Unobserved stochastic nodes: 19
##    Total graph size: 1987
## 
## Initializing model
## 
## Adaptive phase..... 
## Adaptive phase complete 
##  
## 
##  Burn-in phase, 2500 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 2500 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(m1.sim)

## JAGS output for model 'mod1.bug', generated by jagsUI.
## Estimates based on 3 chains of 5000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 2500 iterations and thin rate = 1,
## yielding 7500 total samples from the joint posterior. 
## MCMC ran for 0.121 minutes at time 2019-03-02 13:15:15.
## 
##              mean    sd     2.5%      50%    97.5% overlap0     f  Rhat
## b0[1]      -3.594 0.575   -4.702   -3.596   -2.471    FALSE 1.000 1.002
## b0[2]      -2.833 0.342   -3.503   -2.839   -2.154    FALSE 1.000 1.004
## b0[3]      -4.339 0.602   -5.614   -4.292   -3.327    FALSE 1.000 1.056
## b0[4]      -3.784 0.452   -4.665   -3.794   -2.872    FALSE 1.000 1.004
## b0[5]      -3.972 0.685   -5.317   -3.974   -2.610    FALSE 1.000 1.000
## b0[6]      -4.509 0.408   -5.345   -4.489   -3.756    FALSE 1.000 1.061
## b0[7]      -4.637 0.435   -5.478   -4.634   -3.785    FALSE 1.000 1.011
## b0[8]      -2.986 0.281   -3.559   -2.980   -2.442    FALSE 1.000 1.051
## b0[9]      -2.912 0.348   -3.605   -2.917   -2.248    FALSE 1.000 1.022
## b0[10]     -4.480 0.424   -5.329   -4.478   -3.658    FALSE 1.000 1.029
## b0[11]     -1.549 0.787   -3.024   -1.575    0.024     TRUE 0.973 1.006
## b0[12]     -5.214 0.728   -6.646   -5.218   -3.779    FALSE 1.000 1.003
## b0[13]     -2.665 0.295   -3.295   -2.641   -2.145    FALSE 1.000 1.046
## b0[14]     -3.958 0.474   -4.899   -3.943   -3.064    FALSE 1.000 1.022
## b0[15]     -2.546 0.732   -3.957   -2.551   -1.095    FALSE 0.999 1.001
## b0[16]     -3.069 0.803   -4.641   -3.083   -1.449    FALSE 1.000 1.001
## b1          1.490 0.252    0.988    1.491    1.979    FALSE 1.000 1.001
## mub0       -3.562 0.480   -4.512   -3.565   -2.599    FALSE 1.000 1.004
## sdb0        1.079 0.222    0.734    1.050    1.593    FALSE 1.000 1.011
## deviance 1276.447 5.359 1267.708 1275.836 1288.291    FALSE 1.000 1.005
##          n.eff
## b0[1]     1570
## b0[2]     1567
## b0[3]       43
## b0[4]      691
## b0[5]     7500
## b0[6]       41
## b0[7]      212
## b0[8]       46
## b0[9]      160
## b0[10]     109
## b0[11]     354
## b0[12]     789
## b0[13]      52
## b0[14]     188
## b0[15]    1285
## b0[16]    2270
## b1        1862
## mub0       604
## sdb0       191
## deviance   647
## 
## 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 = 14.3 and DIC = 1290.765 
## DIC is an estimate of expected predictive error (lower is better).

Podemos ver como es la posterior de la pendiente

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(density(m1.sim$sims.list$b1), lwd = 2, main = "", xlab = "b1")

par(op)

Parece que hay un efecto claro de la frecuencia de visitas en la proporción de pistilos rotos. También podemos ver gráficamente las diferencias entre las posteriores de las ordenadas al origen para los distintos sitios b0 y la ordenada al origen para la “población” de sitios mub0.

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")

plot(density(m1.sim$sims.list$mub0), ylim = c(0, 2), xlim = c(-10, 3), lwd = 4, 
    main = "", xlab = "Ordenada al origen")

for (i in 1:S) {
    lines(density(m1.sim$sims.list$b0[, i]), col = "gray", lwd = 2)
}

par(op)

Veamos ahora las curvas estimadas para cada sitio en un gráfico

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")

plot(1, type = "n", ylim = c(0, 1.2), xlim = c(0, 3.5), pch = 20, xlab = "Frecuencia de Visitas", 
    ylab = "Proporción de pistilos rotos")

colores <- rainbow(16, alpha = 0.5)

xf <- seq(min(freq), max(freq), length.out = 100)

for (i in 1:max(Nsitio)) {
    idx <- which(Nsitio == i)
    points((datos$Frecuencia[idx]), jitter(datos$Rotos[idx]/datos$Total[idx]), 
        pch = 21, bg = colores[i])
    
    lines(xf, plogis(m1.sim$mean$b0[i] + m1.sim$mean$b1 * xf), col = colores[i], 
        lwd = 2)
    
}
lines(xf, plogis(m1.sim$mean$mub0 + m1.sim$mean$b1 * xf), lwd = 4)

par(op)

Pero también teníamos plantas dentro de sitios. Podemos tener cuanta la variablidad entre plantas asumiendo que no todas las plantas dentro de un sitio están en la mismas condiciones.

Escribimos el modelo:

cat(file = "mod2.bug", "    
model {
  for(i in 1: n) {
    rotos[i] ~ dbin(p[i], N[i])
    logit(p[i]) <- b0[Nplanta[i]] + b1 * freq[i]
    }
    
  for( j in 1:P){
    b0[j] ~ dnorm(mub0[Ns[j]], taub0)
  }
    
  for(s in 1:S){
    mub0[s] ~ dnorm(mub0s, taub0s)
  }
    
  b1 ~ dnorm(0, 0.1)
    
  taub0 <- 1/(sdb0 * sdb0)
  sdb0 ~ dnorm(0, 1)T(0, )
  mub0s  ~ dnorm(0, 0.1)

  taub0s <- 1/(sdb0s * sdb0s)
  sdb0s ~ dnorm(0, 1)T(0, )
  }
")

Luego definimos los datos, la función para los valores iniciales y los parámetros que queremos guardar.

rotos <- datos$Rotos
freq <- datos$Frecuencia
n <- length(rotos)
N <- datos$Total

Nsitio <- as.numeric(datos$Sitio)
S <- max(Nsitio)

Nplanta <- datos$Planta
ns <- unique(Nsitio)

P <- max(Nplanta)
Ns <- aggregate(Nsitio, by = list(Nplanta), FUN = "max")[, 2]

data <- list("rotos", "freq", "N", "n", "S", "P", "Nplanta", "Ns")
params <- c("b0", "b1", "mub0", "sdb0", "mub0s", "sdb0s")

inits <- function() {
    list(b1 = rnorm(1, 0, 1), mub0 = rnorm(S, 0, 1), sdb0 = runif(1, 0, 10), 
        mub0s = rnorm(1, 0, 1), sdb0s = runif(1, 0, 5))
}

Ahora definimos cuántas iteraciones queremos correr, el “thining”, el “burn in”, número de cadenas:

ni <- 5000  # número de iteraciones
nt <- 1  # si descartamos 'thin' algunos valores
nb <- 2500  # cuantas iteraciones usamos de 'burn in'
nc <- 3  # y cuantas cadenas corremos

Finalmente llamamos a jags

m2.sim <- jags(data, inits, params, model.file = "mod2.bug", n.chains = nc, 
    n.iter = ni, n.burnin = nb, n.thin = nt)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 477
##    Unobserved stochastic nodes: 206
##    Total graph size: 2707
## 
## Initializing model
## 
## Adaptive phase..... 
## Adaptive phase complete 
##  
## 
##  Burn-in phase, 2500 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 2500 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(m2.sim)

## JAGS output for model 'mod2.bug', generated by jagsUI.
## Estimates based on 3 chains of 5000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 2500 iterations and thin rate = 1,
## yielding 7500 total samples from the joint posterior. 
## MCMC ran for 0.165 minutes at time 2019-03-02 13:15:23.
## 
##             mean     sd    2.5%     50%   97.5% overlap0     f  Rhat n.eff
## b0[1]     -7.938  1.467 -10.896  -7.865  -5.395    FALSE 1.000 1.014   170
## b0[2]     -3.730  0.681  -5.081  -3.715  -2.423    FALSE 1.000 1.017   122
## b0[3]     -7.655  1.300 -10.360  -7.619  -5.303    FALSE 1.000 1.012   203
## b0[4]     -7.781  1.512 -10.818  -7.609  -4.969    FALSE 1.000 1.148    21
## b0[5]     -8.050  1.754 -12.014  -7.833  -5.218    FALSE 1.000 1.110    24
## b0[6]     -7.277  1.539 -10.494  -7.168  -4.512    FALSE 1.000 1.039    61
## b0[7]     -7.296  1.543 -10.683  -7.171  -4.707    FALSE 1.000 1.013   222
## b0[8]     -7.919  1.779 -11.502  -7.766  -4.920    FALSE 1.000 1.047    51
## b0[9]     -7.272  1.699 -11.097  -7.078  -4.320    FALSE 1.000 1.015   156
## b0[10]    -7.550  1.477 -10.700  -7.464  -5.083    FALSE 1.000 1.047   234
## b0[11]    -8.295  1.658 -12.035  -8.159  -5.454    FALSE 1.000 1.040    58
## b0[12]    -7.332  1.534 -10.947  -7.130  -4.613    FALSE 1.000 1.127    36
## b0[13]    -7.044  1.340  -9.744  -6.970  -4.602    FALSE 1.000 1.071    49
## b0[14]    -7.181  1.501 -10.649  -7.101  -4.591    FALSE 1.000 1.116    29
## b0[15]    -7.134  1.376  -9.986  -7.002  -4.858    FALSE 1.000 1.137    21
## b0[16]    -7.355  1.863 -11.690  -7.010  -4.434    FALSE 1.000 1.127    24
## b0[17]    -6.930  1.570  -9.915  -6.865  -3.942    FALSE 1.000 1.043    52
## b0[18]    -6.623  1.285  -9.065  -6.620  -4.109    FALSE 1.000 1.039    65
## b0[19]    -3.578  0.685  -5.001  -3.556  -2.328    FALSE 1.000 1.015   179
## b0[20]    -4.255  0.960  -6.174  -4.220  -2.427    FALSE 1.000 1.010   201
## b0[21]    -4.656  0.911  -6.507  -4.617  -3.009    FALSE 1.000 1.009   268
## b0[22]    -6.898  1.715 -10.657  -6.625  -4.134    FALSE 1.000 1.135    20
## b0[23]    -6.765  1.577  -9.932  -6.711  -3.971    FALSE 1.000 1.063    42
## b0[24]    -6.576  1.768 -10.185  -6.470  -3.405    FALSE 1.000 1.064    38
## b0[25]    -7.008  1.644 -10.625  -6.946  -4.210    FALSE 1.000 1.044    73
## b0[26]    -7.260  1.762 -11.204  -7.037  -4.551    FALSE 1.000 1.092    34
## b0[27]    -5.416  0.948  -7.467  -5.384  -3.672    FALSE 1.000 1.055    44
## b0[28]    -6.030  1.020  -8.136  -5.984  -4.079    FALSE 1.000 1.020   104
## b0[29]    -4.927  0.898  -6.759  -4.931  -3.182    FALSE 1.000 1.051    46
## b0[30]    -7.087  1.154  -9.492  -7.052  -4.933    FALSE 1.000 1.150    18
## b0[31]    -2.073  1.392  -4.448  -2.204   0.847     TRUE 0.925 1.046    67
## b0[32]    -3.528  0.904  -5.223  -3.551  -1.686    FALSE 1.000 1.024   103
## b0[33]    -2.275  1.236  -4.562  -2.330   0.294     TRUE 0.962 1.062    37
## b0[34]    -7.692  1.465 -10.879  -7.523  -5.188    FALSE 1.000 1.104    25
## b0[35]    -4.095  0.913  -5.906  -4.086  -2.316    FALSE 1.000 1.041    63
## b0[36]    -7.803  1.325 -10.463  -7.773  -5.338    FALSE 1.000 1.106    25
## b0[37]    -8.142  1.262 -11.007  -8.072  -5.919    FALSE 1.000 1.025   103
## b0[38]    -6.173  1.615  -9.638  -6.142  -3.383    FALSE 1.000 1.130    26
## b0[39]    -6.710  1.759 -10.097  -6.779  -3.344    FALSE 1.000 1.302    11
## b0[40]    -6.326  1.656  -9.549  -6.317  -3.296    FALSE 1.000 1.045    50
## b0[41]    -3.685  1.382  -6.748  -3.539  -1.353    FALSE 1.000 1.010   232
## b0[42]    -6.437  1.801 -10.094  -6.236  -3.398    FALSE 1.000 1.067    43
## b0[43]    -5.760  1.904  -9.915  -5.544  -2.397    FALSE 1.000 1.128    22
## b0[44]    -4.195  1.336  -6.938  -4.082  -1.951    FALSE 1.000 1.047    49
## b0[45]    -6.069  1.852  -9.981  -6.057  -2.609    FALSE 1.000 1.085    31
## b0[46]    -6.062  1.740  -9.573  -5.966  -2.940    FALSE 1.000 1.058    40
## b0[47]    -6.115  1.098  -8.375  -6.151  -3.873    FALSE 1.000 1.056    41
## b0[48]    -9.327  1.427 -12.308  -9.286  -6.652    FALSE 1.000 1.089    29
## b0[49]    -7.528  1.141  -9.809  -7.508  -5.202    FALSE 1.000 1.049    47
## b0[50]    -5.830  1.107  -8.147  -5.804  -3.716    FALSE 1.000 1.032    80
## b0[51]    -7.953  1.181 -10.351  -7.949  -5.721    FALSE 1.000 1.064    41
## b0[52]    -8.528  1.228 -10.982  -8.498  -6.158    FALSE 1.000 1.047    50
## b0[53]    -7.549  1.153  -9.889  -7.524  -5.420    FALSE 1.000 1.051    46
## b0[54]    -7.185  1.107  -9.362  -7.189  -4.906    FALSE 1.000 1.047    50
## b0[55]    -5.972  1.274  -8.412  -5.996  -3.449    FALSE 1.000 1.035    65
## b0[56]    -5.341  1.651  -9.404  -5.153  -2.779    FALSE 1.000 1.081    67
## b0[57]    -5.646  1.730  -9.116  -5.497  -2.605    FALSE 1.000 1.031   103
## b0[58]    -5.809  1.582  -8.889  -5.767  -2.816    FALSE 1.000 1.044    53
## b0[59]    -5.729  1.490  -8.451  -5.862  -2.900    FALSE 1.000 1.009   274
## b0[60]    -5.616  1.538  -8.720  -5.606  -2.616    FALSE 1.000 1.227    13
## b0[61]    -3.921  1.172  -6.189  -3.909  -1.802    FALSE 1.000 1.003  1133
## b0[62]    -5.598  1.766  -8.954  -5.514  -2.587    FALSE 1.000 1.236    13
## b0[63]    -5.880  1.675  -9.258  -5.848  -2.846    FALSE 1.000 1.059    39
## b0[64]    -3.841  1.236  -6.182  -3.777  -1.699    FALSE 1.000 1.013   678
## b0[65]    -5.940  1.687  -9.109  -5.988  -2.620    FALSE 1.000 1.052    46
## b0[66]    -5.477  1.630  -8.588  -5.413  -2.602    FALSE 1.000 1.137    21
## b0[67]    -5.911  1.679  -9.011  -5.927  -2.735    FALSE 1.000 1.027    88
## b0[68]    -5.446  1.446  -8.370  -5.351  -2.983    FALSE 1.000 1.142    20
## b0[69]    -5.537  1.469  -8.753  -5.405  -2.923    FALSE 1.000 1.030    85
## b0[70]    -2.350  0.727  -3.839  -2.329  -0.977    FALSE 0.999 1.006  1229
## b0[71]    -4.672  1.158  -7.306  -4.590  -2.731    FALSE 1.000 1.024    99
## b0[72]    -2.239  0.891  -4.074  -2.230  -0.514    FALSE 0.994 1.008   319
## b0[73]    -2.631  0.640  -3.937  -2.617  -1.436    FALSE 1.000 1.003   949
## b0[74]    -4.622  1.050  -6.757  -4.570  -2.745    FALSE 1.000 1.067    36
## b0[75]    -6.021  1.434  -9.022  -5.916  -3.483    FALSE 1.000 1.008  2784
## b0[76]    -2.177  0.913  -4.001  -2.182  -0.393    FALSE 0.989 1.001  3287
## b0[77]    -5.542  1.492  -8.785  -5.452  -2.909    FALSE 1.000 1.023   123
## b0[78]    -6.058  1.553  -9.378  -5.827  -3.528    FALSE 1.000 1.024   108
## b0[79]    -5.926  1.484  -8.998  -5.844  -3.330    FALSE 1.000 1.063    43
## b0[80]    -3.273  0.666  -4.606  -3.259  -2.050    FALSE 1.000 1.023    95
## b0[81]    -5.600  1.482  -8.833  -5.438  -3.169    FALSE 1.000 1.012   351
## b0[82]    -8.242  1.130 -10.442  -8.259  -6.055    FALSE 1.000 1.073    34
## b0[83]    -5.058  1.040  -7.162  -5.061  -3.035    FALSE 1.000 1.034    88
## b0[84]    -4.380  1.230  -6.722  -4.422  -1.774    FALSE 0.998 1.022    95
## b0[85]    -2.848  1.495  -5.569  -2.902   0.204     TRUE 0.966 1.064    37
## b0[86]    -6.192  1.054  -8.325  -6.163  -4.228    FALSE 1.000 1.026   114
## b0[87]    -5.393  1.073  -7.601  -5.384  -3.332    FALSE 1.000 1.097    26
## b0[88]    -8.748  1.390 -11.687  -8.674  -6.102    FALSE 1.000 1.133    21
## b0[89]    -5.395  1.131  -7.713  -5.383  -3.229    FALSE 1.000 1.067    39
## b0[90]    -3.503  1.124  -5.811  -3.520  -1.259    FALSE 0.999 1.071    33
## b0[91]    -5.647  1.033  -7.740  -5.629  -3.712    FALSE 1.000 1.084    29
## b0[92]    -5.425  1.084  -7.623  -5.425  -3.315    FALSE 1.000 1.047    50
## b0[93]    -8.964  1.465 -12.004  -8.830  -6.384    FALSE 1.000 1.020   183
## b0[94]    -3.130  1.293  -5.528  -3.167  -0.377    FALSE 0.987 1.023   118
## b0[95]    -4.360  1.140  -6.580  -4.279  -2.385    FALSE 1.000 1.086    29
## b0[96]    -5.176  1.510  -9.321  -4.963  -2.853    FALSE 1.000 1.078    75
## b0[97]    -4.957  1.425  -8.053  -4.820  -2.525    FALSE 1.000 1.086    30
## b0[98]    -4.950  1.479  -8.311  -4.869  -2.350    FALSE 1.000 1.056    41
## b0[99]    -5.110  1.574  -9.014  -4.893  -2.480    FALSE 1.000 1.042    66
## b0[100]   -4.580  1.550  -7.940  -4.488  -1.709    FALSE 1.000 1.049    53
## b0[101]   -3.266  1.165  -5.745  -3.257  -1.066    FALSE 1.000 1.014   414
## b0[102]   -3.592  1.192  -6.234  -3.530  -1.530    FALSE 1.000 1.010   541
## b0[103]   -5.071  1.414  -7.795  -5.056  -2.353    FALSE 1.000 1.023   168
## b0[104]   -5.150  1.569  -8.898  -4.841  -2.680    FALSE 1.000 1.107    29
## b0[105]   -5.302  1.413  -8.146  -5.249  -2.724    FALSE 1.000 1.090    27
## b0[106]   -1.914  0.723  -3.339  -1.905  -0.547    FALSE 0.999 1.009   269
## b0[107]   -0.539  0.531  -1.592  -0.524   0.471     TRUE 0.844 1.042    53
## b0[108]   -5.527  1.454  -8.525  -5.524  -2.800    FALSE 1.000 1.013  1192
## b0[109]   -5.726  1.607  -8.982  -5.592  -3.084    FALSE 1.000 1.263    12
## b0[110]   -5.592  1.507  -8.712  -5.576  -2.859    FALSE 1.000 1.023   120
## b0[111]   -5.248  1.485  -8.391  -5.237  -2.312    FALSE 1.000 1.067    36
## b0[112]   -1.545  0.643  -2.813  -1.514  -0.349    FALSE 0.996 1.020   111
## b0[113]   -1.539  0.613  -2.844  -1.486  -0.480    FALSE 0.999 1.008   298
## b0[114]   -5.029  1.435  -7.801  -4.997  -2.350    FALSE 1.000 1.059    39
## b0[115]   -5.391  1.571  -8.396  -5.391  -2.469    FALSE 1.000 1.033    75
## b0[116]   -3.554  1.074  -5.891  -3.458  -1.657    FALSE 1.000 1.087    30
## b0[117]   -5.814  1.840  -9.344  -5.898  -2.580    FALSE 1.000 1.309    11
## b0[118]   -5.642  1.805  -9.212  -5.559  -2.525    FALSE 1.000 1.084    40
## b0[119]   -3.248  0.741  -4.830  -3.218  -1.939    FALSE 1.000 1.030   109
## b0[120]   -4.394  0.956  -6.421  -4.295  -2.754    FALSE 1.000 1.116    24
## b0[121]   -3.059  0.728  -4.627  -3.020  -1.794    FALSE 1.000 1.013   361
## b0[122]   -5.952  1.538  -9.331  -5.822  -3.225    FALSE 1.000 1.071    54
## b0[123]   -5.990  1.442  -8.958  -5.912  -3.572    FALSE 1.000 1.234    13
## b0[124]   -5.380  1.602  -9.021  -5.162  -2.751    FALSE 1.000 1.020   176
## b0[125]   -4.407  0.968  -6.512  -4.355  -2.710    FALSE 1.000 1.019   139
## b0[126]   -3.587  0.939  -5.808  -3.473  -2.143    FALSE 1.000 1.080    56
## b0[127]   -5.341  1.223  -8.139  -5.247  -3.094    FALSE 1.000 1.019   142
## b0[128]   -2.411  0.594  -3.633  -2.388  -1.314    FALSE 1.000 1.005   411
## b0[129]   -5.721  1.199  -8.319  -5.604  -3.716    FALSE 1.000 1.034    93
## b0[130]   -1.453  0.645  -2.732  -1.447  -0.229    FALSE 0.987 1.014   150
## b0[131]   -3.687  0.890  -5.771  -3.635  -2.147    FALSE 1.000 1.067    38
## b0[132]   -3.746  0.838  -5.506  -3.720  -2.220    FALSE 1.000 1.008   528
## b0[133]   -7.215  1.639 -10.418  -7.168  -4.339    FALSE 1.000 1.019   197
## b0[134]   -3.789  0.757  -5.339  -3.768  -2.356    FALSE 1.000 1.009   257
## b0[135]   -7.638  1.707 -11.408  -7.546  -4.787    FALSE 1.000 1.068    70
## b0[136]   -5.423  1.250  -7.984  -5.288  -3.323    FALSE 1.000 1.104    25
## b0[137]   -7.462  1.448 -10.537  -7.376  -5.024    FALSE 1.000 1.132    22
## b0[138]   -7.069  1.517 -10.693  -6.901  -4.537    FALSE 1.000 1.125    22
## b0[139]   -3.864  0.723  -5.293  -3.859  -2.516    FALSE 1.000 1.026    89
## b0[140]   -7.816  1.741 -11.007  -7.675  -4.853    FALSE 1.000 1.411     9
## b0[141]   -7.376  1.550 -10.800  -7.233  -4.716    FALSE 1.000 1.056    48
## b0[142]   -7.530  1.508 -10.765  -7.493  -4.813    FALSE 1.000 1.041   103
## b0[143]   -6.913  1.499 -10.134  -6.870  -4.267    FALSE 1.000 1.056    54
## b0[144]   -7.441  1.420 -10.228  -7.441  -4.767    FALSE 1.000 1.003   864
## b0[145]   -4.523  1.085  -6.658  -4.473  -2.460    FALSE 1.000 1.026   101
## b0[146]   -6.375  1.084  -8.487  -6.375  -4.209    FALSE 1.000 1.044    53
## b0[147]   -4.072  1.206  -6.385  -4.078  -1.812    FALSE 1.000 1.032    71
## b0[148]   -1.955  1.534  -4.672  -2.072   1.545     TRUE 0.898 1.115    23
## b0[149]   -4.033  1.131  -6.271  -4.021  -1.914    FALSE 1.000 1.071    34
## b0[150]   -4.153  1.194  -6.377  -4.195  -1.607    FALSE 1.000 1.025    90
## b0[151]   -3.323  1.391  -6.030  -3.361  -0.584    FALSE 0.988 1.012   218
## b0[152]   -2.491  1.489  -5.316  -2.485   0.688     TRUE 0.952 1.006   698
## b0[153]   -4.052  1.282  -6.500  -4.072  -1.382    FALSE 1.000 1.027    81
## b0[154]   -3.311  1.308  -5.675  -3.410  -0.538    FALSE 0.986 1.019   111
## b0[155]   -5.964  1.172  -8.403  -5.874  -3.867    FALSE 1.000 1.019   112
## b0[156]   -3.047  0.751  -4.502  -3.031  -1.620    FALSE 0.999 1.035   101
## b0[157]   -3.805  0.850  -5.509  -3.798  -2.167    FALSE 1.000 1.030    72
## b0[158]   -7.313  1.333 -10.546  -7.233  -4.937    FALSE 1.000 1.040   129
## b0[159]   -6.549  1.108  -8.772  -6.526  -4.533    FALSE 1.000 1.022   103
## b0[160]   -5.673  0.904  -7.605  -5.637  -3.947    FALSE 1.000 1.032   244
## b0[161]   -4.648  0.869  -6.389  -4.636  -2.881    FALSE 1.000 1.009   221
## b0[162]   -5.790  0.870  -7.653  -5.714  -4.279    FALSE 1.000 1.048    49
## b0[163]   -6.857  1.279  -9.633  -6.784  -4.538    FALSE 1.000 1.008   834
## b0[164]   -3.866  0.740  -5.279  -3.888  -2.384    FALSE 1.000 1.025   119
## b0[165]   -6.742  1.530 -10.358  -6.606  -4.115    FALSE 1.000 1.202    15
## b0[166]   -5.971  1.247  -8.733  -5.848  -3.748    FALSE 1.000 1.019   155
## b0[167]   -1.622  1.736  -4.997  -1.549   1.552     TRUE 0.814 1.147    18
## b0[168]   -2.958  1.453  -5.726  -2.988   0.061     TRUE 0.973 1.099    25
## b0[169]   -2.072  1.787  -5.351  -2.195   1.786     TRUE 0.878 1.080    31
## b0[170]   -5.017  1.154  -7.368  -4.995  -2.832    FALSE 1.000 1.044    53
## b0[171]   -2.038  1.659  -5.065  -2.124   1.371     TRUE 0.881 1.042    55
## b0[172]   -3.354  1.464  -6.044  -3.432  -0.369    FALSE 0.985 1.011   215
## b0[173]   -2.128  1.835  -5.442  -2.149   1.473     TRUE 0.870 1.264    12
## b0[174]   -2.284  1.664  -5.238  -2.400   1.345     TRUE 0.906 1.008   474
## b0[175]   -1.838  1.696  -4.903  -1.958   1.901     TRUE 0.858 1.010   231
## b0[176]   -3.323  1.506  -6.163  -3.382  -0.206    FALSE 0.980 1.166    17
## b0[177]   -2.057  1.797  -5.255  -2.195   1.735     TRUE 0.870 1.036    62
## b0[178]   -2.334  1.982  -6.044  -2.442   1.663     TRUE 0.865 1.034    75
## b0[179]   -5.598  1.177  -7.835  -5.629  -3.211    FALSE 1.000 1.053    47
## b0[180]   -4.767  1.241  -7.200  -4.766  -2.288    FALSE 1.000 1.063    50
## b0[181]   -4.762  1.272  -7.302  -4.748  -2.250    FALSE 1.000 1.032    70
## b0[182]   -5.213  1.283  -7.871  -5.200  -2.732    FALSE 1.000 1.035    71
## b0[183]   -2.741  1.552  -5.577  -2.840   0.387     TRUE 0.953 1.066    36
## b0[184]   -3.286  1.395  -6.000  -3.317  -0.475    FALSE 0.991 1.035    72
## b0[185]   -4.096  1.347  -6.943  -4.064  -1.415    FALSE 1.000 1.040    81
## b0[186]   -6.736  1.229  -9.141  -6.749  -4.329    FALSE 1.000 1.054    45
## b1         2.165  0.341   1.500   2.161   2.854    FALSE 1.000 1.064    38
## mub0[1]   -5.327  0.909  -7.078  -5.321  -3.537    FALSE 1.000 1.060    38
## mub0[2]   -4.467  0.729  -5.919  -4.449  -3.099    FALSE 1.000 1.010   258
## mub0[3]   -5.600  0.892  -7.401  -5.588  -3.889    FALSE 1.000 1.165    17
## mub0[4]   -5.471  0.778  -7.011  -5.462  -3.967    FALSE 1.000 1.035    66
## mub0[5]   -5.547  0.998  -7.573  -5.524  -3.585    FALSE 1.000 1.065    38
## mub0[6]   -6.406  0.743  -7.881  -6.404  -4.961    FALSE 1.000 1.016   158
## mub0[7]   -6.991  0.786  -8.554  -6.986  -5.463    FALSE 1.000 1.028    84
## mub0[8]   -4.311  0.639  -5.586  -4.304  -3.086    FALSE 1.000 1.021   109
## mub0[9]   -4.699  0.725  -6.151  -4.687  -3.327    FALSE 1.000 1.034    78
## mub0[10]  -6.220  0.794  -7.866  -6.209  -4.712    FALSE 1.000 1.068    35
## mub0[11]  -2.992  1.176  -5.165  -3.042  -0.539    FALSE 0.992 1.064    37
## mub0[12]  -6.897  1.040  -9.067  -6.866  -4.886    FALSE 1.000 1.044    53
## mub0[13]  -4.310  0.632  -5.587  -4.295  -3.098    FALSE 1.000 1.017   170
## mub0[14]  -5.357  0.804  -6.884  -5.381  -3.733    FALSE 1.000 1.065    37
## mub0[15]  -4.076  1.056  -6.158  -4.101  -1.963    FALSE 1.000 1.030    73
## mub0[16]  -4.752  1.128  -7.007  -4.775  -2.469    FALSE 1.000 1.038    65
## sdb0       1.838  0.166   1.520   1.837   2.164    FALSE 1.000 1.161    17
## mub0s     -5.156  0.638  -6.425  -5.155  -3.902    FALSE 1.000 1.056    44
## sdb0s      1.300  0.305   0.777   1.270   1.985    FALSE 1.000 1.010   214
## deviance 827.776 17.743 795.278 826.845 864.091    FALSE 1.000 1.038    64
## 
## **WARNING** Rhat values indicate convergence failure. 
## 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 = 152.4 and DIC = 980.218 
## DIC is an estimate of expected predictive error (lower is better).

Podemos hacer un gráfico un poco sobrepoblado con las ordenadas al origen por sitio (en negro) y por planta (en gris)

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(density(m2.sim$sims.list$mub0s), ylim = c(0, 1), xlim = c(-12, 3), lwd = 3, 
    main = "", xlab = "b0")
ns <- unique(Nsitio)
for (i in 1:S) {
    lines(density(m2.sim$sims.list$mub0[, i]), ylim = c(0, 2), xlim = c(-12, 
        3), lwd = 2)
    for (j in which(Ns == ns[i])) {
        lines(density(m2.sim$sims.list$b0[, j]), col = "gray")
    }
}

par(op)

Ahora hagamos una comparación más interesante entre las posteriores de las pendientes en los dos modelos:

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(density(m1.sim$sims.list$b1), lwd = 3, main = "", xlab = "b1", xlim = c(0, 
    4))
lines(density(m2.sim$sims.list$b1), lwd = 3, col = 2)

par(op)

Ejercicio:

Ajusten modelos equivalentes a los dos anteriores pero asumiendo alguna relación entre visitas y rotura de pistilos que no use la transformación logit.

---

Juan Manuel Morales. 30 de Septiembre del 2015. Última actualización: 2019-03-02