Selección y regularización

JPAG1

1unam

Abstract

Ajustes alternativos que pueden producir una mejor precisión de predicción e interpretabilidad en la modelación.

Contents

1 Introducción

Los procedimientos de ajuste alternativos exploran los siguientes procedimientos:

• Precisión de la predicción : siempre que la relación entre la respuesta y sus predictores sea aproximadamente lineal, las estimaciones de mínimos cuadrados tendrán un sesgo bajo. Si \(n>>p\) , lo que significa que el número de observaciones \(n\) es mucho mayor que el número de predictores \(p\), entonces las estimaciones de mínimos cuadrados tienden a tener también baja varianza. A medida que \(p\) se acerca a \(n\), puede haber mucha variabilidad en el ajuste de mínimos cuadrados, lo que podría resultar en un sobreajuste y predicciones deficientes en las observaciones futuras. Si \(p>n\) , ya no hay una estimación única del coeficiente de mínimos cuadrados; el método no funciona. Al restringir o reducir los coeficientes estimados, se reduce significativamente la varianza a costa de un aumento insignificante del sesgo.

• Interpretabilidad del modelo : es común que las variables predictoras utilizadas en un modelo de regresión múltiple no estén asociadas con la respuesta. La inclusión de estas variables irrelevantes conduce a una complejidad innecesaria en el modelo resultante. Si se pueden eliminar estas variables estableciendo sus coeficientes iguales a cero, se puede obtener un modelo más simple e interpretable. La probabilidad de que los mínimos cuadrados produzcan un coeficiente cero es bastante baja.

Hay tres clases importantes de métodos:

1. Selección de subconjuntos. Este enfoque implica identificar un subconjunto de los predictores que se cree están relacionados con la respuesta.

2. Contracción. Este enfoque implica ajustar un modelo de todos los predictores \(p\), pero reducir (también conocido como regularizar) los coeficientes de algunos predictores hacia cero. Esto también puede resultar en una selección de variables cuando los coeficientes se reducen exactamente a cero.

3. Reducción de dimensión . Este enfoque implica proyectar los \(p\) predictores en un subespacio \(M\)-dimensional, donde \(M<p\). Esto se logra calculando \(M\) diferentes combinaciones lineales , o proyecciones , de las variables. Luego, se usan estas proyecciones como predictores.

2 Selección de subconjunto

2.1 Selección del mejor subconjunto

Sa ajusta una regresión de mínimos cuadrados separada para cada combinación posible de los \(p\) predictores . Es decir, se ajustan todos los \(p\) modelos que contienen exactamente un predictor, todos los \(\binom{p}{2}\) que contienen exactamente dos predictores, y así sucesivamente. Una vez ajustados todos, se identifica el que es mejor.

Estos son los pasos:

1. Sea \(M_0\) el modelo nulo sin predictores. Esta es simplemente la media muestral.

2. Para \(k = 1,2,...p\):

1. Ajustar todos los \(\binom{p}{k}\) modelos que contienen exactamente \(k\) predictores.

2. Elegir el mejor entre estos \(\binom{p}{k}\) modelos mediante la\(R^2\) más grande.

3. Seleccionar un único mejor modelo de \(M_0,...M_p\) usando validación cruzada, \(C_p\), AIC, BIC o \(R^2\) ajustada.

Una vez que completado el paso \(2\), se reduce el problema de uno de \(2^p\) modelos posibles a uno de \(p+1\) modelos posibles. Para completar el paso \(3\), ya no se puede usar una métrica como \(R^2\), ya que aumenta monótonamente a medida que aumenta el número de características incluidas en la modelación. Por lo tanto, se debe elegir el modelo con el error de prueba estimado más bajo.

Figure 1: R² vs predictors

Como se dijo anteriormente, \(R^2\) siempre aumentará a medida que se agregan más predictores. En este caso, aumenta a través de tres predictores antes de aplanarse. Se puede aplicar la misma idea a otro tipo de modelos como la regresión logística. En lugar de ordenar por la \(R^2\), se podría ordenar por la devianza , la medida que juega el papel de \(R^2\) para una clase más amplia de modelos; cuanto menor sea la devianza, mejor será el ajuste.

Un problema con la mejor selección de subconjuntos es el costo computacional. Ajustar \(2^p\) modelos posibles rápidamente se vuelve prohibitivamente costoso.

Figure 2: Time vs predictors

Ajustar la selección del mejor subconjunto para los \(20\) predictores de la gráfica anterior tomaría menos de diez segundos, lo cual es razonable. A los \(40\) predictores, ese número es de unos \(100,000\) segundos, más de un día completo. Dado que es común que un conjunto de datos contenga cientos, si no miles, de predictores, se necesitan explorar alternativas más eficientes desde el punto de vista computacional.

3 Selección escalonada

Debido a los costos computacionales, la mejor selección de subconjuntos no se puede aplicar con una \(p\) muy grande . También hay problemas estadísticos cuando \(p\) es grande. Cuanto mayor sea el espacio de búsqueda, mayor será la probabilidad de encontrar modelos que se vean bien en los datos de entrenamiento. Un espacio de búsqueda enorme puede provocar un sobreajuste y una gran variación de las estimaciones de los coeficientes.

Por ambas razones, los métodos escalonados , que exploran un conjunto más restringido de modelos, son alternativas atractivas.

3.1 Selección escalonada hacia adelante

La selección escalonada hacia adelante es una alternativa computacionalmente eficiente a la mejor selección de subconjuntos. Comienza con un modelo que no contiene predictores y luego agrega predictores uno a la vez, hasta que todos los predictores están en el modelo. En cada paso, se agrega al modelo la variable que brinda la mayor mejora adicional al ajuste.

Estos son los pasos:

