[R-sig-ME] MCMCglmm specification

[Jarrod Hadfield]

Hi,

I think the model makes sense, although the priors you use are likely  
to be quite informative with so few individuals. idh structures with  
nu=2 is equivalent to fitting a prior to each individual variance with  
nu=2. (Note this changed in versions >=2.05 see Section 3.6.2 of the  
CourseNotes, because the old specification was creating confusion).  
Generally nu=2 is high, but especially if you have 14 individuals  
whose "effect" will be poorly resolved with 4 observations. Also,  
rather than fitting the id random effect for the zero-inflation but  
fixing the associated variance to zero, you could just leave it it  
out.  random = ~idh(at.level(trait,1)):id  -  it should ru a bit  
quicker.

I also noticed the traces look ed pretty bad for the parameters of the  
zero-inflated model. This seems to be a general problem with ZIP  
models (in MCMCglmm) particularly when zero-inflation is close to zero  
or one. Hurdle models, and zero-altered models fit similar models but  
they tend to mix much better. The zero-altered model in particular, is  
a nice way of assessing whether zero-inflation (or deflation) exists.

Cheers,

Jarrod



Quoting Ray Danner <rdanner at vt.edu>:

> 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"))
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>



-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.