M1U1-METD

suppressMessages(suppressWarnings(library(cluster)))
suppressMessages(suppressWarnings(library(knitr)))
suppressMessages(suppressWarnings(library(rmarkdown)))
suppressMessages(suppressWarnings(library(paqueteMETODOS)))
suppressMessages(suppressWarnings(library(naniar)))
suppressMessages(suppressWarnings(library(mice)))
suppressMessages(suppressWarnings(library(ggmice)))
suppressMessages(suppressWarnings(library(pastecs)))
suppressMessages(suppressWarnings(library(stringr)))
suppressMessages(suppressWarnings(library(stringi)))
suppressMessages(suppressWarnings(library(dplyr)))
suppressMessages(suppressWarnings(library(tidyverse)))
suppressMessages(suppressWarnings(library(kableExtra)))
suppressMessages(suppressWarnings(library(reticulate)))
suppressMessages(suppressWarnings(library(lsr)))
suppressMessages(suppressWarnings(library(tibble)))
suppressMessages(suppressWarnings(library(patchwork)))
suppressMessages(suppressWarnings(library(sf)))
suppressMessages(suppressWarnings(library(osmdata)))
suppressMessages(suppressWarnings(library(mapview)))
suppressMessages(suppressWarnings(library(leaflet)))
suppressMessages(suppressWarnings(library(heatmaply)))
suppressMessages(suppressWarnings(library(ggplot2)))
suppressMessages(suppressWarnings(library(stats)))
suppressMessages(suppressWarnings(library(e1071)))
suppressMessages(suppressWarnings(library(grid)))
suppressMessages(suppressWarnings(library(gridExtra)))
suppressMessages(suppressWarnings(library(moments)))
suppressMessages(suppressWarnings(library(diptest)))
suppressMessages(suppressWarnings(library(LaplacesDemon)))
suppressMessages(suppressWarnings(library(ggExtra)))
suppressMessages(suppressWarnings(library(RColorBrewer)))
suppressMessages(suppressWarnings(library(hexbin)))
suppressMessages(suppressWarnings(library(MASS)))
suppressMessages(suppressWarnings(library(ggfortify)))
suppressMessages(suppressWarnings(library(latex2exp)))
suppressMessages(suppressWarnings(library(lmtest)))
suppressMessages(suppressWarnings(library(interactions)))
suppressMessages(suppressWarnings(library(leafpop)))
suppressMessages(suppressWarnings(library(corrplot)))
suppressMessages(suppressWarnings(library(factoextra)))
suppressMessages(suppressWarnings(library(plotly)))
suppressMessages(suppressWarnings(library(FactoMineR)))
suppressMessages(suppressWarnings(library(gplots)))
data(vivienda)

Problema

Una empresa inmobiliaria líder en una gran ciudad está buscando comprender en profundidad el mercado de viviendas urbanas para tomar decisiones estratégicas más informadas. La empresa posee una base de datos extensa que contiene información detallada sobre diversas propiedades residenciales disponibles en el mercado. Se requiere realizar un análisis holístico de estos datos para identificar patrones, relaciones y segmentaciones relevantes que permitan mejorar la toma de decisiones en cuanto a la compra, venta y valoración de propiedades.

Solución

Análisis Exploratorio

Valores únicos

unique_counts <- sapply(vivienda, function(x) length(unique(x[!is.na(x) & x != "" & !is.nan(x)])))

unique_counts_df <- data.frame('Conteo' = unique_counts)

rownames(unique_counts_df) <- names(unique_counts)

paged_table(unique_counts_df, options = list(rows.print = 15)) %>%
  kable(caption = 'Valores únicos por Variable del Dataset Original')

Valores únicos por Variable del Dataset Original

	Conteo
id	8319
zona	5
piso	12
estrato	4
preciom	539
areaconst	652
parqueaderos	10
banios	11
habitaciones	11
tipo	2
barrio	436
longitud	2928
latitud	3679

 

 

Variables numéricas

df = subset(vivienda, select = -c(id, zona, piso, tipo, barrio))
options(scipen = 999)
summary_df <- stat.desc(df) 
summary_df <- round(summary_df, 2)

paged_table(summary_df, options = list(rows.print = 15)) %>%
  kable(caption = 'Variables Numéricas')

Variables Numéricas

	estrato	preciom	areaconst	parqueaderos	banios	habitaciones	longitud	latitud
nbr.val	8319.00	8320.00	8319.00	6717.00	8319.00	8319.00	8319.00	8319.00
nbr.null	0.00	0.00	0.00	0.00	45.00	66.00	0.00	0.00
nbr.na	3.00	2.00	3.00	1605.00	3.00	3.00	3.00	3.00
min	3.00	58.00	30.00	1.00	0.00	0.00	-76.59	3.33
max	6.00	1999.00	1745.00	10.00	10.00	10.00	-76.46	3.50
range	3.00	1941.00	1715.00	9.00	10.00	10.00	0.13	0.16
sum	38547.00	3609981.00	1455283.75	12327.00	25883.00	29993.00	-636641.47	28431.38
median	5.00	330.00	123.00	2.00	3.00	3.00	-76.53	3.42
mean	4.63	433.89	174.93	1.84	3.11	3.61	-76.53	3.42
SE.mean	0.01	3.60	1.57	0.01	0.02	0.02	0.00	0.00
CI.mean.0.95	0.02	7.06	3.07	0.03	0.03	0.03	0.00	0.00
var	1.06	108009.01	20438.74	1.27	2.04	2.13	0.00	0.00
std.dev	1.03	328.65	142.96	1.12	1.43	1.46	0.02	0.04
coef.var	0.22	0.76	0.82	0.61	0.46	0.40	0.00	0.01

 

 

Valores vacíos

