Как предсказать риск ишемической болезни сердца с помощью логистической регрессии?

Болезни сердца являются основной причиной смерти во всем мире, на них приходится одна треть смертей в 2019 году. Число случаев сердечных заболеваний за этот период почти удвоилось: с 271 миллиона в 1990 году до >523 млн в 2019 году, а число смертей от сердечно-сосудистых заболеваний выросло с 12,1 млн до 18,6 млн.

Ишемическая болезнь сердца (ИБС) включает снижение притока крови к сердечной мышце из-за образования бляшек в артериях сердца.

Одной из самых сложных задач является выявление причин этого заболевания и его профилактика по мере возможности. Медицинские диагностические рассуждения становятся популярным применением машинного обучения, в котором экспертные системы и схемы на основе моделей предоставляют механизмы, используемые для разработки гипотез, которые затем будут проверены с использованием методов моделирования и моделирования или статистического анализа.

Прогнозирование сердечно-сосудистых заболеваний считается одной из самых важных тем в области здравоохранения. Алгоритмы машинного обучения и наличие больших объемов данных позволяют экстраполировать информацию, которая может помочь врачам делать более точные прогнозы.

Прогнозирование ИБС является очень сложной задачей, учитывая уровень опыта и знаний, необходимых для точного результата. Согласно опросу ВОЗ,медицинские работники могут правильно прогнозировать сердечные заболевания с точностью 67%.

В этой статье ряд независимых переменных, таких как пол, возраст, cigsPerDay, totChol, sysBP и глюкоза, будут использоваться вместе с зависимой переменной (TenYearCHD класс) на этапе обучения для построения модели классификации. Цель классификации состоит в том, чтобы предсказать, есть ли у пациента 10-летний риск развития сердечно-сосудистых заболеваний в будущем (ИБС) или нет.

Хорошо, приступим!

Шаг 1. Импорт необходимых библиотек

import pandas as pd
import numpy as np
import statsmodels.api as sm
import scipy.stats as st
import matplotlib.pyplot as plt
import seaborn as sn
from sklearn.metrics import confusion_matrix
import matplotlib.mlab as mlab
%matplotlib inline

Шаг 2. Подготовка данных

Набор данных взят с веб-сайта Kaggle (набор данных исследования Framingham Heart). Набор данных предоставляет информацию о пациентах. Он включает более 4000 записей и 15 атрибутов, как указано ниже.

df = pd.read_csv('framingham.csv')
df.drop(['education'],axis=1,inplace=True)
df.rename(columns={'male':'sex_male'},inplace=True)
df.head()

Входные переменные:

пол : мужской или женский (номинальный)
age : возраст пациента (непрерывно)
currentSmoker: является ли пациент курильщиком в настоящее время (номинальное значение).
cigsPerDay : количество сигарет, выкуриваемых человеком в среднем за один день (непрерывно).
BPMeds : принимал ли пациент лекарства от артериального давления (номинальное значение).
prevalentStroke: был ли у пациента ранее инсульт (номинальное значение).
prevalentHyp: была ли у пациента артериальная гипертензия (номинальное значение).
диабет : был ли у пациента диабет (номинальное значение).
totChol: уровень общего холестерина (непрерывно)
sysBP : систолическое кровяное давление (непрерывно)
diaBP : диастолическое кровяное давление (непрерывно)
ИМТ: индекс массы тела (непрерывный)
heartRate : частота сердечных сокращений (непрерывно).
глюкоза: уровень глюкозы (непрерывно)
TenYearCHD : 10-летний риск ишемической болезни сердца (ИБС) (двоичный код: 1 (да), 0 (нет))

Шаг 3. Работа с отсутствующими значениями

df.isnull().sum()

count = 0
for i in df.isnull().sum(axis=1):
    if i > 0:
        count = count+1
print('Total number of rows with missing values is', count)

Из приведенного выше вывода мы получили общее количество строк с отсутствующим значением 489. В этом случае, поскольку это всего лишь 12 % от всего набора данных, строки с отсутствующими значениями исключаются.

Мы получаем данные из процесса исключения, как показано ниже:

df.dropna(axis=0, inplace=True)
df.info()

Шаг 4. Исследовательский визуальный анализ

Цель этого шага — увидеть распределение всех данных.