1. Selección progresiva hacia adelante La selección progresiva hacia adelante es una alternativa computacionalmente eficiente a la mejor selección de subconjuntos. Comienza con un modelo que no contiene predictores y luego agrega predictores uno a la vez, hasta que todos los predictores están en el modelo. En cada paso, se agrega al modelo la variable que brinda la mayor mejora adicional al ajuste.

Estos son los pasos:

1. Sea \(M_0\) el modelo nulo sin predictores. Esta es simplemente la media muestral.

2. Para \(k=0,1,..,p-1\):

1. Considerara todos los \(p-k\) modelos que aumentan los predictores en \(M_k\) con un predictor adicional.

2. Elegir el mejor entre estos \(p-k\) modelos y llamarlo \(M_{k-1}\). Mejor se define como el que tiene la \(R^2\) más alta.

3. Seleccionar un único mejor modelo de \(M_0,...M_p\) usando validación cruzada, \(C_p\), AIC, BIC o \(R^2\) ajustada.

La selección escalonada hacia adelante se ajusta a los \(1+p(p+1)/2\) modelos , que se pueden aproximar a \(p^2\). Computacionalmente, esto escala mucho mejor a medida que \(p\) crece.

Figure 3: Forward stepwise vs best subset selection

Elegir el mejor modelo de los que aumentan \(M_k\) con un predictor adicional es simple. Simplemente se puede usar el valor de \(R²\). El paso \(3\), en el que se evalúa qué modelo entre los modelos seleccionados es el mejor, requiere usar un método que pueda estimar el error de prueba.

La desventaja de ser computacionalmente económico en comparación con la selección del mejor subconjunto es que no siempre encontrará el mejor modelo entre los \(2^p\) modelos posibles. Suponiendo un conjunto de datos que contiene \(p=3\) predictores. El mejor modelo posible de una variable contiene a \(X_1\) , y el mejor modelo posible de dos variables contiene a \(X_2\) y \(X_3\). La selección escalonada hacia adelante nunca producirá el modelo ideal de dos variables, ya que siempre se ajustará \(X_1\) al modelo de una variable. Por lo tanto, \(M_2\) siempre contendrá a \(X_1\).

3.2 Selección escalonada hacia atrás

La selección por pasos hacia atrás es similar a la selección por pasos hacia adelante, pero comienza con el modelo de mínimos cuadrados completo que contiene todos los predictores \(p\). Luego, elimina iterativamente el predictor menos útil, uno a la vez.

Estos son los pasos:

1. Sea \(M_0\) el modelo completo de todos los predictores.

2. Para \(k=p,p-1,...1\):

1. Considerar todos los \(k\) modelos que contienen todos menos uno de los predictores en \(M_k\), para un total de \(k-1\) predictores.

2. Elegir el mejor entre estos \(k\) modelos y llamrlo \(M_{k-1}\). Mejor se define como el que tiene la \(R^2\) más alta.

3. Seleccionar un único mejor modelo de \(M_0,...M_p\) usando validación cruzada, \(C_p\), AIC, BIC o \(R^2\) ajustada.

Al igual que en la selección hacia adelante, también se busca a través de aproximadamente \(p²\) modelos. Además, al igual que la selección hacia adelante, no se garantiza que se produzca el mejor modelo que contenga un subconjunto de los \(p\) predictores . La selección hacia atrás requiere que \(n\) sea mayor que \(p\) para ajustarse al modelo completo.

3.3 Enfoques híbridos

Otra alternativa es un enfoque híbrido. Las variables se pueden agregar al modelo secuencialmente, como en la selección hacia adelante. Sin embargo, después de agregar cada nueva variable, el método también puede eliminar cualquier variable que ya no proporcione una mejora en el ajuste del modelo. Tal enfoque intenta imitar la mejor selección de subconjuntos mientras que conserva las ventajas computacionales de la selección escalonada.

4 Elección del modelo óptimo

4.1 Eatadística de Mallow, AIC, BIC y R cuadrada ajustada.

En general, el MSE del conjunto de entrenamiento es una subestimación del MSE del conjunto de prueba . Cuando se ajusta un modelo a los datos de entrenamiento usando mínimos cuadrados, se estiman específicamente los coeficientes de regresión de manera que el RSS de entrenamiento sea lo más pequeño posible. El error de entrenamiento siempre disminuirá a medida que agreguemos más variables al modelo, pero es posible que el error de prueba no lo haga. Por lo tanto, no se pueden usar métricas como la \(R^2\) para seleccionar modelos que contienen diferentes números de variables.

Se tienen varias técnicas para ajustar el error de entrenamiento. Se consideran cuatro de estos enfoques: \(C_p\), criterio de información de Akaike (AIC) , criterio de información bayesiano (BIC) y \(R^2\) ajustada.

Para un modelo de mínimos cuadrados ajustado que contiene \(d\) predictores , la estimación de \(C_p\) del MSE de pruebda se calcula usando la ecuación:

\[\begin{equation} C_p = \frac{1}{n}\Big(RSS+2d\widetilde{\sigma}^2\Big), \tag{1} \end{equation}\]

donde \(\widetilde{\sigma}^2\) es una estimación de la varianza del error \(\varepsilon\) asociado con cada medición de respuesta. Normalmente, \(\widetilde{\sigma}^2\) se estima utilizando el modelo completo con todos los predictores. La estadística \(C_p\) esencialmente agrega una penalización de \(2d\widetilde{\sigma}^2\) al RSS de entrenamiento para ajustar el hecho de que el error de entrenamiento tiende a subestimar el error de prueba. Cuantos más predictores, mayor será la penalización.

El criterio AIC se define para una gran clase de modelos ajustados por máxima verosimilitud. En el caso del modelo con errores gaussianos, máxima verosimilitud y mínimos cuadrados son lo mismo.

