Analítica Avanzada de Datos
Caso #3
Introducción
Los ARN no codificantes largos (lncRNAs en inglés) comprenden un vasto repertorio de ARN que desempeñan una amplia variedad de funciones cruciales en la fisiología de los tejidos de manera específica para cada célula. A pesar de estar implicados en miles de mecanismos reguladores, todavía no se ha asignado ninguna función a muchos lncRNAs.
Cada vez hay más pruebas de que los lncRNAs se expresan de forma aberrante en la progresión de la Enfermedad de Alzheimer (EA). En nuestro grupo, hemos analizado la expresión de cientos de lncRNAs en 15 individuos diagnosticados con EA, y 15 voluntarios sin la enfermedad en el marco del proyecto de investigación “Nuevos ARN no codificantes exo- somales y su papel en la patogénesis de la Enfermedad de Alzheimer”, código 121584468097 (844/2019), contract 416-2020, otorgado al Grupo en Genética y Medicina Molecular (Departamento de Medicina), y al Grupo de Productividad y Competitividad (Departamento de Ingeniería Industrial) de la Universidad del Norte, Barranquilla, Colombia.
Las variables relevantes son
id: identificador del participante.dx: diagnóstico del participante (0: no EA;1: EA)x1,x2,…,x29809: niveles de expresión del lncRNA \(j\), con \(j=1,2,\ldots,29809\).
El objeto d tiene 30 filas y 29811 columnas.
Importante
A partir de los datos, el objetivo es identificar los lncRNAs cuya expresión incrementa o reduce la probabilidad de diagnóstico de EA.
Fecha de entrega: Junio 11, 2024.
Como la variable respuesta es binaria, podemos utilizar las funciones que ya hemos empleado con anterioridad en Regresión Logística (RL) para medir el desempeño del modelo.
## load functions
source('https://www.dropbox.com/s/xclvdugfbrf5ryn/logistic-functions.R?dl=1')Preguntas
P1
Use una prueba \(t\) de Student de 2 colas para identificar los lncRNAs diferencialmente expresados en personas con EA. Extraiga el valor \(p\) y compárelo con el nivel de significancia \(\alpha = 0.05\). Cuántos de estos son estadísticamente significativos? Use la función p.adjust() para ajustar los valores \(p\) utilizando el método FDR. Compare.
t_test_results <- sapply(1:(ncol(d)-2), function(j) t.test(d[ ,2 +j] ~ dx, d)$p.value)
significant <- length(which(t_test_results < 0.05))Con esto, se tiene que hay 1718 lncRNAs estadísticamente significativos.
A continuación, se comparan los p valores originales como los p valores ajustados con el método FDR.
p_values_fdr <- p.adjust(t_test_results, method = "fdr")
comparison <- data.frame(
original = t_test_results, fdr = p_values_fdr
)%>%
mutate(dif = fdr - original)
head(comparison, 10)De la tabla, es posible obtener que el mínimo valor que se obtiene de los valores ajustados es de 0.3028446. Es decir, ningún p valor ajustado por el método FDR es significativo. En la siguiente gráfica se puede visualizar la relación entre los p valores originales con los ajustados, en la que se observa cómo estos últimos son considerablemente mayores en todos los casos.
ggplot(comparison, aes(x = original, y = fdr)) +
geom_point() +
geom_abline(slope = 1, intercept = 0, color = "red", ) +
labs(title = "Comparación de p valores originales con p valores ajustados",
x = "p valores originales",
y = "p valores ajustados con FDR") +
theme_minimal() +
coord_cartesian(ylim = c(min(comparison$fdr) - min(comparison$fdr), max(comparison$fdr) + 0))
A través de estos resultados, parece concluirse que todos los hallazgos de significancia de la prueba original son falsos positivos, pues ningún p valor mantuvo su significancia. Sin embargo, el hecho de que no halla resultados significativos en los valores ajustados se puede deber a que la muestra es muy pequeña para realizar la comparación. En este caso, se cuenta con 30 observaciones contra casi 30 mil variables.
P2
Construya una función que, a partir del número del lncRNA, estime el modelo de RL, extraiga el coeficiente estimado, el estadístico \(t\) y el valor \(p\). Por ejemplo, para el lncRNA \(j\), aplicar my_magical_function(j) debe producir \(\hat{\beta}_j\), \(t_j\) y \(p_j\).
Respuesta.
Para esta función, se toma como argumento el número del lncRNa, y se calcula el modelo de regresión logística. Luego, se extraen el coeficiente estimado, el estadístico \(t\) y el valor \(p\).
my_magical_function <- function(j){
model <- glm(dx ~ d[ , 2+j], family = "binomial", data = d)
summary(model)$coefficients[2, ][c(1, 3, 4)]
}Ahora, se utiliza la función para hallar los p valores que midan la significancia de cada una de las moléculas bajo el efecto conjunto del resto de ellas.
test_total <- sapply(1:(ncol(d)-2), my_magical_function)
test_total <- as.data.frame(t(test_total))
head(test_total, 10)De acuerdo a estos resultados, hay 1304 moleculas significativas para el diagnóstico de EA.
P3
Realice el mismo análisis que en P1, pero ahora estandarice los niveles de expresión. Apóyese en la función scale(). Considera usted que la estandarización es necesaria? Por qué si? Por qué no? Compare los resultados cuando usa estandarización y cuando no.
Respuesta.
Para este ejercicio, se utiliza el mismo código que en P1, con la única diferencia de que en este caso los valores de \(x_j\) están estandarizados.
t_test_results_scale <- sapply(1:(ncol(d)-2), function(j) t.test(scale(d[ ,2 +j]) ~ dx, d)$p.value)
significant_scale <- length(which(t_test_results_scale < 0.05))A través de estos cálculos, se obtiene que existen 1718 moléculas significativas para el diagnóstico.
Con los resultados anteriores se puede ver que: escalar los datos no influye en el número de moléculas significativas para el diagnóstico. Por lo tanto, es viable afirmar que la estandarización no es necesaria.
A continuación se puede ver una tabla que resume los valores obtenidos con los datos originales y los datos estandarizados.
comp_scaled_not <- data.frame(Original = t_test_results, Estandarizados = t_test_results_scale)
head(comp_scaled_not, 10)P4
En P1 se usó una prueba \(t\) de Student y en P2 un modelo de RL. Explique por qué hay diferencias en el número de lncRNAs que resultan ser significativos al emplear estos métodos.
Respuesta.
La diferencia principal se encuentra en que aplicando el t.test, se asumen que el resto de condiciones de las observaciones son homogeneas. En este caso, se consideraría que todos los individuos presentan la misma cantidad en cada lncRNA distinto al que se evalua en la prueba. Por el contrario, por medio del modelo de RL se tienen en consideración el efecto del resto de moléculas sobre el lncRNA estudiado. La consideración de estas otras moléculas es importante por el hecho de que al entrar en contacto, estas pueden contrarestar o potenciar el efecto de cada una en la absorción de componentes de las células.
En otras palabras, realizar el t.test permite determinar si hay diferencias significativas entre cada una de las moléculas, pero el modelo RL halla la probabilidad de obtener la diagnóstico teniendo en cuenta la existenciae influencia de otras moléculas dentro de la células.
P5
Ordene los resultados obtenidos de acuerdo con el valor \(p\) y extraiga los 10 lncRNAs más significativos de acuerdo con el valor \(p\). Cuántos de ellos confieren riesgo y cuántos protección contra EA? Concluya.
Respuesta.
Primero, se seleccionan las variables con los menores p valores, las cuales representan los lncRNAs más significativos para el diagnóstico.
colnames(test_total) <- c("coeff", "valor_z", "p_value")
test_total <- test_total%>%
mutate(id = seq(1, (ncol(d)-2), 1))
p_value_sort <- test_total%>%
arrange(p_value)
head(p_value_sort, 10)Ahora, para determinar cuáles confieren riesgo y protección contra EA y cuántos hay de cada uno, se procede de la siguiente manera.
Se crea un data frame con la información de los lncRNAs más significativos:
top_10 <- p_value_sort%>%
head(10)lncRNAs que confieren protección:
length(which(top_10$coeff < 0))[1] 2lncRNAs que confieren riesgo:
length(which(top_10$coeff > 0))[1] 8De los resultados obtenidos se puede observar que el 80% de los lncRNAs más significativos para el diagnóstico confieren un riesgo. Que si bien las moléculas que confieren protección presentan los coeficientes con mayor incidencia con -4.804 y -7.161 para los lncRNAs 14508 y 26435 respectivamente, no es posible contrarestar con el efecto de riesgo que aportan el resto de lncRNAs.
En efecto,
Influencia de protección:
\[ \begin{align} \hat{OR}_{14508} & = e^{-4.804} = 0.00819 = \frac{P(\text{estar enfermo})}{P(\text{estar sano})} \\ & \Longrightarrow 121.9974 = \frac{P(\text{estar sano})}{P(\text{estar enfermo})} \end{align} \] Por cada unidad de cantidad que aumente en el lcnRNA 14508, el paciente tiene 122 veces más de probabilidad de estar sano. El lncRNA 14508 es un factor protector.
\[ \begin{align} \hat{OR}_{26435} & = e^{-7.1605} = 0.000776 = \frac{P(\text{estar enfermo})}{P(\text{estar sano})} \\ & \Longrightarrow 1287.555 = \frac{P(\text{estar sano})}{P(\text{estar enfermo})} \end{align} \] Por cada unidad de cantidad que aumente en el lcnRNA 26435, el paciente tiene 1287 veces más de probabilidad de estar sano. El lncRNA 26435 es un factor protector.
Influencia de riesgo:
\[ \hat{OR}_{4634} = e^{4.595} = 98.98814 = \frac{P(\text{estar enfermo})}{P(\text{estar sano})} \] Por cada unidad de cantidad que aumente en el lcnRNA 4634, el paciente tiene 99 veces más de probabilidad de estar enfermo. El lncRNA 4634 es un factor de riesgo.
\[ \hat{OR}_{6948} = e^{3.789} = 44.2122 = \frac{P(\text{estar enfermo})}{P(\text{estar sano})} \] Por cada unidad de cantidad que aumente en el lcnRNA 6948, el paciente tiene 44 veces más de probabilidad de estar enfermo. El lncRNA 6948 es un factor de riesgo.
\[ \hat{OR}_{23214} = e^{2.135} = 8.4570 = \frac{P(\text{estar enfermo})}{P(\text{estar sano})} \] Por cada unidad de cantidad que aumente en el lcnRNA 23214, el paciente tiene 8 veces más de probabilidad de estar enfermo. El lncRNA 23214 es un factor de riesgo.
\[ \hat{OR}_{28007} = e^{2.405} = 11.078 = \frac{P(\text{estar enfermo})}{P(\text{estar sano})} \] Por cada unidad de cantidad que aumente en el lcnRNA 28007, el paciente tiene 11 veces más de probabilidad de estar enfermo. El lncRNA 28007 es un factor de riesgo.
\[ \hat{OR}_{11576} = e^{2.776} = 16.055 = \frac{P(\text{estar enfermo})}{P(\text{estar sano})} \] Por cada unidad de cantidad que aumente en el lcnRNA 11576, el paciente tiene 16 veces más de probabilidad de estar enfermo. El lncRNA 11576 es un factor de riesgo.
\[ \hat{OR}_{2114} = e^{4.361} = 78.335 = \frac{P(\text{estar enfermo})}{P(\text{estar sano})} \] Por cada unidad de cantidad que aumente en el lcnRNA 2114, el paciente tiene 78 veces más de probabilidad de estar enfermo. El lncRNA 2114 es un factor de riesgo.
\[ \hat{OR}_{27199} = e^{2.343} = 10.412 = \frac{P(\text{estar enfermo})}{P(\text{estar sano})} \] Por cada unidad de cantidad que aumente en el lcnRNA 27199, el paciente tiene 10 veces más de probabilidad de estar enfermo. El lncRNA 27199 es un factor de riesgo.
\[ \hat{OR}_{2237} = e^{1.687} = 5.403 = \frac{P(\text{estar enfermo})}{P(\text{estar sano})} \] Por cada unidad de cantidad que aumente en el lcnRNA 2237, el paciente tiene 5 veces más de probabilidad de estar enfermo. El lncRNA 2237 es un factor de riesgo.
P6
Construya un modelo de RL que incluya los lncRNAs seleccionados en P5. Evalúe la idoneidad. Concluya.
Respuesta.
Se construye el modelo con los 10 lncRNAs seleccionados en el punto 5, los cuales son los más significativos dentro del total.
modelo_p6 <- glm(dx ~ x4634 + x14508 + x6948 + x23214 + x26435 + x28007 + x11576 + x2114 + x27199 + x2237, family = "binomial", data = d)
summary(modelo_p6)
Call:
glm(formula = dx ~ x4634 + x14508 + x6948 + x23214 + x26435 +
x28007 + x11576 + x2114 + x27199 + x2237, family = "binomial",
data = d)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -6.234e+01 3.853e+06 0 1
x4634 1.604e+00 6.266e+05 0 1
x14508 -5.857e+01 2.001e+05 0 1
x6948 2.065e+01 1.444e+05 0 1
x23214 2.485e+01 1.879e+05 0 1
x26435 9.805e+00 5.700e+05 0 1
x28007 -3.074e+01 3.988e+05 0 1
x11576 2.158e+01 2.133e+05 0 1
x2114 2.786e+01 2.074e+05 0 1
x27199 1.370e+01 1.496e+05 0 1
x2237 -5.430e+00 2.308e+05 0 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 4.1589e+01 on 29 degrees of freedom
Residual deviance: 5.2465e-10 on 19 degrees of freedom
AIC: 22
Number of Fisher Scoring iterations: 25Se observa que todos los p valores son 1, lo que ha un indicio de un mal funcionamiento del modelo.
Ahora, al aplicar la prueba de Hosmer-Lemeshow, se obtienen los siguientes resultados:
## paquete ResourceSelection
require(ResourceSelection)Loading required package: ResourceSelectionResourceSelection 0.3-6 2023-06-27## valor p de la prueba de Hosmer & Lemwshow
g <- length(coefficients(modelo_p6)) + 1
hl <- hoslem.test(modelo_p6$y, fitted(modelo_p6), g = g)
hl$p.value[1] 1Al ser el p valor 1 > 0.05, se podría decir que el modelo está muy bien calibrado. Algo que no concuerda con el panorama visto a través de los p valores del modelo. Este mismo comportamiento se da al analizar la gráfica ROC.
require(pROC)Loading required package: pROCType 'citation("pROC")' for a citation.
Attaching package: 'pROC'The following objects are masked from 'package:stats':
cov, smooth, vartheta_hat <- predict(modelo_p6, type = 'response')
ROC <- roc(dx ~ theta_hat, data = d)Setting levels: control = 0, case = 1Setting direction: controls < casesplot(ROC, las = 1, col = "navy")
Se ve que el valor AUC es 1, dando a entender que el modelo explica perfectamente la variabilidad.
En casos como estos, es prudente sospechar que algo no está bien, ya sea con el modelo o con los datos en general. Una posible razón de esto, como se mencionó anteriormente, es la poca cantidad de observaciones en los datos. Además, el contar con valores \(z = 0\) para todas las variables, permite concluir que no hay efecto significativo en la variable respuesta al no tener en cuenta al presencia del resto de lncRNAs en el modelo.
Por lo tanto, el modelo que considera las variables halladas en P5 no está bien calibrado.
P7
Utilizando el criterio \(S_e \approx S_p\), se determinó que el cutoff óptimo como \(\theta_0 = 0.3\) con base en el modelo modelo de RL ajustastado en P5. Está usted de acuerdo? Concluya.
Respuesta.
Primero, se descargan algunas funciones que ayudan con los análisis:
This is a collection of functions created by
Professor Jorge Vélez, PhD <https://jorgeivanvelez.netlify.app/>
for the fitting and validation of Binary Classification models
In case of problems, please an email to
<jvelezv@uninorte.edu.co> with enough details
Last modification: April 23, 2024
Ahora, para verificar que el cutoff sea en efecto 0.3, se deben condiderar todos los posbiles valores que este puede tomar. Eso se hace mediante el siguiente código.
## definimos 500 posibles valores para el cutoff
cutoff <- seq(.1, .9, length = 500)
## cálculos
real <- d$dx
res <- data.frame(cutoff = cutoff, bycutoff2(modelo_p6, real = real, cutoff))Utilizando el criterio \(S_e \approx S_p\), se obtienen los siguientes resultados:
res[which.min(abs(res$sensitivity - res$specificity)), ] De lo anterior se puede observar que el valor del cutoff en este caso es de 0.1, no de 0.3 como se afirmaba en un comienzo. Sin embargo, al momento de realizar la clasificación en la matriz de confusión, la escogencia del cutoff no presenta alguna incidencia en la misma. Esto se debe que a para cualquier valor del cutoff la clasificación es la misma.
tail(res, 15)En la tabla anterior se puede ver cómo incluso usando los valores de cutoff más altos, la clasificación en la matriz de confusión sigue siendo de 15 verdaderos positivos y 15 verdaderos negativos. Esto puede dar explicación a los resultados obtenidos en el resumen del modelo, en el cual es observaba una separación perfecta de las categorías, lo que causaba un valor de \(z = 0\) para todas las variables.
P8
Construya una función que, a partir de los niveles de expresión de lncRNAs en una persona, determine el diagnóstico de EA. Tenga en cuenta en valor de \(\theta_0\) derivado en P7. Si \(\mathbf{x}_0\) representa el vector de dichos niveles de expresión, entonces la función debería reportar 1 si la persona es diagnosticada con EA, y 0 en otro caso.
alzheimer <- function(x0){
theta0 <- 0.1
x0 <- as.data.frame(matrix(x0, nrow = 1))
colnames(x0) <- c("x4634", "x14508", "x6948", "x23214", "x26435", "x28007",
"x11576", "x2114", "x27199", "x2237")
mod <- predict(modelo_p6, newdata = x0,
type = 'response', se.fit = FALSE)
1 * (mod >= theta0)
}
## ejemplo predicción
x0 <- c(0.564, 3.555, 4.694, 0.457, 1.853, 3.459, 2.409, 1.495, 3.023, 0.474)
alzheimer(x0)1
0 P9
Repita P2 pero use un modelo de Regresión Lineal Simple(RLS) en lugar de un modelo de RL. Compare los resultados y concluya. Con base en los resultados obtenidos, cree usted que tiene sentido (o hace alguna diferencia) usar RLS en lugar de RL?. Justifique.
Respuesta.
Se define la función, tal como en el punto, usando la función lm() en lugar de glm().
my_linear_function <- function(j){
model <- lm(dx ~ d[ , 2+j], data = d)
summary(model)$coefficients[2, ][c(1, 3, 4)]
}Luego, se procede a aplicar la función para cada uno de los lncRNAs de los datos.
test_linear <- sapply(1:(ncol(d)-2), my_linear_function)
test_linear <- as.data.frame(t(test_linear))
head(test_linear, 10)Al igual que en el punto 2, se muestran los lncRNAs más significativos para el modelo.
colnames(test_linear) <- c("coeff", "valor_z", "p_value")
test_linear <- test_linear%>%
mutate(id = seq(1, (ncol(d)-2), 1))
p_value_sort_linear <- test_linear%>%
arrange(p_value)
head(p_value_sort_linear, 10)Se puede observar que a pesar de obtener p valores más pequeños para determinar la significancia, el efecto que presentan los lncRNAs para el modelo de regresión lineal simple es muy pequeño. Donde, de hecho, el máximo y mínimo valor de los coeficientes en este caso son:
c(min(p_value_sort_linear$coeff),max(p_value_sort_linear$coeff))[1] -1.378302 1.326080De los cuales, solo el mínimo se ve expresado dentro de los 10 lncRNAs más significativos.
Esto, más el hecho de que una regresión lineal simple requiere de una variable respuesta con naturaleza numérica y continua, lo cual no se cumple en este estudio, hace que un modelo de esta índole no sea conveniente para determinar los resultados que se pretenden.
P10
Suponga que se ajustan los modelos
\(M_1:\) glm(dx ~ x1, family = 'binomial', data = d)
\(M_2:\) glm(dx ~ x1 + x2, family = 'binomial', data = d)
tal que \(M_1 \subseteq M_2\).
Después del proceso de estimación, se obtuvieron los coeficientes \((\hat{\beta}_{0,M_1}, \hat{\beta}_{1,M_1})\) y \((\hat{\beta}_{0,M_2}, \hat{\beta}_{1,M_2}, \hat{\beta}_{2,M_2})\) con base en una muestra de tamaño \(n\) almacenada en el objeto d.
Interprete \(\hat{\beta}_{1,M_1}\) y \(\hat{\beta}_{1,M_2}\). Qué representa \(\hat{\beta}_{2,M_2}\)?.
Respuesta.
Aplicación de \(M_1\)
m1 <- glm(dx ~ x1, family = 'binomial', data = d)
summary(m1)
Call:
glm(formula = dx ~ x1, family = "binomial", data = d)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.9932 2.7555 0.360 0.719
x1 -0.2899 0.7972 -0.364 0.716
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 41.589 on 29 degrees of freedom
Residual deviance: 41.456 on 28 degrees of freedom
AIC: 45.456
Number of Fisher Scoring iterations: 3Aplicación de \(M_1\)
m2 <- glm(dx ~ x1 + x2, family = 'binomial', data = d)
summary(m2)
Call:
glm(formula = dx ~ x1 + x2, family = "binomial", data = d)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.8105 9.7775 0.492 0.623
x1 -0.1422 0.8754 -0.162 0.871
x2 -0.6017 1.4762 -0.408 0.684
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 41.589 on 29 degrees of freedom
Residual deviance: 41.288 on 27 degrees of freedom
AIC: 47.288
Number of Fisher Scoring iterations: 4Interpretaciones:
- Por cada unidad de aumento en la cantidad del lncRNA 1, la probabilidad de ser diagnosticado con EA se reduce en un %33.63
\[ OR_{1, M_1} = \frac{1}{e^{-0.2899}} = 1.3363 = \frac{P(\text{estar sano})}{P(\text{estar enfermo})} \]
- Para un valor fijo en la cantidad del lncRNA 2, por cada unidad de aumento en la cantidad del lncRNA 1, la probabilidad de ser diagnosticado con EA se reduce en un %15.28.
\[ OR_{1, M_2} = \frac{1}{e^{-0.1422}} = 1.1528 = \frac{P(\text{estar sano})}{P(\text{estar enfermo})} \]
- En este modelo, \(\hat{\beta}_{2,M_2}\) representa qué tanto cambia el \(logit(\theta)\) de ser diagnosticado con EA por cada unidad de aumento en la cantidad del lncRNA 2, con una cantidad fija del lncRNA 1. En particular, es posible decir que el \(logit(\theta)\) disminuirá en 0.6017 unidades.