Instalación de paquetes
Estandarizamos nombres, tipamos variables y quitamos duplicados por
id.
data("vivienda") # proviene de paqueteMODELOS
raw <- vivienda %>% as_tibble() %>% clean_names()
df <- raw %>%
rename(
precio_mill = preciom,
area_const = areaconst
) %>%
mutate(
piso = readr::parse_number(piso),
estrato = factor(estrato),
tipo = factor(tipo),
zona = factor(zona),
barrio = factor(barrio)
) %>%
distinct(id, .keep_all = TRUE)
glimpse(df)
## Rows: 8,320
## Columns: 13
## $ id <dbl> 1147, 1169, 1350, 5992, 1212, 1724, 2326, 4386, 1209, 159…
## $ zona <fct> Zona Oriente, Zona Oriente, Zona Oriente, Zona Sur, Zona …
## $ piso <dbl> NA, NA, NA, 2, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, …
## $ estrato <fct> 3, 3, 3, 4, 5, 5, 4, 5, 5, 5, 6, 4, 5, 6, 4, 5, 5, 4, 5, …
## $ precio_mill <dbl> 250, 320, 350, 400, 260, 240, 220, 310, 320, 780, 750, 62…
## $ area_const <dbl> 70, 120, 220, 280, 90, 87, 52, 137, 150, 380, 445, 355, 2…
## $ parqueaderos <dbl> 1, 1, 2, 3, 1, 1, 2, 2, 2, 2, NA, 3, 2, 2, 1, 4, 2, 2, 2,…
## $ banios <dbl> 3, 2, 2, 5, 2, 3, 2, 3, 4, 3, 7, 5, 6, 2, 4, 4, 4, 3, 2, …
## $ habitaciones <dbl> 6, 3, 4, 3, 3, 3, 3, 4, 6, 3, 6, 5, 6, 2, 5, 5, 4, 3, 3, …
## $ tipo <fct> Casa, Casa, Casa, Casa, Apartamento, Apartamento, Apartam…
## $ barrio <fct> 20 de julio, 20 de julio, 20 de julio, 3 de julio, acopi,…
## $ longitud <dbl> -76.51168, -76.51237, -76.51537, -76.54000, -76.51350, -7…
## $ latitud <dbl> 3.43382, 3.43369, 3.43566, 3.43500, 3.45891, 3.36971, 3.4…
colSums(is.na(df))
## id zona piso estrato precio_mill area_const
## 1 1 2636 1 1 1
## parqueaderos banios habitaciones tipo barrio longitud
## 1603 1 1 1 1 1
## latitud
## 1
recipes)Imputamos numéricos con KNN y normalizamos para PCA/cluster.
num_vars <- c("piso","precio_mill","area_const","parqueaderos","banios","habitaciones")
df <- df %>%
filter(!if_all(all_of(num_vars), ~ is.na(.x)))
rec <- recipe(~ ., data = df) %>%
update_role(id, new_role = "id") %>%
step_impute_knn(all_of(num_vars), neighbors = 7) %>% # K=7
step_normalize(all_of(num_vars))
prep_rec <- prep(rec, training = df)
df_proc <- bake(prep_rec, new_data = NULL)
# Guardamos dataset procesado (num vars ya estandarizadas)
df_proc %>% select(all_of(num_vars)) %>% summary()
## piso precio_mill area_const parqueaderos
## Min. :-1.1307 Min. :-1.1437 Min. :-1.0138 Min. :-0.7159
## 1st Qu.:-0.7078 1st Qu.:-0.6508 1st Qu.:-0.6640 1st Qu.:-0.7159
## Median :-0.2849 Median :-0.3161 Median :-0.3633 Median :-0.3098
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.5609 3rd Qu.: 0.3228 3rd Qu.: 0.3782 3rd Qu.: 0.2316
## Max. : 3.5213 Max. : 4.7620 Max. :10.9822 Max. : 7.8114
## banios habitaciones
## Min. :-2.17847 Min. :-2.4702
## 1st Qu.:-0.77811 1st Qu.:-0.4148
## Median :-0.07794 Median :-0.4148
## Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.62224 3rd Qu.: 0.2704
## Max. : 4.82330 Max. : 4.3813
Reducimos dimensionalidad y exploramos estructura.
x_num <- df_proc %>% select(all_of(num_vars))
# Filtrar filas no finitas por seguridad
x_num <- x_num %>% filter(if_all(everything(), ~ is.finite(.)))
# PCA sobre datos ya normalizados
pca_fit <- prcomp(as.matrix(x_num), center = FALSE, scale. = FALSE)
summary(pca_fit)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 1.8107 1.0552 0.8628 0.59964 0.55989 0.4361
## Proportion of Variance 0.5465 0.1856 0.1241 0.05993 0.05225 0.0317
## Cumulative Proportion 0.5465 0.7321 0.8561 0.91605 0.96830 1.0000
fviz_eig(pca_fit)
fviz_pca_biplot(pca_fit, geom.ind = "point", repel = TRUE)
Redefinimos localmente
cat_varspara evitar errores por ejecución fuera de orden, y nos aseguramos de no pasar NAs adaisy().
# Asegurar un id de fila estable en df_proc
df_proc <- df_proc %>% mutate(.row_id = dplyr::row_number())
# Definir variables disponibles desde df_proc
num_vars <- intersect(c("piso","precio_mill","area_const","parqueaderos","banios","habitaciones"),
names(df_proc))
cat_vars <- intersect(c("tipo","zona","estrato"), names(df_proc))
# Conjunto para distancia de Gower (sin NAs y con .row_id para alineación)
mix_df <- df_proc %>%
dplyr::select(dplyr::all_of(c(num_vars, cat_vars, ".row_id"))) %>%
dplyr::mutate(dplyr::across(dplyr::all_of(cat_vars), as.factor)) %>%
tidyr::drop_na()
# Guardar ids de filas que entran a clustering y quitar .row_id para daisy()
mix_ids <- mix_df$.row_id
mix_df_daisy <- dplyr::select(mix_df, - .row_id)
# Chequeo final de NAs
stopifnot(all(colSums(is.na(mix_df_daisy)) == 0))
# Distancia Gower y PAM con k por silueta
gower_dist <- cluster::daisy(mix_df_daisy, metric = "gower")
sil_scores <- sapply(2:6, function(k){
cluster::pam(gower_dist, k = k, diss = TRUE)$silinfo$avg.width
})
k_opt <- which.max(sil_scores) + 1; k_opt
## [1] 2
pam_fit <- cluster::pam(gower_dist, k = k_opt, diss = TRUE)
cluster_lbl <- factor(pam_fit$clustering)
# Ensamble de clusters con datos originales para reporte/mapa
df_cluster <- df %>%
dplyr::mutate(.row_id = dplyr::row_number()) %>%
dplyr::semi_join(tibble::tibble(.row_id = mix_ids), by = ".row_id") %>%
dplyr::mutate(cluster = cluster_lbl)
table(df_cluster$cluster, df_cluster$estrato)
##
## 3 4 5 6
## 1 814 725 984 696
## 2 639 1404 1766 1291
# Proyección a PCA usando exactamente las filas de mix_df
mix_num <- df_proc %>%
dplyr::filter(.row_id %in% mix_ids) %>%
dplyr::select(dplyr::all_of(num_vars)) %>%
as.matrix()
pcs_mix <- predict(pca_fit, newdata = mix_num)
pca_df <- tibble::as_tibble(pcs_mix[,1:2]) %>%
rlang::set_names(c("PC1","PC2")) %>%
dplyr::mutate(cluster = cluster_lbl)
ggplot2::ggplot(pca_df, ggplot2::aes(PC1, PC2, color = cluster)) +
ggplot2::geom_point(alpha = .6) +
ggplot2::labs(title = "Clusters proyectados en PC1–PC2")
pal <- colorFactor("Set2", df_cluster$cluster)
leaflet(df_cluster) %>%
addTiles() %>%
addCircleMarkers(~longitud, ~latitud,
color = ~pal(cluster), radius = 3,
stroke = FALSE, fillOpacity = .8) %>%
addLegend("bottomright", pal = pal, values = ~cluster, title = "Cluster")
cat_df <- df %>% select(tipo, zona, estrato) %>% drop_na()
mca_fit <- MCA(cat_df, graph = FALSE)
fviz_mca_ind(mca_fit, label = "none", habillage = df_cluster$cluster, addEllipses = TRUE)
factoextra::fviz_screeplot(mca_fit)
df_cluster %>%
group_by(cluster) %>%
summarise(
n = n(),
precio_med = median(precio_mill, na.rm = TRUE),
area_med = median(area_const, na.rm = TRUE),
banios_med = median(banios, na.rm = TRUE),
hab_med = median(habitaciones,na.rm = TRUE)
)
## # A tibble: 2 × 6
## cluster n precio_med area_med banios_med hab_med
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 3219 430 240 4 4
## 2 2 5100 279 90 2 3