1 Carga de librerías y datos
Lo primero que voy a hacer es echar un ojo a los datos y configurarlos, como veo que hay 52 observaciones uso como variable identificativa el nombre de provincia. Hecho esto todas las variables son numéricas así que miro las observaciones y las principales medidas sobre las distribuciones de las variables.
library("readxl")
library("lattice")
library("pastecs")
library("corrplot")
library("ggplot2")
library("factoextra")
library("FactoMineR")
library("cluster")
library("NbClust")
library("RColorBrewer")
library("MASS")
library("gridExtra")
Provincias <- read_excel("Provincias.xlsx")
datos <- as.data.frame(Provincias)
rownames(datos)<- datos$Prov
datos$Prov <- NULL
datospsych::describe(Filter(is.numeric, datos))En primer lugar llama la atención la dimension de los datos ya que tenemos (52 filas por 18 columnas), luego el número de observaciones frente a variables es bastante reducido. Sobre las características de las variables mirando las diferencias entre la media y la mediana, o la curtósis y simetría vamos situaciones bastante dispares, así que no podemos llegar a ninguna conclusión previa al análisis.
2 Matriz de correlación
Vamos a representar la matriz de correlación entre las variables, de esta forma descubriremos que variables están más relacionadas.
Podemos observar como las variables: IPC, VS, Industria, TVF, Poblacion, CTH, Infor, Construccion, NumEmpresas, Ocupados, AFS, APT y PIB están positiva y fuertemente relacionadas entre sí, sobre todo a partir de Industria la correlación es superior a 0.75 entre todas las variables del conjunto.
Respecto a las variables más negativamente correlacionadas observamos que Mortalidad tiene una correlación inferior a -0.5 con Natalidad, TasaActividad y TasaParo y una correlación inferior a -0.25 con Industria, TVF, Poblacion, CTH, Infor, Construccion, NumEmpresas, Ocupados, AFS, APT y PIB (al tenerla con una ya podíamos suponer que iba a tenerla con el resto, ya que sabiamos que estas variables presentaban una fuerte correlación entre sí). También observamos una correlación negativa, ya mayor que -0.25 entre TasaParo y CTH, Infor, Construccion, NumEmpresas, Ocupados, AFS, APT y PIB y algunos otros pares como Natalidad y CANE, Natalidad y VS, TasaParo y CANE, CANE con Infor y APT y PIB. Por último hay que resaltar una correlación negativa y fuerte (inferior a -0.5) entre Natalidad e IPC y aún mayor entre TasaParo e IPC. Vistas estas correlaciones tenemos una aproximación razonable de como se agrupan nuestras variables según su comportamiento.
3 Selección de componentes
Ahora voy a realizar directamente un análisis de componentes principales sobre el set de datos. Lo haré sobre los datos escalados a la unidad y miraré todos los autovectores asociados a la matriz de correlaciones, para posteriormente quedarme con aquellos que explican la mayor parte de la variabilidad del set.
par(mfrow=c(1,2), mar=c(2.1, 3.1, 2.1, .5))
fit<- PCA(datos, scale.unit=TRUE, ncp=7, graph=TRUE)
En el mapa de los individuos en las dos dimensiones principales del análisis podemos ver como Madrid, Barcelona, Valencia y Alicante aparecen bastante separados del resto mostrando valores muy altos de la primera componente. Separados por la segunda componente vemos a Ceuta y Melilla por arriba y a Zamora y Ourense por debajo. El resto de provincias aparecen bastante aglutinadas en torno al cero.
Si miramos ahora el mapa de las variables vemos que muchas de ellas están situadas en torno al 1 de la primera dimesión, luego toda su información debe estar recogida de forma bastante precisa en esta componente, por supuesto estas variables son: VS, Industria, TVF, Poblacion, CTH, Infor, Construccion, NumEmpresas, Ocupados, AFS, APT y PIB, aquellas que habiamos visto que presentaban una correlación positiva acentuada en el gráfico de las correlaciones. Situadas muy próximas a la segunda dimensión observamos TasaParo y Natalidad, también lo hace Cane pero muy cerca del cero, lo que implica que su varianza esté poco explicada (por estas dos dimensiones). TasaActividad, IPC y Mortalidad aparecen como una combinación equlibrada de ambas componentes.
par(mfrow=c(1,1), mar=c(3.1, 3.1, 3.1, 3.1))
as.data.frame(fit$eig)fviz_eig(fit, addlabels=TRUE,
main="Variabilidad explicada", ncp=ncol(datos),
ggtheme = theme_light())
Mirando ahora la tabla y el gráfico con los autovalores y el porcentage de variabilidad explicado por cada componente vemos que tomando las 7 primeras componentes queda explicado un 98,7% y que a partir de la séptima componente el resto explican menos del 0.7% de la variabilidad. Así que con 7 componentes conservamos la mayor parte de la información presente en los datos. Siendo un poco menos ambiciosos creo que sería suficiente con emplear 4 componentes, que explican un 92,2% de forma que añadiendo otras 3 tan sólo hemos mejorado un 6,5%, y casi hemos duplicado el número de componentes.
4 PCA - 4
Tomando ya 4 componentes vamos a proceder al análisis de los datos. Durante todo el proceso me referiré a las variables tratándose ya de las variables estandarizadas, es decir, con media 0 y desviación típica 1: \(X_i^* = \tfrac{X_i - \bar{X}}{\sqrt{S(X)}}\).
4.1 Coeficientes de las componentes
En primer lugar miramos la expresión de las componentes en función de las variables originales. Para ello nos fijamos en la matriz \(P_R\), sabiendo que: \(R = P_R\ \Lambda\ P_R^{-1}\) siendo \(\Lambda\) la matriz de autovalores. Luego \(P_R\) es la matriz de cambio de base entre el espacio de autovectores (o componentes) y el de las variables originales.
fit<- PCA(datos, scale.unit=TRUE, ncp=4, graph=FALSE)
P<- as.data.frame(fit$svd$V)
rownames(P)<- colnames(datos)
colnames(P)<- c("dim.1", "dim.2", "dim.3", "dim.4")
PPodemos así decir que la Componente Principal 1 (dim.1) es:
0,29 \(\cdot\) Poblacion + -0,11 \(\cdot\) Mortalidad + 0,04 \(\cdot\) Natalidad + 0,11 \(\cdot\) IPC + 0,29 \(\cdot\) NumEmpresas + 0,29 \(\cdot\) Industria + 0,29 \(\cdot\) Construccion + 0,29 \(\cdot\) CTH + 0,28 \(\cdot\) Infor + 0,29 \(\cdot\) AFS + 0,29 \(\cdot\) APT + 0,11 \(\cdot\) TasaActividad + -0,01 \(\cdot\) TasaParo + 0,29 \(\cdot\) Ocupados + 0,29 \(\cdot\) PIB + 0,02 \(\cdot\) CANE + 0,29 \(\cdot\) TVF + 0,17 \(\cdot\) VS.
Y así podemos ver que las componentes de los autovectores en función de las variables son las columnas de la matriz \(P_R\). Ya llama la atención el elevado número de variables cuya componente para el primer autovector es ~0.29. Todas estás variables deben ser las que veíamos en el mapa de variables muy similares a la primera dimensión y en el siguiente apartado deberían presentar una correlación alta.
4.2 Correlaciones componente - variable
Ahora vamos a ver las correlacones entre las variables y las componentes, estas correlaciones serán proporcionales a la componente del autovector en el espacio de las variables según la expresión. \[ Corr(dim_i, X_j) = \frac{Cov(dim_i, X_j)}{\sqrt{V(dim_i) \cdot V(X_j)}} = e_{ij} \cdot \frac{\sqrt{\lambda_i}}{s_j} \]
Como nuestras variables están estandarizadas (\(s_j = 1\)) la expresión queda: \(Corr(dim_i, X_j) = e_{ij} \cdot \sqrt{\lambda_i}\) luego la correlación entre Población y la primera componente es: 11,47\(^{1/2}\ \cdot\ \) 0,29 = 0,99. Haciendo esto para todas las variables:
A la vista de la tabla podemos decir que las variables con mayor correlación son:
Para la primera dimensión (correlación superior a 0.95): Poblacion, NumEmpresas, Industria, Construccion, CTH, Infor, AFS, APT, Ocupados, PIB y TVF.
Para la segunda dimension (correlación superior a 0.7): Natalidad y TasaParo.
Para la tercera dimensión (correlación superior a 0.4): CANE y TasaParo.
Para la cuarta dimensión (correlación superior a 0.4): IPC, TasaActividad y VS.
Vemos que al ir avanzando en las dimensiones hemos tenido que ir rebajando el criterio sobre la correlación para tener variables que lo superasen, esto se debe a que tomamos las dimensiones ordenadas según la variabilidad que explican y en consecuencia las primeras contienen información de más variables.
4.3 Variables en los planos de las componentes
Vamos ahora a ver como se distribuyen las variables en los planos formados por las nuevas componentes. Como tenemos 4 dimensiones tenemos 6 planos posibles: {1,2}, {1,3}, {1,4}, {2,3}, {2,4}, {3,4}.
colfunc<-colorRampPalette(c("firebrick","goldenrod","darkolivegreen3",
"cornflowerblue", "darkviolet"))
p1<- fviz_pca_var(fit, axes = c(1, 2), col.var = "cos2", repel = TRUE,
gradient.cols = colfunc(6), label="all",
title="Plano 1 - 2") + theme(legend.position="none")
p2<- fviz_pca_var(fit, axes = c(1, 3), col.var = "cos2", repel = TRUE,
gradient.cols = colfunc(6), label="all",
title="Plano 1 - 3") + theme(legend.position="none")
p3<- fviz_pca_var(fit, axes = c(1, 4), col.var = "cos2", repel = TRUE,
gradient.cols = colfunc(6), label="all",
title="Plano 1 - 4") + theme(legend.position="none")
p4<- fviz_pca_var(fit, axes = c(2, 3), col.var = "cos2", repel = TRUE,
gradient.cols = colfunc(6), label="all",
title="Plano 2 - 3") + theme(legend.position="none")
p5<- fviz_pca_var(fit, axes = c(2, 4), col.var = "cos2", repel = TRUE,
gradient.cols = colfunc(6), label="all",
title="Plano 2 - 4") + theme(legend.position="none")
p6<- fviz_pca_var(fit, axes = c(3, 4), col.var = "cos2", repel = TRUE,
gradient.cols = colfunc(6), label="all",
title="Plano 3 - 4") + theme(legend.position="none")
gridExtra::grid.arrange(p1, p2, p3, p4, p5, p6, ncol=2, nrow=3)
A la vista de los diagramas podemos valorar de forma aproximada a que variables equivale cada componente e inferir una descripción de ellas. Como las coordenadas de la variable son la correlación entre la variable y la componente estos resultados se solapan con los del apartado anterior, luego podemos decir que:
La primera componente recoge la información de Poblacion, NumEmpresas, Industria, Construccion, CTH, Infor, AFS, APT, Ocupados, PIB y TVF por lo que podría considerarse un indicador económico.
La segunda variable recoge información de Natalidad y Mortalidad por lo que podemos considerarla un indicador demográfico (obviando la TasaParo que está en una posición intermedia).
La tercera componente recoge la información de CANE, TasaActividad y TasaParo por lo que podríamos valorarla como un indicador de empleo.
La cuarta variable recoge información de VS e IPC por lo que la podemos ver como un indicador de consumo.
Estas interpretaciones de las componentes son bastante libres y están forzadas por la intención de repartir las variables entre las componentes, no obstante las he realizado observando los gráficos e intentando mantener la coherencia con el resultado anterior. Además la interpretación es incompleta por que desconozco (y no he encontrado) el qué marcan algunas variables referidas mediante siglas como APT, TVF, AFS, CANE, …
4.4 Proporción de la varianza explicada por variable
Si la correlación entre una dimensión y una variable era la raiz del autovalor asociado a la dimensión por la proyección de la componente sobre dicha variable, si elevamos este valor al cuadrado tenemos la proporción de la variamza de cada variable explicada por cada dimensión, o el coseno del ángulo formado entre la variable y la dimensión al cuadrado.
var.cos2<- as.data.frame(fit$var$cos2)
var.cos2$suma<- rowSums(var.cos2)
var.cos2<- var.cos2[order(var.cos2$suma, decreasing = T),]
var.cos2corrplot(as.matrix(var.cos2), is.corr=FALSE, cl.ratio = 0.5,
method="color", cl.align.text = "l", cl.lim = c(0,1),
col=brewer.pal(n=10, name="RdYlBu"))
Mirando la columna suma vemos la proporión de variabilidad explicada en total por las cuatro dimensiones. Así podemos ver que todas las variables tienen explicado al menos el 75% de su varianza. Podemos señalar las peor explicadas como aquellas cuya proporción de varianza explicada es menor al 80% que son: IPC, Natalidad y VS.
4.5 Contribución a la varianza de cada variable
Podemos ahora calcular que porccentaje de la varianza de cada componente es debida a la varianza de cada variable. Para ello bastará con dividir el porcentaje de varianza explicado de la variable entre la suma de los pocentajes de varianza explicada para todas las variables (que será el autovalor asociado a la componente). Si seguimos los cálculos a partir de la matriz de correlaciones: \[ Corr(Dim_i, X_j) = e_{ij} \cdot \sqrt{\lambda_i}\\ VarExp(X_j, Dim_i) = e_{ij}^2 \cdot \lambda_i\\ ContrVar(Dim_i, X_j) = \frac{e_{ij}^2 \cdot \lambda_i}{\sum_j (e_{ij}^2 \cdot \lambda_i)} = \frac{e_{ij}^2 \cdot \lambda_i}{\lambda_i \cdot \sum_j e_{ij}^2} = \frac{VarExp(X_j, Dim_i)}{\lambda_i} \] Donde he usado en la última línea que: \(\sum_j e_{ij}^2 = 1\), es decir, que los autovectores están normalizados. De esta forma podemos ver la contribución a la varianza de cada variable como la parte de la variabilidad total explicada por esa dimensión debida a esa variable (ya que la variabilidad total explicada por cada dimensión es proporcional a su autovalor). En nuestro caso los valores son:
as.data.frame(fit$var$contrib)corrplot(as.matrix(fit$var$contrib), is.corr=FALSE, cl.ratio = 0.5,
method="color", cl.align.text = "l", cl.lim = c(0,50))
p1<- fviz_contrib(fit,choice="var",axes=1,top=12,ggtheme=theme_light())
p2<- fviz_contrib(fit,choice="var",axes=2,top=12,ggtheme=theme_light())
p3<- fviz_contrib(fit,choice="var",axes=3,top=12,ggtheme=theme_light())
p4<- fviz_contrib(fit,choice="var",axes=4,top=12,ggtheme=theme_light())
gridExtra::grid.arrange(p1, p2, p3, p4, ncol=2, nrow=2)
Podemos ver como la primera dimensión tiene 11 variables que aportan más de la media (marcada en línea roja discontínua), mientras que la segunda tiene 5, la tercera 6 y la cuarta 4. Debo señalar que la media es igual para todas ya que la contribución de todas las variables suma siempre el 100% (de la variabilidad explicada por esa dimensión).
4.6 Observaciones en los planos de las componentes
Es ya momento de ver como quedan los individuos en los planos de las componentes.
p1<- fviz_pca_ind(fit, axes = c(1, 2), col.ind = 'cos2', repel = TRUE,
gradient.cols = colfunc(6), label="all",
title="Plano 1 - 2") + theme(legend.position="none")
p2<- fviz_pca_ind(fit, axes = c(1, 3), col.ind = "cos2", repel = TRUE,
gradient.cols = colfunc(6), label="all",
title="Plano 1 - 3") + theme(legend.position="none")
p3<- fviz_pca_ind(fit, axes = c(1, 4), col.ind = "cos2", repel = TRUE,
gradient.cols = colfunc(6), label="all",
title="Plano 1 - 4") + theme(legend.position="none")
p4<- fviz_pca_ind(fit, axes = c(2, 3), col.ind = "cos2", repel = TRUE,
gradient.cols = colfunc(6), label="all",
title="Plano 2 - 3") + theme(legend.position="none")
p5<- fviz_pca_ind(fit, axes = c(2, 4), col.ind = "cos2", repel = TRUE,
gradient.cols = colfunc(6), label="all",
title="Plano 2 - 4") + theme(legend.position="none")
p6<- fviz_pca_ind(fit, axes = c(3, 4), col.ind = "cos2", repel = TRUE,
gradient.cols = colfunc(6), label="all",
title="Plano 3 - 4") + theme(legend.position="none")
gridExtra::grid.arrange(p1, p2, p3, p4, p5, p6, ncol=2, nrow=3)
A la vista de los mapas, dado que tenemos un elevado número de individuos a catalogar, he elavorado el siguiente gráfico, que ayuda a hacerse una idea de forma visual de la situación de cada provincia, dividiendo todo el rango de valores que tomaban las coordenadas para cada dimensión en 4 tramos y ya interpretando las dimensiones como propuse en el apartado 4.3.
ind.coord<- as.data.frame(fit$ind$coord)
haz_factor<- function (vect, n, factor) {
aux<- vect
aux2<- c(rep(0, n-1), 0.001)
M<- max(vect)
m<- min(vect)
piece<- (M - m) / n
for (i in seq(1,n)){
aux<- replace(aux, which(m+(i-1)*piece <= vect &
vect < m+(i+aux2[i])*piece), (i-1)/(n-1))
}
if(factor){
aux<- as.factor(aux)
}
return(aux)
}
interpretacion<- data.frame(row.names = row.names(datos))
interpretacion$Economia<- haz_factor(ind.coord$Dim.1, 4, FALSE)
interpretacion$Demografia<- haz_factor(ind.coord$Dim.2, 4, FALSE)
interpretacion$Empleo<- haz_factor(ind.coord$Dim.3, 4, FALSE)
interpretacion$Consumo<- haz_factor(ind.coord$Dim.4, 4, FALSE)
par(mar=c(2,4.1,3,2))
corrplot(t(as.matrix(interpretacion)), is.corr=FALSE, cl.ratio = 0.6,
method="color", cl.align.text = "l", cl.lim = c(0,1),
cl.pos="b", mar=c(1,0.3,0.1,2), col=colfunc(8))
En primer lugar hay que destacar que las dimensiones 2 y 3 están negativamente correlacionadas con Natalidad y TasaActividad de tal forma que valores altos en la dimensión 2 implican una demografía en recesión e igualmente ocurre con el empleo y la dimensión 3. Aclarado esto simplemente echando un ojo sobre nuestro gráfico identíficamos:
Las provincias económicamente más fuertes son Madrid y Barcelona.
Las provincias demográficamente más débiles (con mayor mortalidad) son Ceuta y Melilla.
Las provincias más débiles a nivel de empleo (con mayor tasa de paro) son Jaen y Valencia.
Las provincias más fuertes a nivel de consumo son Alicante, Baleares, Girona, Tarragona y Valencia.
Lo cotejamos con nuestros datos y vemos que, aunque no literalmente, si es cierto que estos individuos toman posiciones elevadas al ordenar los datos según la suma (o la resta) de las variables con mayor correlación indicadas en el apartado 4.3.
4.7 Índice de desarrollo económico
Como ya hemos visto la primera dimensión podría considerarse un indicador económico, aunque también es cierto que destacan en exceso las grandes ciudades, así que he decidido construir 2 índices:
El primer índice se define como la suma pesada (mediante el porcentaje de variabilidad explicada sobre el total explicado) de mis 4 dimensiones, tomando la raíz cuadrada para que el índice cumpla la condición \(\sum_j e_{ij}^2 = 1\) de normalidad e invirtiendo el signo del coeficiente para la segunda y la tercera dimesión para que midan los efectos que considero positivos para el desarrollo económico. Con estos criterios el Índice de Desarrollo Económico sería: \[ IDE_1 = \frac{1}{\sqrt{0.922}} \cdot \left(\sqrt{0.637} \cdot dim1\ -\ \sqrt{0.142} \cdot dim2\ -\ \sqrt{0.091} \cdot dim3\ +\ \sqrt{0.052} \cdot dim4\right) \]
El segundo índice lo he definido de una forma un poco menos rigurosa, simplemente sumando las dimesiones. De nuevo he tomado la raíz para que el índice cumpla la condición de normalidad y he invertido el signo del coeficiente para la segunda y la tercera dimensión. Con estos criterios el Índice de Desarrollo Económico sería: \[ IDE_2 = \frac{1}{\sqrt{4}} \cdot \left(dim1\ -\ dim2\ -\ dim3\ +\ dim4\right) \]
Si sustituimos la expresión de las dimensiones respecto a las variables originales en esa ecuación, la aplicamos sobre las variables estandarizadas, y ordenamos el resultado nos quedarían los siguientes valores:
coef1<- sqrt(fit$eig[1,2]/fit$eig[4,3])
coef2<- -sqrt(fit$eig[2,2]/fit$eig[4,3])
coef3<- -sqrt(fit$eig[3,2]/fit$eig[4,3])
coef4<- sqrt(fit$eig[4,2]/fit$eig[4,3])
coefs<- coef1 * P$dim.1 + coef2 * P$dim.2 + coef3 * P$dim.3 +
coef4 * P$dim.4
escalados<- as.matrix(scale.default(datos))
IDEs<- data.frame(row.names = rownames(datos))
IDEs$IDE<- rowSums(escalados*coefs)
IDEs$aux<- 0
IDEs<- IDEs[order(IDEs$IDE),]
IDEs$aux<- NULL
IDEspar(mar=c(1,6,1,1))
barplot(t(as.matrix(IDEs)), horiz=T, col=colfunc(52),
border = "white", main="IDE por provincia",
cex.names=0.9, las=1, names.arg = rownames(IDEs),
space = 0.02
)
A la vista del gráfico podemos concluir que el IDE\(_1\) no reduce las diferencias de las grandes ciudades con el restovrespecto a usar solo la dimesión 1 pero da unos resultados a priori acordes a la experiencia, con Madrid y Barcelona a la cabeza y Valencia, Alicante, Baleares, Málaga y Sevilla en la zona alta, sin embargo, me llama enormemente la atención lo abajo que aparecen las provincias de Euskadi, además de Navarra y otras provincias del norte en general. Esto me lleva a pensar que como indicador económico tal vez no sea idóneo, pero para indagar más a fondo debería conocer mejor las variables. La expresión asociada es: \[\begin{aligned} IDE_1 &= 0.215Pobl + 0.021Mort - 0.102Nat + 0.420IPC\\ &+ 0.236NumE + 0.246Indu + 0.259Cons + 0.225CTH\\ &+ 0.218Info + 0.215AFS + 0.228APT + 0.189TAct\\ &- 0.367TPar + 0.236Ocup + 0.236PIB - 0.163CANE\\ &+ 0.222TVF + 0.168VS\\ &\text{Todas ellas normalizadas} \end{aligned}\]
Que es una expresión que tiene dos peros, en primer lugar está dando pesos muy similares a varias variables y en segundo penaliza la Natalidad mientras que la Mortalidad suma.
Para este índice los valores pedidos son:
Madrid: 8.700
Melilla: -0.445
Si repetimos el proceso para el segundo IDE\(_2\) resulta:
coef<- sqrt(1/4)
coefs2<- rowSums(coef * P)
IDEs<- data.frame(row.names = rownames(datos))
IDEs$IDE<- rowSums(escalados*coefs2)
IDEs$aux<- 0
IDEs<- IDEs[order(IDEs$IDE),]
IDEs$aux<- NULL
IDEspar(mar=c(1,6,1,1))
barplot(t(as.matrix(IDEs)), horiz=T, col=colfunc(52),
border = "white", main="IDE por provincia",
cex.names=0.9, las=1, names.arg = rownames(IDEs),
space = 0.02
)
Podemos ahora concluir que el IDE\(_2\) tampoco reduce las diferencias de Madrid y Barcelona con el resto pero los resultados a parecen más similares a mi preconcepción, las que estaban a la cabeza se mantienen, salvo tal vez alicante que cae bastante y las provincias vascas ascienden de forma sensible en la lista. Aún así tampoco este me parece un indicador económico fiable. La expresión asociada es: \[\begin{aligned} IDE_2 &= 0.146Pobl - 0.303Mort + 0.078Nat - 0.041IPC\\ &+ 0.103NumE + 0.155Indu + 0.104Cons + 0.151CTH\\ &- 0.023Info + 0.112AFS + 0.046APT + 0.273TAct\\ &- 0.307TPar + 0.110Ocup + 0.042PIB - 0.525CANE\\ &+ 0.217TVF + 0.539VS\\ &\text{Todas ellas normalizadas} \end{aligned}\]
Que también tiene dos peros, pese a haber resuelto la incoherencia del otro índice con la Natalidad y la Mortalidad en este la TasaParo suma y el IPC penaliza lo que parece ser también un contrasentido.
Para este índice los valores pedidos son:
Madrid: 8.070
Melilla: -1.641
5 Técnicas de clustering
Ahora vamos a hacer clustering sobre las provincias, para ello lo primero que hago es normalizar los datos (esta vez manualmente) para que la diferente escala entre las medidas no afecte a la medida de la distancia. Y ver la distancia euclídea entre las provincias en un mapa de calor, así podemos hacernos una idea de como deben agruparse.
datos_esc<- scale(datos)
dist_esc<- dist(datos_esc, method = "euclidean")
colfunc<-colorRampPalette(c("#fdbb2d","#b21f1f","#1a2a6c"))
fviz_dist(dist_esc, gradient = list(low="white",
mid=colfunc(3)[2],
high=colfunc(3)[3]))
El mapa únicamente nos permite corroborar que Madrid y Barcelona están muy alejadas de las demás, aunque tampoco parecen estar próximas entre sí. Por otro lado si da la sensación de que algunas de las provincias que aparecen próximas (en el mapa de calor) si están a menor distancia entre sí.
5.1 Jerárquico
Ahora voy a realizar un agrupamiento jerárquico mediante el método Ward que minimiza la varianza entre los individuos pertenecientes a un grupo. Para elegir el número de grupos voy a ver cómo resultan los agrupamientos de 4 y 6 grupos (ya he ensayado otras opciones) y me acabaré decantando por una de ellos.
res.hc_st <- hclust(dist_esc, method="ward.D2")
grp <- cutree(res.hc_st, k = 6)
p1<- fviz_dend(res.hc_st, k = 6, # ut in four groups
cex = 0.5, # label size
k_colors = colfunc(6),
color_labels_by_k = TRUE, # color labels by groups
rect = TRUE, # Add rectangle around groups
main="Dendograma. 6 grupos",
horiz = TRUE
)
p2<- fviz_cluster(list(data = datos_esc, cluster = grp),
palette = colfunc(6),
ellipse.type = "convex", # Concentration ellipse
repel = TRUE, # Avoid label overplotting (slow)
show.clust.cent = FALSE, ggtheme = theme_minimal(),
main="Cluster-jerárquico. 6 grupos")
grp <- cutree(res.hc_st, k = 4)
p3<- fviz_dend(res.hc_st, k = 4, # ut in four groups
cex = 0.5, # label size
k_colors = colfunc(4),
color_labels_by_k = TRUE, # color labels by groups
rect = TRUE, # Add rectangle around groups
main="Dendograma. 4 grupos",
horiz = TRUE
)
p4<- fviz_cluster(list(data = datos_esc, cluster = grp),
palette = colfunc(4),
ellipse.type = "convex", # Concentration ellipse
repel = TRUE, # Avoid label overplotting (slow)
show.clust.cent = FALSE, ggtheme = theme_minimal(),
main="Cluster-jerárquico. 4 grupos"
)
gridExtra::grid.arrange(p1, p3, p2, p4, ncol=2, nrow=2)
La diferencia consiste esencialmente en separar un bloque inmenso con todas las provincias ‘promedio’ en dos y separar Ceuta y Melilla de otro bloque. Ambas cosas me parecen positivas para el clustering ya que parece lógico que Ceuta y Melilla vayan juntas y por su cuenta y que se divida el grupo más numeroso del clustering de 4 grupos en 2 más o menos igual de poblados. El gráfico de los grupos en las dos dimensiones principales corroboran la bondad del clustering de 6 grupos, Ceuta y Melilla aparecen muy apartadas y aunque hay 3 grupos que parecen solaparse en exceso hay que recordar que es la proyección en dos dimensiones de un set de datos que tiene 18.
5.2 Proximidad
Lo primero que necesito para hacer clustering por proximidad es ver cúal es el número adecuado de clústers, para ello dispongo fundamentalmente de 2 criterios, el criterio de Elbow que consiste en fijarse en el momento en el que el aumento del número de grupos involucrados ya no reduce la variabilidad intra-clúster: \[ W(k) = \sum_k\left(\sum_{j=1}^p \sum_{i\in C_K} (x_{ij} - \bar{x}_{jK}) \right)\\ \text{siendo p el número de variables, i las observaciones,} \\ \text{k el número de clústeres y K refiere un clúster concreto} \]
Bucamos \(k\) tal que \(W(k+1) \approx W(k)\) para ello representamos \(W(k)\) en función de \(k\).
Por otro lado tenemos el criterio Silhouette que consiste en maximizar la media de la diferencia entre la distancia de una observación a las observaciones externas a su clúster y la distancia de esa observación a las observaciones internas, dividido entre el máximo de ambas. \[ \bar{s}(k) = \frac{\sum_{i=1}^n\left(\frac{b_i(k) - a_i(k)}{\max(b_i,\ a_i)} \right)}{n} \] Bucamos \(k\) tal que \(s(k) = \max(s(k))\) y para ello representaré \(s(k)\) en función de \(k\). Si observamos ambos criterios.
fviz_nbclust(datos_esc, kmeans, method = "wss") +
geom_vline(xintercept = 3, linetype = 2) +
geom_vline(xintercept = 6, linetype = 4) +
labs(subtitle = "Elbow method")
fviz_nbclust(datos_esc, kmeans, method = "silhouette")+
labs(subtitle = "Silhouette method")
Vemos que con el criterio Elbow, a priori podríamos pensar en 6 grupos, ya que entre 5 y 6 la mejora es sensible, no obstante si miramos la diferencia entre 4 y 5 es mínima y la diferencia entre 4 y 6 no parece suficiente como para justificar el uso de 2 clústers más. Por su parte el criterio Silhouette indica 3 ó 5. Con la curva Elbow vemos que 4 no representan una gran mejora sobre 3 y tampoco lo hacen 5 sobre 4, es por esto que voy a mirar la calidad de estas 4 posibilidades para poder decantarme.
Para comprobar la calidad de los clústers representamos \(s(i)\) que es el valor que encontramos dentro del sumatorio de la expresión de \(s(k)\), y los clústeres sobre las observaciones proyectadas en las dimensiones principales.
km.res <- kmeans(datos_esc, 3)
sil <- silhouette(km.res$cluster, dist_esc)
rownames(sil) <- rownames(datos)
p1<- fviz_silhouette(sil, legend="none", font.main=9,
print.summary = FALSE)
p2<- fviz_cluster(km.res, datos_esc, labelsize = 0,
main="3-Cluster Plot")
gridExtra::grid.arrange(p1, p2, ncol=2, nrow=1)
km.res.sel <- kmeans(datos_esc, 4)
sil <- silhouette(km.res.sel$cluster, dist_esc)
rownames(sil) <- rownames(datos)
p1<- fviz_silhouette(sil, legend="none", font.main=9,
print.summary = FALSE)
p2<- fviz_cluster(km.res.sel, datos_esc, labelsize = 0,
main="4-Cluster Plot")
gridExtra::grid.arrange(p1, p2, ncol=2, nrow=1)
km.res <- kmeans(datos_esc, 5)
sil <- silhouette(km.res$cluster, dist_esc)
rownames(sil) <- rownames(datos)
p1<- fviz_silhouette(sil, legend="none", font.main=9,
print.summary = FALSE)
p2<- fviz_cluster(km.res, datos_esc, labelsize = 0,
main="5-Cluster Plot")
gridExtra::grid.arrange(p1, p2, ncol=2, nrow=1)
km.res <- kmeans(datos_esc, 6)
sil <- silhouette(km.res$cluster, dist_esc)
p1<- fviz_silhouette(sil, legend="none", font.main=9,
print.summary = FALSE)
p2<- fviz_cluster(km.res, datos_esc, labelsize = 0,
main="6-Cluster Plot")
gridExtra::grid.arrange(p1, p2, ncol=2, nrow=1)
Lo primero que me llama la atención es que el valor promedio del cociente de distancias no se corresponde con el observado en la curva \(s(k)\), los valores mayores se presentan ahora para 3 y 4 clústeres y el valor para 5 disminuye bastante (he revisado el código sin encontrar la posible fuente de discrepancia). Viendo además las 4 opciones notamos que para 5 clústers tan solo hay dos valores negativo, pero más allá de los cocientes de distancias, mirando los grupos, el que más me convence es el de 4 clústers, ya que no existe solapamiento entre clústeres, el criterio Elbow lo señalaba como un valor acertado y creo que es una cantidad de categorías razonable para las 52 observaciones.
5.3 Interpretación socioeconómica de los clústers
Como este apartado ya es puramente personal voy a seleccionar las variable que me sirvan a mi para tener una idea de la situación socioeconómica, ya que hya siglas cuyo significado no he sido capaz de encontrar. por este motivo voy a fijarme en las que comprendo que son: Población, Mortalidad, Natalidad, IPC, NumEmpresas, Industria, Construcción, TasaActividad, TasaParo, Ocupados y PIB. Como 11 es un número desagradable y todo mi análisis va a ser comparativo voy a hacerlo sobre los datos normalizados e incorporaré como variable numero 12 el índice construido en el apartado 4.7 IDE\(_1\), de esta forma además conservo algo de información de las variables que estoy eliminando.
datos_fin<- datos[,c("Poblacion", "Mortalidad", "Natalidad",
"IPC", "NumEmpresas", "Industria",
"Construccion", "TasaActividad",
"TasaParo", "Ocupados", "PIB")]
datos_fin_esc<- as.data.frame(scale(datos_fin))
datos_fin_esc$IDE<- rowSums(escalados*coefs)
# Añado la agrupación en 4 clústers realizada sobre todo el set
datos_fin_esc$grupo<- as.factor(km.res.sel$cluster)
p1<- ggplot(datos_fin_esc, aes(x=Poblacion, fill=grupo)) +
geom_density(alpha=.3, show.legend = FALSE) +
theme(axis.title.x=element_blank()) + ggtitle("Población")
p2<- ggplot(datos_fin_esc, aes(x=Mortalidad, fill=grupo)) +
geom_density(alpha=.3, show.legend = FALSE) +
theme(axis.title.x=element_blank()) + ggtitle("Mortalidad")
p3<- ggplot(datos_fin_esc, aes(x=Natalidad, fill=grupo)) +
geom_density(alpha=.3) +
theme(axis.title.x=element_blank()) + ggtitle("Natalidad")
p4<- ggplot(datos_fin_esc, aes(x=IPC, fill=grupo)) +
geom_density(alpha=.3, show.legend = FALSE) +
theme(axis.title.x=element_blank()) + ggtitle("IPC")
p5<- ggplot(datos_fin_esc, aes(x=NumEmpresas, fill=grupo)) +
geom_density(alpha=.3, show.legend = FALSE) +
theme(axis.title.x=element_blank()) + ggtitle("NumEmpresas")
p6<- ggplot(datos_fin_esc, aes(x=Industria, fill=grupo)) +
geom_density(alpha=.3) +
theme(axis.title.x=element_blank()) + ggtitle("Industria")
p7<- ggplot(datos_fin_esc, aes(x=Construccion, fill=grupo)) +
geom_density(alpha=.3, show.legend = FALSE) +
theme(axis.title.x=element_blank()) + ggtitle("Construccion")
p8<- ggplot(datos_fin_esc, aes(x=TasaActividad, fill=grupo)) +
geom_density(alpha=.3, show.legend = FALSE) +
theme(axis.title.x=element_blank()) + ggtitle("TasaActividad")
p9<- ggplot(datos_fin_esc, aes(x=TasaParo, fill=grupo)) +
geom_density(alpha=.3) +
theme(axis.title.x=element_blank()) + ggtitle("TasaParo")
p10<- ggplot(datos_fin_esc, aes(x=Ocupados, fill=grupo)) +
geom_density(alpha=.3, show.legend = FALSE) +
theme(axis.title.x=element_blank()) + ggtitle("Ocupados")
p11<- ggplot(datos_fin_esc, aes(x=PIB, fill=grupo)) +
geom_density(alpha=.3, show.legend = FALSE) +
theme(axis.title.x=element_blank()) + ggtitle("PIB")
p12<- ggplot(datos_fin_esc, aes(x=IDE, fill=grupo)) +
geom_density(alpha=.3) +
theme(axis.title.x=element_blank()) + ggtitle("IDE")
gridExtra::grid.arrange(p1, p2, p3, p4, p5, p6, p7, p8, p9,
p10, p11, p12, ncol=3, nrow=4)
grupos<- rep(datos_fin_esc$grupo, 12)
valores<- c()
variables<- c()
for (col in names(datos_fin_esc)[seq(1:12)]){
valores<- c(valores, datos_fin_esc[[col]])
variables<- c(variables, rep(col, 52))
}
to_box_plot<- data.frame(grupos, variables, valores)
ggplot(to_box_plot, aes(x=variables, y=valores, fill=grupos)) +
geom_boxplot() +
facet_wrap(~variables, ncol=3, scale="free") +
theme(axis.text.x = element_blank(),
axis.ticks = element_blank(),
legend.title = element_blank())
Para hacer una correcta interpretación de los gráficos primero debemos recordar que el clúster número 1 tiene tan sólo 2, el número dos tiene 16, el número 3 8 y el 4 26, así que estas medidas de dispersión en el primer clúster no tienen mucho sentido aunque si nos son útiles para comparar con el resto. Por otro lado es evidente que la información de los dos gráficos es redundante, no obstante, considero que se apoyan entre sí, y uno puede servir para despejar dudas que despierte el otro.
Sin mirar en detalle ya podemos notar algunas cuestiones: las variables Natalidad e IPC apenas nos permiten discriminar entre clústeres, la variable TasaParo discrimina muy bien entre clústeres, la distribución del primer clúster suele estar alejada de las distribuciones del resto y las distribuciones del segundo y cuarto clúster suelen tener amplias regiones de solapamiento.
Si miramos ahora los clústers individualmente, pero en comparación con el resto:
Las provincias del Clúster 1 (Barcelona y Madrid) tienen Población, NumEmpresas, Industria, Construcción, Ocupados y PIB muy superiores a las demás; TasaActividad, IPC E IDE elevadas; una Natalidad promedio y Mortalidad y TasaParo bajas luego podemos decir que son provincias socioeconómicamente muy desarrolladas.
Las provincias del Clúster 2 (Albacete, Almería, Badajoz, Castellón, Ceuta, CiudadReal, Cádiz, Córdoba, Granada, Guadalajara, Huelva, Jaén, Melilla, Palmas, SantaCruz y Toledo) tienen TasaParo elevada; Natalidad, Mortalidad y TasaActividad promedio y Poblacion, IPC, NumEmpresas, Industria, Construccion, Ocupados, PIB e IDE bajas luego se trata de provincias con bajo desarrollo socioeconómico.
Las provincias del Clúster 3 (Alicante, Balears, Girona, Murcia, Málaga, Sevilla, Tarragona y Valencia) ocupan las regiones intermedias en todas las variables, luego podemos decir que se trata de provincias con un desarrollo socioeconómico moderado.
Las provincias del Clúster 4 (Álava, Asturias, Bizkaia, Burgos, Cantabria, Coruña, Cuenca, Cáceres, Gipuzkoa, Huesca, León, Lleida, Lugo, Navarra, Ourense, Palencia, Pontevedra, Rioja, Salamanca, Segovia, Soria, Teruel, Valladolid, Zamora, Zaragoza y Ávila) tienen Mortalidad alta; TasaActividad e IPC promedio y Poblacion, Natalidad, NumEmpresas, Industria, Construccion, TasaParo, Ocupados, PIB e IDE bajas luego son provincias de bajo desarrollo socioeconómico.
6 Muchas gracias
Ha sido un ejercicio tremendamente útil para asentar los conocimientos adquiridos en clase, muy adecuadamente enfocado.