<jupyter_start><jupyter_text>4 - Otimização de modelosConcluída a fase de modelagem, vimos que dois modelos, regressão logística e XGBoost, preveem adequadamente os resultados de testes de COVID. Vamos agora tentar otimizá-los.Nessa fase, buscaremos os melhores parâmetros para utilizar com estes modelos, a fim de determinar o modelo final a ser utilizado para predição de novos dados.<jupyter_code># para carregar a base de dadas limpa
import pickle

import numpy as np
import pandas as pd
from scipy import stats as spst

from pprint import pprint

# gráficos
import seaborn as sns
from matplotlib import rcParams, pyplot as plt

# parâmetros do matplotlib
# essencialmente, para deixar os gráficos maiores por padrão
rcParams['figure.dpi'] = 120
rcParams['figure.figsize'] = (10, 8)

# warnings
import warnings
warnings.filterwarnings("ignore")

# pacote com funções para análise desse projeto
import os
cwd = os.getcwd()
os.chdir("../")
import scripts.plots as splt, scripts.metrics as smetrics
os.chdir(cwd)<jupyter_output><empty_output><jupyter_text>Importação dos modelos fitadosVamos importar os modelos já fitados porém com parâmetros padrão.<jupyter_code>with open(r'../models/modelo_default.model', 'rb') as modelfile:
    pickler = pickle.Unpickler(file = modelfile)
    modelos_import = pickler.load()
modelos = modelos_import['modelo']

X_train, X_test, y_train, y_test = modelos_import['train_test_split']
covid = modelos_import['base']

covid.info()
covid.head()<jupyter_output><empty_output><jupyter_text>Otimização dos modelosCada modelo de aprendizado de máquina é construído em cima de alguns parâmetros. Por exemplo, no modelo de regressão logística, há um parâmetro regulando a presença ou não do intercepto (`fit_intercept`) e um outro controlando a regularização do modelo, *i.e.* a penalização para coeficientes muito diferentes (`penalty`), entre outros. Estes são chamados ***hiperparâmetros***.A otimização do modelo consiste na escolha dos hiperparâmetros de forma a maximizar/minimizar alguma métrica. Neste caso, para sermos consistentes com o estudo realizado na parte de modelagem, vamos testar várias combinações de hiperparâmetros e ranqueá-las através do *Brier score*. Um modelo se adere mais ao resultado se essa métrica for menor que a de um outro modelo; logo, iremos minimizar o *Brier score* das combinações de hiperparâmetros.<jupyter_code>nome_score = 'roc_auc'  # por que neg? Porque as classes do scikit-learn trabalham com maximização. 
                                # Como a lógica do Brier score é de minimização, basta colocar um sinal negativo para inverter a lógica<jupyter_output><empty_output><jupyter_text>*Cross validation*O ***cross validation*** (validação cruzada) é uma técnica tipicamente utilizada para medir o **poder de generalização** do nosso modelo. São feitos testes na prática para saber o quão preciso é um modelo, ou seja o seu desempenho em um novo conjunto de dados. De início, faremos a separação da mesma forma em treino e teste, em seguida usaremos a seguinte técnica de generalização:O *cross validation* pode ser usado para encontrar hiperparâmetros. Para cada combinação de hiperparâmetros, realizamos o ajuste do modelo utilizando-os, e então rodamos o *cross validation* para verificar se o modelo tem o poder de generalização necessário.  Busca de hiperparâmetros: `RandomizedSearchCV`E como se determinam as opções de hiperparâmetros? Para cada hiperparâmetro, especificamos as alternativas, que passamos para um buscador denominado `RandomizedSearchCV`.  Implementação do algoritmo de otimizaçãoHá dois modelos a serem otimizados e comparados:* Regressão Logística* XGBoostPara cada modelo, há um conjunto distinto de parâmetros a serem tentados.<jupyter_code>modelos.keys()
params = {}

params['logistic'] = {
    'penalty': ['elasticnet', 'l1', 'l2'],
    'fit_intercept': [True, False],
    'solver': [ 'lbfgs', 'saga' ],
}

