SAS ENTERPRISE GUIDE 6.1: ЛИНЕЙНАЯ РЕГРЕССИЯ.

Полный текст

Введение

При выполнении научных исследований может возникать задача предсказания значения какого-либо показателя в зависимости от концентрации того или иного лабораторного маркера [6, 8, 9, 10]. Одним из возможных подходов для решения такого рода задач, может быть построение линейной регрессии с использованием SAS Enterprise Guide 6.1 [4].

Цель работы состояла в разработке алгоритма прогнозирования повреждения миокарда на 3 день после коронарного шунтирования с искусственным кровообращением в зависимости от концентрации тропонина I на 1 день после операции.

Задачи включали разработку алгоритма оценки характера распределения исследуемых переменных, выбор адекватной прогностической модели, оценку точность построенной модели в программном пакете SAS Enterprise Guide 6.1.

Материалы и методы

В проспективное нерандомизированное исследование AMIRI — CABG в ПСПбГМУ им. И.П. Павлова с 2016 по 2019 годы включено 336 пациентов, имеющих показания к коронарному шунтированию. Пациенты распределены в три независимые группы наблюдения. Исходная тяжесть состояния пациентов была сопоставима по группам наблюдения, сравнение групп выполнялось с использованием SAS Enterprise Guide 6.1 [3]. Создана база данных, в которую внесены данные: концентрации тропонина I на 1 и 3 день после коронарного шунтирования. Для разработки алгоритма отобрана группа пациентов после коронарного шунтирования с искусственным кровообращением (n=128). Определен характер распределения исследуемых переменных, для построения прогностической модели выбрана линейная регрессия, оценена ее точность.

Результаты и их обсуждение

В программном комплексе SAS Enterprise Guide 6.1 реализован алгоритм линейной регрессии. Линейная регрессия представляет собой математическую модель линейной взаимосвязи двух переменных, позволяющую предсказывать значения одной переменной по значениям другой.

Ранее было показано, как загрузить базу данных из файла Excel в SAS [5].

На первом этапе с целью выбора построения прогностической модели была выполнена оценка характера распределения исследуемых переменных (TnI1 – тропонин I на 1-й день, TnI3 – тропонин I на 3-й день):

DATA TRIALDATA;

    SET WORK.'20_06_2019 work1'n; 

RUN;

TITLE 'Тест на нормальность';

PROC UNIVARIATE DATA=TRIALDATA normaltest;

    WHERE CPBType=2; /* Выбираем группу коронарного шунтирования c искусственным кровообращением */

    VAR LGTnI1 LGTnI3;

RUN;

Подробное описание алгоритма проверки нормальности распределения исследуемых переменных и интерпретации полученных статистических критериев представлено ранее [5]. Критерии Колмогорова-Смирнова и Шапиро-Уилка показали, что распределение переменных TnI1 и TnI3 отличается от нормального (p<0,05).

Для приведения распределения TnI1 и TnI3 к нормальному использовалось Бокс-Кокс – преобразование:

proc transreg details data=TRIALDATA ss2 plots=(transformation(dependent) obp);

WHERE CPBType=2; /* Выбираем группу коронарного шунтирования c искусственным кровообращением */

model BoxCox(TnI3) = identity(TnI1);

run;

На представленном графике λ=0,25 (рис. 1), что говорит о возможности получения нормального распределения переменных TnI1 и TnI3 с помощью возведения в степень 0,25.

 

Рис. 1.  Бокс-Кокс преобразование.

Fig. 1. Box-Cox transformation

Исходя из правила, представленного в табл. 1 преобразование переменных TnI1 и TnI3 выполнялось следующим образом:

DATA TRIALDATA;

    SET WORK.'20_06_2019 work1'n; 

    LGTnIPreOp=log(TnIPreOp); /* концентрация тропонина перед операцией */

     LGTnIEndOp=log(TnIEndOp); /* концентрация тропонина I после операции */

TRTnI1=((TnI1)**0.25-1)/0.25;/* концентрация тропонина I на 1-й день после операции (Бокс-Кокс преобразование) - возведение в  степень 0,25 */

    LGTnI2=log(TnI2); /* концентрация тропонина I на 2-й день после операции */

TRTnI3=TRTnI1=((TnI2)**0.25-1)/0.25; /* концентрация тропонина I на 3-й день после операции Бокс-Кокс преобразование  - возведение в  степень 0,25 */

    LGTnI4=log(TnI4); /* концентрация тропонина I на 4-й день после операции */

    LGTnI5=log(TnI5); /* концентрация тропонина I на 5-й день после операции */

