Actividad 1 Modelación

Analisis de Componentes Principales PCA

Se procede a observar y analizar la base de datos siguiente. Este análisis tiene como objetivo comprender la estructura, distribución y relaciones entre las variables, así como identificar posibles patrones o anomalías que puedan ser relevantes para el estudio.

La base de datos contiene información relacionada con ciertas caracteristicas de viviendas de la ciudad de cali.

data("vivienda")
head(vivienda)

## # A tibble: 6 × 13
##      id zona    piso  estrato preciom areaconst parqueaderos banios habitaciones
##   <dbl> <chr>   <chr>   <dbl>   <dbl>     <dbl>        <dbl>  <dbl>        <dbl>
## 1  1147 Zona O… <NA>        3     250        70            1      3            6
## 2  1169 Zona O… <NA>        3     320       120            1      2            3
## 3  1350 Zona O… <NA>        3     350       220            2      2            4
## 4  5992 Zona S… 02          4     400       280            3      5            3
## 5  1212 Zona N… 01          5     260        90            1      2            3
## 6  1724 Zona N… 01          5     240        87            1      3            3
## # ℹ 4 more variables: tipo <chr>, barrio <chr>, longitud <dbl>, latitud <dbl>

Estas son las variables numéricas con datos faltantes y se muestran a continuación:

vivienda_new = vivienda[,c(4:9)]
colSums(is.na(vivienda_new))

##      estrato      preciom    areaconst parqueaderos       banios habitaciones 
##            3            2            3         1605            3            3

sapply(vivienda_new, class)

##      estrato      preciom    areaconst parqueaderos       banios habitaciones 
##    "numeric"    "numeric"    "numeric"    "numeric"    "numeric"    "numeric"

dim(vivienda_new)

## [1] 8322    6

Tenemos 8322 filas con 6 columnas (variables) lo cual nos da una idea con que prueba se puede imputar, en tanto que Shapiro tiene un limite no mayor de 5000 registros.

Usamos Kolmogorov-Smirnov para verificar la normalidad de las variables numérica:

Kolmogorov = sapply(vivienda_new, function(x) {
  n = sum(!is.na(x))
  if (n >= 3 && n <= 5000) {
     ks.test(x[!is.na(x)], "pnorm", mean = mean(x, na.rm = TRUE), sd = sd(x, na.rm = TRUE))$p.value
  } else if (n > 5000) {
    ks.test(x[!is.na(x)], "pnorm", mean = mean(x, na.rm = TRUE), sd = sd(x, na.rm = TRUE))$p.value
  } else {
    NA
  }
})

## Warning in ks.test.default(x[!is.na(x)], "pnorm", mean = mean(x, na.rm = TRUE),
## : ties should not be present for the one-sample Kolmogorov-Smirnov test
## Warning in ks.test.default(x[!is.na(x)], "pnorm", mean = mean(x, na.rm = TRUE),
## : ties should not be present for the one-sample Kolmogorov-Smirnov test
## Warning in ks.test.default(x[!is.na(x)], "pnorm", mean = mean(x, na.rm = TRUE),
## : ties should not be present for the one-sample Kolmogorov-Smirnov test
## Warning in ks.test.default(x[!is.na(x)], "pnorm", mean = mean(x, na.rm = TRUE),
## : ties should not be present for the one-sample Kolmogorov-Smirnov test
## Warning in ks.test.default(x[!is.na(x)], "pnorm", mean = mean(x, na.rm = TRUE),
## : ties should not be present for the one-sample Kolmogorov-Smirnov test
## Warning in ks.test.default(x[!is.na(x)], "pnorm", mean = mean(x, na.rm = TRUE),
## : ties should not be present for the one-sample Kolmogorov-Smirnov test

print(Kolmogorov)

##       estrato       preciom     areaconst  parqueaderos        banios 
## 1.571815e-314 9.661993e-206 3.477542e-242  0.000000e+00 4.289098e-292 
##  habitaciones 
##  0.000000e+00

