from scipy.stats import norm
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from __future__ import print_function
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
import squarify    # pip install squarify (algorithm for treemap)
from matplotlib.gridspec import GridSpec

def plot_tree_maps(P_vn, P_fp, P_fn, P_vp, plot1, plot2, plot3, plot4, plot5):

  def stringfy(num):
    if (num>1000000):
      return f"{num/1000000:.1f}MM"
    if (num>1000):
      return f"{num/1000:.0f}M"
    return f"{int(num)}"

  VN = int(P_vn * 200_000_000)
  FP = int(P_fp * 200_000_000)
  FN = int(P_fn * 200_000_000)
  VP = int(P_vp * 200_000_000)

  df_square = pd.DataFrame({
    'size':[VN,FP,VP,FN],
    'label':[f"VN\n{stringfy(VN)}", f"FP\n{stringfy(FP)}", f"VP\n{stringfy(VP)}", f"FN\n{stringfy(FN)}"],
    # 'color':["green","darkorange","indianred", "darkred"],
    'color':["green","blue","orange", "red"],
    'type':['VN','FP','VP','FN']
  })

  alpha = .5


  positive_test_mask = ((df_square['type']=='VP') | (df_square['type']=='FP'))
  negative_test_mask = (~positive_test_mask)
  diseased_mask = ((df_square['type']=='VP') | (df_square['type']=='FN'))
  not_diseased_mask = (~diseased_mask)

  # Todos os grupos
  plot1.axis('off')
  plot1.set_title('População toda')
  squarify.plot(sizes=df_square['size'], label=df_square['label'], color=df_square['color'], alpha=alpha, ax=plot1)

  # Apenas testes positivos
  plot2.axis('off')
  plot2.set_title('Testes positivos')
  df_temp = df_square[positive_test_mask]
  squarify.plot(sizes=df_temp['size'], label=df_temp['label'], color=df_temp['color'], alpha=alpha, ax=plot2)

  # Apenas testes negativos
  plot3.axis('off')
  plot3.set_title('Testes negativos')
  df_temp = df_square[negative_test_mask]
  squarify.plot(sizes=df_temp['size'], label=df_temp['label'], color=df_temp['color'], alpha=alpha, ax=plot3)

  # Apenas pessoas doentes
  plot4.axis('off')
  plot4.set_title('Doentes')
  df_temp = df_square[diseased_mask]
  squarify.plot(sizes=df_temp['size'], label=df_temp['label'], color=df_temp['color'], alpha=alpha, ax=plot4)

  # Apenas pessoas saudáveis
  plot5.axis('off')
  plot5.set_title('Saudáveis')
  df_temp = df_square[not_diseased_mask]
  squarify.plot(sizes=df_temp['size'], label=df_temp['label'], color=df_temp['color'], alpha=alpha, ax=plot5)

# formata float como % com 2 casas decimais
def percentify(val):
    return f'{100*val:.2f}%'

# Discretização do eixo X
xs_glob = list(np.linspace(0,12,100))

