Estatística Multivariada - 5a lista de exercícios - Especialização em Estatística Aplicada, ICEX/UFMG Profa. Sueli Aparecida Mingoti
José Augusto Barros de Oliveira
16 de outubro de 2016
Introdução
O objetivo deste trabalho é refazer a lista de exercícios 5, prevista para o Minitab, usando o R.
Inicialização
A base de dados é carregada e os pacotes necessários são inicializados
library("NbClust", lib.loc="~/R/win-library/3.2")
library("ggplot2", lib.loc="~/R/win-library/3.2")
library("dplyr", lib.loc="~/R/win-library/3.2")
library("reshape2", lib.loc="~/R/win-library/3.2")
library("cluster", lib.loc="~/R/win-library/3.2")
vendedores <- read.csv2("D:/Guto/pessoal/especializacao/Multivariada/vendedores.csv")Escolha da região de referência
Gera a matriz de distâncias (euclidiana) e os agrupamentos pelo método hierárquico e ligação média. Não existe uma saída para esse pacote como a do Minitab então, para decidir sobre a região de referência, vamos usar o dendrograma e a relação de “alturas”, equivalente ao “distance levels”
mdist <- dist(vendedores[,2:8], method="euclidean")
cl1 <- hclust(mdist, method="average")
s1 <- data.frame(seq(44,2, by =-1),cl1$height)
colnames(s1) <- c("Num_clusters","altura")
s1## Num_clusters altura
## 1 44 0.0000000
## 2 43 0.7445804
## 3 42 1.2131364
## 4 41 1.2306096
## 5 40 1.2626955
## 6 39 1.2764404
## 7 38 1.3552859
## 8 37 1.4214546
## 9 36 1.6859715
## 10 35 1.7266767
## 11 34 1.7632709
## 12 33 1.8172782
## 13 32 2.0084853
## 14 31 2.0141499
## 15 30 2.1146867
## 16 29 2.1688737
## 17 28 2.2095022
## 18 27 2.2342336
## 19 26 2.3202172
## 20 25 2.4690064
## 21 24 2.5477087
## 22 23 2.6145363
## 23 22 2.6495127
## 24 21 2.8629977
## 25 20 2.9896298
## 26 19 3.0841398
## 27 18 3.1565503
## 28 17 3.2725220
## 29 16 3.3643315
## 30 15 3.3735687
## 31 14 3.3754290
## 32 13 3.4543007
## 33 12 3.7958243
## 34 11 4.2062027
## 35 10 4.4447759
## 36 9 4.4700781
## 37 8 4.5974823
## 38 7 5.6351266
## 39 6 6.2731284
## 40 5 6.6387747
## 41 4 6.9984098
## 42 3 9.1367681
## 43 2 12.8334707plot(cl1)
rect.hclust(cl1, h=6, border="red")
rect.hclust(cl1, h=9, border="blue")
Com altura = 6, teremos 6 clusters, dois dos quais com uma única observação. Se optarmos pela altura = 7, o número de clusters reduz para 3. Vamos optar por essa região de referência (3 <= k <= 6)
Cálculo do R2 e Pseudo F e escolha de k
sstc<- (nrow(vendedores)-1)*sum(apply(vendedores[,2:8],2,var))
for (i in c(3:6)) {
temp <- cutree(cl1, k=i)
vendedores <- cbind(vendedores, temp)
}
temp <- colnames(vendedores)
temp[9:12] <- c("k3","k4","k5","k6")
colnames(vendedores) <- temp
rm(temp)
# calculo do ssw para k = 3
ssw <- 0
por_k <- group_by(vendedores,k3)
temp <- summarise(por_k, var(Vendas), var(Lucro), var(Clientes), var(Escrita), var(Logica), var(Social), var(Matematica))
temp <- apply(temp[,2:8], 1, sum)
for (i in 1:3) {
temp2 <- as.vector((table(por_k$k3)[i]-1)*temp[i])
ssw <- ssw + temp2
}
sswk3 <- ssw
# calculo do ssw para k = 4
ssw <- 0
por_k <- group_by(vendedores,k4)
temp <- summarise(por_k, var(Vendas), var(Lucro), var(Clientes), var(Escrita), var(Logica), var(Social), var(Matematica))
temp <- apply(temp[,2:8], 1, sum)
temp[3] <-0
for (i in 1:4) {
temp2 <- as.vector((table(por_k$k4)[i]-1)*temp[i])
ssw <- ssw + temp2
}
sswk4 <- ssw
# calculo do ssw para k = 5
ssw <- 0
por_k <- group_by(vendedores,k5)
temp <- summarise(por_k, var(Vendas), var(Lucro), var(Clientes), var(Escrita), var(Logica), var(Social), var(Matematica))
temp <- apply(temp[,2:8], 1, sum)
temp[4] <-0
for (i in 1:5) {
temp2 <- as.vector((table(por_k$k5)[i]-1)*temp[i])
ssw <- ssw + temp2
}
sswk5 <- ssw
# calculo do ssw para k = 6
ssw <- 0
por_k <- group_by(vendedores,k6)
temp <- summarise(por_k, var(Vendas), var(Lucro), var(Clientes), var(Escrita), var(Logica), var(Social), var(Matematica))
temp <- apply(temp[,2:8], 1, sum)
temp[c(4,6)] <-0
for (i in 1:6) {
temp2 <- as.vector((table(por_k$k6)[i]-1)*temp[i])
ssw <- ssw + temp2
}
sswk6 <- ssw
#calculo dos R2 e Pseudo F
k <- seq(3,6)
ssw <- c(sswk3,sswk4,sswk5,sswk6)
r2 <- (sstc-ssw)/sstc
psf3 <- (44-3)/(3-1)*r2[1]/(1-r2[1])
psf4 <- (44-4)/(4-1)*r2[2]/(1-r2[2])
psf5 <- (44-5)/(5-1)*r2[3]/(1-r2[3])
psf6 <- (44-6)/(6-1)*r2[4]/(1-r2[4])
psf <- c(psf3,psf4,psf5,psf6)
saida <- cbind(k,r2,psf)
saida## k r2 psf
## [1,] 3 0.7401371 58.38775
## [2,] 4 0.7580311 41.77016
## [3,] 5 0.8401787 51.25563
## [4,] 6 0.8550736 44.84040Há uma redução de cerca de 10 pontos percentuais entre o R2 para k=5 e k=3. No entanto, o valor do pseudo F para k = 3 é maior. Considerando-se a esses resultados e a avaliação do dendrograma, optamos por k = 3.
O pacote NbClust referenda essa escolha por meio da comparação de 30 índices, entre eles o Pseudo F de Calinski e Harabasz (CH). 14 deles indicam o particionamento em 3 grupos:
NbClust(vendedores[,2:8], method = 'average', index = 'all')$Best.nc
## *** : The Hubert index is a graphical method of determining the number of clusters.
## In the plot of Hubert index, we seek a significant knee that corresponds to a
## significant increase of the value of the measure i.e the significant peak in Hubert
## index second differences plot.
## 
## *** : The D index is a graphical method of determining the number of clusters.
## In the plot of D index, we seek a significant knee (the significant peak in Dindex
## second differences plot) that corresponds to a significant increase of the value of
## the measure.
##
## *******************************************************************
## * Among all indices:
## * 6 proposed 2 as the best number of clusters
## * 14 proposed 3 as the best number of clusters
## * 1 proposed 13 as the best number of clusters
## * 2 proposed 15 as the best number of clusters
##
## ***** Conclusion *****
##
## * According to the majority rule, the best number of clusters is 3
##
##
## *******************************************************************## KL CH Hartigan CCC Scott Marriot
## Number_clusters 3.0000 3.0000 3.0000 2.0000 3.0000 3
## Value_Index 74.5285 58.3878 40.4687 20.5042 63.7147 17740029799
## TrCovW TraceW Friedman Rubin Cindex DB
## Number_clusters 3.0 3.0000 13.00 3.0000 2.0000 2.0000
## Value_Index 48236.8 533.6752 19068.82 -82.0399 0.3548 0.6416
## Silhouette Duda PseudoT2 Beale Ratkowsky Ball
## Number_clusters 2.0000 3.0000 3.0000 3.0000 2.000 3.0000
## Value_Index 0.4357 5.0747 -8.8324 -3.3607 0.435 377.8348
## PtBiserial Frey McClain Dunn Hubert SDindex Dindex
## Number_clusters 3.0000 1 2.0000 15.0000 0 3.0000 0
## Value_Index 0.5995 NA 0.2337 0.3633 0 0.4736 0
## SDbw
## Number_clusters 15.0000
## Value_Index 0.0268Composição dos grupos
p1 <- ggplot(vendedores, aes(x=k3, y=Funcionario, label=Funcionario)) + geom_point() +geom_text(aes(label=Funcionario),hjust=-1, vjust=0) + theme(text = element_text(size=9))
p1
vend <- melt(vendedores, id.vars=c(1,9), measure.vars=c(2:8))
p2 <- ggplot(vend, aes(x=variable, y=value)) + geom_boxplot()+facet_wrap(~ k3, ncol=3)+ theme(text = element_text(size=9), axis.text.x = element_text(angle=90, vjust=0))
p2
A inspeção dos boxplots revela que:
cada cluster apresenta as variáveis de desempenho (vendas, lucro e clientes) em faixas distintas de valor: menores no grupo 2, intermediários no 1 e maiores no grupo 3. Chama a atenção a média da variável lucro no grupo 3, muito destacada em relação aos demais grupos.
as variáveis escrita, lógica e social, associadas às habilidades dos funcionários, apresentam baixos valores e dispersão no grupo de pior desempenho (2), mas parecem não distinguir os outros grupos
os valores de matemática são muito diferentes entre os grupos, acompanhando o padrão das variáveis de desempenho: menores (média aprox = 8) para o grupo 2, intermediários (média aprox = 15) para o grupo 1 e maiores (média aprox = 21) para o grupo 3
Para confirmar essas percepções e caracterizar os clusters, vamos fazer a análise de variância das médias:
Variáveis de desempenho
summarise(group_by(vendedores,k3), mean(Vendas), mean(Lucro), mean(Clientes))## Source: local data frame [3 x 4]
##
## k3 mean(Vendas) mean(Lucro) mean(Clientes)
## (int) (dbl) (dbl) (dbl)
## 1 1 24.75565 26.17217 25.61304
## 2 2 21.73000 23.00750 24.01750
## 3 3 26.61154 29.87846 26.83385summary(aov(vendedores$Vendas ~ as.factor(vendedores$k3)))## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(vendedores$k3) 2 118.01 59.01 73.69 2.66e-14 ***
## Residuals 41 32.83 0.80
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(vendedores$Lucro ~ as.factor(vendedores$k3)))## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(vendedores$k3) 2 246.82 123.41 109.2 <2e-16 ***
## Residuals 41 46.35 1.13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(vendedores$Clientes ~ as.factor(vendedores$k3)))## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(vendedores$k3) 2 39.52 19.761 36.83 6.98e-10 ***
## Residuals 41 22.00 0.536
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ao nível de significância de 95%, as médias das variáveis de desempenho são diferentes entre os clusters
Variáveis de habilidades
summarise(group_by(vendedores,k3), mean(Escrita), mean(Logica), mean(Social), mean(Matematica))## Source: local data frame [3 x 5]
##
## k3 mean(Escrita) mean(Logica) mean(Social) mean(Matematica)
## (int) (dbl) (dbl) (dbl) (dbl)
## 1 1 5.565217 7.000000 5.347826 14.36957
## 2 2 4.375000 4.812500 4.125000 6.68750
## 3 3 6.576923 8.346154 5.846154 21.23077summary(aov(vendedores$Escrita ~ as.factor(vendedores$k3)))## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(vendedores$k3) 2 24.34 12.170 4.437 0.018 *
## Residuals 41 112.45 2.743
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(vendedores$Logica ~ as.factor(vendedores$k3)))## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(vendedores$k3) 2 61.84 30.919 23.19 1.83e-07 ***
## Residuals 41 54.66 1.333
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(vendedores$Social ~ as.factor(vendedores$k3)))## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(vendedores$k3) 2 14.94 7.471 8.806 0.000658 ***
## Residuals 41 34.78 0.848
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(vendedores$Matematica ~ as.factor(vendedores$k3)))## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(vendedores$k3) 2 1066.6 533.3 87.85 1.5e-15 ***
## Residuals 41 248.9 6.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ao nível de significância de 95%, as médias das variáveis de habilidade são diferentes entre os clusters
Componentes principais
Vamos usar as duas primeiras componentes principais para avaliar se a opção por três conglomerados representa bem o perfil dos funcionários
cp <- princomp(~ Vendas + Lucro + Clientes + Escrita + Logica + Social + Matematica, data= vendedores, cor=FALSE)
summary(cp)## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4
## Standard deviation 6.517031 1.85407870 1.00917117 0.98839299
## Proportion of Variance 0.879809 0.07121067 0.02109689 0.02023709
## Cumulative Proportion 0.879809 0.95101966 0.97211655 0.99235364
## Comp.5 Comp.6 Comp.7
## Standard deviation 0.445120128 0.371288522 0.1820212029
## Proportion of Variance 0.004104339 0.002855695 0.0006863296
## Cumulative Proportion 0.996457975 0.999313670 1.0000000000cp$loadings##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## Vendas -0.272 -0.124 -0.128 0.307 0.651 0.410 0.457
## Lucro -0.385 -0.125 0.420 -0.200 0.446 -0.645
## Clientes -0.162 -0.177 -0.279 0.115 -0.416 -0.447 0.692
## Escrita -0.123 -0.801 -0.403 -0.301 -0.286
## Logica -0.181 -0.423 0.655 0.412 -0.361 0.240
## Social -0.101 -0.345 0.751 -0.317 -0.449
## Matematica -0.832 0.340 -0.139 -0.167 -0.267 0.224 -0.152
##
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## SS loadings 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## Proportion Var 0.143 0.143 0.143 0.143 0.143 0.143 0.143
## Cumulative Var 0.143 0.286 0.429 0.571 0.714 0.857 1.000cor(cp$scores, vendedores[,2:8], method = "pearson")## Vendas Lucro Clientes Escrita Logica
## Comp.1 -0.95904047 -0.971902459 -0.89534688 -0.4551699015 -0.725676528
## Comp.2 -0.12374315 -0.089741445 -0.27759211 -0.8426456356 -0.482448177
## Comp.3 -0.06974484 0.164082461 -0.23792881 -0.2304813802 0.406371099
## Comp.4 0.16378759 -0.076645417 0.09589796 -0.1685326329 0.250445209
## Comp.5 0.15645107 0.076932471 -0.15641667 0.0004727412 -0.098810815
## Comp.6 0.08212275 -0.092803752 -0.14023574 0.0192994660 0.054862434
## Comp.7 0.04491331 -0.004834084 0.10657305 -0.0295141606 -0.001896405
## Social Matematica
## Comp.1 -0.61902666 -0.992189277
## Comp.2 0.05450522 0.115179678
## Comp.3 -0.32741500 -0.025653291
## Comp.4 0.69852277 -0.030125010
## Comp.5 0.02359277 -0.021712833
## Comp.6 -0.11055086 0.015207798
## Comp.7 -0.07695722 -0.005074356As duas primeiras componentes respondem por cerca de 95% da variância total. A primeira está fortemente correlacionada às variáveis de desempenho e a Lógica, Social e Matemática, compondo uma síntese do desempenho e esses atributos de habilidade. A segunda destaca os vendedores com boa avaliação em Escrita. Como os sinais são negativos, valores menores indicam melhores avaliações em ambas as componentes.
vendedores <- cbind(vendedores,cp$scores[,1:2])
p3 <- ggplot(vendedores, aes(x=Comp.1, y=Comp.2, color=as.factor(k3))) + geom_text(aes(label=Funcionario))
p3
Os conglomerados formados pelo particionamento em 3 se alinham perfeitamente aos scores da primeira componente, que, como se viu, é um índice do desempenho e de 3 das cinco habilidades. Concluimos, portanto, que a decisão foi adequada.
Ranking dos funcionários
Vamos criar um índice para ordenar os vendedores por desempenho
desempenho<- princomp(~ Vendas + Lucro + Clientes, data=vendedores, cor=TRUE)
summary(desempenho)## Importance of components:
## Comp.1 Comp.2 Comp.3
## Standard deviation 1.6667775 0.39262628 0.26018728
## Proportion of Variance 0.9260491 0.05138513 0.02256581
## Cumulative Proportion 0.9260491 0.97743419 1.00000000desempenho$loadings##
## Loadings:
## Comp.1 Comp.2 Comp.3
## Vendas -0.586 -0.158 0.795
## Lucro -0.576 -0.609 -0.545
## Clientes -0.570 0.778 -0.266
##
## Comp.1 Comp.2 Comp.3
## SS loadings 1.000 1.000 1.000
## Proportion Var 0.333 0.333 0.333
## Cumulative Var 0.333 0.667 1.000cor(desempenho$scores[1:44],vendedores[,2:4])## Vendas Lucro Clientes
## [1,] -0.9764127 -0.9606117 -0.9497319A primeira componente é responsável por 92,6% da variância e está fortemente correlacionada às variáveis originais de interesse. Vamos usar apenas a primeira, portanto. Como os coeficientes têm sinal negativo, baixos valores indicam melhor desempenho.
escore <- desempenho$scores[1:44]
vendedores <- cbind(vendedores, escore)
vendedores <- arrange(vendedores,escore)
#melhores
head(vendedores$Funcionario,5)## [1] 8 35 39 31 25#piores
vendedores <- arrange(vendedores,desc(escore))
head(vendedores$Funcionario,5)## [1] 41 44 16 23 2ranking pela primeira componente que utiliza desempenho e habilidade
vendedores <- arrange(vendedores,Comp.1)
#melhores
head(vendedores$Funcionario,5)## [1] 8 28 36 30 31#piores
vendedores <- arrange(vendedores,desc(Comp.1))
head(vendedores$Funcionario,5)## [1] 41 44 32 16 2Os melhores funcionários estão no cluster 3 e os piores no 2, independentemente do índice usado para o ordenamento. A diferença é que o índice que criamos por último (desempenho), leva em consideração apenas as variáveis Vendas, Lucro e Clientes, enquanto o primeiro (cp) considera também as variáveis de habilidade.
K Médias
O default nesse caso é determinar as sementes aleatoriamente entre os vetores disponíveis.
vendedores <- arrange(vendedores, Funcionario)
kmedia <- kmeans(vendedores[,2:8],3)
kmedia## K-means clustering with 3 clusters of sizes 14, 17, 13
##
## Cluster means:
## Vendas Lucro Clientes Escrita Logica Social Matematica
## 1 26.59071 29.78000 26.79714 6.535714 8.285714 5.857143 21.03571
## 2 25.05706 26.53529 25.79059 5.500000 7.000000 5.529412 15.38235
## 3 22.37923 23.57077 24.34462 4.884615 5.615385 4.307692 8.00000
##
## Clustering vector:
## [1] 3 3 2 2 2 3 2 1 2 1 2 2 2 2 2 3 2 2 2 2 3 1 3 1 1 3 1 1 3 1 1 3 2 3 1
## [36] 1 3 2 1 1 3 2 1 3
##
## Within cluster sum of squares by cluster:
## [1] 161.3543 155.0217 149.5806
## (between_SS / total_SS = 78.1 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"O R2 é apresentado na saída: 78,1% O Pseudo F é (44-3)/(3-1)*0.781/(1-0.781) = 73,10. As duas medidas são superiores àquelas obtidas pelo método hierárquico (R2= 74% e Pseudo F = 58,4)
km3 <- kmedia$cluster
vendedores <- cbind(vendedores,km3)
p4 <- ggplot(vendedores, aes(x=Comp.1, y=Comp.2, color=as.factor(k3), shape=as.factor(km3), size=5)) + geom_point()
p4
É possível observar que o método de K médias formou conglomerados mais homogêneos quanto ao número de elementos, alocando elementos do conglomerado intermediário obtido pelo método hierárquico (1) nos outros dois.
O gráfico seguinte é uma outra forma de se visualizar isso.
clusplot(vendedores[,2:8], kmedia$cluster, color=TRUE, shade=TRUE, labels=2, lines=0)