Simulation study

gs <- expand.grid((0:10)/10, (0:10)/10)
niter <- 100
arr <- array(NA, dim=c(niter, 2, 2), dimnames=list(niter=1:niter, stat=c("beta","beta.var"), dep=c(1,0)))
 
# gau
set.seed(111)
gd <- grf(grid=gs, cov.pars=c(1,0.2), mean=2, nsim=niter)
 
for(i in 1:niter){
cat(i)
	ggdd <- gd
	ggdd$data <- gd$data[,i]
	lf <- likfit(ggdd, ini.cov.pars=c(1,0.2), messages=FALSE)
	arr[i,,1] <- c(lf$beta, lf$beta.var)
}
 
# gau independente
set.seed(111)
gd <- grf(grid=gs, cov.pars=c(0,0), mean=2, nsim=niter, nugget=1)
 
for(i in 1:niter){
cat(i)
	ggdd <- gd
	ggdd$data <- gd$data[,i]
	lf <- likfit(ggdd, ini.cov.pars=c(0.5,0.5), messages=F)
	arr[i,,2] <- c(lf$beta, lf$beta.var)
}
 
pdf("betabias.pdf")
par(mfrow=c(2,2))
hist(arr[,1,1], main="beta for mu=2 s2=1 phi=0.2 t2=0")
hist(arr[,1,2], main="beta for mu=2 s2=0 phi=0 t2=1")
hist(arr[,2,1], main="beta.var")
hist(arr[,2,2], main="beta.var")
dev.off()

## Modelo: variáveis com mu dependente e correlação componente espacial parcialmente comum
##	Y1 = mu1 + S00 + S11 + e1 = mu1 + sig00 R + sig11 R11 + e1
##	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# 

i indexa localizações j indexa idades

• Definir grid
• Definir sig00 (usado para gerir associação espacial entre idades)
• Simular Ij gaussiana
  • set sig00 e phi00 = sig00/5
  • set sig## = 1-sig00 e phi## = sig##/5
  • set mu## = 0.4
  • sim S ~ MVN(0,Sig) Sig=f(sig, phi)
  • sim mu ~ MVN({mu1,mu2,m3}, Sigmu) Sigmu = {indep, dep}
  • set Ij = muj + S00 + Sjj
• Construir Y = prod(exp{Ij})
• Construir Pj = exp{Ij}/sum_j(exp{Ij})
• Construir Ij = Y*Pj para garantir que Y = sum_j(Ij)
• Aleatorizar vector de localizações (100)
• Construir amostras de Y, P e I
• Ajustar modelo geo
• Estimar Y
  • Yhat = exp{beta+beta.var/2}
  • Ybar = sum_i(Yi)/n
• Estimar Ij
  • Ihatj = Yhat * 1/n sum_i(Pij) sendo Pij = Iij/sum_j(Iij)
  • Ibarj = sum_i(Iij)/n

## INI
# libraries
library(geoR)
library(lattice)
library(MASS)
library(RandomFields)
source("funs.R")
 
# definindo grid de simulação e outros parametros da sim
gs <- expand.grid((0:40)/40, (0:40)/40)
n <- nrow(gs)
niter <- 100
spcor <- seq(0,1,len=5)
 
# objectos
Isim <- array(NA, dim=c(n,3,niter,5,2), dimnames=list(loc=1:n, age=1:3, niter=1:niter, spcor=spcor, mucor=c("indep", "dep")))
Isim.ln <- array(NA, dim=c(n,3,niter,5,2), dimnames=list(loc=1:n, age=1:3, niter=1:niter, spcor=spcor, mucor=c("indep", "dep")))
Psim <- array(NA, dim=c(n,3,niter,5,2), dimnames=list(loc=1:n, age=1:3, niter=1:niter, spcor=spcor, mucor=c("indep", "dep")))
Ysim <- array(NA, dim=c(n,1,niter,5,2), dimnames=list(loc=1:n, age="all", niter=1:niter, spcor=spcor, mucor=c("indep", "dep")))
Yres <- array(NA, dim=c(2,niter,5,2), dimnames=list(stat=c("Ybar","Yhat"), niter=1:niter, spcor=spcor, mucor=c("indep", "dep")))
Isim.lnhat <- array(NA, dim=c(3,niter,5,2), dimnames=list(age=1:3, niter=1:niter, spcor=spcor, mucor=c("indep", "dep")))
Ibar <- array(NA, dim=c(3,niter,5,2), dimnames=list(age=1:3, niter=1:niter, spcor=spcor, mucor=c("indep", "dep")))
Ihat <- array(NA, dim=c(3,niter,5,2), dimnames=list(age=1:3, niter=1:niter, spcor=spcor, mucor=c("indep", "dep")))
 