def plot_pdfs(mu_doente, sigma_doente, mu_saldavel, sigma_saldavel, prop, threshold):

    # Truque pra usar menos pontos e mesmo assim hachurar a área bonitinho
    xs = pd.Series([threshold-0.001,threshold+0.001]+xs_glob).sort_values()

    threshold_index = np.where(xs>threshold)[0][0]
    fig = plt.figure(figsize=(22,11))
    gs = GridSpec(ncols=5, nrows=4, hspace=.4, height_ratios=[2,1,1,2])
    plot = fig.add_subplot(gs[0:2,:4])
    matrix = fig.add_subplot(gs[0,4:])
    coefs = fig.add_subplot(gs[2,:4])
    roc_plot = fig.add_subplot(gs[1:3,4:])

    plot1 = fig.add_subplot(gs[3,0])
    plot2 = fig.add_subplot(gs[3,1])
    plot3 = fig.add_subplot(gs[3,2])
    plot4 = fig.add_subplot(gs[3,3])
    plot5 = fig.add_subplot(gs[3,4])
    
    # valores da PDF para os DOENTES
    y1s = norm.pdf(xs, loc=mu_doente, scale=sigma_doente)*prop
    plot.plot(xs, y1s)
    # hachurar os falsos negativos
    error1_vals = np.array((list([0.0]*threshold_index) + list(y1s[threshold_index:])))
    plot.fill_between(xs, error1_vals, alpha=0.5)
    
    # valores da PDF para os NÃO DOENTES
    y2s = norm.pdf(xs, loc=mu_saldavel, scale=sigma_saldavel)*(1-prop)
    plot.plot(xs, y2s)
    # hachurar os falsos positivos
    error2_vals = np.array((list(y2s[:threshold_index]) + list([0.0]*(len(xs)-threshold_index))))
    plot.fill_between(xs, error2_vals, alpha=0.5)
    
    # linha vertical de limiar/threshold
    plot.axvline(threshold, 0, 1, color='black', alpha=200)
    
    # legendas
    plot.legend(labels=[
        "Distribuição doentes",
        "Distribuição saldáveis",
        "Threshold ou Limiar [positivos | negativos]",
        "Prob. falsos negativos",
        "Prob. falsos positivos",
    ])

    # fig.title("Distribuições das respostas ao teste")
    
    # Cálculo das probabilidades
    P_vp = norm.cdf(xs[threshold_index], loc=mu_doente, scale=sigma_doente)*prop      # verdadeiro positivo
    P_fn = norm.sf(xs[threshold_index], loc=mu_doente, scale=sigma_doente)*prop       # falso positivo
    P_vn = norm.sf(xs[threshold_index], loc=mu_saldavel, scale=sigma_saldavel)*(1-prop)   # verdadeiro negativo
    P_fp = norm.cdf(xs[threshold_index], loc=mu_saldavel, scale=sigma_saldavel)*(1-prop)  # falso negativo
    
    # Criando matriz de confusão
    df_conf = pd.DataFrame()
    df_conf['Teste negativo'] = [P_vn,P_fn]
    df_conf['Teste positivo'] = [P_fp,P_vp]
    df_conf['idx'] = ['Não doente','Doente']
    df_conf.set_index('idx',inplace=True)

    sns.heatmap(df_conf, annot=True, ax=matrix, cmap="Blues")

    # Construção da matriz de confusão
    d = [
        {'Real':'Positivo', 'Teste positivo':percentify(P_vp), 'Teste negativo':percentify(P_fn), '  % da pop.':percentify(P_vp+P_fn)},
        {'Real':'Negativo', 'Teste positivo':percentify(P_fp), 'Teste negativo':percentify(P_vn), '  % da pop.':percentify(P_fp+P_vn)},
        {'Real':'% da pop.', 'Teste positivo':percentify(P_vp+P_fp), 'Teste negativo':percentify(P_fn+P_vn), '  % da pop.':percentify(P_vp+P_fn+P_fp+P_vn)}
    ]
    print(pd.DataFrame(d).set_index('Real'))
    print('-'*100)
    
    # Cálculo dos indicadores de qualidade do teste
    sensitivity = P_vp/(P_vp+P_fn)    # prob. de classificar corretamente alguém como positivo
    specificity = P_vn/(P_vn+P_fp)    # prob. de classificar corretamente alguém como negativo
    PPV = P_vp/(P_vp+P_fp)            # prob. de ser positivo dado que o teste deu positivo
    NPV = P_vn/(P_vn+P_fn)            # prob. de ser negativo dado que o teste deu negativo
    print(f"Precision/PPV: {100*PPV:.2f}%   NPV: {100*NPV:.2f}%   sensitivity: {100*sensitivity:.2f}%   specificity: {100*specificity:.2f}%")


    precisao = P_vp/(P_vp+P_fp)
    sensibilidade = P_vp/(P_vp+P_fn)
    especificidade = P_vn/(P_vn+P_fp)
    TFP = P_fp/(P_vn+P_fp)
    df_coefs = pd.DataFrame()
    df_coefs['precisao'] = [precisao]
    df_coefs['especificidade'] = [especificidade]
    df_coefs['sensibilidade / TVP'] = [sensibilidade]
    df_coefs['TFP'] = [TFP]
    sns.heatmap(df_coefs, annot=True, ax=coefs, cmap="Blues")


    df_roc = pd.DataFrame()
    df_roc['thresholds'] = np.linspace(-3,13,100)
    df_roc['P_vp'] = norm.cdf(df_roc['thresholds'], loc=mu_doente, scale=sigma_doente)*prop      # verdadeiro positivo
    df_roc['P_fn'] = norm.sf(df_roc['thresholds'], loc=mu_doente, scale=sigma_doente)*prop       # falso positivo
    df_roc['P_vn'] = norm.sf(df_roc['thresholds'], loc=mu_saldavel, scale=sigma_saldavel)*(1-prop)   # verdadeiro negativo
    df_roc['P_fp'] = norm.cdf(df_roc['thresholds'], loc=mu_saldavel, scale=sigma_saldavel)*(1-prop)  # falso negativo
    df_roc['TVP'] = df_roc['P_vp']/(df_roc['P_vp']+df_roc['P_fn'])
    df_roc['TFP'] = df_roc['P_fp']/(df_roc['P_vn']+df_roc['P_fp'])
    
    # sns.lineplot(x=df_roc['TFP'],y=df_roc['TVP'], ax=roc_plot)
    roc_plot.plot(df_roc['TFP'],df_roc['TVP'])
    roc_plot.set_title('ROC')

    # Labels
    roc_plot.set_xlabel("TFP")
    roc_plot.set_ylabel("TVP")

    # Limites
    lim_inf = -0.1
    lim_sup = 1.1
    roc_plot.set_xlim(lim_inf, lim_sup)
    roc_plot.set_ylim(lim_inf, lim_sup)

    # tracejados da caixa
    if(lim_inf<0 and lim_sup>1):
        roc_plot.plot([0, 0], [lim_inf, lim_sup], linestyle="--", alpha=0.5, color='#888888')
        roc_plot.plot([1, 1], [lim_inf, lim_sup], linestyle="--", alpha=0.5, color='#888888')
        roc_plot.plot([lim_inf, lim_sup], [0, 0], linestyle="--", alpha=0.5, color='#888888')
        roc_plot.plot([lim_inf, lim_sup], [1, 1], linestyle="--", alpha=0.5, color='#888888')

    # ponto
    x_p = df_coefs['TFP'][0]
    y_p = df_coefs['sensibilidade / TVP'][0]
    roc_plot.scatter(x_p, y_p)

    # ticks 
    ticks = []
    roc_plot.set_xticks(ticks+[df_coefs['TFP'][0]])
    roc_plot.set_yticks(ticks+[df_coefs['sensibilidade / TVP'][0]])

    # linhas guia
    roc_plot.plot([x_p, x_p], [lim_inf, y_p], linestyle="--", alpha=0.75, color='#000000')
    roc_plot.plot([lim_inf, x_p], [y_p, y_p], linestyle="--", alpha=0.75, color='#000000')



    # Treemaps
    plot_tree_maps(P_vn, P_fp, P_fn, P_vp, plot1, plot2, plot3, plot4, plot5)

