Установка пакетов

Установим необходимые для занятия пакеты (это действие необходимо, только если пакеты ещё не установлены на ваш компьютер) (для установки выполните следующую строку без символа # вначале).

# install.packages(c('ggplot2', 'reshape2', 'dplyr', 'gplots', 'lsr', 'pwr'))

Подгружаем загруженные пакеты для их использования в текущей сессии работы R.

library(ggplot2)
library(reshape2)
library(dplyr)
library(gplots)
library(lsr)
library(pwr)

Двухфакторный дисперсионный анализ (Two-way ANOVA)

С помощью двухфакторного дисперсионного анализа можно проверить влияние двух факторов на зависимую переменную. Продолжим работать с данными про свинок и их зубы. Попробуем проверить, влияет ли на длину зубов и дозировака витамина C, и тип введения его в организм свинки (supp: апельсиновый сок (OJ) или аскорбиновая кислота (VC)).

Посмотрим на сочетание признаков

data(ToothGrowth)
ToothGrowth$dose <- as.factor(ToothGrowth$dose)
table(ToothGrowth$supp, ToothGrowth$dose)
##     
##      0.5  1  2
##   OJ  10 10 10
##   VC  10 10 10

И на средние значения зависимой переменной в этих группах

aggregate(ToothGrowth$len, by=list(ToothGrowth$supp, ToothGrowth$dose), FUN=mean)
##   Group.1 Group.2     x
## 1      OJ     0.5 13.23
## 2      VC     0.5  7.98
## 3      OJ       1 22.70
## 4      VC       1 16.77
## 5      OJ       2 26.06
## 6      VC       2 26.14

Закажем процедуру двухфакторного дисперсионного анализа

fit3 <- aov(len ~ dose + supp, data=ToothGrowth)
summary(fit3)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## dose         2 2426.4  1213.2   82.81  < 2e-16 ***
## supp         1  205.4   205.4   14.02 0.000429 ***
## Residuals   56  820.4    14.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Дисперсионный анализ (3x2 ANOVA) показал, что и дозировака витамина C влияет на длину зубов, F(2, 56) = 82.81, p < 0.001, и тип введения его в организм свинки, F(1, 56) = 14.02, p < 0.001.

Независимые переменные могут взаимодействовать дуг с другом, и их взаимодействие также может оказывать влияние на зависимую переменную. В предыдущий анализ кроме двух основных факторов можно добавить их взаимодействие. Давайте это сделаем.

fit4 <- aov(len ~ dose + supp + dose:supp, data=ToothGrowth) # взаимодействие факторов добавляется через с помощью добавление в формула нового члена, состоящего из имен двух взаимодействующих факторов, разделённых двоеточием - dose:supp
summary(fit4)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## dose         2 2426.4  1213.2  92.000  < 2e-16 ***
## supp         1  205.4   205.4  15.572 0.000231 ***
## dose:supp    2  108.3    54.2   4.107 0.021860 *  
## Residuals   54  712.1    13.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Можно было сделать короче: **len ~ dose*supp**. Эта формула добавляет сразу и две переменные по-отдельности, и их взаимодействие. Результат получится тот же самый.

fit5 <- aov(len ~ dose*supp, data=ToothGrowth) 
summary(fit5)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## dose         2 2426.4  1213.2  92.000  < 2e-16 ***
## supp         1  205.4   205.4  15.572 0.000231 ***
## dose:supp    2  108.3    54.2   4.107 0.021860 *  
## Residuals   54  712.1    13.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Дисперсионный анализ (3x2 ANOVA) показал, что есть как два основных эффекта дозы и типа введения по-отдельности, F(2, 54) = 92.00, p < 0.001, F(1, 54) = 15.57, p < 0.001, так и эффект их взаимодействия, F(2, 54) = 4.107, p = 0.022.

Посмотрим, какие именно группу отличаются друг от друга с помощью теста Тьюки.

TukeyHSD(fit5)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = len ~ dose * supp, data = ToothGrowth)
## 
## $dose
##         diff       lwr       upr   p adj
## 1-0.5  9.130  6.362488 11.897512 0.0e+00
## 2-0.5 15.495 12.727488 18.262512 0.0e+00
## 2-1    6.365  3.597488  9.132512 2.7e-06
## 
## $supp
##       diff       lwr       upr     p adj
## VC-OJ -3.7 -5.579828 -1.820172 0.0002312
## 
## $`dose:supp`
##                 diff        lwr         upr     p adj
## 1:OJ-0.5:OJ     9.47   4.671876  14.2681238 0.0000046
## 2:OJ-0.5:OJ    12.83   8.031876  17.6281238 0.0000000
## 0.5:VC-0.5:OJ  -5.25 -10.048124  -0.4518762 0.0242521
## 1:VC-0.5:OJ     3.54  -1.258124   8.3381238 0.2640208
## 2:VC-0.5:OJ    12.91   8.111876  17.7081238 0.0000000
## 2:OJ-1:OJ       3.36  -1.438124   8.1581238 0.3187361
## 0.5:VC-1:OJ   -14.72 -19.518124  -9.9218762 0.0000000
## 1:VC-1:OJ      -5.93 -10.728124  -1.1318762 0.0073930
## 2:VC-1:OJ       3.44  -1.358124   8.2381238 0.2936430
## 0.5:VC-2:OJ   -18.08 -22.878124 -13.2818762 0.0000000
## 1:VC-2:OJ      -9.29 -14.088124  -4.4918762 0.0000069
## 2:VC-2:OJ       0.08  -4.718124   4.8781238 1.0000000
## 1:VC-0.5:VC     8.79   3.991876  13.5881238 0.0000210
## 2:VC-0.5:VC    18.16  13.361876  22.9581238 0.0000000
## 2:VC-1:VC       9.37   4.571876  14.1681238 0.0000058

