1 Introducción

Con la intención de generar un análisis reproducible, fácil y entendible, se elabora el siguiente código para realizar análisis de datos. Se recupera la información de Kaggle.

Sin antes, se ha de mencionar que, al ser un análisis sencillo no se profundiza teóricamente en los supuestos de algebra matricial, ni de las librerias a utilizar, por lo que es responsabilidad de la persona lectora de este código buscar la información complementaría.

Se agradece cualquier tipo de comentario y se advierte que el código presente se mantiene en actualización.

1.1 Putos clave para la discusión del análisis

  • Objetivo
  • Supuestos
  • Validaciones
  • Selección de componentes principales

1.1.1 Librerías de trabajo

A continuación, se enlistas las librerías a utilizar, mismas que serán mencionadas a lo largo del código para poder identificar el uso que tienen.

El siguiente código instala las paqueterías necesarias en caso de no contar previamente con ellas y posteriormente las carga.

if(!require('pacman')) install.packages('pacman')
pacman::p_load(tidyverse,#manipulación de datos
               janitor,#limpieza y tablas de frecuencia
               tidymodels,#modelos
               naniar,#NA's
               GGally,#Grafico distribuciones
               corrplot,#Grafico correlacion
               RColorBrewer,#Paletas de colores
               FactoMineR,#PCA
               caret,#PCA
               factoextra,#graficos pca
               fpc,#Cluster
               embed,#UMAP
               highlight,#RMarkDown
               tidytext)#Texto en gráfico
#Notación cientifica
options(scipen = 999)

1.1.2 Descargamos datos

Los datos provienen de las bases de Kaggle, para este caso la data corresponde a información sobre los clientes de un supermercado.

#Cargamos datos
customers <- read.delim("data/raw/marketing_campaign.csv", stringsAsFactors = FALSE)

1.1.3 Objetivos del Modelo

Crear clusters de clientes dadas las caracteristicas que los definen.

1.1.4 Estructura de Datos

Clientes

  • ID: Customer’s unique identifier
  • Year_Birth: Customer’s birth year
  • Education: Customer’s education level
  • Marital_Status: Customer’s marital status
  • Income: Customer’s yearly household income
  • Kidhome: Number of children in customer’s household
  • Teenhome: Number of teenagers in customer’s household
  • Dt_Customer: Date of customer’s enrollment with the company
  • Recency: Number of days since customer’s last purchase
  • Complain: 1 if customer complained in the last 2 years, 0 otherwise

Productos

  • MntWines: Amount spent on wine in last 2 years
  • MntFruits: Amount spent on fruits in last 2 years
  • MntMeatProducts: Amount spent on meat in last 2 years
  • MntFishProducts: Amount spent on fish in last 2 years
  • MntSweetProducts: Amount spent on sweets in last 2 years
  • MntGoldProds: Amount spent on gold in last 2 years

Promoción

  • NumDealsPurchases: Number of purchases made with a discount
  • AcceptedCmp1: 1 if customer accepted the offer in the 1st campaign, 0 otherwise
  • AcceptedCmp2: 1 if customer accepted the offer in the 2nd campaign, 0 otherwise
  • AcceptedCmp3: 1 if customer accepted the offer in the 3rd campaign, 0 otherwise
  • AcceptedCmp4: 1 if customer accepted the offer in the 4th campaign, 0 otherwise
  • AcceptedCmp5: 1 if customer accepted the offer in the 5th campaign, 0 otherwise
  • Response: 1 if customer accepted the offer in the last campaign, 0 otherwise

Lugar

  • NumWebPurchases: Number of purchases made through the company’s web site
  • NumCatalogPurchases: Number of purchases made using a catalogue
  • NumStorePurchases: Number of purchases made directly in stores
  • NumWebVisitsMonth: Number of visits to company’s web site in the last month.