params['xgboost'] = {
    'criterion': ['friedman_mse', 'mae' ],
    'max_features': [ 'log2', 'sqrt' ],
    'learning_rate': [ 0.01, 0.05, 0.1, 0.5, 1, ],
    'max_depth': [ 3, 4, 5, 10 ],
    'n_estimators': [ 5 ,10, 15, 20, 100, 150 ],
}

params2 = {}
for nome_modelo, params_dict in params.items():
    params2[nome_modelo] = { f'{nome_modelo}__{k}': v for k, v in params_dict.items() }

params = params2
params
from sklearn.model_selection import StratifiedKFold, RandomizedSearchCV
from datetime import datetime as dt

searchers = {}

cv = StratifiedKFold(
    n_splits = 3,
)

# deletar a variável t0 caso exista, para que a lógica de timing dos testes não se confunda
try:
    del t0
except NameError:
    pass

for nome_modelo, modelo in modelos.items():

    params_grid = params[nome_modelo]
    
    try:
        tfim = dt.now()
        delta = tfim - t0
        print(f"... OK. Duração da busca: {delta})\n")
    except NameError:   # variavel t0 não existe: primeira iteração do loop
        pass

    print(f"Procurando hiperparâmetros ótimos para o modelo '{nome_modelo}'... ")
    t0 = dt.now()
    random_search = RandomizedSearchCV(
        estimator = modelo, 
        param_distributions = params_grid, 
        scoring = nome_score,
        n_iter = 150,
        cv = cv, 
        verbose = 2,
        n_jobs = -1,
    )
    
    searchers[nome_modelo] = random_search.fit(X_train, y_train)

tfim = dt.now()
delta = tfim - t0
print(f"... OK. Duração da busca: {delta}")<jupyter_output>Procurando hiperparâmetros ótimos para o modelo 'logistic'... 
Fitting 3 folds for each of 12 candidates, totalling 36 fits
... OK. Duração da busca: 0:00:17.782146)

Procurando hiperparâmetros ótimos para o modelo 'xgboost'... 
Fitting 3 folds for each of 150 candidates, totalling 450 fits
[20:59:34] WARNING: C:/Users/Administrator/workspace/xgboost-win64_release_1.5.1/src/learner.cc:576: 
Parameters: { "criterion", "max_features" } might not be used.

  This could be a false alarm, with some parameters getting used by language bindings but
  then being mistakenly passed down to XGBoost core, or some parameter actually being used
  but getting flagged wrongly here. Please open an issue if you find any such cases.


... OK. Duração da busca: 0:04:22.084204<jupyter_text>Avaliação dos modelos otimizadosVamos agora avaliar as otimizações feitas.  Melhores parâmetros encontrados<jupyter_code>from pprint import pprint

for nome_modelo, searcher in searchers.items():
    print(f"Melhores parâmetros para modelo '{nome_modelo}':")
    pprint(searcher.best_params_)
    print('')<jupyter_output>Melhores parâmetros para modelo 'logistic':
{'logistic__fit_intercept': False,
 'logistic__penalty': 'l2',
 'logistic__solver': 'saga'}

Melhores parâmetros para modelo 'xgboost':
{'xgboost__criterion': 'friedman_mse',
 'xgboost__learning_rate': 0.1,
 'xgboost__max_depth': 3,
 'xgboost__max_features': 'sqrt',
 'xgboost__n_estimators': 100}<jupyter_text>Resultado da otimização<jupyter_code>cvresults = {}

for nome_modelo, searcher in searchers.items():
    rsdf = pd.DataFrame(searcher.cv_results_)
    rsdf = rsdf.set_index([ col for col in rsdf.columns if col.startswith('param_') ])
    cvresults[nome_modelo] = rsdf.sort_values(by = 'rank_test_score')