Los valores p de la prueba Kolmogorov-Smirnov son pequeños es decir que todas las variables analizadas no siguen una distribución normal con un nivel de confianza muy alto (es decir, rechazamos la hipótesis nula de normalidad). Sin embargo, el aviso indica que hay empates en los datos, lo cual es un problema para esta prueba porque el test asume que los valores en la muestra son únicos.

Por lo anterior, se procede a manejar los empates añadiendo un ruido pequeño.

Kolmogorov = sapply(vivienda_new, function(x) {
  n <- sum(!is.na(x))
  if (n >= 3 && n <= 5000) {
    ks.test(jitter(x), "pnorm", mean = mean(x, na.rm = TRUE), sd = sd(x, na.rm = TRUE))
  } else if (n > 5000) {
    ks.test(jitter(x), "pnorm", mean = mean(x, na.rm = TRUE), sd = sd(x, na.rm = TRUE))
  } else {
    NA
  }
})
print(Kolmogorov)

##             estrato                                        
## statistic   0.133729                                       
## p.value     1.199315e-129                                  
## alternative "two-sided"                                    
## method      "Asymptotic one-sample Kolmogorov-Smirnov test"
## data.name   "jitter(x)"                                    
## exact       FALSE                                          
##             preciom                                        
## statistic   0.1683143                                      
## p.value     3.732583e-205                                  
## alternative "two-sided"                                    
## method      "Asymptotic one-sample Kolmogorov-Smirnov test"
## data.name   "jitter(x)"                                    
## exact       FALSE                                          
##             areaconst                                      
## statistic   0.1826687                                      
## p.value     1.555851e-241                                  
## alternative "two-sided"                                    
## method      "Asymptotic one-sample Kolmogorov-Smirnov test"
## data.name   "jitter(x)"                                    
## exact       FALSE                                          
##             parqueaderos                                   
## statistic   0.2110675                                      
## p.value     2.428664e-260                                  
## alternative "two-sided"                                    
## method      "Asymptotic one-sample Kolmogorov-Smirnov test"
## data.name   "jitter(x)"                                    
## exact       FALSE                                          
##             banios                                         
## statistic   0.157462                                       
## p.value     1.38897e-179                                   
## alternative "two-sided"                                    
## method      "Asymptotic one-sample Kolmogorov-Smirnov test"
## data.name   "jitter(x)"                                    
## exact       FALSE                                          
##             habitaciones                                   
## statistic   0.228242                                       
## p.value     0                                              
## alternative "two-sided"                                    
## method      "Asymptotic one-sample Kolmogorov-Smirnov test"
## data.name   "jitter(x)"                                    
## exact       FALSE

Verificamos de forma grafica la no - normalidad de las variables.

variables = c("preciom", "areaconst", "banios", "habitaciones", "parqueaderos")
for (var in variables) {
  x = vivienda[[var]]
  x = na.omit(x)
  qqnorm(x, main = paste("Gráfico Q-Q de", var))
  qqline(x, col = "red")
}

Asi pues, procedemos a imputar los valores faltantes con la mediana y posteriormnte a reemplazar las columnas de vivienda_new por las imputadas.

num_vars_new = vivienda_new[sapply(vivienda_new, is.numeric)] ## filtro
for (var in names(num_vars_new)) {
  num_vars_new[[var]][is.na(num_vars_new[[var]])] = median(num_vars_new[[var]], na.rm = TRUE)}
vivienda_new[sapply(vivienda_new, is.numeric)] = num_vars_new

Ningun dato faltante.

colSums(is.na(vivienda_new))

##      estrato      preciom    areaconst parqueaderos       banios habitaciones 
##            0            0            0            0            0            0

Procedemos a estandarizar las variables.

Con el fin de evitar que las variables que tiene una escala con valores más grandes afecten las estimaciones realizadas (sesgos) se realiza la estandarización de las variables antes de proceder a realizar el proceso de estimación de los componentes principales.

vivienda_num = vivienda_new[sapply(vivienda_new, is.numeric)]
viviendaZ = scale(vivienda_num)
head(viviendaZ)