# Eliminar valores faltantes comunes a todas las columnas
vivienda <- vivienda %>% drop_na(tipo)

p1 <- vis_miss(vivienda)  + 
  labs(title = 'Mapa de calor para valores faltantes') + 
  theme(plot.title = element_text(color = 'black', face = 'bold', size = 12, hjust = 0.5))
p2 <- plot_pattern(vivienda, square = TRUE, rotate = TRUE) +
  theme(
    legend.position = "none",
    axis.title = element_blank(),
    axis.title.x.top = element_blank(),
    axis.title.y.right = element_blank()
  )+ 
  labs(title = 'Patrón de valores faltantes') + 
  theme(plot.title = element_text(color = 'black', face = 'bold', size = 12, hjust = 0.5))

p1 | p2

Limpieza de datos

data(vivienda)
# Función para limpieza de datos
clean_vivienda_data <- function(df) {

  df$tipo <- tolower(df$tipo)
  df$barrio <- tolower(df$barrio)
  original_rows <- nrow(df)
  

  df$tipo <- gsub("apto", "apartamento", df$tipo)
  

  df$barrio <- iconv(df$barrio, "UTF-8", "ASCII//TRANSLIT")
  barrio_replacements <- c(
    'urbanizacion la flora' = 'la flora',
    'caney' = 'el caney',
    'cristales' = 'los cristales',
    'alf?crez real' = 'alferez real',
    'parcelaciones pance' = 'pance',
    'juanamb??' = 'juanambu',
    'santa monica residencial' = 'santa monica',
    'mel?cndez' = 'melendez',
    'el aguacatal' = 'aguacatal',
    'bajo aguacatal' = 'aguacatal',
    'miradol del aguacatal' = 'aguacatal',
    'sector aguacatal' = 'aguacatal',
    'arboleda campestre candelaria' = 'arboledas',
    'arboleda' = 'arboledas'
  )


  df$barrio <- sapply(df$barrio, function(b) {
    if (b %in% names(barrio_replacements)) {
      barrio_replacements[[b]]
    } else {
      b
    }
  })
  
  df <- df[df$latitud >= 3.35 & df$latitud <= 3.55 & 
           df$longitud >= -76.60 & df$longitud <= -76.45, ]
  
  df <- subset(df, banios != 0 & habitaciones != 0)
  df <- df[!duplicated(df$id), ]
  df <- df[, !(names(df) %in% "piso")]
  df <- df[!is.na(df$parqueaderos), ]
  df$price_per_sqm <- df$preciom / df$areaconst
  filtered_rows <- nrow(df)
  difference <- original_rows - filtered_rows
  percentage_loss <- (difference / original_rows) * 100
  list(cleaned_data = df, information_loss = percentage_loss)
}


result <- clean_vivienda_data(vivienda)
df <- result$cleaned_data
loss_percentage <- result$information_loss
cat("El porcentaje de pérdida de información respecto a filas, al eliminar la columna piso y los valores NaN de parqueaderos es de:", round(loss_percentage,2), "% respecto al conjunto de datos inicial.\n")

El porcentaje de pérdida de información respecto a filas, al eliminar la columna piso y los valores NaN de parqueaderos es de: 23.84 % respecto al conjunto de datos inicial.

 

 

---

Análisis de Componentes Principales

 

Matriz de covarianza

• estrato (X1) y preciom (X2) tienen una correlación positiva del 0.67.
• estrato (X1) y habitaciones (X6) tienen una leve correlación negativa de -0.39.
• preciom (X2) y parqueaderos (X4) están correlacionados positivamente (0.57).
• areaconst (X3) y price_per_sqm (X9) muestran un correlación negativa débil (-0.38).

num_df = subset(df, select = -c(id, zona, tipo, barrio))
names(num_df) = c("X1","X2","X3","X4", "X5", "X6","X7", "X8", "X9")
M = cor(num_df)
set.seed(0)
corrplot::corrplot.mixed(M, lower="ellipse", upper="number", order="hclust")

 

 

---

PCA

En este caso se ha aplicado análisis de componentes principales (PCA) para evaluar la reducción de dimensión del dataset de variables numéricas obtenido después del proceso de limpieza de datos. Siguiendo la lógica del análisis de componentes principales se busca determinar la cantidad de componentes que permitirían describir la varianza del conjunto de datos de variables numéricas en un 90%. Como podemos ver en la siguiente gráfica, podemos capturar un 90% de la varianza con los primeros 5 componentes principales.

num_df <- subset(df, select = -c(id, zona, tipo, barrio))
scaled_df <- scale(num_df)
res.pca <- prcomp(scaled_df)
scree_plot <- fviz_eig(res.pca, addlabels = TRUE)


cumulative_variance <- cumsum(res.pca$sdev^2) / sum(res.pca$sdev^2) * 100
cum_var_plot <- ggplot(data = data.frame(PC = 1:length(cumulative_variance), 
                                         CumulativeVariance = cumulative_variance),
                       aes(x = PC, y = CumulativeVariance)) +
  geom_line() +
  geom_point() +
  geom_hline(yintercept = 90, col = "red", linetype = "dashed") +
  labs(x = "Principal Components", y = "Cumulative Variance (%)", 
       title = "Cumulative Variance Explained by PCA") +
  theme_minimal()


grid.arrange(scree_plot, cum_var_plot, ncol = 2)

 

 

---

Contribución de variables a los componentes principales (PCs)

• Las variables latitud, longitud, habitaciones y price_per_sqm son las más influyentes sobre varios componentes principales, lo que las convierte en variables importantes para capturar la variación de los datos.

• Existen variables dominantes para el tercer y cuarto componente. En otras palabras los componetes PC3 y PC4 se encuentran casi totalmente compuestos por latitud y longitud respectivamente.

