Medidas de rendimiento diagnóstico

5.1. Medidas de rendimiento diagnóstico#

En esta libreta veremos cómo calcular algunas métricas claves de rendimiento diagnóstico como la sensibilidad y la especificidad entre otras.

Truco

Podrás utilizar las funciones que utilizaremos en esta unidad en tus propios análisis ya que el objetivo es generar código reutilizable.

Antes de comenzar, te recomiendo que revises el tema en statpearls.

5.1.1. Librerías y datos#

Los datos que utilizaremos en esta libreta fueron simulados, el código para generarlos está al final de la libreta.

import pandas as pd
import scipy.stats as stats

url = 'https://raw.githubusercontent.com/chrisdewa/curso_python/refs/heads/main/bases/datos_dx.csv'
datos = pd.read_csv(url)
datos.head()
Estándar Test
0 0 0
1 0 0
2 1 1
3 0 0
4 0 0

Vemos que los datos están compuestos por dos columnas, «Estándar» y «Test», la primera columna representa el estado de la enfermedad de acuerdo al estándar de oro, la segunda, representa la predicción realizada por la prueba diagnóstica.

5.1.2. Tablas de contingencia#

En general, para todos los cálculos que realizaremos, asumimos que contamos con datos en una matriz de confusión de 1 grado de libertad, es decir, 2x2 y se estructura de la siguiente manera:

Estándar de Oro positivo

Estándar de Oro negativo

Test positivo

Verdaderos positivos

Falsos positivos

Test negativo

Falsos negativos

Verdaderos negativos

Los verdaderos positivos (VP) son aquellos sujetos positivos en el test y en el estándar de oro. Los falsos positivos (FP) son aquellos que fueron positivos en el test pero no en el estándar. Los falsos negativos (FN) fueron negativos en el test pero positivos en el estándar. Los verdaderos negativos (VN) fueron negativos en ambos casos.

Veamos cómo generar la tabla con python.

El modo más directo es pd.crosstab que produce una tabulación cruzada entre las columnas.

tab1 = pd.crosstab(
    index=datos['Test'],       # las filas
    columns=datos['Estándar'], # las columnas
)
tab1
Estándar 0 1
Test
0 72 3
1 5 20

Analiza la tabla generada, ¿notas algo fuera de lugar?. Si observas la tabla, está invertida, en el formato:

Estándar de Oro negativo

Estándar de Oro positivo

Test negativo

Verdaderos negativos

Falsos negativos

Test positivo

Falsos positivos

Verdaderos positivos

Esto es muy importante de entender, ya que si no lo contemplamos los cálculos pueden salir al revés. Para muchos cálculos en python, como el OR, no es importante el cambio, dados los valores de los encabezados y el índice, los datos están en las columnas correctas, solo es el orden de las columnas el que no está de forma habitual.

Primero, veamos cómo corregirlo

tab2 = pd.crosstab(
    index=datos['Test'],       # las filas
    columns=datos['Estándar'], # las columnas
).reindex(index=[1,0], columns=[1, 0])
tab2
Estándar 1 0
Test
1 20 5
0 3 72

Lo que hicimos fue reordenar la tabla por medio del método reindex.

Ahora bien, ¿esto importa?, para el ojo humano sí, claro, porque en el contexto científico se espera este orden, pero veamos qué pasa cuando queremos obtener datos de la tabla para realizar los cálculos.

Busquemos los VP:

tab1.loc[1, 1]

np.int64(20)

El atributo loc permite buscar en un dataframe un dato particular, la sintaxis es de indización, por eso los corchetes, similar a listas y diccionarios. El primer argumento es la fila o filas y el segundo es la columna o columnas.

Nota

Si queremos todas las columnas a partir de la 3a para la fila 17, el código sería:

dataframe.loc[17, 3:]

Eso quiere decir que buscamos la fila índice 1 en la columna 1. Es decir, «Estándar» igual a 1 y «Test» igual a 1; o Verdaderos positivos, que es lo mismo. Por lo tanto, el código funcionará exactamente igual si utilizamos la tabulación 1 (automática) o la tabulación 2 (reindexada).

tab1.loc[1, 1] == tab2.loc[1, 1]

np.True_

Por lo tanto, no es necesario reindexar las tabulaciones cruzadas a menos que se busque pasar la tabla a un documento para publicación.

5.1.2.1. Extracción de datos de la tabulación#

La siguiente función toma una tabulación cruzada y extrae sus valores VP, FP, FN, VN.

tab1
Estándar 0 1
Test
0 72 3
1 5 20 
def descomponer_tabulacion(tabulacion: pd.DataFrame):
    """
    Extrae los valores de la tabla de contingencia 2x2 para un test diagnóstico.

    Asume que la tabla está organizada con:
    - Filas: resultado del test (1 = positivo, 0 = negativo)
    - Columnas: estándar diagnóstico (1 = enfermedad presente, 0 = enfermedad ausente)

    Parámetros
    ----------
    tabulacion : pd.DataFrame
        Tabla de contingencia 2x2 con los resultados del test y el estándar diagnóstico.

    Retorna
    -------
    tuple
        Una tupla con cuatro enteros: (VP, FP, FN, VN)
        - VP: Verdaderos Positivos
        - FP: Falsos Positivos
        - FN: Falsos Negativos
        - VN: Verdaderos Negativos
    """
    vp = tabulacion.loc[1, 1] 
    fp = tabulacion.loc[1, 0] 
    fn = tabulacion.loc[0, 1] 
    vn = tabulacion.loc[0, 0]

    return vp, fp, fn, vn
# el print no es necesario, es solo para agilizar la compilación del libro
# en producción.
print('función de extraer datos definida')

función de extraer datos definida

Veamos cómo funciona y si realmente es igual para la tabulación 1 y para la 2.

print(descomponer_tabulacion(tab1))
print(descomponer_tabulacion(tab2))

(np.int64(20), np.int64(5), np.int64(3), np.int64(72))
(np.int64(20), np.int64(5), np.int64(3), np.int64(72))

Podemos corroborar que efectivamente el resultado es el mismo, independientemente de si reindexamos o no.

5.1.3. Rendimiento diagnóstico#

Las principales métricas, o las más utilizadas, se resumen en la siguiente tabla

Métrica

Fórmula

Interpretación breve

Sensibilidad

VP / (VP + FN)

Capacidad de detectar correctamente a los enfermos

Especificidad

VN / (VN + FP)

Capacidad de identificar correctamente a los sanos

VPP

VP / (VP + FP)

Probabilidad de que un positivo sea realmente enfermo

VPN

VN / (VN + FN)

Probabilidad de que un negativo sea realmente sano

LR+

Sensibilidad / (1 - Especificidad)

Cuánto más probable es un resultado positivo en un enfermo que en un sano

LR−

(1 - Sensibilidad) / Especificidad

Cuánto más probable es un resultado negativo en un enfermo que en un sano

Viendo las fórmulas y contando con los valores extraidos de la función que definimos, calcular los valores en python es trivial.

def sensibilidad(tabulacion):
    """
    Calcula la sensibilidad (verdaderos positivos / total de enfermos).

    Parámetro:
    - tabulacion: DataFrame 2x2 con los resultados del test y estándar diagnóstico.

    Retorna:
    - Sensibilidad como flotante.
    """
    vp, fp, fn, vn = descomponer_tabulacion(tabulacion)
    return vp / (vp + fn)

def especificidad(tabulacion):
    """
    Calcula la especificidad (verdaderos negativos / total de sanos).
    """
    vp, fp, fn, vn = descomponer_tabulacion(tabulacion)
    return vn / (vn + fp)

def vpp(tabulacion):
    """
    Calcula el valor predictivo positivo (probabilidad de enfermedad dado test positivo).
    """
    vp, fp, fn, vn = descomponer_tabulacion(tabulacion)
    return vp / (vp + fp)

