ANOVA

База данных со случаями возникновением рака молочной железы. Включает следующие переменные:

• country (страна)
• continent (континент)
• year (год)
• lifeExp (ожидаемая продолдительность жизни)
• pop (население)
• gdpPercap (ВВП на душу населения)
• NewCasesOfBreastCancerIn2002 (новые случаи рака)
• AlcoholComsumption (потребление алкоголя)
• BloodPressure (кровяное давление)
• BodyMassIndex (индекс массы тела)
• Cholestorol (холестерин)
• Smoking (количество курящих).

База данных доступна по ссылке :https://github.com/bioinformatics-core-shared-training/intermediate-stats/blob/master/globalBreastCancerRisk.csv.

Необходимо определить, есть ли различия между континентами в количестве новых случаях рака.

data <- read.csv('globalBreastCancerRisk.csv')
summary(data)

##    country           continent              year         lifeExp     
##  Length:123         Length:123         Min.   :2002   Min.   :39.19  
##  Class :character   Class :character   1st Qu.:2002   1st Qu.:54.87  
##  Mode  :character   Mode  :character   Median :2002   Median :70.75  
##                                        Mean   :2002   Mean   :65.53  
##                                        3rd Qu.:2002   3rd Qu.:75.10  
##                                        Max.   :2002   Max.   :82.00  
##                                                                      
##       pop              gdpPercap       NewCasesOfBreastCancerIn2002
##  Min.   :2.880e+05   Min.   :  446.4   Min.   :  3.90              
##  1st Qu.:4.448e+06   1st Qu.: 1330.2   1st Qu.: 19.50              
##  Median :1.060e+07   Median : 5351.6   Median : 28.10              
##  Mean   :4.554e+07   Mean   :10117.2   Mean   : 37.94              
##  3rd Qu.:3.123e+07   3rd Qu.:12262.0   3rd Qu.: 50.60              
##  Max.   :1.280e+09   Max.   :44684.0   Max.   :101.10              
##                                                                    
##  AlcoholComsumption BloodPressure   BodyMassIndex    Cholestorol   
##  Min.   : 0.020     Min.   :117.7   Min.   :19.72   Min.   :4.089  
##  1st Qu.: 2.070     1st Qu.:123.9   1st Qu.:22.67   1st Qu.:4.393  
##  Median : 5.550     Median :126.5   Median :25.12   Median :4.825  
##  Mean   : 6.041     Mean   :126.8   Mean   :24.69   Mean   :4.780  
##  3rd Qu.: 9.545     3rd Qu.:129.7   3rd Qu.:26.30   3rd Qu.:5.103  
##  Max.   :16.270     Max.   :134.2   Max.   :30.34   Max.   :5.658  
##                                                                    
##     Smoking     
##  Min.   : 0.30  
##  1st Qu.: 3.40  
##  Median : 9.10  
##  Mean   :13.81  
##  3rd Qu.:25.80  
##  Max.   :40.10  
##  NA's   :30

boxplot(data$NewCasesOfBreastCancerIn2002 ~ data$continent, las = 2)

library(car)

## Загрузка требуемого пакета: carData

Делаем тест, показывающий равны ли дисперсии.

leveneTest(data$NewCasesOfBreastCancerIn2002 ~ data$continent)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value   Pr(>F)   
## group   4  4.8828 0.001116 **
##       118                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Дисперсии не равны, так как Pr < 0.05 (а это значит, что мы можем отвергнуть нулевую гипотезу о равенстве дисперсий)

oneway.test(data$NewCasesOfBreastCancerIn2002 ~ data$continent, var.equal = F) 

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  data$NewCasesOfBreastCancerIn2002 and data$continent
## F = 72.151, num df = 4.0000, denom df = 7.7944, p-value = 3.333e-06

Есть разница между средними в группах, так как p-value < 0.05.

Посмотрим, между какими группами наблюдается разница.

Так как дисперсии не равны, используем попарный t-test с поправкой Бонферонни.

pairwise.t.test(data$NewCasesOfBreastCancerIn2002, data$continent, 
                adjust = "bonferroni")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  data$NewCasesOfBreastCancerIn2002 and data$continent 
## 
##          Africa  Americas Asia    Europe 
## Americas 7.5e-05 -        -       -      
## Asia     0.05475 0.10080  -       -      
## Europe   < 2e-16 7.2e-08  8.9e-13 -      
## Oceania  7.5e-07 0.00047  3.2e-05 0.11521
## 
## P value adjustment method: holm

Видно, что значение P value < 0.05 наблюдаются в след. группах: Africa - Americas, Europe, Oceania; Americas - Europe, Oceania; Asia - Europe, Oceania; Именно между этими группами есть значимая разница.

Посмотрим какая именно, визуализировав медианы в группах.

library(sjPlot)

## Learn more about sjPlot with 'browseVignettes("sjPlot")'.

plot_grpfrq(data$NewCasesOfBreastCancerIn2002, data$continent,  type = "box")

## Warning in plot_grpfrq(data$NewCasesOfBreastCancerIn2002, data$continent, : в
## результате преобразования созданы NA

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique

## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique

## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique

## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique

Видно, что, например, в Европе и Океании больше случаев возникновения рака, чем в Африке, Америке и Азии.

library(sjstats)

omega_sq(aov(data$NewCasesOfBreastCancerIn2002 ~ data$continent))

## Warning: 'omega_sq' устарел
## Используйте 'effectsize::omega_sqared()'.
## См. help("Deprecated")

## For one-way between subjects designs, partial omega squared is equivalent to omega squared.
## Returning omega squared.

## # Effect Size for ANOVA
## 
## Parameter      | Omega2 |       95% CI
## --------------------------------------
## data$continent |   0.57 | [0.47, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

eta_sq(aov(data$NewCasesOfBreastCancerIn2002 ~ data$continent))

## Warning: 'eta_sq' устарел
## Используйте 'effectsize::eta_squared()'.
## См. help("Deprecated")

## For one-way between subjects designs, partial eta squared is equivalent to eta squared.
## Returning eta squared.

## # Effect Size for ANOVA
## 
## Parameter      | Eta2 |       95% CI
## ------------------------------------
## data$continent | 0.59 | [0.49, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

Размер эффекта между группами (то есть разница) очень большой, так как Omega2, Eta2 > 0.14.

И давайте повторим все тоже самое, но при помощи одного графика

library(ggstatsplot)

## You can cite this package as:
##      Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
##      Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167

ggbetweenstats(
  data = data,
  x = continent,
  y = NewCasesOfBreastCancerIn2002,
  xlab = "Continent",
  ylab = "New Cases Of Breast Cancer In 2002",
  plot.type = "box")

Пример 2

База данных с обзором вин, содержащая следующие переменные:

• country (страна производства вина)
• description (описание вина)
• designation (происхождение винограда)
• points (количество баллов WineEnthusiast по шкале 1-100)
• price (стоимость бутылки вина)
• province (штат или провинция производства вина) и др.

(Ссылка: https://www.kaggle.com/datasets/zynicide/wine-reviews?resource=download)

data2 <- read.csv('wine.csv')

table(data2$country)

## 
##       Chile     Germany New Zealand 
##        1485         756         487

boxplot(data2$price ~ data2$country)

Делаем тест, показывающий равны ли дисперсии.

leveneTest(data2$price ~ data2$country)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    2   56.65 < 2.2e-16 ***
##       2680                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Дисперсии не равны, так как Pr < 0.05

oneway.test(data2$price ~ data2$country, var.equal = F)

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  data2$price and data2$country
## F = 57.227, num df = 2.0, denom df = 1156.8, p-value < 2.2e-16

Есть разница между средними в группах, так как p-value < 0.05.

Посмотрим, между какими группами наблюдается разница.

Так как дисперсии не равны, используем попарный t-test с поправкой Бонферонни.

pairwise.t.test(data2$price, data2$country, 
                adjust = "bonferroni")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  data2$price and data2$country 
## 
##             Chile   Germany
## Germany     < 2e-16 -      
## New Zealand 0.004   3.7e-14
## 
## P value adjustment method: holm

Видно, что между всеми группами есть разница.

plot_grpfrq(data2$price, data2$country,  type = "box")

## Warning in plot_grpfrq(data2$price, data2$country, type = "box"): в результате
## преобразования созданы NA

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique

## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique

По рисунку можно сказать, что в Германии цены больше, чем в двух других странах. Но в Новой Зеландии цены выше, чем в Чили.

omega_sq(aov(data2$price ~ data2$country))

## Warning: 'omega_sq' устарел
## Используйте 'effectsize::omega_sqared()'.
## См. help("Deprecated")

## For one-way between subjects designs, partial omega squared is equivalent to omega squared.
## Returning omega squared.

## # Effect Size for ANOVA
## 
## Parameter     | Omega2 |       95% CI
## -------------------------------------
## data2$country |   0.06 | [0.05, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

eta_sq(aov(data2$price ~ data2$country))

## Warning: 'eta_sq' устарел
## Используйте 'effectsize::eta_squared()'.
## См. help("Deprecated")

## For one-way between subjects designs, partial eta squared is equivalent to eta squared.
## Returning eta squared.

## # Effect Size for ANOVA
## 
## Parameter     | Eta2 |       95% CI
## -----------------------------------
## data2$country | 0.06 | [0.05, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

Размер эффекта между группами (то есть разница) средний, так как Omega2, Eta2 в районе 0.06 - 0.14.

И опять график при помощи функции ggbetweenstats()

ggbetweenstats(
  data = data2,
  x = country,
  y = price,
  xlab = "country",
  ylab = "price",
  plot.type = "box")

По боксплотам было видно, что распределение по группам далеко от нормального, это значит, что можно попробовать либо преобразовать числовую переменную, либо использовать непараметрический тест. Давайте попробуем оба варианта.

hist(data2$price)

Действительно, как и у многих переменных, касающихся дохода или цены, наблюдается длинный хвост справа.

Вариант 1, логарифмируем.

data2$logprice <- log(data2$price)
hist(data2$logprice)

Так явно лучше.

Делаем тест, показывающий равны ли дисперсии.

leveneTest(data2$logprice ~ data2$country)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    2  43.576 < 2.2e-16 ***
##       2680                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Дисперсии не равны, так как Pr < 0.05

oneway.test(data2$logprice ~ data2$country, var.equal = F)

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  data2$logprice and data2$country
## F = 201.38, num df = 2.0, denom df = 1177.7, p-value < 2.2e-16

Есть разница между средними в группах, так как p-value < 0.05.

Посмотрим, между какими группами наблюдается разница.

Так как дисперсии не равны, используем попарный t-test с поправкой Бонферонни.

pairwise.t.test(data2$logprice, data2$country, 
                adjust = "bonferroni")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  data2$logprice and data2$country 
## 
##             Chile   Germany
## Germany     < 2e-16 -      
## New Zealand < 2e-16 9.1e-13
## 
## P value adjustment method: holm

Видно, что между всеми группами есть разница.

plot_grpfrq(data2$logprice, data2$country,  type = "box")

## Warning in plot_grpfrq(data2$logprice, data2$country, type = "box"): в
## результате преобразования созданы NA

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique

## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique

По рисунку можно сказать, что в Германии цены больше, чем в двух других странах. Но в Новой Зеландии цены выше, чем в Чили.

omega_sq(aov(data2$logprice ~ data2$country))

## Warning: 'omega_sq' устарел
## Используйте 'effectsize::omega_sqared()'.
## См. help("Deprecated")

## For one-way between subjects designs, partial omega squared is equivalent to omega squared.
## Returning omega squared.

## # Effect Size for ANOVA
## 
## Parameter     | Omega2 |       95% CI
## -------------------------------------
## data2$country |   0.15 | [0.13, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

eta_sq(aov(data2$logprice ~ data2$country))

## Warning: 'eta_sq' устарел
## Используйте 'effectsize::eta_squared()'.
## См. help("Deprecated")

## For one-way between subjects designs, partial eta squared is equivalent to eta squared.
## Returning eta squared.

## # Effect Size for ANOVA
## 
## Parameter     | Eta2 |       95% CI
## -----------------------------------
## data2$country | 0.15 | [0.13, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

Размер эффекта между группами (то есть разница) большая, так как \(\omega^2\) и \(\eta^2\) больше 0.14.

И опять график при помощи функции ggbetweenstats()

ggbetweenstats(
  data = data2,
  x = country,
  y = logprice,
  xlab = "country",
  ylab = "price",
  plot.type = "box")

Вариант 2. Непараметрический тест

kruskal.test(data2$logprice ~ data2$country)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  data2$logprice by data2$country
## Kruskal-Wallis chi-squared = 461.48, df = 2, p-value < 2.2e-16

p-value < 0.05, а это значит, что средние ранги различаются значимо.

Делаем post hoc анализ (чтобы сделать вывод, какие группы различаются между собой)

library (FSA)

## Registered S3 methods overwritten by 'FSA':
##   method       from
##   confint.boot car 
##   hist.boot    car

## ## FSA v0.9.3. See citation('FSA') if used in publication.
## ## Run fishR() for related website and fishR('IFAR') for related book.

## 
## Присоединяю пакет: 'FSA'

## Следующий объект скрыт от 'package:sjstats':
## 
##     se

## Следующий объект скрыт от 'package:car':
## 
##     bootCase

dunnTest(data2$logprice, data2$country, method = "bonferroni" )

## Warning: 'g' variable was coerced to a factor.

## Warning: Some rows deleted from 'x' and 'g' because missing data.

## Dunn (1964) Kruskal-Wallis multiple comparison

##   p-values adjusted with the Bonferroni method.

##              Comparison          Z      P.unadj        P.adj
## 1       Chile - Germany -20.259818 2.910347e-91 8.731040e-91
## 2   Chile - New Zealand -12.626244 1.513301e-36 4.539904e-36
## 3 Germany - New Zealand   4.138874 3.490143e-05 1.047043e-04

ggbetweenstats(
  data = data2,
  x = country,
  y = logprice,
  xlab = "Country",
  ylab = "Log(price)",
  plot.type = "box",
  type = "np",
  conf.level = 0.95)

ggbetweenstats() так же подходит и для визуализации непараметрических тестов

ggbetweenstats(
  data = data2,
  x = country,
  y = logprice,
  xlab = "country",
  ylab = "price",
  type = "np",
  plot.type = "box", )