• Los componentes principales PC1, PC2 y PC5 muestran contribuciones más balanceadas entre la contribución de varias variables, por lo que pueden estar capturando relaciones más complejas de los datos.

p1 <- fviz_contrib(res.pca, choice = "var", axes = 1, top = 10, title="Variable contrib PC1")


p2 <- fviz_contrib(res.pca, choice = "var", axes = 2, top = 10, title="Variable contrib PC2")


p3 <- fviz_contrib(res.pca, choice = "var", axes = 3, top = 10, title="Variable contrib PC3")


p4 <- fviz_contrib(res.pca, choice = "var", axes = 4, top = 10, title="Variable contrib PC4")


p5 <- fviz_contrib(res.pca, choice = "var", axes = 5, top = 10, title="Variable contrib PC5")


p6 <- fviz_contrib(res.pca, choice = "var", axes = 1:5, top = 10, title="Overall PC1-PC5")


blankPlot <- ggplot() + geom_blank(aes(1, 1)) + theme_void()


grid.arrange(p1, p2, p3,
             p4, p5, p6,
             blankPlot, blankPlot, blankPlot,  # 3 placeholders
             ncol = 3)

 

 

---

Calidad de representación (\(Cos^{2}\))

• \(Cos^{2}\): Mide la calidad de la representación de una variable en el espacio definido por PC1-PC5. Es el cuadrado del coseno del ángulo entre el vector de la variable y el vector del componente principal. Lo mostrado por esta gráfica muestra cuanta varianza por variable puede explicarse por un componente particular.

  • Valores de \(Cos^{2}\) cercanos a 1: Las variables longitud, latitud, price_per_sqm y parqueaderos están bien representadas por los primeros 5 componentes. Del mismo modo que en la gráfica de contribuciones, si una de estas variables se muestra cercana a la dirección del eje del componente principal, el componente se encuentra capturará la mayor parte de la varianza asociada a esta variable.

  • Valores de \(Cos^{2}\) cercanos a 0: Aunque para los primeros 5 componentes no hay variables con \(Cos^{2}\) realmente cercanos a 0, se tienen valores intermedios (0.26-0.36) donde los componentes tomados no capturan la mayor parte de varianza asociada a la variable. Esto es especialmente importante para las variables de estrato (0.294) y preciom (0.286).

loadings <- res.pca$rotation[, 1:5]


cos2 <- rowSums(loadings^2)


plot <- plot_ly(x = ~loadings[, 1], y = ~loadings[, 2], z = ~loadings[, 5], 
                text = rownames(loadings), 
                color = ~cos2, # Color by cos2
                colors = c("#00AFBB", "#E7B800", "#FC4E07"), 
                type = "scatter3d", 
                mode = "markers", 
                marker = list(size = 5)) %>%
  layout(title = "3D PCA Variables Plot (Cos2)",
         scene = list(xaxis = list(title = 'PC1'),
                      yaxis = list(title = 'PC2'),
                      zaxis = list(title = 'PC5')))


plot

loadings <- res.pca$rotation[, 1:5]


contrib <- rowSums(loadings^2) / sum(res.pca$sdev^2) * 100


plot <- plot_ly(x = ~loadings[, 1], y = ~loadings[, 2], z = ~loadings[, 5], 
                text = rownames(loadings), 
                color = ~contrib, 
                colors = c("#FF7F00",  "#034D94"), 
                type = "scatter3d", 
                mode = "markers", 
                marker = list(size = 5)) %>%
  layout(title = "3D PCA Variables Plot (Contribution)",
         scene = list(xaxis = list(title = 'PC1'),
                      yaxis = list(title = 'PC2'),
                      zaxis = list(title = 'PC3')))


plot

 

 

---

Biplot de individuos y variables

p1 <- fviz_pca_biplot(res.pca,
                labelsize=3,
                addEllipses = T,
                repel=TRUE,
                # Individuals
                geom.ind = "point",
                geom.var = c("point", "text"),
                fill.ind = df$tipo,
                pointshape = 21 ,
                pointsize = 2,
                alpha.ind=0.1,
                # Variables
                col.var="black",
                alpha.var =1, 
                title = "PCs vs Tipo",
                legend.title = list(fill = "Tipo")+
  
  guides(color=guide_legend("Tipo"),fill= T))


p2 <- fviz_pca_biplot(res.pca,
                labelsize=3,
                addEllipses = T,
                repel=TRUE,
                # Individuals
                geom.ind = "point",
                geom.var = c("point", "text"),
                fill.ind = df$zona,
                pointshape = 21 ,
                pointsize = 2,
                alpha.ind=0.1,
                # Variables
                col.var="black",
                alpha.var =1, 
                title = "PCs vs Zona",
                legend.title = list(fill = "Zona")+
                guides(color=guide_legend("Zona"),fill= T))

p1 | p2

• Representación de los datos: Los puntos mostrados toman su color con respecto a si el tipo y zona de las categorías de cada variable.

• Componetes Principales (PCs): El eje x y el eje y corresponden con el primer y segundo componente principal respectivamente (Dim1 y Dim2). Estos dos componentes estarían explicando un 62% de la varianza en el conjunto de variables numéricas.

• Variables: Los nombres de las variables usadas representan su contribución a los dos primeros componentes principales (Dim1 y Dim2). La dirección y espacio del origen hasta el punto donde aparece el nombre de la variable indican la contribución de cada variable al componente principal. Si tomamos como referencia a Dim1 encontraremos que parqueadero, preciom, banios y areaconst se encuentran fuertemente alineados con este primer componente, lo que sugiere que estas variables contribuyen de forma significativa a la varianza capturada por el primer componente principal Dim1. Por otro lado, Dim2 se encuentra principalmente asociado a price_per_sqm o precio por metro cuadrado y en menor medida por habitaciones y estrato.