for nome_modelo, rsdf in cvresults.items():
    print(f"Resultados da busca de hiperparâmetros do modelo '{nome_modelo}':")
    display(rsdf)<jupyter_output>Resultados da busca de hiperparâmetros do modelo 'logistic':<jupyter_text>A melhora no desempenho de cada modelo é:<jupyter_code>from sklearn.metrics import SCORERS

fscore = SCORERS[nome_score]

print(f"Usando score = '{nome_score}'\n")

best_scores = { nome_modelo: searcher.best_score_ for nome_modelo, seacher in searchers.items() }

for nome_modelo, modelo in modelos.items():

    print(f"Resultados para modelo '{nome_modelo}':")
    score_original = fscore(
        estimator = modelos[nome_modelo],
        X = X_test,
        y_true = y_test,
    )
    score_otimizado = cvresults[nome_modelo].iloc[0]['mean_test_score']

    if nome_score.startswith('neg'):
        melhoria = score_original / score_otimizado - 1
    else:
        melhoria = score_otimizado / score_original - 1

    frase_score_orig = f"Score '{nome_score}' original"
    frase_score_opt = f"Score '{nome_score}' otimizado"
    frase_melhoria = 'Melhoria'

    ncaract = len(frase_score_opt) + 4

    print(f"{frase_score_orig:>{ncaract}s}: {score_original:}")
    print(f"{frase_score_opt:>{ncaract}s}: {score_otimizado:}")
    print(f"{frase_melhoria:>{ncaract}s}: {melhoria:.3%}")
    print('')<jupyter_output>Usando score = 'roc_auc'

Resultados para modelo 'logistic':
     Score 'roc_auc' original: 0.6527829092680578
    Score 'roc_auc' otimizado: 0.6469761085350275
                     Melhoria: -0.890%

Resultados para modelo 'xgboost':
     Score 'roc_auc' original: 0.6489462964321899
    Score 'roc_auc' otimizado: 0.6536565903083834
                     Melhoria: 0.726%<jupyter_text>O modelo XGBoost teve desempenho pior com parâmetros *default*, mas com parâmetros otimizados parece ser o melhor modelo.<jupyter_code>from copy import deepcopy

# vamos fazer um outro dicionario de modelos
opt_modelos = {}

# para cada item no dicionario de modelos original...

for nome_modelo, modelo in modelos.items():
    # obtem os hiperparametros otimizados
    opt_params = searchers[nome_modelo].best_params_

    # faz uma cópia do modelo original, para não dar problemas de passagem por referencia
    opt_modelo = deepcopy(modelo)

    # seta os hiperparametros
    opt_modelo.set_params(**opt_params)

    # retreina com os novos parâmetros
    opt_modelo.fit(
        X_train, y_train
    )

    # e armazena o novo modelo
    opt_modelos[f'{nome_modelo}_opt'] = opt_modelo<jupyter_output>[21:04:16] WARNING: C:/Users/Administrator/workspace/xgboost-win64_release_1.5.1/src/learner.cc:576: 
Parameters: { "criterion", "max_features" } might not be used.

  This could be a false alarm, with some parameters getting used by language bindings but
  then being mistakenly passed down to XGBoost core, or some parameter actually being used
  but getting flagged wrongly here. Please open an issue if you find any such cases.


[21:04:16] WARNING: C:/Users/Administrator/workspace/xgboost-win64_release_1.5.1/src/learner.cc:576: 
Parameters: { "criterion", "max_features" } might not be used.

  This could be a false alarm, with some parameters getting used by language bindings but
  then being mistakenly passed down to XGBoost core, or some parameter actually being used
  but getting flagged wrongly here. Please open an issue if you find any such cases.<jupyter_text>Vamos ainda plotar a curva ROC e calcular a área sob a mesma (métrica AUROC).<jupyter_code>modelos_auroc = {}
modelos_auroc.update(modelos)
modelos_auroc.update(opt_modelos)

