ANOVA, ANCOVA, MANOVA
Julianna Waluś
31 01 2022
Wczytaj zbiór danych Salaries z biblioteki carData i zapoznaj się z nim. Potraktuj zmienną salary jako odpowiedź a zmienne rank, discipline, sex jako zmienne wyjaśniające. Przy pomocy funkcji table sprawdź czy grupy są równoliczne. Wykonaj dwa wykresy skrzynkowe, jeden dla mężczyzn a drugi dla kobiet.
Porównaj ze sobą wyniki trójczynnikowych analiz wariancji wykonanych bez interakcji oraz z interakcjami. Zinterpretuj wyniki. Czy lata od doktoratu (yrs.since.phd), staż pracy (yrs.service) mogą być istotnymi zmiennymi współ-oddziaływującymi na płacę (salary)? Czy istnieje istotna różnica w latach od doktoratu (yrs.since.phd) i stażu pracy (yrs.service) różnej rangi profesorów?
ANOVA (dla ‘salary’ względem ‘rank’, ‘sex’, ‘discipline’)
Hipotezy
H0- płeć, dyscyplina oraz tytuł nie mają znaczący wpływ na pensje osoby. H1- płeć, dyscyplina oraz tytuł mają znaczący wpływ na pensje osoby. # Równoliczność grup Pierwszym krokiem jest sprawdzenie, czy wszysykie zmienne posiadają taką samą ilość danych
describe(Salaries)## vars n mean sd median trimmed mad min
## rank* 1 397 2.50 0.77 3 2.62 0.00 1
## discipline* 2 397 1.54 0.50 2 1.55 0.00 1
## yrs.since.phd 3 397 22.31 12.89 21 21.83 14.83 1
## yrs.service 4 397 17.61 13.01 16 16.51 14.83 0
## sex* 5 397 1.90 0.30 2 2.00 0.00 1
## salary 6 397 113706.46 30289.04 107300 111401.61 29355.48 57800
## max range skew kurtosis se
## rank* 3 2 -1.12 -0.38 0.04
## discipline* 2 1 -0.18 -1.97 0.03
## yrs.since.phd 56 55 0.30 -0.81 0.65
## yrs.service 60 60 0.65 -0.34 0.65
## sex* 2 1 -2.69 5.25 0.01
## salary 231545 173745 0.71 0.18 1520.16Jak możemy zauważyć wszysykie grupy są równoliczne.
Poniżej przedstawione są wykresy skrzynkowe obrazujące wynagrodzenie kobiet oraz mężczyzn.
Porządkujemy dane, aby łatwiej było nam zauważyć istotne różnice.
indf1 <- Salaries
v <- c(1:nrow(indf1))
indf1 <- cbind(v,indf1 )
small <- indf1 %>%
group_by(rank, discipline, sex) %>%
get_summary_stats(salary, type = "mean")
small## # A tibble: 12 x 6
## rank discipline sex variable n mean
## <fct> <fct> <fct> <chr> <dbl> <dbl>
## 1 AsstProf A Female salary 6 72933.
## 2 AsstProf A Male salary 18 74270.
## 3 AsstProf B Female salary 5 84190.
## 4 AsstProf B Male salary 38 84647.
## 5 AssocProf A Female salary 4 72128.
## 6 AssocProf A Male salary 22 85049.
## 7 AssocProf B Female salary 6 99436.
## 8 AssocProf B Male salary 32 101622.
## 9 Prof A Female salary 8 109632.
## 10 Prof A Male salary 123 120619.
## 11 Prof B Female salary 10 131836.
## 12 Prof B Male salary 125 133518.Poniżej przedstawione są wykresy skrzynkowe obrazujące wynagrodzenie kobiet oraz mężczyzn.
pudelkowy<-ggboxplot(Salaries, x="rank",
y="salary",color="discipline",palette="jco",facet.by="sex")
pudelkowy
ggplot(data = indf1, aes(y = salary, x = discipline, color = sex)) + geom_boxplot() +
facet_wrap(~rank) + xlab("") + ylab("Salary") 
### Założenia użycia Anovy Żeby mieć możliwość użycia testu Anovy dane muszą spełnić następujące założenia.
Badanie normlaności
W tym kroku sprawdzamy normalność rozkładu zmiennych
indf1 %>%
group_by(rank, discipline, sex) %>%
shapiro_test(salary)## # A tibble: 12 x 6
## rank discipline sex variable statistic p
## <fct> <fct> <fct> <chr> <dbl> <dbl>
## 1 AsstProf A Female salary 0.870 0.226
## 2 AsstProf A Male salary 0.941 0.300
## 3 AsstProf B Female salary 0.889 0.354
## 4 AsstProf B Male salary 0.941 0.0458
## 5 AssocProf A Female salary 0.863 0.269
## 6 AssocProf A Male salary 0.878 0.0113
## 7 AssocProf B Female salary 0.635 0.00117
## 8 AssocProf B Male salary 0.967 0.416
## 9 Prof A Female salary 0.934 0.549
## 10 Prof A Male salary 0.952 0.000259
## 11 Prof B Female salary 0.974 0.923
## 12 Prof B Male salary 0.978 0.0435Z powyższego testu wynika, że dla większości grup nie ma podstaw do odrzucenia hipotezy zerowej o normalności rozkładu. Ponieważ Anova jest mocnym testem, możemy przyjąć normalność rozkładu zmiennych.
Zobaczmy, jak wyglądają wykresy normalności wg podgrup:
ggqqplot(indf1, "salary") +
facet_grid(sex~discipline~rank)
Wykresy podgrup nie pokazują nam idealnej normalności rozkładu zmiennych.
Jednorodność wariancji
Założenie jednorodności wariancji
indf1 %>%
levene_test(salary ~ rank*sex*discipline)## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 11 385 9.05 2.06e-14P-value < 0.05, a więc stwierdzamy występowanie istotnej różnicy wariancji. Zatem wariancje nie są jednorodne.
Obydwa warunki nie zostały spełnione. Bez modyfikacji danych nie jesteśmy wstanie przeprowadzić testu Anovy.
Modyfikacja Danych
Na obecnych danych nie jesteśmy wstanie przeprowadzić testu Anovy. Jedną z popularniejszym możliwości jest zlogarytmowanie zestawu danych Salaries. Dodajemy z zlogarytmizowaną zmienna Salary o nazwie salary_log.
indf1$salary_log <- log(indf1$salary)
indf1 %>%
group_by(rank, discipline, sex) %>%
shapiro_test(salary_log)## # A tibble: 12 x 6
## rank discipline sex variable statistic p
## <fct> <fct> <fct> <chr> <dbl> <dbl>
## 1 AsstProf A Female salary_log 0.852 0.164
## 2 AsstProf A Male salary_log 0.941 0.302
## 3 AsstProf B Female salary_log 0.896 0.387
## 4 AsstProf B Male salary_log 0.931 0.0218
## 5 AssocProf A Female salary_log 0.847 0.216
## 6 AssocProf A Male salary_log 0.904 0.0363
## 7 AssocProf B Female salary_log 0.611 0.000616
## 8 AssocProf B Male salary_log 0.980 0.789
## 9 Prof A Female salary_log 0.936 0.575
## 10 Prof A Male salary_log 0.986 0.235
## 11 Prof B Female salary_log 0.976 0.939
## 12 Prof B Male salary_log 0.991 0.639Możemy zaobserować, że więkoszość danych posiada normalną dystrybucje. Teraz sprawdzamy homogeniczność wairancji.
indf1 %>%
levene_test(salary_log ~ rank*sex*discipline)## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 11 385 7.09 5.41e-11Warunek jednorodności waraincji cały czas nie został spełniony. Dzięki zlogarytmowaniu zmiennej Salaries jesteśmy wstanie przyjąć założenie o normalności jednak warunek jednorodności wariancji blokuje nam możliwość przeprowadzenia Anovy. Przperowadziemy nieparametryczną wersję Anovy, czyli test Kruskala-Wallisa.
indf1grouped <- indf1 %>% unite(col = "group", rank, discipline, sex)
kruskal.test(data =indf1grouped , salary ~ group)##
## Kruskal-Wallis rank sum test
##
## data: salary by group
## Kruskal-Wallis chi-squared = 220.78, df = 11, p-value < 2.2e-16p<0.05 więc odrzucamy H0. Co świadczy o tym, że płeć, dyscyplina oraz tytuł mają znaczący wpływ na pensje osoby.
ANCOVA
Sprawdźmy czy lata od doktoratu (yrs.since.phd), staż pracy (yrs.service) mogą być istotnymi zmiennymi współ-oddziaływującymi na płacę (salary)?
Hipotezy
H0- czas od zdobycia doktoratu oraz staż pracy nie mają znaczący wpływ na pensje osoby. H1- czas od zdobycia doktoratu oraz staż pracy mają znaczący wpływ na pensje osoby.
Zobaczmy jak wyglądają wykresy zależności czasu od doktoratu, a wynagrodzenia.
ggplot(data = indf1, aes(y = yrs.since.phd, x = salary, colour= rank)) + geom_point() +
xlab("") + ylab("Lata od dokroratu") +
ggtitle('Czas od doktoratu, a wynagrodzenie')
Następnie przedstawmy wykres dotyczący drugiej części pytania czyli lata stażu , a pensja przy jednoczesnym uwzględnieniu stanowiska.
ggplot(data = indf1, aes(y = yrs.service, x = salary,colour= rank)) + geom_point() +
xlab("") + ylab("Wysokość stażu") +
ggtitle('Lata stażu, a wynagrodzenie')
Na obu wykresach możemy zauważyć, że Profesorowie są najbardziej rozproszeni.Doktorzy habilitowani (AsstProf) zarabiają najmniej lecz jednocześnie posiadają najniższy staż.
Zależność czasu od doktoratu, a wysokości pensji
Sprawdzamy czy zmienne czasu od doktoratu, a wynagrodzenia są niezależne.
anova(aov(salary ~ group * yrs.since.phd, data = indf1grouped))## Analysis of Variance Table
##
## Response: salary
## Df Sum Sq Mean Sq F value Pr(>F)
## group 11 1.6365e+11 1.4878e+10 28.1805 <2e-16 ***
## yrs.since.phd 1 1.3298e+08 1.3298e+08 0.2519 0.6160
## group:yrs.since.phd 11 2.5915e+09 2.3559e+08 0.4462 0.9343
## Residuals 373 1.9692e+11 5.2794e+08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Zmienne są niezależne, dzięki czemu możemy przejść do wyliczeniu testu ANCOVA.
ANCOVA
Test dla zmiennej czas od doktoratu, a wysokość pensji.
ancova <- aov(salary ~ yrs.since.phd + rank * discipline * sex, data = indf1)
anova(ancova)## Analysis of Variance Table
##
## Response: salary
## Df Sum Sq Mean Sq F value Pr(>F)
## yrs.since.phd 1 6.3852e+10 6.3852e+10 122.8943 < 2.2e-16 ***
## rank 2 7.9596e+10 3.9798e+10 76.5980 < 2.2e-16 ***
## discipline 1 1.8380e+10 1.8380e+10 35.3756 6.12e-09 ***
## sex 1 6.4694e+08 6.4694e+08 1.2452 0.2652
## rank:discipline 2 5.4634e+08 2.7317e+08 0.5258 0.5915
## rank:sex 2 1.6548e+08 8.2739e+07 0.1592 0.8528
## discipline:sex 1 4.6068e+08 4.6068e+08 0.8867 0.3470
## rank:discipline:sex 2 1.4019e+08 7.0097e+07 0.1349 0.8738
## Residuals 384 1.9951e+11 5.1957e+08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Czas od doktoratu ma duży wpływ na wysokość pensji. Aby móc w pełni odrzucić hipoteze 0 należy zbadać wpływ czsu stażu.
Zależność czasu stażu, a wysokości pensji
Sprawdzamy czy zmienne czasu stażu, a wynagrodzenia są niezależne.
anova(aov(salary ~ group * yrs.service, data = indf1grouped))## Analysis of Variance Table
##
## Response: salary
## Df Sum Sq Mean Sq F value Pr(>F)
## group 11 1.6365e+11 1.4878e+10 28.7513 <2e-16 ***
## yrs.service 1 3.4976e+08 3.4976e+08 0.6759 0.4115
## group:yrs.service 11 6.2842e+09 5.7129e+08 1.1040 0.3566
## Residuals 373 1.9301e+11 5.1746e+08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Zmienne są niezależne, dzięki czemu możemy przejść do wyliczeniu testu ANCOVA.
ANCOVA
Test dla zmiennej czas stażu, a wysokość pensji.
ancova <- aov(salary ~ yrs.service + rank * discipline * sex, data = indf1)
anova(ancova)## Analysis of Variance Table
##
## Response: salary
## Df Sum Sq Mean Sq F value Pr(>F)
## yrs.service 1 4.0709e+10 4.0709e+10 78.4376 < 2.2e-16 ***
## rank 2 1.0358e+11 5.1789e+10 99.7848 < 2.2e-16 ***
## discipline 1 1.7617e+10 1.7617e+10 33.9443 1.201e-08 ***
## sex 1 7.7669e+08 7.7669e+08 1.4965 0.2220
## rank:discipline 2 5.0702e+08 2.5351e+08 0.4885 0.6140
## rank:sex 2 2.0277e+08 1.0139e+08 0.1953 0.8226
## discipline:sex 1 4.7892e+08 4.7892e+08 0.9228 0.3374
## rank:discipline:sex 2 1.3479e+08 6.7393e+07 0.1299 0.8783
## Residuals 384 1.9930e+11 5.1900e+08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Czas stażu ma duży wpływ na wysokość pensji. Odżucamy hipoteze 0.
Test KRUSKAL-WALLIS
Zbadajmy czy istnieje istotna różnica w latach od doktoratu (yrs.since.phd) i stażu pracy (yrs.service) różnej rangi profesorów? # Hipotezy H0-nie istnieje różnica w latach od doktoratu i stażu pracy dla różnych rangą profesorów. H1-istnieje różnica w latach od doktoratu i stażu pracy dla różnych rangą profesorów.
ggplot(data = indf1, aes(x = yrs.since.phd, y = yrs.service)) + geom_point() +
facet_grid(~rank) + xlab('czas od doktoratu') + ylab('Czas stażu') +
ggtitle('Czas od doktoratu i stażu pracy różnej rangi profesorów')
kruskal.test(data = indf1, yrs.since.phd - yrs.service ~ rank)##
## Kruskal-Wallis rank sum test
##
## data: yrs.since.phd - yrs.service by rank
## Kruskal-Wallis chi-squared = 7.1521, df = 2, p-value = 0.02799p-value<0,05 Odżucamy hipoteze 0 W związku z tym istnieje różnica w różnicach między latami, które upłynęły od doktoratu, a stażem pracy dla różnych rangą profesorów.