##         estrato    preciom  areaconst parqueaderos     banios habitaciones
## [1,] -1.5876059 -0.5595266 -0.7339788   -0.8561040 -0.0779236    1.6410785
## [2,] -1.5876059 -0.3465092 -0.3841860   -0.8561040 -0.7782261   -0.4146749
## [3,] -1.5876059 -0.2552161  0.3153997    0.1313523 -0.7782261    0.2705762
## [4,] -0.6158454 -0.1030608  0.7351511    1.1188087  1.3226815   -0.4146749
## [5,]  0.3559152 -0.5290955 -0.5940617   -0.8561040 -0.7782261   -0.4146749
## [6,]  0.3559152 -0.5899576 -0.6150492   -0.8561040 -0.0779236   -0.4146749

Elección del número de componentes principales

res.pca = prcomp(viviendaZ)
fviz_eig(res.pca, addlabels = TRUE)

En este caso el primer componente principal explica el 56.1% de la variabilidad contenida en la base de datos y entre los dos primeros el 76.1%, lo cual indicaría que con solo una variable (CP1) que se obtiene mediante una combinación lineal de las variables se puede resumir gran parte de la variabilidad que contiene la base de datos.

fviz_pca_var(res.pca,
             col.var = "contrib", 
             gradient.cols = c("#FF7F00",  "#034D94"),
             repel = TRUE   ) 

Para entender mejor el grafico, tenemos que mirar las contribuciones, en tanto que tenemos 2 componenetes con unos pesos si se quiere, en PC1, la variable con mayor contribución es preciom con 0.23, seguido por banios con 0.221 y areaconst con 0.20. Esto significa que PC1 está más relacionado con el precio, los baños y el área construida. En PC2, la variable más importante es habitaciones 0.4629, lo que indica que esta componente refleja la variabilidad en el número de habitaciones.

contrib = res.pca$rotation^2
contrib

##                     PC1        PC2          PC3        PC4        PC5
## estrato      0.09321396 0.39877772 0.2168467376 0.03761240 0.20874592
## preciom      0.23093492 0.06323765 0.0003038743 0.06911761 0.03211104
## areaconst    0.20866747 0.03621701 0.0128800344 0.52659336 0.03383652
## parqueaderos 0.15509069 0.01585556 0.6070929374 0.19337113 0.01238919
## banios       0.22139252 0.02302822 0.1047302168 0.08107465 0.48091830
## habitaciones 0.09070044 0.46288384 0.0581461995 0.09223085 0.23199902
##                     PC6
## estrato      0.04480326
## preciom      0.60429491
## areaconst    0.18180561
## parqueaderos 0.01620048
## banios       0.08885608
## habitaciones 0.06403966

fviz_pca_ind(res.pca, col.ind = "gray")

