Este documento contiene las instrucciones y enunciados asociados a la Prueba Especial Programada 1 - Parte práctica (PEP1) del curso de Estadística Computacional. Con respecto a ello:
La inferencia filogenética es un área de la Bioinformática que busca proponer una hipótesis para explicar las relaciones evolutivas entre organismos. Estas relaciones generalmente son representadas por medio de un árbol filogenético, que puede ser construido usando las características morfológicas de los organismos o sus secuencias moleculares: ADN, ARN o aminoácidos.
El trabajo desarrollado por (Villalobos-Cid et al. 2020) encontró que usando un conjunto específico de 22 genes de los más de 6,015 que componen el genoma de la levadura Saccharomyces cerevisiae, clave en la industria productora de bebidas alcohólicas, se puede llegar a construir más del 30% de su historia evolutiva. A raíz de ello, el equipo de microbiología ha planteado la siguiente pregunta: ¿existirán otros grupos de 22 genes que alcancen igual o mayor porcentaje de reconstrucción de la historia evolutiva de esta levadura? Bajo este contexto:
Primero que todo, se cargarán las bibliotecas a usar en la resolución de las preguntas.
#Bibliotecas
library("ggplot2") #Biblioteca para graficar
library("plotly") #Biblioteca para graficar
#library("nortest") #Biblioteca
library("psych") #Biblioteca par estadísticosCon respecto a la parte A del enunciado, el número de combinaciones se puede calcular con la siguiente fórmula \(C_n=\frac{n!}{x!(n-x)!}\),en que \(n\) es igual a 100 y \(x\) es 22. Esto es:
#===============
# Respuesta - Pregunta 1 - Parte A
#===============
# Fórmula Cn = n! / (x! * (n-x)!)
x=22
n=100
combinaciones = factorial(n)/(factorial(x)*factorial(n-x))Lo que implica que se deberán evaluar 7.3320669^{21} combinaciones de genes. Haga conciencia de esa cifra ¡Es enorme!
Aquí se nos solicita calcular el tiempo que tardaría un algoritmo exhaustivo en recorrer todas las combinaciones, dado que evaluar una de ellas requiere 1 nanosegundo. Esto se puede calcular como:
#===============
# Respuesta - Pregunta 1 - Parte B
#===============
tiempo_segundos = (combinaciones * 10^-9)
tiempo_horas = tiempo_segundos/3600
tiempo_dias = tiempo_horas/24
tiempo_anos = tiempo_dias/365Esto implica que el algoritmos demorará 2.3249832^{5} años en evaluar todas las combinaciones. Esto equivale a aproximadamente 232.5 milenios. Sí, esperemos sentaditos(as) a que termine.
Los marcapasos cardíacos tienen una vida media de 10 años libres de fallas, por lo que deben ser reemplazados quirúrgicamente para evitar problemas en la salud de los y las pacientes (Dutta and Barman 2021). Dado esto:
Al hablar de vida media asociada a años como unidad de tiempo, inmediatamente podríamos considerar una distribución Exponencial o de Poisson, de \(\lambda\) igual a 10. En este caso, asumiremos una distribución discreta, por lo tanto, las tres primeras preguntas pueden ser fácilmente respondidas con la función ppois(), que calcula la probabilidad acumulada hasta una entrada de valor determinado (\(prob(x<z)\)). A raíz de ello, la primera parte se puede calcular como:
#===============
# Respuesta - Pregunta 2 - Parte A
#===============
# Se asume distribución de Poisson
lambda=10
anos=3
respuesta = ppois(anos,lambda,lower.tail = T)Resultando que la probabilidad de que un marcapasos falle antes de los 3 años es 0.01. ¿Medicamente esta probabilidad será alta o baja? Desde mi punto de vista es baja aún. Es por ello que algunos(as) expertos(as) recomiendan cambiar el equipo antes de los 5 años de uso.
Se puede seguir la misma lógica de la respuesta anterior para responder la parte siguiente. En este caso:
#===============
# Respuesta - Pregunta 2 - Parte B
#===============
# Se asume distribución de Poisson
lambda=10
anos=10
respuesta = ppois(anos,lambda,lower.tail = T)La probabilidad de que un marcapasos falle antes de los 10 años es 0.583, que es extremadamente alta.
Este caso es similar al anterior, sin embargo, se debe ajustar la colas (lower.tail = F) de la función ppois() para calcular \(prob(x>Z)\), o usar su complemento: \(1-prob(x<z)\).
#===============
# Respuesta - Pregunta 2 - Parte C
#===============
# Se asume distribución de Poisson
lambda=10
anos=15
respuesta = ppois(anos,lambda,lower.tail = F)La probabilidad de que un marcapasos falle después de los 15 años es 0.049. ¿Esto es bueno? Claro que no, ya que la probabilidad de que falle antes de loas 15 años es \(1-prob(x<15)=0.951\).
El gráfico que representa la distribución es el siguiente:
#===============
# Respuesta - Pregunta 2 - Parte D
#===============
#Datos
anos=seq(1:20)
lambda=10
distribucion = dpois(anos,lambda)
datos=data.frame(anos,distribucion)
#Gráfico
grafico = ggplot(data=datos,aes(x=anos,y=distribucion))
grafico = grafico + geom_bar(stat="identity",fill="lightblue3")
grafico = grafico + theme_bw() + ggtitle("Distribución de probabilidades para cambio de marcapasos")
grafico = grafico + xlab("Años") + ylab("Probabilidad")
ggplotly(grafico)La tasa de incidencia de la tuberculosis en Chile es de 14 casos por cada 100,000 personas (Herrera M. 2020). Dado esto:
Al preguntarnos por probabilidades sin reposición, podemos inmediamente asociar una distribución hipergeométrica al problema, con probabilidad base \(p(x)=14/100000=0.00014\) y acumulada calculable con la función phyper().
En este caso podemos usar la probabilidad base para obtener los casos asociados a la categoría “tener tuberculosis” (A), considerando como total las 30,000 personas. Su complemento serán las personas que no poseen esta patología (B).
#===============
# Respuesta - Pregunta 3 - Parte A
#===============
A=round(14/100000*30000,0) #Total de personas con tuberculosis. Se aproxima a un valor entero.
B=30000-A #Total de personas sin tuberculosis.
pac_tub = 1 #Total de personas que esperamos que tengan tuberculosis (éxito).
pac_muestreados = 6 #Total de pacientes muestreados (intentos).
respuesta=dhyper(x=pac_tub, m=A, k=pac_muestreados, n=B)Como resultado, se obtiene que la probabilidad de que al seleccionar seis pacientes sin reposición, uno de ellos tenga tuberculosis es 7.9960004^{-4}. Este valor tan pequeño se debe al pequeño tamaño de B en la población/muestra.
El mismo procedimiento se puede aplicar para responder la segunda parte:
#===============
# Respuesta - Pregunta 3 - Parte B
#===============
A=round(14/100000*30000,0) #Total de personas con tuberculosis. Se aproxima a un valor entero.
B=30000-A #Total de personas sin tuberculosis.
pac_tub = 3 #Total de personas que esperamos que tengan tuberculosis (éxitos).
pac_muestreados = 100 #Total de pacientes muestreados (intentos).
respuesta=dhyper(x=pac_tub, m=A, k=pac_muestreados, n=B)En este caso, la probabilidad de que al seleccionar 100 pacientes sin reposición y que tres de ellos tenga tuberculosis es 1.4328288^{-7}.
En ambos casos se usa una distribución hypergeométrica. Primero se graficará la probabilidad de obtener pacientes con tuberculosis considerando una muestra de 6 personas, seleccionadas sin reposición.
#===============
# Respuesta - Pregunta 3 - Parte C
#===============
#Datos
pac_tub = seq(0,5)
pac_muestreados=6 #Pacientes muestrados
distribucion = dhyper(x=pac_tub, m=A, k=pac_muestreados, n=B)
datos=data.frame(pac_tub,distribucion)
#Gráfico
grafico = ggplot(data=datos,aes(x=pac_tub,y=distribucion))
grafico = grafico + geom_bar(stat="identity",fill="lightblue3")
grafico = grafico + theme_bw() + ggtitle("Distribución de probabilidades")
grafico = grafico + xlab("Prob. de obtener entre 0 a 5 enferm@s en al seleccionar 6 sin rep.") +
ylab("Probabilidad")
ggplotly(grafico)Ahora se graficará la probabilidad de obtener pacientes con tuberculosis considerando una muestra de 100 personas seleccionadas sin reposición.
#===============
# Respuesta - Pregunta 3 - Parte C
#===============
#Datos
pac_tub = seq(0,5)
pac_muestreados=100 #Pacientes muestrados
distribucion = dhyper(x=pac_tub, m=A, k=pac_muestreados, n=B)
datos=data.frame(pac_tub,distribucion)
#Gráfico
grafico = ggplot(data=datos,aes(x=pac_tub,y=distribucion))
grafico = grafico + geom_bar(stat="identity",fill="lightblue3")
grafico = grafico + theme_bw() + ggtitle("Distribución de probabilidades")
grafico = grafico + xlab("Prob. de obtener entre 0 a 10 enferm@s en al seleccionar 100 sin rep.") +
ylab("Probabilidad")
ggplotly(grafico)En ambos casos la probabilidad de tener pacientes con tuberculosis se acerca a cero, ya que es mucho más probable seleccionar a pacientes saludables \((x=0)\).
También se puede hacer el ejercicio de dejar fijo el número de casos de éxito y variar el número de pacientes seleccionados. Por ejemplo, veamos el comportamiento al considerar sólo un paciente con tuberculosis.
#===============
# Respuesta - Pregunta 3 - Parte C
#===============
#Datos
pac_tub = 1
pac_muestreados=seq(0,30000) #Pacientes muestrados
distribucion = dhyper(x=pac_tub, m=A, k=pac_muestreados, n=B)
datos=data.frame(pac_muestreados,distribucion)
#Gráfico
grafico = ggplot(data=datos,aes(x=pac_muestreados,y=distribucion))
grafico = grafico + geom_bar(stat="identity",fill="lightblue3")
grafico = grafico + theme_bw() + ggtitle("Distribución de probabilidades")
grafico = grafico + xlab("Prob. de obtener 1 enfermo al seleccionar de 1 a 30,000 pacientes") +
ylab("Probabilidad")
ggplotly(grafico)Ahora veámoslo con 3 pacientes con tuberculosis
#===============
# Respuesta - Pregunta 3 - Parte C
#===============
#Datos
pac_tub = 3
pac_muestreados=seq(0,100000) #Pacientes muestrados
distribucion = dhyper(x=pac_tub, m=A, k=pac_muestreados, n=B)
datos=data.frame(pac_muestreados,distribucion)
#Gráfico
#Gráfico
grafico = ggplot(data=datos,aes(x=pac_muestreados,y=distribucion))
grafico = grafico + geom_bar(stat="identity",fill="lightblue3")
grafico = grafico + theme_bw() + ggtitle("Distribución de probabilidades")
grafico = grafico + xlab("Prob. de obtener 3 enfermos al seleccionar de 1 a 30,000 pacientes") +
ylab("Probabilidad")
ggplotly(grafico)El periodo de incubación del virus Sars-CoV-2 se puede aproximar a una distribución normal de media 7 y desviación estándar 2 (Paul and Lorin 2021). Dado esto:
En este caso la respuesta es 0, ya que se trata de una distribución continua. Sin embargo, se puede efectuar una aproximación usando la función dnorm(), que trabaja con cuantiles (vér ayuda de la función).
#===============
# Respuesta - Pregunta 3 - Parte A
#===============
x=4.5
mu=7
sd=2
respuesta=dnorm(x,mu,sd)En este caso, la probabilidad sería 0.0913245.
En este caso tenemos una probabilidad binomial considerando \(x\) éxitos (20) para \(n\) (100) ensayos de Bernoullí. El éxito corresponde a tener un periodo de incubación viral menor a 7 días. Esto se traduce en:
#===============
# Respuesta - Pregunta 4 - Parte B
#===============
x=7
mu=7
sd=2
p=pnorm(x,mu,sd) #Probabilidad p(x<7)
#Distribución binomial
x_exitos = 20
n_ensayos = 100
respuesta = dbinom(x_exitos, size = n_ensayos,prob = p)En este caso, la probabilidad de que entre 100 pacientes, 20 tengan una incubación menor a 7 días es 4.2281633^{-10}.
La distribución normal inicial se puede graficar de la siguiente manera:
#===============
# Respuesta - Pregunta 4 - Parte C
#===============
T_incubacion = seq(0,15,by=0.2)
distribucion = dnorm(T_incubacion, mean=7,sd = 2)
datos=data.frame(T_incubacion,distribucion)
#Gráfico
library("ggplot2")
grafico = ggplot(data=datos,aes(x=T_incubacion,y=distribucion))
grafico = grafico + geom_bar(stat="identity",fill="lightblue3")
grafico = grafico + theme_bw() + ggtitle("Distribución de probabilidades")
grafico = grafico + xlab("Periodo de incubación") + ylab("Probabilidad")
ggplotly(grafico)En cambio, la distribución binomial se puede graficar como:
#===============
# Respuesta - Pregunta 4 - Parte C
#===============
exitos = seq(0,100)
distribucion = dbinom(exitos, size = n_ensayos,prob = p)
datos=data.frame(exitos,distribucion)
#Gráfico
library("ggplot2")
grafico = ggplot(data=datos,aes(x=exitos,y=distribucion))
grafico = grafico + geom_bar(stat="identity",fill="lightblue3")
grafico = grafico + theme_bw() + ggtitle("Distribución de probabilidades")
grafico = grafico + xlab("Éxitos") + ylab("Probabilidad")
ggplotly(grafico)Como la probabilidad base de que el periodo de incubación sea menor a 7 días es 0.5, la probabilidad de que de 100 casos, 20 cumplan con la condición de éxito es casi 0.