Simulation study

#################################################################################################
#
#	PAPER COMPANION: SIMULATION STUDY
#	Authors:	Ernesto Jardim <ernesto@ipimar.pt> 
#			Paulo Ribeiro Jr. <paulojus@ufpr.br>
#	Date: 20081013
#	Objective: 	Test paper methodology against design-based estimators regarding
# 			association between age compositions and the age aggregated abundance.
#			Tree options are simulated: (i) no association, (ii) week association, 
#			(iii) strong association.
#	Note:		The association is induced by considering correlated multivariate gaussian 
#			processes as the basis for building the age aggregated abundance and age 
#			compositions. 
#
#################################################################################################
 
# R libraries
library(geoR)
library(lattice)
library(MASS)
library(RandomFields)
library(compositions)
 
#==============================================================================================
# SIMULATING THE PROCESSES
#==============================================================================================
# grid
gs <- expand.grid((0:40)/40, (0:40)/40)
# size of processes
n <- nrow(gs)
# number of replicates
nsim <- 250
 
#----------------------------------------------------------------------------------------------
# STEP 1 : gaussian processes to build abundance and compositions
#----------------------------------------------------------------------------------------------
 
# object
arr0 <- array(NA, dim=c(n,7,nsim), dimnames=list(loc=1:n, age=c("all","1i","2i","1w","2w","1s","2s"), nsim=1:nsim))
# variance-covariance matrix
s2 <- 0.5
Sig <- diag(c(1,1,1,1,1,1,1))
Sig[1,4] <- 0.45; Sig[4,1] <- 0.45; Sig[1,6] <- 0.9; Sig[6,1] <- 0.9; Sig[6,4] <- 0.5; Sig[4,6] <- 0.5
Sig <- Sig*s2
# simulation
set.seed(111)
for(i in 1:nsim){
	arr0[,,i] <- mvrnorm(nrow(gs), c(1,0,0,0,0,0,0), Sig)
}
 
#----------------------------------------------------------------------------------------------
# STEP 2: generate 250 replicates of a log-gaussian spatial process (Diggle and Ribeiro Jr, 2007)
#         with  mu=1, phi=0.2, sigma2=0.5, tau2=0.5
#----------------------------------------------------------------------------------------------
 
# Generate a spatial Gaussian process Z
phi <- 0.2
sigmasq <- 0.5
set.seed(222)
Z <- grf(grid=gs, cov.pars=c(sigmasq, phi), nsim=nsim)
# build a log gausian process Y
Ysim <- Z
Ysim$data <- exp(Z$data+arr0[,1,])
# Y characteristics
Ysim.lmean <- apply(log(Ysim$data), 2, mean)
Ysim.lvar <- apply(log(Ysim$data), 2, var)
Ysim.lnhat <- exp(Ysim.lmean+Ysim.lvar/2)
 
#----------------------------------------------------------------------------------------------
# STEP 3: build compositions
#----------------------------------------------------------------------------------------------
 
# objects
arr00 <- array(1, dim=c(n,3,nsim), dimnames=list(loc=1:n, age=1:3, nsim=1:nsim))
Psim <- array(NA, dim=c(n,3,nsim,3), dimnames=list(loc=1:n, age=1:3, nsim=1:nsim, mucor=c("indep", "dep045","dep090")))
# option 1: no association between compositions and age aggregated abundance
arr00[,1:2,] <- exp(arr0[,2:3,])
Psim[,,,1] <- aperm(apply(arr00,c(1,3),function(x) x/sum(x)),c(2,1,3))
# option 2: week association between compositions and age aggregated abundance
arr00[,1:2,] <- exp(arr0[,4:5,])
Psim[,,,2] <- aperm(apply(arr00,c(1,3),function(x) x/sum(x)),c(2,1,3))
# option 3: strong association between compositions and age aggregated abundance
arr00[,1:2,] <- exp(arr0[,6:7,])
Psim[,,,3] <- aperm(apply(arr00,c(1,3),function(x) x/sum(x)),c(2,1,3))
# P characteristics
 
#----------------------------------------------------------------------------------------------
# STEP 4: build abundance at age 
#----------------------------------------------------------------------------------------------
 
# objects
Isim.ln <- array(NA, dim=c(n,3,nsim,3), dimnames=list(loc=1:n, age=1:3, nsim=1:nsim, mucor=c("indep", "dep045","dep090")))
# sim
for(i in 1:3){
	for(j in 1:3){
		Isim.ln[,i,,j] <- Psim[,i,,j]*Ysim$data
	}
}
# I characteristics
Isim.lnmean <- apply(log(Isim.ln), c(2,3,4), mean)
Isim.lnvar <- apply(log(Isim.ln), c(2,3,4), var)
Isim.lnhat <- exp(Isim.lnmean+Isim.lnvar/2)
 
#==============================================================================================
# ESTIMATION
#==============================================================================================
# number of samples to be drawn from each of the 250 replicate
ns <- 2
 