RUN;

Таблица 1 Преобразование переменной в зависимости от λ

λ	Преобразование
≠ 0
0

 

С целью оценки наличия линейной связи между исследуемыми переменными строили скатерограмму (рис. 2).

 

Рис. 2.  Скаттерограмма переменных

Fig. 2. Scatter plot for TRTnI1 и TRTnI3

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

PROC REG DATA=TRIALDATA outest=parms plots=all;

    WHERE CPBType=2; /* Выбираем группу коронарного шунтирования без искусственного кровообращения */

MODEL TRTnI3=TRTnI1/CLB STB CLI CLM COLLIN VIF TOL DW SPEC HCC WHITE HCCMETHOD=0 INFLUENCE R LACKFIT;

    output out=reg_res r=res; /* сохраним остатки для их последующего анализа */

RUN;

Графический результат построения регрессии представлен на рис. 3

Рис. 3.  Линейная регрессия зависимости концентрации тропонина I на 3 - й день после коронарного шунтирования (TRTnI3) от концентрации тропонина I на 1 – й день после операции (TRTnI1)

Примечание. Observations – количество наблюдений, которые использовались для построения линейной регрессии. Parameters – количество коэффициентов регрессии, Error DF – степень свободы ошибки, MSE – среднеквадратическая ошибка, R-square – R² (коэффициент детерминации), Adj R-Square – скорректированный коэффициент детерминации.

Fig. 3. Linear regression for relation of troponin I 3^d day after coronary artery bypass grafting (TRTnI3) and troponin I 1^d day after coronary artery bypass grafting (TRTnI1)

Observations – number of observations. Parameters – number of coefficients, Error DF – error freedom degree, MSE – mean squared error, R-square – coefficient of determination, Adj R-Square – adjusted coefficient of determination.

На рис. 3 представлены характеристики регрессии: коэффициент детерминации R² = 0.66, что является хорошим показателем. Для построения линейной регрессии использовались 35 наблюдений, остальные наблюдения не использовались, так как в случае отсутствия хотя бы одного наблюдения в паре, SAS исключает его из расчетов.

Также оценивался квадратный корень из среднеквадратической ошибки (RMSE), чем меньше его значение, тем лучше (рис. 4).

 

Рис. 4.  Параметры линейной регрессии зависимости концентрации тропонина I на 3 - й день после коронарного шунтирования от концентрации тропонина I на 1 – й день после операции

Fig. 4. Linear regression parameters for relation of troponin I 3^d day after coronary artery bypass grafting (TRTnI3) and troponin I 1^d day after coronary artery bypass grafting (TRTnI1)

Также оценивалась гомоскедастичность остатков (постоянство дисперсии остатков). На рис. 4 в разделе Heteroscedasticity Consistent p<0,05 (нулевую гипотезу о гомоскедастичности остатков следует отвергнуть) для коэффициентов регрессии, что говорит о наличии гетероскедастичности (дисперсия остатков непостоянна). Таким образом, полученная модель линейной регрессии с рассчитанными коэффициентами может быть недостаточно точной. SAS позволяет рассчитать коэффициенты регрессии с поправкой на гетероскедастичность. Далее будет представлен код, позволяющий это сделать.

Еще одним показателем качества регрессии является критерий который приближается к оптимальному значению 2 (рис. 5).

 

Рис. 5.  Критерий Durbin-Watson для построенной линейной регрессии

Fig. 5. Durbin-Watson statistic for linear regression

Расстояние Кука (Cook’s D, показывает наличие выбросов) должно быть примерно одинаковым и менее 2, что будет говорить о хорошем качестве построенной модели (рис. 6).

Рис. 6.  Расстояние Кука для построенной линейной регрессии

Fig. 6. Cook’s distance for linear regression

В данном случае расстояние Кука неодинаково и превышает 2, что говорит о погрешности рассчитанных коэффициентов регрессии.

Распределение остатков соответствовало характеру нормального (рис. 7).

Рис. 7.  Распределение остатков

Fig. 7. Residuals distribution

Помимо визуальной оценки характера распределения остатков, были рассчитаны критерии Колмогорова-Смирнова и Шапиро-Уилка:

PROC UNIVARIATE DATA=reg_res NORMAL PLOT;

    VAR res;

