##################################### ## Carregando pacotes necessários ##################################### require(sp) ## pacote que define classes para representação de dados espaciais no R require(rgdal) ## importação/exportação de mapas, imagens e objetos geográficos em geral require(spdep) ## funções para análises de dados de áreas incluindo econometria espacial require(classInt) ## rotinas para faclitar a divisão de dados em classes por vários critérios require(RColorBrewer) ## usada aqui para criar palhetas de cores nas visualizações em mapas ### guarda parametros originais da janela grafica par.ori <- par(no.readonly=TRUE) ## Carregando os dados (atributos e geometria de polygonos) disponíveis no pacote spdep ## lendo dados tipo shapefiles loc <- system.file("etc/shapes", package="spdep")[1] loc auckland <- readOGR(loc, 'auckland') ## entendendo a classe/tipo de dado plot(auckland) class(auckland) slotNames(auckland) ## visualizando os atributos head(auckland@data) ## alguns resumos dos atributos summary(auckland) summary(auckland$M77_85) summary(auckland$Und5_81) ## taxas anuais de mortalidade summary(auckland$M77_85/(9*auckland$Und5_81)) ## idem por 1000 habitantes summary(1000*auckland$M77_85/(9*auckland$Und5_81)) taxa1000 <- 1000*auckland$M77_85/(9*auckland$Und5_81) hist(taxa1000, freq=F) rug(taxa1000) lines(density(taxa1000)) ## calculos de taxas etc res <- probmap(auckland$M77_85, 9*auckland$Und5_81) str(res) head(res, n=10) (auckland$M77_85/(9*auckland$Und5_81))[1:10] options(width=100, digits=2) head(cbind(auckland$M77_85, 9*auckland$Und5_81, res), n=10) ## dividindo dados em classes (usando pacote classInt) ## veja ?classIntervals kmeans7_raw <- classIntervals(res$raw*1000, n=7, style="kmeans") kmeans7_raw ############################################################### ## Preparando visualização gráfica das taxas brutas por classes ############################################################### ## arredondando limites das classes brks_raw_rate <- round(kmeans7_raw$brks, digits=2) brks_raw_rate ## recalculando as frequencias com os valores arredondados das classes rkmeans7_raw <- classIntervals(res$raw*1000, style="fixed", fixedBreaks=brks_raw_rate) rkmeans7_raw ## taxa global sum(auckland$M77_85)/sum(9*auckland$Und5_81)*1000 ## criando uma palheta de cores (usando pacote RColorBrewer) de azul para vermelho pal_raw <- c(rev(brewer.pal(3, "Blues")), brewer.pal(4, "Reds")) ## vendo as cores criadas pal_raw plot(1:7, 1:7, pch=19, cex=10, col=pal_raw, xlim=c(0, 8), ylim=c(0, 8)) ## atribuindo as cores adequadas a cada grupo cols <- findColours(rkmeans7_raw, pal_raw) plot(auckland, col=cols) title(main="Taxas brutas anuais (por 1000)") ## criando uma legenda table <- attr(cols, "table") legtext <- paste(names(table), " (", table, ")", sep="") legend("bottomleft", fill=attr(cols, "palette"), legend=legtext, bty="n", cex=0.9, y.inter=0.8) ################################################### ## idem anterior, agora para riscos relativos ################################################### kmeans7_smr <- classIntervals(res$relRisk, n=7, style="kmeans") brks_smr <- round(kmeans7_smr$brks) rkmeans7_smr <- classIntervals(res$relRisk, style="fixed", fixedBreaks=brks_smr) rkmeans7_smr pal_smr <- pal_raw cols <- findColours(rkmeans7_smr, pal_smr) plot(auckland, col=cols) title(main="Taxas de mortalidade padronizadas (SMR)") table <- attr(cols, "table") legtext <- paste(names(table), " (", table, ")", sep="") legend("bottomleft", fill=attr(cols, "palette"), legend=legtext, bty="n", cex=0.7, y.inter=0.8) ################################################### ### Mapa de probabilidade (Poisson) ################################################### brks_prob <- c(0,0.05,0.1,0.2,0.8,0.9,0.95,1) fixed_prob <- classIntervals(res$pmap, style="fixed", fixedBreaks=brks_prob) pal_prob <- rev(brewer.pal(5, "RdBu")) cols <- findColours(fixed_prob, pal_prob) plot(auckland, col=cols) title(main="Mapa de probabilidades Poisson") table <- attr(cols, "table") legtext <- paste(names(table), " (", table, ")", sep="") legend("bottomleft", fill=attr(cols, "palette"), legend=legtext, bty="n", cex=0.7, y.inter=0.8) ###################################################################################################### ### Outro tipo de mapa de probabilidades indicando distritos "significativamente" atípicos ###################################################################################################### resCh <- choynowski(auckland$M77_85, 9*auckland$Und5_81) resCh ## definindo cores para taxas "normais" e contagens mais altas e baixas que esperado cols <- rep("white", length(resCh$pmap)) cols[(resCh$pmap < 0.05) & (resCh$type)] <- "grey35" cols[(resCh$pmap < 0.05) & (!resCh$type)] <- "grey75" plot(auckland, col=cols) title(main="Mapa de probabilidade de Choynowski") legend("bottomleft", fill=c("grey35", "white", "grey75"), legend=c("low", "N/S", "high"), bty="n")