#----------------------------------------------------------------------------------------------
# STEP 5: estimation with methodology proposed by the paper
#----------------------------------------------------------------------------------------------
 
# objects
Yres <- array(NA, dim=c(ns,4,nsim), dimnames=list(samp=1:ns, stat=c("Ybar", "Ybarvar", "Yhat", "Yhatvar"), nsim=1:nsim))
Ihat <- array(NA, dim=c(ns,3,nsim,3), dimnames=list(samp=1:ns, age=1:3, nsim=1:nsim, mucor=c("indep", "dep045","dep090")))
xd <- expand.grid(dimnames(Ihat)[-c(1,2)])
# estimation
set.seed(222)
for(j in 1:ns){
cat("\n",j ,":", date(),"-", sep="")
	samp <- sample(1:n, 100)
	for(i in 1:nrow(xd)){
		x <- xd[i,]
		cat(".", sep="")
		# building the sample
		Isamp <- Isim.ln[samp,,x[[1]],x[[2]]]
		locsamp <- gs[samp,]
		Ysamp <- apply(Isamp,1,sum)
		# estimation of age aggregated abundance
		lf <- likfit(as.geodata(cbind(locsamp, Ysamp)), lambda=0, ini.cov.pars=c(1,0.2), messages=FALSE)
		Yhat <- exp(lf$beta+lf$tausq/2+lf$sigmasq/2)
		Yhatvar <- exp(2*lf$beta+lf$tausq+lf$sigmasq)*(exp(lf$tausq+lf$sigmasq)-1)
		Yres[j,c("Yhat", "Yhatvar"),x[[1]]] <- c(Yhat, Yhatvar)	
		# estimation of abundance-at-age
		prop <- acomp(Isamp)
		Ihat[j,,x[[1]],x[[2]]] <- Yhat*c(mean(prop))
	}
}
 
#----------------------------------------------------------------------------------------------
# STEP 6: estimation with design-based estimators
#----------------------------------------------------------------------------------------------
 
# objects
Ibar <- array(NA, dim=c(ns,3,nsim,3), dimnames=list(samp=1:ns, page=1:3, nsim=1:nsim, mucor=c("indep", "dep045","dep090")))
# estimation
set.seed(222)
for(j in 1:ns){
cat("\n",j ,":", date(),"-", sep="")
	samp <- sample(1:n, 100)
	for(i in 1:nrow(xd)){
		x <- xd[i,]
		cat(".", sep="")
		# building the sample
		Isamp <- Isim.ln[samp,,x[[1]],x[[2]]]
		Ysamp <- apply(Isamp,1,sum)
		# store means
		Yres[j,c("Ybar", "Ybarvar"),x[[1]]] <- c(mean(Ysamp), var(Ysamp))	
		# estimation of abundance-at-age
		Ibar[j,,x[[1]],x[[2]]] <- apply(Isamp,2,mean)
	}
}
 
#==============================================================================================
# RESULTS
#==============================================================================================
 
#----------------------------------------------------------------------------------------------
# STEP 7: summary statistics
#----------------------------------------------------------------------------------------------
 
# objects
sim.res <- array(NA, dim=c(nsim,3,3,2,3), dimnames=list(iter=1:nsim, age=1:3, mucor=c("indep", "dep045","dep090"), meth=c("geo","samp"), stats=c("bias", "mse", "ICcov")))
# computation
for(i in 1:nsim){
	lmusim <- Isim.lnhat[,i,]
	# geo 
	m <- Ihat[,,i,]
	q025 <- apply(m,c(2,3), quantile, probs=0.025, na.rm=T)
	q975 <- apply(m,c(2,3), quantile, probs=0.975, na.rm=T)
	sim.res[i,,,"geo","ICcov"] <- q025<=lmusim & q975>=lmusim
	for(j in 1:ns){
		m[j,,] <- m[j,,]-lmusim 
	}
	sim.res[i,,,"geo","bias"] <- apply(m, c(2,3), mean)
	sim.res[i,,,"geo","mse"] <- apply(m, c(2,3), var)
 
	# samp
	m <- Ibar[,,i,]
	q025 <- apply(m,c(2,3), quantile, probs=0.025, na.rm=T)
	q975 <- apply(m,c(2,3), quantile, probs=0.975, na.rm=T)
	sim.res[i,,,"samp","ICcov"] <- q025<=lmusim & q975>=lmusim
	for(j in 1:ns){
		m[j,,] <- m[j,,]-lmusim 
	}
	sim.res[i,,,"samp","bias"] <- apply(m, c(2,3), mean)
	sim.res[i,,,"samp","mse"] <- apply(m, c(2,3), var)
}
# table
tab.res <-  apply(sim.res, c(2,3,4,5), mean)
 
#################################################################################################
# THE END (or the beginning ...) 
#################################################################################################