RUN;

Остатки имели нормальный характер распределения, так как значимость критериев Колмогорова-Смирнова и Шапиро-Уилка p>0,05 (рис. 8).

Рис. 8.  Тест нормальности распределения остатков

Fig. 8. Normality test for residuals

 

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

proc report data=parms;

column TRTnI1 Intercept;

   define TRTnI1/ 'b' FORMAT=8.2  ;

   define Intercept/ 'a' FORMAT=8.2  ;

run;

Результаты представлены в табл. 2.

Таблица 2 Коэффициенты линейной регрессии

b	a
0.72	-0.80

Полученные коэффициенты необходимо подставить в уравнение линейной регрессии . Учитывая преобразование Бокс-Кокса, выполненное ранее, уравнение линейной регрессии примет вид

b + a

Отсюда зависимость концентрации тропонина I на 3-й день после коронарного шунтирования от концентрации тропонина I на 1-й день после операции будет представлена формулой .

Для построения графика зависимости TnI3 от TnI1 в соответствие с вышеприведенной формулой, использовался следующий код:

data function;

    set parms;

    do TnI1 = 0 to 30 by 0.1;

    TnI3 = (TRTnI1*(TnI1**0.25)+Intercept/4)**4;

    output;

    end;

run;

proc sgplot data=function;

series x=TnI1 y=TnI3;

xaxis min=0; /* не показывать ось х менее 0 */

run;

Рис. 9.  График зависимости концентрации тропонина I на 3-й день после коронарного шунтирования от концентрации тропонина I на 1-й день после операции

Fig. 9. Linear regression parameters for relation of troponin I 3^d day after coronary artery bypass grafting (TnI3) and troponin I 1^d day after coronary artery bypass grafting (TnI1)

Как уже упоминалось выше, полученные коэффициенты регрессии могут иметь существенную погрешность, т. к. имеется гетероскедастичность остатков (рис. 4), различное расстояние Кука, превышающее 2 (рис. 6, рис. 10)

Рис. 10.  Расстояние Кука для модели линейной регрессии

Fig. 10. Cook’s distance for linear regression model

 

В таком случае для расчета коэффициентов линейной регрессии следует использовать более сложные методы вычисления коэффициентов регрессии. Следующий алгоритм реализует расчет коэффициентов устойчивой (робастной) регрессии:

PROC ROBUSTREG DATA=TRIALDATA  outest=parms_robust plots=all FWLS METHOD=M (WF=HUBER) ;

    WHERE CPBType=2;

    model TRTnI3=TRTnI1;

    output out = robust_out stdp=ST_DP weight=wgt residual=robust  stdi=ST_DI;

RUN;

R²=0.61 говорит о хорошем качестве регрессионной модели (рис. 11).

 

Рис. 11. Устойчивая (робастная) регрессия. TRTnI1 – концентрация тропонина I на 1 – й день после операции (переменная после Бокс-Кокс преобразования), TRTnI3 – концентрация тропонина I на 3 – й день после операции (переменная после Бокс-Кокс преобразования).

Fig. 11. Robust regression. TRTnI1 – troponin I level 1^st day after coronary artery bypass grafting (variable after Box-Cox transformation), TRTnI3 – troponin I level 3^d day after coronary artery bypass grafting (variable after Box-Cox transformation)

 

Остатки были проверены на нормальность распределения с помощью следующего кода:

PROC UNIVARIATE DATA=robust_out NORMAL PLOT;

    VAR robust;

RUN;

Значимость критериев Колмогорова-Смирнова и Шапиро-Уилка p>0,05, что говорит о нормальном распределении остатков и адекватности расчета коэффициентов регрессии (рис. 12).

Рис. 12.  Проверка на нормальность распределения остатков устойчивой (робастной) линейной регрессии

Fig. 12. Normality test for residuals of robust regression

Для построения графиков линейной и устойчивой регрессии, а также вывода коэффициентов регрессии на экран использовался следующий код:

DATA PARMS_ROBUST; /* объединение в одном наборе данных коэффициентов линейной и устойчивой регрессии */

    SET PARMS_ROBUST;

    TRTnI1_RB=TRTnI1;

    Intercept_RB=Intercept;

RUN;

DATA PARMS;

    MERGE PARMS_ROBUST PARMS;

RUN;