• Elipses: Las elipses representan los intervalos del 95% de los puntos de datos que pertenecen a cada categoría de tipo. Estas elipses proveen una indicación visual de cómo los grupos se sobrelapan mostrando cierta similaridad entre apartamentos y casas. En relación con la categoría zona vemos tres subgrupos:

  • Grupo 1: Zona Sur y Norte con sobrelapes importantes. Recordando las medias y medianas de price_per_sqm en estas dos zonas son muy similares para la categoría de apartamento.
  • Grupo 2: Zona Oriente y Zona Centro, en este caso las medias de price_per_sqm no parecen distinguir sobre la variable tipo. Estas presentan precios por metro cuadrado similares para apartamentos de zona oriente y casas de zona centro.
  • Grupo 3: Zona Oeste, presenta algunas similitudes con propiedades de tipo casa para la zona sur.

df %>%
  group_by(tipo, zona) %>%
  mutate(Q1 = quantile(price_per_sqm, 0.25),
         Q3 = quantile(price_per_sqm, 0.75),
         IQR = Q3 - Q1,
         lower_bound = Q1 - 1.5 * IQR,
         upper_bound = Q3 + 1.5 * IQR) %>%
  filter(price_per_sqm >= lower_bound, price_per_sqm <= upper_bound) -> df_filtered


results <- df_filtered %>%
  group_by(tipo, zona) %>%
  summarise(mean = mean(price_per_sqm),
            median = median(price_per_sqm),
            std_dev = sd(price_per_sqm), .groups = 'drop')

paged_table(results, options = list(rows.print = 10)) %>%
  kable(caption = 'Valores de price_per_sqm agrupados por Tipo y Zona')

Valores de price_per_sqm agrupados por Tipo y Zona

tipo	zona	mean	median	std_dev
apartamento	Zona Centro	2.224967	2.192982	0.6725103
apartamento	Zona Norte	2.968890	2.919872	0.8786715
apartamento	Zona Oeste	3.878913	3.882353	1.0140589
apartamento	Zona Oriente	1.519318	1.586431	0.4760633
apartamento	Zona Sur	2.948351	2.926829	0.6402131
casa	Zona Centro	1.610017	1.500000	0.4713779
casa	Zona Norte	1.740763	1.623377	0.5455837
casa	Zona Oeste	2.239288	1.995223	0.9553809
casa	Zona Oriente	1.354955	1.287059	0.5018179
casa	Zona Sur	2.254948	2.133333	0.8138701

 

 

---

Extra: Distribución de Zona y Tipo en componentes PC1, PC2 y PC5

pc1 <- res.pca$x[, 1]
pc2 <- res.pca$x[, 2]
pc5 <- res.pca$x[, 5]


plot_ly(x = ~pc1, y = ~pc2, z = ~pc5, 
        color = ~df$tipo, 
        colors = "Set2", 
        marker = list(size = 3),
        type = "scatter3d", 
        mode = "markers") %>%
  layout(title = "3D PCA Plot: PCs 1,2,5 vs Tipo",
         scene = list(xaxis = list(title = 'PC1'),
                      yaxis = list(title = 'PC2'),
                      zaxis = list(title = 'PC5')))

plot_ly(x = ~pc1, y = ~pc2, z = ~pc5, 
        color = ~df$zona, 
        colors = "Set2",
        marker = list(size = 3),
        type = "scatter3d", 
        mode = "markers") %>%
  layout(title = "3D PCA Plot: PCs 1,2,5 vs Zona",
         scene = list(xaxis = list(title = 'PC1'),
                      yaxis = list(title = 'PC2'),
                      zaxis = list(title = 'PC5')))

Análisis de Correspondencia

 

 

---

CA: Tipo vs Zona

A partir del análisis de correspondencia se determina que las variables tipo y zona se encuentran relacionadas. Como se puede ver en la siguiente gráfica, las zonas norte, oeste y sur están asociadas principalmente con propiedades de tipo apartamento. Para la zona sur no hay una división entre las frecuencias por tipo de propiedad.

cat_df <- subset(df, select = c(zona, tipo, barrio))
contingency_table <- table(cat_df$tipo, cat_df$zona)
chisq <- chisq.test(contingency_table)
numb <- formatC(chisq$p.value, format = "e", digits = 2)
cat("Las variables `tipo` y `zona` se encuentran significativamente asociadas de acuerdo con ($\\chi^{2}$=", chisq$statistic, " p-value = ", numb,").\n")

Las variables tipo y zona se encuentran significativamente asociadas de acuerdo con (\(\chi^{2}\)= 562.5876 p-value = 1.93e-120 ).

# Convert the data as a table
dt <- as.table(as.matrix(contingency_table))

balloonplot(t(dt), main ="zona vs tipo", xlab ="", ylab="",
            label = FALSE, show.margins = FALSE)

res.ca <- CA(contingency_table, graph = TRUE)
eig.val <- get_eigenvalue(res.ca)
cat("La asociación de variables `tipo` y `zona` pueden representarse con una sóla dimensión.\n")

La asociación de variables tipo y zona pueden representarse con una sóla dimensión.

eign_df <- as.data.frame(eig.val)
paged_table(eign_df, options = list(rows.print = 10)) %>%
  kable(caption = 'CA para variables de Tipo y Zona')

CA para variables de Tipo y Zona

	eigenvalue	variance.percent	cumulative.variance.percent
Dim.1	0.0887642	100	100

 

 

---

CA: Barrio vs Tipo

# categories
df <- result$cleaned_data
cat_df <- subset(df, select = c(zona, tipo, barrio))
contingency_table <- table(cat_df$barrio, cat_df$tipo)
chisq <- chisq.test(contingency_table, simulate.p.value = TRUE)
numb <- formatC(chisq$p.value, format = "e", digits = 2)
cat("Las variables `barrio` y `zona` se encuentran significativamente asociadas de acuerdo con  ($\\chi^{2}$=", chisq$statistic," p-value = ", numb,").\n")

