Кореляційний та регресійний аналіз – це статистичні дослідження залежності між випадковими величинами. Різниця між цими двома методами у меті: для кореляційного аналізу важливо виявити, чи існує зв’язок між двома величинами, а для регресійного – аналітичний вираз (функцію регресії) цієї залежності. Таким чином, регресійний аналіз створює математичну модель зв’язку, а кореляційний аналіз формально відношення до математичного моделювання не має. Однак, наприклад, коефіцієнт кореляції Пірсона спроможний оцінити силу лише лінійної залежності між величинами. Загалом, для дослідження залежностей між двома величинами використовують наступні методи:

Лінійна залежність між двома величинами має наступний вид: \[y=a_0+bx\], де \(y\) – залежна змінна, \(a_0\) – зсув – значення функції при значенні аргумента рівному нулю, \(b\) – нахил, або кутовий коефіцієнт, або коефіцієнт пропорційності, \(х\) – аргумент, або незалежна змінна (у випадку регресійного аналізу також прийнято \(y\) називати “відгук”, а \(x\) – “предиктор” або “регресор”). При \(b>0\) значення функції \(y\) зростає тоді, коли зростає значення аргументу – пряма пропорційність. При \(b<0\) значення функції зростає тоді, коли значення аргументу зменшується – обернена пропорційність.

Модель множинної лінійної регресії має вигляд \(y=a_0+a_1x_1+a_2x_2+…+a_nx_n\), де \(а_1…а_n\) – коефіцієнти пропорційності, \(x_1…x_n\) – предиктори (регресори).

Для побудови моделі множинної лінійної регресії у R використовується функція lm(), аргументом якої є формула залежності. Крім того, за допомогою нелінійних перетворень як відгуку, так і регресорів, можна нелінійні моделі виразити через лінійні і таким чином знайти їх коефіцієнти та обчислити коефіцієнт кореляції. Якщо кількість незалежних між собою регресорів дорівнює кількості спостережень у залежній змінній, то навіть при відсутності будь-якого зв’язку між відгуком і регресорами, можна розв’язати систему лінійних алгебраїчних рівнянь і отримати точну залежність y від множини x. Причому, якщо кількість незалежних між собою регресорів є більшою від кількості спостережень у залежній змінній, то таких точних розв’язків є безліч. Тому із множини потенційних регресорів часто потрібно вибрати невелику підмножину кращих, які і описують поведінку відгуку. Існує багато методів для здійснення такого вибору, одним з яких є метод покрокової регресії. Даний метод на кожному кроці приймає поточну модель (для першого кроку як модель приймають середнє значення відгуку), і аналізує ефективність нових моделей, утворених шляхом додавання до поточної моделі одного із регресорів. Таким чином знаходиться один кращий на даному етапі регресор, він вноситься у модель. Також перед кожним кроком поточна модель перевіряється щодо можливості виключення із неї зайвого регресора. Варто пам’ятати, що при внесенні у модель будь-якої випадкової незалежної змінної кореляційний коефіцієнт даної регресії трошки зросте. Тому при регресійному аналізі слід дотримуватись принципу “леза Оккама” - “не приумножай сутностей без конче на те необхідності”. Також варто пам’ятати, що наявність значної кореляції, хоч часто і означає наявність деякого зв’язку, однак далеко не завжди цей зв’язок є причинно-наслідковим. Так, існує відомий приклад про пожежників (Елисеева И.И., Юзбашев М.М. Общая теория статистики: Учебник / Под ред. И.И. Елисеевой — 4-е издание, переработанное и дополненное. — Москва: Финансы и Статистика, 2002. — 480 с): Дослідивши розміри збитків при різних пожежах та кількості пожежників, що були задіяні для ліквідації пожеж, можна виявити сильний кореляційну залежність: чим більше пожежників гасять пожежу, тим більшими є збитки від неї. Але це не значить, що звільнивши пожежників, збитки від пожеж будуть меншими. Справа в тому, що існує ще одна величина – сила пожежі. Чим сильніша пожежа, тим більше пожежників залучається до її ліквідації, і одночасно більшими є і збитки. Також сильний кореляційний зв’язок може виявитись випадковістю, наприклад при аналіз великих об’ємів даних. Так, було знайдено, що рівень смертності на дорогах у США добре корелює із кількістю лимонів, імпортованих до США із Мексики.