\[\begin{equation} AIC = \frac{1}{n\widetilde{\sigma}^2}\Big(RSS+2d\widetilde{\sigma}^2\Big), \tag{2} \end{equation}\]

Para mínimos cuadrados, AIC y \(C_p\) son proporcionales entre sí.

BIC parece similar pero puede imponer una penalización más alta en modelos con muchas variables.

\[\begin{equation} BIC = \frac{1}{n\widetilde{\sigma}^2}\Big(RSS+\log(n) \widetilde{\sigma}^2\Big), \tag{3} \end{equation}\]

BIC reemplaza \(3D\widetilde{\sigma}^2\) usado por \(C_p\) con \(\log(n)\). Dado que \(\log(n)>2\) para cualquier \(n>7\), BIC impone una penalización mayor a los modelos con muchas variables y tiende a seleccionar modelos más pequeños que \(C_p\).

La \(R^2\) ajustada es otro enfoque para seleccionar entre un conjunto de modelos. La \(R^2\) se define como \(1-RSS/TSS\), dado que RSS sólo puede disminuir a medida que se agregan más variables al modelo, la \(R^2\) siempre aumenta a medida que se agregan más variables. Para un modelo de mínimos cuadrados con \(d\) variables, la \(R^2\) ajustada se calcula como:

\[\begin{equation} R^2_{adj}=1−\Big[\frac{(1−R^2)(n−1)}{n−d−1}\Big] \tag{4} \end{equation}\]

Maximizar la \(R^2\) ajustada es equivalente a minimizar el numerador, \(RSS/(n-d-1)\). La idea detrás de la \(R^2\) ajustada es que las variables de ruido adicionales aumentarán pero también aumentará \(d\), lo que en última instancia puede aumentar \(RSS/(n-d-1)\). La \(R^2\) ajustada paga un precio por la inclusión de variables innecesarias en el modelo.

\(C_p\), AIC y BIC tienen justificaciones teóricas rigurosasy además, AIC y BIC también se pueden definir para tipos de modelos más generales.

4.2 Validación cruzada

Se puede estimar directamente el error de prueba utilizando el conjunto de validación y los métodos de validación cruzada.

Este procedimiento tiene una ventaja sobre AIC, BIC y \(C_p\) en que proporciona una estimación directa del error de prueba y hace menos suposiciones sobre el modelo subyacente. También se puede utilizar en casos en los que sea difícil estimar \({\sigma}^2\) y / o no se conozca el número de grados de libertad.

La validación cruzada se ha convertido en un enfoque más atractivo a medida que ha aumentado la potencia computacional.

5 Métodos de reducción

Los métodos de subconjunto anteriores implican ajustar un modelo de mínimos cuadrados que contiene un subconjunto de predictores. También se puede ajustar un modelo que contenga todos los \(p\) predictores que restringen o regularizan las estimaciones de los coeficientes. Reducir las estimaciones de los coeficientes puede reducir significativamente su varianza. Las dos técnicas más conocidas son la regresión de ridge y la de lasso.

5.1 Regresión de ridge

La regresión de ridge es muy similar a los mínimos cuadrados, excepto que estamos minimizando los coeficientes en una cantidad diferente.

\[\begin{equation} \sum_{i=1}^n \Big(y_i-{\beta}_0-\sum_{j=1}^p{\beta}_j x_{ij}\Big)^2 +{\lambda}\sum_{j=1}^p {\beta}_j^2 = RSS + {\lambda}\sum_{j=1}^p{\beta}_j^2 \tag{5} \end{equation}\]

Donde \({\lambda}\) es un parámetro de ajuste o hiperparámetro que se determinará por separado. De manera similar a los mínimos cuadrados, se buscan estimaciones de coeficientes que reduzcan el valor de RSS. Sin embargo, se agrega un segundo término, \({\lambda}\sum_{j=1}^p{\beta}_j^2\), llamado penalización por contracción . Esta penalización por contracción es pequeña cuando \({\beta}_1,...{\beta}_p\) están cerca de cero y, por lo tanto, tiene el efecto de reducir las estimaciones de \({\beta}_j\) hacia cero. La penalización, \({\beta}_1^2\), también se conoce como penalización L2. \({\lambda}\) controla el impacto relativo de estos dos términos (en \(λ=0\), esta es la estimación de mínimos cuadrados original). La regresión de ridge producirá un conjunto diferente de estimaciones de coeficientes para cada valor de \({\lambda}\). Seleccionar un buen valor de \({\lambda}\) es fundamental.

Se observe que la contracción no se aplica a la intersección \(β_0\). Si las variables se han centrado para tener una media de cero antes de realizar la regresión ridge, entonces la intersección estimada tomará el valor de la media de la muestra.

5.1.1 ¿Por qué mejora la regresión de ridge con respecto a los mínimos cuadrados?

Se vuelve nuevamente a la compensación de sesgo-varianza . A medida que \(λ\) aumenta y los coeficientes se reducen, la flexibilidad del ajuste disminuye, lo que reduce la varianza pero aumenta el sesgo. Si \(λ\) se vuelve demasiado grande y los coeficientes se reducen aún más, tienden a subestimarse y la reducción de la varianza se ralentiza mientras que el aumento del sesgo se acelera.

Cuando \(n\) está cerca de \(p\), la estimación de mínimos cuadrados será extremadamente variable. La regresión de ridge funciona mejor cuando las estimaciones de mínimos cuadrados tienen una alta varianza.

Computacionalmente, la regresión ridge es mucho más rápida que la mejor selección de subconjuntos, que busca a través de \(2^p\) modelos . Por el contrario, la regresión de ridge se ajusta a un solo modelo para cada valor de \({\lambda}\) que verifica.

5.2 El lasso

Una desventaja de la regresión de ridge es que nunca reduce completamente los coeficientes hacia cero. Esto puede no ser un problema para la precisión de la predicción, pero ciertamente podría afectar la interpretabilidad del modelo.

Si se quiere construir un modelo disperso que involucre menos predictores, entonces se necesita explorar un método de contracción diferente.

El lasso es una alternativa a la regresión ridge. Los coeficientes del lasso, \(\widehat{\beta}_j^L\) , minimizan la cantidad:

\[\begin{equation} \sum_{i=1}^n \Big(y_i-{\beta}_0-\sum_{j=1}^p{\beta}_j x_{ij}\Big)^2 +{\lambda}\sum_{j=1}^p |{\beta}_j| = RSS + {\lambda}\sum_{j=1}^p{\beta}_j^2 \tag{6} \end{equation}\]

La única diferencia entre esto y la regresión de ridge es el coeficiente de penalización. El lasso usa una penalización L1 en lugar de una penalización L2.

El lasso también reduce los coeficientes hacia cero. Pero debido a la penalización L1, el lazo puede forzar que algunos coeficientes sean exactamente cero cuando \(λ\) es suficientemente grande. Por tanto, el lasso tiene la capacidad de realizar una selección de variables , produciendo modelos dispersos que involucran un subconjunto de los predictores originales.

Otra forma de pensar en las penalizaciones es imaginarlas como una restricción o presupuesto. Se puede pensar en la regresión ridge y el lasso como un ajuste de mínimos cuadrados con un componente de restricción adicional , definido por la regularización L1 / L2.

6 Estudio de caso: subtipo de enfermedad a partir de datos genómicos

Se empleará un conjunto de datos reales de biopsias de tumores. Se usarán los datos de expresión génica de muestras de tumores de glioblastoma del proyecto (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga), The Cancer Genome Atlas. Se intentará predecir el subtipo de esta enfermedad mediante marcadores moleculares. Este subtipo se caracteriza por alteraciones epigenéticas a gran escala llamadas fenotipo metilador de isla CpG o CIMP; la mitad de la base de datos de pacientes tienen este subtipo y el resto no, y se trata de predecir cuáles tienen el subtipo CIMP. La base de datos se organiza con pacientes en cada fila y cada columna es un valor de expresión génica. Hay \(184\) muestras de tumores. La información se encuentra en el objeto tgexp3, y se inspecciona con la función glimpse().

glimpse(tgexp3[,1:5])

## Rows: 184
## Columns: 5
## $ subtype <fct> noCIMP, CIMP, CIMP, CIMP, noCIMP, CIMP, noCIMP, noCIMP, noCIMP…
## $ XIRP2   <dbl> 0.7110307, 0.5467647, 0.0000000, 1.3373135, 1.0695719, 0.74831…
## $ RPS4X   <dbl> 3.851703, 4.354194, 4.309313, 4.113485, 3.911905, 4.013907, 4.…
## $ DDX3X   <dbl> 3.864235, 3.700598, 3.777171, 3.856144, 3.849256, 3.841601, 3.…
## $ KDM5C   <dbl> 3.335923, 3.175486, 3.356665, 3.387141, 3.307022, 3.509791, 3.…

En este punto, se opta por dividir los datos en las particiones de prueba y de entrenamiento. La razón de esto es que es necesaria una prueba independiente en la que no se entrena, pero sin tener un conjunto de prueba por separado, no se puede evaluar el rendimiento de cierta modelo o ajustarla correctamente.

Para empezar, se dividirá al \(30\%\) de los datos como prueba. Este método es el estándar de oro para probar el rendimiento de un modelo. Al hacer esto, se tiene un conjunto de datos separado que el modelo nunca ha visto.

set.seed(3031)
tgexp_split <- initial_split(tgexp3, prop = 0.7, strata = subtype)
tgexp_train <- training(tgexp_split)
tgexp_test <- testing(tgexp_split)
glimpse(tgexp_train[,1:5]) 

## Rows: 128
## Columns: 5
## $ subtype <fct> CIMP, CIMP, CIMP, CIMP, CIMP, CIMP, CIMP, CIMP, CIMP, CIMP, CI…
## $ XIRP2   <dbl> 0.5467647, 0.0000000, 1.3373135, 0.7483121, 0.5423897, 0.70096…
## $ RPS4X   <dbl> 4.354194, 4.309313, 4.113485, 4.013907, 4.206883, 4.870087, 4.…
## $ DDX3X   <dbl> 3.700598, 3.777171, 3.856144, 3.841601, 3.617547, 3.533291, 3.…
## $ KDM5C   <dbl> 3.175486, 3.356665, 3.387141, 3.509791, 3.285774, 3.071785, 3.…

glimpse(tgexp_test[,1:5])

## Rows: 56
## Columns: 5
## $ subtype <fct> noCIMP, noCIMP, CIMP, CIMP, CIMP, noCIMP, noCIMP, noCIMP, noCI…
## $ XIRP2   <dbl> 0.3699761, 0.0000000, 0.8625845, 0.1862215, 0.3607638, 0.49096…
## $ RPS4X   <dbl> 4.118965, 3.794960, 4.016565, 3.961572, 4.195175, 3.845546, 4.…
## $ DDX3X   <dbl> 3.853480, 3.673748, 3.664307, 3.789902, 3.533150, 3.951146, 3.…
## $ KDM5C   <dbl> 3.414320, 3.399652, 3.244334, 3.420281, 3.220843, 3.376112, 3.…

Paso siguiente, se tendrán que pre-procesar los datos antes de comenzar a entrenar. Esto podría incluir un análisis de datos exploratorio para ver cómo las variables y las muestras se relacionan entre sí. Por ejemplo, es posible verificar la correlación entre las variables predictoras y mantener solo una variable de ese grupo. Además, algunos algoritmos de entrenamiento pueden ser sensibles a escalas de datos o valores atípicos. Esos problemas se llevan a cabo en este paso. En algunos casos, los datos pueden tener valores perdidos, se puede optar por eliminar las muestras a las que les faltan valores o intentar imputarlas ya que muchas algorítmicas de aprendizaje automático no podrán lidiar con los valores perdidos.

En genómica, es posible que no se tengan valores para ciertos genes o ubicaciones genómicas debido a problemas técnicos durante los experimentos. check_missing() crea una especificación de recipe que verificará si las variables contienen valores perdidos.

recipe(tgexp_train) %>%
  check_missing(all_predictors()) %>%
  prep()

## Data Recipe
## 
## Inputs:
## 
##   51 variables (no declared roles)
## 
## Training data contained 128 data points and no missing data.
## 
## Operations:
## 
## Check missing values for <none> [trained]

También se filtrarán las variables predictoras que tienen poca variación. No es probable que tengan importancia predictiva, ya que no hay mucha variación y solo ralentizarán la modelación. Cuantas más variables, más lentos serán los algoritmos en general. step_nzv() se encarga de esto. Centrar y escalar los datos se utiliza generalmente para mejorar la estabilidad numérica de algunos cálculos. Las variables predictoras que están altamente correlacionadas también se filtran, y step_corr() lo hará dentro de la pipeline, especificando el argumento threshold. Dicho esto se crea la siguiente recipe:

tgexp_recipe2 <- recipe(subtype ~ ., data = tgexp_train) %>%
  step_nzv(all_predictors()) %>%
  step_center(all_predictors()) %>%
  step_corr(all_numeric_predictors(), threshold = .9)

Ahora se aplica la regresión logística, pues es un método estadístico que se utiliza para modelar una variable de respuesta binaria basada en variables predictoras,, ya que la variable obejtivo CMIP de muestras de tumores es una respuesta binaria o un problema de dos clases. La regresión logística se puede considerar como un problema de estimación de máxima verosimilitud en el que se tratan de encontrar parámetros estadísticos que maximicen la probabilidad de que los datos observados sean muestreados a partir de la distribución estadística de interés. Esto también está muy relacionado con el enfoque de la función de costo / pérdida general en los algoritmos de aprendizaje automático supervisados. En el caso de las variables de respuesta binaria, el modelo de regresión lineal simple, sería una mala elección porque puede generar fácilmente valores fuera del límite de \(0\) a \(1\).

En la regresión logística, la variable de respuesta se modela con una distribución binomial o bernoulli. El valor de cada variable de respuesta es \(0\) o \(1\), y se necesitan averiguar los valores del parámetro que podrían generar tal distribución de \(0\) y \(1\). Si se encuentran los mejores valores para cada muestra de tumor, se estará maximizando la función de log-verosimilitud del modelo sobre los datos observados.

Entonces, se usará una variable predictora, la expresión de un gen para clasificar las muestras de tumores en los subtipos CIMP y noCIMP usando la interfaz de metapaqueterías en tidymodels, donde se proporcionan los nombres de la respuesta y las variables predictoras en un workflow.

set.seed(3031)

log_reg <- logistic_reg() %>%
  set_engine("glm") %>% 
  set_mode("classification")

log_fit3 <- workflow() %>%
  add_recipe(tgexp_recipe2) %>%
  add_model(log_reg) %>%
  fit(data = tgexp_train)

## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

log_fit3

## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: logistic_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 3 Recipe Steps
## 
## • step_nzv()
## • step_center()
## • step_corr()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## 
## Call:  stats::glm(formula = ..y ~ ., family = stats::binomial, data = data)
## 
## Coefficients:
## (Intercept)        XIRP2        RPS4X        DDX3X        KDM5C        USP9X  
##      4.2548     -15.9294     -57.6326      77.2900      40.1340    -129.6705  
##      EIF1AX        LTC4S          POR         UTS2       TTTY14        GSTK1  
##     51.6012      -1.9226     -67.1621      36.1571       9.8250     101.2314  
##        CHGB       ZFYVE9        GSTP1        CYYR1        TSHZ3        HOXA6  
##      4.5469      24.4609     -35.8569      40.1073     -26.6271      12.3161  
##       MEOX1         OPA3      SLC17A6         LGR5        SHOC2        MOV10  
##     -1.0115      27.5331      -2.5232     -14.9773      13.5160      60.1655  
##        WIBG         NEFH         GZMM         IRX1        SFRP1         IRS4  
##     37.3805       8.0981       5.8313      -0.8489       8.7699       4.1113  
##      HOTAIR        HOXB9         GJB5        GRIK5         PRKY         PDXP  
##     13.1294     -11.5801     -19.6575     -10.3203       2.2611      13.8764  
##       KLRC1      TMSB15B       NLGN4X       DMRTA1        ABCC2       GABPB2  
##      6.3006     -25.2256      19.6514      21.1458     -20.9282      -7.2582  
##      GABRG1       SLC6A9        IGFL4        HOXB2        HOXA3       MOBKL3  
##     -2.2053      37.6044     -29.1260      -0.5449       2.7801     -12.1784  
##       MAGT1  
##    -66.5697  
## 
## Degrees of Freedom: 127 Total (i.e. Null);  79 Residual
## Null Deviance:       177.4 
## Residual Deviance: 2.219e-09     AIC: 98

Ahora se predicen las probabilidades para los datos de prueba.

tgexp_test_predictions2 <- predict(log_fit3, tgexp_test) %>%
  bind_cols(tgexp_test)
tgexp_test_predictions2

## # A tibble: 56 x 52
##    .pred_class subtype XIRP2 RPS4X DDX3X KDM5C USP9X EIF1AX LTC4S   POR  UTS2
##    <fct>       <fct>   <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
##  1 CIMP        noCIMP  0.370  4.12  3.85  3.41  3.71   3.41 1.32   3.34 0    
##  2 noCIMP      noCIMP  0      3.79  3.67  3.40  3.60   3.38 0.572  3.08 0.281
##  3 CIMP        CIMP    0.863  4.02  3.66  3.24  3.61   3.25 1.67   3.17 0    
##  4 CIMP        CIMP    0.186  3.96  3.79  3.42  3.55   3.32 1.42   3.15 0    
##  5 CIMP        CIMP    0.361  4.20  3.53  3.22  3.34   3.24 1.85   3.25 0    
##  6 noCIMP      noCIMP  0.491  3.85  3.95  3.38  3.51   3.29 1.67   3.33 0    
##  7 noCIMP      noCIMP  0.624  4.09  3.69  3.36  3.45   3.10 1.78   3.41 0.216
##  8 noCIMP      noCIMP  0.543  4.03  3.94  3.52  3.58   3.34 1.48   3.30 0    
##  9 noCIMP      noCIMP  0.725  3.99  3.95  3.51  3.66   3.34 1.35   3.20 0    
## 10 CIMP        CIMP    0.823  4.25  3.75  3.35  3.59   3.38 1.22   3.30 0.222
## # … with 46 more rows, and 41 more variables: TTTY14 <dbl>, GSTK1 <dbl>,
## #   CHGB <dbl>, ZFYVE9 <dbl>, GSTP1 <dbl>, CYYR1 <dbl>, TSHZ3 <dbl>,
## #   HOXA6 <dbl>, MEOX1 <dbl>, OPA3 <dbl>, CYorf15A <dbl>, SLC17A6 <dbl>,
## #   LGR5 <dbl>, SHOC2 <dbl>, MOV10 <dbl>, WIBG <dbl>, NEFH <dbl>, GZMM <dbl>,
## #   IRX1 <dbl>, DAOA <dbl>, SFRP1 <dbl>, IRS4 <dbl>, HOTAIR <dbl>, HOXB9 <dbl>,
## #   GJB5 <dbl>, GRIK5 <dbl>, PRKY <dbl>, PDXP <dbl>, KLRC1 <dbl>,
## #   TMSB15B <dbl>, NLGN4X <dbl>, DMRTA1 <dbl>, ABCC2 <dbl>, GABPB2 <dbl>,
## #   GABRG1 <dbl>, SLC6A9 <dbl>, IGFL4 <dbl>, HOXB2 <dbl>, HOXA3 <dbl>,
## #   MOBKL3 <dbl>, MAGT1 <dbl>

Con estas predicciones se evalúan las metricas de rendimiento

tgexp_test_predictions2 %>%
  metrics(truth = subtype, estimate = .pred_class)

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.821
## 2 kap      binary         0.643

Y se observa la matriz de confusión.

tgexp_test_predictions2 %>%
  conf_mat(truth = subtype, estimate = .pred_class)

##           Truth
## Prediction CIMP noCIMP
##     CIMP     22      4
##     noCIMP    6     24

La precisión de la prueba puede mejorar, por lo que puede resulta útil una técnica basada en remuestreo, como la validación cruzada.

tgexp_vfold <- vfold_cv(tgexp_train, v = 10, strata = subtype)
tgexp_vfold

## #  10-fold cross-validation using stratification 
## # A tibble: 10 x 2
##    splits           id    
##    <list>           <chr> 
##  1 <split [114/14]> Fold01
##  2 <split [114/14]> Fold02
##  3 <split [114/14]> Fold03
##  4 <split [114/14]> Fold04
##  5 <split [116/12]> Fold05
##  6 <split [116/12]> Fold06
##  7 <split [116/12]> Fold07
##  8 <split [116/12]> Fold08
##  9 <split [116/12]> Fold09
## 10 <split [116/12]> Fold10

Y desarrollando nuevamente un workflow para este remuestreo, se debe emplear fit_resamples() para ajustar.

set.seed(3031)

log_fit4 <- workflow() %>%
  add_recipe(tgexp_recipe2) %>%
  add_model(log_reg) %>%
  fit_resamples(resamples = tgexp_vfold)

log_fit4

## Warning: This tuning result has notes. Example notes on model fitting include:
## preprocessor 1/1, model 1/1: glm.fit: algorithm did not converge, glm.fit: fitted probabilities numerically 0 or 1 occurred
## preprocessor 1/1, model 1/1: glm.fit: algorithm did not converge, glm.fit: fitted probabilities numerically 0 or 1 occurred
## preprocessor 1/1, model 1/1: glm.fit: algorithm did not converge, glm.fit: fitted probabilities numerically 0 or 1 occurred

## # Resampling results
## # 10-fold cross-validation using stratification 
## # A tibble: 10 x 4
##    splits           id     .metrics         .notes          
##    <list>           <chr>  <list>           <list>          
##  1 <split [114/14]> Fold01 <tibble [2 × 4]> <tibble [1 × 1]>
##  2 <split [114/14]> Fold02 <tibble [2 × 4]> <tibble [1 × 1]>
##  3 <split [114/14]> Fold03 <tibble [2 × 4]> <tibble [1 × 1]>
##  4 <split [114/14]> Fold04 <tibble [2 × 4]> <tibble [1 × 1]>
##  5 <split [116/12]> Fold05 <tibble [2 × 4]> <tibble [1 × 1]>
##  6 <split [116/12]> Fold06 <tibble [2 × 4]> <tibble [1 × 1]>
##  7 <split [116/12]> Fold07 <tibble [2 × 4]> <tibble [1 × 1]>
##  8 <split [116/12]> Fold08 <tibble [2 × 4]> <tibble [1 × 1]>
##  9 <split [116/12]> Fold09 <tibble [2 × 4]> <tibble [1 × 1]>
## 10 <split [116/12]> Fold10 <tibble [2 × 4]> <tibble [1 × 1]>