def draw_histograms(dataframe, features, rows, cols):
    fig=plt.figure(figsize=(20,20))
    for i, feature in enumerate(features):
        ax=fig.add_subplot(rows,cols,i+1)
        dataframe[feature].hist(bins=20,ax=ax,facecolor='red')
        ax.set_title(feature+"Distribution", color='blue')
    fig.tight_layout()
    plt.show()
draw_histograms(df, df.columns, 6, 3)

Из приведенной выше гистограммы мы можем лучше понять распределение данных.

Следующим шагом является определение того, сколько пациентов подвержены риску развития ишемической болезни сердца или нет.

sn.countplot(x='TenYearCHD', data=df)

Из приведенного выше рисунка можно сделать вывод о наличии 3179пациентов с отсутствием ишемической болезни сердца и 572 пациентов с риском ишемической болезни сердца. Болезнь сердца.

Шаг 5. Давайте моделировать!

Логистическая регрессия – это алгоритм классификации машинного обучения, который используется для прогнозирования вероятности категориальной зависимой переменной. Это расширение модели линейной регрессии для задач классификации. В отличие от линейной регрессии, которая выводит непрерывные числовые значения, логистическая регрессия преобразует свои выходные данные с помощью функции логистической сигмоиды, чтобы вернуть значение вероятности, которое затем может быть сопоставлено с двумя или более дискретные классы.

Первый шаг, который мы собираемся сделать с точки зрения процесса моделирования, — добавить постоянное значение в очищенный набор данных.

from statsmodels.tools import add_constant as add_constant
df_constant = add_constant(df)
df_constant.head()

И на следующем этапе мы пытаемся увидеть параметр и результат нашей модели логистической регрессии:

st.chisqprob = lambda chisq, df: st.chi2.sf(chisq, df)
cols = df_constant.columns[:-1]
model = sm.Logit(df.TenYearCHD, df_constant[cols])
result = model.fit()
result.summary()

Приведенный выше результат показывает некоторые атрибуты со значением P выше предпочтительного альфа (5%), что предполагает, что атрибуты имеют низкую статистически значимую связь с вероятностью ишемической болезни сердца (ИБС).

Чтобы решить эту проблему, мы используем метод обратного исключения для удаления атрибутов с самым высоким значением P по одному за раз, после чего регрессия выполняется повторно, пока все атрибуты не получат P -значения менее 0,05.

Выбор функций: обратное исключение (подход P-value)

Зависимая переменная и список имен столбцов многократно запускают регрессию, удаляя признаки с P-значением выше альфа (5%) по одному и возвращая сводку регрессии с все p-значения ниже альфа (5%).

def back_feature_elem (data_frame, dep_var, col_list):
    while len(col_list)>0 :
        model = sm.Logit(dep_var,data_frame[col_list])
        result = model.fit(disp=0)
        largest_pvalue = round(result.pvalues,3).nlargest(1)
        if largest_pvalue[0]<(0.05):
            return result
            break
        else:
            col_list = col_list.drop(largest_pvalue.index)
result = back_feature_elem(df_constant, df.TenYearCHD, cols)
result.summary()

В конце мы получаем атрибуты, p-значения которых меньше альфа (5%):

секс
возраст
сигарет в день
тотхол
системный БП
глюкоза

и мы получаем это уравнение логистической регрессии:

Шаг 6. Интерпретация: отношение шансов, доверительные интервалы и P-значения

Следующим шагом мы собираемся интерпретировать отношение шансов, доверительные интервалы и P-значения, используя этот код:

params = np.exp(result.params)
conf = np.exp(result.conf_int())
conf['OR'] = params
pvalue = round(result.pvalues,3)
conf['pvalue'] = pvalue
conf.columns = ['CI 95%(2.5%)','CI 95%(97.5%)', 'Odds Ratio', 'pvalue']
print((conf))

Интерпретация результатов выше:

