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 ! <=                                    #
############################################################################