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}<d_{6,10}$, 
embora a taxa brutra do município 12 seja 1.03 por mil e 
do município 10 seja 0,91 por mil, esta sendo mais próxima 
da taxa bruta do município 6, 0,92 por mil. 
Esse caso é curioso porque observamos aqui que a medida 
de $d_{i,j}$ não é simétrica. Podemos observar também 
que o município 12 tem população quase o dobro maior que 
a população dos municípios 6 e 10, e ocorreu o dobro do 
número de casos, 8 casos. 
Neste caso, a variância da taxa bruta do município 12 
é menor que dos municípios 6 e 10 e esta pode ser uma 
justificativa plausível para o fato de que $d_{6,12}<d_{6,10}$.

A seguir, nós redefinimos a função do cálculo de 
distância para que retorne apenas $d_{i,j}$, 
calculamos a lista de custos e criamos o 
objeto da classe \code{listw} que contém a 
lista de vizinhança e a lista de custos. 
\begin{Schunk}
\begin{Sinput}
> 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<q)=1-\alpha$ com $\alpha$ o nível de 
significância e $P(\chi)$ a distribuição qui-quadrado com 
um grau de liberdade. 
A distribuição qui-quadrado justifica-se pelo fato de que 
sob-hipótese núla a diferença da soma sequencial de 
homogeneidades é uma diferença de verossimihanças, 
com diferença de um parâmetro. 
\begin{Schunk}
\begin{Sinput}
> 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}
