====== 8.e Intervalos de confianza ======
Cuando calculamos una métrica realmente no tenemos el valor real de la métrica sino un valor (llamado estimador) que nos ayuda para saber el verdadero valor de la métrica.


El valor real de la métrica no lo podemos obtener pero si al menos in intervalo en el que es probable que esté ese valor real. A ese intervalo se le llama el intervalo de confianza.

En este prospecto médico para un test de COVID se muestran los intervalos de confianza la para la sensibilidad y la especificidad: {{ :clase:iabd:pia:2eval:sars-cov-2_rapid_antigen_test_es.pdf |}}.

Los datos que se muestran son los siguientes:

$$
\begin{array}
\\
E&=&115&=&Nº \; de \; enfermos
\\
S&=&311&=&Nº \; de \; sanos
\\
TP&=&111&=&Nº \; predichos \; positivos
\\
TN&=&310&=&Nº \; predichos \; negativos
\end{array}
$$

$$
\begin{array}
\\
Sensibilidad&=&\frac{111}{115}&=&0,9652&=&96,52\%
\\
Especificidad&=&\frac{310}{311}&=&0,9968&=&99,68\%
\end{array}
$$

La sensibilidad real y la especificidad real no la conocemos ya que con otros paciente habría dado unos resultados distintos, por ello se calcula el intervalo de confianza que nos dice que ese intervalo con un 95% de probabilidad tendrá dentro el valor buscado.

Siguiendo con el test de covid los intervalos de confianza mostrados son:

$$
\begin{array}
\\
Sensibilidad&=&[91,33\% \; ‑ \; 99,04\%]
\\
Especificidad&=&[98,22\% \; ‑ \; 99,99\%]
\end{array}
$$

Es decir que:
  * El intervalo de confianza de $[91,33\% \; ‑ \; 99,04\%]$ tiene una probabilidad del 95% de contener al valor real de la **Sensibilidad**
  * El intervalo de confianza de $[98,22\% \; ‑ \; 99,99\%]$ tiene una probabilidad del 95% de contener al valor real de la **Especificidad**


===== Cálculo del intervalo de confianza =====
Hay varias formas de calcular el intervalo de confianza, en el ejemplo de COVID se usa ''Clopper-Pearson''.

Para calcular esos intervalos se necesitan los siguientes datos:
  * Tamaño de la muestra o número de intentos: Realmente es el denominador de nuestra métrica
  * Número de éxitos o de positivos: Realmente es el numerador de nuestra métrica
  * Nivel de confianza: Suele ser el 95% pero podría ser cualquier valor , aunque si se cambia debería ser un valor mayor aun para estar más seguro.

