Manova

Czasem mamy do czynienia z więcej niż jedną zmienną odpowiedzi równocześnie. Jeśli podejrzewamy (na podstawie innych testów), że zmienne te są niezależne od siebie to możemy wykonać zwykłą analizę wariancji dla sumy tych zmiennych wynikowych. Jeśli jednak te zmienne są ze sobą powiązane lub po prostu zmienna wynikowa ma charakter wielowymiarowy to należy skorzystać z wielowymiarowej analizy wariancji.

Data structure

Sprawdzimy czy istnieje istotna różnica w długościach płatków i działek pomiędzy różnymi gatunkami irysów.

irysy <- iris
cor(irysy[,-5])

##              Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length    1.0000000  -0.1175698    0.8717538   0.8179411
## Sepal.Width    -0.1175698   1.0000000   -0.4284401  -0.3661259
## Petal.Length    0.8717538  -0.4284401    1.0000000   0.9628654
## Petal.Width     0.8179411  -0.3661259    0.9628654   1.0000000

Sample size assumption

irysy %>%
  group_by(Species) %>%
  summarise(N = n())

## # A tibble: 3 × 2
##   Species        N
##   <fct>      <int>
## 1 setosa        50
## 2 versicolor    50
## 3 virginica     50

Ponieważ powyższa tabela zawiera 50 obserwacji na grupę, założenie o odpowiedniej liczebności próby jest spełnione.

Identyfikacja wartości odstających

Odchylenia mogą być łatwo zidentyfikowane przy użyciu metod wykresów pudełkowych, zaimplementowanych w funkcji R identify_outliers() [pakiet rstatix].

Pogrupuj dane według gatunków, a następnie zidentyfikuj wartości odstające w zmiennej Sepal.Length:

irysy %>%
  group_by(Species) %>%
  identify_outliers(Sepal.Length)

## # A tibble: 1 × 7
##   Species   Sepal.Length Sepal.Width Petal.Length Petal.Width is.outlier
##   <fct>            <dbl>       <dbl>        <dbl>       <dbl> <lgl>     
## 1 virginica          4.9         2.5          4.5         1.7 TRUE      
## # … with 1 more variable: is.extreme <lgl>

Pogrupuj dane według gatunków, a następnie zidentyfikuj wartości odstające w zmiennej Petal.Length:

irysy %>%
  group_by(Species) %>%
  identify_outliers(Petal.Length)

## # A tibble: 5 × 7
##   Species    Sepal.Length Sepal.Width Petal.Length Petal.Width is.outlier
##   <fct>             <dbl>       <dbl>        <dbl>       <dbl> <lgl>     
## 1 setosa              4.3         3            1.1         0.1 TRUE      
## 2 setosa              4.6         3.6          1           0.2 TRUE      
## 3 setosa              4.8         3.4          1.9         0.2 TRUE      
## 4 setosa              5.1         3.8          1.9         0.4 TRUE      
## 5 versicolor          5.1         2.5          3           1.1 TRUE      
## # … with 1 more variable: is.extreme <lgl>

W zmiennej Sepal.Length i Petal.length nie było żadnych skrajnych wartości jednoczynnikowych, co oceniono za pomocą metod wykresów pudełkowych.

Zauważ, że w sytuacji występowania skrajnych wartości odstających może to być spowodowane: 1) błędów we wprowadzaniu danych, błędów pomiarowych lub nietypowych wartości.

Możesz mimo wszystko włączyć wartość odstającą do analizy, jeśli uważasz, że nie wpłynie ona w znaczący sposób na wynik. Można to ocenić poprzez porównanie wyniku MANOVA z i bez wartości odstającej.

Pamiętaj, aby w części poświęconej wynikom pisemnym zamieścić informacje o wszelkich decyzjach podjętych w związku ze znalezionymi wartościami odstającymi.

Wykrywanie wielowymiarowych wartości odstających