Визуализируем результаты №1 : Простой, но не очень хороший вариант

interaction.plot(ToothGrowth$dose, ToothGrowth$supp, ToothGrowth$len, type='b', col=c('red','blue'), pch=c(16, 18))

№2: Более сложный, но хороший вариант

pd = position_dodge(0.2)
ggplot(ToothGrowth, aes(as.numeric(dose), len, color = supp)) + 
    stat_summary(fun.data = mean_cl_boot, geom = 'errorbar', width = 0.2, lwd = 0.8, position = pd)+ 
    stat_summary(fun.data = mean_cl_boot, geom = 'line', size = 1.5, position = pd) +
    stat_summary(fun.data = mean_cl_boot, geom = 'point', size = 5, position = pd, pch=15) +
    theme_bw() + 
    xlab('Доза витамина C')+
    ylab('Длина зубов')

Size effect

Рассчитаем размер эффекта типа введение Витамина C на длину зубов. Воспользуемся пакетом lsr, в котором есть функция etaSquared, рассчитывающая размер эффекта eta-Squared.

etaSquared(fit5)
##               eta.sq eta.sq.part
## dose      0.70286419   0.7731092
## supp      0.05948365   0.2238254
## dose:supp 0.03137672   0.1320279

Pазмер эффекта дозы eta-Squared = 0.773. Это эффект cредней величины.

Анализ мощности

Теперь проверим, сколько необходимо людей, чтобы зафиксировать эффект такого размера. Воспользуется пакетом pwr, в котором есть функция pwr.anova.test (есть и для других дизайнов).

pwr.anova.test(k = 2, f= 0.773, sig.level = 0.05, power = 0.95)
## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 2
##               n = 11.92686
##               f = 0.773
##       sig.level = 0.05
##           power = 0.95
## 
## NOTE: n is number in each group

Анализ мощности показал, что для того, чтобы с 95% мощностью зафиксировать эффект дозы (eta-Squared = 0.773), необходимо всего по 12 человек на группу.

Дисперсионный анализ с повторными измерениями (ANOVA with repeated measures)

Пример из книги R в действии. Набор данных CO2 содержит результаты исследования холодоустойчивости северных и южных популяций злака ежовника. Сравнивали интенсивность фотосинтеза охлажденных и неохлажденных растений при разных концентрациях углекислого газа в окружающей среде. Половина растений происходила из Квебека, а половина – из штата Миссисипи.

Plant - a unique identifier for each plant.

Type - the origin of the plant / штат происхождения: Квебек или Миссисипи

Treatment - nonchilled chilled / охлаждённый неохлаждённый

conc - ambient carbon dioxide concentrations / концентрация углекислого газа в окружающей среде

uptake - carbon dioxide uptake rates / уровень потребления углекислого газа

Посмотрим на структуру данных

str(CO2)
## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame':   84 obs. of  5 variables:
##  $ Plant    : Ord.factor w/ 12 levels "Qn1"<"Qn2"<"Qn3"<..: 1 1 1 1 1 1 1 2 2 2 ...
##  $ Type     : Factor w/ 2 levels "Quebec","Mississippi": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Treatment: Factor w/ 2 levels "nonchilled","chilled": 1 1 1 1 1 1 1 1 1 1 ...
##  $ conc     : num  95 175 250 350 500 675 1000 95 175 250 ...
##  $ uptake   : num  16 30.4 34.8 37.2 35.3 39.2 39.7 13.6 27.3 37.1 ...
##  - attr(*, "formula")=Class 'formula'  language uptake ~ conc | Plant
##   .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv> 
##  - attr(*, "outer")=Class 'formula'  language ~Treatment * Type
##   .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv> 
##  - attr(*, "labels")=List of 2
##   ..$ x: chr "Ambient carbon dioxide concentration"
##   ..$ y: chr "CO2 uptake rate"
##  - attr(*, "units")=List of 2
##   ..$ x: chr "(uL/L)"
##   ..$ y: chr "(umol/m^2 s)"

Отбираем только охлажденные растения.

w1b1 <- subset(CO2, Treatment=='chilled')
fit6 <- aov(uptake ~ conc*Type + Error(Plant/conc), data = w1b1)
summary(fit6)
## 
## Error: Plant
##           Df Sum Sq Mean Sq F value  Pr(>F)   
## Type       1 2667.2  2667.2   60.41 0.00148 **
## Residuals  4  176.6    44.1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Plant:conc
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## conc       1  888.6   888.6  215.46 0.000125 ***
## conc:Type  1  239.2   239.2   58.01 0.001595 ** 
## Residuals  4   16.5     4.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 30    869   28.97

Результаты дисперсионного анализа, представленные в таблице, свидетельствуют о том, что главные эффекты (штат и концентрация), а также взаимодействие между ними (все p-value < 0.05).

Визуализируем

boxplot(uptake ~ Type*conc, data=w1b1, col=(c('gold', 'green')), las = 2, cex.axis=0.6)