#Estructura de datos
glimpse(customers)
## Rows: 2,240
## Columns: 29
## $ ID                  <int> 5524, 2174, 4141, 6182, 5324, 7446, 965, 6177, 485…
## $ Year_Birth          <int> 1957, 1954, 1965, 1984, 1981, 1967, 1971, 1985, 19…
## $ Education           <chr> "Graduation", "Graduation", "Graduation", "Graduat…
## $ Marital_Status      <chr> "Single", "Single", "Together", "Together", "Marri…
## $ Income              <int> 58138, 46344, 71613, 26646, 58293, 62513, 55635, 3…
## $ Kidhome             <int> 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1,…
## $ Teenhome            <int> 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1,…
## $ Dt_Customer         <chr> "04-09-2012", "08-03-2014", "21-08-2013", "10-02-2…
## $ Recency             <int> 58, 38, 26, 26, 94, 16, 34, 32, 19, 68, 11, 59, 82…
## $ MntWines            <int> 635, 11, 426, 11, 173, 520, 235, 76, 14, 28, 5, 6,…
## $ MntFruits           <int> 88, 1, 49, 4, 43, 42, 65, 10, 0, 0, 5, 16, 61, 2, …
## $ MntMeatProducts     <int> 546, 6, 127, 20, 118, 98, 164, 56, 24, 6, 6, 11, 4…
## $ MntFishProducts     <int> 172, 2, 111, 10, 46, 0, 50, 3, 3, 1, 0, 11, 225, 3…
## $ MntSweetProducts    <int> 88, 1, 21, 3, 27, 42, 49, 1, 3, 1, 2, 1, 112, 5, 1…
## $ MntGoldProds        <int> 88, 6, 42, 5, 15, 14, 27, 23, 2, 13, 1, 16, 30, 14…
## $ NumDealsPurchases   <int> 3, 2, 1, 2, 5, 2, 4, 2, 1, 1, 1, 1, 1, 3, 1, 1, 3,…
## $ NumWebPurchases     <int> 8, 1, 8, 2, 5, 6, 7, 4, 3, 1, 1, 2, 3, 6, 1, 7, 3,…
## $ NumCatalogPurchases <int> 10, 1, 2, 0, 3, 4, 3, 0, 0, 0, 0, 0, 4, 1, 0, 6, 0…
## $ NumStorePurchases   <int> 4, 2, 10, 4, 6, 10, 7, 4, 2, 0, 2, 3, 8, 5, 3, 12,…
## $ NumWebVisitsMonth   <int> 7, 5, 4, 6, 5, 6, 6, 8, 9, 20, 7, 8, 2, 6, 8, 3, 8…
## $ AcceptedCmp3        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,…
## $ AcceptedCmp4        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ AcceptedCmp5        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,…
## $ AcceptedCmp1        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,…
## $ AcceptedCmp2        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ Complain            <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ Z_CostContact       <int> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,…
## $ Z_Revenue           <int> 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11…
## $ Response            <int> 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0,…
#NA's
customers %>% 
  summarise(across(everything(), ~ sum(is.na(.)))) %>% 
  pivot_longer(everything(), names_to = "variables", values_to = "na_count")
#library(naniar)
miss_var_summary(customers)
  • Las fechas no estan en el formato adecuado
  • Valores perdidos en la variable income
  • Codificar variables categoricas
  • Reclasificar para reducir variables categoricas

1.1.5 Limpieza

dta <- customers

#Fecha
dta$Dt_Customer <- dta %>% mutate(Dt_Customer = as.Date(dmy(Dt_Customer)))

#Grupos de Edad
  # Genera los puntos de corte (breaks)
breaks <- seq(min(dta$Year_Birth), max(dta$Year_Birth) + 5, by = 5)
  # Genera las etiquetas para los intervalos
labels <- paste(breaks[-length(breaks)], breaks[-1] - 1, sep = "-")
  # Crea una nueva columna "Age_Group" con los grupos de 5 años
dta$Age_Group <- cut(
  dta$Year_Birth, 
  breaks = breaks, 
  right = FALSE, 
  labels = labels)

#Reducir categorias en estado civil
dta <- dta %>% mutate(Marital_Status = replace(Marital_Status, Marital_Status == "Divorced" | Marital_Status == "Widow" | Marital_Status == "Alone" | Marital_Status == "Absurd" | Marital_Status == "YOLO", "Single"))
dta <- dta %>% mutate(Marital_Status = replace(Marital_Status, Marital_Status == "Together" | Marital_Status == "Married", "Taken"))
dta$Marital_Status <- factor(dta$Marital_Status)