##Приклад 1. Розглянемо ще раз таблицю оптових та роздрібних цін на лікарські препарати компанії “Фармак” (див. приклад для розділу “АВС та XYZ–аналіз”). Визначимо, наскільки сильний кореляційний зв’язок між оптовою ціною на лікарський засіб та його середньою націнкою. Так, якщо націнка на лікарські препарати формується як певний відсоток від оптової ціни, то очікується сильний прямий кореляційний зв’язок.

adress <- "http://stat.org.ua/data/abc.csv"
dani <- read.csv(adress, sep = ";", fileEncoding = "CP1251")
# Якщо у вас Windows української чи російської локалізації - аргумент
# fileEncoding можна пропустити
nacinka <- dani$aptcina - dani$optcina  #обчислили націнку.

Коефіцієнти кореляції обчислюються за допомогою функції cor.test(). Ця функція по замовчуванню дає нам кореляційний коефіцієнт Пірсона.

cor.test(dani$optcina, nacinka)

## 
## 	Pearson's product-moment correlation
## 
## data:  dani$optcina and nacinka
## t = 4.695, df = 50, p-value = 2.111e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3301 0.7177
## sample estimates:
##    cor 
## 0.5532

Тобто кореляційний коефіцієнт Пірсона знаходиться внизу, \(cor≈0.55\), що означає наявність слабкого зв’язку. Крім того, p-величина для гіпотези про рівність коефіцієнта кореляції нулю достатньо мала, тобто можна стверджувати про невипадковий характер зв’язку. Однак варто пам’ятати, що коректне визначення цієї p-величини та обчислення довірчих інтервалів для коефіцієнта кореляції Пірсона можливе лише при нормальному розподілі даних. Перевіримо оптову ціну щодо нормальності розподілу:

shapiro.test(dani$optcina)

## 
## 	Shapiro-Wilk normality test
## 
## data:  dani$optcina
## W = 0.5914, p-value = 8.627e-11

Сукупність оптових цін не підкоряється закону нормального розподілу. Cлід обережно ставитись до отриманої вище р-величини. Тим не менше, саме значення коефіцієнта кореляції Пірсона від нормальності не залежить.

Обчислимо також кореляційний коефіцієнт Спірмена (ранговий):

cor.test(dani$optcina, nacinka, method = "spearman")

## Warning: Cannot compute exact p-value with ties

## 
## 	Spearman's rank correlation rho
## 
## data:  dani$optcina and nacinka
## S = 3089, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.8682

Тобто кореляційний коефіцієнт Спірмена \(\rho≈0.87\). Бачимо, що коефіцієнт кореляції Спірмена значно більший, ніж Пірсона. Це є індикатором нелінійної залежності націнки від оптової ціни.

Побудуємо лінійну регресійну модель. Оскільки націнку зазвичай утворюють як певний відсоток від оптової ціни, то у лінійній кореляційній функції зсув \(a_0\) буде близький до нуля (гіпотетична границя “0 грн ціна – 0 грн націнка”) , а коефіцієнт \(b\) – близький за значенням до цього відсотка (у формі частки від цілого).

model <- lm(nacinka ~ dani$optcina)
summary(model)

## 
## Call:
## lm(formula = nacinka ~ dani$optcina)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -46.22  -5.91  -4.07   0.64 114.61 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    6.7757     3.2240     2.1    0.041 *  
## dani$optcina   0.0591     0.0126     4.7  2.1e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.9 on 50 degrees of freedom
## Multiple R-squared:  0.306,	Adjusted R-squared:  0.292 
## F-statistic:   22 on 1 and 50 DF,  p-value: 2.11e-05

