Simulation study

library(geoR)
library(lattice)
library(MASS)
 
## definindo grid de simulação
gs <- expand.grid((0:50)/50, (0:50)/50)
n <- nrow(gs)
 
## 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# 
 
nsim <- 100
set.seed(475)
 
## parâmetros do componente comum (S00) 
sig00 <- 0.9
phi00 <- sig00
 
## parâmetros dos componentes individuais 
sig11 <- sig22 <- sig33 <- 1-sig00
phi11 <- phi22 <- phi33 <- 1-phi00
mu1 <- mu2 <- mu3 <- 1
 
## simulando S
S00 <- grf(grid=gs, cov.pars=c(sig00, phi00), nsim=nsim)
S11 <- grf(grid=gs, cov.pars=c(sig11, phi11), nsim=nsim)
S22 <- grf(grid=gs, cov.pars=c(sig22, phi22), nsim=nsim)
S33 <- grf(grid=gs, cov.pars=c(sig33, phi33), nsim=nsim)
 
## simulando mu
## ruídos
tsq0 <- tsq1 <- tsq2 <- tsq3 <- 1
 
# independente
set.seed(111)
Cor <- diag(c(1,1,1))
m1 <- mvrnorm(n, c(mu1, mu2, mu3), Cor*tsq0)
# dependent
Cor[2,1] <- -0.75
Cor[1,2] <- -0.75
m2 <- mvrnorm(n, c(mu1, mu2, mu3), Cor*tsq0)
 
m <- m1
 
## construindo Y = log(S)+e
Y <- Y1 <- Y2 <- Y3 <- S00
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
 
## use plots to check dependency
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))))))
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))
 
df0.mean <- tapply(df0$prop, list(age=df0$age, Yfac=df0$Yfac, samp=df0$samp), mean)
df0.mean <- data.frame(expand.grid(dimnames(df0.mean)), prop=c(df0.mean))
bwplot(prop~age|Yfac, data=df0.mean, auto.key=list(space="right"))