data function; /* построение графика в наборе данных */

    set parms;

         exp=exp(1);

             do x = 0 to 20 by 0.1;

                 y = (TRTnI1*(x**0.25)+Intercept/4)**4;

                 y_robust=(TRTnI1_RB*(x**0.25)+Intercept_RB/4)**4;

    output;

    end;

run;

 

proc sgplot data=function;

series x=x y=y/ curvelabel='линейная регрессия';

series x=x y=y_robust/ curvelabel='устойчивая регрессия';

xaxis min=0;

label y='тропонин I на 3-й день после КШ, нг/мл';

label y_robust='тропонин I на 3-й день после КШ, нг/мл';

label x='тропонин I на 1-й день после КШ, нг/мл';

run;

TITLE 'Коэффициенты линейной и устойчивой (робастной) регрессии';

proc report data=parms;

column TRTnI1  Intercept TRTnI1_RB  Intercept_RB;

   define TRTnI1/ 'b' FORMAT=8.2  ;

   define Intercept/ 'a' FORMAT=8.2  ;

   define TRTnI1_RB/ 'b (Robust)' FORMAT=8.2  ;

   define Intercept_RB/ 'a (Robust)' FORMAT=8.2  ;

run;

Таким образом, графики зависимости TnI3 от TnI1 согласно уравнениям линейной и устойчивой регрессии, представлены на рис. 13.

 

Рис. 13. Линейная и устойчивая (робастная) регрессия модели предсказания концентрации тропонина I на 3 - й день после коронарного шунтирования (КШ)

в зависимости от концентрации тропонина I на 1 - й день после операции

Fig. 13. Linear and robust regression for troponin I level on 3^d day after operation coronary artery bypass grafting prediction depending on troponin I level on 1^st day after operation

Примечание. КШ – коронарное шунтирование.

Коэффициенты регрессии представлены в табл. 3

Таблица 3 Коэффициенты линейной и устойчивой (робастной) регрессии

b	a	b (Robust)	a (Robust)
0.72	-0.80	0.78	-0.96

 

Подставив в уравнение  коэффициенты устойчивой регрессии получим . Это более точное уравнение регрессии, так как внесена поправка на гетероскедастичность остатков и выбросы (рис. 4, рис. 5, рис. 10). На рис. 12 можно вдеть, что в случае отсутствии такой поправки, построенная модель занижает прогнозируемую концентрацию тропонина I.

Таким образом для построения модели линейной регрессии необходимо, чтобы исследуемые переменные были количественные, а наблюдения независимые. Минимальное количество наблюдений по данным литературы колеблется от 10 до 20 [1]. Также необходимо предпринять следующие шаги:

1. Определить характер распределения исследуемых переменных (количественные). Ранее рассмотрено, как это сделать [5].
2. Если распределение отличается от нормального, то провести преобразование к нормальному виду (Бокс-Кокс преобразование).
3. Если распределение исследуемых переменных не удалось привести к нормальному, следует отказаться от построения линейной регрессии и использовать другие методы статистической обработки данных.
4. Если распределение исследуемых переменных имеет характер нормального или удалось преобразовать к нормальному виду, возможно построение линейной регрессии.
5. Построить скатерограмму для оценки наличия возможной линейной связи [2, 7].
6. После построения линейной регрессии и вычисления коэффициентов регрессии необходимо оценить возможность ошибки (табл. 4) [7].
7. Если какие-либо требования к линейной регрессии из табл. 4 не соблюдаются, то следует вносить поправку при расчете коэффициентов регрессии, в рассмотренном примере использовалась устойчивая (робастная) регрессия, однако, существуют и другие методы расчета коэффициентов.
8. После выполнения всех необходимых тестов и преобразований необходимо подставить соответствующие коэффициенты в уравнение

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Таблица 4 Ключевые параметры и тесты для оценки линейной регрессии

Параметр	Чем проверить	Характеристика
		Хорошо	Плохо
