
######## Análisis de clasificación con "hclust" de la "Enquesta Metropolitana"

# Matriz de datos del ejemplo (RMB1986F2.rda) con 4 factores rotados: RC1 a RC4
load(url("http://pagines.uab.cat/plopez/sites/pagines.uab.cat.plopez/files/RMB1986F2.rda"))
attach(RMB1986F2)
# Variable añadida en el análisis de clasificación
ZSM=(seccion_municipio-mean(seccion_municipio))/sd(seccion_municipio)
RMB1986F2=data.frame(RMB1986F2,ZSM)

# Matriz de distancias
dist=dist(RMB1986F2[,c("RC1","RC2","RC3","RC4","ZSM")],method="euclidean") 
# Clasificación jerárquica ascendente con ward
ward=hclust(d=dist,method="ward.D2")
# Dendrograma
plot(ward)
rect.hclust(ward,k=8,border="red")
# Proceso de agregación
ward$merge[1:20,]
ward$merge[3490:3508,]
# Coeficiente de distancia a la que se forman los grupos
d=round(ward$height[3490:3508],3)
d
# Incremento o diferencias primeras (resta cada valor del anterior) 
d1=c(0,d[-length(d)])
round(d1,3)
round(d-d1,3)

# Determinación automàtica del número de grupos: Nbclust tiene 30 procedimientos distintos
# install.packages("NbClust")
library(NbClust)
attach(RMB1986F2)
DCL=data.frame(RC1,RC2,RC3,RC4,ZSM)
# Método 1 de Krzanowski and Lai (1988): kl
ngruposkl=NbClust(data=DCL,diss=NULL, distance="euclidean", min.nc=2, max.nc=10, method="ward.D", index="kl")
ngruposkl$All.index
ngruposkl$Best.nc
# Todos los métodos: all
# Ojo: La ejecución tarda unos 7 minutos
ngrupos<-NbClust(data=DCL,diss=NULL, distance="euclidean", min.nc=2, max.nc=10, method="ward.D", index="all")
round(res$Best.nc, digits=0)




fviz_nbclust(DCL, hcut, method = c("silhouette", "wss", "gap_stat"))
fviz_nbclust(DCL, hcut, method = "wss")
fviz_nbclust(DCL, hcut, method = "silhouette")
fviz_nbclust(DCL, hcut, method = "gap_stat")

# Corte del dendrograma en k grupos
G8=cutree(ward,k=8)
Grupos8=as.factor(G8)

# Matriz de datos con las clasificaciones
RMB1986F2C=data.frame(RMB1986F2,Grupos8)
save(RMB1986F2C,file="C:/Datos/RMB1986F2C.rda")

# Tabla de frecuencias y cruces
# install.packages(Deducer)
library(Deducer)
frequencies(RMB1986F2C[c("Grupos8")],r.digits=1)
round(t(aggregate(RMB1986F2C[,5:34],by=list(G8),FUN=mean)),1)
descriptive.table(vars=d(RMB1986F2C[,5:34]),data=RMB1986F2C,func.names=c("Mean"))

# Unidades según el grupo en el gráfico factorial
ggplot(data=RMB1986F2C) + geom_point(aes(x=RC1,y=RC2,colour=Grupos8)) +
  geom_vline(xintercept=0.0) +  geom_hline(yintercept=0.0)
ggplot(data=RMB1986F2C) + geom_point(aes(x=RC1,y=RC3,colour=Grupos8)) +
  geom_vline(xintercept=0.0) +  geom_hline(yintercept=0.0)
ggplot(data=RMB1986F2C) + geom_point(aes(x=RC3,y=RC4,colour=Grupos8)) +
  geom_vline(xintercept=0.0) +  geom_hline(yintercept=0.0)





# Análisis de correspondencias múltiples
# install.packages("FactoMineR")
library(FactoMineR)
# install.packages("psych")
library(psych) # Para cor.plot, cortest.bartlett, principal

load(url("http://pagines.uab.cat/plopez/sites/pagines.uab.cat.plopez/files/RMB1986.rda"))
attach(RMB1986)
DF=data.frame(menores,mayores,inmigracion,llegados,analfabetos,univesitarios,escolarizacion14_24,
              parados_ocupados,parados_1erempleo,mujeresactivas,profesion_alta,profesion_baja,
              terciario_medio, terciario_alto,agropecuario,vehiculotrabajo,seccion_municipio)
ACP2=principal(DF, nfactors= 4, rotate="varimax", covar=FALSE, scores=TRUE)
print(ACP2,sort=T)
hcpc=HCPC(ACP2,min=2,method="ward.D2",nb.clus=-1)    # Elección automática de la partición


ACP=PCA(DF,ncp=2)
summary(ACP)
# Análisis de clasificación jerárquica
hcpc=HCPC(ACP,min=2,method="ward.D2",nb.clus=-1)    # Elección automática de la partición
hcpc=HCPC(ACP,min=2,method="ward.D2",nb.clus=0)     # Para elegir la partición
# Elección del número óptimo de grupos o tipos de la partición
hcpc$call$t$nb.clust
# Descripción de las clases de la clasificación (tipos de la tipología)
hcpc$desc.var$test.chi2   # Variables
hcpc$desc.var$category    # Categorías
hcpc$desc.axes            # Factores
hcpc$desc.ind             # individuos

