library(rstatix) library(dplyr) library(car)

data(“Salaries”, package = “carData”) sal <- Salaries

v <- c(1:nrow(sal)) sal <- cbind(v, sal)

a <- sal %>% group_by(rank, discipline, sex) %>% get_summary_stats(salary, type = “mean”)

#plots ggplot(data = sal, aes(y = salary, x = discipline, color = sex)) + geom_boxplot() + facet_wrap(~rank) + xlab(““) + ylab(”Salary”)

#Assumptions outliers_data <- sal %>% group_by(sex, rank, discipline) %>% identify_outliers(salary)

x <- outliers_data %>% filter(is.outlier == ‘TRUE’) x <- x[, 4] x <- as.numeric(unlist(x)) sal1 <- sal[-which(sal$v %in% x),]

#normality

a <- sal1 %>% group_by(rank, discipline, sex) %>% shapiro_test(salary)

library(ggqqplot)

ggqqplot(Salaries, “salary”, ggtheme = theme_light()) + facet_grid(sex + rank ~ discipline) + labs(title = “Normality of salary in relation to other factors”)

#Homogeneity of variances a <- leveneTest(salary ~ sex * rank * discipline, data = sal1)

Ancova

ggplot(data = sal1, aes(y = salary, x = yrs.since.phd, color = sex)) + geom_point() + facet_wrap(_disciplinerank) + xlab(““) + ylab(”Salary”) + ggtitle(‘Years since phd vs salary for different groups of Professors’)

ggplot(data = sal1, aes(y = salary, x = yrs.service, color = sex)) + geom_point() + facet_wrap(_disciplinerank) + xlab(““) + ylab(”Salary”) + ggtitle(‘Years of service vs for different groups of Professors’)

#mancova manova_model <- lm(cbind(yrs.since.phd, yrs.service) ~ salary, data = sal1) library(MASS) manova_result <- manova(manova_model) summary(manova_result, test = “Pillai”)