[R-sig-ME] MCMCglmm specification

[Ray Danner]

Thanks Jarrod for the suggestions!  I am now testing if the data are
zero-inflated with the use of a zero-altered model.  It seems to me
that zero-inflation is low per level of predator.type.  My example is
a bit more complicated than that presented in the MCMCglmm
CourseNotes, as it has a random effect.  If anyone would be willing to
critique model specification or share rules of thumb for the degree of
zero-inflation that would warrant accounting for zero-inflation
(through a zero-inflated or hurdle model), I would appreciate it.

Thanks in advance. Ray

Zero-altered model based on data set d:

#Create data set
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)

#Set prior
prior<-list(R = list(V = diag(c(1, 1)), nu = 0.002, fix = 2), G =
list(G1 = list(V = diag(1), nu = 0.002, fix = 1)))

#Zero-altered model
fitza<-MCMCglmm(num.backflips ~ trait * predator.type, random =
~idh(at.level(trait,1)):id, rcov = ~idh(trait):units, data =
avoidance, prior = prior, family = "zapoisson", nitt = 130000, thin =
10, burnin = 3000, verbose = FALSE, pl = TRUE) #This takes ~1 minute
to compute.
#plot(fitza)
summary(fitza)

On Sat, Jan 29, 2011 at 9:58 AM, Jarrod Hadfield <j.hadfield at ed.ac.uk> wrote:
> 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.
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>



-- 
Raymond M. Danner
PhD Candidate, Virginia Tech and the Smithsonian Migratory Bird Center
PO Box 37012 MRC 5503, Washington, DC 20013-7012
(Fed-Ex and UPS deliveries: 3001 Connecticut Ave NW, Washington, DC 20008)
P 202-633-9444 | F 202-673-0040 | E rdanner at vt.edu
https://filebox.vt.edu/users/rdanner/
http://nationalzoo.si.edu/ConservationAndScience/Scientific_Staff/staff_scientists.cfm?id=290