Imagine que você precisa criar um teste para uma doença...

Notebook interativo pra testar o efeito de diversos fatores como:

• Proporção entre pessas doentes e não doentes: prop
  • ex.: se a doença afeta 10% da população, deixar em 0.1
• Médias dos resultados dos testes: mu_doente, mu_saldavel
• Desvios dos resultados dos testes: sigma_doente, sigma_saldavel
• Limiar para resultado do teste: threshold

Exemplos de efeitos observados:

• Separar as médias ou reduzir os desvios melhora a qualidade dos testes e faz a ROC tender ao ideal
• Aproximar as médias ou aumentar os desvios piora os testes e aproxima a ROC da linha diagonal
• Diminuir o threshold leva o ponto sob a ROC em direção ao extremo (0,0) e aumentar leva em direção ao extremo (1,1)
• O valor de threshold controla sua relação de falsos negativos e falsos positivos
• bla, bla, bla...

# Configurando sliders interativos
interact(
    plot_pdfs,
    mu_doente=widgets.IntSlider(min=0, max=10, step=1, value=4),
    sigma_doente=widgets.FloatSlider(min=0.1, max=5, step=.1, value=1.0),
    mu_saldavel=widgets.IntSlider(min=0, max=10, step=1, value=6),
    sigma_saldavel=widgets.FloatSlider(min=0.1, max=5, step=.1, value=1.0),
    prop=widgets.FloatSlider(min=0.01, max=0.99, step=.01, value=0.5),
    threshold=widgets.FloatSlider(min=0, max=10, step=.25, value=5),
);

• Precision / PPV (Positive Predictive Value) is the percentage of patients with a positive test who actually have the disease $$Precision = PPV = \frac{TP}{TP + FP}$$
• NPV (Negative Predictive Value) is the percentage of patients with a negative test who do not have the disease $$NPV = \frac{TN}{FN + TN}$$
• Specificity is the ability of a test to correctly classify an individual as disease-free $$Specificity = \frac{TN}{TN + FP}$$
• Sensitivity / TPR (True Positive Rate) is the ability of a test to correctly classify an individual as ′diseased′ $$Sensitivity = TPR = \frac{TP}{TP + FN} = [TVP (Taxa de Verdadeiros Positivos)]$$
• FPR (False Positive Rate) / False alarm ratio $$FPR = \frac{FP}{FP + TN} = [TFP (Taxa de Falsos Positivos)]$$

  fonte: Understanding and using sensitivity, specificity and predictive values

Obs:

• Mudar a proporção entre DOENTES e NÃO DOENTES afeta PPV e NPV, mas não afeta Specificity e Sensitivity
• Mesmo deixando tudo simétrico, pode ser que os valores fiquem com uma pequena diferença causada pela discretização grosseira do eixo X

---

Minha dúvida: O balanço (prop) não tem nenhum efeito sobre ROC mesmo?