[R-sig-ME] Zipoisson MCMCglmm for abundance data

Daniel Sol dsolrueda at gmail.com
Tue Nov 27 15:13:29 CET 2012

Dear Jarrod,

many thanks for the rapid answer and very useful advice. Following
your advice, I'm gonna use the following code.

zi.prior <-  list(R = list(V = diag(2), n = 0.002, fix = 2),
	      G = list(G1 = list(V = 1, n = 0.002),
	      G2 = list(V = 1, n = 0.002),
                   G3 = list(V = 1, n = 0.002),
                   G4 = list(V = 1, n = 0.002)))

m2 <- MCMCglmm(abund.urb ~ trait-1 + at.level(trait,1):dens.surr,
			 random = ~idh(at.level(trait,1)):location	+						
              idh(at.level(trait,1)):animal +
				       idh(at.level(trait,1)):sp2 +
			                    rcov = ~ idh(trait):units,
				       prior = zi.prior,
				       data = dat0.phyl, family = "zipoisson", verbose = TRUE,
				       pr = FALSE, pl = FALSE)

However, I'm still wondering how to run the same model allowing the
relationship between abundance in the city (abund.urb) and density in
the surrounding habitats (dens.surr) to vary across locations
(location) in both intercepts and slopes.

Many thanks again,


2012/11/27 Jarrod Hadfield <j.hadfield at ed.ac.uk>:
> 1/3b You need to drop at.level(trait,1):location from the fixed model as you
> have it in the random part of the model (although this may just be a typo).
> I would also have trait-1 as you do not want the intercept for the Poisson
> process and the zero-inflation to be the same.
> 2. If you have many observations per species then I would put a
> non-phylogenetic species effect in too. If there are few (at the limit, only
> one) then it may be hard to separate the phylogenetic from the
> non-phylogenetic.
> 3a. this looks fine but make sure to put mesd in dat0.phyl (I presume this
> is the case otherwise MCMcglmm should spit an error, if it did not please
> tell me). Not sure how effort is measured, but you may not expect a linear
> relationship between 1/(effort) and the measurement error variance of the
> counts on the log scale. (I presume ** should be ^ in your code)
> 4. With  the 2x2 "idh" structure on the residuals I would use nu=0.002
> rather than nu=1.002. Only with a 2x2 "us" structure
> is the degree of belief for the marginal distribution of a single variance
> 0.002 when specifying nu=1.002. Parameter expanded priors might also be
> entertained for the random effect variances. They will also improve mixing
> if the varinaces are close to zero.

Daniel Sol
CREAF (Centre for Ecological Research and Applied Forestries)
CSIC (Centre for Advanced Studies of Blanes-Spanish National Research Council)
Autonomous University of Barcelona, Bellaterra, Catalonia E-08193, Spain
TEL: +34 93-5814678
FAX: +34 93-5814151
E-MAIL: d.sol at creaf.uab.es

More information about the R-sig-mixed-models mailing list