Residuals – відхилення, або залишки. Досить серйозна відмінність медіани залишків від нуля та значна диспропорція між абсолютними значеннями максимуму та мінімуму свідчать про те, що препарат, якого стосується цей максимум, для даної моделі є викидом.

Coefficients – таблиця коефіцієнтів. Перший рядок стосується зсуву \(а_0\), другий – коефіцієнт пропорційності регресора optcina \(b\). Перший стовпчик – встановлені значення коефіцієнтів. Із значень коефіцієнтів випливає, що найбільш ближчий метод формування націнки полягає у додаванні фіксованої додаткової вартості (6.78 грн) та 5.9% від оптової ціни препарату. Другий стовпчик – стандартні відхилення коефіцієнтів – дозволяють оцінити наскільки сильно може відхилятись істинне значення коефіцієнту від обчисленого (досить сильно може відхиляєтися зсув). Третій стовпчик – значення t-статистики для нульової гіпотези про рівність коефіцієнта нулю. Четвертий стовпчик – p-величина для цієї гіпотези. Оскільки розподіл отриманої t-статистики описується розподілом Стьюдента тільки при нормально розподілених відгуці та регресорі (а попередньо нами було встановлено невідповідність регресора умовам нормальності), то отримана p-величина є не зовсім коректною. Тобто стверджувати статистично значиму відмінність \(а_0\) від нуля не можна. Однак, оскільки p-величина для коефіцієнта регресора є достатньо малою, то можна стверджувати, що наявність кореляційного зв’язку між націнкою та оптовою ціною є статистично значимою. Нижче також подано значення стандартної похибки відхилень (Residual standard error) – 19.9 - досить велике, множинний коефіцієнт кореляції (коефіцієнт детермінації, Multiple R-squared) – 0.306 та скоректований коефіцієнт детермінації (Adjusted R-squared) – 0.292. Скоректований коефіцієнт детермінації враховує поправку на кількість регресорів і є придатним для порівняння регресійних моделей із різною кількістю регресорів. Також подано значення F-критерію (критерію Фішера для порівняння дисперсій) та його p-величина (при цьому порівнюється дисперсія відгуку із дисперсією залишків, якщо гіпотеза про рівність дисперсій не відхиляється, то отримана модель не відрізняється статистично значимо від “нульової” моделі - середнього арифметичного). Але так як ні розподіл відгуку, ні залишків не відповідає нормальному (це можна перевірити за допомогою відповідного критерію), а критерій Фішера належить до параметричних методів аналізу, то отримані значення також є не зовсім коректними. Побудуємо графік залежності націнки від оптової ціни та зобразимо на ньому отриману лінійну модель:

plot(dani$optcina, nacinka)  #графік залежності націнки від оптової ціни
abline(model$coefficients)  #лінійна модель залежності

На графіку помітно викид (препарат із націнкою біля 140грн). Крім того, при зростанні оптової ціни націнка спочатку збільшується швидко, а потім повільно. Тобто більш підходящою є нелінійна залежність. Описана поведінка може бути апроксимована (наближена), наприклад, за допомогою функції квадратного кореня. Визначимо регресійну функцію наступною формулою: \(y=a_0+b\sqrt{x}\). Видно, що після добування кореня формула ідентична лінійній функції. Отож:

model2 <- lm(nacinka ~ sqrt(dani$optcina))
summary(model2)

## 
## Call:
## lm(formula = nacinka ~ sqrt(dani$optcina))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -46.25  -2.75  -1.65  -0.31 106.97 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -4.769      4.102   -1.16     0.25    
## sqrt(dani$optcina)    2.141      0.358    5.99  2.3e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.3 on 50 degrees of freedom
## Multiple R-squared:  0.418,	Adjusted R-squared:  0.406 
## F-statistic: 35.9 on 1 and 50 DF,  p-value: 2.27e-07

Бачимо, що коефіцієнт детермінації зріс. Однак значення його все одно є невисоким. Вилучення викиду значно покращить ситуацію. Тому спробуйте це зробити самостійно і повторіть побудову моделей. Ця задача є навчальною, оптові ціни і роздрібні ціни стосуються різних часових періодів, тому отримані результати не відображають дійсності.

