# Sample posterior | |
post <- extract.samples(eggsFMod) | |
# PI of P(no clutch at all) | |
dens(logistic(post$ap), show.HPDI = T, xlab = "ZIP Bernoulli(p)") | |
# Run simulations w/ averages of all predictors, except parasite 0 / 1 | |
lambdaNoP <- exp(post$a + 0*post$bP + 0*post$bA + | |
0*post$bGS + 0*post$bES + 0*0*post$bPA) | |
simFledgeNoPar <- rpois(n = length(lambdaNoP), lambda = lambdaNoP) | |
lambdaP <- exp(post$a + 1*post$bP + 0*post$bA + | |
0*post$bGS + 0*post$bES + 1*0*post$bPA) | |
simFledgePar <- rpois(n = length(lambdaP), lambda = lambdaP) | |
table(simFledgeNoPar) | |
table(simFledgePar) | |
# Simulate with varying age | |
rangeA <- seq(-3, 3, length.out = 100) | |
# No parasite | |
predictions <- sapply(rangeA, function(x){ | |
exp(post$a + 0*post$bP + x*post$bA + 0*post$bGS + | |
0*post$bES + 0*x*post$bPA) | |
}) | |
hdpiPois <- apply(predictions, 2, HPDI, prob = .95) | |
meanPois <- colMeans(predictions) | |
plot(rangeA, meanPois, type = "l", ylim = c(0, 3), yaxp = c(0, 3, 3), | |
xlab = "Min Age (standardized)", ylab = expression(lambda)) | |
shade(hdpiPois, rangeA) | |
# Parasite | |
predictionsP <- sapply(rangeA, function(x){ | |
exp(post$a + 1*post$bP + x*post$bA + 0*post$bGS + | |
0*post$bES + x*post$bPA) | |
}) | |
hdpiPoisP <- apply(predictionsP, 2, HPDI, prob = .95) | |
meanPoisP <- colMeans(predictionsP) | |
lines(rangeA, meanPoisP, lty = 2, col = "red") | |
shade(hdpiPoisP, rangeA, col = rgb(1,0,0,.25)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment