##########################################################################################
##                                        V2                                            ##
##                         Ana Flávia Delbem Vidigal Nazareth                           ##
##                                                                                      ##
##                   ANÁLISE DE AGRUPAMENTOS COM VARIÁVEIS MISTAS                       ##
##                       ***        K-PROTÓTIPOS         ***                            ##
##                       BANCO DE DADOS GEOLÓGICO-GEOTÉCNICO                            ##
##                                                                                      ##
##########################################################################################

options(max.print = 5.5E5)
par(mar = c(4, 10, 1, 1)) #bottom, left, top, right

load <- function(pkg) {
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg))
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)}

packages <- c("devtools",
              "fastDummies",
              "NbClust",
              "clustMixType",
              "cluster",
              "factoextra",
              "FactoMineR",
              "xda",
              "corrplot",
              "RColorBrewer",
              "plyr",
              "Rcpp",
              "stats")

load(packages)


##########################################################################################
##               FUNÇÕES MODIFICADAS A SEREM UTILIZADAS PARA AS ANÁLISES                ##
##########################################################################################

source("clustMixType modified functions.R")

##########################################################################################
##                    OBTENÇÃO E MANIPULAÇÃO DA MATRIZ DE DADOS                         ##
##########################################################################################

##  As etapas do pré-tratamento do banco de dados de transformação de variáveis qualitativas nominais em variáveis binárias e normalização das variáveis numéricas são as únicas etapas do pré-tratamento realizadas com o auxílio do R. As demais etapas foram realizadas utilizando o Excel e macros em linguagem VBA.

#1# Obtenção das variáveis nominais originais para a transformação em variáveis binárias
#dados = read.table("01_01_HIGIENIZACAO8.DATA", header=T, sep = ';',na.strings='?',
#  colClasses = c("factor","factor","factor"))
                  #ESTRUTURA FAMÍLIA	LITOTIPO

## Criação de variáveis Dummy para as variáveis qualitativas nominais
#dados.dummy <- fastDummies::dummy_cols(dados)
#write.csv(dados.dummy, file = "01_01_HIGIENIZACAO8_DUMMY.csv") #Transformação salva em uma planilha externa


#2# Obtenção dos dados para a normalização das variáveis numéricas e realização da análise de agrupamentos
##TESTE
#dados = read.table("TESTE.DATA", header = T, sep = ';', na.strings = '?',
#  colClasses = c("factor", "factor", "factor", "factor", "numeric", "numeric", "numeric", "numeric"))

dados = read.table("DADOS_FINAL.DATA", header = T, sep = ';', na.strings = '?',                 
   colClasses = c("factor", #PONTO
    "factor","factor","factor","factor","factor","factor","factor","factor","factor","factor",
    "factor","factor","factor","factor","factor","factor","factor","factor","factor","factor", #Litotipos 1 a 20
    "factor","factor","factor","factor", #Tipos das estruturas 1 a 4
    "factor","factor","factor","factor","factor","factor","factor","factor", #Famílias das estruturas 1 a 8
    "numeric","numeric","numeric","numeric","numeric","numeric","numeric","numeric","numeric","numeric","numeric","numeric")) #X,Y,Z,R,W,ESPAÇAMENTO,PERSISTÊNCIA,ABERTURA,RUGOSIDADE,PREENCHIMENTO,ALTERAÇÃO_PAREDES,JW

## Separação das variáveis numéricas
var.num <- sapply(dados, is.numeric) #Retorna um vetor de TRUE ou FALSE para a condição.
dados.num <- dados[,var.num]

## Separação das variáveis categóricas nominais
var.cat <- sapply(dados, is.factor) #Retorna um vetor de TRUE ou FALSE para a condição.
dados.cat <- dados[,var.cat] #Separa a matriz de dados de acordo com o tipo das variáveis.
excluir <- c("PONTOS") #Excluir a variável contendo os nomes das observações do banco de dados
#excluir <- c("SETORES")
dados.cat <- dados.cat[,!(names(dados.cat)%in% excluir)]

#dados.setores <- dados$SETORES

## Normalização das variáveis numéricas: (x - min(x))/(max(x) - min(x))
j <- ncol(dados.num)
dados.num.n <- as.data.frame(apply(dados.num[, 1:j], 2, function(x) (x - min(x))/(max(x)-min(x))))

## Junção das matrizes, criando um novo banco de dados onde as variáveis numéricas estão normalizadas

## Converter as variáveis nominais de character para factor
#for (i in 1:ncol(dados.cat)) {
#  if (is.character(dados.cat[, i])) {
#    dados.cat[, i] = factor(dados.cat[, i])}}