##Приклад 2. При пошуку нових лікарських засобів інтенсивно використовуються сучасні комп’ютерні технології та обчислення. У даному прикладі ми проведемо дослідження зв’язку між структурою сполук і їх біологічною активністю. Такі дослідження називаються QSAR (Quantitative structure-activity relationship). У результаті ми отримаємо математичну модель, за допомогою якої можна прогнозувати активність нових, ще не синтезованих сполук, і, таким чином, вибирати для синтезу лише ті, що мають найкращий потенціал. Візьмемо як приклад реальне дослідження Девіняка О.Т. та співавторів. У цій публікації об’єктом аналізу є логарифми констант інгібування монокарбоксилат транспортеру 1 двадцяти двома похідними тієно[2,3-d]піримідин-1,4-діону, узяті зі статті S. Guile та ін.. Монокарбоксилат транспортер 1 є новою мішенню для розробки імуномодулюючих лікарських засобів.

Виконаємо побудову лінійної моделі імуномодулюючої активності похідних тієно[2,3-d]піримідин-1,4-діону за допомогою методу покрокової регресії. Відгуком є від’ємний логарифм констант інгібування монокарбоксилат транспортеру 1 (чим більший, тим краща активність), регресорами – числові характеристики структури молекули (дескриптори). Завантажимо дані:

adress1 <- "http://stat.org.ua/data/stepwise-immuno.csv"
data <- read.csv(adress1, sep = ";")
fix(data)  #перевірка

Перший стовпчик activity містить значення біологічної активності, інші стовпчики - молекулярні дескриптори - числа, що описують різні структури тих хімічних речовин. Для проведення покрокової регресії слід буде задати максимально насичену формулу (у якій відгук залежить від усіх регресорів). Створимо таку формулу за допомогою функції paste(), яка об’єднує фрагменти тексту. Назви стовпчиків можна отримати через функцію colnames(). Із тих стовпчиків перший є відгуком, а всі інші є регресорами. У формулу слід передати тільки назви регресорів.

form = as.formula(paste("~", paste(colnames(data)[-1], collapse = "+")))

Додатковий аргумент collapse задає розділювач для створених різних значень. Нам потрібно, щоб між назвами регресорів знаходився плюс. Створюємо лінійну модель за допомогою покрокової регресії із максимальним числом кроків рівним 3 (слідуєм принципу “леза Оккама”). Оскільки при цьому здійснюється значна кількість обчислень, доведеться трохи почекати результату.

model <- step(lm(activity ~ 1, data = data), form, steps = 3)

## Start:  AIC=-15.57
## activity ~ 1
## 
##                  Df Sum of Sq  RSS   AIC
## + PEOE_RPC..1     1      6.15 3.74 -35.0
## + VDistEq         1      5.49 4.40 -31.4
## + Q_RPC..1        1      5.14 4.76 -29.7
## ....
## 
## Step:  AIC=-34.96
## activity ~ PEOE_RPC..1
## 
##                  Df Sum of Sq  RSS   AIC
## + AM1_dipole      1      0.95 2.79 -39.4
## + SlogP_VSA5      1      0.89 2.86 -38.9
## + VDistEq         1      0.86 2.88 -38.7
## ....
## 
## Step:  AIC=-39.39
## activity ~ PEOE_RPC..1 + AM1_dipole
## 
##                  Df Sum of Sq  RSS   AIC
## + weinerPath      1      1.00 1.79 -47.2
## + E_ang           1      1.00 1.80 -47.1
## + VDistEq         1      0.94 1.85 -46.5
## ....
## Step:  AIC=-47.17
## activity ~ PEOE_RPC..1 + AM1_dipole + weinerPath

Таким чином, у модель було внесено 3 кращих регресори. Подивимося на статистичні параметри отриманої моделі:

summary(model)