## SIM
for(i in 1:length(spcor)){
	Isim[,,,i,] <- [[artigos:ernesto3:sim:isim|isim]](gs, spcor[i], niter, 0.2)
}
 
Isim.mean <- apply(Isim, c(2,3,4,5), mean)
Isim.var <- apply(Isim, c(2,3,4,5), var)
 
# Ysim como produto de lognormais
Ysim[,1,,,] <- apply(exp(Isim), c(1,3,4,5), prod)
# cractarerísticas de Y
Ysim.smean <- apply(Ysim, c(3,4,5), mean)
Ysim.svar <- apply(Ysim, c(3,4,5), var)
Ysim.lmean <- apply(log(Ysim), c(3,4,5), mean)
Ysim.lvar <- apply(log(Ysim), c(3,4,5), var)
# composições
Psim[] <- aperm(apply(exp(Isim), c(1,3,4,5), function(x) x/sum(x)), c(2,1,3,4,5))
# observações das marginais
Isim.ln[,1,,,] <- Ysim*Psim[,1,,,,drop=FALSE]
Isim.ln[,2,,,] <- Ysim*Psim[,2,,,,drop=FALSE]
Isim.ln[,3,,,] <- Ysim*Psim[,3,,,,drop=FALSE]
Isim.lnmean <- apply(log(Isim.ln), c(2,3,4,5), mean)
Isim.lnvar <- apply(log(Isim.ln), c(2,3,4,5), var)
Isim.lnhat <- exp(Isim.lnmean+Isim.lnvar/2)
 
## running our model
xd <- expand.grid(dimnames(Ihat)[-1])
samp <- sample(1:n, 100)
 
for(i in 1:nrow(xd)){
	x <- xd[i,]
	cat(as.character(unlist(x)),"\n")
	Isamp <- Isim.ln[samp,,x[[1]],x[[2]],x[[3]]]
	locsamp <- gs[samp,]
	Ysamp <- Ysim[samp,,x[[1]],x[[2]],x[[3]]]
	# geodata
	gd <- as.geodata(cbind(locsamp, Ysamp))	
	lf <- likfit(gd, lambda=0, ini.cov.pars=c(1,0.2), messages=FALSE)
	Yhat <- exp(lf$beta+lf$beta.var/2)
	# store means
	Yres[,x[[1]],x[[2]],x[[3]]] <- c(mean(Ysamp), Yhat)	
	# compositions
	prop <- apply(Isamp,1,function(x) x/sum(x))
	prop[is.na(prop)]<-0
	Ihat[,x[[1]],x[[2]],x[[3]]] <- Yhat*apply(prop,1,mean)
	Ibar[,x[[1]],x[[2]],x[[3]]] <- apply(Isamp,2,mean)
}
 
Ihat.bias <- apply(Ihat - exp(1.2+log(0.33)), c(1,3,4), mean)
Ihat.mse <- apply(Ihat - exp(1.2+log(0.33)+0.7/2), c(1,3,4), var)
Ibar.bias <- apply(Ibar - exp(1.2+log(0.33)+0.7/2), c(1,3,4), mean)
Ibar.mse <- apply(Ibar - exp(1.2+log(0.33)+0.7/2), c(1,3,4), var)