Não foi possível enviar o arquivo. Será algum problema com as permissões?
Essa é uma revisão anterior do documento!
Simulation study
Test beta bias
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()
O problema está na variância do beta que não é um estimador da variância do processo. Logo quando fazemos a backtrans há uma subestimação da média do processo.
Algorítmo
## 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 (função isim abaixo) for(i in 1:length(spcor)){ Isim[,,,i,] <- 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) }
isim <- function(gs, sig00, niter, tsq00, seed=333){ set.seed(seed) arr <- array(NA, dim=c(n=nrow(gs),3,niter,2), dimnames=list(loc=1:n, age=1:3, niter=1:niter, mucor=c("indep", "dep"))) m.arr <- array(NA, dim=c(n=nrow(gs),3,niter), dimnames=list(loc=1:n, age=1:3, niter=1:niter)) # parâmetros do componente comum (S00) phi00 <- sig00/5 # parâmetros dos componentes individuais sig11 <- sig22 <- sig33 <- 1-sig00 phi11 <- phi22 <- phi33 <- (1-sig00)/5 mu1 <- mu2 <- mu3 <- 0.4 ## simulando S S00 <- grf(grid=gs, cov.pars=c(sig00/2, phi00), nsim=niter) S11 <- grf(grid=gs, cov.pars=c(sig11/2, phi11), nsim=niter) S22 <- grf(grid=gs, cov.pars=c(sig22/2, phi22), nsim=niter) S33 <- grf(grid=gs, cov.pars=c(sig33/2, phi33), nsim=niter) ## simulando mu # independent set.seed(111) Cor <- diag(c(1,1,1)) for(i in 1:niter){ m.arr[,,i] <- mvrnorm(n, c(mu1, mu2, mu3), Cor*tsq00) } ## construindo Y = S+e arr[,1,,1] <- m.arr[,1,] + S00$data + S11$data arr[,2,,1] <- m.arr[,2,] + S00$data + S22$data arr[,3,,1] <- m.arr[,3,] + S00$data + S33$data # dependent set.seed(111) Cor[2,1] <- -0.9 Cor[1,2] <- -0.9 for(i in 1:niter){ m.arr[,,i] <- mvrnorm(n, c(mu1, mu2, mu3), Cor*tsq00) } ## construindo Y = S+e arr[,1,,2] <- m.arr[,1,] + S00$data + S11$data arr[,2,,2] <- m.arr[,2,] + S00$data + S22$data arr[,3,,2] <- m.arr[,3,] + S00$data + S33$data # out arr }
resultados
</code R> Ihat.bias ← apply(Ihat - exp(1.2+log(0.33)+0.7/2), 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)
Ihat.bias
, , mucor = indep
spcor
age 0 0.25 0.5 0.75 1
1 -0.2310553 -0.2861606 -0.3006835 -0.1239977 0.3697025 2 -0.2123689 -0.2996022 -0.2987973 -0.1059902 0.3663539 3 -0.2199671 -0.2902498 -0.2877976 -0.1249653 0.3581888
, , mucor = dep
spcor
age 0 0.25 0.5 0.75 1
1 -0.1963665 -0.2576779 -0.2853226 -0.12670738 0.3707011 2 -0.1859164 -0.2815529 -0.2817175 -0.09679016 0.3764652 3 -0.2268756 -0.3127838 -0.3202839 -0.17392410 0.2859450
Ibar.bias
, , mucor = indep
spcor
age 0 0.25 0.5 0.75 1
1 1.347212 3.053397 4.841803 6.413725 10.24748 2 1.460882 2.975048 4.877596 6.232822 10.83875 3 1.507800 2.930321 4.902475 6.425564 10.72584
, , mucor = dep
spcor
age 0 0.25 0.5 0.75 1
1 0.8780813 2.174725 3.593623 5.212714 8.202600 2 0.8637434 2.198764 3.656170 4.997954 8.087312 3 1.2216741 2.651468 4.023795 5.665337 9.817139
</code>