#Reducir categorias de educación
dta <- dta %>% mutate(Education = replace(Education, Education == "2n Cycle", "non-graduate"))
dta$Education <- factor(dta$Education)

#Renombrar
dta <- dta %>% rename(wines = MntWines, fruits = MntFruits, meat = MntMeatProducts, fish = MntFishProducts, sweet = MntSweetProducts, gold = MntGoldProds )

#Crear una nueva variable
dta <- dta %>% mutate(Total_spent = wines + fruits + meat + fish + sweet + gold)

#Quitar variables dedundantes 
dta <- dta %>% select(- ID, - Year_Birth, - Dt_Customer, - Z_CostContact, - Z_Revenue)

#Formateo de base, quitar NA's y atipicos
dta <- dta %>%
  janitor::clean_names() %>% 
  na.omit() %>% 
  filter(income < 200000)
# Revisamos base de datos limpia
glimpse(dta)
## Rows: 2,215
## Columns: 26
## $ education             <fct> Graduation, Graduation, Graduation, Graduation, …
## $ marital_status        <fct> Single, Single, Taken, Taken, Taken, Taken, Sing…
## $ income                <int> 58138, 46344, 71613, 26646, 58293, 62513, 55635,…
## $ kidhome               <int> 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, …
## $ teenhome              <int> 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, …
## $ recency               <int> 58, 38, 26, 26, 94, 16, 34, 32, 19, 68, 59, 82, …
## $ wines                 <int> 635, 11, 426, 11, 173, 520, 235, 76, 14, 28, 6, …
## $ fruits                <int> 88, 1, 49, 4, 43, 42, 65, 10, 0, 0, 16, 61, 2, 1…
## $ meat                  <int> 546, 6, 127, 20, 118, 98, 164, 56, 24, 6, 11, 48…
## $ fish                  <int> 172, 2, 111, 10, 46, 0, 50, 3, 3, 1, 11, 225, 3,…
## $ sweet                 <int> 88, 1, 21, 3, 27, 42, 49, 1, 3, 1, 1, 112, 5, 1,…
## $ gold                  <int> 88, 6, 42, 5, 15, 14, 27, 23, 2, 13, 16, 30, 14,…
## $ num_deals_purchases   <int> 3, 2, 1, 2, 5, 2, 4, 2, 1, 1, 1, 1, 3, 1, 1, 3, …
## $ num_web_purchases     <int> 8, 1, 8, 2, 5, 6, 7, 4, 3, 1, 2, 3, 6, 1, 7, 3, …
## $ num_catalog_purchases <int> 10, 1, 2, 0, 3, 4, 3, 0, 0, 0, 0, 4, 1, 0, 6, 0,…
## $ num_store_purchases   <int> 4, 2, 10, 4, 6, 10, 7, 4, 2, 0, 3, 8, 5, 3, 12, …
## $ num_web_visits_month  <int> 7, 5, 4, 6, 5, 6, 6, 8, 9, 20, 8, 2, 6, 8, 3, 8,…
## $ accepted_cmp3         <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, …
## $ accepted_cmp4         <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ accepted_cmp5         <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, …
## $ accepted_cmp1         <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, …
## $ accepted_cmp2         <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ complain              <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ response              <int> 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, …
## $ age_group             <fct> 1953-1957, 1953-1957, 1963-1967, 1983-1987, 1978…
## $ total_spent           <int> 1617, 27, 776, 53, 422, 716, 590, 169, 46, 49, 6…
#library(corrplot)
#library(RColorBrewer)
#Variables Númericas
col.numeric <- dta %>% select(where(is.numeric))
#Matriz de Correlación
cor_matrix <- cor(col.numeric, use = "complete.obs")
#Gráfico de correlación
corrplot(
  cor_matrix, 
  tl.col = "black",       # Color de las etiquetas de las variables
  col = colorRampPalette(brewer.pal(11, "PuOr"))(200)  # Paleta de colores
)

#Gráficos de distribución de las variables
col.numeric %>% 
  # Convertir la base de datos en formato largo
  pivot_longer(cols = everything(),#Mover todas las variables
               names_to = "variable",#Pasarlas a una sola
               values_to = "valor") %>% #Pasar las observaciones a una sola
  ggplot(aes(x = valor)) +
  geom_histogram(bins = 25) +
  facet_wrap(~ variable, scales = "free") +
  theme_bw()