dados.num.n <- as.data.frame.matrix(dados.num.n) #Matriz para dataframe

dados1 <- cbind(dados.cat, dados.num.n) #Data frame a ser utilizado no agrupamento
dados2 <- cbind(dados.cat, dados.num)

## As primeiras variáveis serão sempre as nominais seguidas das variáveis numéricas.
## Esta ordem é necessária para o correto funcionamento da função kproto.modif.


##########################################################################################
##                          ESTATÍSTICA DESCRITIVA DOS DADOS                            ##
##########################################################################################

## Medidas descritivas das variáveis
summary(dados1)

## Medidas descritivas das variáveis numéricas
dados.num.summary <- numSummary(dados.num)
dados.num.summary
write.csv(dados.num.summary, file = "C:/Users/anafl/Desktop/dados.num.summary.csv", col.names = NA)

## Medidas descritivas das variáveis numéricas normalizadas
dados.num.n.summary <- as.data.frame(dados.num.n)
dados.num.n.summary <- numSummary(dados.num.n.summary)
dados.num.n.summary
write.csv(dados.num.n.summary, file = "C:/Users/anafl/Desktop/dados.num.n.summary.csv", col.names = NA)

## Desvio padrão médio das variáveis numéricas
mean(x = dados.num.n.summary$sd)

## Boxplot dos dados numéricos em escala logarítmica
dados.num_log = dados.num
dados.num_log[(dados.num_log == 0)] <- 1

par(mfrow=c(1,1), mar = c(4, 7.7, 0.1, 0.1),cex=1.5)

boxplot(dados.num_log,
  #main = "BoxPlot - Variáveis Numéricas",
  horizontal = TRUE,
  log = "x",
  xlab = "Variáveis Numéricas (escala logarítmica)",
  #ylab = "Variáveis",
  col = "#FF1F1F",
  las = 1)

## Boxplot das variáveis numéricas normalizadas em escala aritmética
boxplot(dados.num.n,
  #main = "BoxPlot - Variáveis Numéricas Normalizadas",
  horizontal = TRUE,
  xlab = "Variáveis Numéricas Normalizadas",
  #ylab = "Variáveis",
  col = "#FF1F1F",
  las = 1)

## Matriz de correlações das variáveis numéricas normalizadas
#Igual a matriz de covariância dos dados padronizados R(dados) = S(dados_escala)
R = cor(dados.num.n, method = 'spearman')
R

## Função para cálculo do p-valor das correlações
cor.mteste <- function(mat, ...) {
  mat <- as.matrix(mat)
  n <- ncol(mat)
  p.mat <- matrix(NA, n, n)
  diag(p.mat) <- 0
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      tmp <- cor.test(mat[, i], mat[, j], ...)
      p.mat[i, j] <- p.mat[j, i] <- tmp$p.value}}
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat}

## Matriz dos p-valores das correlações
p.mat <- cor.mteste(R)
p.mat

## Correlograma das variáveis normalizadas - correlações com nível de significância p<0,05 (5%) ou intervalo de confiança > 95%
par(mfrow=c(1,1), mar = c(1, 1, 1, 1),cex=1.5)
corrplot(R, method = "color",
         type = "upper", order = "hclust", number.cex = .5,
         addCoef.col = "gray20",
         tl.col = "black", tl.srt = 45,
         tl.cex=0.7, cl.cex=0.6,
         p.mat = p.mat, sig.level = 0.05, insig = "blank",
         pch.cex = .3, diag = FALSE)

## Histogramas de densidade das variáveis numéricas, originais e normalizadas
par(mfrow=c(3,4), mar = c(2, 3, 2, 0.1), cex=0.45, cex.main=2.12)
#dados.num
hist(dados.num$X, freq = F, ylab = "", xlab="", main = "X",col="dodgerblue4")
abline(v=median(dados.num$X), col="green",lwd=1)
abline(v = mean(dados.num$X), col = "red", lwd = 1)
hist(dados.num$Y, freq = F, ylab = "", xlab="", main = "Y",col="dodgerblue4")
abline(v=median(dados.num$Y), col="green",lwd=1)
abline(v = mean(dados.num$Y), col = "red", lwd = 1)
hist(dados.num$Z, freq = F, ylab = "", xlab="", main = "Z",col="dodgerblue4")
abline(v=median(dados.num$Z), col="green",lwd=1)
abline(v = mean(dados.num$Z), col = "red", lwd = 1)
hist(dados.num$R, freq = F, ylab = "", xlab="", main = "R",col="dodgerblue4")
abline(v=median(dados.num$R), col="green",lwd=1)
abline(v = mean(dados.num$R), col = "red", lwd = 1)
hist(dados.num$W, freq = F, ylab = "", xlab="", main = "W",col="dodgerblue4")
abline(v=median(dados.num$W), col="green",lwd=1)
abline(v = mean(dados.num$W), col = "red", lwd = 1)
hist(dados.num$Espaçamento, freq = F, ylab = "", xlab="", main = "Espaçamento",col="dodgerblue4")
abline(v=median(dados.num$Espaçamento), col="green",lwd=1)
abline(v = mean(dados.num$Espaçamento), col = "red", lwd = 1)
hist(dados.num$Persistência, freq = F, ylab = "", xlab="", main = "Persistência",col="dodgerblue4")
abline(v=median(dados.num$Persistência), col="green",lwd=1)
abline(v = mean(dados.num$Persistência), col = "red", lwd = 1)
hist(dados.num$Abertura, freq = F, ylab = "", xlab="", main = "Abertura",col="dodgerblue4")
abline(v=median(dados.num$Abertura), col="green",lwd=1)
abline(v = mean(dados.num$Abertura), col = "red", lwd = 1)
hist(dados.num$Rugosidade, freq = F, ylab = "", xlab="", main = "Rugosidade",col="dodgerblue4")
abline(v=median(dados.num$Rugosidade), col="green",lwd=1)
abline(v = mean(dados.num$Rugosidade), col = "red", lwd = 1)
hist(dados.num$Preenchimento, freq = F, ylab = "", xlab="", main = "Preenchimento",col="dodgerblue4")
abline(v=median(dados.num$Preenchimento), col="green",lwd=1)
abline(v = mean(dados.num$Preenchimento), col = "red", lwd = 1)
hist(dados.num$AlteraçãoParedes, freq = F, ylab = "", xlab="", main = "Alteração das Paredes",col="dodgerblue4")
abline(v=median(dados.num$AlteraçãoParedes), col="green",lwd=1)
abline(v = mean(dados.num$AlteraçãoParedes), col = "red", lwd = 1)
hist(dados.num$JW, freq = F, ylab = "", xlab="", main = "JW",col="dodgerblue4")
abline(v=median(dados.num$JW), col="green",lwd=1)
abline(v = mean(dados.num$JW), col = "red", lwd = 1)
#legend("topleft", legend=c("Mediana", "Média"), fill=c("green", "red"), bty="n")

par(mfrow=c(3,4), mar = c(2, 3, 2, 0.1), cex=0.45, cex.main=2.12)
#dados.num.n
hist(dados.num.n$X, freq = F, ylab = "", xlab="", main = "X",col="dodgerblue4")
abline(v=median(dados.num.n$X), col="green",lwd=1)
abline(v = mean(dados.num.n$X), col = "red", lwd = 1)
hist(dados.num.n$Y, freq = F, ylab = "", xlab="", main = "Y",col="dodgerblue4")
abline(v=median(dados.num.n$Y), col="green",lwd=1)
abline(v = mean(dados.num.n$Y), col = "red", lwd = 1)
hist(dados.num.n$Z, freq = F, ylab = "", xlab="", main = "Z",col="dodgerblue4")
abline(v=median(dados.num.n$Z), col="green",lwd=1)
abline(v = mean(dados.num.n$Z), col = "red", lwd = 1)
hist(dados.num.n$R, freq = F, ylab = "", xlab="", main = "R",col="dodgerblue4")
abline(v=median(dados.num.n$R), col="green",lwd=1)
abline(v = mean(dados.num.n$R), col = "red", lwd = 1)
hist(dados.num.n$W, freq = F, ylab = "", xlab="", main = "W",col="dodgerblue4")
abline(v=median(dados.num.n$W), col="green",lwd=1)
abline(v = mean(dados.num.n$W), col = "red", lwd = 1)
hist(dados.num.n$Espaçamento, freq = F, ylab = "", xlab="", main = "Espaçamento",col="dodgerblue4")
abline(v=median(dados.num.n$Espaçamento), col="green",lwd=1)
abline(v = mean(dados.num.n$Espaçamento), col = "red", lwd = 1)
hist(dados.num.n$Persistência, freq = F, ylab = "", xlab="", main = "Persistência",col="dodgerblue4")
abline(v=median(dados.num.n$Persistência), col="green",lwd=1)
abline(v = mean(dados.num.n$Persistência), col = "red", lwd = 1)
hist(dados.num.n$Abertura, freq = F, ylab = "", xlab="", main = "Abertura",col="dodgerblue4")
abline(v=median(dados.num.n$Abertura), col="green",lwd=1)
abline(v = mean(dados.num.n$Abertura), col = "red", lwd = 1)
hist(dados.num.n$Rugosidade, freq = F, ylab = "", xlab="", main = "Rugosidade",col="dodgerblue4")
abline(v=median(dados.num.n$Rugosidade), col="green",lwd=1)
abline(v = mean(dados.num.n$Rugosidade), col = "red", lwd = 1)
hist(dados.num.n$Preenchimento, freq = F, ylab = "", xlab="", main = "Preenchimento",col="dodgerblue4")
abline(v=median(dados.num.n$Preenchimento), col="green",lwd=1)
abline(v = mean(dados.num.n$Preenchimento), col = "red", lwd = 1)
hist(dados.num.n$AlteraçãoParedes, freq = F, ylab = "", xlab="", main = "Alteração das Paredes",col="dodgerblue4")
abline(v=median(dados.num.n$AlteraçãoParedes), col="green",lwd=1)
abline(v = mean(dados.num.n$AlteraçãoParedes), col = "red", lwd = 1)
hist(dados.num.n$JW, freq = F, ylab = "", xlab="", main = "JW",col="dodgerblue4")
abline(v=median(dados.num.n$JW), col="green",lwd=1)
abline(v = mean(dados.num.n$JW), col = "red", lwd = 1)
#legend("topleft", legend=c("Mediana", "Média"), fill=c("green", "red"), bty="n")


##########################################################################################
##           ÍNDICES ESTATÍSTICOS PARA DETERMINAÇÃO DO MELHOR NÚMERO DE GRUPOS          ##
##########################################################################################

## Determinando o numero ótimo de grupos a partir das variáveis numéricas apenas

# Método do cotovelo
fviz_nbclust(dados.num.n, kmeans, nstart = 30, method = "wss") + geom_vline(xintercept = 4, linetype = 2) + labs(subtitle = "Elbow method")

# Várias estatísticas - Regra da maioria
nb <- NbClust(dados.num.n, distance = "euclidean", min.nc = 2, max.nc = 10, method = "kmeans")
fviz_nbclust(nb)

# Observar a tendência de agrupamentos graficamente (estatística Hopkins)
get_clust_tendency(dados.num.n, 394, graph = TRUE, gradient = list(low = "red", mid ="white", high = "blue"))

# Método da silhueta - k que maximiza a silhueta média
fviz_nbclust(dados.num.n, kmeans, nstart = 30, method = "silhouette") + labs(subtitle = "Silhouette method")

#Gráfico da silhueta das observações
par(mar = c(4, 3, 3, 3)) #bottom, left, top, right
pam.res <- pam(dados.num.n, 3)
fviz_silhouette(silhouette(pam.res))
par(mar = c(4, 10, 1, 1)) #bottom, left, top, right

# Análise das cargas nas duas primeiras componentes principais (CPs)
#A partir dos resultados da função kproto.modif, mas apenas as variáveis numéricas
dados1.kproto <- kproto.modif(
  dados1,
  k = 3,
  lambda = 1,
  iter.max = 100000,
  nstart = 30, #Número de repetições com protótipos iniciais diferentes. O resultado com distâncias totais mínimas é armazenado.
  na.rm = TRUE,
  keep.data = FALSE,
  verbose = TRUE)
#Gráfico das cargas de cada observação nas duas primeiras componentes principais
#fviz_cluster()
fviz_cluster.modif(dados1.kproto) #Componentes principais apenas com as variáveis numéricas
# Cargas das variáveis (2 primeiras CPs)
PCA.num <- PCA(dados.num.n)

## Determinando o numero ótimo de grupos a partir de todas as variáveis

#Utilizando a função "kproto" modificada com o coeficiente de similaridade desejado
Es <- numeric(20) #for(i in 1:k) e plot(1:k...) : usar o mesmo valor de k abaixo
for (i in 1:20) { #for(i in 1:k) e plot(1:k...) = Limite do número de grupos do gráfico scree plot
  kpres <- kproto.modif(dados1,
                        k = i,
                        lambda = 1, #Alterar os pesos se desejado.
                        iter.max = 10000,
                        nstart = 30)
  Es[i] <- kpres$tot.withinss}

#Plotagem do gráfico Scree Plot
par(mfrow=c(1,1), mar = c(4, 4, 1, 1), cex=1.5, cex.axis=0.8)

plot(1:20, #for(i in 1:k) e plot(1:k...) : usar o mesmo valor de k acima
     Es,
     pch = 16,
     col = "black",
     ylab = "Função Objetivo",
     xlab = "Número de Grupos",
     xaxt='n',yaxt='n',
     type = "o",
     main = "Scree Plot - k-Protótipos")
axis(side=1, at=c(0:20))
axis(side=2, at=seq(0, 900, by=100))
abline(v = 5, col = "black", lty = 2)
abline(v = 8, col = "black", lty = 2)


##########################################################################################
##         PARÂMETROS DO ALGORITMO K-PROTÓTIPOS A SEREM UTILIZADOS NA ANÁLISE           ##
##########################################################################################

## Número de Grupos (k):
K = 3

## Peso lambda:
## Estimativa do peso das variáveis numéricas e qualitativas (binárias) a partir da função lambdaest (compara as variâncias para especificar os pesos)
#lambdaest(dados1, num.method = 1, fac.method = 1, outtype = "numeric") #Peso único às variáveis qualitativas.
#lambdaest(dados1, num.method = 1, fac.method = 1, outtype = "vector") #Diferentes pesos para cada variável.
#num.method e fac.method heurística: 1 para variância e 2 para desvio padrão.  
#"Ratio of averages over all numeric/factor variables is returned. In case of outtype = "vector" the separate lambda for all variables is returned as the inverse of the single variables' variation as specified by the num.method and fac.method argument. outtype = "variation" directly returns these quantities and is not ment to be passed directly to kproto()."  

## Peso às variáveis qualitativas (binárias) (LAMB=0: só considera as variáveis numéricas)
LAMB = 1

## Pesos a variáveis específicas
#Vetor de pesos específicos de cada variável, na ordem em que são apresentadas.
#Neste caso as distâncias de cada variável serão multiplicadas pelo seu peso correspondente.
LAMB <- c( #O primeiro termo refere-se ao peso das variáveis categóricas, aplicado à distância entre as variáveis categóricas em conjunto. Os demais termos são os pesos dados individualmente a cada variável quantitativa. O tamanho do vetor deve ser o número de variáveis numéricas mais um.
  #Peso das variáveis qualitativas
  1, #Estrutura (1 a 4)
         #Família da estrutura (1 a 8)
         #Litotipo (1 a 18)
  #Pesos das variáveis numéricas
  0, 0, 0, #Coordenadas espaciais (X;Y;Z)
  1, 1, 1, 1, 1, 1, 1, 1, 1) #Qualitativas ordinais (R;W;ESPAÇAMENTO;PERSISTÊNCIA;ABERTURA;RUGOSIDADE;PREENCHIMENTO;ALTERAÇÃO_PAREDES;JW)

##Vetor LAMB para a função kproto()
#LAMB <- c(
#  #Peso das variáveis qualitativas
#  1,1,1,1, #Estrutura (1 a 4)
#  1,1,1,1,1,1,1,1,#Família da estrutura (1 a 8)
#  1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,#Litotipo (1 a 18)
#  #Pesos das variáveis numéricas
#  4, 4, 4, #Coordenadas espaciais (X;Y;Z)
#  1, 1, 1, 1, 1, 1, 1, 1, 1) #Qualitativas ordinais (R;W;ESPAÇAMENTO;PERSISTÊNCIA;ABERTURA;RUGOSIDADE;PREENCHIMENTO;ALTERAÇÃO_PAREDES;JW)

##TESTE
#LAMB <- c(0.5, 0.5, 1, 1, 1)

##########################################################################################
##                 AGRUPAMENTOS PELO MÉTODO NÃO HIERÁRQUICO K-PROTÓTIPOS                ##
##########################################################################################

#Função para o algoritmo de agrupamentos k-Protótipos, modificada dos pacotes "ClustMixtype", função "kproto", e "ade4", função "dist.binary".
dados1.kproto <- kproto.modif(
  dados1,
  k = K,
  lambda = LAMB,
  iter.max = 100000,
  nstart = 100, #Número de repetições com protótipos iniciais diferentes. O resultado com distâncias totais mínimas é armazenado.
  na.rm = TRUE,
  keep.data = FALSE,
  verbose = TRUE)

dados1.kproto[["withinss"]]
dados1.kproto[["dists"]]
dados1.kproto[["tot.withinss"]]
mean(dados1.kproto[["withinss"]])
dados1.kproto$centers

#Função original do pacote "ClustMixtype"
#dados2.kproto <- kproto( 
#  dados1,
#  k = K,
#  lambda = LAMB,
#  iter.max = 100000,
#  nstart = 100,
#  na.rm = TRUE,
#  keep.data = FALSE,
#  verbose = FALSE)

#dados2.kproto[["withinss"]]
#dados2.kproto[["tot.withinss"]]
#mean(dados2.kproto[["withinss"]])

## Sumarização dos resultados obtidos no agrupamento e visualização gráfica dos mesmos
summary.kp <- summary.kproto.modif(dados1.kproto, pct.dig = 3)
clprofiles.modif(dados1.kproto, dados1) # Boxplots para variáveis numéricas e Barplots para variáveis nominais de cada grupo são gerados.

## Extração dos grupos a que cada observação pertence
RC1 = as.data.frame(dados1.kproto$cluster) #Extrai um resultado (objeto) da lista resultante do Kproto.modif.
write.csv(RC1, file = "C:/Users/anafl/Desktop/RES29.csv", col.names = NA)


##########################################################################################
##                 VALIDAÇÃO ESTATÍSTICA DOS AGRUPAMENTOS GERADOS                       ##
##########################################################################################
load("RES00.RData")

## MANOVA
# Utilizando apenas as variáveis numéricas

y=cbind(dados.num.n[,1],dados.num.n[,2],dados.num.n[,3], #X, Y e Z - não entram nos resultados RES 6a11 15a17 21a23 e 27a29
        dados.num.n[,4], dados.num.n[,5],dados.num.n[,6],
        dados.num.n[,7],dados.num.n[,8],dados.num.n[,9],
        dados.num.n[,10],dados.num.n[,11],dados.num.n[,12])
#y12 <- y + 1e12 #Para resultados que apresentarem deficiência de posto (rank deficiency)

trat=factor(RC1[,1])
trat

m1=manova(y~trat)
#m1=manova(y12~trat)
m1

summary(m1)
rowSums(y)

#summary(m1)$SS
H = summary(m1)$SS$trat
#H
E = summary(m1)$SS$Res
#E

det(E)/det(E+H)

summary.manova(m1,test="Wilks")
summary(m1,test="Hotelling-Lawley")
summary(m1,test="Pillai")
summary(m1,test="Roy")
#?summary.manova


##########################################################################################
##                                FINALIZAÇÃO DO SCRIPT                                 ##
##########################################################################################

save.image(file = "RES30.RData", version = NULL, ascii = FALSE, compress = FALSE, safe = TRUE)
#load("RES30.RData")

rm(K, LAMB, dados1.kproto, summary.kp, RC1, y, trat, m1, H, E)

#Para obter novo resultado executar novamente a partir da linha de "PARÂMETROS DO ALGORITMO K-PROTÓTIPOS A SEREM UTILIZADOS NA ANÁLISE"

#rm(list=ls())




##########################################################################################
##                    EXEMPLO UTILIZANDO O BANCO DE DADOS IRIS                          ##
##########################################################################################

