### Objetivo: ajustar modelo de poisson para dados de area ### considerando efeito aleatorio com priori CAR espacial ### como exemplo vamos considera dados de mortalidade infantil ### no estado da Carolina do Norte, EUA ### Esses dados estao no pacote spdep ### carregando o pacote spdep e os dados require(spdep) data(nc.sids) ### objetos carregados ls() ### "ncCC89.nb" "ncCR85.nb" "nc.sids" "sidscents" "sidspolys" ### ver help com ###?nc.sids ### temos a variavel numero de obitos infantis de 1979-1984 ### e a variavel numero de nascimentos em 1979-1984 head(nc.sids) ### vamos criar o numero esperado de obitos infantis oesp <- nc.sids$BIR79 * sum(nc.sids$SID79)/sum(nc.sids$BIR79) ### sumario do observado e esperado (medias sao iguais) summary(nc.sids$SID79) summary(oesp) ### plotando o SMR smr <- nc.sids$SID79/oesp summary(smr) cl.smr <- findInterval(smr, c(0,.7,1,1.5,10)) table(cl.smr) cores <- c("yellow", "orange", "red", "brown") plot(sidspolys, col=cores[cl.smr]) ### para ajustar um modelo poisson com efeito aleatorio com ### priori CAR espacial, vamos utilizar a funcao MCMCmetrop1R() ### do pacote MCMCpack para amostrar da distribuicao a postiori. ### Para fazer isso, temos que escrever a funcao a posteriori em R. ### priori para sigma2 s2.prior <- function(sigma2, a=1, b=0.01) dgamma(1/sigma2, a, b, log=TRUE) ### priori independente para os efeitos aleatorios siid.prior <- function(s, sigma2) sum(dnorm(s, 0, sqrt(sigma2), log=TRUE)) ### verossimilhanca pois.vero <- function(y, lambda) sum(dpois(y, lambda, log=TRUE)) ### posteriori post.piid <- function(pars, y, E) { n <- length(y) lam <- exp(pars[1:n])*E s2.prior(1/exp(pars[n+1])) + siid.prior(pars[1:n], exp(pars[n+1])) + pois.vero(y, lam) } ### carregando pacote MCMCpack require(MCMCpack) ### simulando 1000 amostras a posteriori para calibrar variancia post0 <- MCMCmetrop1R(post.piid, rep(0,101), 1000, 1000, 1, verb=0, V=diag(rep(0.003,101)), y=nc.sids$SID79, E=oesp) dim(post0) ### simulando da posteriori tun <- c(rep(.23,100),.5) vcov <- diag(tun)%*%cov(post0)%*%diag(tun) post <- MCMCmetrop1R(post.piid, colMeans(post0), 50000, 50000, 100, verb=10000, V=vcov, y=nc.sids$SID79, E=oesp) dim(post) par(mfrow=c(10,10), mar=c(0,0,0,0), mgp=c(1.5,.5,0)) for (i in 1:100) plot.ts(post[,i]) par(mar=c(3,3,2,2)) plot(sqrt(exp(post[,101]))) b.m <- colMeans(post[,1:100]) summary(b.m) smrp <- exp(b.m) summary(smrp) cor(smr,smrp) par(mfrow=c(1,1), mar=c(3,3,1,1), mgp=c(2,1,0)) plot(smr,smrp, xlim=range(smr,smrp), ylim=range(smr,smrp)) abline(c(0,1)) resc <- function(x) (x-min(x))/diff(range(x)) par(mfrow=c(2,1), mar=c(1,1,0,0), mgp=c(1,.1,0)) plot(sidspolys, col=gray(1-resc(smr))) plot(sidspolys, col=gray(1-resc(smrp))) cl.smrp <- findInterval(smrp, c(0,.7,1,1.5,10)) table(cl.smr) table(cl.smrp) table(cl.smr,cl.smrp) par(mfrow=c(2,1), mar=c(1,1,1,0), mgp=c(1,.1,0)) plot(sidspolys, col=cores[cl.smr]) plot(sidspolys, col=cores[cl.smrp]) legend("bottomleft", leglabs(c(0,.7,1,1.5,10)), fill=cores, bty="n") ############################################################### ### priori CAR para os efeitos aleatorios scar.prior <- function(s, sigma2, w, nviz) sum(dnorm(s, w%*%s, sqrt(sigma2/nviz), log=TRUE)) ### posteriori post.pcar <- function(pars, y, E, w, nviz) { n <- length(y) lam <- exp(pars[1:n])*E s2.prior(1/exp(pars[n+1])) + scar.prior(pars[1:n], exp(pars[n+1]), w, nviz) + pois.vero(y, lam) } ### definindo a matriz de vizinhanca w <- nb2mat(ncCR85.nb) nviz <- sapply(ncCR85.nb, length) ### simulando 1000 amostras a posteriori para calibrar variancia post0b <- MCMCmetrop1R(post.pcar, rep(0,101), 1000, 1000, 1, verb=0, V=diag(rep(0.001,101)), y=nc.sids$SID79, E=oesp, w=w, nviz=nviz) dim(post0b) ### simulando da posteriori tunb <- c(rep(.23,100),.5) vcovb <- diag(tunb)%*%cov(post0b)%*%diag(tunb) postb <- MCMCmetrop1R(post.pcar, colMeans(post0b), 50000, 50000, 100, verb=10000, V=vcovb, y=nc.sids$SID79, E=oesp, w=w, nviz=nviz) dim(postb) par(mfrow=c(10,10), mar=c(0,0,0,0), mgp=c(1.5,.5,0)) for (i in 1:100) plot.ts(postb[,i]) par(mar=c(3,3,2,2)) plot(sqrt(exp(postb[,101]))) bb.m <- colMeans(postb[,1:100]) summary(bb.m) smrpb <- exp(bb.m) summary(smrpb) cor(cbind(smr,smrp,smrpb)) par(mfrow=c(1,1), mar=c(3,3,1,1), mgp=c(2,1,0)) plot(smr,smrpb, xlim=range(smr,smrpb), ylim=range(smr,smrpb)) abline(c(0,1)) resc <- function(x) (x-min(x))/diff(range(x)) par(mfrow=c(2,1), mar=c(1,1,0,0), mgp=c(1,.1,0)) plot(sidspolys, col=gray(1-resc(smr))) plot(sidspolys, col=gray(1-resc(smrpb))) cl.smrpb <- findInterval(smrpb, c(0,.7,1,1.5,10)) table(cl.smr) table(cl.smrpb) table(cl.smr,cl.smrpb) par(mfrow=c(3,1), mar=c(1,1,1,0), mgp=c(1,.1,0)) plot(sidspolys, col=cores[cl.smr]) plot(sidspolys, col=cores[cl.smrp]) plot(sidspolys, col=cores[cl.smrpb]) legend("bottomleft", leglabs(c(0,.7,1,1.5,10)), fill=cores, bty="n") lam1 <- smrp*oesp lam2 <- smrpb*oesp summary(cbind(oesp,lam1,lam2)) cor(cbind(oesp,lam1,lam2)) cl0 <- findInterval(oesp, c(0,2.5,5,10,70)) cl1 <- findInterval(lam1, c(0,2.5,5,10,70)) cl2 <- findInterval(lam2, c(0,2.5,5,10,70)) table(cl0) table(cl1) table(cl2) par(mfrow=c(3,1), mar=c(1,1,1,0), mgp=c(1,.1,0)) plot(sidspolys, col=cores[cl0]) plot(sidspolys, col=cores[cl1]) plot(sidspolys, col=cores[cl2]) legend("bottomleft", leglabs(c(0,2.5,5,10,70)), fill=cores, bty="n")