Вероятность получения диагноза ишемической болезни сердца (ИБС)для мужчин (пол = 1) выше, чем у женщин(sex=0) равно exp(0,5815) = 1,788687. С точки зрения процентного изменения, мы можем сказать, что шансы для мужчин получить диагноз 78,8% выше, чем шансы для женщин.
Коэффициент возраста говорит, что если мы хотим посмотреть на вероятность диагноза ИБС, взглянув на параметр увеличения возраста на один год, он будет примерно равен exp(0,0655) = 1,067644. или около 7%.
Каждая выкуренная дополнительная сигарета увеличивает 2%(exp(0,0197)= 1,019895) вероятность постановки диагноза CDH.
Коэффициент увеличивается на 1,7% (exp(0,0174)=1,017552) на каждую единицу увеличения систолического артериального давления.
Существует 0,7% увеличение (exp(0,0076)=1,007628) шансов поставить диагноз CDH на каждую единицу увеличения уровня глюкозы.

Шаг 7. Разделение данных: обучение и тестирование

import sklearn
new_features =  df[['age','sex','cigsPerDay','totChol','sysBP','glucose','TenYearCHD']]
x = new_features.iloc[:,:-1]
y = new_features.iloc[:,-1]
from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x,y,test_size=.20, random_state=5)
from sklearn.linear_model import LogisticRegression
logreg = LogisticRegression()
logreg.fit(x_train, y_train)
y_pred = logreg.predict(x_test)

Шаг 8. Оценка модели

Оценка модели направлена на оценку точности обобщения модели на будущих (невидимых/вне выборки) данных.

1. Точность

Точность описывает, насколько точно модель может правильно классифицировать. Это количество правильных прогнозов, сделанных как отношение всех сделанных прогнозов. Мы используем модуль sklearn для вычисления точности задачи классификации.

sklearn.metrics.accuracy_score(y_test,y_pred)

Исходя из этого вывода, мы можем интерпретировать, составляет ли точность полученной модели около 87,5 % в проверочном наборе.

2. Матрица путаницы

Матрица путаницы определяет количество TP (истинно положительный), TN (истинно отрицательный), FP (ложноположительный), FN (False Negative) в прогнозах классификатора.

Из матрицы путаницы мы можем получить точность, которая определяется как сумма скорректированных прогнозов, деленная на общее количество прогнозов:

Точность = TP+TN/TP+FP+FN+TN

Также другие стандартные показатели производительности:

Чувствительность или Отзыв = TP / TP + FN или True-Positive Rate (TPR)
Специфичность = TN / TN + FP или Истинный отрицательный показатель (TNR)
Точность = TP / TP + FPили Положительное прогнозируемое значение
Коэффициент ложноположительных результатов (FPR) = FP / FP + TN
Коэффициент ложноотрицательных результатов (FNR) = FN / FN + TP
Оценка F1 = 2*(напоминание * точность) / (напоминание + точность)

Для хороших классификаторов TPR и TNR должны быть ближе к 100 %. То же самое и с параметрами precision и accuracy. Напротив, FPR и FNR должны быть как можно ближе к 0%.

from sklearn.metrics import confusion_matrix
cm = confusion_matrix(y_test, y_pred)
conf_matrix = pd.DataFrame(data=cm, columns=['Predicted:0','Predicted:1'],index=['Actual:0','Actual:1'])
plt.figure(figsize = (8,5))
sn.heatmap(conf_matrix, annot=True, fmt='d', cmap='YlGnBu')

Приведенная выше матрица путаницы показывает 652 + 5 = 657 правильных предсказаний и 87 + 5 = 92 неправильных.

Истинно положительные результаты (TP): 5
Истинно отрицательные результаты (TN):652
Ложные срабатывания (FP): 7 (ошибка типа I)
Ложноотрицательные результаты (FN): 87 (ошибка типа II)

TN = cm[0,0]
TP = cm[1,1]
FN = cm[1,0]
FP = cm[0,1]
sensitivity = TP/float(TP+FN)
specificity = TN/float(TN+FP)
print('The acuuracy of the model = TP+TN/(TP+TN+FP+FN) = ',(TP+TN)/float(TP+TN+FP+FN),'\n',
'The Missclassification = 1-Accuracy = ',1-((TP+TN)/float(TP+TN+FP+FN)),'\n',
'Sensitivity or True Positive Rate = TP/(TP+FN) = ',TP/float(TP+FN),'\n',
'Specificity or True Negative Rate = TN/(TN+FP) = ',TN/float(TN+FP),'\n',
'Positive Predictive value = TP/(TP+FP) = ',TP/float(TP+FP),'\n',
'Negative Predictive Value = TN/(TN+FN) = ',TN/float(TN+FN),'\n',
'Positive Likelihood Ratio = Sensitivity/(1-Specificity) = ',sensitivity/(1-specificity),'\n',
'Negative Likelihood Ratio = (1-Sensitivity)/Specificity = ',(1-sensitivity)/specificity)