splt.plot_roc_auc_curve(
    models = modelos_auroc,
    X_test = X_test, y_test = y_test
)<jupyter_output><empty_output><jupyter_text>Vemos que o modelo XGBoost parece ser o que mais acuradamente representa as probabilidades previstas.  Otimização de *threshold* de decisãoApós a análise das probabilidades que o modelo prediz, resta ainda uma pergunta: dada uma probabilidade de o teste ser positivo, qual deve ser o limiar a partir do qual a predição deve ser "teste positivo"?Para tal, há que se fazer uma análise do "custo" de ocorrer um falso positivo vs um falso negativo.<jupyter_code>modelo_final = opt_modelos['xgboost_opt']

y_prob = modelo_final.predict_proba(X_test)[:, 1]<jupyter_output><empty_output><jupyter_text>Neste caso, **é mais prejudicial um falso negativo**, *i.e.* prevermos que um paciente receberá um teste negativo quando ele estiver infectado. A consequência para a comunidade é altamente negativa, pois esse paciente pode infectar outras pessoas sem se dar conta que está potencialmente acometido por COVID.No caso de um falso positivo (prevermos que o paciente receberá um teste positivo quando ele não estiver infectado), o paciente incorre no prejuízo de se isolar desnecessariamente.<jupyter_code>def calc_threshold_optimization(y_true, y_prob, beta = 1, calcs = dict()):

    # generate threshold dataframe
    thresholds = pd.DataFrame(np.linspace(0, 1, 101), columns = ['threshold'])

    # calc tn, fp, fn, tp for each threshold
    threshold_conf_mat = thresholds.apply(
            smetrics.calc_confusion_matrix, 
            y_true = y_true, y_prob = y_prob, 
            axis = 1
        )

    threshold_opt = pd.concat([
        thresholds,
        threshold_conf_mat
    ], axis = 1)

    # calc a few other metrics
    threshold_opt['precision'] = (
        (threshold_opt.tp / (threshold_opt.tp + threshold_opt.fp))
            .where(threshold_opt.tp + threshold_opt.fp > 0, np.nan)
    )

    threshold_opt['recall'] = (
        (threshold_opt.tp / (threshold_opt.tp + threshold_opt.fn))
            .where(threshold_opt.tp + threshold_opt.fn > 0, np.nan)
    )

    threshold_opt['accuracy'] = (
        (threshold_opt.tp + threshold_opt.tn) /
        (threshold_opt.tn + threshold_opt.fp + threshold_opt.fn + threshold_opt.tp)
    )

    threshold_opt[f'f{beta:n}_score'] = (
        (
            (1+beta**2) * threshold_opt.precision * threshold_opt.recall /
            (threshold_opt.precision * beta**2 + threshold_opt.recall)
        )
            .where(
                threshold_opt.precision * beta**2 + threshold_opt.recall > 0,
                np.nan
            )
    )

    threshold_opt['cohen_kappa_score'] = np.vectorize(
        lambda threshold: smetrics.calc_cohen_kappa_score(
            y_true = y_true,
            y_prob = y_prob,
            threshold = threshold)
    )(threshold_opt['threshold'])

    for name, fun_config in calcs.items():
        fcalc = fun_config['fcalc']
        kwds = fun_config['kwds']

        fapply = lambda row: fcalc(row, **kwds)

        threshold_opt[name] = threshold_opt.apply(fapply, axis = 1)

        try:
            if fun_config['normalize']:
                colvals = threshold_opt[name]
                threshold_opt[name] = (colvals - colvals.min()) / (colvals.max() - colvals.min())
        except KeyError:
            pass

    # return the dataframe (index is threshold level)
    return threshold_opt.set_index('threshold')
# chamar a função de análise de threshold
# inclui uma coluna com o custo
# custo de cada um [[custo_tn, custo_fp], [custo_fn, custo_tp]]
# caso a função de otimização seja maximização, basta dar o sinal negativo

calcs = {
    'cost': {
        'fcalc': smetrics.calc_confusion_matrix_cost,
        'kwds': {
            'cost_matrix': [[1, -1], [-5, 1]],    # custo de fn é 5 vezes mais alto que custo de fp
        },
        'normalize': True     # custo normalizado para ficar na escala de 0 a 1
    }
}