## 
## Call:
## lm(formula = activity ~ PEOE_RPC..1 + AM1_dipole + weinerPath, 
##     data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4572 -0.2298  0.0208  0.2405  0.5121 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.474706   1.043785    1.41   0.1748    
## PEOE_RPC..1 41.313243   4.677977    8.83  5.8e-08 ***
## AM1_dipole  -0.119618   0.031572   -3.79   0.0013 ** 
## weinerPath   0.000685   0.000216    3.17   0.0053 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.316 on 18 degrees of freedom
## Multiple R-squared:  0.819,	Adjusted R-squared:  0.789 
## F-statistic: 27.1 on 3 and 18 DF,  p-value: 6.74e-07

Коефіцієнт детермінації R2≈0.82, що вказує на досить високу точність отриманої моделі. Значення p-величин для гіпотези про рівність коефіцієнтів нулю не повинні вводити в оману: так як кожен з регресорів вибирався із великої кількості, імовірність появи випадково такого регресора значно більша і залежить від кількості доступних регресорів. Встановити правильний набір регресорів можливо завдяки більш складним методам вибору змінних та проведенню крос-валідації. Для візуалізації точності моделі побудуємо графік залежності прогнозованих активностей від істинних

plot(model$fitted.values ~ data$activity)
abline(0, 1)

Дані лежать близько до прямої, що підтверджує хорошу точність. Викидів, тобто речовин, активність яких модель не спромоглась пояснити, - немає. Слід зауважити, що при невеликих кількостях спостережень (сполук) порівняно із великою кількістю характеристик (дескрипторів), можна легко знайти моделі високої точності завдяки простій випадковості. Щоб оцінити чи дана модель має чисто випадковий характер можна завдяки тесту перестановки відгуку (Response Permutation Test). При цьому справжній відгук заміняється випадковою перестановкою його значень і далі тотожним алгоритмом будується модель. Якщо статистичні показники побудованих таким чином моделей досягають показників нашої моделі, отриманої із справжнім відгуком, то наша модель по суті є випадковою, і тому не спроможна відображати реальну залежність (якщо така існує). При реальних дослідженнях проводять від сотні до тисяч та десятків тисяч таких перестановок, однак з метою економії часу ми проведемо всього двадцять.

Спочатку ініціалізуємо за допомогою функції set.seed() генератор випадкових чисел за допомогою явно вказаного числа (у даному випадку 2901). Це проводиться для того, щоб при повторному виконанні цього прикладу читачем результат був той самий.

set.seed(2901)
fikt <- NULL  #пуста змінна для наступних 'фіктивних' значень
for (i in 1:20) {
    model2 <- step(lm(sample(activity) ~ 1, data = data), form, steps = 3)
    fikt[i] <- summary(model2)$r.squared
}

summary(fikt)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0054  0.0142  0.0525  0.0544  0.4640

Бачимо, що медіана коефіцієнтів детермінації випадкових моделей становить 0.014, а максимальне значення \(R^2\), яке отримали при симуляціях - 0.464. Тобто випадково отримати такий великий коефіцієнт детермінації \(R^2≈0.82\) практично неможливо.

Виконаємо прогноз активності сполук із ряду нових похідних тієно[2,3-d]піримідин-1,4-діону (за допомогою функції predict):

adress2 <- "http://stat.org.ua/data/predict-immuno.csv"
newdata <- read.csv(adress2, sep = ";")
fix(newdata)  #перевіряєм успішність завантаження
prognoz <- predict(model, newdata)  #прогноз
summary(prognoz)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    7.88    8.16    9.76    9.47   10.00   11.30

Таким чином, для більшої частини досліджуваних сполук прогнозується наявність імуномодулюючої активності \(pK_i>9\). Щоб подивитись номери сполук із \(pK_i>10\) виконаємо команду:

which(prognoz > 10)

## 16 17 18 19 20 21 
## 16 17 18 19 20 21

Якраз ці сполуки слід вибрати для проведення спрямованого синтезу, адже в них є високі шанси виявитись високоактивними імуномодуляторами (блокаторами монокарбоксилат транспортеру 1).