[R-sig-ME] MCMCglmm specification

[Ray Danner]

Jarrod, Thanks very much for the help and again, for creating this package!
Sincerely,
Ray

On Tue, Feb 1, 2011 at 11:21 AM, Jarrod Hadfield <j.hadfield at ed.ac.uk> wrote:
> 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.
>
> _______________________________________________
> 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