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

Próxima revisão		Revisão anterior
publications:papercompanions:ejpjmstaa [2008/05/09 05:55] – criada 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 ==== 
 + 
 +=== 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)