Tus modelos tienen delulu

Hay gente que ven conejos en la luna, o ven un rostro humano en aquella montaña de Marte que la NASA fotografió en los 70, e incluso quien quiere ve fenómenos paranormales en todos lados.

Tenemos un cerebro tan bueno encontrando patrones de caras que las encuentra incluso donde no las hay. A este fenómeno se le llama pareidolia: el detector de patrones funciona dónde no debería, generando falsas alarmas. Esto, que parece una curiosidad de sobremesa, es uno de los problemas más comunes y graves por desgracia en el campo del aprendizaje no supervisado, y el motivo por el que nunca he sido muy fan de esta disciplina.

Cualquier método que busque “grupos ocultos” o “perfiles de individuos” en unos datos acabará por alucinar. Estos métodos encuentran grupos. Siempre. Le des lo que le des, los modelos siempre devuelven una respuesta. Pero no pueden especificar si esa respuesta es válida o no. Entonces: si un método siempre encuentra grupos, ¿cómo sabemos si los que encuentra son de verdad o son humo?

No tengo para nada en claro la respuesta, este post va más de dejar constancia de mi ignorancia que de mis conocimientos. La única forma honesta que se me ocurre de encontrar límites a este dilema es… haciendo trampa: voy a coger unos datos en los que ya sé la respuesta correcta y comprobar si el método la descubre por su cuenta. Si la encuentra, oye, que bien, tenemos una herramienta más el repertorio de data scientist. Si no, hemos aprendido algo sobre sus límites. Para esa trampa tenemos un clásico perfecto: las tres especies de flores más famosas de toda la historia de la ciencia de datos. Pero antes de ponerme a hacer pruebas quiero sentar las bases de con que probalemas me estoy encontrando.

¿Qué es un “perfil latente”?

Imagina que tienes una base de datos de individuos y, de cada uno, varias medidas numéricas continuas: cuánto mide, cuánto pesa, su puntuación en algún test. Siempre variables numéricas, no categóricas (eso para más adelante)

La sospecha es que esos individuos no son todos iguales: que hay tipos. Por ejemplo, en una base de datos de deportistas tendrías un tipo grande y pesado (posiblemente especializado en fuerza), y otro pequeño y ligero (especializado en velocidad). A cada uno de esos tipos lo vamos a llamar un perfil.

La palabra latente significa simplemente que no se observa directamente. En la base de datos no hay una columna que diga “este individuo es del tipo 1”. Esa etiqueta está escondida; la única pista de que existe son las medidas continuas. El método tiene que deducir los tipos a partir de las medidas, no al revés. Es un problema de análisis no supervisado clásico.

Un modelo de perfiles latentes (en inglés, Latent Profile Analysis, abreviado LPA) hace justo eso:

Supone que cada individuo pertenece a uno de unos pocos perfiles ocultos.
Supone que cada perfil tiene unos valores medios característicos en las medidas.
A partir de los datos, estima cuántos perfiles hay, cuáles son sus medias y con qué probabilidad pertenece cada individuo a cada perfil.

Quédate con la palabra probabilidad. El modelo no dice “este individuo es del tipo 1 y punto”. Dice “es del tipo 1 con un 80% de probabilidad y del tipo 2 con un 20%”. Esa incertidumbre será importante.

Técnicamente, el LPA asume que cada perfil genera sus datos siguiendo una distribución normal con su media y su anchura propias. Los perfiles, además, pueden tener tamaños distintos: uno puede ser la mayoría de la población y otro una minoría. Todos juntos forman una mezcla de normales, es decir, varias campanas de distinto tamaño superpuestas. El trabajo del modelo es separar las campanas que están mezcladas.

La Figura 1 dibuja exactamente esa idea. A la izquierda, lo que tú mides: una nube de individuos que, vista de golpe, parece una sola población. A la derecha, lo que en realidad hay debajo: esa misma nube son tres perfiles distintos, de distinto tamaño y forma, apilados.

Figura 1: La misma población vista de dos formas. A la izquierda, lo que mides: una única distribución. A la derecha, los tres perfiles de los datos que la componen, cada uno con su tamaño (peso), su centro y su anchura. Las tres campanas ponderadas suman exactamente la curva de la izquierda.

Fíjate en el truco visual de la Figura 1: la curva gris de la izquierda no es una campana de verdad, es la suma de las tres de la derecha. Cada perfil contribuye según su peso, por eso ninguno sobresale por encima de la población: encajan debajo de ella. El reto del LPA es exactamente este: solo le damos lo de la izquierda y tiene que reconstruir lo de la derecha. Y la dificultad que va a marcar todo el artículo ya se intuye: cuanto más se pisan las campanas, más difícil es separarlas.

El reto: tres especies que el modelo no conoce

Para poner el LPA a prueba con la respuesta correcta en la mano usamos iris, una base de datos incluida en R: 150 flores de tres especies: setosa, versicolor y virginica; hay 50 de cada tipo, con cuatro medidas continuas (largo y ancho del pétalo y del sépalo).

Data summary
Name	iris
Number of rows	150
Number of columns	5
_______________________
Column type frequency:
factor	1
numeric	4
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
Species	0	1	FALSE	3	set: 50, ver: 50, vir: 50

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Sepal.Length	1	5.84	0.83	4.3	5.1	5.80	6.4	7.9	▆▇▇▅▂
Sepal.Width	1	3.06	0.44	2.0	2.8	3.00	3.3	4.4	▁▆▇▂▁
Petal.Length	1	3.76	1.77	1.0	1.6	4.35	5.1	6.9	▇▁▆▇▂
Petal.Width	1	1.20	0.76	0.1	0.3	1.30	1.8	2.5	▇▁▇▅▃

Nosotros sí sabemos la especie de cada flor, pero al modelo se la esconderemos.

Antes de modelar nada ¿Cómo se sospecha que hay perfiles?

Aquí está la primera lección de mundo real, no de juguete. Nadie mira un solo gráfico y decide que hay grupos. El trabajo previo es buscar pistas.

Un mapa de correlaciones no sirve para esto. La correlación mide qué variables se mueven juntas, no qué individuos se agrupan; es ciega a las relaciones no lineales y, peor aún, la propia mezcla de grupos fabrica correlaciones que no existen dentro de cada grupo. No es una herramienta de detección de perfiles.

La primera herramienta que si nos vale la pena: mirar la distribución de cada variable por separado. La pista que delata grupos ocultos es la multimodalidad, es decir, que una variable no tenga una sola joroba sino dos o más. Recuerda la Figura 1: si dos perfiles están lo bastante separados, su suma deja de parecer una campana y le salen bultos. Esta vista se hace con un único gráfico faceteado, y funciona igual con 4 variables que con 15:

Código

iris |>
  pivot_longer(-Species, names_to = "variable", values_to = "valor") |>
  ggplot(aes(valor)) +
  geom_density(fill = col_general, color = col_general, alpha = 0.4) +
  facet_wrap(~ variable, scales = "free") +
  labs(x = NULL, y = "Densidad") +
  theme_minimal(base_size = 13)

Figura 2: Distribución de cada variable por separado. Las medidas del pétalo muestran dos jorobas claras (pista de grupos ocultos); las del sépalo, una sola. Esta misma receta sirve para escanear 15 variables de golpe.

En la Figura 2, las dos medidas del pétalo tienen claramente dos jorobas: eso grita “aquí hay al menos dos grupos”. Las del sépalo parecen una sola campana: solas, no te harían sospechar nada. Moraleja: una sola variable puede engañarte en cualquier dirección, por eso se miran todas.

La segunda herramienta cruza las variables de dos en dos, para ver si la nube de puntos se apelmaza en grumos. Para un puñado de variables, la matriz de dispersión lo enseña todo a la vez (en la diagonal, cada variable a solas; fuera de ella, cada par):

Código

library(GGally)

ggpairs(
  iris, columns = 1:4,
  lower = list(continuous = wrap("points", alpha = 0.5, color = "#243A6E", size = 0.8)),
  diag  = list(continuous = wrap("densityDiag", color = "#243A6E", fill = "#243A6E", alpha = 0.3)),
  upper = list(continuous = wrap("cor", color = "gray20"))
) +
  theme_minimal(base_size = 10)

Figura 3: Matriz de dispersión de las cuatro medidas, sin etiquetas de especie. En varios paneles se ve un grumo claramente despegado del resto: la sospecha de que hay grupos nace aquí, antes de ajustar ningún modelo.

En la Figura 3, varios paneles muestran un grumo despegado del resto de la nube. Sin saber todavía qué es, ya tienes una sospecha fundada: hay estructura, probablemente al menos dos grupos.

¿Y si tuvieras quince variables? La matriz de la Figura 3 se vuelve ilegible enseguida, pero el escaneo de densidades de la Figura 2 escala sin problema y sigue siendo el caballo de batalla: la multimodalidad delata grupos haya las variables que haya.

Aplastar muchas dimensiones en dos: UMAP y t-SNE

Las dos herramientas anteriores enseñan los datos a trozos: una variable cada vez, o un par cada vez. Pero a nuestro ojo se le da de maravilla detectar grumos en una nube de puntos, y eso solo sabemos hacerlo en dos dimensiones.

Eso es exactamente una proyección a dos dimensiones: aplastar el espacio de muchas variables sobre un plano intentando respetar las distancias. Piensa en la sombra de un objeto tridimensional sobre una pared. Pierdes información, pero si la luz incide bien, la silueta sigue siendo reconocible. Toda la cuestión está en desde qué ángulo iluminar.

La forma más simple, el análisis de componentes principales (PCA), ilumina con luz recta: busca el ángulo que conserva la mayor variación posible y proyecta. Es rápido y fiel a la estructura global, pero al ser lineal no sabe “desenrollar” grupos separados por fronteras curvas; si dos perfiles se distinguen de forma no lineal, la sombra de PCA los puede dejar pegados.

UMAP y t-SNE son las versiones no lineales del mismo truco. Miran quién es vecino de quién en el espacio original y construyen un mapa en 2D que respete esas vecindades. Por eso sacan a la luz grupos con formas retorcidas que PCA aplastaría. Ambas tienen un mando que regula cuán local es la mirada: en t-SNE se llama perplexity; en UMAP, n_neighbors; y conviene tenerlo presente, porque ese mando cambia el dibujo.

Probémoslo con iris. Estandarizamos las cuatro medidas (para que ninguna domine solo por su unidad, algo que repetiremos al modelar) y pedimos las dos proyecciones:

Código

library(uwot)
library(Rtsne)

# Debemos pasar los datos a formato matriz y normalizarlos para que estas técnicas funcionen
X_mat <- iris |> 
  dplyr::select(-Species) |>
   scale() |> 
   as.matrix()

set.seed(20242026)
emb_umap <- umap(X_mat, n_neighbors = 15, min_dist = 0.1)
emb_tsne <- Rtsne(X_mat, perplexity = 30, check_duplicates = FALSE)$Y

proyecciones <- bind_rows(
  tibble(metodo = "UMAP",  d1 = emb_umap[, 1], d2 = emb_umap[, 2], Especie = iris$Species),
  tibble(metodo = "t-SNE", d1 = emb_tsne[, 1], d2 = emb_tsne[, 2], Especie = iris$Species)
)


proyecciones_2x2 <- bind_rows(
  proyecciones |> mutate(vista = "Lo que verías", grupo = "Sin etiqueta"),
  proyecciones |> mutate(vista = "La solución real",  grupo = as.character(Especie))
) |>
  mutate(vista = factor(vista, levels = c("Lo que verías",
                                          "La solución real")))

colores_2x2 <- c("Sin etiqueta" = col_general,
                 setNames(col_grupos, levels(iris$Species)))

ggplot(proyecciones_2x2, aes(d1, d2, color = grupo)) +
  geom_point(size = 2, alpha = 0.45) +
  facet_grid(vista ~ metodo, scales = "free") +
  scale_color_manual(values = colores_2x2,
                     breaks = levels(iris$Species),
                     name = "Especie real") +
  labs(x = "Dimensión 1", y = "Dimensión 2") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = "bottom",
    panel.border  = element_rect(color = "gray70", fill = NA, linewidth = 0.5),
    panel.spacing = unit(1, "lines"),
    strip.background = element_rect(fill = "gray90", color = "gray70", linewidth = 0.5)
  )

Figura 4: Mapa con al proyección en 2 variables de Las cuatro medidas de cada flor, tanto con UMAP y con t-SNE. El color (especie real) no estaría disponible en una exploración de verdad; lo añadimos para comprobar. Ambas técnicas cuentan lo mismo: setosa forma un grumo aislado, y versicolor y virginica comparten otro con una frontera difusa.

En la Figura 4 las dos técnicas cuentan la misma historia: un grumo perfectamente aislado y otro, mayor, más enmarañado. Coloreando por la especie real (ese lujo que en un problema de verdad no tendrías, y que aquí usamos solo para comprobar) se ve que el núcleo de puntos solitario es setosa y que el bulto grande puntos son versicolor y virginica rozándose, sin una frontera limpia. Otra vez lo mismo: dos grupos clarísimos y un tercero que pelea por despegarse del segundo. La proyección no te da el número de perfiles, pero te deja la sospecha afinada: dos seguros, quizá tres.