La forma de calcular los intervalos de confianza es:
  * Mediante la siguiente página web: [[https://epitools.ausvet.com.au/ciproportion|Calcular límites de confianza para una proporción muestral]]
  * Mediante la función ''proportion_confint'' de python: [[https://www.statsmodels.org/v0.12.2/generated/statsmodels.stats.proportion.proportion_confint.html|statsmodels.stats.proportion.proportion_confint]]


Desgraciadamente hay varios formas de calcular el intervalo de confianza, los vamos a nombrar ahora y tambien como se //llaman// en la función de ''proportion_confint''.

^ Nombre ^ Nombre en ''proportion_confint'' ^
| Aproximación normal (Wald) |  ''normal''  |
| Clopper-Pearson exacto basado en la distribución Beta |  ''beta''  |
| Wilson |  ''wilson''  |
| Jeffreys con pior Beta(0.5, 0.5)  |  ''jeffreys''  |
| Agresti-Coull |  ''agresti_coull''  |

Veamos un ejemplo en Python que calcula todos los intervalos de confianza para todos los tipos:

<sxh python>
from statsmodels.stats.proportion import proportion_confint

# Número de éxitos
success = 111
# Número total de ensayos
total = 115
# Nivel de confianza deseado
confianza=95

methods=['normal','beta','wilson','jeffreys','agresti_coull']
for method in methods:
    lower, upper = proportion_confint(success, total, alpha=1-(confianza/100), method=method)

    print(f"Intervalo de confianza al {confianza}% con el método \"{method}\"= [{lower:.4f}, {upper:.4f}]")
</sxh>

Y el resultado es:
<sxh base>
Intervalo de confianza al 95% con el método "normal"= [0.9317, 0.9987]
Intervalo de confianza al 95% con el método "beta"= [0.9133, 0.9904]
Intervalo de confianza al 95% con el método "wilson"= [0.9140, 0.9864]
Intervalo de confianza al 95% con el método "jeffreys"= [0.9194, 0.9881]
Intervalo de confianza al 95% con el método "agresti_coull"= [0.9111, 0.9893]
</sxh>

Y podemos observar que el cálculo del intervalo del test de COVID se hizo con ''Clopper-Pearson''

===== Elección del método =====
De forma general y sin querer entrar en detalles el método que deberemos usar en las métricas será el Jeffreys (''jeffreys'')

¿Porque Jeffreys? Porque es el que usa estadística bayesiana para calcular el intervalo de confianza. Así que realmente no es un intervalo de confianza sino un intervalo de credibilidad. ¿cual es la diferencia? Vamos a explicarlo con el ejemplo que estamos usando: 
  * **Intervalo de confianza**  de ''Clopper-Pearson'' de ''[0.9133, 0.9904]'' significa que hay un 95% de probabilidad de que ese intervalo contenga al valor real que estamos buscando.
  * **Intervalo de credibilidad**  de ''Jeffreys'' de ''[0.9194, 0.9881]'' significa que hay un 95% de probabilidad de que el valor real que estamos buscando esté en ese intervalo. Es decir, que el valor real es una variable aleatorio y cada valor posible tiene una probabilidad. Y el intervalo de credibilidad es el conjunto de valores reales cuya probabilidad es el 95%. 

En el siguiente código se puede ver como da el mismo resultado usando ''proportion_confint(method="jeffreys")'' que calculándolo mediante el teorema de bayes y usando como prior una $Beta(0.5,0.5)$. Y por lo tanto ''Jeffreys'' es un **Intervalo de credibilidad** ya que es lo que se calcula en la estadística bayesiana.

**NOTA:La explicación del cálculo bayesiano va más allá del ámbito de este curso de especialización.**

<sxh python>
from scipy.stats import beta
from statsmodels.stats.proportion import proportion_confint

success=111
total=115
confidence=0.95

lower,upper=proportion_confint(success,total,method="jeffreys",alpha=1-confidence)
print(f"Intervalo de confianza usando \"Jeffreys\"\t\t\t=[{lower:.6f} - {upper:.6f}]")
beta_a=0.5+success
beta_b=0.5+total-success
lower,upper=beta.interval(confidence,beta_a,beta_b)
print(f"Intervalo de confianza usando Bayes con prior Beta(0.5,0.5)\t=[{lower:.6f} - {upper:.6f}]")

</sxh>

<sxh base>
Intervalo de confianza usando "Jeffreys"			=[0.919431 - 0.988148]
Intervalo de confianza usando Bayes con prior Beta(0.5,0.5)	=[0.919431 - 0.988148]
</sxh>

¿Y cuando se podrían usar otros? El resto no son bayesianos así que dan "peor" información pero en lineas generales, podemos decir que:
  * Clopper-Pearson o Jeffreys: Tamaño de muestra pequeño
  * Wilson o Agresti-Coull: Tamaño de muestra grande y proporción no extrema
  * Normal (Wald): No es recomendable



===== Comparación de métricas =====

La mayor utilidad del cálculo del intervalo de confianza es la comparación de las métricas entre distintos modelos y así saber si los modelos realmente tienen las mismas métricas.

Para ello seguiremos la siguiente regla:
  * Dos modelos tienen la **misma** métrica si sus intervalos **se solapan** en algún punto
  * Dos modelos tienen **distinta** métrica si sus intervalos **no se solapan** en ningún punto
  * En caso de **"empate"**, se elige el modelo con la menor variabilidad

Veamos un ejemplo.

Supongamos que tenemos la sensibilidad de 5 modelos y queremos saber si realmente son iguales o no. Los intervalos de confianza de la sensibilidad de los 4 modelos son:

$$
\begin{array}
\\
Modelo \; 1&=&[0.9678,0.9856]
\\
Modelo \; 2&=&[0.9464,0.9789]
\\
Modelo \; 3&=&[0.8567,0.9126]
\\
Modelo \; 4&=&[0.7986,0.8934]
\\
Modelo \; 5&=&[0.9245,0.9403]
\\
\end{array}
$$

Podemos comparar los valores pero es mejor mostrarlos en una gráfica para ayudarnos en el trabajo.

{{:clase:iabd:pia:2eval:comparacion_intervalos_confianza.png?direct|}}

Se puede apreciar que:
  * Los modelos 1 y 2 son iguales aunque es mejor el 1 ya que tiene menor variabilidad.
  * Los modelos 3 y 4 son iguales aunque es mejor el 3 ya que tiene menor variabilidad.

Por lo que el orden de los modelos de mejor a peor es: **1, 2, 5, 3 y 4**


La siguiente función en Python crea un gráfico para poder comparar intervalos de confianza:
<sxh python>
def plot_intervals(axes,intervalos,title,xlabel,ylabels):
    for index, (inferior, superior) in enumerate(intervalos):
        position=len(intervalos)-index-1
        media=((superior-inferior)/2)+inferior
        axes.plot([inferior, superior], [position, position], marker='|', markersize=10, linestyle='-', color='#003B80')
        axes.scatter(media, position, color='#003B80')        

    axes.set_yticks(np.arange(len(ylabels)))
    axes.set_yticklabels(ylabels[::-1])
    axes.set_xlabel(xlabel)
    axes.set_title(title)

    axes.set_facecolor("#F0F7FF")
    axes.grid(visible=True, which='major', axis='both',color="#FFFFFF",linewidth=1)

</sxh>


Y para crear la gráfica de nuestro ejemplo se llamaría así
<sxh python>
figure=plt.figure()
axes = figure.add_subplot()

intervalos=[
    [0.9678,0.9856],
    [0.9464,0.9789],
    [0.8567,0.9126],
    [0.7986,0.8934],
    [0.9245,0.9403]
]
ylabels=['Modelo 1', 'Modelo 2', 'Modelo 3', 'Modelo 4', 'Modelo 5']

plot_intervals(axes,intervalos,"Comparación de modelos","Sensibilidad",ylabels)

</sxh> 


Otra forma de hacer el gráfico es poner lineas en los extremos de los intervalos, que aunque se ve mejor, queda un gráfico muy recargado.

<sxh python>
def plot_intervals_with_guides(axes,intervalos,title,xlabel,ylabels):
    for index, (inferior, superior) in enumerate(intervalos):
        position=len(intervalos)-index-1
        media=((superior-inferior)/2)+inferior
        color=next(axes._get_lines.prop_cycler)['color']
        axes.plot([inferior, superior], [position, position], marker='|', markersize=10, linestyle='-', color=color)
        axes.scatter(media, position, color=color) 
        axes.plot([inferior, inferior], [0, len(intervalos)-1], linewidth=0.5, linestyle='dashed', color=color)   
        axes.plot([superior, superior], [0, len(intervalos)-1], linewidth=0.5, linestyle='dashed', color=color)
        
    axes.set_yticks(np.arange(len(ylabels)))
    axes.set_yticklabels(ylabels[::-1])
    axes.set_xlabel(xlabel)
    axes.set_title(title)

    axes.set_facecolor("#F0F7FF")
    axes.grid(visible=True, which='major', axis='both',color="#FFFFFF",linewidth=1)


</sxh>

{{:clase:iabd:pia:2eval:comparacion_intervalos_confianza_lineas.png?direct|}}

===== Comparación avanzada de métricas =====
En caso de que los intervalos se solapen, no significa exactamente que tengan la misma métrica. Para saber si realmente la tienen es necesario restar las distribuciones y ver si el cero (o algo cercano al cero) está en dicha resta.

Para hacer los cálculos vamos a crear una serie de funciones en python.

<sxh python>
def get_hdi(sample, hdi_prob):
    """
            Calculation of 100(1-a)% HDI for sample
            Chen-Shao HPD Estimation algorithm
            http://y-okamoto-psy1949.la.coocan.jp/Python/en/calc_map_est/
    """
    smpl = np.sort(sample)
    n = len(smpl)
    ckW = max(smpl) - min(smpl)
    cki = 0
    RB = int(n * (hdi_prob))
    L0 = 0
    U0 = n-1
    while (RB + cki < n):
        if smpl[RB + cki] - smpl[cki] < ckW:
            ckW = smpl[RB + cki] - smpl[cki]
            L0 = cki
            U0 = RB + cki
        cki += 1
    ckW = max(smpl) - min(smpl)
    cki = 1
    LB = n - int(n * (hdi_prob))
    L1 = 0
    U1 = n-1
    while (LB - cki >= 0):
        if smpl[n - cki] - smpl[LB - cki] < ckW:
            ckW = smpl[n - cki] - smpl[LB - cki]
            U1 = n - cki
            L1 = LB - cki
        cki += 1
    return [smpl[int((L1+L0)/2)], smpl[int((U1+U0)/2)]]



def get_interval(success,total,confidence=0.95):
    beta_a=0.5+success
    beta_b=0.5+total-success
    actual_beta=beta(beta_a,beta_b)
    hdi=get_hdi(actual_beta.rvs(100000),confidence)
    return hdi





def get_hdi_compare(success1,total1,success2,total2,confidence=0.95):
    beta_a1=0.5+success1
    beta_b1=0.5+total1-success1
    beta1=beta(beta_a1,beta_b1)
    beta_a2=0.5+success2
    beta_b2=0.5+total1-success2
    beta2=beta(beta_a2,beta_b2)

    n_samples = 100000
    datos1=beta1.rvs(n_samples)
    datos2=beta2.rvs(n_samples)

    resultado=datos1-datos2
    
    lower,upper=get_hdi(resultado,confidence)

    return [lower,upper]
</sxh>

Lo primero es que hemos creado la función ''get_interval'' que nos retorna el [[https://easystats.github.io/bayestestR/articles/credible_interval.html|intervalo de credibilidad]] (en este caso el Highest Density Interval  o HDI ). Calculamos ambos intervalos y los mostramos:

<sxh python>
TP1=131
FN1=16

TP2=124
FN2=23


interval1=get_interval(TP1,TP1+FN1)
interval2=get_interval(TP2,TP2+FN2)


figure=plt.figure(figsize=(5.5, 3))
axes = figure.add_subplot()

intervalos=[
    interval1,
    interval2
]
ylabels=['Modelo 1', 'Modelo 2']
 
plot_intervals(axes,intervalos,"Comparación de modelos","Sensibilidad",ylabels)
</sxh>

La gráfica resultante es: 
{{:clase:iabd:pia:2eval:comparacion_2_metricas.png?direct|}}


Y para saber si realmente tienen la misma métrica simplemente ejecutamos el siguiente código:
<sxh python>
print(get_hdi_compare(TP1,TP1+FN1,TP2,TP2+FN2))
</sxh>
Y el resultado es:
<sxh base>
[-0.029769860656823077, 0.12498680076438595]
</sxh>

**Si el intervalo incluye al cero es que tienen la misma métrica.** Y en el ejemplo como el intervalo contiene el cero es que tienen la misma métrica.

Veamos ahora otro ejemplo:
<sxh python>
TP1=131
FN1=16

TP2=116
FN2=31


interval1=get_interval(TP1,TP1+FN1)
interval2=get_interval(TP2,TP2+FN2)


figure=plt.figure(figsize=(5.5, 3))
axes = figure.add_subplot()

intervalos=[
    interval1,
    interval2
]
ylabels=['Modelo 1', 'Modelo 2']
 
plot_intervals_with_guides(axes,intervalos,"Comparación de modelos","Sensibilidad",ylabels)
</sxh>
{{:clase:iabd:pia:2eval:comparacion_2_metricas_b.png?direct|}}

<sxh python>
print(get_hdi_compare(TP1,TP1+FN1,TP2,TP2+FN2))
</sxh>

<sxh base>
[0.018992455698845223, 0.18505496723415582]
</sxh>

Por lo que en este caso se solapan pero no tiene la misma métrica. 

De todas en este caso, el intervalo está rozando el cero, hay un concepto llamado [[https://easystats.github.io/bayestestR/articles/region_of_practical_equivalence.html|Region of Practical Equivalence (ROPE)]] que dice que aunque no esté el cero en el intervalo, nos vale una región cercana al cero ya que de forma práctica es como si fuera el cero.

Mas información:
  * {{ :clase:iabd:pia:2eval:a_bayesian_interpretation_of_the_confusion_matrix.pdf |A bayesian interpretation of the confusion matrix}}
  * {{ :clase:iabd:pia:2eval:classifier_uncertainty.evidence_potential_impact_and_probabilistic_treatment.pdf |Classifier uncertainty:evidence, potential impact and probabilistic treatment}}
  * {{ :clase:iabd:pia:2eval:the_balanced_accuracy_and_its_posterior_distribution.pdf |The balanced accuracy and its posterior distribution}}
  * [[https://doingbayesiandataanalysis.blogspot.com/2018/05/just-published-rejecting-or-accepting.html|Just published: "Rejecting or Accepting Parameter Values in Bayesian Estimation"]]
  * [[https://doingbayesiandataanalysis.blogspot.com/2013/08/how-much-of-bayesian-posterior.html|How much of a Bayesian posterior distribution falls inside a region of practical equivalence (ROPE)]]

===== Ejercicios =====

==== Ejercicio 1 ====
Mirando el test de COVID de {{ :clase:iabd:pia:2eval:sars-cov-2_rapid_antigen_test_es.pdf |}}. Usando la web [[https://epitools.ausvet.com.au/ciproportion|Calcular límites de confianza para una proporción muestral]] calcula el intervalo de confianza de la especificidad y comprueba que coincide con el del prospecto.

==== Ejercicio 2 ====
Mirando el test de COVID de {{ :clase:iabd:pia:2eval:sars-cov-2_rapid_antigen_test_es.pdf |}}. Calcula en Python el intervalo de confianza de la especificidad y comprueba que coincide con el del prospecto.

==== Ejercicio 3 ====
Crear 5 modelos distintos para el problema del cancer y para cada uno de ellos:
  * Calcula las métricas de sensibilidad y especificidad.
  * Calcula el intervalo de confianza de las métricas de sensibilidad y especificidad. ¿Que método se va a usar?

Tambien:
  * Muestra en una gráfica el intervalo de confianza de la métrica de sensibilidad de todos los modelos
  * Muestra en una gráfica el intervalo de confianza de la métrica de especificidad de todos los modelos

Ahora:
  * Indica según la sensibilidad el orden de los bueno que es cada modelo (de mejor a peor).
  * Indica según la especificidad el orden de los bueno que es cada modelo (de mejor a peor)
  * ¿Que modelo elegirías?

{{:clase:iabd:pia:2eval:comparacion_modelos_metricas_intervalo.png?direct|}}


==== Ejercicio 4 ====
Elige el mejor de los 5 modelos anteriores y muestra ahora una gráfica, en la que se mostrará:

  * El valor de la sensibilidad y la especificidad según el valor del umbral. La gráfica tendrá un aspecto simular a la siguiente:

{{:clase:iabd:pia:2eval:threshold-metricas.png?direct|}}
¿Que umbral elegirías?



  * Repite ahora la gráfica, pero incluyendo los intervalos de credibilidad de cada métrica. La gráfica tendrá un aspecto simular a la siguiente:

{{:clase:iabd:pia:2eval:threshold-metricas-intervalo-confianza.png?direct|}}
¿Que umbral elegirías?

<note tip>Para hacer este típo de gráficas usa el método:

''axes.fill_between(rango_threshold, rango_sensibilidad_lower, rango_sensibilidad_upper, color=color, alpha=0.5)''</note>