Дисперсійний аналіз - (ANalysis Of VAriances, ANOVA) використовується для дослідження зв’язку між якісними (номінальними) ознаками та кількісною (неперервною) змінною. По суті, дисперсійний аналіз здійснює тестування гіпотези про рівність середніх арифметичних декількох вибірок. Таким чином його можна розглядати як параметричний критерій для порівняння центрів відразу кількох вибірок. Якщо використовувати дисперсійний аналіз для двох вибірок, то результати будуть ідентичні до результатів t-критерію Стьюдента. Однак, на відміну від інших статистичних критеріїв, дисперсійний аналіз дозволяє дослідити проблему більш детально. Дисперсійний аналіз базується на законі, що сума квадратів відхилень об’єднаної вибірки дорівнює сумі квадратів внутрішньогрупових відхилень та сумі квадратів міжгрупових відхилень. Дисперсійний аналіз використовує критерій Фішера для встановлення статистичної значимості відмінності міжгрупових дисперсій від внутрішньогрупових. Однак для цього необхідними передумовами є нормальність розподілу та гомоскедастичність (рівність дисперсій) вибірок. Розрізняють одновимірний (однофакторний) дисперсійний аналіз та багатовимірний (багатофакторний). Одновимірний досліджує залежність досліджуваної величини від однієї ознаки, багатовимірний – відразу від багатьох, а також дозволяє виявити зв’язок між ознаками. Непараметричним аналогом одновимірного дисперсійного аналізу є критерій Краскела-Уолліса.

##Приклад 1 Необхідно виявити чи залежить тривалість надання фармацевтичної опіки від категорії та досвіду фармацевтичного фахівця. Дані для цього прикладу вигадані. Завантажимо їх з інтернету:

adress <- "http://stat.org.ua/data/anova.csv"
data <- read.csv(adress, sep = "\t")
head(data)  #виводимо декілька перших рядочків на екран

##   N time       edu    exp
## 1 1    3 farmacevt   malo
## 2 2    4  provisor bagato
## 3 3    5 farmacevt bagato
## 4 4    7  provisor   malo
## 5 5    6  provisor bagato
## 6 6    7 farmacevt   malo

Залежною змінною є час (time). Незалежними факторами (дискретними, номінальними величинами) є освіта (edu) та досвід (exp). Таким чином досліджується залежність відносної величини від двох номінальних. Отже слід застосувати двофакторний дисперсійний аналіз. Перед цим перевіряємо нормальний розподіл. Так як у нас відразу два фактори утворюють групи, простіше перевіряти відповідність залишків нормальному розподілу, отримуючи ці залишки з готової моделі дисперсійного аналізу (отримується шляхом виклику aov()) за допомогою функції residuals():

result <- aov(time ~ edu + exp, data = data)
shapiro.test(residuals(result))

## 
## 	Shapiro-Wilk normality test
## 
## data:  residuals(result)
## W = 0.9318, p-value = 0.1675

Оскільки p-величина > 0.05, то припускаємо нормальність розподілу та переходимо до тестування відмінностей дисперсій вибірок. Для цього використаємо параметричний критерій множинного порівняння дисперсій, напр. критерій Бартлетта. Оскільки ділити потрібно відразу по двом ознакам, об’єднаємо їх у одну чотирирівневу ознаку за допомогою функції interaction().

bartlett.test(data$time ~ interaction(data$exp, data$edu))

## 
## 	Bartlett test of homogeneity of variances
## 
## data:  data$time by interaction(data$exp, data$edu)
## Bartlett's K-squared = 0.8593, df = 3, p-value = 0.8352

Cуттєвих відмінностей не виявлено, тому можемо проводити дисперсійний аналіз. Власне, саму модель ми вже побудували, залишилось проглянути результати:

summary(result)

##             Df Sum Sq Mean Sq F value Pr(>F)   
## edu          1   46.9    46.9   12.18 0.0028 **
## exp          1   23.7    23.7    6.16 0.0238 * 
## Residuals   17   65.4     3.8                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Рядками цієї таблиці є незалежні фактори, а стовпцями – деякі статистичні показники. Нас цікавить останній – p-величина. Як бачимо, на тривалість опіки впливає і кваліфікація, і стаж. Щоб виявити напрям і розмір цього впливу, варто провести додатковий післяаналіз (post-hoc). Для цього можна застосувати критерій Тьюкі.

TukeyHSD(result)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = time ~ edu + exp, data = data)
## 
## $edu
##                     diff   lwr   upr  p adj
## provisor-farmacevt 3.125 1.236 5.014 0.0028
## 
## $exp
##              diff    lwr   upr  p adj
## malo-bagato 2.096 0.2356 3.956 0.0295

Критерій вказує нам різницю між вибірками, сформованими за допомогою поділу за рівнями факторів, її довірчий інтервал (із ймовірністю 95% різниця між вибірками знаходиться у цьому інтервалі), та скоректовані значення p-величини гіпотези про рівність центрів вибірок. Отож, в середньому провізор витрачає на 3хв більше часу на фармацевтичну опіку, ніж фармацевт, а також мало досвідчені аптекарі витрачають на 2хв більше, ніж досвідчені.

