Dear MCMCglmm users,
I recently started using MCMCglmm (very cool!) for zero-inflated
poisson models and have been reading Hadfield's Jstat paper,
CourseNotes, and Tutorial. However, I'm still uncertain about some of
the syntax and details regarding prior specification and would love
suggestions if possible.
I've been working with a data set of insect behavioral response
(backflips) to simulated predator attacks. Each of 14 individual
insects (id) was exposed to 4 different predators (predator.type).
The "excess" zeros come from individuals that did not respond to any
predators and a general negative gradient of response from predator 1
to predator 4. As each individual was tested 4 times, I would like id
to be a random effect. My hypothesis is that predator #1 would elicit
the highest number of backflips from the insects. I would like to
compare response between predator.types based on overlap of 95%
credible intervals (unless there is a better method).
A prior sensitivity analysis showed that the conclusions are robust to
R-structure variance values from at least 0.01 to 5. I fixed the
G-structure variance to 0. I'm uncertain how to determine the degree
of freedom parameter (n) and am not certain that my syntax is correct
overall.
For this example, I dialed the number of iterations down, though it
would be nice to have an experienced eye examine the trace and density
plots.
Again, any advice regarding this analysis would be helpful.
Thanks in advance.
Ray
#Create dataset
id<-c(11, 4, 9, 5, 6, 12, 3,11, 4, 2, 6, 8, 7, 13, 12, 3, 9, 11, 4,
12, 6, 2, 5, 7, 13, 8, 6, 8, 7, 9, 13, 7, 11, 3, 2, 8, 5, 13, 4, 9, 1,
3, 10, 1, 2, 12, 1, 10, 5, 10, 1, 10)
predator.type<-factor(c(4, 4, 2, 4, 3, 3, 3, 1, 2, 2, 2, 4, 3, 4, 4,
1, 4, 3, 3, 2, 1, 4, 3, 4, 3, 3, 4, 2, 1, 3, 2, 2, 2, 4, 3, 1, 1, 1,
1, 1, 4, 2, 2, 3, 1, 1, 2, 3, 2, 1, 1, 4))
num.backflips<-c(0, 0, 0, 0, 0, 0, 0, 13, 0, 4, 0, 0, 0, 0, 0, 2, 0,
0, 0, 1, 1, 1, 1, 0, 6, 2, 0, 1, 4, 0, 0, 1, 0, 0, 2, 4, 9, 3, 2, 0,
0, 0, 1, 1, 9, 2, 0, 4, 3, 6, 0, 2)
avoidance<-data.frame(id, I(predator.type), num.backflips)
#ZIP model
library("MCMCglmm")
prior.1<-list(R = list(V = diag(2), n = 2, fix = 2), G = list(G1 =
list(V = diag(c(1, 1e-6)), n = 2, fix = 2)))
mod.1<-MCMCglmm(num.backflips ~ trait + trait:predator.type - 1,
random = ~idh(trait):id, rcov = ~idh(trait):units, data = avoidance,
prior = prior.1, family = "zipoisson", nitt = 130000, thin = 10,
burnin = 3000, verbose = FALSE, pl = TRUE)
plot(mod.1)
#summary(mod.1)
#Plot predicted values with 95% credible intervals
head(mod.1$Sol)
means<-c(mean(exp(mod.1$Sol[,1])), mean(exp(mod.1$Sol[,3])),
mean(exp(mod.1$Sol[,5])), mean(exp(mod.1$Sol[,7])))
ci1<-HPDinterval(exp(mod.1$Sol[,1]))
ci2<-HPDinterval(exp(mod.1$Sol[,3]))
ci3<-HPDinterval(exp(mod.1$Sol[,5]))
ci4<-HPDinterval(exp(mod.1$Sol[,7]))
l<-c(ci1[1], ci2[1], ci3[1], ci4[1])
u<-c(ci1[2], ci2[2], ci3[2], ci4[2])
library(gplots)
barplot2(means, beside = TRUE, plot.ci = TRUE, ci.l = l, ci.u=u,
border=c("black"), col = c("white", "white"), xpd = FALSE, axes =
TRUE, axisnames = TRUE, axis.lty = 1, ylim = c(0, 6), ylab =
"Predicted Backflips +/- 95% CRIs", xlab = "Predator Type", names.arg
= c("1", "2", "3", "4"), lwd = 2, ci.lwd = 2, ci.col = c("black"))