Nesta seção mostramos como utilizar a função \code{skater()} para fazer a análise de agrupamentos quando tempos uma variável de contagem. Nós vamos utilizar neste exemplo, a distância função de distância DJ definida. Incialmente vamos carregar o mapa da meso-região de Belo Horizonte: \begin{Schunk} \begin{Sinput} > require(spdep) > bh <- readShapePoly(system.file("etc/shapes/bhicv.shp", package = "spdep")[1]) > (n <- nrow(coordinates(bh))) \end{Sinput} \begin{Soutput} [1] 98 \end{Soutput} \end{Schunk} O objeto \code{bh} contém o mapa das 98 cidades que fazem parte da meso região de Belo Horizonte e 7 variáveis, uma delas sendo a população. Nós vamos dividir o mapa em cinco regiões, mostradas na figura~\ref{fig:bh5} e atribuir uma taxa diferente para cada grupo. Nós simulamos um conjunto de dados a taxa definida para fazer a análise de agrupamentos. Na figura~\ref{fig:bh5} observamos a SMR dos dados simulados. \setkeys{Gin}{width=0.99\textwidth} \begin{figure} \centering \begin{Schunk} \begin{Sinput} > grs <- list() > grs$g1 <- which(coordinates(bh)[, 2] > -19) > grs$g2 <- which(coordinates(bh)[, 2] < (-20.2)) > grs$g3 <- setdiff(which(coordinates(bh)[, 1] > -43.5), unlist(grs)) > grs$g4 <- setdiff(which(coordinates(bh)[, 1] < (-44.3)), unlist(grs)) > grs$g5 <- setdiff(1:n, unlist(grs)) > gr <- integer(n) > for (i in 1:length(grs)) gr[grs[[i]]] <- i > lambda <- c(0.001, 0.001, 5e-04, 5e-04, 1e-04)[gr] > esp0 <- lambda * bh$Population > set.seed(123) > obs <- rpois(n, esp0) > sum(obs == 0) \end{Sinput} \begin{Soutput} [1] 16 \end{Soutput} \begin{Sinput} > tx.obs <- sum(obs)/sum(bh$Population) > esp.obs <- tx.obs * bh$Population > summary(smr <- obs/esp.obs) \end{Sinput} \begin{Soutput} Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0000 0.4212 2.1660 2.7000 4.5950 9.7570 \end{Soutput} \begin{Sinput} > brea.smr <- c(0, 0.2, 0.5, 2, 5, Inf) > table(cl.smr <- findInterval(smr, brea.smr)) \end{Sinput} \begin{Soutput} 1 2 3 4 5 17 13 15 35 18 \end{Soutput} \begin{Sinput} > par(mfrow = c(1, 2), mar = c(0, 0, 0, 0)) > plot(bh, col = gray(1 - lambda * 100)) > text(-43.5, c(-18.2, -18.1), c("0.5%", "G1")) > text(-43.3, c(-20.9, -20.8), c("0.5%", "G2")) > text(-42.7, c(-19.6, -19.5), c("0.2%", "G3")) > text(-45, c(-19.3, -19.2), c("0.2%", "G4")) > text(-44.4, c(-18.9, -18.8), c("0.05%", "G5")) > plot(bh, col = gray(1 - 1:5/7)[cl.smr]) > legend("topleft", leglabs(brea.smr, "<", ">"), fill = gray(1 - + 1:5/7), bty = "n") \end{Sinput} \end{Schunk} \includegraphics{figs/fig-bh5} \caption{Map of Belo Hozionte macro-region with five groups defined for example (left) and the map of SMR obtained with one simulated data} \label{fig:bh5} \end{figure} Nós precisamos definir uma função que calcule $d_{i,j}$ e utiliza-la no argumento \code{otherfun} da função \code{skater()}. Essa função deve ter três argumentos, um \code{data.frame}, um inteiro indicando a área da na qual a função será avaliada e o terceiro, um vetor de inteiros para indicar as áreas vizinhas. Nós criamos um \code{data.frame} com a população e os casos simulados obtemos lista de vizinhança: \begin{Schunk} \begin{Sinput} > dat <- data.frame(pop = bh$Population, obs = obs) > nb.bh <- poly2nb(bh) \end{Sinput} \end{Schunk} Com isso, temos que a primeira coluna do \code{data.frame} contém a população e a segunda contém o número de casos. Então, nós definimos a função para cálculo de DJ considerando essa estrutura de dados. Nós trabalhamos na escala logaritma para melhor precisão nos resultados. Além disso, nós consideramos duas correções para evitar o efeito de estar-mos trabalhando com taxas pequenas. 1) se a taxa estimada é zero, fazemos com que seja igual a 1e-10. 2) se o número esperado de casos é menor que 1 e o número observado é zero, fazemos com que o valor máximo que a função de densidade de probabilidade pode assumir seja 1. Inicialmente, nós vamos escrever uma função para o cálculo de DJ que retorne os valores intermediários das quantidades envolvidas nos cálculos e aplicar essa função as dados do município 6: \begin{Schunk} \begin{Sinput} > dj0.pois <- function(data, i, j) { + tx <- c(data[i, 2] + data[j, 2])/c(data[i, 1] + data[j, 1]) + tx[tx == 0] <- 1e-10 + ei <- tx * dat[i, 1] + ej <- tx * dat[j, 1] + fi <- dpois(dat[i, 2], ei, log = TRUE) + fj <- dpois(dat[j, 2], ej, log = TRUE) + maxi <- ifelse(ei < 0.5, ei/10, ei - 0.5) * log(ei) - ei - + lgamma(ifelse(ei < 0.5, 1 + ei/10, ei + 0.5)) + maxj <- ifelse(ej < 0.5, ej/10, ej - 0.5) * log(ej) - ej - + lgamma(ifelse(ej < 0.5, 1 + ej/10, ej + 0.5)) + maxi[dat[i, 2] == 0 & ei < 1] <- 0 + maxj[dat[j, 2] == 0 & ej < 1] <- 0 + data.frame(pop = dat[j, 1], obs = dat[j, 2], tx0 = 1000 * + dat[j, 2]/dat[j, 1], tx.j = 1000 * tx, ei, ej, fi = exp(fi), + fj = exp(fj), mi = exp(maxi), mj = exp(maxj), dij = 1 - + exp(fi - maxi + fj - maxj), row.names = paste(j)) + } > data.frame(dat[6, ], tx = 1000 * dat[6, 2]/dat[6, 1]) \end{Sinput} \begin{Soutput} pop obs tx 6 4320 4 0.925926 \end{Soutput} \begin{Sinput} > round(dj0.pois(dat, 6, nb.bh[[6]]), 3) \end{Sinput} \begin{Soutput} pop obs tx0 tx.j ei ej fi fj mi mj dij 7 6797 10 1.471 1.259 5.440 8.560 0.158 0.112 0.172 0.137 0.252 9 5389 10 1.856 1.442 6.229 7.771 0.124 0.093 0.161 0.144 0.501 10 4393 4 0.911 0.918 3.966 4.034 0.195 0.195 0.202 0.201 0.061 12 7760 8 1.031 0.993 4.291 7.709 0.193 0.139 0.194 0.144 0.044 13 6593 3 0.455 0.641 2.771 4.229 0.154 0.184 0.243 0.196 0.407 \end{Soutput} \begin{Sinput} > 1000 * dat[6, 2]/dat[6, 1] - 1000 * dat[nb.bh[[6]], 2]/dat[nb.bh[[6]], + 1] \end{Sinput} \begin{Soutput} [1] -0.54531138 -0.92970592 0.01538643 -0.10500191 0.47089787 \end{Soutput} \end{Schunk} Observamos que o municípios 9 e 13 são os mais distantes do município 6. Esse resultado é razoável pois a diferença entre a taxa bruta observada em ambos e a taxa bruta do município 6 são as maiores entre os 5 vizinhos. Nós observamos que o município 10 é muito parecido com o município 6, pois tem população com tamanho muito próximo e ocorreu o mesmo número de casos. Porém, $d_{6,12} dj.pois <- function(data, i, j) { + tx <- c(data[i, 2] + data[j, 2])/c(data[i, 1] + data[j, 1]) + tx[tx == 0] <- 1e-10 + ei <- tx * dat[i, 1] + ej <- tx * dat[j, 1] + fi <- dpois(dat[i, 2], ei, log = TRUE) + fj <- dpois(dat[j, 2], ej, log = TRUE) + maxi <- ifelse(ei < 0.5, ei/10, ei - 0.5) * log(ei) - ei - + lgamma(ifelse(ei < 0.5, 1 + ei/10, ei + 0.5)) + maxj <- ifelse(ej < 0.5, ej/10, ej - 0.5) * log(ej) - ej - + lgamma(ifelse(ej < 0.5, 1 + ej/10, ej + 0.5)) + maxi[dat[i, 2] == 0 & ei < 1] <- 0 + maxj[dat[j, 2] == 0 & ej < 1] <- 0 + 1 - exp(fi - maxi + fj - maxj) + } > nb.dj <- nbcosts(nb.bh, dat, meth = "oth", oth = dj.pois) > nb.wj <- nb2listw(nb.bh, nb.dj) \end{Sinput} \end{Schunk} A partir do objeto \code{listw} nós obtemos a árvore geradora mínima: \begin{Schunk} \begin{Sinput} > mstj.bh <- mstree(nb.wj) \end{Sinput} \end{Schunk} \includegraphics{figs/fig-mst} Agora, nós podamos alguns galhos da AGM utilizando o algoritmo SKATER implementado na função \code{skater}. Para isso, precisamos definir a função para o cálculo da medida de homogeneidade. Essa função deve possuir dois argumentos, o primeiro deve ser um \code{data.frame} e o segundo um vetor inteiro com os índices dos elementos do grupo. Nos definimos uma função para o cálculo da verossimilhança avaliada no ponto da estimativa de máxima verossimilhança da taxa estimada com os dados do grupo. \begin{Schunk} \begin{Sinput} > nllpois <- function(data, id) { + tx <- sum(data[id, 2])/sum(data[id, 1]) + -sum(dpois(data[id, 2], data[id, 1] * tx, log = TRUE)) + } \end{Sinput} \end{Schunk} Inicialmente vamos obter 2 grupos, fazendo uma poda na árvore. Atribuimos o resultado a um objeto e em seguida obtemos 3, 4 e 5 grupos sequencialmente, atribuindo cada resultado a um objeto diferente. Ao final, fazemos mais 5 podas para gerar 10 grupos, com objetivo de estudar o decaimento da função de homogeneidade. \begin{Schunk} \begin{Sinput} > sk2 <- skater(mstj.bh[, 1:2], dat, 1, method = "other", other = nllpois) > sk3 <- skater(sk2, dat, 1, method = "other", other = nllpois) > sk4 <- skater(sk3, dat, 1, method = "other", other = nllpois) > sk5 <- skater(sk4, dat, 1, method = "other", other = nllpois) > sk10 <- skater(sk5, dat, 5, method = "other", other = nllpois) \end{Sinput} \end{Schunk} Podemos observar o decaimento da soma de verossimilhanças e utilizar a diferença como critério de parada. Podemos decidir parar quando essa diferença for menor que $q$ tal que $P(\chi sk10$ssw \end{Sinput} \begin{Soutput} [1] 981.6896 563.4223 320.9714 224.4314 207.8563 201.6048 198.3827 196.0463 [9] 193.9180 191.7909 \end{Soutput} \begin{Sinput} > diff(sk10$ssw) \end{Sinput} \begin{Soutput} [1] -418.267315 -242.450828 -96.540053 -16.575047 -6.251574 -3.222115 [7] -2.336325 -2.128342 -2.127078 \end{Soutput} \end{Schunk} Observamos que para esse conjunto de dados simulados, o cinco é o número adequado de grupos, se consideramos $\alpha=0.05$ e $q=3.84$ Na figura~\ref{fig:poisgrupos} temos os mapas dos grupos e as sub-arvores obtidas a cada passo, desde dois grupos até cinco grupos. No gráfico superior esquerdo, a aresta amarela é o a aresta cortada no primeiro passo do algoritmo e as arestas vermelhas e azuis definem os dois grupos gerados. Olhando para cada gráfico, podemos seguir os passos do SKATER. O resultado obtido para cinco grupos, gráfico inferior direito, é bastante satisfatório, pois um grupo foi perfeitamente identificado e os demais não são muito diferentes dos verdadeiros grupos. \begin{figure} \centering \begin{Schunk} \begin{Sinput} > par(mfrow = c(2, 2), mar = c(0, 0, 0, 0)) > plot(bh, col = gray(1 - sqrt(lambda) * 10)) > plot(mstj.bh, coordinates(bh), label.areas = "", cex.cir = 0.01, + cex.lab = 0.001, add = T, col = "yellow") > plot(sk2, coordinates(bh), cex.cir = 0.01, cex.lab = 0.001, add = T, + lwd = 2) > plot(bh, col = gray(1 - sqrt(lambda) * 10)) > plot(sk3, coordinates(bh), cex.cir = 0.01, cex.lab = 0.001, add = T, + lwd = 2) > plot(bh, col = gray(1 - sqrt(lambda) * 10)) > plot(sk4, coordinates(bh), cex.cir = 0.01, cex.lab = 0.001, add = T, + lwd = 2) > plot(bh, col = gray(1 - sqrt(lambda) * 10)) > plot(sk5, coordinates(bh), cex.cir = 0.01, cex.lab = 0.001, add = T, + lwd = 2) \end{Sinput} \end{Schunk} \label{fig:poisgrupos} \caption{The sub-trees obtained by SKATER in results for two groups (top left), three groups (top right), four groups (bottom left) and five groups (bottom right).} \end{figure}