Las variables barrio y zona se encuentran significativamente asociadas de acuerdo con (\(\chi^{2}\)= 1818.37 p-value = 5.00e-04 ).

res.ca <- CA(contingency_table, graph = TRUE)
eig.val <- get_eigenvalue(res.ca)
cat("La asociación de variables `barrio` y `zona` pueden representarse con una sóla dimensión.\n")

La asociación de variables barrio y zona pueden representarse con una sóla dimensión.

eign_df <- as.data.frame(eig.val)
paged_table(eign_df, options = list(rows.print = 10)) %>%
  kable(caption = 'CA para variables de Tipo y Zona')

CA para variables de Tipo y Zona

	eigenvalue	variance.percent	cumulative.variance.percent
Dim.1	0.2868997	100	100

 

 

---

MCA

Del mismo modo que en la sección de PCA, en este caso se ha aplicado análisis de componentes múltiples (MCA) a las variables categóricas para evaluar la reducción de dimensión del dataset obtenido después del proceso de limpieza de datos. Siguiendo la lógica del análisis de componentes múltiples se busca determinar la cantidad de componentes que permitirían describir la varianza del conjunto de datos de variables numéricas en un 90%. Como podemos ver en la siguiente gráfica, podemos capturar un 90% de la varianza con los primeros 285 componentes principales. A diferencia de PCA, en MCA el primer componente está aportando 0.7% por lo que se requiere una gran cantidad de componentes para capturar el 90% de varianza.

res.mca <- MCA(cat_df, graph = FALSE)

scree_plot <- fviz_screeplot(res.mca, addlabels = TRUE, 
                             title = "Scree Plot of MCA")


cumulative_variance <- cumsum(res.mca$eig[, 2])


cum_var_plot <- ggplot(data = data.frame(Dim = 1:length(cumulative_variance), 
                                         CumulativeVariance = cumulative_variance),
                       aes(x = Dim, y = CumulativeVariance)) +
  geom_line() +
  geom_point() +
  geom_hline(yintercept = 90, col = "red", linetype = "dashed") +
  geom_vline(xintercept = 285, col = "red", linetype = "dashed") +
  labs(x = "Dimensions", y = "Cumulative Variance (%)", 
       title = "Cumulative Variance") +
  theme_minimal()

grid.arrange(scree_plot, cum_var_plot, ncol = 2)

cumulative_variance <- cumsum(res.mca$eig[, 2])


num_components <- which(cumulative_variance >= 90)[1]

 

 

---

Contribución de variables a los componentes principales (PCs)

• Como se aprecia en la anterior gráfica y en la gráfica MCA, las primeras dos dimensiones sólo explican un 1.3% de la varianza del conjunto de datos. Sin embargo, podemos ver que al ubicar los individuos u observaciones del dataset en el espacio de componentes MCA se crean subconjuntos de datos con características similares. Por otro lado, los grupos de individuos cercanos a cierta variable (▲) comparten características definidas por ella misma.

• La posición de una variable en el gráfica indica la asociación de dicha variable con la dimensión correspondiente. Del mismo modo que en PCA, las variables cercanas al origen tienen menos influencia en las dimensiones Dim1 y Dim2. Por el contrario, las categorías alejadas al origen como Zona Norte, Zona Oeste, Zona Oriente y Zona Sur tiene una alta influencia o contribución.

• Respecto a los valores de calidad de representación (\(Cos^{2}\)), las zonas Zona Norte, Zona Oeste, Zona Oriente, Zona Sur y en menor medida, los tipos de propiedad casa y apartamento parecen capturar la mayor variabilidad. Esto reafirma lo encontrado en la sección de PCA, donde el tipo y zona asociada con el inmueble muestran valores promedio distintos.

p0 <- fviz_mca_biplot(res.mca,
                labelsize=3,
                addEllipses = F,
                repel=TRUE,
                # Individuals
                geom.ind = "point",
                geom.var = c("point"),
                pointsize = 2,
                alpha.ind=0.07,
                alpha.var=0.8,
                # Variables
                col.var="black",
                #alpha.var =1, 
                title = "MCA")


grp1 <- as.factor(cat_df$tipo)
p1 <- fviz_mca_ind(res.mca, label="none", habillage=grp1,
                   title="PCs vs Tipo",
       addEllipses=TRUE, ellipse.level=0.95)
      
grp2 <- as.factor(cat_df$zona)
p2 <- fviz_mca_ind(res.mca, label="none", habillage=grp2,
                   title="PCs vs Zona",
       addEllipses=TRUE, ellipse.level=0.95)


total_contrib <- rowSums(res.mca$var$contrib)

top_contributors <- names(sort(total_contrib, decreasing = TRUE)[1:15])

p3 <- fviz_mca_var(res.mca, 
             col.var = "contrib", 
             repel = T,
             title = "Top 15 Vars per Contrib",
             gradient.cols =  c("#FF7F00",  "#034D94"),
             select.var = list(name = top_contributors))



total_contrib <- rowSums(res.mca$var$cos2)

top_contributors <- names(sort(total_contrib, decreasing = TRUE)[1:15])
suppressMessages(suppressWarnings(p4 <- fviz_mca(res.mca,
               # Variables
               select.var = list(name = top_contributors),
               col.var = "cos2",
               gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
               # Individuals
               geom.ind = "point",
               pointsize = 1,
               alpha.ind = 0.07,
               repel = TRUE,
               title = "Top 15 Vars per Cos2")))


p0

p4

p3

 

 

---

Biplot de individuos y variables