#Copia de la base
datos <- dta
#Codificamos variables categoricas a numericas
# Encoding the categorical features to numeric
# Convert education levels to numeric values using case_when
datos$education_numeric <- case_when(
  datos$education == "Basic" ~ 1,
  datos$education == "Graduation" ~ 2,
  datos$education == "Master" ~ 3,
  datos$education == "non-graduate" ~ 4,
  datos$education == "PhD" ~ 5,
  TRUE ~ NA # Handle missing values
)


datos$marital_status_numeric <- case_when(
 datos$marital_status == "Single" ~ 6,
 datos$marital_status == "Taken" ~ 7,
  TRUE ~ NA # Manejar valores faltantes
)

# Assuming "datos" is your data frame
datos$age_group_numeric <- case_when(
  datos$age_group == "1893-1897" ~ 8,
  datos$age_group == "1898-1902" ~ 9,
  datos$age_group == "1903-1907" ~ 10,
  datos$age_group == "1908-1912" ~ 11,
  datos$age_group == "1913-1917" ~ 12,
  datos$age_group == "1918-1922" ~ 13,
  datos$age_group == "1923-1927" ~ 14,
  datos$age_group == "1928-1932" ~ 15,
  datos$age_group == "1933-1937" ~ 16,
  datos$age_group == "1938-1942" ~ 17,
  datos$age_group == "1943-1947" ~ 18,
  datos$age_group == "1948-1952" ~ 19,
  datos$age_group == "1953-1957" ~ 20,
  datos$age_group == "1958-1962" ~ 21,
  datos$age_group == "1963-1967" ~ 22,
  datos$age_group == "1968-1972" ~ 23,
  datos$age_group == "1973-1977" ~ 24,
  datos$age_group == "1978-1982" ~ 25,
  datos$age_group == "1983-1987" ~ 26,
  datos$age_group == "1988-1992" ~ 27,
  datos$age_group == "1993-1997" ~ 28,
  TRUE ~ NA # Handle invalid age_group values
)
#Con los datos filtrados, creamos una base para trabajar
datos <- datos %>%
  select(-c("education", "marital_status", "age_group"))

2 Modelos

2.1 Tidymodels

Usamos el componente basico de tidymodels para realizare el análisis PCA, posteriormente usamos otras librerias. Hasta el momento desconozco como aplicar PCA de otras librerias dentro de tidymodels.

2.1.1 1. Receta

Primero, debemos decirle recipe() que está pasando con nuestro modelo (observe la fórmula sin resultado ) y qué datos estamos usando.

#library(tidymodels)
pca_rec <- recipe(~.,data = datos) %>%
  #Normalizamos datos
  step_normalize(all_predictors()) %>%
  #Convertinos en Componentes Principales
  step_pca(all_predictors())

Antes de su uso estos pasos se han definido pero no se han ejecutado ni implementado. La prep()función es donde todo se evalúa.

2.1.2 2. Preparación

Para una receta con al menos una operación de preprocesamiento, estime los parámetros requeridos a partir de un conjunto de entrenamiento que luego se puede aplicar a otros conjuntos de datos.

pca_prep <- prep(pca_rec)
pca_prep

2.1.3 3. Datos procesado ~ lo que cocinamos

#Resultados del cocinamiento
juice(pca_prep)

2.1.4 4.Visualización de Resultados

#Visualización de Resultados
tidied_pca <- tidy(pca_prep, 2)

tidied_pca %>%
  filter(component %in% paste0("PC", 1:5)) %>%
  mutate(component = fct_inorder(component)) %>%
  ggplot(aes(value, terms, fill = terms)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~component, nrow = 1) +
  labs(y = NULL)

#library(tidytext)