Wielowymiarowe wartości odstające to punkty danych, które mają nietypową kombinację wartości zmiennych wynikowych (lub zależnych).

W MANOVA, odległość Mahalanobisa jest zwykle używana do wykrywania odstających punktów wieloczynnikowych. Odległość ta mówi nam, jak daleko obserwacja znajduje się od środka chmury, biorąc pod uwagę kształt (kowariancję) chmury.

Funkcja mahalanobis_distance() [pakiet rstatix] może być łatwo użyta do obliczenia odległości Mahalanobisa i oznaczenia wielowymiarowych wartości odstających. Więcej informacji można znaleźć w dokumentacji funkcji.

Metryka ta musi być obliczana przez grupy:

# Compute distance by groups and filter outliers
# Użyj -id aby pominąć kolumnę id w obliczeniach
irysy %>%
 group_by(Species) %>%
 mahalanobis_distance() %>%
 filter(is.outlier == TRUE) %>%
  as.data.frame()

## [1] Sepal.Length Sepal.Width  Petal.Length Petal.Width  mahal.dist  
## [6] is.outlier  
## <0 rows> (or 0-length row.names)

W danych nie występowały wielowymiarowe wartości odstające, co oceniono na podstawie odległości Mahalanobisa (p ≤ 0,001).

Jeśli masz wielowymiarowe wartości odstające, możesz rozważyć przeprowadzenie MANOVA przed i po usunięciu wartości odstających, aby sprawdzić, czy ich obecność zmienia wyniki.

Normalność 1-wymiarowa

Założenie normalności może być sprawdzone przez obliczenie testu Shapiro-Wilka dla każdej zmiennej zależnej na każdym poziomie zmiennej grupującej. Jeżeli dane mają rozkład normalny, wartość p powinna być większa niż 0,05.

irysy %>%
  group_by(Species) %>%
  shapiro_test(Sepal.Length, Petal.Length) %>%
  arrange(variable)

## # A tibble: 6 × 4
##   Species    variable     statistic      p
##   <fct>      <chr>            <dbl>  <dbl>
## 1 setosa     Petal.Length     0.955 0.0548
## 2 versicolor Petal.Length     0.966 0.158 
## 3 virginica  Petal.Length     0.962 0.110 
## 4 setosa     Sepal.Length     0.978 0.460 
## 5 versicolor Sepal.Length     0.978 0.465 
## 6 virginica  Sepal.Length     0.971 0.258

Sepal.Length i Petal.length były normalnie rozłożone dla każdej grupy gatunków, jak oceniono testem Shapiro-Wilka (p ≤ 0.05).

Można również utworzyć wykres QQ dla każdej grupy. Wykres QQ rysuje korelację pomiędzy danymi a rozkładem normalnym.

# QQ plot of Sepal.Length
ggqqplot(irysy, "Sepal.Length", facet.by = "Species",
         ylab = "Sepal Length", ggtheme = theme_bw())

# QQ plot of Petal.Length
ggqqplot(irysy, "Petal.Length", facet.by = "Species",
         ylab = "Petal Length", ggtheme = theme_bw())

Wszystkie punkty wypadają w przybliżeniu wzdłuż linii odniesienia, dla każdej grupy. Możemy więc założyć normalność danych.

Zauważ, że jeśli liczebność próby jest większa niż 50, preferowany jest normalny wykres QQ, ponieważ przy większych liczebnościach test Shapiro-Wilka staje się bardzo wrażliwy nawet na niewielkie odchylenia od normalności.

W sytuacji, gdy założenia nie są spełnione, można rozważyć przeprowadzenie MANOVA na danych po przekształceniu zmiennych wynikowych. Można również wykonać test niezależnie od tego, ponieważ MANOVA jest dość odporna na odchylenia od normalności.

Wielowymiarowa normalność:

irysy %>%
  select(Sepal.Length, Petal.Length) %>%
  mshapiro_test()

## # A tibble: 1 × 2
##   statistic p.value
##       <dbl>   <dbl>
## 1     0.995   0.855

Test nie jest istotny (p ≤ 0,05), więc możemy założyć normalność wielowymiarową.

Określenie współliniowości

Idealnie korelacja pomiędzy zmiennymi wynikowymi powinna być umiarkowana, nie za wysoka. Korelacja powyżej 0,9 jest wskaźnikiem wieloliniowości, która jest problematyczna dla MANOVA.

Z drugiej strony, jeśli korelacja jest zbyt niska, należy rozważyć przeprowadzenie oddzielnej jednokierunkowej ANOVA dla każdej zmiennej wynikowej.

Oblicz parami współczynniki korelacji Pearsona pomiędzy zmiennymi wynikowymi. W poniższym kodzie R użyjemy funkcji cor_test() [pakiet rstatix]. Jeśli masz więcej niż dwie zmienne wynikowe, rozważ użycie funkcji cor_mat():

irysy %>% cor_test(Sepal.Length, Petal.Length)

## # A tibble: 1 × 8
##   var1         var2           cor statistic        p conf.low conf.high method 
##   <chr>        <chr>        <dbl>     <dbl>    <dbl>    <dbl>     <dbl> <chr>  
## 1 Sepal.Length Petal.Length  0.87      21.6 1.04e-47    0.827     0.906 Pearson

Nie stwierdzono współliniowości, co oceniono na podstawie korelacji Pearsona (r = 0,87, p = 0,0001).

W sytuacji, gdy mamy do czynienia z współliniowością, można rozważyć usunięcie jednej ze zmiennych zależnych, która jest silnie skorelowana.

Liniowość

Zależność między parami zmiennych wynikowych powinna być liniowa dla każdej grupy. Można to sprawdzić wizualnie, tworząc macierz wykresu rozrzutu za pomocą funkcji R ggpairs() [pakiet GGally]. W naszym przykładzie mamy tylko jedną parę:

# Tworzenie macierzy rozrzutu według grup
library(GGally)
results <- irysy %>%
  select(Sepal.Length, Petal.Length, Species) %>%
  group_by(Species) %>%
  doo(~ggpairs(.) + theme_bw(), result = "plots")
results$plots

## [[1]]

## 
## [[2]]

## 
## [[3]]

Istnieje liniowa zależność pomiędzy długością działek i długością płatków w każdej grupie gatunków, co zostało ocenione na wykresie rozrzutu.

W sytuacji, gdy wykryjesz nieliniowe zależności, możesz: - przekształcić lub usunąć zmienne wynikowe, których to dotyczy; - przeprowadzić analizę mimo wszystko. W ten sposób stracisz trochę mocy.

Wielowymiarowa jednorodność wariancji

Można to ocenić za pomocą testu Box’s M-test zaimplementowanego w pakiecie rstatix.

box_m(irysy[, c("Sepal.Length", "Petal.Length")], irysy$Species)

## # A tibble: 1 × 4
##   statistic  p.value parameter method                                           
##       <dbl>    <dbl>     <dbl> <chr>                                            
## 1      58.4 9.62e-11         6 Box's M-test for Homogeneity of Covariance Matri…

Test jest statystycznie istotny (tzn. p = 0,001), więc dane naruszyły założenie o jednorodności macierzy wariancji-kowariancji.

Zauważ, że jeśli masz zrównoważony projekt (tzn. grupy o podobnych rozmiarach), nie musisz się zbytnio przejmować naruszeniem założenia o jednorodności macierzy wariancji-kowariancji i możesz kontynuować swoją analizę.

Jednorodność wariancji

Dla każdej ze zmiennych wynikowych, jednokierunkowa MANOVA zakłada, że istnieją równe wariancje między grupami. Można to sprawdzić za pomocą testu równości wariancji Levene’a. Kluczowa funkcja R: levene_test() [pakiet rstatix].

Procedura:

1. Zbierz zmienne zależne w pary klucz-wartość.
2. Pogrupuj według zmiennych
3. Oblicz test Levene’a

irysy %>% 
  gather(key = "variable", value = "value", Sepal.Length, Petal.Length) %>%
  group_by(variable) %>%
  levene_test(value ~ Species)

## # A tibble: 2 × 5
##   variable       df1   df2 statistic            p
##   <chr>        <int> <int>     <dbl>        <dbl>
## 1 Petal.Length     2   147     19.5  0.0000000313
## 2 Sepal.Length     2   147      6.35 0.00226

Test Levene’a jest istotny (p ≤ 0,05), więc nie było jednorodności wariancji.

Zauważ, że jeśli nie ma homogeniczności wariancji, możesz spróbować przekształcić zmienną wynikową (zależną), aby skorygować nierówne wariancje.

MANOVA

Najczęściej zalecaną statystyką wielowymiarową jest Lambda Wilksa.

Jednakże, statystyka Pillai’s jest bardziej odporna i jest zalecana, gdy masz niezrównoważoną konstrukcję i statystycznie istotny wynik Box’s M.

Zauważ, że “Pillai” jest wartością domyślną w funkcji R Manova() [pakiet ‘car’].

model <- lm(cbind(Sepal.Length, Petal.Length) ~ Species, irysy)
Manova(model, test.statistic = "Pillai")

## 
## Type II MANOVA Tests: Pillai test statistic
##         Df test stat approx F num Df den Df    Pr(>F)    
## Species  2    0.9885   71.829      4    294 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Stwierdzono statystycznie istotną różnicę pomiędzy gatunkami na połączonych zmiennych zależnych (długości płatków i działek irysów), F(4, 294) = 71.829, p ≤ 0.0001. Istotnie, tak sformułowana zmienna dwuwymiarowa posiada różne wartości średnie w zależności od gatunku kwiatu.

Testy post-hoc

Statystycznie znacząca jednokierunkowa MANOVA może być kontynuowana przez jednokierunkową ANOVA badającą, oddzielnie, każdą zmienną zależną. Celem jest zidentyfikowanie specyficznych zmiennych zależnych, które przyczyniły się do znaczącego efektu globalnego.

Procedura:

1. Zebrać zmienne zależne w pary klucz-wartość.
2. Pogrupuj według zmiennych
3. Oblicz jednokierunkowy test ANOVA

Poniższe kody R pokazują, jak używać każdej z powszechnie używanych funkcji:

# pogrupuj po zmiennej
grouped.data <- irysy %>%
  gather(key = "variable", value = "value", Sepal.Length, Petal.Length) %>%
  group_by(variable)

# Do welch one way anova test
grouped.data %>% welch_anova_test(value ~ Species)

## # A tibble: 2 × 8
##   variable     .y.       n statistic   DFn   DFd        p method     
## * <chr>        <chr> <int>     <dbl> <dbl> <dbl>    <dbl> <chr>      
## 1 Petal.Length value   150     1828.     2  78.1 2.69e-66 Welch ANOVA
## 2 Sepal.Length value   150      139.     2  92.2 1.51e-28 Welch ANOVA

# or do Kruskal-Wallis test
grouped.data %>% kruskal_test(value ~ Species)

## # A tibble: 2 × 7
##   variable     .y.       n statistic    df        p method        
## * <chr>        <chr> <int>     <dbl> <int>    <dbl> <chr>         
## 1 Petal.Length value   150     130.      2 4.8 e-29 Kruskal-Wallis
## 2 Sepal.Length value   150      96.9     2 8.92e-22 Kruskal-Wallis

# or use aov()
grouped.data %>% anova_test(value ~ Species)

## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()

## # A tibble: 2 × 8
##   variable     Effect    DFn   DFd     F        p `p<.05`   ges
## * <chr>        <chr>   <dbl> <dbl> <dbl>    <dbl> <chr>   <dbl>
## 1 Petal.Length Species     2   147 1180. 2.86e-91 *       0.941
## 2 Sepal.Length Species     2   147  119. 1.67e-31 *       0.619

Stwierdzono statystycznie istotną różnicę w długości płatków (F(2, 147) = 119, p ≤ 0.0001 ) i długości płatków (F(2, 147) = 1180, p ≤ 0.0001 ) pomiędzy gatunkami tęczówek.

Porównania wielokrotne

Statystycznie istotna jednoczynnikowa ANOVA może być kontynuowana przez wielokrotne porównania parami w celu określenia, które grupy są różne.

Funkcje R tukey_hsd() [pakiet rstatix] mogą być użyte do obliczenia testów post-hoc Tukeya, jeśli założenie o jednorodności wariancji jest spełnione.

Jeśli założenie o jednorodności wariancji zostało naruszone, tak jak w naszym przykładzie, można użyć testu post-hoc Games-Howella. Możliwe jest również użycie funkcji pairwise_test() [rstatix] z opcją pool.sd = FALSE i var.equal = FALSE.

pwc <- irysy %>%
  gather(key = "variables", value = "value", Sepal.Length, Petal.Length) %>%
  group_by(variables) %>%
  games_howell_test(value ~ Species) %>%
  select(-estimate, -conf.low, -conf.high) # Remove details
pwc

## # A tibble: 6 × 6
##   variables    .y.   group1     group2        p.adj p.adj.signif
##   <chr>        <chr> <chr>      <chr>         <dbl> <chr>       
## 1 Petal.Length value setosa     versicolor 1.85e-11 ****        
## 2 Petal.Length value setosa     virginica  1.68e-11 ****        
## 3 Petal.Length value versicolor virginica  4.45e-10 ****        
## 4 Sepal.Length value setosa     versicolor 2.86e-10 ****        
## 5 Sepal.Length value setosa     virginica  0        ****        
## 6 Sepal.Length value versicolor virginica  5.58e- 7 ****

Wszystkie porównania parami były istotne dla każdej ze zmiennych wynikowych (długość i szerokość płatka).

Raport końcowy

Jednoczynnikowa analiza wariancji została przeprowadzona w celu określenia wpływu gatunku irysa na długość działek i płatków. Istnieją trzy różne gatunki: setosa, versicolor i virginica.

Wystąpiła statystycznie istotna różnica pomiędzy gatunkami na połączonych zmiennych zależnych (Sepal.Length i Petal.Length), F(4, 294) = 71.829, p = 0.0001.

Kolejne jednoczynnikowe ANOVA, używając skorygowanego poziomu alfa Bonferroniego 0.025, wykazały, że istnieje statystycznie istotna różnica w Sepal.Length (F(2, 147) = 119, p ≤ 0,0001 ) i Petal.Length (F(2, 147) = 1180, p ≤ 0,0001 ) między gatunkami irysów.

Wszystkie porównania parami między grupami były znaczące dla każdej zmiennej zależnej (Sepal.Length i Petal.Length).

# Wizualizacja: wykresy pudełkowe z wartościami p
pwc <- pwc %>% add_xy_position(x = "Species")
test.label <- create_test_label(
  description = "MANOVA", statistic.text = quote(italic("F")),
  statistic = 71.83, p= "<0.0001", parameter = "4,294",
  type = "expression", detailed = TRUE
  )
ggboxplot(
  irysy, x = "Species", y = c("Sepal.Length", "Petal.Length"), 
  merge = TRUE, palette = "jco"
  ) + 
  stat_pvalue_manual(
    pwc, hide.ns = TRUE, tip.length = 0, 
    step.increase = 0.1, step.group.by = "variables",
    color = "variables"
    ) +
  labs(
    subtitle = test.label,
    caption = get_pwc_label(pwc, type = "expression")
  )