Você está aqui: start » publications » papercompanions » ejpjmstaa

Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

Link para esta página de comparações

Ambos lados da revisão anteriorRevisão anterior
Próxima revisão		Revisão anterior
publications:papercompanions:ejpjmstaa [2008/05/09 05:56] – ernestopublications:papercompanions:ejpjmstaa [2008/05/09 06:39] (atual) – ernesto
Linha 1: Linha 1:
 ==== Jardim, E. and Ribeiro Jr., P.J. Submitted. Modelling Spatio-temporal Abundance at Age with Bayesian Geostatistics and Compositional Data Analysis. Canadian Journal of Fisheries and Aquatic Science ==== ==== Jardim, E. and Ribeiro Jr., P.J. Submitted. Modelling Spatio-temporal Abundance at Age with Bayesian Geostatistics and Compositional Data Analysis. Canadian Journal of Fisheries and Aquatic Science ====
  
-=== Paper companion ===+==== Paper companion ====
  
-== R code ==+=== R code ===
  
 +<code R>
 +require(geoR)
 +require(MASS)
 +require(FLCore)
  
 +############################################################################
 +#              => DESIGN ESTIMATES <=                                      #
 +############################################################################
 +
 +inds.nhat <- matrix(inds.orig$nhat[inds.orig$age<6], ncol=18, nrow=6)
 +inds.nhat <- FLQuant(inds.nhat, dimnames=list(age=0:5, year=1989:2006, unit="unique", season="unique", area="unique"))
 +inds.varnhat <- matrix(inds.orig$varnhat[inds.orig$age<6], ncol=18, nrow=6)
 +inds.varnhat <- FLQuant(inds.varnhat, dimnames=list(age=0:5, year=1989:2006, unit="unique", season="unique", area="unique"))
 +inds.cv <- sqrt(inds.varnhat)/inds.nhat
 +
 +############################################################################
 +#              => COMPOSITIONAL DATA ANALYSIS <=                           #
 +############################################################################
 +require(compositions)
 +
 +# data (age 2 last)
 +data.comp <- split(cd8906[,c(19,20,22,23,24,21)], cd8906$y)
 +data.acomp <- lapply(data.comp, "+", 0.1)
 +
 +# compute compositions
 +data.acomp <- lapply(data.acomp, acomp)
 +
 +# parametric bootstrap
 +data.comps <- lapply(data.acomp, function(x){
 + xt <- alr(x)[] 
 + mu <- apply(xt,2,mean)
 + sig <- sqrt(apply(xt,2,var)/nrow(xt))
 + rho <- cor(xt)
 + sigmu <- rho*sig%*%t(sig)
 + mat <- mvrnorm(1000, mu, sigmu, empirical=T)
 + alr.inv(mat)
 +})
 +
 +detach()
 +
 +############################################################################
 +#              => GEOSTATISTICS <=                                         #
 +############################################################################
 +
 +#===========================================================================
 +# GLM
 +#===========================================================================
 +
 +# data object
 +cd8906.glm <- transform(cd8906, 
 + n05h=a0+a1+a2+a3+a4+a5, 
 + n05o=(a0+a1+a2+a3+a4+a5)*dur/60, 
 + n05d=((a0+a1+a2+a3+a4+a5)*dur/60)/dist)
 +
 +cd8906.glm <- cd8906.glm[,c(4,5,7,8,44:48)]
 +
 +# glm
 +n05h.glm <- glm(round(n05h,0)~factor(y)*factor(dcode), family=negative.binomial(0.5), data=cd8906.glm)
 +
 +# residuals and predictions
 +n05h.res <- residuals(n05h.glm)
 +newdata <- expand.grid(y=factor(1989:2006), dcode=factor("0"))
 +n05h.pred <- predict(n05h.glm, newdata)
 +names(n05h.pred) <- 1989:2006
 +n05h.pred <- apply(cd8906.glm,1,function(x) n05h.pred[x["y"]]) + n05h.res
 +
 +# check normality
 +sum(unlist(lapply(tapply(log(cd8906.glm$n05h+0.1), cd8906.glm$y, shapiro.test),"[","p.value"))>0.01)
 +sum(unlist(lapply(tapply(n05h.pred, cd8906.glm$y, shapiro.test),"[","p.value"))>0.01)
 +
 +#===========================================================================
 +# GEOSTATS
 +#===========================================================================
 +
 +# geodata
 +n05hp.gd <- mkgd(cbind(cd8906.glm[,c("y","Xit","Yit")], n05h.pred))
 +
 +# likfit
 +n05hp.lf <- lapply(n05hp.gd, likfit, messages=F, ini.cov.pars=expand.grid(seq(0,2,len=3),seq(0,40,len=5)))
 +
 +# bayes 
 +require(VGAM)
 +PC <- prior.control(
 + phi.prior = "exponential", phi=20, phi.discrete= seq(0,100, l=20),
 + tausq.rel.prior = dzipois(0:8,1.25,0.25)/sum(dzipois(0:8,1.25,0.25)), tausq.rel.discrete = seq(0,2, l=9)
 +)
 +detach()
 +
 +MC <- model.control()
 +OC <- output.control(n.predictive = 1000, messages=FALSE)
 +kbtg <- split(as.character(1989:2006), as.character(1989:2006))
 +
 +n05hp.kb <- lapply(kbtg, function(x){
 + v <- date()
 + cat(x,":", v, "\n")
 + gd <- n05hp.gd[[x]]
 + krige.bayes(gd, prior=PC, model=MC, output=OC)
 +})
 +
 +#===========================================================================
 +# SPATIAL PREDICTIONS: VERY TIME DEMANDING
 +#===========================================================================
 +
 +# non rotated data
 +gd.lst <- mkgd(cbind(y=cd8906.glm$y, cd8906[,c(30,31)], n05h.pred))
 +OC <- output.control(n.predictive=1, messages=F)
 +
 +n05hp.kr <- split(as.character(1989:2006), as.character(1989:2006))
 +n05hp.kr <- lapply(n05hp.kr, function(x){
 + cat(x, sep=",")
 + gd <- gd.lst[[x]]
 + pars <- n05hp.kb[[x]]$posterior$sample
 + v <- apply(pars, 1, function(y){
 + cat(".")
 + KC <- krige.control(cov.model = "matern", cov.pars = c(y[2],y[3]), kappa = 0.5, 
 + nugget = y[2]*y[4], beta=y[1])
 + krige.conv(gd, locations=gr, krige=KC, output=OC)
 + })
 + sims <- lapply(v,"[[","simulations")
 + sims <- do.call("cbind", sims)
 + list(locations=gr, pars=pars, simulations=sims)
 +})
 +
 +############################################################################
 +#              => ALL TOGETHER NOW ! ESTIMATE ABUNDANCE AT AGE <=          #
 +############################################################################
 +
 +#===========================================================================
 +# SIMULATIONS
 +#===========================================================================
 +
 +# initialize objects
 +iay.arr <- array(dim=c(1,18,1000), dimnames=list(age="all", year=1989:2006,sim=1:1000))
 +iaya.arr <- array(dim=c(6,18,1000), dimnames=list(age=0:5, year=1989:2006,sim=1:1000))
 +
 +# beta simulations
 +n05hp.iay <- lapply(n05hp.kb, function(x) exp(x$posterior$sample$beta))
 +n05.iay <- iay.arr
 + for(i in names(n05hp.iay)) n05.iay[,i,] <- n05hp.iay[[i]]
 +
 +# abundance at age simulations
 +lst <- split(as.character(1989:2006), as.character(1989:2006))
 +n05hp.iaya <- lapply(lst, function(x) n05hp.iay[[x]]*data.comps[[x]][,c(1,2,6,3,4,5)])
 +
 +n05.iaya <- iaya.arr
 + for(i in names(n05hp.iaya)) n05.iaya[,i,] <- t(n05hp.iaya[[i]])
 +
 +#===========================================================================
 +# SUMMARY STATS & ABUNDANCE INDEX
 +#===========================================================================
 +
 +n05.iay.med <- apply(n05.iay, c(1,2), median)
 +n05.iay.mad <- apply(n05.iay, c(1,2), mad, constant=1)
 +n05.iay.madmed <- n05.iay.mad/n05.iay.med
 +
 +n05.iaya.med <- apply(n05.iaya, c(1,2), median)
 +n05.iaya.mad <- apply(n05.iaya, c(1,2), mad, constant=1)
 +n05.iaya.madmed <- n05.iaya.mad/n05.iaya.med
 +
 +############################################################################
 +#              => THAT'S ALL FOLKS ! <=                                    #
 +############################################################################
 +
 +</code>
  
  
  

QR Code
QR Code publications:papercompanions:ejpjmstaa (generated for current page)