library(tidyverse)
library(fpc) #para o calculo do indice de dunn
library(cluster) #para o calculo do indice de silhuetaAnálise de Agrupamento
Introdução
Esta análise tem como objetivo realizar uma análise de agrupamentos para o conjunto de dados apresentados na Tabela 12.4 do livro texto. Esta tabela apresenta para 22 empresas dos EUA os valores das seguintes variáveis:
X1: Índice de cobertura de encargos fixos (renda/dívida).
X2: Taxa de retorno sobre o capital.
X3: Custo por kW de capacidade instalada.
X4: Fator de carga anual.
X5: Crescimento da demanda de pico de kWh de 1974 para 1975.
X6: Vendas (consumo anual de kWh).
X7: Percentual de energia nuclear.
X8: Custo total de combustível (centavos por kWh).
Base de Dados
O código abaixo realiza a leitura da base, verifica que não há na, retira a coluna de identificação da linha e padroniza as variáveis numéricas.
base_original = read_csv("tabela_12_4.csv")
#contando o num de na
sum(is.na(base_original))[1] 0glimpse(base_original)Rows: 22
Columns: 9
$ Company <chr> "Arizona Public Service", "Boston Edison Co.", "Central Louisi…
$ X1 <dbl> 1.06, 0.89, 1.43, 1.02, 1.49, 1.32, 1.22, 1.10, 1.34, 1.12, 0.…
$ X2 <dbl> 9.2, 10.3, 15.4, 11.2, 8.8, 13.5, 12.2, 9.2, 13.0, 12.4, 7.5, …
$ X3 <dbl> 151, 202, 113, 168, 192, 111, 175, 245, 168, 197, 173, 178, 19…
$ X4 <dbl> 54.4, 57.9, 53.0, 56.0, 51.2, 60.0, 67.6, 57.0, 60.4, 53.0, 51…
$ X5 <dbl> 1.6, 2.2, 3.4, 0.3, 1.0, -2.2, 2.2, 3.3, 7.2, 2.7, 6.5, 3.7, 6…
$ X6 <dbl> 9077, 5088, 9212, 6423, 3300, 11127, 7642, 13082, 8406, 6455, …
$ X7 <dbl> 0.0, 25.3, 0.0, 34.3, 15.6, 22.5, 0.0, 0.0, 0.0, 39.2, 0.0, 0.…
$ X8 <dbl> 628.000, 1.555, 1.058, 700.000, 2.044, 1.241, 1.652, 390.000, …#1a coluna e um indicador da empresa e nao deve ser usada no metodo
base = base_original[,-1]
#padronizando a base (somente as var numercias, mas nesse exemplo todas sao numericas)
base = scale(base)Método de Agrupamento Hierárquico
O agrupamento hierárquico é realizado com a função hclust() que recebe como argumento de entrada a matriz de distância dos dados. Esta é obtida no R pela função dist().
dist_matrix = dist(base)Escolha da medida de distância entre grupos
O argumento method da função hclust() determina qual a medida de distância entre os grupos que será utilizada. O padrão é method = "complete", mas existem outras alternativas. Algumas serão verificadas aqui, para mais informações consute o help da função hclust().
cluster_h_complete = hclust(dist_matrix)
cluster_h_centroid = hclust(dist_matrix, method = "centroid")
cluster_h_single = hclust(dist_matrix, method = "single")
cluster_h_ward = hclust(dist_matrix, method = "ward.D2")Escolha do número de grupos
Nesse momento ainda não temos um agrupamento. Precisamos ainda definir o número de grupos.
Pela inspesão visual do dendograma
Uma das formas de escolher o número de grupos é pelo dendograma.
plot(cluster_h_complete, hang = -1)
plot(cluster_h_centroid, hang = -1)
plot(cluster_h_single, hang = -1)
plot(cluster_h_ward, hang = -1)
Os dendogramas para as medidas “centroid” e “single” apresenta agrupamentos tais que a última fusão junta um grupo com um único elemento aos demais. O que não parece um bom agrupamento. Já o dendograma para as medidas “complete” e “ward.D2” não têm essa característica. Vamos seguir somente com esses dois últimos.
No agrupamento hierárquico, o vetor cluster_h$heightguarda as distâncias de fusão entre os clusters a cada etapa do processo. O vetor é crescente: ele começa com as fusões mais próximas e termina com a última fusão (os dois maiores grupos).
Podemos escolher o número de grupos por um critério visual: Corte onde há o maior salto de altura (“maior gap” entre fusões consecutivas). A ideia é encontrar um ponto onde há um grande aumento repentino nas alturas de fusão, sugerindo que unir mais grupos a partir dali resultaria em clusters muito diferentes. O gráfico com a marcação nos locais de fusão dos grupos pode ajudar nessa decisão.
plot(cluster_h_complete, hang = -1)
abline(h=cluster_h_complete$height, lty=2, col="gray")
Analisando o gráfico pensamos em usar 3 ou 4 grupos para o agrupamento hierárquico com medida entre grupos “complete”.
plot(cluster_h_ward, hang = -1)
abline(h=cluster_h_ward$height, lty=2, col="gray")
Analisando o gráfico pensamos em usar 4 grupos para o agrupamento hierárquico com medida entre grupos “ward.D2”.
Índice de Silhueta
Outra forma de escolher o número de grupos é pelo índice de Silhueta. O índice de silhueta mede o quão bem uma observação está agrupada dentro do seu próprio cluster em comparação com outros clusters.
Para cada observação \(i\):
\(a(i)\) representa a distância média entre o ponto \(i\) e todos os demais pontos do mesmo cluster (coerência interna).
\(b(i)\) representa a menor distância média entre o ponto \(i\) e os pontos de qualquer outro cluster (separação externa).
A fórmula da silhueta para a observação \(i\) é:
\[ s(i) = \frac{b(i) - a(i)}{\max\{a(i), b(i)\}} \]
Onde:
- \(s(i)\) varia de \(-1\) a \(1\).
- Valores próximos de \(1\) indicam que a observação está bem alocada em seu cluster.
- Valores próximos de \(0\) indicam sobreposição entre clusters.
- Valores negativos sugerem possível má alocação.
Entre dois agrupamentos vamos preferir aquele que tiver maior valor médio entre os ídnices de silhueta de todas as observações.
Vejamos como encontrar esse valor no R. Vamos comparar o índice de silhueta para um dado agrupamento hierárquicom considerando 3 e 4 grupos.
cluster_h_3 = cutree(cluster_h_complete, 3)
cluster_h_4 = cutree(cluster_h_complete, 4)Agora, para cada um desses agrupamentos, será definido o índice de silhueta médio.
si_3 <- silhouette(cluster_h_3, dist_matrix)
mean(si_3[, 3])[1] 0.1857484si_4 <- silhouette(cluster_h_4, dist_matrix)
mean(si_4[, 3])[1] 0.2164058Podemos fazer isso para vários números de grupos e para os dois agrupamentos que temos até o momento.
Primeiro para o agrupamento considerando a distância “complete”.
num_grupos = 2:10
si = NULL
for(k in num_grupos){
cluster_h_k = cutree(cluster_h_complete, k)
si_k = silhouette(cluster_h_k,dist_matrix)
si = c(si,mean(si_k[, 3]))
}
plot(x = num_grupos, y=si, type="b")
O gráfico sugere realizar o agrupamento com 4 grupos.
Agora para a medida de distância entre grupos “ward.D2”.
num_grupos = 2:10
si = NULL
for(k in num_grupos){
cluster_h_k = cutree(cluster_h_ward, k)
si_k = silhouette(cluster_h_k,dist_matrix)
si = c(si,mean(si_k[, 3]))
}
plot(x = num_grupos, y=si, type="b")
O método da silhueta já sugere 5 grupos e não 4, como foi escolhido pela inspesão visual. Veja pelo dendograma que o agrupamento pela medida “ward.D2” considerando 5 grupos deixa um grupo com um único elemento, o de número 5. Isto não é um agrupamento interessante, um grupo unitário. Vamos então ficar com 4 grupos como sugeriu a isnpesão visual.
Índice de Dunn
O índicer de DUNN pode também ser usado como critério de escolha do número de grupos. Trata-se de uma medida que avalia a qualidade de um agrupamento, considerando simultaneamente a separação entre os grupos e o espalhamento de cada grupo. Quanto mais separados os grupos estiverem entre si, melhor. Quanto menos espalhados os grupos forem, melhor.
O índice de Dunn é calculado como:
\[ \text{Dunn} = \frac{\text{Distância mínima entre clusters}}{\text{Maior dispersão interna (diâmetro) dentro dos clusters}} \]
Buscamos o número de grupos que maximiza o índice de Dunn, garantindo grupos bem separados e com pouca dispersão interna.
D = NULL
for(k in num_grupos){
cluster_h_k = cutree(cluster_h_complete, k)
stats = cluster.stats(dist_matrix,cluster_h_k)
D = c(D,stats$dunn)
}
plot(x=num_grupos,y=D)
D = NULL
for(k in num_grupos){
cluster_h_k = cutree(cluster_h_ward, k)
stats = cluster.stats(dist_matrix,cluster_h_k)
D = c(D,stats$dunn)
}
plot(x=num_grupos,y=D)
O índice de Dunn não apresentou resultados interessantes para a escolha do número de grupos neste exemplo.
Até agora temos dois agrupamentos que parecem interessante: método hierárquico com 4 grupos e medidas de distância entre grupos “complete” e “ward.D2”.
Método K-means
O método K-means é uma técnica de agrupamento que parte de um número pré-definido de grupos (\(k\)), com centroides iniciais geralmente escolhidos de forma aleatória. Como há um elemento de aleatoriedade no processo de inicialização, recomenda-se fixar uma semente aleatória antes da execução do algoritmo. Isso garante que os resultados sejam reproduzíveis e que o agrupamento retorne sempre os mesmos grupos quando o código for reexecutado.
set.seed(123456789)Escolha do número de grupos
A esolha do número de grupos para o k-means será feita pelos índices de Silhueta e Dunn.
Índice de Silhueta
Entre dois agrupamentos vamos preferir aquele que tiver maior valor médio entre os índices de silhueta de todas as observações.
num_grupos = 2:10
si = NULL
for(k in num_grupos){
cluster_k_means = kmeans(base,centers = k)
cluster_k = cluster_k_means$cluster
si_k = silhouette(cluster_k,dist_matrix)
si = c(si,mean(si_k[, 3]))
}
plot(x = num_grupos, y=si, type="b")
O método da silhueta sugeriu a criação de 4 grupos a partir do método k-means.
Índice de Dunn
Entre dois agrupamentos vamos preferir aquele que tiver maior valor médio entre os índices de Dunn.
D = NULL
for(k in num_grupos){
cluster_k_means = kmeans(base, k)
cluster_k = cluster_k_means$cluster
stats = cluster.stats(dist_matrix,cluster_k)
D = c(D,stats$dunn)
}
plot(x=num_grupos,y=D)
Novamente o índice de Dunn não apresentou resultados interessantes para a escolha do número de grupos.
Comparação entre os Agrupamentos
A análise até o momento sugere a criação de 4 grupos pelos três métodos validados: hierárquico com distância “complete”, hierárquico com distância “ward.D2” e k-means.
cluster_h_complete = hclust(dist_matrix)
cluster_h_complete_4 = cutree(cluster_h_complete,k = 4)
cluster_h_ward = hclust(dist_matrix,method = "ward.D2")
cluster_h_ward_4 = cutree(cluster_h_ward,k = 4)
cluster_k_means = kmeans(base,centers = 4)
cluster_k_means_4 = cluster_k_means$clusterUma vez definido os agrupamentos, podemos calcular o índice de silhueta e/ou o índice de Dunn para escolher qual deles será o agrupamento final.
Índice de Silhueta
si_h_complete = silhouette(cluster_h_complete_4,dist_matrix)
mean(si_h_complete[, 3])[1] 0.2164058si_h_ward = silhouette(cluster_h_ward_4,dist_matrix)
mean(si_h_ward[, 3])[1] 0.2270425si_h_kmeans = silhouette(cluster_k_means_4,dist_matrix)
mean(si_h_kmeans[, 3])[1] 0.224333Os valores do índice de silhueta ficaram bem próximos, mas o maior deles foi do agrupamento hierárquico com medida de distância “ward.D2”.
Índice de Dunn
stats_h_complete = cluster.stats(dist_matrix,cluster_h_complete_4)
stats_h_complete$dunn[1] 0.4393243stats_h_ward = cluster.stats(dist_matrix,cluster_h_ward_4)
stats_h_ward$dunn[1] 0.4251652stats_k_means = cluster.stats(dist_matrix,cluster_k_means_4)
stats_k_means$dunn[1] 0.4251652Os índices de Dunn também ficaram bem parecidos uns com os outros, mas o melhor deles foi o método hierárquico com medida de distância “complete”.
Tabela de Contingência
A semelhança entre os índices de Silhueta e índices de Dunn para os três agrupamentos nos leva a crer que os três são agrupamentos semelhantes. Isso pode ser verificado pela tabela de contingência, que para cada par de agrupamentos mostra a semelhança entre eles.
table(
Hierárquico_complete = cluster_h_complete_4,
Hierárquico_ward = cluster_h_ward_4) Hierárquico_ward
Hierárquico_complete 1 2 3 4
1 6 0 1 0
2 0 2 5 0
3 0 5 0 0
4 0 0 0 3Para comparar os dois agrupamentos organize as colunas da tabela de forma que os maiores valores fiquem na diagonal principal e depois veja quantas observação ficaram fora da diagonal principal. Para a comparação entre os dois métodos hierárquicos temos 3 observações agrupadas de forma diferente.
table(
Hierárquico_complete = cluster_h_complete_4,
Hierárquico_kmeans = cluster_k_means_4) Hierárquico_kmeans
Hierárquico_complete 1 2 3 4
1 1 0 0 6
2 6 0 1 0
3 0 0 5 0
4 0 3 0 0Comparando o método hierárquico completo com o k-means, são apenas duas observações que foram agrupadas de forma diferente.
table(
Hierárquico_ward = cluster_h_ward_4,
Hierárquico_kmeans = cluster_k_means_4) Hierárquico_kmeans
Hierárquico_ward 1 2 3 4
1 0 0 0 6
2 1 0 6 0
3 6 0 0 0
4 0 3 0 0A tabela acima nos mostra que a comparação entre o método hierárquico ward e o k-means, apenas uma observação foi agrupada de forma diferente.
Os agrupamentos finais foram muito semelhantes. Vamos seguir para a análise exploratória com o agrupamento hierárquico usando a medida “ward.D2”.
Análise Exploratória
Agora, voltemos para a base original a fim de investigar as diferenças entre os grupos criados. Para facilitar essa análise será incluída uma coluna na base original que define o grupo da observação.
base_original = base_original |> cbind(grupo = cluster_h_ward_4 )Primeiro vejamos as empresas em cada grupo.
base_original |> filter(grupo == 1) |> select(Company) Company
1 Arizona Public Service
2 Central Louisiana Electric Co.
3 Florida Power & Light Co.
4 Oklahoma Gas & Electric Co.
5 The Southern Co.
6 Texas Utilities Co.base_original |> filter(grupo == 2) |> select(Company) Company
1 Boston Edison Co.
2 Hawaiian Electric Co.
3 New England Electric Co.
4 Pacific Gas & Electric Co.
5 San Diego Gas & Electric Co.
6 United Illuminating Co.
7 Virginia Electric & Power Co.base_original |> filter(grupo == 3) |> select(Company) Company
1 Commonwealth Edison Co.
2 Consolidated Edison Co. (N.Y.)
3 Kentucky Utilities Co.
4 Madison Gas & Electric Co.
5 Northern States Power Co.
6 Wisconsin Electric Power Co.base_original |> filter(grupo == 4) |> select(Company) Company
1 Idaho Power Co.
2 Nevada Power Co.
3 Puget Sound Power & Light Co.X1: Índice de cobertura de encargos fixos (renda/dívida).
boxplot(base_original$X1 ~ base_original$grupo)
O gráfico do boxplot da variável \(X_1\) mostra que os grupos 1 e 3 são aqueles com maiores índice de cobertura de encargos fixos.
X2: Taxa de retorno sobre o capital
boxplot(base_original$X2 ~ base_original$grupo)
O gráfico do boxplot da variável \(X_2\) mostra que os grupos 1 e 3 são aqueles com maiores taxas de retorno sobre capital.
X3: Custo por kW de capacidade instalada
boxplot(base_original$X3 ~ base_original$grupo)
O gráfico do boxplot da variável \(X_3\) mostra que as empresas do grupo 4 tem maior custo por kW de capacidade instalada e as empresas do grupo 1 tem o menor valor para esta variável.
X4: Fator de carga anual
boxplot(base_original$X4 ~ base_original$grupo)
O gráfico do boxplot da variável \(X_4\) mostra que as empresas do grupo 2 apresentam maior fator de carga anual quando comparadas com as empresas dos outros grupos.
X5: Crescimento da demanda de pico de kWh de 1974 para 1975
boxplot(base_original$X5 ~ base_original$grupo)
As empresas do grupo 4 apresentam maior Crescimento da demanda de pico de kWh de 1974 para 1975, e as empresas do grupo 1 os menores crescimentos.
X6: Vendas (consumo anual de kWh)
boxplot(base_original$X6 ~ base_original$grupo)
As empresas do grupo 4 apresentam maior consumo anual de kWh e as empresas dos grupos 2 e 3 os menores consumos.
X7: Percentual de energia nuclear
boxplot(base_original$X7 ~ base_original$grupo)
As empresas do grupo 3 apresentam maior percentual de energia nuclear, e as empresas dos grupos 1 e 4 os menores percentuais.
X8: Custo total de combustível
boxplot(base_original$X8~ base_original$grupo)
As empresas dos grupos 3 e 4 apresentam maior custo total de combustível e as empresas do grupo 2 os menores custos. Vale destacar que no grupo 1 esse custo tem maior variabilidade que nos demais grupos.
Conclusão
A análise comparativa mostrou que o método hierárquico com medida de distância entre grupos “woard.D2” e quatro grupos apresentou melhor desempenho em termos de separação entre os grupos (índice de silhueta). A escolha entre os métodos deve considerar tanto os critérios quantitativos quanto a interpretabilidade dos agrupamentos.
Uma análise exploratório permitiu entender as diferenças entre os grupos e com isso concluir as seguintes características de cada um deles.
Quando comparadas com as empresas dos outros grupos, as empresas do grupo 1 apresentam índice de cobertura de encargos fixos altos, taxa de retorno sobre o capita altas, custo por kW de capacidade instalada baixo, fator de carga anual baixo, crescimento da demanda de pico de kWh de 1974 para 1975 baixo, percentual de energia nuclear baixo e uma grande variabilidade no custo total de combustível.
Quando comparadas com as empresas dos outros grupos, as empresas do grupo 2 apresentam índice de cobertura de encargos fixos baixos, fator de carga anual alto e consumo anual de kWh baixo.
Quando comparadas com as empresas dos outros grupos, as empresas do grupo 3 apresentam índice de cobertura de encargos fixos alto, taxa de retorno sobre o capital alta, fator de carga anual com grande variabilidade, consumo anual de kWh baixo e percentual de energia nuclear baixo.
Quando comparadas com as empresas dos outros grupos, as empresas do grupo 4 apresentam índice de cobertura de encargos fixos baixo, taxa de retorno sobre o capital baixa, custo por kW de capacidade instalada alto, crescimento da demanda de pico de kWh de 1974 para 1975 alto, consumo anual de kWh alto e percentual de energia nuclear baixo.
De fato o agrupamento apresentou uma partição dos dados de forma que os grupos de empresas apresentam diferentes características.