Y se observan las nuevas métricas de rendimiento.

log_fit4 %>% collect_metrics()

## # A tibble: 2 x 6
##   .metric  .estimator  mean     n std_err .config             
##   <chr>    <chr>      <dbl> <int>   <dbl> <chr>               
## 1 accuracy binary     0.861    10  0.0295 Preprocessor1_Model1
## 2 roc_auc  binary     0.907    10  0.0234 Preprocessor1_Model1

El objeto workflow completo quedaría resumido en la siguiente pipeline.

workflow() %>%
  add_recipe(tgexp_recipe2) %>%
  add_model(log_reg) %>%
  fit_resamples(resamples = tgexp_vfold) %>%
  collect_metrics()

## # A tibble: 2 x 6
##   .metric  .estimator  mean     n std_err .config             
##   <chr>    <chr>      <dbl> <int>   <dbl> <chr>               
## 1 accuracy binary     0.861    10  0.0295 Preprocessor1_Model1
## 2 roc_auc  binary     0.907    10  0.0234 Preprocessor1_Model1

6.1 Regresión ridge

Ahora se analizarán los modelos de regularización y el ajuste de hiperparámetros en la penalización. Se usará la librería glmnet como engine para realizar la regresión de ridge pues parsnip no tiene una función dedicada para crear una especificación en la penalización a través de glm. Se debe usar logistic_reg() y configurar mixture = 0 para especificar un modelo de ridge. El argumento mixture especifica la cantidad de diferentes tipos de regularización, mixture = 0 especifica solo la regularización de ridge y mixture = 1 especifica solo la regularización de lasso. Establecer mixture un valor entre \(0\) y \(1\) permite usar ambas. Al usar el engine de glmnet, también se vuelve necesario configurar el argumento penalty para poder ajustar el modelo.

Se usará la función validation_split() para asignar el \(20\) por ciento de los diagnósticos en tgexp_train al conjunto de validación y el resto al conjunto de entrenamiento.

set.seed(3031)

val_set <- validation_split(tgexp_train, 
                            strata = subtype, 
                            prop = 0.80)
val_set

## # Validation Set Split (0.8/0.2)  using stratification 
## # A tibble: 1 x 2
##   splits           id        
##   <list>           <chr>     
## 1 <split [102/26]> validation

Así, se especifica mixture = 0.

lr_mod <- 
  logistic_reg(penalty = tune(), mixture = 0) %>% 
  set_engine("glmnet")

lr_workflow <- 
  workflow() %>% 
  add_model(lr_mod) %>% 
  add_recipe(tgexp_recipe2)

Se crea la cuadrícula para alimentar a la función de penalización y se pasa al objeto tipo workflow.

lr_reg_grid <- tibble(penalty = 10^seq(-4, -1, length.out = 30))

lr_res <- 
  lr_workflow %>% 
  tune_grid(val_set,
            grid = lr_reg_grid,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))
lr_res 

## # Tuning results
## # Validation Set Split (0.8/0.2)  using stratification 
## # A tibble: 1 x 5
##   splits           id         .metrics         .notes          .predictions     
##   <list>           <chr>      <list>           <list>          <list>           
## 1 <split [102/26]> validation <tibble [30 × 5… <tibble [0 × 1… <tibble [780 × 6…

Y por último se crean los objetos necesarios que guardan las métricas de rendimiento y permiten visualizarlas.

lr_plot <- 
  lr_res %>% 
  collect_metrics() %>% 
  ggplot(aes(x = penalty, y = mean)) + 
  geom_point() + 
  geom_line() + 
  ylab("Area under the ROC curve") +
  scale_x_log10(labels = scales::label_number())
lr_plot 

lr_best <- 
  lr_res %>% 
  collect_metrics() %>% 
  arrange(penalty)
lr_best

## # A tibble: 30 x 7
##     penalty .metric .estimator  mean     n std_err .config              
##       <dbl> <chr>   <chr>      <dbl> <int>   <dbl> <chr>                
##  1 0.0001   roc_auc binary     0.988     1      NA Preprocessor1_Model01
##  2 0.000127 roc_auc binary     0.988     1      NA Preprocessor1_Model02
##  3 0.000161 roc_auc binary     0.988     1      NA Preprocessor1_Model03
##  4 0.000204 roc_auc binary     0.988     1      NA Preprocessor1_Model04
##  5 0.000259 roc_auc binary     0.988     1      NA Preprocessor1_Model05
##  6 0.000329 roc_auc binary     0.988     1      NA Preprocessor1_Model06
##  7 0.000418 roc_auc binary     0.988     1      NA Preprocessor1_Model07
##  8 0.000530 roc_auc binary     0.988     1      NA Preprocessor1_Model08
##  9 0.000672 roc_auc binary     0.988     1      NA Preprocessor1_Model09
## 10 0.000853 roc_auc binary     0.988     1      NA Preprocessor1_Model10
## # … with 20 more rows

lr_auc <- 
  lr_res %>% 
  collect_predictions(parameters = lr_best) 
lr_auc