Los puntos mostrados toman su color con respecto a si el tipo y zona de las categorías de cada variable. A diferencia de PCA, MCA permite distinguir más efectivamente entre tipo casa-apartamento y las diferentes zonas: Zona Norte, Zona Oeste, Zona Oriente y Zona Sur. Al graficar los intervalos del 95% de los puntos de los datos que pertenecen a cada categoría de tipo se encuentra una división difusa para los distintos inmuebles. Por el contrario, las zonas están bien definidas para los individuos más cercanos entre sí que corresponderían con una u otra zona. Los individuos que no hacen parte de los intervalos podrían presentar características mixtas.

p1 | p2 

Análisis de Conglomerados 1

 

 

---

num_df y cat_df

Definimos a nuestros conjuntos iniciales como num_df y cat_df.

• num_df posee las variables: estrato, preciom, areaconst, parqueaderos, banios, habitaciones, longitud, latitud y price_per_sqm.
• cat_df posee las variables: zona, tipo y barrio

 

 

---

Extracción de componentes PCA

result <- clean_vivienda_data(vivienda)
df <- result$cleaned_data
cat_df <- subset(df, select = c(zona, tipo, barrio))
num_df <- subset(df, select = -c(id, zona, tipo, barrio))
scaled_df <- scale(num_df)
res.pca <- prcomp(scaled_df)


select_pca_components <- function(res.pca, variance_threshold = 90) {
  
  cumulative_variance <- cumsum(res.pca$sdev^2) / sum(res.pca$sdev^2) * 100
  
  
  num_components <- which(cumulative_variance >= variance_threshold)[1]
  
  
  selected_components <- res.pca$x[, 1:num_components]
  
  
  return(list(num_components = num_components, 
              cumulative_variance = cumulative_variance[1:num_components],
              selected_components = selected_components))
}


result <- select_pca_components(res.pca)
cat("Número de componentes necesarios (PCA) para 90% de varianza:", result$num_components, "\n")

Número de componentes necesarios (PCA) para 90% de varianza: 6

cat("Varianza acumulada:", round(result$cumulative_variance))

Varianza acumulada: 40 62 74 82 88 93

selected_pca_components <- result$selected_components

 

 

---

Extracción de componentes MCA

res.mca <- MCA(cat_df, ncp=300, graph = FALSE)

select_mca_components <- function(res.mca, variance_threshold = 90) {

  cumulative_variance <- cumsum(res.mca$eig[, 2]) 
  num_components <- which(cumulative_variance >= variance_threshold)[1]
  selected_var_components <- res.mca$var$coord[, 1:num_components]
  selected_ind_components <- res.mca$ind$coord[, 1:num_components]
  return(list(num_components = num_components, 
              cumulative_variance = cumulative_variance[1:num_components],
              selected_var_components = selected_var_components,
              selected_ind_components = selected_ind_components))
}


result_mca <- select_mca_components(res.mca)
cat("Número de componentes necesarios (MCA) para 90% de varianza:", result_mca$num_components)

Número de componentes necesarios (MCA) para 90% de varianza: 285

selected_mca_components <- result_mca$selected_ind_components

combined_components <- cbind(selected_pca_components, selected_mca_components)
combined_components <- as.data.frame(combined_components)

 

 

---

KMeans clustering

set.seed(2024) 
cluster_range <- 2:10
avg_silhouette_scores <- numeric(length(cluster_range))

for (i in cluster_range) {
  
  kmeans_result <- kmeans(combined_components, centers = i)
  silhouette_scores <- silhouette(kmeans_result$cluster, dist(combined_components))
  avg_silhouette_scores[i - 1] <- mean(silhouette_scores[, 3])
}


optimal_silhouette_score <- max(avg_silhouette_scores)
optimal_clusters <- cluster_range[which.max(avg_silhouette_scores)]
cat("Número óptimo de clusters:", optimal_clusters,"\n")

Número óptimo de clusters: 4

cat("Silhouette score obtenido:", round(optimal_silhouette_score, 4))

Silhouette score obtenido: 0.0637

plot(cluster_range, avg_silhouette_scores, type = "b", 
     xlab = "Nro de clusters", ylab = "Silhouette Score promedio",
     main = "Silhouette Score vs Nro Clusters")

final_kmeans <- kmeans(combined_components, centers = optimal_clusters)
combined_components$cluster <- as.factor(final_kmeans$cluster)

 

 

---

Descripción de clusters obtenidos

• Cluster 1 muestra valores más altos de price_per_sqm, preciom y areconst, en simultáneo tiene valores bajos para parqueaderos y habitaciones. Este cluster está constituído principalmente por propiedades de tipo apartamento en la zona sur.
• Cluster 2 muestra valores más moderados y poco dispersos para la mayoría de variables presentando únicamente más variación en estrato. Este cluster está constituído principalmente por propiedades de tipo apartamento en la zona norte.
• Cluster 3 muestra una media menor que cluster 1 en la variable price_per_sqm por lo que su precio por metro cuadrado sería atractivo en ventas. Del mismo modo muestra variación en valores de parqueaderos, banios y habitaciones. Llama la atención que la mayoría de sus valores se encuentran asociados al estrato 5. Este cluster está constituído principalmente por propiedades de tipo casa en la zona sur.
• Cluster 4 muestra valores bajos para la mayoría de variables comparables al cluster 2, tiene más variación en las variables banios y especialmente en habitaciones, donde supera las medias del resto de clusters. Este cluster está constituído principalmente por propiedades de tipo casa en la zona oriente.

En resumen, el cluster 1 representa propiedades más caras y con mayor variación en sus atributos mientras que el cluster 4 muestra las propiedades de menor precio por metro cuadrado. Los cluster 2 y 3 presentan características intermedias, disintiguiendo que las propiedades del cluster 3 tienden a pertenecer a estratos más altos y a tener áreas mayores

