Analítica Avanzada de Datos

Caso #3

Autor/a

Afiliación

Jorge I. Vélez, PhD

Universidad del Norte, Barranquilla

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.

Clasificación de las funciones de los lncRNAs. Imagen tomada de aquí.

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.

Una muestra de los datos puede obtenerse haciendo:

## disponibilidad del paquete data.table
if(!require(data.table)) install.packages('data.table')
require(data.table)

## lectura de datos con fread() de data.table
url <- 'https://www.dropbox.com/scl/fi/osljigmhxwcriewhpxq1s/AD_lncRNAs.txt?rlkey=j8cfuqsx51uz9hp3a7qvbzhst&dl=1'
d <- fread(url, header = TRUE)
d <- as.data.frame(d)

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')

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 
 

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.

resultados_t <- sapply(1:(ncol(d)-2), function(j) t.test(d[, 2 + j] ~ d$dx)$p.value)
significativos <- sum(resultados_t < 0.05)
p_ajustados <- p.adjust(resultados_t, method = "fdr")
comparacion <- data.frame(
  originales = resultados_t, 
  ajustados = p_ajustados
)

tabla_comparativa <- kable(head(comparacion), "html", caption = "p-valor original y ajustados  lncRNAs") %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%", height = "500px")

tabla_comparativa

p-valor original y ajustados lncRNAs

originales	ajustados
0.7273179	0.9731852
0.6163452	0.9690406
0.0008775	0.4359370
0.8337504	0.9850402
0.4156953	0.9547774
0.0133456	0.6887817

La acumulación de valores p ajustados cerca de 1 implica que ninguno de los lncRNAs serían considerados estadísticamente significativos después del ajuste FDR. Esto indica que podría ser necesario revisar la potencia estadística del estudio, posiblemente aumentando el tamaño de la muestra o revisando el diseño experimental.

Tipo de Valor p	Significativos
Original	1718
Ajustado	0

La diferencia dramática en el número de lncRNAs significativos antes y después del ajuste FDR subraya la importancia de este último en estudios de alta dimensionalidad. Sin el ajuste, es probable que muchos de los lncRNAs identificados como significativos simplemente sean falsos positivos causados por el azar debido al gran número de comparaciones.

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.

my_magical_function <- function(j) {
  modelo_rl <- glm(dx ~ d[, j + 2], family = "binomial", data = d)
  resultados <- summary(modelo_rl)$coefficients[2, c(1, 3, 4)]
  resultados<-c(paste("x",j),resultados)
  names(resultados) <- c("lncRNA","Coeficiente", "Estadístico_t", "Valor_p")
  return(resultados)
}

Resultados de la Regresión Logística para lncRNAs

lncRNA	Coeficiente	Estadístico_t	Valor_p
x 3	2.41602535206397	2.68721898136842	0.00720497004292235
x 4	0.137360252624823	0.219066724716692	0.826598071027626
x 5	0.591416878442149	0.827795094814097	0.407786553109876
x 6	2.9116293970772	2.19046044140059	0.0284908607216054
x 7	-0.702452979400386	-1.46588922458273	0.142678467983934
x 8	0.181979799139377	0.0978954230405511	0.922015334415466

Resumen de lncRNAs Significativos

Descripción	Cantidad
lncRNAs Significativos utilizando GLM	1304

La utilización de la regresión logística en este contexto ha revelado un número sustancial de lncRNAs (1304) que podrían estar relacionados significativamente con la enfermedad de Alzheimer, ofreciendo un enfoque más robusto y específico comparado con los resultados obtenidos por la prueba t de Student, especialmente después del ajuste por FDR.

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. Utilizamos la función scale() para estandarizar los niveles de expresión antes de realizar la prueba t, comparando los resultados con y sin

resultados_t_estandarizados <- sapply(3:ncol(d), function(j) {
  t.test(scale(d[, j]) ~ d$dx)$p.value
})

significativos_estandarizados <- sum(resultados_t_estandarizados < 0.05)

Número de lncRNAs Significativos con y sin Estandarización

Tipo	Significativos
Original	1718
Estandarizado	1718

La prueba t que es utilizada aquí para comparar las medias entre dos grupos (casos de enfermedad de Alzheimer vs. controles), es invariante a la transformación de estandarización. Esto se debe a que la estandarización no cambia la relación subyacente entre los grupos en términos de diferencias proporcionales

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.

Número de LncRNAs Significativos por Técnica

Técnica	LncRNAs Significativos
Prueba t de Student	1718
Regresión Logística	1304