tidied_pca %>%
  filter(component %in% paste0("PC", 1:4)) %>%
  group_by(component) %>%
  top_n(8, abs(value)) %>%
  ungroup() %>%
  mutate(terms = reorder_within(terms, abs(value), component)) %>%
  ggplot(aes(abs(value), terms, fill = value > 0)) +
  geom_col() +
  facet_wrap(~component, scales = "free_y") +
  scale_y_reordered() +
  labs(
    x = "Absolute value of contribution",
    y = NULL, fill = "Positive?"
  )

  • El componente uno (C1) se explica mayormente por las compras de productos por catalogo y en tienda, así como el ingreso y el total de compra.
  • El C2 se explica mayormente por el tipo de compra, el grupo de edad y los años de educación.
  • El C3 se explica por el número de respuesta a la oferta en las campañas, así como la compra de frutas y pescado.

2.2 FactoMiner

2.2.1 1. Normalizar los datos

#library(FactoMineR)
#library(caret)
#Restamos media y dividimos entre varianza = Valores Z para normalizar
#train_data.standard <- preProcess(datos, method = c("center", "scale"))
#train_data.norm <- predict(train_data.standard, train_data)

#Alternativamente puedo escalar los datos con scale
dat <- scale(datos)

#Corremos PCA.
pca1 <- PCA(dat, graph = FALSE)

#Corremos PCA nativo
acp <- prcomp(dat)

Numero de componentes a retener: ScreePlot

El Screenplot nos ayuda visualmente, a través de los eigenvalores, la varianza explicada por cada componente.

# Scree plot para determinar el número de componentes a retener
fviz_eig(pca1, addlabels = TRUE, ylim = c(0, 50))

  • El 29.7% de la varianza de los datos esta explicada en el CP1

  • El CP2 y CP3, explican casi la misma varianzación del modelo.

  • Criterio del Codo con Eigenvalores

  • Criterio del eigenvalor mayor a 1

#library(factoextra)
#Criterio valores > 1
fviz_eig(pca1, choice = "eigenvalue")

  • El CP 1 y 2, explican en buena proporcion el modelo; ligeramente el cp3 lo hace.

¿Cuántos grupos conservamos?

#Donde estan representados nuestros puntos en las dimensiones del PCA
fviz_pca_ind(pca1,
             col.ind = "cos2", #Calidad de la representación 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07","#008f39"))

#Scores después de la rotación
#Nos dice como contribuyen de forma positiva o negativa las variables
#En donde estan las variables respecto al PCA
fviz_pca_var(pca1, 
             axes = c(1, 2),
             col.var="contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07", "#008f39"),
             repel = TRUE # Avoid text overlapping
             )

#library(rgl)
#Visualizar la contribución de los score
#plot3d(pca1$var$contrib[,1:3])
fviz_pca_biplot(pca1,
                col.var = "#00AFBB",
                col.ind = "#E7B800")

#Contribucaciones del Componente Principal 1,  por Variable
fviz_contrib(pca1, choice = "var", axes = 1, top = 10)

  • Las variables total_spent, income y num_catalog_purch, influyen significativamente en el PC1.
#Contribuciones por dimensión a los CP
evi.val <- get_eigenvalue(acp)
  • Hasta 7 dimensiones se tienen por PCA con un eigenvalor mayor a uno los cuales tienen un nivel de explicación del 62.6%, con una varianza de 3.

2.2.2 2. Conclusiones

Dados los resultados podemos concluir que:

  1. Se pueden observar dos grupos a partir de las pruebas realizadas.

2.3 Cluster con K-means

El análisis de clusters funciona como criterios adicionales para la formación de grupos tras el análisis de PCA.

#Análisis de Clusters
library(fpc)
set.seed(123)
#Intentamos con dos grupos, dada la conclusión obtenida en PCA.
k1 <- kmeansruns(dat, #datos centrados
                 k =2,#2 CP
                 runs = 100)#Al menos 100 iteraciones

#Visualizamos
fviz_cluster(k1, data = dat)

También podemos trabajar con otra función

# Calcule k-medias con k = 2
set.seed(123)
km.2 <- kmeans(dat, 2, nstart = 25)
# Imprime los resultados

#Agregar los clusters a la base de datos inicial
aggregate(datos, by=list(cluster=km.2$cluster), mean)
#Visualizamos
fviz_cluster(
   km.2,
   data = dat,
   palette = c("#2E9FDF", "#00AFBB"),
   ellipse.type = "euclid",
   # Concentration ellipse
   star.plot = TRUE,
   # Avoid label overplotting (slow)
   ggtheme = theme_minimal()
)