Распределение исследуемых переменных нормальное	Тест на нормальность (тест Колмогорова-Смирнова или др.)	p>0,05 - нормальное	p<0,05 - отличное от нормального
Рекомендации	Преобразование к нормальному виду (Бокс-Кокс преобразование)
Автокорреляция (независимость наблюдений)	Критерий Дарбина-Уотсона (DW)	DW приближается к 2,0	DW приближается к 0 или 4
Рекомендации	Использование процедуры Autoreg
Гомоскедастичность (постоянная дисперсия остатков)	Тест Уайта, Тест Бройша — Пагана	p>0,05 гомоскедастичность	p<0,05 гетероскедастичность
Рекомендации	Использование устойчивой (робастной) регрессии
Распределение остатков нормальное	Тест на нормальность (тест Колмогорова-Смирнова или др.)	p>0,05 нормальное	p<0,05 отличное от нормального
Расстояние Кука	Cook’s D distance	Примерно одинаковое для всех наблюдений и не превышает 2	Разное, >2
Рекомендации	Использование устойчивой (робастной) регрессии
Квадратный корень из среднеквадратической ошибки	RMSE	Маленькое значение	Большое значение
Коэффициент детерминации/ Скорректированный коэффициент детерминации	R2/R2 adj	Приближается к 1	Менее 0,6
Рекомендации	Попробовать применить методы непараметрической статистики

Параметры линейной и устойчивой (робастной) регрессии могут быть представлены в виде таблиц в стиле APA (табл. 5, табл. 6).

Таблица 5 Параметры линейной регрессии

 

	B	95% CI	β	t	p
a	-0,79	[-1,22; -0,37]	0	-3,83	0,0005
b	0,72	[0, 54; 0, 91]	0,81	7,98	0,0001

Таблица 6 Параметры устойчивой (робастной) регрессии

 

	B	95% CI	χ2	p
a	-0,96	[-1,3; -0,62]	30,53	0,0001
b	0,78	[0, 64; 0, 93]	107,6	0,0001

 

Выводы

Алгоритм построения линейной регрессии, разработанный в SAS Enterprise Guide 6.1, может использоваться для обработки данных научных и клинических исследований.

Финансирование

Работа выполнена в рамках государственного задания по теме «Оценка регенеративного потенциала пациента при операциях на сердце».

Соблюдение этических норм

Выполнение исследования одобрено протоколом локального этического комитета ФГБОУ ВО ПСПбГМУ им. И.П. Павлова.

 

Конфликт интересов

Авторы заявляют об отсутствии конфликта интересов.

 

 

 

 

 

 

Об авторах

Николай Сергеевич Буненков

ФГБОУ ВО «Первый Санкт-Петербургский государственный медицинский университет им. акад. И.П. Павлова»

Автор, ответственный за переписку.
Email: bunenkov2006@gmail.com
ORCID iD: 0000-0003-4331-028X

аспирант кафедры факультетской хирургии

Россия, Санкт-Петербург

Владимир Владимирович Комок

ФГБОУ ВО «Первый Санкт-Петербургский государственный медицинский университет им. акад. И.П. Павлова»

Email: vladimir_komok@mail.ru
ORCID iD: 0000-0002-3834-7566
SPIN-код: 3572-5180

канд. мед. наук, кардиохирург отделения кардиохирургии № 2

Россия, Санкт-Петербург

Василий Андреевич Бобыльков

ФГБОУ ВО «Первый Санкт-Петербургский медицинский университет им. акад. И.П. Павлова» Минздрава России, Санкт-Петербург

Email: basilbob@yandex.ru

клинический ординатор

Россия

Никита Владимирович Грудинин

ФГБУ НМИЦ Трансплантологии и искусственных органов имени академика В.И. Шумакова, Минздрава России, Москва

Email: nikita.grudinin@mail.ru
ORCID iD: 0000-0003-1100-9513

кардиохирург отделения трансплантации легких

Москва

Гульнара Фазулиевна Буненкова

ФГБОУ ВО «Первый Санкт-Петербургский государственный медицинский университет им. акад. И.П. Павлова»

Email: gulnara533@gmail.com

клинический ординатор кафедры госпитальной терапии, 

Россия, Санкт-Петербург

Олег Александрович Гриненко

ФГБОУ ВО «Первый Санкт-Петербургский государственный медицинский университет им. акад. И.П. Павлова»

Email: klinika@spb-gmu.ru

д-р мед. наук, проректор по лечебной работе

Россия, Санкт-Петербург

Александр Сергеевич Немков

ФГБОУ ВО «Первый Санкт-Петербургский государственный медицинский университет им. акад. И.П. Павлова»

Email: nemk_as@mail.ru
ORCID iD: 0000-0002-5152-0001
SPIN-код: 2853-4634

д-р мед. наук, профессор, кардиохирург отделения кардиохирургии № 2

Россия, Санкт-Петербург

Список литературы

Дополнительные файлы