data_with_clusters <- cbind(df, cluster = combined_components$cluster)

suppressMessages(suppressWarnings(
p1 <- ggplot(data_with_clusters, aes_string(x = "cluster", y = "price_per_sqm", fill = "cluster")) +
    geom_boxplot() +
    labs(title = paste("price_per_sqm"),
         x = "Cluster",
         y = "price_per_sqm") +
    theme_minimal() +
    scale_fill_brewer(palette = "Set3")))
suppressMessages(suppressWarnings(
p2 <- ggplot(data_with_clusters, aes_string(x = "cluster", y = "preciom", fill = "cluster")) +
    geom_boxplot() +
    labs(title = paste("preciom"),
         x = "Cluster",
         y = "preciom") +
    theme_minimal() +
    scale_fill_brewer(palette = "Set3")))

suppressMessages(suppressWarnings(
p3 <- ggplot(data_with_clusters, aes_string(x = "cluster", y = "areaconst", fill = "cluster")) +
    geom_boxplot() +
    labs(title = paste("areaconst"),
         x = "Cluster",
         y = "areaconst") +
    theme_minimal() +
    scale_fill_brewer(palette = "Set3")))

suppressMessages(suppressWarnings(
p4 <- ggplot(data_with_clusters, aes_string(x = "cluster", y = "parqueaderos", fill = "cluster")) +
    geom_boxplot() +
    labs(title = paste("parqueaderos"),
         x = "Cluster",
         y = "parqueaderos") +
    theme_minimal() +
    scale_fill_brewer(palette = "Set3")))

suppressMessages(suppressWarnings(
p5 <- ggplot(data_with_clusters, aes_string(x = "cluster", y = "banios", fill = "cluster")) +
    geom_boxplot() +
    labs(title = paste("banios"),
         x = "Cluster",
         y = "banios") +
    theme_minimal() +
    scale_fill_brewer(palette = "Set3")))

suppressMessages(suppressWarnings(
p6 <- ggplot(data_with_clusters, aes_string(x = "cluster", y = "habitaciones", fill = "cluster")) +
    geom_boxplot() +
    labs(title = paste("habitaciones"),
         x = "Cluster",
         y = "habitaciones") +
    theme_minimal() +
    scale_fill_brewer(palette = "Set3")))

suppressMessages(suppressWarnings(
p7 <- ggplot(data_with_clusters, aes_string(x = "cluster", y = "estrato", fill = "cluster")) +
    geom_boxplot() +
    labs(title = paste("estrato"),
         x = "Cluster",
         y = "estrato") +
    theme_minimal() +
    scale_fill_brewer(palette = "Set3")))


blankPlot <- ggplot() + geom_blank(aes(1, 1)) + theme_void()

grid.arrange(p1, p2, p3, p4, 
             p5, p6, p7,
             blankPlot, blankPlot, blankPlot, blankPlot, # 4 placeholders
             ncol = 4)

