Diferenças

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

--- artigos:ernesto3:sim [2008/09/23 10:18] – criada ernesto
+++ artigos:ernesto3:sim [2008/10/13 13:00] (atual) – ernesto
@@ Linha 1: / Linha 1: @@
 ======Simulation study ======
+===== Algorithm =====
 <code R>
+#################################################################################################
+#
+#	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)
-## definindo grid de simulação
+#==============================================================================================
-gs <- expand.grid((0:50)/50, (0:50)/50)
+# SIMULATING THE PROCESSES
+#==============================================================================================
+# grid
+gs <- expand.grid((0:40)/40, (0:40)/40)
+# size of processes
 n <- nrow(gs)
+# number of replicates
+nsim <- 250
-## Modelo: variáveis com mu dependente e correlação componente espacial parcialmente comum
+#----------------------------------------------------------------------------------------------
-##	Y1 = mu1 + S00 + S11 + e1 = mu1 + sig00 R + sig11 R11 + e1
+# STEP 1 : gaussian processes to build abundance and compositions
-##	Y2 = mu2 + S00 + S22 + e2 = mu1 + sig00 R + sig22 R22 + e2
+#----------------------------------------------------------------------------------------------
-##	Y3 = mu3 + S00 + S33 + e3 = mu1 + sig00 R + sig33 R33 + e3
-##	Constraint by:
-##		sig00 + sig## = 1; e# = 1
-##		mu ~ MVG(c(mu1, mu2, mu3), Sigma)
-##		mu1 = mu2 = mu3 = 1
-##		diag(Sigma) = e#
-nsim <- 100
+# object
-set.seed(475)
+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
+#----------------------------------------------------------------------------------------------
-## parâmetros do componente comum (S00)
+# Generate a spatial Gaussian process Z
-sig00 <- 0.9
+phi <- 0.2
-phi00 <- sig00
+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)
-## parâmetros dos componentes individuais
+#----------------------------------------------------------------------------------------------
-sig11 <- sig22 <- sig33 <- 1-sig00
+# STEP 3: build compositions
-phi11 <- phi22 <- phi33 <- 1-phi00
+#----------------------------------------------------------------------------------------------
-mu1 <- mu2 <- mu3 <- 1
-## simulando S
+# objects
-S00 <- grf(grid=gs, cov.pars=c(sig00, phi00), nsim=nsim)
+arr00 <- array(1, dim=c(n,3,nsim), dimnames=list(loc=1:n, age=1:3, nsim=1:nsim))
-S11 <- grf(grid=gs, cov.pars=c(sig11, phi11), nsim=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")))
-S22 <- grf(grid=gs, cov.pars=c(sig22, phi22), nsim=nsim)
+# option 1: no association between compositions and age aggregated abundance
-S33 <- grf(grid=gs, cov.pars=c(sig33, phi33), nsim=nsim)
+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
-## simulando mu
+#----------------------------------------------------------------------------------------------
-## ruídos
+# STEP 4: build abundance at age
-tsq0 <- tsq1 <- tsq2 <- tsq3 <- 1
+#----------------------------------------------------------------------------------------------
-# independente
+# objects
-set.seed(111)
+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")))
-Cor <- diag(c(1,1,1))
+# sim
-m1 <- mvrnorm(n, c(mu1, mu2, mu3), Cor*tsq0)
+for(i in 1:3){
-# dependent
+	for(j in 1:3){
-Cor[2,1] <- -0.75
+		Isim.ln[,i,,j] <- Psim[,i,,j]*Ysim$data
-Cor[1,2] <- -0.75
+	}
-m2 <- mvrnorm(n, c(mu1, mu2, mu3), Cor*tsq0)
+}
+# 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
+#----------------------------------------------------------------------------------------------
-m <- m1
+# 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)
+	}
+}
-## construindo Y = log(S)+e
+#==============================================================================================
-Y <- Y1 <- Y2 <- Y3 <- S00
+# RESULTS
-Y1$data <- exp(m[,1] + S00$data + S11$data)
+#==============================================================================================
-Y2$data <- exp(m[,2] + S00$data + S22$data)
-Y3$data <- exp(m[,3] + S00$data + S33$data)
-# abundância total
+#----------------------------------------------------------------------------------------------
-Y$data <- Y1$data + Y2$data + Y3$data
+# STEP 7: summary statistics
+#----------------------------------------------------------------------------------------------
-## use plots to check dependency
+# objects
-Yfac <- apply(Y$data, 2, function(x) as.numeric(cut(x, c(-1,quantile(x, probs=c(0.2,0.4,0.6,0.8,1))))))
+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")))
-df0 <- data.frame(expand.grid(list(loc=1:n, samp=1:100, age=1:3)), I=c(Y1$data, Y2$data, Y3$data), prop=c(Y1$data, Y2$data, Y3$data)/c(Y$data), Yfac=factor(Yfac))
+# 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)
-df0.mean <- tapply(df0$prop, list(age=df0$age, Yfac=df0$Yfac, samp=df0$samp), mean)
+	# samp
-df0.mean <- data.frame(expand.grid(dimnames(df0.mean)), prop=c(df0.mean))
+	m <- Ibar[,,i,]
-bwplot(prop~age|Yfac, data=df0.mean, auto.key=list(space="right"))
+	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 ...)
+#################################################################################################
 </code>
+===== Results =====

Diferenças

Navegação

Busca

Ferramentas

QR Code