def vpn(tabulacion):
    """
    Calcula el valor predictivo negativo (probabilidad de no tener la enfermedad dado test negativo).
    """
    vp, fp, fn, vn = descomponer_tabulacion(tabulacion)
    return vn / (vn + fn)

def verosimilitud_positiva(tabulacion):
    """
    Calcula la razón de verosimilitud positiva: sensibilidad / (1 - especificidad).
    """
    sens = sensibilidad(tabulacion)
    espe = especificidad(tabulacion)
    return sens / (1 - espe)

def verosimilitud_negativa(tabulacion):
    """
    Calcula la razón de verosimilitud negativa: (1 - sensibilidad) / especificidad.
    """
    sens = sensibilidad(tabulacion)
    espe = especificidad(tabulacion)
    return (1 - sens) / espe

print('Funciones de rendimiento definidas')

Funciones de rendimiento definidas

Ya puedes utilizar estas fórmulas en tu código, con facilidad, sin embargo, podemos empaquetar toda la lógica en una sola función que genere una tabla de resumen con solo la tabulación.

def rendimiento_diagnostico(tabulacion):
    """
    Calcula las principales métricas de rendimiento diagnóstico y las devuelve en un DataFrame.

    Cada métrica se coloca como una fila con su valor en la columna 'valor'.

    Parámetro:
    - tabulacion: DataFrame 2x2 con los resultados del test y el estándar diagnóstico.

    Retorna:
    - DataFrame con las métricas en el índice y sus valores en la columna 'valor'.
    """
    # creamos primero un diccionario con los resultados.
    valores = {
        'Sensibilidad': sensibilidad(tabulacion),
        'Especificidad': especificidad(tabulacion),
        'VPP': vpp(tabulacion),
        'VPN': vpn(tabulacion),
        'Verosimilitud positiva': verosimilitud_positiva(tabulacion),
        'Verosimilitud negativa': verosimilitud_negativa(tabulacion),
    }
    # convertimos el resultado en un dataframe
    resultado = (
        pd.DataFrame  # es un dataframe
        .from_dict(   # creado desde un diccionario (from_dict)
            valores,  # pasamos el diccionario de valores
            orient='index',   # esto hace que las llaves se pasen al índice
            columns=['Valor'] # esto define el nombre de la columna.
        )
        .round(3) # redondeamos 3 puntos decimales para que sea legible.
    )
    return resultado

print('Función definida')

Función definida

Veamos cómo funciona con nuestras tabulaciones

rendimiento_diagnostico(tab1)
Valor
Sensibilidad 0.870
Especificidad 0.935
VPP 0.800
VPN 0.960
Verosimilitud positiva 13.391
Verosimilitud negativa 0.139

5.1.4. Intervalos de confianza#

Para cualquier proporción es posible calcular el intervalo de confianza. Hagámoslo para la sensibilidad. Recordando \(Sensibilidad = VP/(VP+FN)\) Para calcular utilizaremos el método de Wilson, puede revisar el siguiente paper para conocer cómo funcionan los demás.

Erdoğan S, Gülhan OT. Alternative Confidence Interval Methods Used in the Diagnostic Accuracy Studies. Comput Math Methods Med. 2016;2016:7141050. doi: 10.1155/2016/7141050. Epub 2016 Jul 11. PMID: 27478491; PMCID: PMC4958484.

from statsmodels.stats.proportion import proportion_confint

sens = sensibilidad(tab1)
vp, fp, fn, vn = descomponer_tabulacion(tab1)
bajo, alto = proportion_confint(
    count=vp,   # número de "éxitos"
    nobs=vp+fn, # número total de intentos
    method='wilson', # Revisa la cita. Este método ofrece un buen balance en la estimación.
)
print(f'Sensibilidad: {sens:.2%} (IC 95% {bajo:.2%} - {alto:.2%})')

Sensibilidad: 86.96% (IC 95% 67.87% - 95.46%)