## Aquisição e tratamento dos dados
#Exemplo <- iris
#str(Exemplo)
## Padronização
#Exemplo.p <- as.data.frame(scale(Exemplo[-match(c('Species'), names(Exemplo))]))
#Exemplo.p1 <- cbind(Exemplo$Species, Exemplo.p)
#colnames(Exemplo.p1) = c("Species", "Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")
## Normalização
#Exemplo.n <- as.data.frame(apply(Exemplo[, 1:4], 2, function(x) (x - min(x))/(max(x)-min(x))))
#Exemplo.n1 <- cbind(Exemplo$Species, Exemplo.n)
#colnames(Exemplo.n1) = c("Species", "Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")
## Preparo do banco de dados normalizado para a função kproto.modif
#Exemplo.n2 <- fastDummies::dummy_cols(Exemplo.n1$Species) #Transformação das variáveil nominal em binárias
#colnames(Exemplo.n2) = c("Species", "Setosa", "Versicolor","Virginica")
#Exemplo.n2 <- cbind(Exemplo.n2[-match(c('Species'), names(Exemplo.n2))], Exemplo.n)
#str(Exemplo.n2)
#cols <- c("Setosa", "Versicolor","Virginica")
#Exemplo.n2[cols] <- lapply(Exemplo.n2[cols], factor) 
#str(Exemplo.n2)
#
## Descrição e vizualização dos dados
## scatterplot comprimentos vs. especies vs.largura petala (bubble)
#ggplot(Exemplo.n1, aes(x = Sepal.Length, y = Petal.Length, color = Species, size = Petal.Width), alpha = I(0.7)) + # alpha: reduz overplotting
#  geom_point() + xlab("Comprimento da sépala (cm)") +
#  ylab("Comprimento da pétala (cm)") + ggtitle("Conjunto de Dados Íris") +
#  scale_color_discrete(name ="Espécies") + scale_size_continuous(name = "Largura pétala") +
#  theme(legend.position = c(1, 0), legend.justification = c(1,0))
## scatterplot matrix com correlacoes e p-valor
## funcao para personalizacao do painel
#painel.cor <- function(x, y, digits = 2, cex.cor, ...){
#  usr <- par("usr"); on.exit(par(usr))
#  par(usr = c(0, 1, 0, 1))
## coeficiente de correlacao
#  r <- cor(x, y)
#  txt <- format(c(r, 0.123456789), digits = digits)[1]
#  txt <- paste("r= ", txt, sep = "")
#  text(0.5, 0.6, txt, cex = 1.2)
## calculo do p-valor
#  p <- cor.test(x, y)$p.value
#  txt2 <- format(c(p, 0.123456789), digits = digits)[1]
#  txt2 <- paste("p= ", txt2, sep = "")
#  if(p < 0.01) txt2 <- paste("p ", "<0.01", sep = "")
#  text(0.5, 0.4, txt2, cex = 1.2)}
#pairs(Exemplo.n1[,2:5], pch = 21,
#      bg = c("red", "green3", "blue")[unclass(Exemplo.n1$Species)],
#      upper.panel = painel.cor,
#      labels = c("Comprimento\nsepala", "Largura\nsepala",
#                 "Comprimento\npetala", "Largura\npetala"))
#legend('bottomright', legend = levels(Exemplo.n1$Species), col = 1:3, cex = 0.8, fill = unique(species.cor))
#
#par(mfrow=c(1,3), mar = c(4, 4, 1, 1), cex=0.5)
#plot(Exemplo$Sepal.Length, Exemplo$Sepal.Width, col=Exemplo$Species)
#plot(Exemplo.n1$Sepal.Length, Exemplo.n1$Sepal.Width, col=Exemplo.n1$Species)
#plot(Exemplo.p1$Sepal.Length, Exemplo.p1$Sepal.Width, col=Exemplo.p1$Species)
#par(mfrow=c(1,1))
#
## Determinação do melhor número de grupos
#fviz_nbclust(Exemplo.n, kmeans, nstart = 30, method = "wss") + geom_vline(xintercept = 3, linetype = 2) + labs(subtitle = "Elbow method")
#fviz_nbclust(Exemplo.n, kmeans, nstart = 30, method = "silhouette") + labs(subtitle = "Silhouette method")
#nb.mouse <- NbClust(Exemplo.n, distance = "euclidean", min.nc = 2, max.nc = 10, method = "kmeans")
#fviz_nbclust(nb.mouse)
#get_clust_tendency(Exemplo.n, 50, graph = TRUE, gradient = list(low = "red", mid ="white", high = "blue"))
#
## Análise pelo algoritmo de agrupamentos k-médias e pelo kproto.modif
#K.iris = 3
#LAMB.iris = c(0.5,1,1,1,1) #O primeiro é o peso único dado às variáveis nominais 
#
#kmeans.iris <- kmeans(Exemplo.n, centers = K.iris, iter.max = 10000, nstart = 30)
#
#kproto.iris <- kproto.modif(
#  Exemplo.n2,
#  k = K.iris,
#  lambda = LAMB.iris,
#  iter.max = 100000,
#  nstart = 30, #Número de repetições com protótipos iniciais diferentes. O resultado com distâncias totais mínimas é armazenado.
#  na.rm = TRUE,
#  keep.data = FALSE,
#  verbose = TRUE)
#
## Comparação visual dos grupos gerados
#par(mfrow=c(3,1), mar = c(4, 4, 1, 1), cex=1)
#plot(Exemplo.n1$Sepal.Length, Exemplo.n1$Petal.Length, col=Exemplo.n1$Species, pch = 19, xlab = "Comprimento Sepala", ylab = "Comprimento Petala")
#legend('topright', title = "Espécies", legend = levels(as.factor(Exemplo.n1$Species)), col = 1:3, cex = 0.8, pch = 19)
#plot(Exemplo.n$Sepal.Length, Exemplo.n$Petal.Length, col=kmeans.iris$cluster, pch = 19, xlab = "Comprimento Sepala", ylab = "Comprimento Petala")
#legend('topright', title = "Grupos Kmeans", legend = levels(as.factor(kmeans.iris$cluster)), col = 1:3, cex = 0.8, pch = 19)
#plot(Exemplo.n2$Sepal.Length, Exemplo.n2$Petal.Length, col=kproto.iris$cluster, pch = 19, xlab = "Comprimento Sepala", ylab = "Comprimento Petala")
#legend('topright', title = "Grupos Kproto.modif", legend = levels(as.factor(kproto.iris$cluster)), col = 1:3, cex = 0.8, pch = 19)