Mientras que la prueba t podría sugerir que varios lncRNAs son significativos simplemente por diferencias de medias, la regresión logística puede refinar estos resultados al mostrar qué lncRNAs tienen un efecto significativo cuando se consideran todos los factores relevantes. por lo tantola Regresión Logística ofrece una visión más integral, considerando cómo cada lncRNA contribuye al diagnóstico de Alzheimer.

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.

resultados_lncRNA$Coeficiente <- as.numeric(resultados_lncRNA$Coeficiente)

lncRNA_sig <- resultados_lncRNA %>%
  arrange(Valor_p) %>%
  slice_head(n = 10) %>%
  mutate(Odds_Ratio = exp(Coeficiente))

Top 10 lncRNAs más Significativos con sus Odds Ratios

lncRNA	Coeficiente	Estadístico t	Valor p	Odds Ratio
x 4634	4.595273	3.00871221123096	0.00262357505403248	99.0151322
x 14508	-4.803780	-2.96482160667883	0.00302858496710436	0.0081987
x 6948	3.478922	2.86753787238015	0.00413679289685028	32.4247446
x 23214	2.134990	2.86444760984818	0.0041773725038525	8.4569666
x 26435	-7.160553	-2.86399449413214	0.00418335283013234	0.0007766
x 28007	2.404967	2.8514921169776	0.00435145652154728	11.0780629
x 11576	2.775605	2.8508258396312	0.00436058465323587	16.0483376
x 2114	4.361054	2.84968394380929	0.00437626922787505	78.3396976
x 27199	2.342667	2.84703679927126	0.00441282609773261	10.4089622
x 2237	1.686658	2.82929003618742	0.00466513997362706	5.4014005

El coeficiente positivo y el OR extremadamente alto para el lncRNA x4634 indican que un aumento en la expresión de este lncRNA está fuertemente asociado con un aumento en la probabilidad de tener EA. Específicamente, el OR de aproximadamente 99 sugiere que cada incremento unitario en la expresión de x4634 aumenta la probabilidad de tener EA en 99 veces, en comparación con los que tienen menores niveles de expresión de este lncRNA. Este lncRNA podría ser un potencial biomarcador de riesgo para la EA.

En contraste, el lncRNA x14508 muestra un coeficiente negativo y un OR muy bajo, lo que indica que un aumento en su expresión está asociado con una disminución en la probabilidad de tener EA. El OR de aproximadamente 0.0082 implica que cada incremento unitario en la expresión de x 14508 reduce la probabilidad de tener EA en más de 99%, lo que lo convierte en un potencial factor protector o biomarcador de protección contra la enfermedad.

lncRNAs_riesgo <- sum(lncRNA_sig$Coeficiente > 0)
lncRNAs_proteccion<- sum(lncRNA_sig$Coeficiente < 0)
riesgo_proteccion_df <- data.frame(
  Tipo = c("Confieren Riesgo", "Confieren Protección"),
  Cantidad = c(lncRNAs_riesgo, lncRNAs_proteccion)
)

kable(riesgo_proteccion_df, caption = "Clasificación de los Top 10 lncRNAs Significativos") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F, font_size = 14)%>%
  column_spec(1, bold = TRUE, width = "70%", color = "blue")

Clasificación de los Top 10 lncRNAs Significativos

Tipo	Cantidad
Confieren Riesgo	8
Confieren Protección	2

De los 10 lncRNAs mas significativos, 8 se clasifican como factores de riesgo y 2 como protectores. La mayoría de los lncRNAs identificados en este análisis son clasificados como factores de riesgo, Esto indica que estos lncRNAs, cuando están sobreexpresados o presentan cambios en su expresión, están asociados con un aumento en la probabilidad de desarrollar la enfermedad de Alzheimer.

P6

Construya un modelo de RL que incluya los lncRNAs seleccionados en P5. Evalúe la idoneidad. Concluya.

Respuesta.

# Crear un modelo incluyendo solo los 10 lncRNAs más significativos
modelo_final <- glm(dx ~ x4634+x14508+x6948+x23214+x26435+x28007+x11576+x2114+x27199+x2237, family = binomial(link = "logit"), data = d)
summary(modelo_final)

