[R-sig-ME] MCMCglmm specification

Ray Danner rdanner at vt.edu
Thu Jan 27 02:01:55 CET 2011

Dear MCMCglmm users,

I recently started using MCMCglmm (very cool!) for zero-inflated
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.

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
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, ci2, ci3, ci4)
u<-c(ci1, ci2, ci3, ci4)
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"))