##Приклад 2 Проаналізуємо результати вивчення кількісного вмісту амінокислот досліджуваних зразків трави квасолі [1]. Дані наведено у відсотках від сухої ваги. Сформулюємо дві нульові гіпотези:

1. Вміст амінокислот у досліджуваних видах квасолі не залежить від природи амінокислот (однаковий для всіх амінокислот)
2. Вміст амінокислот не залежить від виду квасолі (однаковий для всіх досліджуваних видів).

Тобто існує одна залежна змінна (неперервна величина) – вміст речовини, і два незалежні фактори (дискретні, номінальні величини) – природа (назва) амінокислоти та вид квасолі. Таким чином дослідження слід проводити за допомогою двофакторного дисперсійного аналізу. Завантажимо дані з інтернету. Оскільки таблиця складена українською мовою, на всякий випадок зазначимо кодування файлу - CP1251 (також відоме під назвою Windows-1251)

adress <- "http://pharm.uz.ua/files/kvasolya.csv"
data <- read.csv(adress, sep = ";", fileEncoding = "CP1251")
fix(data)  #перевіримо успішність завантаження

Якщо із розшифруванням файлу виникнуть труднощі - спробуйте видалити аргумент разом з комою перед ним і повторити. Звернемо увагу на структуру датафрейму. Значення факторів формують стовпці (види квасолі) та рядки (назви амінокислот) отриманої таблиці, а у кожній клітинці – вміст відповідної амінокислоти у відповідному виду квасолі. Дані у такому вигляді називають крос-табульованими. Щоб проводити дисперсійний аналіз, потрібно перевести дані до звичайної, “довгої” таблиці, де кожен рядок відповідає одиничному спостереженню із характеристиками “вміст”, “амінокислота”, “вид”. Для цього існує функція reshape()

data <- reshape(data, direction = "long", varying = list(2:6), v.names = "vmist", 
    ids = data$Речовина, times = names(data)[-1], timevar = "species")

Перетворюмо data у “довгу” (direction="long") таблицю, по черзі записуючи значення 2-6 стовпців (varying=list(2:6)), залежну змінну назвемо “vmist” (v.names="vmist"), об’єктами спостереження приймемо амінокислоти (ids=data$Речовина), групуюча змінна – назви видів (times=names(data)[-1]), і назвемо цю групуючу змінну “species" (timevar="species")

fix(data)  #подивимось на результат перетворення:

Назвами рядочків стали унікальні комбінації факторів амінокислота.вид", другий стовпчик містить назви амінокислот. Ми також продублювали назви амінокислот в останньому стовпці (id). Третій стовпчик – фактор species - види квасолі, і четвертий стовчик – залежна змінна vmist. Перевіримо нормальність розподілу вибірок для цієї задачі. Для цього відразу побудуєм модель дисперсійного аналізу:

res <- aov(vmist ~ species + id, data = data)
shapiro.test(residuals(res))

## 
## 	Shapiro-Wilk normality test
## 
## data:  residuals(res)
## W = 0.9802, p-value = 0.2518

Результат: \(p>0.05\), немає статистично значимих відмінностей розподілу від нормального. Перевіримо рівність дисперсій.

bartlett.test(data$vmist ~ data$species)

## 
## 	Bartlett test of homogeneity of variances
## 
## data:  data$vmist by data$species
## Bartlett's K-squared = 0.3824, df = 4, p-value = 0.9839

Дисперсії рівні, досліджуєм результати дисперсійного аналізу.

summary(res)

##             Df Sum Sq Mean Sq F value Pr(>F)    
## species      4    0.5    0.13    1.68   0.17    
## id          15   96.8    6.45   86.18 <2e-16 ***
## Residuals   60    4.5    0.07                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

Твердження: якщо потрібно перевірити тільки гіпотезу про рівність вмісту амінокислот у видах квасолі, то слід виконати однофакторний дисперсійний аналіз виду

res <- aov(vmist~species,data=data)

Спростування: в такому разі усі значення вмісту амінокислот для одного виду розглядаються як результати визначення однієї і тієї ж речовини, тобто є одне істинне значення, а відхилення від нього – випадкові (дисперсія). І ми протестуємо рівність центрів незалежних вибірок. А насправді – вибірки залежні, кожне спостереження вмісту одної амінокислоти пов’язане із вмістом цієї ж амінокислоти в інших видах. Таким чином ми отримаємо суттєво інше (більше) значення p-величини. Подивимося:

summary(res)

##            Df Sum Sq Mean Sq F value Pr(>F)
##species      4    0.5   0.126   0.093  0.984
##Residuals   75  101.2   1.350    

Вихід: слід вказати при проведенні аналізу, що vmist може залежати і від id (амінокислоти), однак ця залежність – похибка зумовлена експериментом (тобто нас не цікавить). Такий метод аналізу називається “дисперсійний аналіз з повторюваними вимірюваннями”.

res <- aov(vmist~species+Error(id), data=data)
summary(res)

##Error: id
##          Df Sum Sq Mean Sq F value Pr(>F)
##Residuals 15  96.75    6.45               
##
##Error: Within
##          Df Sum Sq Mean Sq F value Pr(>F)
##species    4  0.504 0.12599   1.683  0.166
##Residuals 60  4.491 0.07484  

Результат для рядка species той же, що при res <- aov(vmist ~ species + id, data = data).