threshold_opt = calc_threshold_optimization(
    y_true = y_test, 
    y_prob = y_prob,
    calcs = calcs,
    beta = 1.5 # priorizacao do recall sobre o precision
)

threshold_opt
fig = splt.plot_precision_recall_curve(y_true = y_test, y_prob = y_prob)
fig = splt.plot_threshold_curve(
    y_true = y_test, y_prob = y_prob,
    threshold_opt = threshold_opt,
    metrics = {
        'f1.5_score': r'F$\beta$-score ($\beta$ = 1.5)',
        'cost': 'Est. Custo (sinal invertido)'
    },
    best_threshold = True,
    minimize = False
)<jupyter_output><empty_output><jupyter_text>O gráfico acima evidencia que o *threshold* $\theta$ padrão não é adequado para tomada de decisão.Segundo a estimativa de custo (em que tomamos por hipótese que o custo de falsos negativos é 5 vezes maior que o custo de falsos positivos), o limiar ideal em que o resultado do modelo passa de "teste negativo" para "teste positivo" é $\theta = 24\%$. Em outras palavras, se a probabilidade de um teste positivo for menor que 24%, predizemos que o teste dará negativo. Já se a probabilidade for maior que 24%, predizemos que o teste dará positivo. Calculamos também o ***F-beta score*** em função do limiar de decisão. O *F-beta score* é similar ao *F1 score*, ou seja é uma média harmônica entre a precisão e o *recall*. A diferença é que essa média é ponderada, com um fator $\dfrac{1}{\beta _2}$ para a precisão. Logo, se $\beta  1$, o *F-beta* passa a priorizar o *recall* (igualmente sem ignorar completamente a precisão).Como julgamos que os falsos negativos são mais custosos, devemos priorizar o *recall*. Com um valor de $\beta = 1.5$, o ponto em que *F-beta* é máximo corresponde ao *threshold* $\theta = 24\%$, concordando com a análise similar do custo estimado. Este resultado igual é devido ao patamar que aparece nas curvas, bem como na curva de ***Precision-Recall Threshold***. Com o limiar $0\% < \theta <  \sim 30\%$, a precisão, *recall*, *F-beta* e custo estimado variam pouco. Em outras palavras, temos evidências de que o limiar ótimo está entre $0\%$ e $30\%$, mas não há insumos fortes o suficientes para definirmos categoricamente o *threshold* otimizado. Se a aplicação desse modelo for insumo para tomada de decisão de médicos, é apropriado apenas citar as probabilidades de um teste positivo de COVID.No entanto, para aplicações em que o consumidor final é leigo, é mais prudente que a saída do modelo seja uma decisão. Na falta de informações mais definitivas, **estabelecemos o limiar ótimo em $\theta = 24\%$**.<jupyter_code>pprint({ k: v for k, v in modelo_final.get_params().items() if k.startswith('xgboost__') })
from sklearn.metrics import f1_score, accuracy_score, recall_score, precision_score

y_pred_threshold = (y_prob >= 0.24) * (1)

metricas_args = [y_test, y_pred_threshold]

metricas_threshold = {
    'f1': f1_score(*metricas_args),
    'precision': precision_score(*metricas_args),
    'recall': recall_score(*metricas_args),
    'accuracy': accuracy_score(*metricas_args)
}

pprint(metricas_threshold)<jupyter_output>{'accuracy': 0.4549,
 'f1': 0.6162348634187553,
 'precision': 0.4472433702927801,
 'recall': 0.9904945117121196}
{'fraction_non_alphanumeric': 0.09209276334171557, 'fraction_numerical': 0.018720312936574463, 'mean_word_length': 4.021324354657688, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '49259', 'license': 'MIT', 'n_tokens_mistral': 6422, 'n_tokens_neox': 6090, 'n_words': 1881, 'path': 'notebooks_models/4-otimizacao.ipynb', 'repo_name': 'flimao/modelos-covid', 'size': '17895'}