# Sample posterior | |
post <- extract.samples(eggsLMod) | |
# 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) | |
# Sim with varying age | |
# 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, 7), yaxp = c(0, 7, 7), | |
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