#Evaluamos el número de cluster óptimo con la base de datos dat, es es la que se trabaja y se encuentra estandarizada.
set.seed(123)

#Método de la suma de cuadrados dentro de los clusters
fviz_nbclust(dat,
             kmeans, #Con KMedias
             method = "wss")+#método de separación de medias al cuadrado
  labs(subtitle = "Elbow method")#Gráfico de codo

#Método de interpretación y validación de la coherencia dentro del análisis de grupos
fviz_nbclust(dat,
             kmeans, #Con KMedias
             method = "silhouette")+#meétodo de separación de medias al cuadrado
  labs(subtitle = "Silhouette method")

#Método para estimar el número óptimo de clústeres comparando la disperción de datos originales con muestras creadas
#Gab Statistic
fviz_nbclust(dat, kmeans, nstart = 25, method = "gap_stat", nboot = 50)+
  labs(subtitle = "Gap Statistic Method")

#Metodo Completo
#library(NbClust)
#res.nbclust <- NbClust(dat,distance = "euclidean",min.nc = 2,max.nc = 9, method = "complete",index="all")
  • El método de codo y de silueta arrojan un K de 2
  • El método gap con toda la realidad bien alterada arroja uno de 10
#Dado que son aleatorios, veamos que son estables
cf1 <- clusterboot(dat, 
                   #numero de remuestreos
                   B=100,
                   #metodo de remuestreo
                   bootmethod = c("jitter", "boot"),
                   #Metodo de clustering
                   clustermethod = kmeansCBI,
                   #Numero de clusters a evaluar
                   krange=2, seed = 123)

#cf1
cf1
## * Cluster stability assessment *
## Cluster method:  kmeans 
## Full clustering results are given as parameter result
## of the clusterboot object, which also provides further statistics
## of the resampling results.
## Number of resampling runs:  100 
## 
## Number of clusters found in data:  2 
## 
##  Clusterwise Jaccard jittering mean:
## [1] 0.9878500 0.9920876
## dissolved:
## [1] 2 0
## recovered:
## [1] 98 98
##  Clusterwise Jaccard bootstrap (omitting multiple points) mean:
## [1] 0.9546201 0.9700675
## dissolved:
## [1] 6 0
## recovered:
## [1] 94 94
  • Con 2 grupos, se arrojan las estadisticas e indices más estables. Por lo que se puede tomar la decisión de agrupación de datos en 2 clusters o dos dimensiones para PCA.
#Con los datos sin transformar evaluamos los dos clusters o PCA
dta %>% 
  ggplot(aes(x=factor(k1$cluster), y=income, fill = factor(k1$cluster)))+
  geom_boxplot()+geom_point()+
  xlab("Clusters")+labs(fill = "Cluster")

2.3.1 Agrupamos Clusters

# Base inicial sin estandarizar
cluster.dta <- datos

# Data frame con segmentación de clusters
# Recuperamos la clasificaicón en nuestro modelo k1 
cluster.dta$cluster <- km.2$cluster

#Tabla de frecuencias de cada cluster correspondiente al PCA
tabyl(cluster.dta$cluster)
#Factor tipo de cluster
cluster.dta$cluster <- factor(cluster.dta$cluster)

head(cluster.dta, n = 4)

2.4 Tidymodels 2.0

Uno de los beneficios del ecosistema tidymodels es la flexibilidad y facilidad de probar diferentes enfoques para el mismo tipo de tarea. Por ejemplo, podemos cambiar PCA por UMAP , un algoritmo completamente diferente para la reducción de dimensionalidad basado en ideas del análisis de datos topológicos. El paquete de inserción proporciona pasos de recetas para crear incrustaciones, incluido UMAP. Cambiemos el paso PCA por el paso UMAP.

2.4.1 1. Dividimos datos de entrenamiento y datos de prueba

# Plantamos semilla para reproducibilidad 
set.seed(123)

# seleccionamos el 80% para el entrenamiento
data_split <- initial_split(cluster.dta, 
                            prop = .8,
                            strata = cluster)

# Create data frames for the two sets:
train_data <- training(data_split)
test_data  <- testing(data_split)

