SAS ENTERPRISE GUIDE 6.1: LINEAR REGRESSION.

Full Text

Введение

При выполнении научных исследований может возникать задача предсказания значения какого-либо показателя в зависимости от концентрации того или иного лабораторного маркера [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, может использоваться для обработки данных научных и клинических исследований.

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

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

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

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

 

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

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

 

 

 

 

 

 

About the authors

Nikolay S. Bunenkov

Pavlov First Saint Petersburg State Medical University

Author for correspondence.
Email: bunenkov2006@gmail.com
ORCID iD: 0000-0003-4331-028X

aspirant, Department of Faculty Surgery

Russian Federation, Saint Petersburg

Vladimir V. Komok

Pavlov First Saint Petersburg State Medical University

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

PhD, cardiac surgeon, Department of Cardiac Surgery No. 2

Russian Federation, Saint Petersburg

Vailiy A. Bobylkov

1Pavlov First Saint Petersburg State Medical University, Saint Petersburg, Russia

Email: basilbob@yandex.ru

resident

Russian Federation

Nikita V. Grudinin

Federal State Budgetary Institution “Academician V.I. Shumakov Federal Research Center of Transplantology and Artificial Organs”, Ministry of Health of the Russian Federation, Moscow

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

caridac surgeon, Department of Lung Transplantation

Москва

Gulnara F. Bunenkova

Pavlov First Saint Petersburg State Medical University

Email: gulnara533@gmail.com

resident, Department of Hospital Therapy

Russian Federation, Saint Petersburg

Oleg A. Grinenko

Pavlov First Saint Petersburg State Medical University

Email: klinika@spb-gmu.ru

Doctor of Medical Science, Vice-Rector

Russian Federation, Saint Petersburg

Alexander S. Nemkov

Pavlov First Saint Petersburg State Medical University

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

Doctor of Medical Science, Professor, cardiac surgeon, Department of Cardiac Surgery No. 2

Russian Federation, Saint Petersburg

References

Supplementary files