¿Y si el solape se deshace en tres dimensiones?

Para cerrar esta sección me queda una duda legítima: ¿y si en realidad las dos especies se separasen en una tercera que el plano ha pisado? La clásica de los físicos: “mi sistema funciona en 27 dimensiones”. Podemos hacer la proyección en 3D, UMAP por ejemplo entrega un mapa en 3D con solo pedirle una dimensión más (n_components = 3), y con plotly lo volvemos interactivo: gíralo con el ratón y pasa el cursor por encima de cualquier flor para ver sus medidas.

Código

library(plotly)

set.seed(2024)
emb_3d <- umap(X_mat, n_neighbors = 15, min_dist = 0.1, n_components = 3)

datos_3d <- iris |>
  mutate(
    U1 = emb_3d[, 1], U2 = emb_3d[, 2], U3 = emb_3d[, 3],
    etiqueta = paste0(
      "Especie: ", Species,
      "<br>Sépalo: ", Sepal.Length, " x ", Sepal.Width, " cm",
      "<br>Pétalo: ", Petal.Length, " x ", Petal.Width, " cm"
    )
  )

plot_ly(
  datos_3d, x = ~U1, y = ~U2, z = ~U3,
  color = ~Species, colors = col_grupos,
  text = ~etiqueta, hoverinfo = "text",
  type = "scatter3d", mode = "markers",
  marker = list(size = 4, opacity = 0.85)
) |>
  layout(
    scene = list(
      xaxis = list(title = "UMAP 1"),
      yaxis = list(title = "UMAP 2"),
      zaxis = list(title = "UMAP 3")
    ),
    legend = list(title = list(text = "Especie real"))
  )

Figura 5: UMAP en tres dimensiones. Gira la nube con el ratón y pasa el cursor sobre cualquier punto para ver su especie y sus cuatro medidas. setosa se mantiene aislada desde cualquier ángulo; versicolor y virginica siguen interpenetradas también en la tercera dimensión.

Puedes mover la Figura 5 a tu gusto con ambos clicks de ratón para hacerte la idea: setosa sigue siendo una isla desde cualquier ángulo, pero versicolor y virginica continúan pegadas también en la tercera dimensión. El solapamiento sigue pasando a pesar de la proyección de a una dimensión adicional. Ni el mapa 2D ni el 3D son la verdad, son resúmenes deformados de ella, en los que mires como mires, estas dos especies son bastante cercanas si no tenemos más variables por las que separarlas. La proyección 3D aporta una confirmación visual más, valiosísima cuando tienes quince variables y ningún diagrama capaz de mostrarlas todas a la vez, porque girar la nube es lo más cerca que vas a estar de “ver” el espacio entero.

No quiero extender más esta sección con un caso de”y si hiciermos…” iterativo, en algún punto hay que acabar. Entonces, me toca repetir el aviso: UMAP y t-SNE siempre dibujan grumos bonitos: están diseñadas para agrupar, así que agrupan. Ni el tamaño de cada mancha ni la distancia entre manchas tienen que significar lo que parecen —dos grumos pegados pueden estar lejísimos en el espacio real. Por eso una proyección es un generador de sospechas, nunca una prueba.

La conclusión vendrá de otras cosas la dará el modelo y sobretodo al final, el expertisement que se tenga del problema que se está tratando.

AL menos sabemos que hay grupos, sabemos que hay algo, ya tiene sentido llamar al LPA.

Calibrando un modelo de LPA

Buscando el mejor modelo

Vamos a usar tidyLPA, una interfaz que ordena todo el flujo del LPA en pasos limpios. Por dentro se apoya en el motor de mclust, pero nos deja una sintaxis cómoda y comparable.

Antes de ajustar, estandarizamos las cuatro medidas (media 0, anchura 1) para que ninguna pese más que otra solo por su unidad.

Código

library(tidyLPA)
library(knitr)
library(DT)

indicadores <- iris |> 
  dplyr::select(-Species) |> 
  scale() |> 
  as.data.frame()

datatable(round(indicadores, 2), options = list(pageLength = 5))

Entonces toca tomar una decisión que mucha gente pasa por alto: qué forma le permitimos tener a cada campana. Como dijimos antes, un LPA es una mezcla de gaussianas: cada perfil es una distribución normal. Pero… ¿podrían las campanas no ser normales, ser asimétricas o tener colas pesadas? En un LPA hasta el nivel que controlo, no.

Que cada componente sea normal es precisamente lo que define al método; existen otras familias de mixturas: de distribuciones t para colas pesadas o normales sesgadas para asimetría; pero aquí la campana es siempre normal.

Cabe recordar que “normal” no es una sola forma, es toda una familia. En una dimensión, una normal queda descrita por dos números: su media (dónde está el centro) y su varianza (cómo de ancha es). En varias dimensiones la media pasa a ser un punto en el espacio de datos, y la varianza se convierte en algo más rico, la matriz de covarianzas, este parámetro de modelo decide cuánta de esa libertad concedemos a cada distribución, y lo hace respondiendo a tres preguntas. - ¿Tienen todas las variables la misma anchura dentro de un perfil, o puede cada una tener la suya? - ¿Permitimos que las variables se correlacionen (que la campana se incline), o la obligamos a quedarse alineada con los ejes? - ¿Comparten todos los perfiles la misma forma, o puede cada perfil tener la suya propia?

En un extremo, el modelo más rígido responde “todo igual y sin inclinar”;Una forma rígida es austera y estable, pero, como vamos a ver enseguida, cuando le prohíbes a la campana inclinarse, el modelo se las apaña para capturar esa correlación que no puede modelar… inventándose perfiles de más. En el otro, el modelo más flexible responde “cada cual a su aire”: cada perfil es su propio elipsoide, con su tamaño, su anchura por variable y su inclinación. Una forma más flexible se ajusta mejor a grupos reales que son elípticos o están inclinados, pero cuesta muchos más parámetros (la matriz de covarianzas crece deprisa con el número de variables) y pide más datos para estimarse bien.

tidyLPA lo controla la forma de cada campana con argumento que es el “número de modelo”. para este caso, me interesan dos extremos:

Modelo 1 (rígido): cada perfil es una campana “redonda”, y se prohíbe que las variables estén correlacionadas dentro de un perfil.
Modelo 6 (flexible): cada perfil puede tener su propia forma y permitir que sus variables se correlacionen libremente.

Ajustamos de 2 a 4 perfiles con ambos modelos, que eran el número de opciones plausibles que hemos detectado previamente, y miramos todos los indicadores de ajuste:

Código

modelos_LPA <- indicadores |>
  estimate_profiles(2:4, models = c(1, 6))

datatable(
    get_fit(modelos_LPA) |> mutate(across(everything(),~round(.x,2))) , 
    options = list(pageLength = 10)
  )

El diagnóstico, criterio por criterio

Elegir modelo no es mirar un número y sentenciar. Existe toda una familia de criterios, y casi nunca coinciden; el oficio está en entender qué pregunta responde cada uno y por qué discrepan. La tabla de ajuste mezcla además dos decisiones distintas, así que las separamos: primero qué forma de campana, y después cuántos perfiles.

Empecemos por verlo todo junto. Giramos la tabla para tener cada modelo en una columna y cada métrica en una fila, que con pocas etiquetas y muchas métricas se lee mejor así:

Código

ajuste <- get_fit(modelos_LPA)

ajuste |>
  filter(Classes >= 2) |>
  transmute(
    modelo = paste0("Modelo ", Model, " · ", Classes, "p"),
    `Log-verosimilitud` = LogLik,
    `Parámetros`        = parameters,
    AIC, BIC, SABIC, CAIC, AWE,
    Entropía            = Entropy,
    `Prob. mínima`      = prob_min,
    `N grupo menor`     = n_min * n,
    `BLRT (p)`          = BLRT_p
  ) |>
  pivot_longer(-modelo, names_to = "Métrica", values_to = "valor") |>
  pivot_wider(names_from = modelo, values_from = valor) |>
  mutate(across(where(is.numeric), \(x) round(x, 2))) |>
  kable(caption = "Todos los modelos ajustados, con las métricas en filas. En los criterios de información más bajo es mejor; en entropía y probabilidad mínima, más alto.")

Tabla 1: Todos los modelos ajustados, con las métricas en filas. En los criterios de información más bajo es mejor; en entropía y probabilidad mínima, más alto.

Métrica	Modelo 1 · 2p	Modelo 1 · 3p	Modelo 1 · 4p	Modelo 6 · 2p	Modelo 6 · 3p	Modelo 6 · 4p
Log-verosimilitud	-597.25	-469.77	-464.43	-322.69	-288.53	-275.83
Parámetros	13.00	18.00	23.00	41.00	62.00	83.00
AIC	1220.51	975.54	974.87	727.39	701.05	717.65
BIC	1259.65	1029.73	1044.11	850.82	887.71	967.53
SABIC	1218.50	972.77	971.32	721.07	691.49	704.85
CAIC	1272.65	1047.73	1067.11	891.82	949.71	1050.53
AWE	1361.80	1172.01	1226.59	1177.26	1382.43	1630.50
Entropía	0.99	0.96	0.88	1.00	0.97	0.96
Prob. mínima	0.99	0.97	0.66	1.00	0.98	0.96
N grupo menor	50.00	45.00	14.00	50.00	45.00	12.00
BLRT (p)	0.01	0.01	0.20	0.01	0.02	0.67

Primera decisión: ¿qué forma de campana?

Para empezar debemos elegir un modelo de los 6 que hemos calibrado. La métrica más habitual por la que empezar a mirar esto es el BIC (Bayesian information criteria): resume en un solo número lo bien que un modelo ajusta los datos, poniéndole un coste por cada parámetro que usa. Se suele empezar el análisis por aquí porque premiaa un buen ajuste sin recompensar poner complejidad por las buenas. viendo su fórmula se entiende el criterio:

\[\text{BIC} = \underbrace{-2\ln(\hat{L})}_{\text{lo mal que ajusta}} \;+\; \underbrace{k \cdot \ln(n)}_{\text{peaje por complejidad}}\]

Los dos sumandos son costes, así que más bajo es mejor. El primero baja cuando el modelo ajusta mejor; el segundo sube $\ln(n)$ por cada parámetro $k$, ajuste o no ajuste. El BIC es una factura y queremos pagar poco.

Entonces el modelo flexible (6) gana al rígido (1) a cualquier número de perfiles en BIC. Con dos perfiles, 851 frente a 1260. Y eso tiene miga, porque el flexible es mucho más complejo: 41 parámetros frente a 13. Paga un peaje bastante grande, unos $28 \times \ln(150) \approx 140$ puntos más, y aun así el ajuste mejora 549 puntos, muchísimo más de lo que cuesta esa complejidad.

¿Por qué los modelos flexibles ajustan mejor? Por que deja que las cuatro medidas se correlacionen dentro de cada grupo, cosa que el ajuste rígido lo prohíbe. Significa que esas correlaciones existen y son fuertes. Y eso es justo lo que buscábamos: unas variables que se mueven juntas y forman grupos con forma y orientación propias son la huella de que hay perfiles de verdad. El rígido, al obligar a cada perfil a ser una campana redonda, se niega a verlos como son.

Decisión cerrada: nos quedamos con el flexible. Todo lo que sigue ocurre dentro de él.

nota: el flexible no gana solo en BIC; gana en los cinco criterios de información, así que la decisión es unánime y no depende de cuál miremos. Pa ra este análisis los miré todos, pero pensé que enseñar los cinco para que digan lo mismo cinco veces era paja, de modo que me centré en uno.

Segunda decisión: ¿cuántos perfiles?

LLegamos al quid de la cuestión. Aquí no vale hacer trampa diciendo “pongo 3 porque sé que son 3”. En un problema real no se conoces el número, así que hay que tomar una decisión final entre criterios que muy posiblemente que se contradicen.

Un criterio de información resume el modelo en un solo número que combina dos cosas: cómo de bien ajusta y cuántos parámetros gasta, penalizando lo segundo para no premiar la complejidad porque sí. Más bajo, mejor. Lo que distingue a unos de otros es la dureza de esa penalización. De más blando (tolera más perfiles) a más duro (los recorta):

AIC — Akaike Information Criterion (criterio de información de Akaike, por Hirotugu Akaike). El más permisivo: su peaje es fijo y pequeño (2 por parámetro) y no crece con el tamaño de la muestra. Tiende a quedarse con soluciones más grandes, así que cuidado, porque en muestras decentes suele sobreestimar el número de perfiles.
SABIC — Sample-size Adjusted BIC (BIC ajustado por el tamaño muestral; también ssBIC o aBIC). El BIC, pero con el peaje suavizado mediante una corrección por tamaño muestral. Más indulgente que el BIC, más prudente que el AIC. Suele ser la referencia recomendada cuando esos dos discrepan.
BIC — Bayesian Information Criterion (criterio de información bayesiano; también Schwarz Information Criterion o SIC, por Gideon Schwarz). El de referencia, el que usamos en el diagnóstico del tipo de modelo. Su peaje por parámetro crece con la muestra: cuantos más datos, más severo con la complejidad. Por construcción, es conservador.
CAIC — Consistent AIC (AIC consistente, de Bozdogan). Lleva “AIC” en el nombre, pero es en realidad el BIC con un puntito más de penalización: siempre vota lo mismo que el BIC, o algo más austero.
AWE — Approximate Weight of Evidence (peso aproximado de la evidencia, de Banfield y Raftery). El más duro de todos. Penaliza no solo la complejidad sino también la mala clasificación (lo difusos que quedan los grupos), así que castiga doblemente las soluciones con muchos perfiles solapados. Es el que más recorta y el más reacio a añadir grupos.

Código

ic_largo <- ajuste |>
  filter(Model == 6) |>
  dplyr::select(Classes, AIC, BIC, SABIC, CAIC, AWE) |>
  pivot_longer(-Classes, names_to = "criterio", values_to = "valor")

minimos <- ic_largo |> group_by(criterio) |> slice_min(valor, n = 1)

ggplot(ic_largo, aes(Classes, valor)) +
  geom_line(color = "gray60", linewidth = 0.8) +
  geom_point(color = "gray60", size = 2) +
  geom_point(data = minimos, color = "#2DD4CF", size = 4) +
  facet_wrap(~ criterio, scales = "free_y") +
  scale_x_continuous(breaks = 2:4) +
  labs(x = "Número de perfiles", y = "Valor (más bajo = mejor)") +
  theme_minimal(base_size = 12) +
  theme(panel.border     = element_rect(color = "gray70", fill = NA, linewidth = 0.5),
        panel.spacing    = unit(0.8, "lines"),
        strip.background = element_rect(fill = "gray90", color = "gray70", linewidth = 0.5))

Figura 6: Cinco criterios de información para el modelo flexible. El punto turquesa marca el mínimo (la elección de ese criterio).

Los criterios no coinciden, y eso es el hallazgo.

Los de penalización dura (BIC, CAIC y AWE) tocan fondo en 2 perfiles.
Los más suaves (AIC y SABIC) lo hacen en 3.

Todo el cisma está en el salto de 2 a 3 perfiles: cuesta 21 parámetros y mejora el ajuste en unos 68 puntos.

¿De quién me fío? El $\ln(n)$ del BIC sale de una aproximación pensada para muestras grandes, y con 150 flores y un modelo de 62 parámetros estamos lejos de ese ideal: en ese régimen el BIC penaliza de más y tiende a quedarse corto de perfiles. El SABIC nació precisamente para corregir un problema como este, ajustando el peaje al tamaño real de la muestra, y por eso recupera mejor el número de grupos en muestras pequeñas o medianas, que es justo nuestro caso.

Conclusión: 3 perfiles Toma algo de ventaja como respuesta definitiva.

El BLRT (Bootstrap Likelihood Ratio Test, test de razón de verosimilitud por bootstrap) es una prueba estadística formal a una pregunta muy concreta: ¿el perfil que acabo de añadir mejora el ajuste más de lo que mejoraría por puro azar?

Para entenderlo hay que ver el problema que resuelve. Cuando pasas de 2 a 3 perfiles, el ajuste siempre mejora un poco, aunque el tercer perfil sea basura (overfitting). La pregunta no es si mejora, sino cuánto, y sobre todo: ¿esa mejora es de las que solo da una estructura real, o es del tamaño de la que saldría por casualidad partiendo en dos un grupo que era uno solo?

El BLRT responde a ese “cuánto” mediante bootstrap. La mecánica, en cuatro pasos:

Mide la mejora de ajuste real al pasar de $k$ a $k+1$ perfiles. La medida es el estadístico de razón de verosimilitudes, $LR = -2\ln(\hat{L}_k / \hat{L}_{k+1})$, que crece cuanto más mejora el ajuste al añadir el perfil.
Adopta la hipótesis nula; $H_0$: el modelo verdadero es el de menos perfiles, el de $k$ (frente a la alternativa de que haga falta uno más, $k+1$). Asumiéndola cierta, genera su distribución correspondiente simulando muchos conjuntos de datos a partir del modelo de $k$ perfiles ya estimado.
A cada conjunto simulado le ajusta los modelos de k y k+1 perfiles y recalcula su $LR$. Esto reconstruye la distribución de $LR$ bajo $H_0$: el catálogo de mejoras que produce el azar cuando no hay un perfil más que encontrar.
El valor p es la proporción de réplicas cuyo $LR$ iguala o supera al observado, es decir, la probabilidad de ver una mejora como la nuestra si $H_0$ fuera cierta.

Así que si, la regla de del estadístico es la clásica: p pequeño, la mejora es demasiado grande para ser casualidad, el perfil extra es real; p grande, la mejora cabe holgadamente dentro de lo que da el azar, el perfil no aporta.

Código

ajuste |>
  filter(Model == 6, !is.na(BLRT_p)) |>
  ggplot(aes(Classes, BLRT_p)) +
  geom_hline(yintercept = 0.05, linetype = "dashed", color = "gray50") +
  geom_line(color = "gray60", linewidth = 0.8) +
  geom_point(aes(color = BLRT_p < 0.05), size = 4) +
  scale_color_manual(values = c("TRUE" = "#2DD4CF", "FALSE" = "#8C2C7C"),
                     labels = c("TRUE" = "Aporta (p < 0,05)", "FALSE" = "No aporta"),
                     name = NULL) +
  scale_x_continuous(breaks = 2:4) +
  labs(x = "Número de perfiles", y = "Valor p del BLRT (k frente a k−1)") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Figura 7: Valor p del BLRT al añadir cada perfil secuencialmente. Por debajo de la línea, la mejora supera al azar y el perfil aporta; por encima, no.

El BLRT no penaliza a ojo: contrasta la mejora real contra su propia distribución por azar. No es casualidad que, en los estudios de simulación que enfrentan a todos estos indicadores, sea el que más veces acierta el número real de grupos, por delante del BIC. Cuando ambos discrepan, hay buenas razones para hacerle caso. Mirando la Figura 7 hay una premisa interesante, Pasar de 1 a 2 perfiles supera al azar (p ≈ 0,01) y pasar de 2 a 3 también (p ≈ 0,01), pero pasar de 3 a 4 no (p = 0,65): esa mejora es esperable por casualidad, así que el cuarto perfil no añade estructura, solo ruido.

El BLRT señala que lo ideal son 3 perfiles

El LPA no asigna cada flor a un perfil directamente sino que como todo modelo probabilístico le da una probabilidad de pertenencia a cada uno (70 % al perfil 1, 20 % al 2, 10%al 3 por ejemplo). Cuando esas probabilidades son nítidas (casi 100 % a un perfil y casi 0 % al resto), los grupos están bien separados; cuando quedan repartidas (50 % y 50 %), el modelo está dudando. Dos métricas resumen esa nitidez desde ángulos distintos:

La entropía mira el conjunto: resume, en un solo número de 0 a 1, cómo de seguras son todas las asignaciones a la vez. Vale 1 cuando cada flor pertenece a su perfil sin ambigüedad y cae hacia 0 cuando las pertenencias están, en promedio, repartidas. Es la separación global de la solución. El umbral habitual es ≥ 0,80 para considerar que la clasificación es fiable. Por encima de 0,80 se habla de separación buena; por encima de 0,90, excelente
La probabilidad mínima mira el eslabón más débil: de entre todos los perfiles, toma el que peor se clasifica y mira con qué probabilidad media se asigna correctamente a sus miembros. Una entropía alta puede esconder un perfil concreto mal definido; esta métrica lo destapa, porque se fija justo en el peor. La referencia es ≥ 0,70 por clase; es decir, ningún perfil debería clasificar a sus miembros con una certeza media inferior a 0,70.

Código

ajuste |>
  filter(Model == 6, !is.na(Entropy)) |>
  dplyr::select(Classes, Entropía = Entropy, `Prob. mínima` = prob_min) |>
  pivot_longer(-Classes, names_to = "métrica", values_to = "valor") |>
  ggplot(aes(Classes, valor, color = métrica)) +
  geom_line(linewidth = 1) +
  geom_point(size = 3) +
  scale_color_manual(values = c("Entropía" = "#6F4197", "Prob. mínima" = "#2DD4CF"),
                     name = NULL) +
  scale_x_continuous(breaks = 2:4) +
  coord_cartesian(ylim = c(0.85, 1)) +
  labs(x = "Número de perfiles", y = "Valor (más alto = mejor)") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Figura 8: Separación entre perfiles (entropía) y seguridad de la clasificación (probabilidad mínima) para el modelo flexible.

Ambas métricas coinciden: en 3 perfiles la separación global sigue siendo casi perfecta (entropía 0,97) y, mirando el perfil peor definido, la clasificación sigue siendo segurísima (probabilidad mínima 0,98). El tercer perfil es un grupo nítido, y además ninguno de los tres queda mal delimitado. Solo al llegar a 4 empiezan ambas a ceder, señal de que ahí ya se está forzando una división que los datos no sostienen.

Conclusión: los criterios de “separación y confianza”” indica que 3 perfiles son lo ideal.

La última comprobación yo la he llamado “de cordura”: ¿qué tamaño tiene el grupo más pequeño de cada solución? Un perfil que tenga cuatro gatos suele ser un artefacto, no un tipo real.

No he sido capaz de encontrar reglas universales con una base teórica fundamentada (desconfía de quien te venda un número mágico), pero sí pistas bastante coherentes.

En términos relativos, el umbral más citado es que cada clase reúna al menos un 5 % de la muestra (algunos lo suben al 10 % para ser más exigentes).
En términos absolutos, la referencia habitual es un mínimo de 30 casos por clase. Y este segundo manda sobre el primero, porque el porcentaje engaña con el tamaño del dataset: un 5 % de 5.000 son 250 individuos, de sobra para estimar un grupo con estabilidad; pero un 5 % de 200 son 10, un grupo frágil aunque cumpla el porcentaje. El motivo de fondo no es estético sino de estimación: cada perfil calcula sus propias medias y covarianzas a partir de los casos que le tocan, y con el modelo flexible, que gasta muchos parámetros por grupo, un puñado de casos no da para fijarlos; el grupo resultante es inestable, cambia con la semilla y no replica. Esa inestabilidad es la alucinación.

Con todo, un grupo pequeño no es automáticamente falso: las subpoblaciones raras existen y a veces son lo más interesante. La regla no es “bórralo”, sino “un grupo pequeño necesita más pruebas para creértelo”. El de 12 flores las suspende todas: supera por poco el 5 %, pero se queda por debajo de los 25 casos absolutos, el BLRT ya lo daba por no significativo y aparece justo cuando la predicción empeora. No lo descartamos por un umbral aislado, sino porque tres señales independientes apuntan a lo mismo.

Código

ajuste |>
  filter(Model == 6, !is.na(n_min)) |>
  mutate(n_menor = round(n_min * n)) |>
  ggplot(aes(factor(Classes), n_menor)) +
  geom_col(aes(fill = n_menor < 20), width = 0.6) +
  geom_text(aes(label = n_menor), vjust = -0.4, size = 4) +
  scale_fill_manual(values = c("FALSE" = "#6F4197", "TRUE" = "#8C2C7C"), guide = "none") +
  labs(x = "Número de perfiles", y = "Tamaño del grupo más pequeño (nº de flores)") +
  theme_minimal(base_size = 12)

Figura 9: Tamaño del grupo más pequeño en cada solución. El grupúsculo de 12 flores con 4 perfiles es la firma del sobreajuste.

Con 2 y 3 perfiles, el grupo menor tiene 50 y 45 flores: grupos grandes y respetables. Con 4 se desploma a 12 (un 8 % de la muestra): el grupúsculo anecdótico que delata que el modelo ha empezado a inventarse divisiones. Encaja con el BLRT, que en 4 perfiles ya decía “no aporta”.

Todos los criterios anteriores miran lo mismo: cómo de bien encaja el modelo en los datos que ya tiene. La validación cruzada mide algo distinto y complementario: cómo de bien predice datos que no ha visto. La idea es sencilla. Partimos la muestra en cinco trozos; ajustamos el modelo con cuatro de ellos y medimos cómo de verosímiles le resultan las flores del quinto, que se quedó fuera. Rotamos el trozo apartado, y repetimos varias veces con particiones distintas para promediar el azar. Un modelo que se ha inventado perfiles de más ajustará de maravilla sus propios datos, pero predecirá peor los ajenos: el sobreajuste se paga fuera de casa.

Código

library(rsample)
library(mclust)

set.seed(2024)
folds <- vfold_cv(indicadores, v = 5, repeats = 5)

# Log-verosimilitud media por flor en el trozo apartado
loglik_fuera <- function(split, G) {
  entrena <- as.matrix(analysis(split))
  prueba  <- as.matrix(assessment(split))
  m <- Mclust(entrena, G = G, modelNames = "VVV", verbose = FALSE)
  if (is.null(m)) return(NA_real_)
  mean(dens(modelName = m$modelName, data = prueba,
            parameters = m$parameters, logarithm = TRUE))
}

cv <- tidyr::expand_grid(G = 2:4, split = folds$splits) |>
  mutate(loglik = purrr::map2_dbl(split, G, loglik_fuera))

cv_resumen <- cv |>
  group_by(G) |>
  summarise(media = mean(loglik, na.rm = TRUE),
            ee    = sd(loglik, na.rm = TRUE) / sqrt(sum(!is.na(loglik))),
            .groups = "drop")

ggplot(cv_resumen, aes(G, media)) +
  geom_line(color = "gray60", linewidth = 0.8) +
  geom_pointrange(aes(ymin = media - ee, ymax = media + ee, color = G == 3),
                  size = 0.9) +
  scale_color_manual(values = c("TRUE" = "#2DD4CF", "FALSE" = "gray40"), guide = "none") +
  scale_x_continuous(breaks = 2:4) +
  labs(x = "Número de perfiles", y = "Log-verosimilitud media fuera de muestra") +
  theme_minimal(base_size = 13)

Figura 10: Capacidad de predecir flores no vistas, por número de perfiles (validación cruzada repetida, modelo flexible). Más arriba es mejor. Las barras son el error estándar entre repeticiones.

El resultado de la Figura 10 es sobrio, y conviene leerlo con honestidad. Las soluciones de 2 y 3 perfiles predicen prácticamente igual de bien (sus intervalos se solapan), mientras que la de 4 predice claramente peor. De ahí salen dos lecturas. La primera: pasar a 4 perfiles sí sobreajusta, y la validación cruzada lo castiga, exactamente como avisaba el BLRT. La segunda, más sutil: elegir 3 en lugar de 2 no cuesta nada en capacidad predictiva; el tercer perfil no es un adorno que estropee la generalización, generaliza tan bien como la opción conservadora. Lo que la validación cruzada no hace es romper el empate entre 2 y 3, y no debemos fingir que sí: para distinguirlos hacía falta el BLRT. Las dos pruebas se reparten el trabajo: el BLRT establece que el tercer perfil es una mejora real y no ruido, y la validación cruzada confirma que esa mejora no se paga prediciendo peor.

Resultados

Reuniendo las cinco vistas del diagnóstico (siempre dentro del modelo flexible, porque el rígido quedó descartado por ajustar mucho peor en todos los criterios), los puntos que sostienen la decisión son:

Forma de campana. El modelo flexible bate al rígido en todos los criterios de información y a cualquier número de perfiles: las cuatro medidas se correlacionan dentro de cada grupo, y modelar esa correlación importa.
Criterios de información, divididos. Los de penalización dura (BIC, CAIC, AWE) tocan fondo en 2 perfiles; los más suaves (AIC, SABIC) en 3. SABIC es la métrica más idonea para esta situación, de modo que conservamos en 3.
BLRT, decisivo. El tercer perfil mejora el ajuste de forma significativa (p ≈ 0,01); el cuarto, no (p = 0,65). El indicador más fiable para esta pregunta traza la frontera en 3.
Separación y confianza, intactas en 3. Entropía 0,97 y probabilidad mínima 0,98 para 3 perfiles: grupos nítidos y clasificación segura. A partir de 4 empiezan a degradarse.
Tamaño, la firma del sobreajuste. El grupo más pequeño tiene 45 flores con 3 perfiles, cuando se va a 4 perfiles, el grupo minoritario se desploma a 12: el grupúsculo anecdótico que delata divisiones inventadas.
Predicción. 2 y 3 perfiles generalizan igual de bien; 4 generaliza peor. Elegir 3 no cuesta capacidad predictiva, y 4 la pierde.

El veredicto que sostienen estas vistas es 3 perfiles: la única solución que el BLRT respalda como mejora real, que mantiene separación y tamaños sólidos, y que predice tan bien como la alternativa conservadora. Los criterios que preferían 2 lo hacían solo por penalizar con dureza la complejidad del modelo flexible; ninguno señala que el tercer perfil sea espurio. Y el 4 queda descartado por triple vía: BLRT no significativo, grupo residual minúsculo y peor predicción.

La prueba del algodón: contra la hoja de respuestas

Hemos decidido ajustar un modelos final de LPA flexible de 3 perfiles, solo con lo que sabemos de los datos. Empecemos por ajustar el modelo final:

# Ajustamos el modelo elegido: 3 perfiles, modelo flexible (6)
modelo_final <- indicadores |>
  estimate_profiles(3, models = 6)

# get_data() devuelve, para cada individuo, sus valores y el resultado de la
# clasificación: la columna 'Class' (el perfil asignado) y las 'CPROB_*'
# (probabilidad de pertenecer a cada perfil). Miramos la tabla completa.
datatable(
  get_data(modelo_final) |> mutate(across(where(is.numeric), \(x) round(x, 2))),
  options = list(pageLength = 5),
  caption = "Salida del modelo: valores, perfil asignado (Class) y probabilidades de pertenencia (CPROB)"
)

Como en iris sí conocemos la verdad, podemos permitirnos comprobar si esa decisión, tomada a ciegas, era la acertada. Empecemos por ajustar la matriz de confusión entre el modelo y la realidad:

Código

# La columna 'Class' es el perfil asignado a cada flor. La pegamos a los datos
# originales (get_data mantiene el mismo orden de filas) como factor.
iris_res <- iris |>
  mutate(Perfil = factor(get_data(modelo_final)$Class))

# Matriz de confusión: cruzamos la especie real (que el modelo NO vio) con el
# perfil que ha estimado. La diagonal son los aciertos; fuera, las confusiones.
library(gt)

table(Especie_real = iris$Species, Perfil = iris_res$Perfil) |>
  as.data.frame.matrix() |>
  tibble::rownames_to_column("Especie real \\ Perfil") |>
  gt(rowname_col = "Especie real \\ Perfil") |>
  tab_header(title = "Matriz de confusión: especie real frente a perfil estimado") |>
  tab_spanner(label = "Perfil estimado", columns = everything()) |>
  data_color(
    columns = where(is.numeric),
    fn = scales::col_numeric(c("white", "#2DD4CF"), domain = NULL)
  ) |>
  tab_options(table.align = "center", table.width = pct(60)) |>
  cols_align(align = "center", columns = where(is.numeric))

	Perfil estimado
Matriz de confusión: especie real frente a perfil estimado
	1	2	3
setosa	50	0	0
versicolor	0	45	5
virginica	0	0	50

Los tres perfiles se corresponden casi uno a uno con las tres especies. Para ponerle número a ese “casi” usamos el índice de Rand ajustado (ARI): un valor entre 0 y 1 que mide cómo de bien coincide la agrupación del modelo con la verdadera, sin que importe cómo se numeren los grupos. 1 es coincidencia perfecta; 0 es como tirar los dados. En este caso, el valor ronda 0.9: una recuperación notable.

Ahora, para ver dónde falla, la Figura 11 dibuja cada variable con su perfil estimado; fíjate en cómo los desacuerdos se concentran justo en la frontera entre las dos especies que ya sabíamos solapadas:

Código

library(GGally)
library(dplyr)

# Correspondencia perfil <-> especie (la moda de cada perfil)
correspondencia <- iris_res |>
  dplyr::count(Perfil, Species) |>
  dplyr::group_by(Perfil) |>
  dplyr::slice_max(n, n = 1) |>
  dplyr::ungroup() |>
  dplyr::select(Perfil, especie_esperada = Species)

iris_plot <- iris_res |>
  dplyr::left_join(correspondencia, by = "Perfil") |>
  dplyr::mutate(mal_clasif = Species != especie_esperada)

# Panel de dispersión: puntos por perfil + disonantes con borde negro
lower_fn <- function(data, mapping, ...) {
  ggplot(data, mapping) +
    geom_point(aes(fill = Perfil), shape = 21, color = "gray75",
               size = 1.8, stroke = 0.25, alpha = 0.85) +
    geom_point(data = function(x) dplyr::filter(x, mal_clasif),
               aes(fill = Perfil), shape = 21, color = "black",
               size = 2.8, stroke = 1.1) +
    scale_fill_manual(values = col_grupos)
}

# Diagonal: densidad por perfil
diag_fn <- function(data, mapping, ...) {
  ggplot(data, mapping) +
    geom_density(aes(fill = Perfil), color = NA, alpha = 0.5) +
    scale_fill_manual(values = col_grupos)
}

ggpairs(
  iris_plot,
  columns = c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width"),
  lower  = list(continuous = lower_fn),
  diag   = list(continuous = diag_fn),
  upper  = list(continuous = wrap("cor", color = "gray30")),
  legend = c(3, 1),
  progress = FALSE
) +
  theme_minimal(base_size = 11) +
  theme(
    legend.position  = "bottom",
    panel.border     = element_rect(color = "gray70", fill = NA, linewidth = 0.4),
    panel.spacing    = unit(0.5, "lines"),
    strip.background = element_rect(fill = "gray90", color = "gray70", linewidth = 0.4)
  )

Figura 11: Las cuatro medidas cruzadas dos a dos. Color = perfil estimado por el modelo. Los puntos con borde negro grueso son los mal clasificados. Sin conocer la naturaleza de las variables, esta vista te dice de un vistazo cómo de limpios son los perfiles y dónde se confunden.

Código

library(tidyr)

bind_rows(
  iris_plot |> transmute(grupo = as.character(Species),
                         Sepal.Length, Sepal.Width, Petal.Length, Petal.Width,
                         fuente = "Especie real"),
  iris_plot |> transmute(grupo = as.character(especie_esperada),
                         Sepal.Length, Sepal.Width, Petal.Length, Petal.Width,
                         fuente = "Perfil estimado")
) |>
  pivot_longer(Sepal.Length:Petal.Width, names_to = "Variable", values_to = "valor") |>
  group_by(Grupo = grupo, Variable, fuente) |>
  summarise(resumen = sprintf("%.2f ± %.2f", mean(valor), sd(valor)), .groups = "drop") |>
  pivot_wider(names_from = fuente, values_from = resumen) |>
  arrange(Grupo, Variable) |>
  knitr::kable(caption = "Media ± desviación típica (cm) por variable: especie real frente al perfil estimado que le corresponde. Casi idénticas: el modelo no solo agrupa bien, reconstruye los estadísticos de cada grupo.")

Tabla 2: Media ± desviación típica (cm) por variable: especie real frente al perfil estimado que le corresponde. Casi idénticas: el modelo no solo agrupa bien, reconstruye los estadísticos de cada grupo.

Grupo	Variable	Especie real	Perfil estimado
setosa	Petal.Length	1.46 ± 0.17	1.46 ± 0.17
setosa	Petal.Width	0.25 ± 0.11	0.25 ± 0.11
setosa	Sepal.Length	5.01 ± 0.35	5.01 ± 0.35
setosa	Sepal.Width	3.43 ± 0.38	3.43 ± 0.38
versicolor	Petal.Length	4.26 ± 0.47	4.19 ± 0.44
versicolor	Petal.Width	1.33 ± 0.20	1.29 ± 0.18
versicolor	Sepal.Length	5.94 ± 0.52	5.90 ± 0.53
versicolor	Sepal.Width	2.77 ± 0.31	2.78 ± 0.31
virginica	Petal.Length	5.55 ± 0.55	5.49 ± 0.57
virginica	Petal.Width	2.03 ± 0.27	1.99 ± 0.29
virginica	Sepal.Length	6.59 ± 0.64	6.55 ± 0.62
virginica	Sepal.Width	2.97 ± 0.32	2.95 ± 0.33

¿Hemos aprendido algo?

Hemos respondido al título. Toca hacer sumario de los descubrimientos y aprendizajes.

Sí: el LPA puede redescubrir las tres especies sin que nadie se las diga. Los modelos de la tienen aplicaciones factibles.
No son modelos muy atuomáticos: lo que hace fácil o difícil este problema no es la cantidad de datos, sino cómo de separados estén los perfiles. Aquí, Setosa, estaba bien apartada, se detectaba sola casi. En cambio, versicolor y virginica, estaba solapadas, solo se separaron una vez guiamos al modelo, razonando sobre toda una batería de métricas para la decisión.
Los modelos no son perfectos. La distancia entre perfiles es el factor que de verdad importa.
En un problema real no hay hoja de respuestas. Adjudicamos entre criterios lo mejor que pudimos, pero nadie nos confirmó por detrás que hubiéramos acertado. Con este dataset hemos hecho trampa: teníamos la respuesta correcta, así que pudimos comprobar que elegir 3 era acertado. ¿Hasta dónde podemos fiarnos de esa adjudicación cuando no hay red de seguridad? ¿Qué grupos somos capaces de detectar de verdad, y a partir de qué punto el modelo empieza, otra vez, a ver caras en la tostada? Creo que en otro post voy a expandir esta idea: poner a prueba un LPA cuando no tienes la solución delante, jugando con la separación y con el tamaño de los grupos para medir hasta dónde llega lo que el modelo puede ver.

Tu estadístico de confianza y, cada día, el de menos gente.

Adiós. Os estimo.

--- title: "¿Puede un modelo descubrir 3 especies de flores sin que nadie se lo diga?" description: "Una introducción a los modelos de perfiles latentes con tidyLPA, contada desde cero y puesta a prueba con datos reales." author: "Aitor" date: today lang: es format: html: toc: true toc-depth: 3 toc-location: right code-fold: true code-tools: true smooth-scroll: true theme: darkly fig-width: 8 fig-height: 5 fig-align: center embed-resources: true execute: echo: true warning: false message: false --- <style> .dataTables_wrapper, .dataTables_wrapper label, .dataTables_wrapper .dataTables_info, .dataTables_wrapper .dataTables_paginate { color: #f8f9fa !important; } .dataTables_wrapper .dataTables_filter input, .dataTables_wrapper .dataTables_length select { background-color: #303030 !important; color: white !important; border: 1px solid #6c757d !important; } .dataTables_wrapper .paginate_button { color: white !important; } </style> ## Tus modelos tienen delulu Hay gente que ven [conejos en la luna](https://es.wikipedia.org/wiki/La_cara_de_la_Luna), o ven un [rostro humano](https://es.wikipedia.org/wiki/Cara_de_Marte) en aquella montaña de Marte que la NASA fotografió en los 70, e incluso quien quiere ve [fenómenos paranormales en todos lados](https://www.tiktok.com/@antonikez/video/7291301271703751968). Tenemos un cerebro tan bueno encontrando patrones de caras que las encuentra incluso donde no las hay. A este fenómeno se le llama **pareidolia**: el detector de patrones funciona dónde no debería, generando falsas alarmas. Esto, que parece una curiosidad de sobremesa, es **uno de los problemas más comunes y graves por desgracia en el campo del aprendizaje no supervisado**, y el motivo por el que nunca he sido muy fan de esta disciplina. Cualquier método que busque "grupos ocultos" o "perfiles de individuos" en unos datos acabará por alucinar. [Estos métodos encuentran grupos. Siempre]{.underline}. Le des lo que le des, los modelos siempre devuelven una respuesta. Pero no pueden especificar si esa respuesta es válida o no. Entonces: si un método siempre encuentra grupos, **¿cómo sabemos si los que encuentra son de verdad o son humo?** No tengo para nada en claro la respuesta, este post va más de dejar constancia de mi ignorancia que de mis conocimientos. La única forma honesta que se me ocurre de encontrar límites a este dilema es... haciendo trampa: voy a coger unos datos en los que ya sé la respuesta correcta y comprobar si el método la descubre por su cuenta. Si la encuentra, oye, que bien, tenemos una herramienta más el repertorio de data scientist. Si no, hemos aprendido algo sobre sus límites. Para esa trampa tenemos un clásico perfecto: las tres especies de flores más famosas de toda la historia de la ciencia de datos. Pero antes de ponerme a hacer pruebas quiero sentar las bases de con que probalemas me estoy encontrando. ## ¿Qué es un "perfil latente"? Imagina que tienes una base de datos de individuos y, de cada uno, varias **medidas numéricas continuas**: cuánto mide, cuánto pesa, su puntuación en algún test. Siempre variables numéricas, no categóricas (eso para más adelante) La sospecha es que esos individuos no son todos iguales: que hay tipos. Por ejemplo, en una base de datos de deportistas tendrías un tipo grande y pesado (posiblemente especializado en fuerza), y otro pequeño y ligero (especializado en velocidad). A cada uno de esos tipos lo vamos a llamar un **perfil**. La palabra **latente** significa simplemente *que no se observa directamente*. En la base de datos no hay una columna que diga "este individuo es del tipo 1". Esa etiqueta está escondida; la única pista de que existe son las medidas continuas. El método tiene que **deducir los tipos a partir de las medidas**, no al revés. Es un problema de análisis no supervisado clásico. Un **modelo de perfiles latentes** (en inglés, *Latent Profile Analysis*, abreviado **LPA**) hace justo eso: 1. Supone que cada individuo pertenece a uno de unos pocos perfiles ocultos. 2. Supone que cada perfil tiene unos valores medios característicos en las medidas. 3. A partir de los datos, estima cuántos perfiles hay, cuáles son sus medias y con qué **probabilidad** pertenece cada individuo a cada perfil. Quédate con la palabra **probabilidad**. El modelo no dice "este individuo es del tipo 1 y punto". Dice "es del tipo 1 con un 80% de probabilidad y del tipo 2 con un 20%". Esa incertidumbre será importante. Técnicamente, el LPA asume que cada perfil genera sus datos siguiendo una **[distribución normal](https://www.uv.es/ceaces/pdf/normal.pdf)** con su media y su anchura propias. Los perfiles, además, pueden tener **tamaños distintos**: uno puede ser la mayoría de la población y otro una minoría. Todos juntos forman una *mezcla de normales*, es decir, varias campanas de distinto tamaño superpuestas. El trabajo del modelo es **separar las campanas que están mezcladas**. La @fig-concepto dibuja exactamente esa idea. A la izquierda, lo que tú mides: una nube de individuos que, vista de golpe, parece una sola población. A la derecha, lo que en realidad hay debajo: esa misma nube son tres perfiles distintos, de distinto tamaño y forma, apilados. ```{r} #| label: fig-concepto #| echo: false #| fig-cap: "La misma población vista de dos formas. A la izquierda, lo que mides: una única distribución. A la derecha, los tres perfiles de los datos que la componen, cada uno con su tamaño (peso), su centro y su anchura. Las tres campanas ponderadas suman exactamente la curva de la izquierda." library(tidyverse) # --- Códigos de color del blog --------------------------------------------- col_general <- "gray20" # para la población / lo agregado col_grupos <- c("#6F4197", # morado "#2DD4CF", # turquesa "#8C2C7C") # magenta # --- Tres perfiles DISTINTOS: distinto centro, anchura y peso --------------- componentes <- tribble( ~perfil, ~media, ~sd, ~peso, "Perfil 1", -0.3, 1.2, 0.55, # mayoría, ancho, centrado "Perfil 2", 1.8, 0.9, 0.30, # mediano, a la derecha "Perfil 3", 3.8, 0.5, 0.15 # minoría, estrecho, lejos ) # Densidad de cada perfil, PONDERADA por su peso (así suman la población) dens_perfiles <- componentes |> crossing(x = seq(-4, 6, length.out = 800)) |> mutate(densidad = peso * dnorm(x, media, sd)) # La población es la suma de las tres densidades ponderadas dens_poblacion <- dens_perfiles |> group_by(x) |> summarise(densidad = sum(densidad), .groups = "drop") |> mutate(perfil = "Población") fig_datos <- bind_rows( dens_poblacion |> mutate(panel = "Lo que mides: una población"), dens_perfiles |> mutate(panel = "Lo que hay debajo: tres perfiles") ) |> mutate(panel = factor(panel, levels = c("Lo que mides: una población", "Lo que hay debajo: tres perfiles"))) ggplot(fig_datos, aes(x, densidad, fill = perfil, color = perfil)) + geom_area(alpha = 0.5, position = "identity", linewidth = 0.7) + facet_wrap(~ panel) + scale_fill_manual(values = c("Población" = col_general, "Perfil 1" = col_grupos[1], "Perfil 2" = col_grupos[2], "Perfil 3" = col_grupos[3])) + scale_color_manual(values = c("Población" = col_general, "Perfil 1" = col_grupos[1], "Perfil 2" = col_grupos[2], "Perfil 3" = col_grupos[3])) + labs(x = "Una medida cualquiera", y = "Densidad", fill = NULL, color = NULL) + theme_minimal(base_size = 13) + theme( legend.position="none" ) ``` Fíjate en el truco visual de la @fig-concepto: la curva gris de la izquierda **no es** una campana de verdad, es la suma de las tres de la derecha. Cada perfil contribuye según su peso, por eso ninguno sobresale por encima de la población: encajan debajo de ella. El reto del LPA es exactamente este: solo le damos lo de la izquierda y tiene que reconstruir lo de la derecha. Y la dificultad que va a marcar todo el artículo ya se intuye: [cuanto más se pisan las campanas, más difícil es separarlas.]{.underline} ## El reto: tres especies que el modelo no conoce Para poner el LPA a prueba con la respuesta correcta en la mano usamos `iris`, una base de datos incluida en R: 150 flores de **tres especies**: *setosa*, *versicolor* y *virginica*; hay 50 de cada tipo, con cuatro medidas continuas (largo y ancho del pétalo y del sépalo). ```{r} #| label: iris-vistazo #| echo: false data(iris) skimr::skim(iris) ``` [Nosotros sí sabemos la especie de cada flor, pero al modelo se la esconderemos.](.underline) ### Antes de modelar nada ¿Cómo se *sospecha* que hay perfiles? Aquí está la primera lección de mundo real, no de juguete. **Nadie mira un solo gráfico y decide que hay grupos.** El trabajo previo es buscar pistas. Un mapa de correlaciones **no** sirve para esto. La correlación mide qué variables se mueven juntas, no qué individuos se agrupan; es ciega a las relaciones no lineales y, peor aún, la propia mezcla de grupos fabrica correlaciones que no existen dentro de cada grupo. No es una herramienta de detección de perfiles. La primera herramienta que si nos vale la pena: **mirar la distribución de cada variable por separado**. La pista que delata grupos ocultos es la *multimodalidad*, es decir, que una variable no tenga una sola joroba sino dos o más. Recuerda la @fig-concepto: si dos perfiles están lo bastante separados, su suma deja de parecer una campana y le salen bultos. Esta vista se hace con un único gráfico faceteado, y funciona igual con 4 variables que con 15: ```{r} #| label: fig-marginales #| #| fig-cap: "Distribución de cada variable por separado. Las medidas del pétalo muestran dos jorobas claras (pista de grupos ocultos); las del sépalo, una sola. Esta misma receta sirve para escanear 15 variables de golpe." iris |> pivot_longer(-Species, names_to = "variable", values_to = "valor") |> ggplot(aes(valor)) + geom_density(fill = col_general, color = col_general, alpha = 0.4) + facet_wrap(~ variable, scales = "free") + labs(x = NULL, y = "Densidad") + theme_minimal(base_size = 13) ``` En la @fig-marginales, las dos medidas del pétalo tienen claramente **dos jorobas**: eso grita "aquí hay al menos dos grupos". Las del sépalo parecen una sola campana: solas, no te harían sospechar nada. Moraleja: una sola variable puede engañarte en cualquier dirección, por eso se miran todas. La segunda herramienta cruza las variables **de dos en dos**, para ver si la nube de puntos se apelmaza en grumos. Para un puñado de variables, la matriz de dispersión lo enseña todo a la vez (en la diagonal, cada variable a solas; fuera de ella, cada par): ```{r} #| label: fig-pares #| fig-cap: "Matriz de dispersión de las cuatro medidas, sin etiquetas de especie. En varios paneles se ve un grumo claramente despegado del resto: la sospecha de que hay grupos nace aquí, antes de ajustar ningún modelo." library(GGally) ggpairs( iris, columns = 1:4, lower = list(continuous = wrap("points", alpha = 0.5, color = "#243A6E", size = 0.8)), diag = list(continuous = wrap("densityDiag", color = "#243A6E", fill = "#243A6E", alpha = 0.3)), upper = list(continuous = wrap("cor", color = "gray20")) ) + theme_minimal(base_size = 10) ``` En la @fig-pares, varios paneles muestran un grumo despegado del resto de la nube. Sin saber todavía qué es, ya tienes una sospecha fundada: **hay estructura, probablemente al menos dos grupos**. ¿Y si tuvieras quince variables? La matriz de la @fig-pares se vuelve ilegible enseguida, pero el escaneo de densidades de la @fig-marginales escala sin problema y sigue siendo el caballo de batalla: la multimodalidad delata grupos haya las variables que haya. ### Aplastar muchas dimensiones en dos: UMAP y t-SNE Las dos herramientas anteriores enseñan los datos a trozos: una variable cada vez, o un par cada vez. Pero a nuestro ojo se le da de maravilla detectar grumos en una nube de puntos, y eso solo sabemos hacerlo en dos dimensiones. Eso es exactamente una **proyección a dos dimensiones**: aplastar el espacio de muchas variables sobre un plano intentando respetar las distancias. Piensa en la sombra de un objeto tridimensional sobre una pared. Pierdes información, pero si la luz incide bien, la silueta sigue siendo reconocible. Toda la cuestión está en desde qué ángulo iluminar. La forma más simple, el **análisis de componentes principales** (PCA), ilumina con luz recta: busca el ángulo que conserva la mayor variación posible y proyecta. Es rápido y fiel a la estructura global, pero al ser lineal no sabe "desenrollar" grupos separados por fronteras curvas; si dos perfiles se distinguen de forma no lineal, la sombra de PCA los puede dejar pegados. **UMAP** y **t-SNE** son las versiones no lineales del mismo truco. Miran *quién es vecino de quién* en el espacio original y construyen un mapa en 2D que respete esas vecindades. Por eso sacan a la luz grupos con formas retorcidas que PCA aplastaría. Ambas tienen un mando que regula cuán local es la mirada: en t-SNE se llama *perplexity*; en UMAP, *n_neighbors*; y conviene tenerlo presente, porque ese mando cambia el dibujo. Probémoslo con iris. Estandarizamos las cuatro medidas (para que ninguna domine solo por su unidad, algo que repetiremos al modelar) y pedimos las dos proyecciones: ```{r} #| label: fig-embedding #| fig-cap: "Mapa con al proyección en 2 variables de Las cuatro medidas de cada flor, tanto con UMAP y con t-SNE. El color (especie real) no estaría disponible en una exploración de verdad; lo añadimos para comprobar. Ambas técnicas cuentan lo mismo: setosa forma un grumo aislado, y versicolor y virginica comparten otro con una frontera difusa." library(uwot) library(Rtsne) # Debemos pasar los datos a formato matriz y normalizarlos para que estas técnicas funcionen X_mat <- iris |> dplyr::select(-Species) |> scale() |> as.matrix() set.seed(20242026) emb_umap <- umap(X_mat, n_neighbors = 15, min_dist = 0.1) emb_tsne <- Rtsne(X_mat, perplexity = 30, check_duplicates = FALSE)$Y proyecciones <- bind_rows( tibble(metodo = "UMAP", d1 = emb_umap[, 1], d2 = emb_umap[, 2], Especie = iris$Species), tibble(metodo = "t-SNE", d1 = emb_tsne[, 1], d2 = emb_tsne[, 2], Especie = iris$Species) ) proyecciones_2x2 <- bind_rows( proyecciones |> mutate(vista = "Lo que verías", grupo = "Sin etiqueta"), proyecciones |> mutate(vista = "La solución real", grupo = as.character(Especie)) ) |> mutate(vista = factor(vista, levels = c("Lo que verías", "La solución real"))) colores_2x2 <- c("Sin etiqueta" = col_general, setNames(col_grupos, levels(iris$Species))) ggplot(proyecciones_2x2, aes(d1, d2, color = grupo)) + geom_point(size = 2, alpha = 0.45) + facet_grid(vista ~ metodo, scales = "free") + scale_color_manual(values = colores_2x2, breaks = levels(iris$Species), name = "Especie real") + labs(x = "Dimensión 1", y = "Dimensión 2") + theme_minimal(base_size = 13) + theme( legend.position = "bottom", panel.border = element_rect(color = "gray70", fill = NA, linewidth = 0.5), panel.spacing = unit(1, "lines"), strip.background = element_rect(fill = "gray90", color = "gray70", linewidth = 0.5) ) ``` En la @fig-embedding las dos técnicas cuentan la misma historia: un grumo perfectamente aislado y otro, mayor, más enmarañado. Coloreando por la especie real (ese lujo que en un problema de verdad no tendrías, y que aquí usamos solo para comprobar) se ve que el núcleo de puntos solitario es *setosa* y que el bulto grande puntos son *versicolor* y *virginica* rozándose, sin una frontera limpia. Otra vez lo mismo: dos grupos clarísimos y un tercero que pelea por despegarse del segundo. La proyección no te da el número de perfiles, pero te deja la sospecha afinada: dos seguros, quizá tres. ### ¿Y si el solape se deshace en tres dimensiones? Para cerrar esta sección me queda una duda legítima: ¿y si en realidad las dos especies se separasen en una tercera que el plano ha pisado? La clásica de los físicos: "mi sistema funciona en 27 dimensiones". Podemos hacer la proyección en 3D, UMAP por ejemplo entrega un mapa en 3D con solo pedirle una dimensión más (`n_components = 3`), y con `plotly` lo volvemos interactivo: gíralo con el ratón y pasa el cursor por encima de cualquier flor para ver sus medidas. ```{r} #| label: fig-umap3d #| fig-cap: "UMAP en tres dimensiones. Gira la nube con el ratón y pasa el cursor sobre cualquier punto para ver su especie y sus cuatro medidas. setosa se mantiene aislada desde cualquier ángulo; versicolor y virginica siguen interpenetradas también en la tercera dimensión." library(plotly) set.seed(2024) emb_3d <- umap(X_mat, n_neighbors = 15, min_dist = 0.1, n_components = 3) datos_3d <- iris |> mutate( U1 = emb_3d[, 1], U2 = emb_3d[, 2], U3 = emb_3d[, 3], etiqueta = paste0( "Especie: ", Species, "<br>Sépalo: ", Sepal.Length, " x ", Sepal.Width, " cm", "<br>Pétalo: ", Petal.Length, " x ", Petal.Width, " cm" ) ) plot_ly( datos_3d, x = ~U1, y = ~U2, z = ~U3, color = ~Species, colors = col_grupos, text = ~etiqueta, hoverinfo = "text", type = "scatter3d", mode = "markers", marker = list(size = 4, opacity = 0.85) ) |> layout( scene = list( xaxis = list(title = "UMAP 1"), yaxis = list(title = "UMAP 2"), zaxis = list(title = "UMAP 3") ), legend = list(title = list(text = "Especie real")) ) ``` Puedes mover la @fig-umap3d a tu gusto con ambos clicks de ratón para hacerte la idea: *setosa* sigue siendo una isla desde cualquier ángulo, pero *versicolor* y *virginica* continúan pegadas también en la tercera dimensión. El solapamiento sigue pasando a pesar de la proyección de a una dimensión adicional. Ni el mapa 2D ni el 3D son *la verdad*, son resúmenes deformados de ella, en los que mires como mires, estas dos especies son bastante cercanas **si no tenemos más variables por las que separarlas**. La proyección 3D aporta una confirmación visual más, valiosísima cuando tienes quince variables y ningún diagrama capaz de mostrarlas todas a la vez, porque girar la nube es lo más cerca que vas a estar de "ver" el espacio entero. No quiero extender más esta sección con un caso de"y si hiciermos..." iterativo, en algún punto hay que acabar. Entonces, me toca repetir el aviso: UMAP y t-SNE **siempre** dibujan grumos bonitos: están diseñadas para agrupar, así que agrupan. Ni el tamaño de cada mancha ni la distancia entre manchas tienen que significar lo que parecen —dos grumos pegados pueden estar lejísimos en el espacio real. Por eso **una proyección es un generador de sospechas, nunca una prueba**. La conclusión vendrá de otras cosas la dará el modelo y sobretodo al final, el expertisement que se tenga del problema que se está tratando. AL menos sabemos que hay grupos, sabemos que hay algo, ya tiene sentido llamar al LPA. ## Calibrando un modelo de LPA ### Buscando el mejor modelo Vamos a usar **`tidyLPA`**, una interfaz que ordena todo el flujo del LPA en pasos limpios. Por dentro se apoya en el motor de `mclust`, pero nos deja una sintaxis cómoda y comparable. Antes de ajustar, **estandarizamos** las cuatro medidas (media 0, anchura 1) para que ninguna pese más que otra solo por su unidad. ```{r} #| label: tidylpa-preparacion library(tidyLPA) library(knitr) library(DT) indicadores <- iris |> dplyr::select(-Species) |> scale() |> as.data.frame() datatable(round(indicadores, 2), options = list(pageLength = 5)) ``` Entonces toca tomar una decisión que mucha gente pasa por alto: **qué forma le permitimos tener a cada campana**. Como dijimos antes, un LPA es una mezcla de gaussianas: cada perfil es una distribución normal. Pero... ¿podrían las campanas no ser normales, ser asimétricas o tener colas pesadas? En un LPA hasta el nivel que controlo, no. Que cada componente sea normal es precisamente lo que define al método; existen otras familias de mixturas: de distribuciones *t* para colas pesadas o normales sesgadas para asimetría; pero aquí la campana es siempre normal. Cabe recordar que "normal" no es una sola forma, es toda una familia. En una dimensión, una normal queda descrita por dos números: su media (dónde está el centro) y su varianza (cómo de ancha es). En varias dimensiones la media pasa a ser un punto en el espacio de datos, y la varianza se convierte en algo más rico, la **matriz de covarianzas**, este parámetro de modelo decide cuánta de esa libertad concedemos a cada distribución, y lo hace respondiendo a tres preguntas. - ¿Tienen todas las variables la misma anchura dentro de un perfil, o puede cada una tener la suya? - ¿Permitimos que las variables se correlacionen (que la campana se incline), o la obligamos a quedarse alineada con los ejes? - ¿Comparten todos los perfiles la misma forma, o puede cada perfil tener la suya propia? En un extremo, el modelo más **rígido** responde "todo igual y sin inclinar";Una forma rígida es austera y estable, pero, como vamos a ver enseguida, cuando le prohíbes a la campana inclinarse, el modelo se las apaña para capturar esa correlación que no puede modelar... inventándose perfiles de más. En el otro, el modelo más **flexible** responde "cada cual a su aire": cada perfil es su propio elipsoide, con su tamaño, su anchura por variable y su inclinación. Una forma más flexible se ajusta mejor a grupos reales que son elípticos o están inclinados, pero cuesta muchos más parámetros (la matriz de covarianzas crece deprisa con el número de variables) y pide más datos para estimarse bien. `tidyLPA` lo controla la forma de cada campana con argumento que es el "número de modelo". para este caso, me interesan dos extremos: - **Modelo 1 (rígido):** cada perfil es una campana "redonda", y se prohíbe que las variables estén correlacionadas dentro de un perfil. - **Modelo 6 (flexible):** cada perfil puede tener su propia forma y permitir que sus variables se correlacionen libremente. Ajustamos de 2 a 4 perfiles con ambos modelos, que eran el número de opciones plausibles que hemos detectado previamente, y miramos todos los indicadores de ajuste: ```{r} #| label: tidylpa-ajuste modelos_LPA <- indicadores |> estimate_profiles(2:4, models = c(1, 6)) datatable( get_fit(modelos_LPA) |> mutate(across(everything(),~round(.x,2))) , options = list(pageLength = 10) ) ``` ## El diagnóstico, criterio por criterio Elegir modelo no es mirar un número y sentenciar. Existe toda una familia de criterios, y casi nunca coinciden; el oficio está en entender qué pregunta responde cada uno y por qué discrepan. La tabla de ajuste mezcla además dos decisiones distintas, así que las separamos: primero **qué forma** de campana, y después **cuántos** perfiles. Empecemos por verlo todo junto. Giramos la tabla para tener cada modelo en una columna y cada métrica en una fila, que con pocas etiquetas y muchas métricas se lee mejor así: ```{r} #| label: tbl-diagnostico ajuste <- get_fit(modelos_LPA) ajuste |> filter(Classes >= 2) |> transmute( modelo = paste0("Modelo ", Model, " · ", Classes, "p"), `Log-verosimilitud` = LogLik, `Parámetros` = parameters, AIC, BIC, SABIC, CAIC, AWE, Entropía = Entropy, `Prob. mínima` = prob_min, `N grupo menor` = n_min * n, `BLRT (p)` = BLRT_p ) |> pivot_longer(-modelo, names_to = "Métrica", values_to = "valor") |> pivot_wider(names_from = modelo, values_from = valor) |> mutate(across(where(is.numeric), \(x) round(x, 2))) |> kable(caption = "Todos los modelos ajustados, con las métricas en filas. En los criterios de información más bajo es mejor; en entropía y probabilidad mínima, más alto.") ``` ### Primera decisión: ¿qué forma de campana? Para empezar debemos elegir un modelo de los 6 que hemos calibrado. La métrica más habitual por la que empezar a mirar esto es el **BIC** (Bayesian information criteria): resume en un solo número lo bien que un modelo ajusta los datos, poniéndole un coste por cada parámetro que usa. Se suele empezar el análisis por aquí porque premiaa un buen ajuste sin recompensar poner complejidad por las buenas. viendo su fórmula se entiende el criterio: $$\text{BIC} = \underbrace{-2\ln(\hat{L})}_{\text{lo mal que ajusta}} \;+\; \underbrace{k \cdot \ln(n)}_{\text{peaje por complejidad}}$$ Los dos sumandos son costes, así que **más bajo es mejor**. El primero baja cuando el modelo ajusta mejor; el segundo sube $\ln(n)$ por cada parámetro $k$, ajuste o no ajuste. El BIC es una factura y queremos pagar poco. Entonces el modelo flexible (6) gana al rígido (1) a cualquier número de perfiles en BIC. Con dos perfiles, 851 frente a 1260. Y eso tiene miga, porque el flexible es mucho más complejo: 41 parámetros frente a 13. Paga un peaje bastante grande, unos $28 \times \ln(150) \approx 140$ puntos más, y aun así el ajuste mejora 549 puntos, muchísimo más de lo que cuesta esa complejidad. **¿Por qué los modelos flexibles ajustan mejor?** Por que **deja que las cuatro medidas se correlacionen dentro de cada grupo**, cosa que el ajuste rígido lo prohíbe. Significa que esas correlaciones existen y son fuertes. Y eso es justo lo que buscábamos: unas variables que se mueven juntas y forman grupos con forma y orientación propias son la huella de que hay **perfiles de verdad**. El rígido, al obligar a cada perfil a ser una campana redonda, se niega a verlos como son. Decisión cerrada: **nos quedamos con el flexible**. Todo lo que sigue ocurre dentro de él. nota: el flexible no gana solo en BIC; gana en los cinco criterios de información, así que la decisión es unánime y no depende de cuál miremos. Pa ra este análisis los miré todos, pero pensé que enseñar los cinco para que digan lo mismo cinco veces era paja, de modo que me centré en uno. ### Segunda decisión: ¿cuántos perfiles? LLegamos al *quid* de la cuestión. Aquí no vale hacer trampa diciendo "pongo 3 porque sé que son 3". En un problema real no se conoces el número, así que hay que tomar una decisión final entre criterios que muy posiblemente que se contradicen. ::: {.panel-tabset} ### Criterios de información Un *criterio de información* resume el modelo en un solo número que combina dos cosas: cómo de bien ajusta y cuántos parámetros gasta, penalizando lo segundo para no premiar la complejidad porque sí. Más bajo, mejor. Lo que distingue a unos de otros es la **dureza** de esa penalización. De más blando (tolera más perfiles) a más duro (los recorta): - **AIC** — *Akaike Information Criterion* (criterio de información de Akaike, por Hirotugu Akaike). El más permisivo: su peaje es fijo y pequeño (2 por parámetro) y no crece con el tamaño de la muestra. Tiende a quedarse con soluciones más grandes, así que cuidado, porque en muestras decentes suele sobreestimar el número de perfiles. - **SABIC** — *Sample-size Adjusted BIC* (BIC ajustado por el tamaño muestral; también ssBIC o aBIC). El BIC, pero con el peaje suavizado mediante una corrección por tamaño muestral. Más indulgente que el BIC, más prudente que el AIC. Suele ser la referencia recomendada cuando esos dos discrepan. - **BIC** — *Bayesian Information Criterion* (criterio de información bayesiano; también Schwarz Information Criterion o SIC, por Gideon Schwarz). El de referencia, el que usamos en el diagnóstico del tipo de modelo. Su peaje por parámetro crece con la muestra: cuantos más datos, más severo con la complejidad. Por construcción, es conservador. - **CAIC** — *Consistent AIC* (AIC consistente, de Bozdogan). Lleva "AIC" en el nombre, pero es en realidad el BIC con un puntito más de penalización: siempre vota lo mismo que el BIC, o algo más austero. - **AWE** — *Approximate Weight of Evidence* (peso aproximado de la evidencia, de Banfield y Raftery). El más duro de todos. Penaliza no solo la complejidad sino también la mala clasificación (lo difusos que quedan los grupos), así que castiga doblemente las soluciones con muchos perfiles solapados. Es el que más recorta y el más reacio a añadir grupos. ```{r} #| label: fig-ic #| fig-cap: "Cinco criterios de información para el modelo flexible. El punto turquesa marca el mínimo (la elección de ese criterio)." ic_largo <- ajuste |> filter(Model == 6) |> dplyr::select(Classes, AIC, BIC, SABIC, CAIC, AWE) |> pivot_longer(-Classes, names_to = "criterio", values_to = "valor") minimos <- ic_largo |> group_by(criterio) |> slice_min(valor, n = 1) ggplot(ic_largo, aes(Classes, valor)) + geom_line(color = "gray60", linewidth = 0.8) + geom_point(color = "gray60", size = 2) + geom_point(data = minimos, color = "#2DD4CF", size = 4) + facet_wrap(~ criterio, scales = "free_y") + scale_x_continuous(breaks = 2:4) + labs(x = "Número de perfiles", y = "Valor (más bajo = mejor)") + theme_minimal(base_size = 12) + theme(panel.border = element_rect(color = "gray70", fill = NA, linewidth = 0.5), panel.spacing = unit(0.8, "lines"), strip.background = element_rect(fill = "gray90", color = "gray70", linewidth = 0.5)) ``` Los criterios **no coinciden**, y eso es el hallazgo. - Los de penalización dura (BIC, CAIC y AWE) tocan fondo en **2 perfiles**. - Los más suaves (AIC y SABIC) lo hacen en **3**. Todo el cisma está en el salto de 2 a 3 perfiles: cuesta 21 parámetros y mejora el ajuste en unos 68 puntos. ¿De quién me fío? El $\ln(n)$ del BIC sale de una aproximación pensada para muestras grandes, y con 150 flores y un modelo de 62 parámetros estamos lejos de ese ideal: en ese régimen el BIC penaliza de más y tiende a quedarse corto de perfiles. El SABIC nació precisamente para corregir un problema como este, ajustando el peaje al tamaño real de la muestra, y por eso recupera mejor el número de grupos en muestras pequeñas o medianas, que es justo nuestro caso. Conclusión: 3 perfiles Toma algo de ventaja como respuesta definitiva. ### BLRT El **BLRT** (*Bootstrap Likelihood Ratio Test*, test de razón de verosimilitud por bootstrap) es una prueba estadística formal a una pregunta muy concreta: **¿el perfil que acabo de añadir mejora el ajuste más de lo que mejoraría por puro azar?** Para entenderlo hay que ver el problema que resuelve. Cuando pasas de 2 a 3 perfiles, el ajuste *siempre* mejora un poco, aunque el tercer perfil sea basura (*overfitting*). La pregunta no es *si* mejora, sino *cuánto*, y sobre todo: ¿esa mejora es de las que solo da una estructura real, o es del tamaño de la que saldría por casualidad partiendo en dos un grupo que era uno solo? El BLRT responde a ese "*cuánto*" mediante **bootstrap**. La mecánica, en cuatro pasos: 1. [Mide la mejora de ajuste real al pasar]{.underline} de $k$ a $k+1$ perfiles. La medida es el estadístico de razón de verosimilitudes, $LR = -2\ln(\hat{L}_k / \hat{L}_{k+1})$, que crece cuanto más mejora el ajuste al añadir el perfil. 2. [Adopta la **hipótesis nula**]{.underline}; $H_0$: el modelo verdadero es el de menos perfiles, el de $k$ (frente a la alternativa de que haga falta uno más, $k+1$). Asumiéndola cierta, genera su distribución correspondiente simulando muchos conjuntos de datos a partir del modelo de $k$ perfiles ya estimado. 3. [A cada conjunto simulado le ajusta los modelos de k y k+1 perfiles]{.underline} y recalcula su $LR$. Esto reconstruye la **distribución de $LR$ bajo $H_0$**: el catálogo de mejoras que produce el azar cuando no hay un perfil más que encontrar. 4. [El **valor p** es la proporción de réplicas cuyo $LR$ iguala o supera al observado]{.underline}, es decir, la probabilidad de ver una mejora como la nuestra si $H_0$ fuera cierta. Así que si, la regla de del estadístico es la clásica: **p pequeño**, la mejora es demasiado grande para ser casualidad, el perfil extra es real; **p grande**, la mejora cabe holgadamente dentro de lo que da el azar, el perfil no aporta. ```{r} #| label: fig-blrt #| fig-cap: "Valor p del BLRT al añadir cada perfil secuencialmente. Por debajo de la línea, la mejora supera al azar y el perfil aporta; por encima, no." ajuste |> filter(Model == 6, !is.na(BLRT_p)) |> ggplot(aes(Classes, BLRT_p)) + geom_hline(yintercept = 0.05, linetype = "dashed", color = "gray50") + geom_line(color = "gray60", linewidth = 0.8) + geom_point(aes(color = BLRT_p < 0.05), size = 4) + scale_color_manual(values = c("TRUE" = "#2DD4CF", "FALSE" = "#8C2C7C"), labels = c("TRUE" = "Aporta (p < 0,05)", "FALSE" = "No aporta"), name = NULL) + scale_x_continuous(breaks = 2:4) + labs(x = "Número de perfiles", y = "Valor p del BLRT (k frente a k−1)") + theme_minimal(base_size = 12) + theme(legend.position = "bottom") ``` El BLRT no penaliza a ojo: contrasta la mejora real contra su propia distribución por azar. No es casualidad que, en los estudios de simulación que enfrentan a todos estos indicadores, sea el que más veces acierta el número real de grupos, por delante del BIC. Cuando ambos discrepan, hay buenas razones para hacerle caso. Mirando la @fig-blrt hay una premisa interesante, Pasar de 1 a 2 perfiles supera al azar (p ≈ 0,01) y pasar de 2 a 3 **también** (p ≈ 0,01), pero pasar de 3 a 4 **no** (p = 0,65): esa mejora es esperable por casualidad, así que el cuarto perfil no añade estructura, solo ruido. **El BLRT señala que lo ideal son 3 perfiles** ### Separación y confianza El LPA no asigna cada flor a un perfil directamente sino que como todo modelo probabilístico le da una **probabilidad de pertenencia** a cada uno (70 % al perfil 1, 20 % al 2, 10%al 3 por ejemplo). Cuando esas probabilidades son nítidas (casi 100 % a un perfil y casi 0 % al resto), los grupos están bien separados; cuando quedan repartidas (50 % y 50 %), el modelo está dudando. Dos métricas resumen esa nitidez desde ángulos distintos: - La **entropía** mira el conjunto: resume, en un solo número de 0 a 1, cómo de seguras son *todas* las asignaciones a la vez. Vale 1 cuando cada flor pertenece a su perfil sin ambigüedad y cae hacia 0 cuando las pertenencias están, en promedio, repartidas. Es la separación global de la solución. El umbral habitual es ≥ 0,80 para considerar que la clasificación es fiable. Por encima de 0,80 se habla de separación buena; por encima de 0,90, excelente - La **probabilidad mínima** mira el eslabón más débil: de entre todos los perfiles, toma el que peor se clasifica y mira con qué probabilidad media se asigna correctamente a sus miembros. Una entropía alta puede esconder un perfil concreto mal definido; esta métrica lo destapa, porque se fija justo en el peor. La referencia es ≥ 0,70 por clase; es decir, ningún perfil debería clasificar a sus miembros con una certeza media inferior a 0,70. ```{r} #| label: fig-separacion #| fig-cap: "Separación entre perfiles (entropía) y seguridad de la clasificación (probabilidad mínima) para el modelo flexible." ajuste |> filter(Model == 6, !is.na(Entropy)) |> dplyr::select(Classes, Entropía = Entropy, `Prob. mínima` = prob_min) |> pivot_longer(-Classes, names_to = "métrica", values_to = "valor") |> ggplot(aes(Classes, valor, color = métrica)) + geom_line(linewidth = 1) + geom_point(size = 3) + scale_color_manual(values = c("Entropía" = "#6F4197", "Prob. mínima" = "#2DD4CF"), name = NULL) + scale_x_continuous(breaks = 2:4) + coord_cartesian(ylim = c(0.85, 1)) + labs(x = "Número de perfiles", y = "Valor (más alto = mejor)") + theme_minimal(base_size = 12) + theme(legend.position = "bottom") ``` Ambas métricas coinciden: en 3 perfiles la separación global sigue siendo casi perfecta (entropía 0,97) y, mirando el perfil peor definido, la clasificación sigue siendo segurísima (probabilidad mínima 0,98). El tercer perfil es un grupo nítido, y además ninguno de los tres queda mal delimitado. Solo al llegar a 4 empiezan ambas a ceder, señal de que ahí ya se está forzando una división que los datos no sostienen. Conclusión: los criterios de "*separación y confianza*"" indica que 3 perfiles son lo ideal. ### Tamaño de grupos La última comprobación yo la he llamado "de cordura": ¿qué tamaño tiene el grupo más pequeño de cada solución? Un perfil que tenga cuatro gatos suele ser un artefacto, no un tipo real. No he sido capaz de encontrar reglas universales con una base teórica fundamentada (desconfía de quien te venda un número mágico), pero sí pistas bastante coherentes. - En **términos relativos**, el umbral más citado es que cada clase reúna al menos un **5 %** de la muestra (algunos lo suben al 10 % para ser más exigentes). - En **términos absolutos**, la referencia habitual es un mínimo de **30 casos** por clase. Y este segundo manda sobre el primero, porque el porcentaje engaña con el **tamaño del dataset**: un 5 % de 5.000 son 250 individuos, de sobra para estimar un grupo con estabilidad; pero un 5 % de 200 son 10, un grupo frágil aunque cumpla el porcentaje. El motivo de fondo no es estético sino de estimación: cada perfil calcula sus propias medias y covarianzas a partir de los casos que le tocan, y con el modelo flexible, que gasta muchos parámetros por grupo, un puñado de casos no da para fijarlos; el grupo resultante es inestable, cambia con la semilla y no replica. Esa inestabilidad *es* la alucinación. Con todo, un grupo pequeño no es automáticamente falso: las subpoblaciones raras existen y a veces son lo más interesante. La regla no es "bórralo", sino "un grupo pequeño necesita más pruebas para creértelo". El de 12 flores las suspende todas: supera por poco el 5 %, pero se queda por debajo de los 25 casos absolutos, el BLRT ya lo daba por no significativo y aparece justo cuando la predicción empeora. No lo descartamos por un umbral aislado, sino porque tres señales independientes apuntan a lo mismo. ```{r} #| label: fig-tamano #| fig-cap: "Tamaño del grupo más pequeño en cada solución. El grupúsculo de 12 flores con 4 perfiles es la firma del sobreajuste." ajuste |> filter(Model == 6, !is.na(n_min)) |> mutate(n_menor = round(n_min * n)) |> ggplot(aes(factor(Classes), n_menor)) + geom_col(aes(fill = n_menor < 20), width = 0.6) + geom_text(aes(label = n_menor), vjust = -0.4, size = 4) + scale_fill_manual(values = c("FALSE" = "#6F4197", "TRUE" = "#8C2C7C"), guide = "none") + labs(x = "Número de perfiles", y = "Tamaño del grupo más pequeño (nº de flores)") + theme_minimal(base_size = 12) ``` Con 2 y 3 perfiles, el grupo menor tiene 50 y 45 flores: grupos grandes y respetables. Con 4 se desploma a **12** (un 8 % de la muestra): el grupúsculo anecdótico que delata que el modelo ha empezado a inventarse divisiones. Encaja con el BLRT, que en 4 perfiles ya decía "no aporta". ### Predicción por validación cruzada Todos los criterios anteriores miran lo mismo: cómo de bien encaja el modelo en los datos que ya tiene. La validación cruzada mide algo distinto y complementario: **cómo de bien predice datos que no ha visto**. La idea es sencilla. Partimos la muestra en cinco trozos; ajustamos el modelo con cuatro de ellos y medimos cómo de verosímiles le resultan las flores del quinto, que se quedó fuera. Rotamos el trozo apartado, y repetimos varias veces con particiones distintas para promediar el azar. Un modelo que se ha inventado perfiles de más ajustará de maravilla sus propios datos, pero predecirá peor los ajenos: el sobreajuste se paga fuera de casa. ```{r} #| label: fig-cv #| fig-cap: "Capacidad de predecir flores no vistas, por número de perfiles (validación cruzada repetida, modelo flexible). Más arriba es mejor. Las barras son el error estándar entre repeticiones." library(rsample) library(mclust) set.seed(2024) folds <- vfold_cv(indicadores, v = 5, repeats = 5) # Log-verosimilitud media por flor en el trozo apartado loglik_fuera <- function(split, G) { entrena <- as.matrix(analysis(split)) prueba <- as.matrix(assessment(split)) m <- Mclust(entrena, G = G, modelNames = "VVV", verbose = FALSE) if (is.null(m)) return(NA_real_) mean(dens(modelName = m$modelName, data = prueba, parameters = m$parameters, logarithm = TRUE)) } cv <- tidyr::expand_grid(G = 2:4, split = folds$splits) |> mutate(loglik = purrr::map2_dbl(split, G, loglik_fuera)) cv_resumen <- cv |> group_by(G) |> summarise(media = mean(loglik, na.rm = TRUE), ee = sd(loglik, na.rm = TRUE) / sqrt(sum(!is.na(loglik))), .groups = "drop") ggplot(cv_resumen, aes(G, media)) + geom_line(color = "gray60", linewidth = 0.8) + geom_pointrange(aes(ymin = media - ee, ymax = media + ee, color = G == 3), size = 0.9) + scale_color_manual(values = c("TRUE" = "#2DD4CF", "FALSE" = "gray40"), guide = "none") + scale_x_continuous(breaks = 2:4) + labs(x = "Número de perfiles", y = "Log-verosimilitud media fuera de muestra") + theme_minimal(base_size = 13) ``` El resultado de la @fig-cv es sobrio, y conviene leerlo con honestidad. Las soluciones de 2 y 3 perfiles predicen **prácticamente igual de bien** (sus intervalos se solapan), mientras que la de 4 predice claramente peor. De ahí salen dos lecturas. La primera: pasar a 4 perfiles **sí sobreajusta**, y la validación cruzada lo castiga, exactamente como avisaba el BLRT. La segunda, más sutil: elegir 3 en lugar de 2 **no cuesta nada en capacidad predictiva**; el tercer perfil no es un adorno que estropee la generalización, generaliza tan bien como la opción conservadora. Lo que la validación cruzada no hace es romper el empate entre 2 y 3, y no debemos fingir que sí: para distinguirlos hacía falta el BLRT. Las dos pruebas se reparten el trabajo: el BLRT establece que el tercer perfil es una mejora real y no ruido, y la validación cruzada confirma que esa mejora no se paga prediciendo peor. ::: ## Resultados Reuniendo las cinco vistas del diagnóstico (siempre dentro del modelo flexible, porque el rígido quedó descartado por ajustar mucho peor en todos los criterios), los puntos que sostienen la decisión son: - **Forma de campana.** El modelo flexible bate al rígido en todos los criterios de información y a cualquier número de perfiles: las cuatro medidas se correlacionan dentro de cada grupo, y modelar esa correlación importa. - **Criterios de información, divididos.** Los de penalización dura (BIC, CAIC, AWE) tocan fondo en 2 perfiles; los más suaves (AIC, SABIC) en 3. SABIC es la métrica más idonea para esta situación, de modo que conservamos en 3. - **BLRT, decisivo.** El tercer perfil mejora el ajuste de forma significativa (p ≈ 0,01); el cuarto, no (p = 0,65). El indicador más fiable para esta pregunta traza la frontera en 3. - **Separación y confianza, intactas en 3.** Entropía 0,97 y probabilidad mínima 0,98 para 3 perfiles: grupos nítidos y clasificación segura. A partir de 4 empiezan a degradarse. - **Tamaño, la firma del sobreajuste.** El grupo más pequeño tiene 45 flores con 3 perfiles, cuando se va a 4 perfiles, el grupo minoritario se desploma a 12: el grupúsculo anecdótico que delata divisiones inventadas. - **Predicción.** 2 y 3 perfiles generalizan igual de bien; 4 generaliza peor. Elegir 3 no cuesta capacidad predictiva, y 4 la pierde. El veredicto que sostienen estas vistas es **3 perfiles**: la única solución que el BLRT respalda como mejora real, que mantiene separación y tamaños sólidos, y que predice tan bien como la alternativa conservadora. Los criterios que preferían 2 lo hacían solo por penalizar con dureza la complejidad del modelo flexible; ninguno señala que el tercer perfil sea espurio. Y el 4 queda descartado por triple vía: BLRT no significativo, grupo residual minúsculo y peor predicción. ## La prueba del algodón: contra la hoja de respuestas Hemos decidido ajustar un **modelos final de LPA flexible de 3 perfiles**, solo con lo que sabemos de los datos. Empecemos por ajustar el modelo final: ```{r} #| label: tidylpa-tres-ajuste #| code-fold: false # Ajustamos el modelo elegido: 3 perfiles, modelo flexible (6) modelo_final <- indicadores |> estimate_profiles(3, models = 6) # get_data() devuelve, para cada individuo, sus valores y el resultado de la # clasificación: la columna 'Class' (el perfil asignado) y las 'CPROB_*' # (probabilidad de pertenecer a cada perfil). Miramos la tabla completa. datatable( get_data(modelo_final) |> mutate(across(where(is.numeric), \(x) round(x, 2))), options = list(pageLength = 5), caption = "Salida del modelo: valores, perfil asignado (Class) y probabilidades de pertenencia (CPROB)" ) ``` Como en `iris` sí conocemos la verdad, podemos permitirnos comprobar si esa decisión, tomada a ciegas, era la acertada. Empecemos por ajustar la matriz de confusión entre el modelo y la realidad: ```{r} #| label: tidylpa-tres-confusion # La columna 'Class' es el perfil asignado a cada flor. La pegamos a los datos # originales (get_data mantiene el mismo orden de filas) como factor. iris_res <- iris |> mutate(Perfil = factor(get_data(modelo_final)$Class)) # Matriz de confusión: cruzamos la especie real (que el modelo NO vio) con el # perfil que ha estimado. La diagonal son los aciertos; fuera, las confusiones. library(gt) table(Especie_real = iris$Species, Perfil = iris_res$Perfil) |> as.data.frame.matrix() |> tibble::rownames_to_column("Especie real \\ Perfil") |> gt(rowname_col = "Especie real \\ Perfil") |> tab_header(title = "Matriz de confusión: especie real frente a perfil estimado") |> tab_spanner(label = "Perfil estimado", columns = everything()) |> data_color( columns = where(is.numeric), fn = scales::col_numeric(c("white", "#2DD4CF"), domain = NULL) ) |> tab_options(table.align = "center", table.width = pct(60)) |> cols_align(align = "center", columns = where(is.numeric)) ``` Los tres perfiles se corresponden casi uno a uno con las tres especies. Para ponerle número a ese "casi" usamos el **índice de Rand ajustado** (ARI): un valor entre 0 y 1 que mide cómo de bien coincide la agrupación del modelo con la verdadera, sin que importe cómo se numeren los grupos. 1 es coincidencia perfecta; 0 es como tirar los dados. En este caso, el valor ronda `r round(mclust::adjustedRandIndex(iris_res$Perfil, iris$Species), 2)`: una recuperación notable. Ahora, para ver dónde falla, la @fig-pares-diagnostico dibuja cada variable con su perfil estimado; fíjate en cómo los desacuerdos se concentran justo en la frontera entre las dos especies que ya sabíamos solapadas: ```{r} #| label: fig-pares-diagnostico #| fig-cap: 'Las cuatro medidas cruzadas dos a dos. Color = perfil estimado por el modelo. Los puntos con borde negro grueso son los mal clasificados. Sin conocer la naturaleza de las variables, esta vista te dice de un vistazo cómo de limpios son los perfiles y dónde se confunden.' library(GGally) library(dplyr) # Correspondencia perfil <-> especie (la moda de cada perfil) correspondencia <- iris_res |> dplyr::count(Perfil, Species) |> dplyr::group_by(Perfil) |> dplyr::slice_max(n, n = 1) |> dplyr::ungroup() |> dplyr::select(Perfil, especie_esperada = Species) iris_plot <- iris_res |> dplyr::left_join(correspondencia, by = "Perfil") |> dplyr::mutate(mal_clasif = Species != especie_esperada) # Panel de dispersión: puntos por perfil + disonantes con borde negro lower_fn <- function(data, mapping, ...) { ggplot(data, mapping) + geom_point(aes(fill = Perfil), shape = 21, color = "gray75", size = 1.8, stroke = 0.25, alpha = 0.85) + geom_point(data = function(x) dplyr::filter(x, mal_clasif), aes(fill = Perfil), shape = 21, color = "black", size = 2.8, stroke = 1.1) + scale_fill_manual(values = col_grupos) } # Diagonal: densidad por perfil diag_fn <- function(data, mapping, ...) { ggplot(data, mapping) + geom_density(aes(fill = Perfil), color = NA, alpha = 0.5) + scale_fill_manual(values = col_grupos) } ggpairs( iris_plot, columns = c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width"), lower = list(continuous = lower_fn), diag = list(continuous = diag_fn), upper = list(continuous = wrap("cor", color = "gray30")), legend = c(3, 1), progress = FALSE ) + theme_minimal(base_size = 11) + theme( legend.position = "bottom", panel.border = element_rect(color = "gray70", fill = NA, linewidth = 0.4), panel.spacing = unit(0.5, "lines"), strip.background = element_rect(fill = "gray90", color = "gray70", linewidth = 0.4) ) ``` ```{r} #| label: tbl-perfiles-comparados library(tidyr) bind_rows( iris_plot |> transmute(grupo = as.character(Species), Sepal.Length, Sepal.Width, Petal.Length, Petal.Width, fuente = "Especie real"), iris_plot |> transmute(grupo = as.character(especie_esperada), Sepal.Length, Sepal.Width, Petal.Length, Petal.Width, fuente = "Perfil estimado") ) |> pivot_longer(Sepal.Length:Petal.Width, names_to = "Variable", values_to = "valor") |> group_by(Grupo = grupo, Variable, fuente) |> summarise(resumen = sprintf("%.2f ± %.2f", mean(valor), sd(valor)), .groups = "drop") |> pivot_wider(names_from = fuente, values_from = resumen) |> arrange(Grupo, Variable) |> knitr::kable(caption = "Media ± desviación típica (cm) por variable: especie real frente al perfil estimado que le corresponde. Casi idénticas: el modelo no solo agrupa bien, reconstruye los estadísticos de cada grupo.") ``` ## ¿Hemos aprendido algo? Hemos respondido al título. Toca hacer sumario de los descubrimientos y aprendizajes. - **Sí: el LPA puede redescubrir las tres especies sin que nadie se las diga.** Los modelos de la tienen aplicaciones factibles. - No son modelos muy atuomáticos: lo que hace fácil o difícil este problema **no es la cantidad de datos, sino cómo de separados estén los perfiles**. Aquí, *Setosa*, estaba bien apartada, se detectaba sola casi. En cambio, *versicolor* y *virginica*, estaba solapadas, solo se separaron una vez guiamos al modelo, razonando sobre toda una batería de métricas para la decisión. - Los modelos no son perfectos. La distancia entre perfiles es el factor que de verdad importa. - **En un problema real no hay hoja de respuestas**. Adjudicamos entre criterios lo mejor que pudimos, pero nadie nos confirmó por detrás que hubiéramos acertado. Con este dataset hemos hecho trampa: teníamos la respuesta correcta, así que pudimos comprobar que elegir 3 era acertado. ¿Hasta dónde podemos fiarnos de esa adjudicación cuando no hay red de seguridad? ¿Qué grupos somos capaces de detectar de verdad, y a partir de qué punto el modelo empieza, otra vez, a ver caras en la tostada? Creo que en otro post voy a expandir esta idea: poner a prueba un LPA cuando no tienes la solución delante, jugando con la separación y con el tamaño de los grupos para medir hasta dónde llega lo que el modelo puede ver. --- *Tu estadístico de confianza y, cada día, el de menos gente.* *Adiós. Os estimo.*