Call:
glm(formula = dx ~ x4634 + x14508 + x6948 + x23214 + x26435 + 
    x28007 + x11576 + x2114 + x27199 + x2237, family = binomial(link = "logit"), 
    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: 25

\[ \text{logit}(\hat{\theta}) = -62.34 + 1.604x_{4634} - 5.857x_{14508} + 2.065x_{6948} + 2.485x_{23214} + 9.805x_{26435} - 3.074x_{28007} + 2.158x_{11576} + 2.786x_{2114} + 1.370x_{27199} - 5.430x_{2237} \]

El modelo de regresión logística muestra signos de inestabilidad numérica con coeficientes y errores estándar extremadamente grandes, además de valores p y z constantes de 0 y 1 respectivamente para todos los predictores. Estos resultados sugieren problemas graves como multicolinealidad severa o separación perfecta en los datos, lo cual podría estar afectando la capacidad del modelo para estimar significancias estadísticas de manera adecuada.

Sin embargo, La devianza residual es extremadamente pequeña comparada con la devianza nula, lo que sugiere que el modelo tiene un buen ajuste a los datos.

## prueba de significancia global
Deviance <- with(modelo_final, null.deviance - deviance)
df <- with(modelo_final, df.null - df.residual)
(p_value <- 1 - pchisq(Deviance, df))

[1] 8.872984e-06

En este caso, la devianza reducida entre el modelo y el modelo nulo es significativa con un p-valor muy cercano a cero. Esto indica que el modelo que incluye los lncRNAs es significativamente mejor que el modelo nulo, lo que sugiere que al menos uno de los lncRNAs incluidos en el modelo tiene un efecto importante sobre la probabilidad de diagnóstico de enfermedad de Alzheimer.

Curva ROC y Área Bajo la Curva (AUC): La curva ROC es una herramienta crítica para evaluar la capacidad del modelo para clasificar correctamente los casos positivos y negativos.

## cargamos el paquete pROC
require(pROC)

Loading required package: pROC

Type 'citation("pROC")' for a citation.

Attaching package: 'pROC'

The following objects are masked from 'package:epicalc':

    auc, ci

The following objects are masked from 'package:stats':

    cov, smooth, var

## curva ROC
theta_hat <- predict(modelo_final, type = 'response')
ROC <- roc(dx ~ theta_hat, data = d)

Setting levels: control = 0, case = 1

Setting direction: controls < cases

plot(ROC, las = 1, col = "purple")

## cálculo del AUC
ROC$auc

Area under the curve: 1

Para este caso la AUC es de 1.0, lo que indica un modelo perfecto que distingue sin error entre los estados de la enfermedad.

Pseudo \(R^2\) El pseudo R^2 proporciona una medida de la varianza explicada por el modelo, similar al R^2 en la regresión lineal.

## paquete pscl
require(pscl)

Loading required package: pscl

Classes and Methods for R originally developed in the
Political Science Computational Laboratory
Department of Political Science
Stanford University (2002-2015),
by and under the direction of Simon Jackman.
hurdle and zeroinfl functions by Achim Zeileis.

## pseudo R^2
pR2(modelo_final)

fitting null model for pseudo-r2

          llh       llhNull            G2      McFadden          r2ML 
-2.623235e-10 -2.079442e+01  4.158883e+01  1.000000e+00  7.500000e-01 
         r2CU 
 1.000000e+00 

Un pseudo R^2 de McFadden de 1.0 indica un excelente ajuste del modelo, con el modelo explicando toda la variabilidad en la respuesta observada. Prueba de Hosmer-Lemeshow: Esta prueba evalúa la calibración del modelo, comparando los eventos observados con los predichos por el modelo en subgrupos de probabilidad.

## hosmer & lemeshow GoF test
if(!require(ResourceSelection)) install.packages('ResourceSelection')

Loading required package: ResourceSelection

ResourceSelection 0.3-6      2023-06-27

require(ResourceSelection)

## valor p de la prueba de H & L
hl <- hoslem.test(modelo_final$y, fitted(modelo_final), g = 11)
hl$p.value 

[1] 1

hl$observed

                                       
cutyhat                                 y0 y1
  [2.22044604925e-16,7.83813846239e-12] 11  0
  (7.83813846239e-12,2.9527725501e-11]   3  0
  (2.9527725501e-11,0.99999999997]       1  1
  (0.99999999997,0.999999999988]         0  3
  (0.999999999988,0.999999999995]        0  3
  (0.999999999995,1]                     0  8

Un p-valor de 1.0 indica una calibración perfecta, lo que sugiere que el modelo predice muy bien las probabilidades a lo largo de diferentes rangos.

El modelo muestra un ajuste excepcional con una capacidad excepcional para diferenciar entre las condiciones de la enfermedad, según lo evidenciado por la prueba de significancia global, la prueba de Hosmer-Lemeshow, la curva ROC y los valores de pseudo R^2. Estos resultados indican que el modelo no solo es estadísticamente significativo sino también prácticamente perfecto en términos de predicción y calibración.

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.

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 
 

cutoff <- seq(.1, .9, length = 500) 

##  cálculos
real <- d$dx
res <- data.frame(cutoff = cutoff, bycutoff2(modelo_final, real = real, cutoff))

##  gráfico
with(res, matplot(cutoff, cbind(sensitivity, specificity, class_rate), 
                  type = 'l', las = 1, ylim = c(0, 1), lty = 1, ylab = 'Value (%)'))
legend('bottomright', c('Se', 'Sp', 'Acc'), lty = 1, 
       col = 1:3, bty = 'n', ncol = 3)

la gráfica muestras que la sensibilidad (Se), la especificidad (Sp), y la precisión (Acc) son constantes y perfectas (1.0 o 100%) a lo largo de todos los posibles valores de cutoff desde 0.1 hasta 0.9 Esto implica que el modelo es excepcionalmente preciso en diferenciar los casos positivos de los negativos en todo el rango de valores de corte evaluados, no solo en el cutoff de 0.1.

res[which.min(abs(res$sensitivity - res$specificity)), ] 

  cutoff  a b c  d sensitivity specificity ppv npv fdr fpr class_rate lift
1    0.1 15 0 0 15           1           1   1   1   0   0          1    2
  delta
1     0

Al usar \(\theta_0 = 0.3\)

La matriz de confusión obtenida es:

## matriz de confusión para theta_0 = 0.5
y_pred <- 1*(theta_hat >= 0.3) 
y_real <- d$dx
confusion <- table(modelo = y_pred, real = y_real)[2:1, 2:1]
confusion

      real
modelo  1  0
     1 15  0
     0  0 15

Finalmente las medidas de desempeño, incluyendo el AUC, pueden calcularse como

## medidas de desempeño 
measureslog(confusion)

sensitivity specificity         ppv         npv         fdr         fpr 
          1           1           1           1           0           0 
 class_rate         auc 
          1           1 

vemos como desde diferente tecnicas estamos llegando a los mismos resultados, en donde el modelo puede identificar con precisión absoluta los casos con y sin enfermedad según los criterios establecidos por el cutoff de 0.3. Este resultado puede deberse a datos muy bien diferenciados o que el modelo esta sobreajustados.

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:

## ejemplo predicción
alzheimer(x0)

debería reportar 1 si la persona es diagnosticada con EA, y 0 en otro caso.

alzheimer <- function(x0) {
  if (is.vector(x0)) {
    x0 <- as.data.frame(t(x0))
    names(x0) <- c("x4634", "x14508", "x6948", "x23214", "x26435", 
                   "x28007", "x11576", "x2114", "x27199", "x2237")
  }
  probabilidad <- predict(modelo_final, newdata = x0, type = "response")
  cutoff <- 0.3  
  resultado <- ifelse(probabilidad >= cutoff, 1, 0)
  if (resultado == 1) {
    message("Diagnóstico: Positivo para la enfermedad de Alzheimer.")
  } else {
    message("Diagnóstico: Negativo para la enfermedad de Alzheimer.")
  }
  
  # Devolver el resultado
  return(resultado)
}

vector_x0 <- c(0.5, 98, 0.2, 0.1, 0.4, 0.3, 450, 0.2, 0.1, 10)
alzheimer(vector_x0)

Diagnóstico: Positivo para la enfermedad de Alzheimer.

1 
1 

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.

Primero, definimos una función que toma como argumento el número de la variable lncRNA que se quiere analizar utilizando un modelo de Regresión Lineal Simple (RLS) en lugar de Regresión Logística (RL):

my_linear_function <- function(j) {
  model <- lm(dx ~ d[, j], data = d)
  s_model <- summary(model)$coefficients[2, ]
  return(s_model[c(1, 3, 4)])  
}

A continuación, aplicamos esta función para cada uno de los lncRNAs en nuestros datos usando sapply:

linear_results <- t(sapply(3:(ncol(d)-2), my_linear_function))
linear_results_df <- as.data.frame(linear_results, row.names = NULL)
names(linear_results_df) <- c("Estimate", "t value", "Pr(>|t|)")

Evaluamos cuántos de estos lncRNAs resultan estadísticamente significativos bajo el modelo lineal:

significant_linear <- sum(linear_results_df$`Pr(>|t|)` < 0.05)
print(significant_linear)

[1] 1756

head(linear_results_df[order(linear_results_df$`Pr(>|t|)`), ], 10)

        Estimate   t value     Pr(>|t|)
4634   0.6220172  5.470868 7.668768e-06
24198  0.5999949  5.024129 2.596272e-05
3994   0.8053454  4.875338 3.899779e-05
12278 -1.3783024 -4.752761 5.451890e-05
22359  0.7057086  4.746792 5.541532e-05
2114   0.6877228  4.666934 6.892112e-05
6101   0.2524093  4.602605 8.214754e-05
1457   0.3652637  4.582242 8.683890e-05
6948   0.4283546  4.581873 8.692619e-05
9628   0.3492110  4.581048 8.712203e-05

Al comparar con los resultados de la regresión logística realizada en P2, observamos que aunque el modelo lineal identifica un número diferente de lncRNAs significativos, el uso de RLS en este contexto no es recomendable. Aunque técnicamente es posible implementar una RLS y obtener resultados, estos pueden ser engañosos o incorrectos para la interpretación en el contexto de una variable de respuesta dicotómica (presencia o ausencia de la enfermedad de Alzheimer). La RL es más apropiada porque modela la probabilidad de la presencia o ausencia de la enfermedad, lo que está directamente alineado con la naturaleza de nuestros datos. Por lo tanto, la RL no solo es más adecuada sino necesaria para obtener interpretaciones válidas y acciones basadas en los resultados del análisis.

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 Primero, ajustamos los modelos para obtener las estimaciones de los coeficientes y proceder con su interpretación:

# Modelo 1: incluye solo x1
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: 3

# Modelo 2: incluye x1 y x2
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: 4

Interpretación de los coeficientes:

\(\hat{\beta}_{1,M_1}\) en \(M_1\) indica el cambio en el logaritmo de las odds de dx por cada unidad incrementada en x1, manteniendo todas las otras variables constantes (en este caso, no hay más variables).

\(\hat{\beta}{1,M_2}\) en \(M_2\) ajusta la interpretación de x1, controlando adicionalmente por x2. Esto puede cambiar la magnitud o incluso el signo de \(\hat{\beta}{1,M_2}\) comparado con \(\hat{\beta}_{1,M_1}\), dependiendo de la correlación entre x1 y x2.

\(\hat{\beta}_{2,M_2}\) representa el cambio en el logaritmo de las odds de dx por cada unidad que aumenta en x2, manteniendo x1 constante. Este coeficiente es crucial para entender el efecto independiente de x2 sobre dx, ajustado por la presencia de x1.

coef_m1 <- coef(summary(m1))
coef_m2 <- coef(summary(m2))

or_m1 <- exp(coef_m1["x1", "Estimate"])
or_m2 <- exp(coef_m2[c("x1", "x2"), "Estimate"])


print(c("Odds ratio for x1 in M1:", or_m1))

[1] "Odds ratio for x1 in M1:" "0.748333313882021"       

print(c("Odds ratio for x1 in M2:", or_m2["x1"]))

                                                   x1 
"Odds ratio for x1 in M2:"        "0.867409948372411" 

print(c("Odds ratio for x2 in M2:", or_m2["x2"]))

                                                   x2 
"Odds ratio for x2 in M2:"        "0.547854357148312" 

Es esencial comparar los modelos \(M_1\) y \(M_2\) para entender cómo la inclusión de x2 afecta la relación entre x1 y dx. Si \(\hat{\beta}{1,M_1}\) y \(\hat{\beta}{1,M_2}\) difieren significativamente, esto podría indicar problemas de omisión de variable en \(M_1\) o de multicolinealidad en \(M_2\). Además, \(\hat{\beta}_{2,M_2}\) ayuda a evaluar si x2 tiene un efecto significativo cuando se ajusta por x1, lo cual es fundamental para modelar correctamente la influencia de los predictores sobre la respuesta.

Conclusión Este caso ha demostrado la utilidad de los análisis estadísticos avanzados en la identificación de lncRNAs relacionados con la Enfermedad de Alzheimer (EA), utilizando pruebas t y regresión logística para descubrir potenciales biomarcadores de riesgo y protección. Aunque las pruebas t iniciales sugirieron numerosos lncRNAs significativos, el ajuste por FDR reveló la necesidad de más pruebas debido a los altos falsos positivos.

Los modelos de regresión logística demostraron ser herramientas valiosas, ofreciendo una comprensión más profunda y modelos predictivos precisos para la EA, validados a través de métricas como ROC y AUC, y la prueba de Hosmer-Lemeshow. Los resultados confirmaron que la regresión logística es más adecuada que la regresión lineal para datos dicotómicos, enfatizando la importancia de elegir el modelo estadístico correcto basado en la naturaleza de los datos.