An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20130207/0c0ff7fd/attachment.pl>
MCMCglmm - ZIP model for mackerel egg data
2 messages · Gastauer, Sven, Jarrod Hadfield
1 day later
Hi Sven, Quoting "Gastauer, Sven" <sven.gastauer at wur.nl> on Thu, 7 Feb 2013 16:11:49 +0000:
Dear all, I am trying to model mackerel egg abundance in the North Sea. The egg abundance information comes from triennial surveys and contains a lot of 0 values (e.g. samples without any mackerel eggs). The idea was hence to build a zero-inflated or a zero-altered model using MCMCglmm. Fixed terms should be Sea Surface Temperature (sst), Longitude, Latitude and years as random factor. The first problem or question I ran into is that of course the number of eggs within a sample is a simple integer, although normally I would define an offset : Egg abundance = Counted Eggs x Correction Factor x Sampling Depth / Filtered Volume (Ce*R*SD/VolFilt), but as far as I understand it this is not an option in MCMCglmm? My easy solution for now was to simply round volume the corrected egg counts, which introduces a minor error.
You can fit an offset by fixing the regression coefficient at 1 in the prior. For example if the offset variable is associated with the 2nd fixed effect out of 3 then: prior<-list(B=list(V=diag(3)*1e+8, mu=c(0,1,0)), ...) prior$B$V[2,2]<-1e-8 or something similar, achieves this.
Based on some online readings etc. I came up with a model and a
prior, but I am not feeling very confident about it, could you
please indicate if this model makes sense?
prior0 <- list(B=list(mu = matrix(0,10,2),
V = diag(10)*1e+6),
R = list(V=diag(2),n=2,fix=2),
G = list(G1=list(V=diag(c(1,
1e-6)),n=2, fix=2)))
m.zip <- MCMCglmm(AE ~ trait + trait:(LONGITUDE*LATITUDE) + trait:sst,
random = ~idh(trait):YEAR,
family = "zipoisson",
rcov = ~ idh(trait):units,
prior = prior0,
nitt = 100000,
burnin = 10000,
thin = 25,
data = eggs)
m.zap <-MCMCglmm(AE ~ trait + trait:(LONGITUDE*LATITUDE) + trait:sst,
random = ~idh(trait):YEAR,
family = "zapoisson",
rcov = ~ idh(trait):units,
prior = prior0,
nitt = 100000,
burnin = 10000,
thin = 25, data = eggs)
If you only want to fit random effects for the count process random = ~idh(at.level(trait,1)):YEAR is probably a better way of doing it. This just fits a single variance rather than a 2x2 covariance matrix. Your prior degrees of freedom are large. With idh structures you have nu=2 on a single variance. I would use something smaller (e.g. 0.002) or preferably parameter expanded priors.
Plotting and looking at the summary both models works fine, but I cannot figure out how to predict the data. I tried the following, but get an Error message (which I do not understand): predict(m.zip, marginal = ~idh(trait):YEAR, type = "response", interval = "confidence") Error in M[, which(rm.v), ] <- 0 : incorrect number of subscripts
the predict method does not work with ZIP models yet, so you will have to do it by hand I'm afraid. The expected value is (1-pi)*lambda where pi is the probability of being zero from the zero-inflation part of the model, and lambda is the Poisson rate (in the CourseNotes notation: (1-plogis(l_2))*exp(l_2)). I am not sure what the expectation would be after marginalising the random effects/overdispersion, particularly if the two processes are correlated (which they are not in your model). Cheers, Jarrod
Any help is greatly appreciated and many thanks in advance, Sven Gastauer IMARES Fisheries Acoustician Email: sven.gastauer at wur.nl<mailto:sven.gastauer at wur.nl> Mobile: +31 (0)61 005 71 03 P.O. Box 68 Haringkade, 1 (0.23) 1970AB IJMUIDEN [[alternative HTML version deleted]]
_______________________________________________ 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.