p1 <- ggplot(data_with_clusters, aes(x = as.factor(cluster), fill = tipo)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(title = "tipo",
       x = "Cluster",
       y = "Percentage",
       fill = "Tipo") +
  theme_minimal() +
  scale_fill_brewer(palette = "RdYlBu")


p2 <- ggplot(data_with_clusters, aes(x = as.factor(cluster), fill = zona)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(title = "zona",
       x = "Cluster",
       y = "Percentage",
       fill = "Zona") +
  theme_minimal() +
  scale_fill_brewer(palette = "RdYlBu")

p1 | p2

Distribución geográfica

df_sf <- st_as_sf(data_with_clusters, coords = c("longitud", "latitud"), crs = 4326)
df$price_per_sqm[is.na(data_with_clusters$price_per_sqm) | is.infinite(data_with_clusters$price_per_sqm)] <- NA
map <- mapview(df_sf, 
               zcol = c("tipo", "estrato", "price_per_sqm", "cluster"), 
               alpha=0, 
               legend = TRUE, 
               cex = 5,
               addLayersControl = TRUE, 
               hide = TRUE,
               popup = popupTable(df_sf, 
                                  zcol = c("estrato", "tipo", "price_per_sqm", "cluster")))

map

 

 

Análisis de Conglomerados 2

 

 

---

num_df, cat_df y estrato como categoría

Definimos a nuestros conjuntos iniciales como num_df y cat_df.

• num_df posee las variables: preciom, areaconst, parqueaderos, banios, habitaciones, longitud, latitud y price_per_sqm.
• cat_df posee las variables: estrato, zona y tipo

result <- clean_vivienda_data(vivienda)
df <- result$cleaned_data
cat_df <- subset(df, select = c(zona, tipo, estrato))
cat_df$estrato <- as.character(cat_df$estrato)
num_df <- subset(df, select = -c(id, zona, tipo, barrio, estrato))
scaled_df <- scale(num_df)
res.pca <- prcomp(scaled_df)


result_pca <- select_pca_components(res.pca)
cat("Número de componentes necesarios (PCA) para 90% de varianza:", result_pca$num_components, "\n")

Número de componentes necesarios (PCA) para 90% de varianza: 5

cat("Varianza acumulada:", round(result_pca$cumulative_variance), "\n")

Varianza acumulada: 42 62 76 85 91

selected_pca_components <- result_pca$selected_components

res.mca <- MCA(cat_df, ncp = 300, graph = FALSE)
result_mca <- select_mca_components(res.mca)
cat("Número de componentes necesarios (MCA) para 90% de varianza:", result_mca$num_components)

Número de componentes necesarios (MCA) para 90% de varianza: 7

selected_mca_components <- result_mca$selected_ind_components
combined_components <- cbind(selected_pca_components, selected_mca_components)
combined_components <- as.data.frame(combined_components)

 

 

---

KMeans clustering

set.seed(2024)
cluster_range <- 2:10
avg_silhouette_scores <- numeric(length(cluster_range))

for (i in cluster_range) {

  kmeans_result <- kmeans(combined_components, centers = i)
  silhouette_scores <- silhouette(kmeans_result$cluster, dist(combined_components))
  avg_silhouette_scores[i - 1] <- mean(silhouette_scores[, 3])
}


optimal_silhouette_score <- max(avg_silhouette_scores)
optimal_clusters <- cluster_range[which.max(avg_silhouette_scores)]
cat("Número óptimo de clusters:", optimal_clusters, "\n")

Número óptimo de clusters: 2

cat("Silhouette score obtenido:", round(optimal_silhouette_score, 4))

Silhouette score obtenido: 0.2713

plot(cluster_range, avg_silhouette_scores, type = "b", 
     xlab = "Nro de clusters", ylab = "Silhouette Score promedio",
     main = "Silhouette Score vs Nro Clusters")

final_kmeans <- kmeans(combined_components, centers = optimal_clusters)
combined_components$cluster <- as.factor(final_kmeans$cluster)

 

 

---

Descripción de clusters obtenidos

• Cluster 1 muestra valores más altos de preciom y areconst, sin embargo los valores de price_per_sqm o precio por metro cuadrado muestran mayor variación que el cluster 2. En general este cluster tiene valores más altos para todas las variables numéricas. Respecto a las variables categóricas, el clúster 1 sólo incluye a las zonas Sur, Oeste y una proporción pequeña de la zona Norte, zonas que aparecen asociadas en alta proporción a los estratos 6, 5 y 4.
• Cluster 2 muestra valores más bajos y poco dispersos para la mayoría de variables. Este segmento tiene principalemente propiedades de la Zona Norte y Zona Centro asociadas también con estratos 5, 4 y 3.

result <- clean_vivienda_data(vivienda)
df <- result$cleaned_data
data_with_clusters <- cbind(df, cluster = combined_components$cluster)

p1 <- ggplot(data_with_clusters, aes_string(x = "cluster", y = "price_per_sqm", fill = "cluster")) +
    geom_boxplot() +
    labs(title = paste("price_per_sqm"),
         x = "Cluster",
         y = "price_per_sqm") +
    theme_minimal() +
    scale_fill_brewer(palette = "Set3")

p2 <- ggplot(data_with_clusters, aes_string(x = "cluster", y = "preciom", fill = "cluster")) +
    geom_boxplot() +
    labs(title = paste("preciom"),
         x = "Cluster",
         y = "preciom") +
    theme_minimal() +
    scale_fill_brewer(palette = "Set3")


p3 <- ggplot(data_with_clusters, aes_string(x = "cluster", y = "areaconst", fill = "cluster")) +
    geom_boxplot() +
    labs(title = paste("areaconst"),
         x = "Cluster",
         y = "areaconst") +
    theme_minimal() +
    scale_fill_brewer(palette = "Set3")


p4 <- ggplot(data_with_clusters, aes_string(x = "cluster", y = "parqueaderos", fill = "cluster")) +
    geom_boxplot() +
    labs(title = paste("parqueaderos"),
         x = "Cluster",
         y = "parqueaderos") +
    theme_minimal() +
    scale_fill_brewer(palette = "Set3")


p5 <- ggplot(data_with_clusters, aes_string(x = "cluster", y = "banios", fill = "cluster")) +
    geom_boxplot() +
    labs(title = paste("banios"),
         x = "Cluster",
         y = "banios") +
    theme_minimal() +
    scale_fill_brewer(palette = "Set3")

data_with_clusters$estrato <- as.character(data_with_clusters$estrato)
p6 <- ggplot(data_with_clusters, aes_string(x = "cluster", y = "habitaciones", fill = "cluster")) +
    geom_boxplot() +
    labs(title = paste("habitaciones"),
         x = "Cluster",
         y = "habitaciones") +
    theme_minimal() +
    scale_fill_brewer(palette = "Set3")

  
blankPlot <- ggplot() + geom_blank(aes(1, 1)) + theme_void()

grid.arrange(p1, p2, p3, 
             p4, p5, p6, 
             blankPlot, blankPlot, blankPlot, 
             ncol = 3)

p1 <- ggplot(data_with_clusters, aes(x = as.factor(cluster), fill = tipo)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(title = "tipo",
       x = "Cluster",
       y = "Percentage",
       fill = "Tipo") +
  theme_minimal() +
  scale_fill_brewer(palette = "RdYlBu")

p2 <- ggplot(data_with_clusters, aes(x = as.factor(cluster), fill = zona)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(title = "zona",
       x = "Cluster",
       y = "Percentage",
       fill = "Zona") +
  theme_minimal() +
  scale_fill_brewer(palette = "RdYlBu")

p3 <- ggplot(data_with_clusters, aes(x = as.factor(cluster), fill = estrato)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(title = "estrato",
       x = "Cluster",
       y = "Percentage",
       fill = "Estrato") +
  theme_minimal() +
  scale_fill_brewer(palette = "RdYlBu")


p1 | p2 | p3

 

 

---

Distribución geográfica

df_sf <- st_as_sf(data_with_clusters, coords = c("longitud", "latitud"), crs = 4326)
df$price_per_sqm[is.na(data_with_clusters$price_per_sqm) | is.infinite(data_with_clusters$price_per_sqm)] <- NA
map <- mapview(df_sf, 
               zcol = c("tipo", "estrato", "price_per_sqm", "cluster"), 
               alpha=0, 
               legend = TRUE, 
               cex = 5,
               addLayersControl = TRUE, 
               hide = TRUE,
               popup = popupTable(df_sf, 
                                  zcol = c("estrato", "tipo", "price_per_sqm", "cluster")))

map

 

 

---