casos1 = rbind(res.pca$x[426,1:2],res.pca$x[1258,1:2]) # CP1
rownames(casos1) = c("426","1258")
casos1 = as.data.frame(casos1)
casos2 = rbind(res.pca$x[5342,1:2], res.pca$x[6645,1:2]) # CP2
rownames(casos2) = c("5342","6645")
casos2 = as.data.frame(casos2)
fviz_pca_ind(res.pca, col.ind = "#DEDEDE", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")) +
  geom_point(data = casos1, aes(x = PC1, y = PC2), color = "red", size = 3) +
  geom_text(data = casos1, aes(x = PC1, y = PC2, label = rownames(casos1)), color = "red", vjust = -1) +
  geom_point(data = casos2, aes(x = PC1, y = PC2), color = "blue", size = 3) +
  geom_text(data = casos2, aes(x = PC1, y = PC2, label = rownames(casos2)), color = "blue", vjust = -1)

datos = rbind(vivienda_new[426,], # ok
              vivienda_new[1258,],
              vivienda_new[5342,],
              vivienda_new[6645,])

datos = as.data.frame(datos)
rownames(datos) = c("Vivienda 0426","Vivienda 1258","Vivienda 5342","Vivienda 6645")
datos

##               estrato preciom areaconst parqueaderos banios habitaciones
## Vivienda 0426       3     113        58            2      2            3
## Vivienda 1258       3     169       160            2      4            6
## Vivienda 5342       6    1200       296            2      0            0
## Vivienda 6645       6    1600       825            8      6            8

Los datos muestran cuatro viviendas con diferencias significativas en estrato socioeconómico, precio, área construida y características internas. Las viviendas 0426 y 1258 pertenecen al estrato 3, mientras que las viviendas 5342 y 6645 están en el estrato 6, lo que se refleja en sus valores de precio y tamaño. Las viviendas de estrato 6 tienen precios considerablemente más altos (1200M y 1600M) en comparación con las de estrato 3 (113M y 169M). Además, la diferencia en el área construida es notable, ya que las viviendas del estrato 6 cuentan con espacios significativamente más amplios, especialmente la vivienda 6645 con 825 m². En relación con el análisis de componentes principales (PCA), es probable que las viviendas de estrato 6 se sitúen en posiciones extremas dentro del espacio de los componentes principales debido a su alto precio y tamaño, mientras que la vivienda 5342 podría actuar como un punto atípico por la ausencia de habitaciones y baños.

grafica2 = fviz_pca_biplot(res.pca, repel = TRUE, col.var = "#034A94", col.ind = "#DEDEDE")
grafica2

Analisis de Conglomerados

head (viviendaZ)

##         estrato    preciom  areaconst parqueaderos     banios habitaciones
## [1,] -1.5876059 -0.5595266 -0.7339788   -0.8561040 -0.0779236    1.6410785
## [2,] -1.5876059 -0.3465092 -0.3841860   -0.8561040 -0.7782261   -0.4146749
## [3,] -1.5876059 -0.2552161  0.3153997    0.1313523 -0.7782261    0.2705762
## [4,] -0.6158454 -0.1030608  0.7351511    1.1188087  1.3226815   -0.4146749
## [5,]  0.3559152 -0.5290955 -0.5940617   -0.8561040 -0.7782261   -0.4146749
## [6,]  0.3559152 -0.5899576 -0.6150492   -0.8561040 -0.0779236   -0.4146749

euclides = dist(viviendaZ, method = "euclidean")
head(euclides)

## [1] 2.210040 2.130157 3.660481 2.917931 2.831690 2.580722

manhatan = dist(viviendaZ, method = "manhattan")
head(manhatan)

## [1] 3.318866 4.411950 8.328627 4.869925 4.148635 4.932491

Ahora procedemos a esclar la matriz con el fin de mejorar la manipulación.

viviendaZDATA = as.data.frame(viviendaZ)
head(viviendaZDATA)

##      estrato    preciom  areaconst parqueaderos     banios habitaciones
## 1 -1.5876059 -0.5595266 -0.7339788   -0.8561040 -0.0779236    1.6410785
## 2 -1.5876059 -0.3465092 -0.3841860   -0.8561040 -0.7782261   -0.4146749
## 3 -1.5876059 -0.2552161  0.3153997    0.1313523 -0.7782261    0.2705762
## 4 -0.6158454 -0.1030608  0.7351511    1.1188087  1.3226815   -0.4146749
## 5  0.3559152 -0.5290955 -0.5940617   -0.8561040 -0.7782261   -0.4146749
## 6  0.3559152 -0.5899576 -0.6150492   -0.8561040 -0.0779236   -0.4146749

Metodo del codo

codo = sapply(1:10, function(k) { 
  kmeans(res.pca$x, centers = k, nstart = 25)$tot.withinss 
})

plot(1:10, codo, type = "b", pch = 19, frame = FALSE, 
     xlab = "Número de clusters K", 
     ylab = "Suma total de cuadrados")

Bajo lo anterior grafica procedemos a determinar el número óptimo de clusters en un análisis bajo el algortimo de Hierarchical Clustering on Principal Components. Donde aparanetmete se muestra una gran caida en k=4, es decir procedemos a trabajor sobre 4 clusters. Importantes mencionar que utilizaremos 4 cluester para cualquier algoritmo usado de ahora en adlenate.

HCPC (Hierarchical Clustering on Principal Components)

res.pca = PCA(vivienda_new, scale.unit = TRUE, ncp = 5, graph = FALSE)
res.HCPC = HCPC(res.pca, nb.clust = 4, graph = TRUE)

# Cluster
res.HCPC$desc.var

## 
## Link between the cluster variable and the quantitative variables
## ================================================================
##                   Eta2 P-value
## estrato      0.5288389       0
## preciom      0.7069952       0
## areaconst    0.5672215       0
## parqueaderos 0.5093225       0
## banios       0.5910463       0
## habitaciones 0.5737458       0
## 
## Description of each cluster by quantitative variables
## =====================================================
## $`1`
##                 v.test Mean in category Overall mean sd in category Overall sd
## parqueaderos -42.45273         1.375094     1.866979      0.4939019   1.012642
## habitaciones -42.58930         2.894050     3.605143      0.7218287   1.459231
## estrato      -49.98811         4.045192     4.633742      0.7690775   1.028998
## areaconst    -50.96314        91.569008   174.916216     47.4280786 142.933164
## preciom      -56.76058       220.461461   433.866979     85.1068266 328.591948
## banios       -63.16287         2.079337     3.111271      0.6375420   1.427868
##              p.value
## parqueaderos       0
## habitaciones       0
## estrato            0
## areaconst          0
## preciom            0
## banios             0
## 
## $`2`
##                 v.test Mean in category Overall mean sd in category Overall sd
## estrato      48.273175         5.466881     4.633742      0.5786298   1.028998
## banios       17.212217         3.523485     3.111271      0.8877995   1.427868
## preciom      15.775368       520.809715   433.866979    193.0064447 328.591948
## parqueaderos  6.082836         1.970293     1.866979      0.5130969   1.012642
## habitaciones -8.896801         3.387395     3.605143      0.7585097   1.459231
##                   p.value
## estrato      0.000000e+00
## banios       2.150387e-66
## preciom      4.597135e-56
## parqueaderos 1.180749e-09
## habitaciones 5.747923e-19
## 
## $`3`
##                 v.test Mean in category Overall mean sd in category Overall sd
## habitaciones  62.99430         6.405660     3.605143      1.5847840   1.459231
## banios        30.69502         4.446541     3.111271      1.4354648   1.427868
## areaconst     28.18009       297.628564   174.916216    126.3538828 142.933164
## estrato      -23.29170         3.903564     4.633742      0.8511479   1.028998
##                    p.value
## habitaciones  0.000000e+00
## banios       6.633339e-207
## areaconst    1.025691e-174
## estrato      5.380717e-120
## 
## $`4`
##                v.test Mean in category Overall mean sd in category Overall sd
## preciom      67.17567      1131.373602   433.866979    353.7091052 328.591948
## parqueaderos 60.54176         3.804251     1.866979      1.5101239   1.012642
## areaconst    54.61098       421.572629   174.916216    200.9327827 142.933164
## banios       44.85998         5.135347     3.111271      1.2471996   1.427868
## estrato      33.21172         5.713647     4.633742      0.5502628   1.028998
## habitaciones 17.05358         4.391499     3.605143      1.3581206   1.459231
##                    p.value
## preciom       0.000000e+00
## parqueaderos  0.000000e+00
## areaconst     0.000000e+00
## banios        0.000000e+00
## estrato      7.292937e-242
## habitaciones  3.287603e-65

Interpretación de Clusters en el Análisis de Viviendas

El análisis de clusters ha permitido segmentar las viviendas en cuatro grupos diferenciados según sus características estructurales y socioeconómicas. A continuación, se presenta una interpretación de cada cluster:

Cluster 1: Viviendas de menor tamaño y menor costo

Este cluster agrupa viviendas con un menor número de habitaciones, baños y parqueaderos. Además, presentan áreas construidas reducidas y un precio significativamente menor al promedio. El estrato socioeconómico predominante en este grupo es inferior al promedio general.

Cluster 2: Viviendas de estrato alto con precios elevados

Las viviendas en este grupo pertenecen a un estrato socioeconómico superior al promedio. Tienen un mayor número de baños y un precio considerablemente más alto que el promedio. Sin embargo, la cantidad de habitaciones y parqueaderos no es significativamente mayor.

Cluster 3: Viviendas amplias con múltiples habitaciones

Las viviendas de este cluster se caracterizan por tener un mayor número de habitaciones y baños en comparación con el promedio. Además, poseen un área construida significativamente mayor. Sin embargo, el estrato socioeconómico de estas viviendas es inferior al promedio.

Cluster 4: Viviendas exclusivas de alto costo

Este cluster agrupa viviendas con los precios más elevados, así como la mayor cantidad de parqueaderos y baños. También poseen áreas construidas amplias y pertenecen a los estratos socioeconómicos más altos.

En conclusión, el análisis de clusters permite identificar patrones en la distribución y características de las viviendas, facilitando una mejor toma de decisiones en materia de planeación urbana y desarrollo inmobiliario.

Se observa la cantidad de registros que quedaron en cada cluster HCPC:

cluster_asignaciones = res.HCPC$data.clust$clust
table(cluster_asignaciones)

## cluster_asignaciones
##    1    2    3    4 
## 3983 2491  954  894

<p>Por otro lado el cluster con mayor elementos es el cluster 1 con 3983 registros. </p>

Modelo KMEANS

set.seed(0)
modelo_kmeans = kmeans(viviendaZ, centers = 4, nstart = 10)  
viviendaZ = data.frame(viviendaZ,
                             modelo_kmeans$cluster) 
aggregate(vivienda_new,
          by = list(viviendaZ$modelo_kmeans.cluster),
          FUN = median)  

##   Group.1 estrato preciom areaconst parqueaderos banios habitaciones
## 1       1       6    1150       367            4      5            4
## 2       2       4     216        78            1      2            3
## 3       3       4     415       282            2      4            6
## 4       4       6     470       153            2      3            3

Cluster 1 representa viviendas de alto estrato (6), con precios elevados (1150 millones), áreas grandes (367 m²) y más parqueaderos y baños en comparación con los otros grupos.

Cluster 2 corresponde a viviendas más económicas (216 millones), con menor área (78 m²) y menos comodidades, como parqueaderos y baños.

Cluster 3 parece ser un segmento intermedio, con precios y áreas moderadas, que se alejan de los extremos de los otros grupos.

Cluster 4 también pertenece al estrato 6, pero con precios y áreas menores en comparación con el Cluster 1, lo que indica que, aunque es de alto estrato, ofrece viviendas más accesibles o con menos lujos.

Se observa la cantidad de registros que quedaron en cada cluster KMEANS:

table(viviendaZ$modelo_kmeans.cluster)

## 
##    1    2    3    4 
##  898 3965  949 2510

El cluster 2 en el modelo Kmeas es el que tiene una mayor cantidad de datos.

Creamos una visualización de los 4 clústeres generados por el modelo K-Means en el conjunto de datos viviendaZ.

fviz_cluster(list(data = viviendaZ[,1:6], 
                  cluster = viviendaZ$modelo_kmeans.cluster),
             palette = c("#cd1076",  "#8deeee", "#00FF00", "#FFD700"),  # 4 colores para 4 clusters
             ellipse.type = "convex", 
             repel = FALSE, 
             show.clust.cent = FALSE, 
             ggtheme = theme_minimal())

Analisis de Correspondencia

vivienda_nov = vivienda[,c(2,4,10,11)]
head(vivienda_nov)

## # A tibble: 6 × 4
##   zona         estrato tipo        barrio     
##   <chr>          <dbl> <chr>       <chr>      
## 1 Zona Oriente       3 Casa        20 de julio
## 2 Zona Oriente       3 Casa        20 de julio
## 3 Zona Oriente       3 Casa        20 de julio
## 4 Zona Sur           4 Casa        3 de julio 
## 5 Zona Norte         5 Apartamento acopi      
## 6 Zona Norte         5 Apartamento acopi

Para desarrollar analsis de correspondencia es necesario transformar las variables categóricas en factores para que puedan ser utilizadas correctamente en el análisis. A

vivienda_AC = vivienda_nov
vivienda_AC$estrato = as.factor(vivienda_AC$estrato)
vivienda_AC$zona = as.factor(vivienda_AC$zona)
vivienda_AC$tipo = as.factor(vivienda_AC$tipo)
vivienda_AC$barrio  = as.factor(vivienda_AC$barrio)

Antes que nada debemos contar los valores faltantes (NA) por columna en las categoricas, este proceso lo hicimos inicalemnte pero con las numericas.

colSums(is.na(vivienda_AC))

##    zona estrato    tipo  barrio 
##       3       3       3       3

Como tebemos datos falatntes y son variables categoricas remplzamos los faltantes con la moda.

get_mode = function(x) {
  uniqx = unique(x)
  uniqx[which.max(tabulate(match(x, uniqx)))]
}
vivienda_AC$estrato[is.na(vivienda_AC$estrato)] = get_mode(vivienda_AC$estrato)
vivienda_AC$zona[is.na(vivienda_AC$zona)] = get_mode(vivienda_AC$zona)
vivienda_AC$tipo[is.na(vivienda_AC$tipo)] = get_mode(vivienda_AC$tipo)
vivienda_AC$barrio [is.na(vivienda_AC$barrio )] = get_mode(vivienda_AC$barrio )

Se ejecuta la prueba Chi-Cuadrado para validar el supuesto de independencia entre los datos

# Independencia entre estrato y zona
tabla1 = table(vivienda_AC$estrato, vivienda_AC$zona)
resultado_chicuadrado_estrato_zona = chisq.test(tabla1)
resultado_chicuadrado_estrato_zona

## 
##  Pearson's Chi-squared test
## 
## data:  tabla1
## X-squared = 3831.8, df = 12, p-value < 2.2e-16

Dado que el valor p es extremadamente pequeño (mucho menor que 0.05), se rechaza la hipótesis nula, que establece que no hay asociación entre las variables.

# Independencia entre estrato y tipo
tabla2 = table(vivienda_AC$estrato, vivienda_AC$tipo)
resultado_chicuadrado_estrato_tipo = chisq.test(tabla2)
resultado_chicuadrado_estrato_tipo

## 
##  Pearson's Chi-squared test
## 
## data:  tabla2
## X-squared = 224.64, df = 3, p-value < 2.2e-16

Existe una relación significativa entre las variables. Sin embargo, si las frecuencias esperadas son muy pequeñas en algunas celdas, puede ser más adecuado considerar otras pruebas estadísticas, como la prueba exacta de Fisher si las tablas son pequeñas o si las celdas tienen frecuencias bajas.

# Independencia entre zona y tipo
tabla3 = table(vivienda_AC$zona, vivienda_AC$tipo)
resultado_chicuadrado_zona_tipo = chisq.test(tabla3)
resultado_chicuadrado_zona_tipo

## 
##  Pearson's Chi-squared test
## 
## data:  tabla3
## X-squared = 690.79, df = 4, p-value < 2.2e-16

Relación estadísticamente significativa entre las variables analizadas, ya que el valor p es mucho menor que el umbral común de 0.05.

Análisis de Correspondencias Múltiples (MCA)

Ahora bien, incilamente existia un problema de baja dimencioanlidad, lo que imposibilitaba ver de fomra grafica las relaciones, por lo anterior para poder visualizar las relaciones entre varias variables categóricas se utiliza el MCA.

mca_data = vivienda_AC %>% select("estrato", "zona", "tipo") %>% mutate(across(everything(),as.factor))
resultado_mca = MCA(mca_data, graph = FALSE)
fviz_mca_var(
  resultado_mca,
  repel= TRUE)

Ahora graficamos el porcentaje de varianza explicado por cada componente principal en un análisis de componentes principales (PCA), en este caso empezamos con las relación estrato y zona (tabla1), esto en funcion a que fue la unica dos variables que se dejaron graficar con la función CA, tal como se iliustra a continuación.

tabla1a = CA(tabla1) # Grafica 1

fviz_screeplot(tabla1a, addlabels = TRUE, ylim = c(0, 80))+ggtitle("")+ # Grafica 2
ylab("Porcentaje de varianza explicado") + xlab("Ejes")

Tanyo PC1 y PC2 explican 97.7% de la variabilidad en los datos, lo que significa que se puede reducir el número de dimensiones en el análisis a estos dos componentes sin perder mucha información.

Ahora graficaremos la varianza explicada pero bajo el MCA, como se ilusra a continuación:

fviz_screeplot(resultado_mca, addlabels = TRUE, ylim = c(0, 80))+ggtitle("")+
ylab("Porcentaje de varianza explicado") + xlab("Ejes")

En este caso no hay muchi que decir,las primeras dimensiones (especialmente las dimensiones 1 y 2) capturan una gran parte de la información. Luego, a medida que avanzas a dimensiones superiores, la cantidad de varianza explicada disminuye.

Visualización de resultados

head(vivienda)

## # A tibble: 6 × 13
##      id zona    piso  estrato preciom areaconst parqueaderos banios habitaciones
##   <dbl> <chr>   <chr>   <dbl>   <dbl>     <dbl>        <dbl>  <dbl>        <dbl>
## 1  1147 Zona O… <NA>        3     250        70            1      3            6
## 2  1169 Zona O… <NA>        3     320       120            1      2            3
## 3  1350 Zona O… <NA>        3     350       220            2      2            4
## 4  5992 Zona S… 02          4     400       280            3      5            3
## 5  1212 Zona N… 01          5     260        90            1      2            3
## 6  1724 Zona N… 01          5     240        87            1      3            3
## # ℹ 4 more variables: tipo <chr>, barrio <chr>, longitud <dbl>, latitud <dbl>

Combinamos alos nuevos dataframes las variables de las corrdenada spara poder ilutrasr a tarves de leaflet los resultados.

vivienda_new$latitud = vivienda$latitud[match(rownames(vivienda_new), rownames(vivienda))]
vivienda_new$longitud = vivienda$longitud[match(rownames(vivienda_new), rownames(vivienda))]
vivienda_AC$latitud = vivienda$latitud[match(rownames(vivienda_new), rownames(vivienda))]
vivienda_AC$longitud = vivienda$longitud[match(rownames(vivienda_new), rownames(vivienda))]
vivienda_new$cluster = res.HCPC$data.clust$clust

Mapas

#install.packages("sf")
#install.packages("leaflet")
library(leaflet)

## Warning: package 'leaflet' was built under R version 4.4.2

library(sf)

## Warning: package 'sf' was built under R version 4.4.2

## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE

pal = colorFactor(palette = "Set1", domain = vivienda_new$estrato)
leaflet(data = na.omit(vivienda_new)) %>%
  addTiles() %>%  
  addCircleMarkers(
    lng = ~longitud, 
    lat = ~latitud, 
    color = ~pal(estrato),  
    radius = 3, 
    popup = ~paste("Estrato:", estrato),
    stroke = FALSE, 
    fillOpacity = 0.8
  ) %>%
  addLegend(
    position = "bottomright",
    title = "Estratos de vivienda",
    pal = pal, 
    values = vivienda_new$estrato,
    opacity = 1
  )

library(leaflet)
pal = colorFactor(palette = "Set2", domain = vivienda_new$cluster)
leaflet(data = na.omit(vivienda_new)) %>%
  addTiles() %>%  
  addCircleMarkers(
    lng = ~longitud, 
    lat = ~latitud, 
    color = ~pal(cluster),  
    radius = 3, 
    popup = ~paste("Cluster:", cluster),
    stroke = FALSE, 
    fillOpacity = 0.8
  ) %>%
  addLegend(
    position = "bottomright",
    title = "Clusters de Vivienda",
    pal = pal,  
    values = vivienda_new$cluster,
    opacity = 1
  )

library(leaflet)
pal <- colorFactor(palette = "Set1", domain = vivienda_AC$tipo)
leaflet(data = na.omit(vivienda_AC)) %>%
  addTiles() %>%
  addCircleMarkers(
    lng = ~longitud, 
    lat = ~latitud, 
    color = ~pal(tipo), 
    radius = 3, 
    popup = ~paste("Tipo de vivienda:", tipo),
    stroke = FALSE, 
    fillOpacity = 0.8
  ) %>%
  addLegend(
    position = "bottomright",
    title = "Tipo de Vivienda",
    pal = pal,  
    values = vivienda_AC$tipo,
    opacity = 1
  )