##########################################################################################
##                                         OUTROS                                       ##
##########################################################################################

## Plotagem da Figura 3
#dados = read.table("Mouse.data", header = T, sep = ';', na.strings = '?',                 
#                   colClasses = c("numeric","numeric","factor"))
#color_easy = c("red", "blue", "green")[dados$OBS]
#ggplot(dados)+geom_point(aes(x=X,y=Y,colour=OBS))+theme_bw()+labs(colour = "Grupos")


##########################################################################################
##                          CARACTERIZAÇÃO DOS SETORES PROPOSTOS                        ##
##########################################################################################
#Dados:
#Setores_Final_Pre-Trat.DATA : BANCO DE DADOS APÓS O PRÉ-TRATAMENTO, COM SETORES PROPOSTOS
#Setores_Final_HIG8.DATA : BANCO DE DADOS ANTERIOR À ETAPA DE UNIÃO DAS LINHAS DE DIFERENTES OBSERVAÇÕES EM UM ÚNICO PONTO, COM SETORES PROPOSTOS

setores = read.table("Setores_Final_Pre-Trat.DATA", header = T, sep = ';', na.strings = '?',        #Banco de dados com os setores propostos a que cada observação pertence         
                   colClasses = c("factor", #SETORES
                                  "factor","factor","factor","factor","factor","factor","factor","factor","factor","factor",
                                  "factor","factor","factor","factor","factor","factor","factor","factor","factor","factor", #Litotipos 1 a 20
                                  "factor","factor","factor","factor", #Tipos das estruturas 1 a 4
                                  "factor","factor","factor","factor","factor","factor","factor","factor", #Famílias das estruturas 1 a 8
                                  "numeric","numeric","numeric","numeric","numeric","numeric","numeric","numeric","numeric","numeric","numeric","numeric"))  #X,Y,Z,R,W,ESPAÇAMENTO,PERSISTÊNCIA,ABERTURA,RUGOSIDADE,PREENCHIMENTO,ALTERAÇÃO_PAREDES,JW


##Boxplots e barras dos setores propostos
vars <- 1:ncol(setores)
dados.setores <- setores$SETOR
clusids <- sort(unique(dados.setores))
nr <- nrow(count(levels(dados.setores)))
col <- brewer.pal(nr, "Set3")
size <- table(factor(dados.setores, levels=1:nr))

## Função para visualização dos dos perfis das variáveis nos setores propostos
## modificada a partir da função "clprofiles" do pacote "clustMixType"
#clprofiles.sectors <- function(dados.setores, setores, vars = NULL, col = NULL){
  if(length(dados.setores) != nrow(setores)) stop("Size of setores does not match cluster result!")
  if(is.null(vars)) vars <- 1:ncol(setores)
  if(!is.numeric(vars)) vars <- sapply(vars, function(z) return(which(colnames(setores)==z)))
  if(length(vars) < 1) stop("Specified variable names do not match setores!")
  if(is.null(col)){
    k <- max(unique(dados.setores))
    if(k > 2)  col <- brewer.pal(k, "Set3")
    if(k == 2) col <- c("lightblue","orange")
    if(k == 1) col <- "lightblue"}
  #clusids <- as.numeric(names(size)) # object size not named for kmeans
  clusids <- sort(unique(dados.setores)) 
  if(length(col) != nr) warning("Length of col should match number of clusters!")
  
  par(ask=TRUE)
  for(i in vars){
    if(is.numeric(setores[,i])){
      boxplot(setores[,i]~dados.setores, col = col, main = colnames(setores)[i])
      legend("topright", legend=clusids, fill = col)} 
    if(is.factor(setores[,i])){
      tab <- table(setores[,i], dados.setores)
      for(j in 1:length(size)) tab[,j] <- tab[,j]/size[j]
      barplot(t(tab), beside = TRUE, main = colnames(setores)[i], col = col, ylim=c(0,1))}} #Adição dos limites do eixo Y.
  par(ask="FALSE")