Como pudes ver, con base en el tamaño de muestra, el intervalo de confianza es elevado.

5.1.4.1. Ejercicio#

  1. Implementa el cálculo del intervalo de confianza para las demás estimaciones.

  2. Intenta incorporarlo dentro de la función rendimiento_diagnostico.

Con estas herramientas podemos analizar el rendimiento diagnóstico de pruebas no solo propias de nuestra línea investigativa sino también de estudios que publiquen este tipo de datos.

5.1.5. VPP Bayesianos#

Ahora, no podemos cerrar sin hablar de que los valores predictivos son directamente influenciados por la prevalencia de la enfermedad, por lo que se recomienda, en la mayoría de los casos, utilizar el cálculo ajustado a prevalencia (forma bayesiana).

Recordemos que el teorema de bayes se define de la siguiente forma:

\(P(Hipótesis|Evidencia) = \frac{P(Evidencia|Hipótesis)P(Hipótesis)}{P(Evidencia)}\)

En este contexto la evidencia es el resultado del test y la hipótesis es que la persona esté enferma (Valor predictivo positivo). Luego entonces:

\(P(Hipótesis|Evidencia) = VPP\)

\(P(Evidencia|Hipótesis) = Sensibilidad\)

\(P(Hipótesis) = Prevalencia\)

\(P(Evidencia) = Probabilidad de test positivo\)

La probabilidad de un test positivo está data en los escenarios donde el test es positivo y el paciente tiene la enfermedad por la probabilidad de tener la enfermedad y cuando la enfermedad no está presente pero el test es positivo por la probabilidad de no tener la enfermedad, dicho de otra forma:

\(P(Evidencia) = Sensibilidad * Prevalencia + (1-Especificidad) * (1-Prevalencia)\)

Implementemos este cálculo en python:

def vpp_bayes(tabulacion, prevalencia):
    """Calcula el valor predictivo positivo para la prevalencia"""
    sens = sensibilidad(tabulacion)
    esp = especificidad(tabulacion)
    numerador = sens * prevalencia
    denominador = numerador + (1-esp) * (1-prevalencia)
    return numerador / denominador

print('función vpp_bayes definida')

función vpp_bayes definida

vpp(tab1), vpp_bayes(tab1, 0.15)

(np.float64(0.8), np.float64(0.7026615969581749))

Nota

Como podemos ver para una prevalencia del 15% nuestro test tiene un rendimiendo medio con un VPP de solo ~70.3%, importante a tener en cuenta con enfermedades raras o estacionales.

5.1.6. Ejercicio#

  1. Implementa la función vpn_bayes.

  2. Busca algún artículo que reporte la sensibilidad y especificidad de un test y utiliza las funciones implementadas para corroborar sus datos.

  3. Intenta graficar las diferentes medidas de rendimiento diagnóstico utilizando matplotlib y/o seaborn.

5.1.7. Nota sobre los datos#

La base de datos utilizada en esta libreta fue generada con el siguiente código:

import pandas as pd
import scipy.stats as stats
n = 100     # tamaño de muestra
prev = 0.15 # prevalencia
sens = 0.87 # sensibilidad
esp = 0.91  # especificidad
# El gold standar define la presencia de la enfermedad con 100% de precisión
# por lo que utilizamos la distribución bernoulli
gold = stats.bernoulli.rvs(prev, size=n)
test = [
    stats.bernoulli.rvs(sens) # en caso positivo usamos la sensibilidad
    if subj == 1 else # comparación si el sujeto fue positivo
    stats.bernoulli.rvs(1-esp)# en caso negativo usamos 1 - especificidad
    for subj in gold # por cada sujeto en "gold" (estándar de oro)
]
data = pd.DataFrame({'Estándar': gold, 'Test': test}) # convertimos ambas variables a un dataset
data.to_csv('datos_dx.csv', index=False)

Puedes generar tus propios datos imulados con estos datos.

El código usa una comprension de lista para generar la variable test, puedes leer al respecto aquí.

Contenido