2.4.2 2. UMAP: Receta para UMAP

#Para reproducibilidad
set.seed(123)

#Preparamos datos
umap_rec <- recipe(cluster ~.,data = train_data) %>%
  #Normalizamos datos
  step_normalize(all_predictors(), -all_nominal()) %>%
  #Convertinos en Componentes Principales
  step_umap(all_numeric_predictors(),outcome = vars(cluster))

2.4.3 2.1 UMAP: Gráfico para UMAP

#Gráfico para UMAP
umap_prep <- prep(umap_rec)

#Gráfico UMAP
juice(umap_prep) %>%
  ggplot(aes(UMAP1, UMAP2, label = cluster)) +
  geom_point(aes(color = factor(cluster)), alpha = 0.5, size = 1) +
  scale_color_manual(values = c("#D6B3EB", "#C4DBA6")) +
  labs(color = NULL)

#Si tuviera más niveles podria ayudar:
#geom_text(check_overlap = TRUE, hjust = "inward", family = "IBMPlexSans")+

2.4.4 3. UMAP: Modelo para clasificación

#Modelo
#Creamos el modelo de arbol de decisiones para clasificación
dt_model <- rand_forest(trees = 1000) %>% 
  #Queremos clasificación
  set_mode("classification") %>% 
  #Multinomial
  set_engine("ranger")

2.4.5 4. UMAP: Flujo de trabajo

#Flujo de trabajo
dt_wf <- workflow() %>% 
  add_recipe(umap_rec) %>%
  add_model(dt_model)

dt_wf
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: rand_forest()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 2 Recipe Steps
## 
## • step_normalize()
## • step_umap()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Random Forest Model Specification (classification)
## 
## Main Arguments:
##   trees = 1000
## 
## Computational engine: ranger

2.4.6 5. UMAP: Entrenar Modelo

#Entrenar modelo
dt_fit <- dt_wf %>% 
  fit(data = train_data)

2.4.7 6. UMAP: Estimar

#Predecimos con test
dt_predict <- test_data %>%
   bind_cols(predict = predict(dt_fit, test_data))

2.4.8 7. UMAP: Evaluación

#Evaluamos el Modelos
f_meas(dt_predict, cluster, .pred_class)

2.5 Random Forest: Clasificación.

2.5.1 1. Receta de datos: Tareas de prepración

dt.recipe <- recipe(cluster ~ .,
                    data = train_data)

#Como no aplicamos cambios, no es necesario la función "juice()"

2.5.2 2. Preparamos Modelo

#Creamos el modelo de arbol de decisiones para clasificación
set.model <- rand_forest(trees = 1000) %>% 
  #Queremos clasificación
  set_mode("classification") %>% 
  #Multinomial
  set_engine("ranger")

2.5.3 3. Colapsamos en flujo de trabajo

#Creamos flujo de trabajo para agregar la receta
dt.workfow <- workflow() %>% 
  add_recipe(dt.recipe) %>% 
  add_model(set.model)

2.5.4 4. Entrenamos modelo con datos de entrenamiento

#Entrenamos el modelo
dt.fit <- dt.workfow %>% 
  fit(data = train_data)

2.5.5 5. Estimamos

#Predecimos la base de prueba con el modelo entrenado
dt.predict <- test_data %>%
  bind_cols(predict = predict(dt.fit, test_data))

2.5.6 6. Prueba de precisión del modelo

#Prueba F para el test del modelo
f_meas(dt_predict, cluster, .pred_class)

2.6 Regresión

CODIGO PENDIENTE…

Porque las PC \(x_1, . . . ,x_n\) son ortogonales, se elimina el problema de multicolinealidad y la ecuación de regresión siempre contendrá todas las variables en X (porque cada PC es una combinación lineal de las variables en X producida por un vector propio de \(Z′Z\)).

Debido al uso de PC ortogonales, se cree que la PCR aumenta la precisión numérica de las estimaciones de regresión. (Hadi y Ling, 1998)

Ahora, necesitamos eliminar la popularidad de las variables y ejecutar nuevamente la PCA:

Para realizar la regresión, utilizamos la matriz Z que consta de r o p componentes principales. De esa manera obtenemos coeficientes al hacer una regresión sobre componentes principales: