1 Planteamiento revisado de la pregunta predictiva

La exploración de las diferencias entre las estructuras de las flores de tres especies del género Iris, y el posterior análisis confirmatorio a través de pruebas de hipótesis (análisis de varianza y contrastes post hoc), pusieron en evidencia que la distribución de las dimensiones de la especie setosa se colocaba, en todos los casos, lejos de las respectivas distribuciones de sus hermanas de género (versicolor y virginica), a través de una brecha importante que permite inferir una clara diferenciación, sobre todo de los pétalos, entre las tres especies.

En efecto, discriminar a los especímenes de I. setosa por el tamaño de sus pétalos constituye una tarea bastante segura, lo cual libera al analista de mayor procedimiento, pues no hay traslape alguno en los valores de estas dimensiones (longitud y anchura), con una brecha ante su pariente más cercana de varios milímetros (entre dos y siete mm, en términos generales).

Sin embargo, también se hizo evidente el traslape de las distribuciones de las otras dos especies (I. versicolor e I. virginica), por lo que la pregunta predictiva se plantea de la manera siguiente: ¿es posible predecir de manera razonable la especie de una planta del género Iris (ya sea I. versicolor o I. virginica) a partir de las dimensiones de sus flores?

2 Dataset y variables seleccionadas

2.1 Carga de datos y paquetes pertinentes

Para este análisis se emplearon varios paquetes de R:

  1. ggplot2 (Wickham 2016),
  2. dplyr (Wickham et al. 2023) y
  3. caret (Kuhn 2008).
library(ggplot2)
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(caret)
## Cargando paquete requerido: lattice

Asimismo, se cargaron los datos generados por Anderson (Anderson 1936), y ampliamente referidos por Fisher (Fisher 1936), los cuales se incluyen en la versión básica de R.

iris_aae <- read.csv("iris_aae.csv", header = TRUE, row.names = "id")

2.2 Exploración del lote de datos

2.2.1 Casos iniciales

head(iris_aae)

2.2.2 Casos finales

tail(iris_aae)

3 Exploración de varios modelos

3.1 Creación de una variable dicotómica para separar a I. setosa

iris_aae <- iris_aae %>%
  mutate(setosa_vs_not = ifelse(variety == "Setosa", 1, 0))

3.1.1 Primeros casos

head(iris_aae)

3.1.2 Últimos casos

tail(iris_aae)

3.1.3 Todas las variables predictoras

modelo_log_setosaVsNot.Todas <- glm(setosa_vs_not ~ petal_length + petal_width +
                                  sepal_length + sepal_width, data = iris_aae,
                                  family = binomial)
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary.glm(modelo_log_setosaVsNot.Todas)
## 
## Call:
## glm(formula = setosa_vs_not ~ petal_length + petal_width + sepal_length + 
##     sepal_width, family = binomial, data = iris_aae)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)
## (Intercept)     -16.946 457457.096       0        1
## petal_length    -20.088 107724.592       0        1
## petal_width     -21.608 154350.613       0        1
## sepal_length     11.759 130504.039       0        1
## sepal_width       7.842  59415.383       0        1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1.9095e+02  on 149  degrees of freedom
## Residual deviance: 3.2940e-09  on 145  degrees of freedom
## AIC: 10
## 
## Number of Fisher Scoring iterations: 25

3.1.4 Sólo las variables que se diferencian en el análisis exploratorio

modelo_log_setosaVsNot.SoloDiferentes <- glm(setosa_vs_not ~ petal_length +
                                            petal_width, data = iris_aae,
                                            family = binomial)
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary.glm(modelo_log_setosaVsNot.SoloDiferentes)
## 
## Call:
## glm(formula = setosa_vs_not ~ petal_length + petal_width, family = binomial, 
##     data = iris_aae)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)      69.45   43042.94   0.002    0.999
## petal_length    -17.60   43449.07   0.000    1.000
## petal_width     -33.89  115850.84   0.000    1.000
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1.9095e+02  on 149  degrees of freedom
## Residual deviance: 5.1700e-09  on 147  degrees of freedom
## AIC: 6
## 
## Number of Fisher Scoring iterations: 25

Comentario: en ambos modelos, tanto con las cuatro variables disponibles como con las seleccionadas por su poder de discriminación, se lanza una advertencia de que las probabilidades de acertar son perfectas. En realidad, el procedimiento no hace más que la exploración previa: constatar que la distribución de I. setosa no se traslapa con las de sus parientes de género, y que, inclusive, hay una brecha que permite al taxónomo separar con las solas medidas del pétalo a los ejemplares.

Optaremos ahora por trabajar con la especie de mayor tamaño (I. virginica), cuyas distribuciones se traslapan, en efecto, con las de I. versicolor. Sin embargo, conservaremos en el esquema analítico los datos de la especie pequeña I. setosa, con fines de control y de apreciación del efecto de su presencia comparativa.

3.2 Creación de variable dummy con respecto a la variedad virginica

iris_aae <- iris_aae %>%
  mutate(virgi_vs_not = ifelse(variety == "Virginica", 1, 0))

3.2.1 Primeros casos

head(iris_aae)

3.2.2 Últimos casos

tail(iris_aae)

3.2.3 Todas las variables predictoras para I. virginica

modelo_log_virgiVsNot.Todas <- glm(virgi_vs_not ~ petal_length + petal_width +
                                      sepal_length + sepal_width,
                                    data = iris_aae, family = binomial)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary.glm(modelo_log_virgiVsNot.Todas) # resumen del modelo
## 
## Call:
## glm(formula = virgi_vs_not ~ petal_length + petal_width + sepal_length + 
##     sepal_width, family = binomial, data = iris_aae)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   -42.638     25.708  -1.659   0.0972 .
## petal_length    9.429      4.737   1.990   0.0465 *
## petal_width    18.286      9.743   1.877   0.0605 .
## sepal_length   -2.465      2.394  -1.030   0.3032  
## sepal_width    -6.681      4.480  -1.491   0.1359  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 190.954  on 149  degrees of freedom
## Residual deviance:  11.899  on 145  degrees of freedom
## AIC: 21.899
## 
## Number of Fisher Scoring iterations: 12

Comentario: Como se constata en el modelo resumido, los coeficientes asociados con las dimensiones de los sépalos no son diferentes de cero, pero sí los de las dimensiones de los pétalos, por lo que procedemos a modelar con base en las solas variables del pétalo.

3.2.4 Sólo las variables significativas para I. virginica

modelo_log_virgiVsNot.SoloDiferentes <- glm(virgi_vs_not ~ petal_length +
                                              petal_width, data = iris_aae,
                                            family = binomial)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary.glm(modelo_log_virgiVsNot.SoloDiferentes) # resumen del modelo
## 
## Call:
## glm(formula = virgi_vs_not ~ petal_length + petal_width, family = binomial, 
##     data = iris_aae)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -45.272     13.611  -3.326 0.000881 ***
## petal_length    5.755      2.306   2.496 0.012575 *  
## petal_width    10.447      3.756   2.782 0.005409 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 190.954  on 149  degrees of freedom
## Residual deviance:  20.564  on 147  degrees of freedom
## AIC: 26.564
## 
## Number of Fisher Scoring iterations: 10

Comentario: en ausencia de las variables con coeficientes de cero, las variables restantes, las del pétalo, separan con claridad a los ejemplares de la I. virginica; inclusive, con cierta sorpresa se tiene que la anchura del pétalo lo hace con mayor fuerza que la longitud.

4 Predicciones

predicciones <- ifelse(predict(modelo_log_virgiVsNot.SoloDiferentes,
                               type = "response") > 0.5, 1, 0)
predicciones
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
##   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
##  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
##   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
##  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80 
##   0   0   0   0   0   0   0   0   0   0   1   0   0   0   0   0   0   1   0   0 
##  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 
##   0   0   0   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 
##   1   1   1   1   1   1   0   1   1   1   1   1   1   1   1   1   1   1   1   0 
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   0   1   1   1   1   1   1 
## 141 142 143 144 145 146 147 148 149 150 
##   1   1   1   1   1   1   1   1   1   1

Comentario: Si el caso se predice con un 0, la planta NO pertenence a la especie I. virginica; si se predice con un 1, la planta SÍ pertenece a esa especie.

Los primeros 100 casos NO provienen de la especie I. virginica. Sin embargo, hay tres ejemplares que erróneamente el modelo predice que sí lo son: los casos 71, 78 y 84. Mientras que hay otros tres casos en los que plantas efectivamente de la especie en cuestión se predicen como ajenas; los casos 107, 120 y 134.

4.1 Matriz de confusión

library(caret)
confusionMatrix(as.factor(predicciones), as.factor(iris_aae$virgi_vs_not))
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  0  1
##          0 97  3
##          1  3 47
##                                          
##                Accuracy : 0.96           
##                  95% CI : (0.915, 0.9852)
##     No Information Rate : 0.6667         
##     P-Value [Acc > NIR] : <2e-16         
##                                          
##                   Kappa : 0.91           
##                                          
##  Mcnemar's Test P-Value : 1              
##                                          
##             Sensitivity : 0.9700         
##             Specificity : 0.9400         
##          Pos Pred Value : 0.9700         
##          Neg Pred Value : 0.9400         
##              Prevalence : 0.6667         
##          Detection Rate : 0.6467         
##    Detection Prevalence : 0.6667         
##       Balanced Accuracy : 0.9550         
##                                          
##        'Positive' Class : 0              
##

5 Visualizaciones clave

5.1 Tamaño del sépalo

ggplot(iris_aae, aes(x = sepal_length, y = sepal_width,
                     color = variety, shape = variety)) +
  geom_point() + 
  scale_color_manual(values = c("orange", "limegreen", "darkred")) +
  labs(title = "Tamaño del sépalo en el género Iris",
       x = "Longitud del sépalo (cm)", y = "Anchura del sépalo (cm)") +
  theme(legend.position = c(.65, .85))

5.2 Tamaño del pétalo

ggplot(iris_aae, aes(x = petal_length, y = petal_width, color = variety,
                     shape = variety)) +
  geom_point() +
  scale_color_manual(values = c("orange", "limegreen", "darkred")) +
  labs(title = "Tamaño del pétalo en el género Iris",
       x = "Longitud del pétalo (cm)", y = "Anchura del pétalo (cm)") +
  theme(legend.position = c(.9, .2))

5.3 Tamaño del sépalo vs. tamaño del pétalo

ggplot(iris_aae, aes(x = petal_length, fill = variety)) +
  geom_density(alpha = 0.3) +
  labs(title = "Pétalos del género Iris", subtitle = "Longitud del pétalo (cm)",
       y = "Densidad", x = "") +
  theme(legend.position = c(.9, .8)) # longitud pétalo

ggplot(iris_aae, aes(x = petal_width, fill = variety)) +
  geom_density(alpha = 0.3) +
  labs(title = "Pétalos del género Iris", subtitle = "Anchura del pétalo (cm)",
       y = "Densidad", x = "") +
  theme(legend.position = c(.9, .8)) # anchura pétalo

ggplot(iris_aae, aes(x = sepal_length, fill = variety)) +
  geom_density(alpha = 0.3) +
  labs(title = "Sépalos del género Iris", subtitle = "Longitud del sépalo (cm)",
       y = "Densidad", x = "") +
  theme(legend.position = c(.9, .8)) # longitud del sépalo

ggplot(iris_aae, aes(x = sepal_width, fill = variety)) +
  geom_density(alpha = 0.3) +
  labs(title = "Sépalos del género Iris", subtitle = "Anchura del sépalo (cm)",
       y = "Densidad", x = "") +
  theme(legend.position = c(.9, .8)) # anchura del sépalo

6 Modelo final

6.1 Fórmula

\(I. virginica = -45.27 + 5.75(LongPetalo) + 10.40(AnchPetalo)\)

6.2 Interpretación de coeficientes

Como ambos coeficientes relacionados con las variables son positivos, tanto la longitud como la anchura aumentan las probabilidades de que la planta pertenezca a la especie I. virginica cuando sus valores son muy altos.

6.3 Significancia

Ambos coeficientes son significativos a más de 95 %. Sin embargo, la anchura del pétalo es más del doble de potente que la longitud del pétalo para separar a la especie I. virginica de las otras.

6.4 Métricas globales

La sensibilidad del modelo es de 97% – sólo aparecen tres plantas de 100 predichas como virginica pero que, en realidad, no lo son.

La especificidad alcanza 94% – sólo tres plantas que, siendo efectivamente virginica, no las reconoce como tales.

La exactitud se eleva a 0.96 – sólo seis de 150 casos se clasificaron erróneamente (intervalo de confianza a 95%: 0.915, 0.985)

6.5 Discusión del modelo:

Basado en las dimensiones de los pétalos, el modelo predice la probabilidad de que una planta provenga de la especie Iris virginica. Las medidas de los sépalos aportan poca información para hacer esta distinción.

Atendiendo a las probabilidades asociadas con la distribución de cada variable, y asociándolas con las de dos o más variables, se establece una regla de predicción probabilística.

Limitaciones: cuando las distribuciones de las variables no están superpuestas, el modelo logit no hace nada en realidad, por lo que se hace necesario realizar un análisis exploratorio exhaustivo que permita al analista vislumbrar las variables pertinentes para la clasificación.

Fortalezas: si las variables se eligen bien, con bases teóricas sólidas y un muestreo cuidadoso, el modelo de predicción logit puede proporcionar poderosas predicciones de una variable dicotómica basada en los valores de uno o varias variables cuantitativas o cualitativas.

7 Conclusiones

  • Las especies presentadas del género Iris pueden clasificarse razonablemente bien a partir de sus flores.
  • De las dos estructuras consideradas (sépalos y pétalos), son estos últimos los que permiten una mejor distinción.
  • Los pétalos de menor tamaño que la caracterizan pertenecen a la especie I. setosa.
  • El traslape de las distribuciones de las medidas de los pétalos entre I. versicorlor e I. virginica puede abordarse adecuadamente a través de una regresión logit, la cual resalta la importancia de los pétalos en la diferenciación de ambas especies.

Referencias

Anderson, Edgar. 1936. “The Species Problem in Iris.” Annals of the Missouri Botanical Garden 23 (3): 457–509. https://doi.org/10.2307/2394164.
Fisher, Ronald A. 1936. “The Use of Multiple Measurements in Taxonomic Problems.” Annals of Eugenics 7: 179–88. https://doi.org/doi:10.1111/j.1469-1809.1936.tb02137.x.
Kuhn, Max. 2008. “Building Predictive Models in r Using the Caret Package.” Journal of Statistical Software 28 (5): 1–26. https://doi.org/10.18637/jss.v028.i05.
Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Segunda edición. Use r! Nueva York: Springer-Verlag New York. link-springer-com.pbidi.unam.mx:2443/book/10.1007/978-3-319-24277-4.
Wickham, Hadley, Romain François, Lionel Henry, Kirill Müller, and Davis Vaughan. 2023. “Dplyr: A Grammar of Data Manipulation.” https://CRAN.R-project.org/package=dplyr.