## # A tibble: 780 x 7
##    id         .pred_CIMP .pred_noCIMP  .row penalty subtype .config             
##    <chr>           <dbl>        <dbl> <int>   <dbl> <fct>   <chr>               
##  1 validation      0.597      0.403       1  0.0001 CIMP    Preprocessor1_Model…
##  2 validation      0.984      0.0165      2  0.0001 CIMP    Preprocessor1_Model…
##  3 validation      0.996      0.00352     3  0.0001 CIMP    Preprocessor1_Model…
##  4 validation      0.984      0.0163      6  0.0001 CIMP    Preprocessor1_Model…
##  5 validation      0.907      0.0931     26  0.0001 CIMP    Preprocessor1_Model…
##  6 validation      0.985      0.0148     27  0.0001 CIMP    Preprocessor1_Model…
##  7 validation      0.869      0.131      28  0.0001 CIMP    Preprocessor1_Model…
##  8 validation      0.985      0.0148     35  0.0001 CIMP    Preprocessor1_Model…
##  9 validation      0.960      0.0403     42  0.0001 CIMP    Preprocessor1_Model…
## 10 validation      0.991      0.00921    48  0.0001 CIMP    Preprocessor1_Model…
## # … with 770 more rows

lr_auc <- 
  lr_res %>% 
  collect_predictions(parameters = lr_best) %>% 
  roc_curve(subtype, .pred_CIMP) %>% 
  mutate(model = "Logistic regression")

autoplot(lr_auc)

La regresión ridge resulta un tanto estática, aunque ofrece una buena métrica de rendimiento. Sin embargo, observar la penalziación del lasso puede brindarle mayor aprendizaje a la modelación.

6.2 Regresión lasso

Se crea el objeto en glmnet.

lr_mod <- 
  logistic_reg(penalty = tune(), mixture = 1) %>% 
  set_engine("glmnet")

Y se pasa bajo un objeto tipo workflow.

lr_workflow <- 
  workflow() %>% 
  add_model(lr_mod) %>% 
  add_recipe(tgexp_recipe2) 

Ahora se especificó una mixture = 1. También se vuelve a pasar la cuadrícula al workflow

lr_res <- 
  lr_workflow %>% 
  tune_grid(val_set,
            grid = lr_reg_grid,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))
lr_res 

## # Tuning results
## # Validation Set Split (0.8/0.2)  using stratification 
## # A tibble: 1 x 5
##   splits           id         .metrics         .notes          .predictions     
##   <list>           <chr>      <list>           <list>          <list>           
## 1 <split [102/26]> validation <tibble [30 × 5… <tibble [0 × 1… <tibble [780 × 6…

Y para finalizar, nuevamente se construyen los objetos respectivos al lasso para visualizar las métricas de rendimiento.

lr_plot <- 
  lr_res %>% 
  collect_metrics() %>% 
  ggplot(aes(x = penalty, y = mean)) + 
  geom_point() + 
  geom_line() + 
  ylab("Area under the ROC curve") +
  scale_x_log10(labels = scales::label_number())
lr_plot 

lr_best <- 
  lr_res %>% 
  collect_metrics() %>% 
  arrange(penalty)
lr_best

## # A tibble: 30 x 7
##     penalty .metric .estimator  mean     n std_err .config              
##       <dbl> <chr>   <chr>      <dbl> <int>   <dbl> <chr>                
##  1 0.0001   roc_auc binary     0.976     1      NA Preprocessor1_Model01
##  2 0.000127 roc_auc binary     0.976     1      NA Preprocessor1_Model02
##  3 0.000161 roc_auc binary     0.976     1      NA Preprocessor1_Model03
##  4 0.000204 roc_auc binary     0.976     1      NA Preprocessor1_Model04
##  5 0.000259 roc_auc binary     0.976     1      NA Preprocessor1_Model05
##  6 0.000329 roc_auc binary     0.976     1      NA Preprocessor1_Model06
##  7 0.000418 roc_auc binary     0.976     1      NA Preprocessor1_Model07
##  8 0.000530 roc_auc binary     0.976     1      NA Preprocessor1_Model08
##  9 0.000672 roc_auc binary     0.970     1      NA Preprocessor1_Model09
## 10 0.000853 roc_auc binary     0.970     1      NA Preprocessor1_Model10
## # … with 20 more rows

lr_auc <- 
  lr_res %>% 
  collect_predictions(parameters = lr_best) 
lr_auc

## # A tibble: 780 x 7
##    id         .pred_CIMP .pred_noCIMP  .row penalty subtype .config             
##    <chr>           <dbl>        <dbl> <int>   <dbl> <fct>   <chr>               
##  1 validation    0.00891 0.991            1  0.0001 CIMP    Preprocessor1_Model…
##  2 validation    1.00    0.000000141      2  0.0001 CIMP    Preprocessor1_Model…
##  3 validation    1.00    0.0000000876     3  0.0001 CIMP    Preprocessor1_Model…
##  4 validation    1.00    0.00000121       6  0.0001 CIMP    Preprocessor1_Model…
##  5 validation    1.00    0.00000662      26  0.0001 CIMP    Preprocessor1_Model…
##  6 validation    1.00    0.000000771     27  0.0001 CIMP    Preprocessor1_Model…
##  7 validation    0.984   0.0157          28  0.0001 CIMP    Preprocessor1_Model…
##  8 validation    1.00    0.000000420     35  0.0001 CIMP    Preprocessor1_Model…
##  9 validation    1.00    0.00000643      42  0.0001 CIMP    Preprocessor1_Model…
## 10 validation    1.00    0.000000263     48  0.0001 CIMP    Preprocessor1_Model…
## # … with 770 more rows

lr_auc <- 
  lr_res %>% 
  collect_predictions(parameters = lr_best) %>% 
  roc_curve(subtype, .pred_CIMP) %>% 
  mutate(model = "Logistic regression")

autoplot(lr_auc)

7 Referencias

• Akalin A. Computational Genomics with R. CRC Press, 2021: https://compgenomr.github.io/book/

• ISLR tidymodels Labs: https://emilhvitfeldt.github.io/ISLR-tidymodels-labs/index.html

• Lucas, B. A Tidy Introduction To Statistical Learning. 2020: https://beaulucas.github.io/tidy_islr/