### importa o shapefile library(rgdal) map <- readOGR("~/maps/", "41mu2500gc") ### importa os dados idhm <- read.csv2("dados/idhmPR.csv") dim(idhm) head(idhm,2) ### primeiras linhas dos dados ### primeiras linhas dos atributos do mapa head(map@data,2) ### ordena dados de acordo com ordem do mapa idhm <- idhm[pmatch(map$GEOCODIGO, idhm$Código), ] table(idhm$Código==map$GEOCODIGO) ### confere ## visualiza no mapa plot(map, col=gray(idhm$IDHM.2010)) q5 <- quantile(idhm$IDHM.2010, 0:5/5) legend('topright', format(q5, dig=2), fill=gray(q5)) ### visualização alternativa map$IDHM <- idhm$IDHM.2010 spplot(map, 'IDHM') ### lista de vizinhança library(spdep) viz <- poly2nb(map) vizw <- nb2listw(viz) ### teste de dependência espacial usual: Indice de Moran moran.test(map$IDHM, vizw) ### pacote que contem funções de auto-regressão (mesmas do spdep) library(spatialreg) hist(idhm$IDHM.2010) m1 <- errorsarlm(IDHM.2010 ~ 1, idhm, vizw) m2 <- lagsarlm(IDHM.2010 ~ 1, idhm, vizw) ### igual anterior no caso sem covariaveis (reparametrizacao da media) m3 <- sacsarlm(IDHM.2010 ~ 1, idhm, vizw, vizw) ### não faz sentido sem covariaveis c(AIC(m1), AIC(m2), AIC(m3)) m1 m2 m3 ### um positivo e outro negativo (redundância, não faz sentido sem covariáveis!!!) names(idhm) idhm$probs <- idhm$'Probabilidade.de.sobrevivência.até.60.anos.2010' idhm$rpc <- (idhm$Renda.per.capita.2010/1000) idhm$rpcc <- (idhm$Renda.per.capita.2010/1000-mean(idhm$Renda.per.capita.2010/1000)) summary(idhm$probs) map$probs <- idhm$probs map$rpc <- idhm$Renda.per.capita.2010/1000 library(gridExtra) grid.arrange(spplot(map, c('rpc')), spplot(map, c('probs'))) ### Curitiba se destaca com valor mais alto de Renda! names(idhm) head(idhm,2) ### Define variável indicadora para Curitiba idhm$cwb <- (idhm$Espacialidade=='Curitiba')+0 ### Diagrama de dispersão: Prob. Sobrevivência VS Renda Per Capita plot(idhm$rpc, idhm$probs, col=idhm$cwb+1, pch=4) plot(idhm$rpc, idhm$probs, log='xy', ### escala log em ambos col=idhm$cwb+1, pch=4) ### Curitiba não se destaca na probabilidade de ### sobrevivência até 60 anos, mas na Renda Per Capita ### desconto no efeito de renda em -3.11827 para Curitiba m00 <- lm(probs ~ rpcc, idhm) m0 <- lm(probs ~ cwb + rpcc, idhm) anova(m00, m0) ### não precisa CWB para explicar Sobrevida (pois o destaque é na RPC, a covariável) coef(summary(m0)) ### Modelo com autocorrelação nos resíduos m1s0 <- errorsarlm(probs ~ rpcc, idhm, vizw) plot(predict(m1s0), idhm$probs, col=idhm$cwb+1, pch=4) ### nota-se que o valor de predito de Prob. Sobr. é afetado pelo valor discrepande da covariável (O outilier é na covariável, não na resposta!) ### Inclusão de Curitiba para corrigir o efeito da renda m1s <- errorsarlm(probs ~ cwb + rpcc, idhm, vizw) anova(m1s, m1s0) ### testando as outras formas de incluir autocorrelação ### autocorrelação na resposta m2s <- lagsarlm(probs ~ cwb + rpcc, idhm, vizw) ### Autocorrelação nos resíduos e na resposta m3s <- sacsarlm(probs ~ cwb + rpcc, idhm, vizw, vizw) c(AIC(m1s), AIC(m2s), AIC(m3s)) ### sumário dos modelos summary(m1s) summary(m2s) summary(m3s) ### componentes do modelo nos dois casos de autocorrelação pred1 <- as.data.frame(predict(m1s)) ### erros correlacionados pred2 <- as.data.frame(predict(m2s)) ### resposta autocorrelacionada head(cbind(true=idhm$probs, pred1), 5) head(cbind(true=idhm$probs, pred2), 5) plot(pred1$fit, idhm$probs, pch=4, asp=1); abline(0:1) points(pred2$fit, idhm$probs, pch=3, col=2) ### trend: efeito de covariável(eis) ### signal: efeito espacial suave (do erro ou da autoregressão) ### fit: soma de trend e signal ### visualiza o signal (termo espacial) de cada modelo map$signal1 <- pred1$signal map$signal2 <- pred2$signal grid.arrange(spplot(map, c('signal1')), spplot(map, c('signal2'))) ### Modelo 2: Checando manualmente ### probs = rho * W * probs + beta_0 + beta_1*cwb + beta_2*rpcc + erro (est <- coef(m2s)) mysignal2 <- est[1]*drop(nb2mat(viz)%*%idhm$probs) ## rho * W * probs mytrend2 <- est[2] + est[3]*idhm$cwb + est[4]*idhm$rpcc ## b0 + b1*cwb + b2*rpcc myfit2 <- mysignal2 + mytrend2 cor(cbind(myfit2, mytrend2, mysignal2), pred2) ## compara com retornado pela função predict