Точность модели = 87,48 %
Неправильная классификация = 12,52 %
Чувствительность или доля истинно положительных результатов (TPR) =54,35 %
Специфичность или доля истинно отрицательных результатов (TNR) =98,94 %
Точность или положительное прогнозируемое значение =41,67 %
Отрицательное прогнозируемое значение = 88,23 %
Положительный коэффициент правдоподобия = 5,12
Отрицательный коэффициент правдоподобия = 0,96

Из приведенной выше статистики видно, что модель является высокоспецифичной, чем чувствительной. Отрицательные значения предсказываются более точночем положительные.

Шаг 9. Прогнозируемые вероятности

0 (ишемическая болезнь сердца: нет) и 1 (ишемическая болезнь сердца: да) для тестовых данных с порогом классификации по умолчанию 0,5.

y_pred_prob = logreg.predict_proba(x_test)[:,:]
y_pred_prob_df = pd.DataFrame(data=y_pred_prob, columns=['Prob of no hearts disease(0)','Prob of Heart Disease (1)'])
y_pred_prob_df.head()

Нижний порог

Из матрицы путаницы мы можем сделать вывод, что если количество ложноотрицательных результатов (FN) (ошибка типа II) очень велико, его можно классифицировать как умеренно опасный, потому что, это означает игнорирование вероятности заболевания, когда на самом деле существует одно ИЛИ Верно.

Следовательно, чтобы повысить чувствительность, порог можно снизить.

from sklearn.preprocessing import binarize
for i in range(1,5):
    cm2=0
    y_pred_prob_yes=logreg.predict_proba(x_test)
    y_pred2=binarize(y_pred_prob_yes,i/10)[:,1]
    cm2=confusion_matrix(y_test,y_pred2)
    print ('With',i/10,'threshold the Confusion Matrix is ','\n',cm2,'\n',
            'with',cm2[0,0]+cm2[1,1],'correct predictions and',cm2[1,0],'Type II errors( False Negatives)','\n\n',
          'Sensitivity: ',cm2[1,1]/(float(cm2[1,1]+cm2[1,0])),'Specificity: ',cm2[0,0]/(float(cm2[0,0]+cm2[0,1])),'\n\n\n')

3. ROC-кривая

Площадь под ROC-кривой — это показатель производительности для измерения способности бинарного классификатора различать положительные и отрицательные классы. ROC — это кривая вероятности, а AUC представляет собой меру разделимости. Он говорит, насколько модель способна различать классы.

from sklearn.metrics import roc_curve
fpr, tpr, thresholds = roc_curve(y_test, y_pred_prob_yes[:,1])
plt.plot(fpr,tpr)
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.0])
plt.title('ROC curve for Heart disease classifier')
plt.xlabel('False positive rate (1-Specificity)')
plt.ylabel('True positive rate (Sensitivity)')
plt.grid(True)

Площадь под кривой (AUC)

Площадь под кривой ROC количественно определяет точность классификации модели; чем выше площадь, тем больше разница между истинными и ложными положительными результатами и тем сильнее модель классифицирует элементы обучающего набора данных. Площадь 0,5 соответствует модели, которая работает не лучше, чем случайная классификация, и хороший классификатор остается как можно дальше от нее. Площадь 1 идеальна. Чем ближе AUC к 1, тем лучше.

sklearn.metrics.roc_auc_score(y_test,y_pred_prob_yes[:,1])

Вывод:

Все атрибуты, выбранные после процесса исключения, показывают P-значения ниже 5%, что позволяет предположить, что выбранные атрибуты играют важную роль в прогнозировании ишемической болезни сердца (ИБС).
Мужчины, по-видимому, более подвержены сердечным заболеваниям, чем женщины. Увеличение возраста, количества сигарет, выкуриваемых в день, и систолического артериального давления также указывают на увеличение вероятности сердечно-сосудистых заболеваний.
Модель предсказала с точностью 87,5%. Модель более специфична, чем чувствительна.
Площадь под кривой ROC составляет 73,5, что в некоторой степени удовлетворительно.
Общая модель может быть улучшена за счет большего количества данных.