[R-sig-ME] MCMCglmm specification

[Jarrod Hadfield]

Hi Ray,

I think you want to estimate the id level variance rather than fixing  
it at V so remove fix=1 in the prior for the G-structure.  Also, it is  
usual in zero-altered models to have the zero bit and the truncated  
poisson bit have the same over-dispersion. You do this by fitting the  
interaction rcov=~traits:units. I think it makes more sense to have  
the random term like this too (i.e. ~trait:id) rather than just  
fitting id effects for the truncated Poisson.

Something that does need to be clarified is whether the trait  
interactions are testing for zero inflation/deflation when the 2x2  
residual covariance matrix (for the binary and the truncated Poisson  
across observations) has the form:

v 0
0 v

(i.e. ~trait:units) where v is the estimated variance

or whether it actually it needs to be:

v v
v v

Perhaps someone is more familiar with zero-altered models than me?

Tentatively, I would say there's not much evidence for zero-infaltion  
in these data and the number of zeros can be explained by the  
treatment effects and/or the additional heterogeneity at the id and  
residual level. It seems like it is hard to assess whether the source  
of the additional heterogeneity is due to both id and observation,  or  
just one alone.

Code is below.

Cheers,

Jarrod

prior<-list(R = list(V = diag(1), nu = 0.002), G =
list(G1 = list(V = diag(1), nu = 0.002)))


#Zero-altered model
fitza1<-MCMCglmm(num.backflips ~ trait * predator.type, random =
~trait:id, rcov = ~trait:units, data = avoidance, prior = prior,  
family = "zapoisson", nitt = 130000, thin = 10, burnin = 3000, verbose  
= FALSE, pl = TRUE)

# Non of the traitza terms are "significant" suggesting that the  
standard overdispersed Poisson may be good enough.

fitza2<-MCMCglmm(num.backflips ~ trait * predator.type, rcov =  
~trait:units, data = avoidance, prior = prior, family = "zapoisson",  
nitt = 130000, thin = 10, burnin = 3000, verbose = FALSE, pl = TRUE)

# Qualitative similar results without id


#Poisson model with id as random
fitza3<-MCMCglmm(num.backflips ~ predator.type, random =
~id, data = avoidance, prior = prior, family = "poisson", nitt =  
130000, thin = 10, burnin = 3000, verbose = FALSE, pl = TRUE)

# the trace of the variances hit zero a bit, so perhaps run for longer  
and/or use parameter expansion.


#Poisson model
fitza4<-MCMCglmm(num.backflips ~ predator.type, data = avoidance,  
prior = prior, family = "poisson", nitt = 130000, thin = 10, burnin =  
3000, verbose = FALSE, pl = TRUE)

Quoting Ray Danner <rdanner at vt.edu>:

> 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
>
>



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