Understanding the data

At the beginning we will add an identifier for each observation - it will be useful later, when dealing with outliers.

sal <- Salaries

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

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

We create a summary table using get_summary_stats to observe mean values of observations. We can see that the increasing trend in salaries is visible when looking at rank of teachers - Assistant Professors earn less than Professors. There is also a difference in discipline - in each case professors from group B earn more than those in group A. Another trend may be seen in gender of professors - men tend to earn more than women. This will be three factors we will consider when performing 3-way Anova: rank, discipline, sex.
Visualization of those variables is visible on a plot below:

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

Assumptions

To be able to perform 3-way Anova, we must check following assumptions: - not significant outliers
- samples come from populations that are normally distributed - normality check
- homogeneity of variances in groups - check of homoscedasticity

1. Outliers

We will try to identify outliers in our data.

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),]

We chose to delete all outliers to make the data easier for later calculations.

2. Normality

To check normality, we will perform shapiro test and visualize the QQ-plots.

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

ggqqplot(residuals(lm(salary ~ rank*discipline*sex, data = sal1))) + theme_minimal() + 
  ggtitle("Residuals for checking normality assumption")

In Shapiro-Wilk test we observe two groups where normality assumption is not met. The QQ plot also does not allow us to see perfectly normal distribution.

3. Homogeneity of variances

To check homogeneity of variance we use Levene’s Test.

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

As we can observe, we cannot accept two main assumptions - normality and homogeneity of variance. P-value from Levene’s Test is basically 0. Hence, we will try making some transformations to the data.

Data transformation

Firstly, we will try to change the data by taking log of salaries. We want to try to find a solution when a 3-way ANOVA can be performed.

Transformed data will again go through assumptions checking.

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

ggqqplot(residuals(lm(salary ~ rank*discipline*sex, data = sal22))) + theme_minimal() + 
  ggtitle("Residuals for checking normality assumption")

a <- leveneTest(data = sal22, salary ~ rank * discipline * sex)

By taking log of salaries we can observe that most of connections rank * discipline * sex turns out to have normal distribution. The only problem occurs with one subgroup. However, we did not manage to increase p-value for homogeneity of variance to be able to accept this assumption.

Hence, we will perform the calculations in two ways. First, we will use log(salaries) to perform a 3-way ANOVA as was asked, even though assumptions are not met. We will assume small sample with not normal distribution does not reflects real population parameters. In second part, we will use non parametric test.

3-way Anova

First, performing 3-way Anova with interactions. We will begin by stating null hypothesis:

• H_0A : Factor A = rank has no effect on the salary of employees in the university.

• H_0B : Factor B = discipline the employees work in has no effect on the salary.

• H_0C : Factor C = sex, gender of employees has no effect on the salary amount.

• H_0AxB : There is no interaction between rank and discipline that has an affect on the salary they get.

• H_0AxC : There is no interaction between rank and sex of employees that has an affect on the salary they get.

• H_0BxC : There is no interaction between discipline and sex of employees performed on participants that has an affect on the salary they get.

• H_0AxBxC : There is no interaction between rank, discipline and sex that has an affect on the salary the employees get.

Now we can create an Anova model and check the null hypothesis.

model1 <- aov(data = sal22, salary ~ rank * discipline * sex)
a <- anova(model1)

From performing 3-way Anova with interactions we can observe these results:
- Rank and discipline have very significant influence on salary value. Hence, we can conclude that there is a significant difference between salaries of Professors, Assistant Professors and Associate Professors as well as between all employees in department A and department B.

• Sex does not turn out to have significant impact on salary. Hence, the null hypothesis cannot be rejected and we do not assume any difference between males and females employees.

• There is no significant interaction between any pair of variables: rank-discipline, rank-sex, discipline-sex, nor rank-discipline-sex. In each of these cases p-value is higher than our significance level. Hence, null hypothesis about existing significant interactions between variables are all not rejected.

model2 <- aov(data = sal22, salary ~ rank + discipline + sex)
a <- anova(model2)

If we wanted to perform 3-way ANOVA without interactions, we will get only three first values from previous calculations. We only check statements whether rank, discipline or sex has any significant impact on the salary of employees - but without mixing the factors. We can again observe that significant values are reported when taking about rank and discipline, however not with gender.

Post-hoc test

As we managed to reject null hypothesis, we will now perform a post-hoc test to check which pairs of possible ranks, disciplines and genders had significant influence on reject H₀.

We should check:

• Simple two-way interaction: run two-way interaction at each level of third variable,

• Simple simple main effect: run one-way model at each level of second variable,

• Simple simple pairwise comparisons: run pairwise or other post-hoc comparisons if necessary.

As we had non significant 3-way interaction nor any significant 2-way interaction, we will only focus on finding main effects.

sal2grouped <- sal2 %>% unite(col = "group", rank, discipline, sex)
pair <- pairwise_t_test(data = sal2grouped, salary ~ group, 
                        p.adjust.method = "bonferroni")

From performing pair-wise comparison between grouped three variables we can observe significant interaction between most of the pairs. Taking from first ten for example between AssocProf_A_Female and AssocProf_B_Female -> what supports the fact that null hypothesis stating that discipline has significant impact on salary score. As there are as many as 66 combinations (12 * 11 / 2) we will only show first 10 for observing the p-value.

Non Parametric Anova

Second possibility we can use is to perform a non parametric version of Anova test. As there exist only an option to do 1-way non parametric Anova - Kruskal-Wallis Anova we will use grouped variable corresponding to rank, discipline and sex at once.

sal1grouped <- sal1 %>% unite(col = "group", rank, discipline, sex)
kruskal.test(data = sal1grouped, salary ~ group)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  salary by group
## Kruskal-Wallis chi-squared = 210.93, df = 11, p-value < 2.2e-16

Using this test, we managed to again reject null hypothesis. However, the only think we know thanks to that is that the rank, discipline a person works in and his/her gender have significant impact on their salary value.

Plots

We would like to visualize possible connections between covariates, salary and other variables values.

ggplot(data = sal22, aes(y = salary, x = yrs.since.phd, color = sex)) + geom_point() + 
  facet_wrap(~discipline~rank) + xlab("") + ylab("Salary") +
  ggtitle('Years since phd vs salary for different groups of Professors')

ggplot(data = sal22, aes(y = salary, x = yrs.service, color = sex)) + geom_point() + 
  facet_wrap(~discipline~rank) + xlab("") + ylab("Salary") +
  ggtitle('Years of service vs for different groups of Professors')

On both plots we can see that the Professors are mostly spread - in both years of service / since phd and salary. Assistant Professors earn less and also have the lowest seniority.

Years since phd

In this test we should also take a look at Ancova assumptions. Most of them were checked in previous exercise, however we will also want to take a look at few others. We are using the log(salary) to be closer to our normality and homogeneity assumptions.

Here we want to check whether the normality assumption of covariate (yrs.since.phd) is met.

ggscatter(sal22, x = "yrs.since.phd", y = "salary", 
          facet.by = c("rank", "discipline")) + stat_smooth(method = "loess")

On this plot we can also observe how years since phd and log(salary) looks like in terms of different ranks and disciplines of employees.

We also must take a look if our covariate variable is independent. To do that, we will perform a Anova test.

a <- anova(aov(salary ~ group * yrs.since.phd, data = sal2grouped))

Here we can observe that interaction between grouped variable and our covariate is not significant (0.85). Hence, this assumption is fulfiled.

Ancova Test

ancova <- aov(salary ~ yrs.since.phd + rank * discipline * sex, data = sal22)
a <- anova(ancova)

Having performed this Ancova test, we can observe that years since phd does have strong impact on salary value. We could have observed that on also previous plots, hence we reject the null hypothesis.

Years of service

Now we will take a look at other suspicious variable - years of service.

First, ploting it agains salary.

ggscatter(sal22, x = "yrs.service", y = "salary", 
          facet.by = c("rank", "discipline")) + stat_smooth(method = "loess")

a <- anova(aov(salary ~ group * yrs.service, data = sal2grouped))

We can observe that our variable is independent on group. Hence, we can fulfill that assumption.

ancova2 <- aov(salary ~ yrs.service + rank * discipline * sex, data = sal22)
a <- anova(ancova2)

This Ancova test also shows significant impact of years of service on salary value.

Hence, both our suspictious variables turned out to be good covariates as they should be taken into account when talking about salary of professors.

Plots & Data

In the last task, we were asked to find whether there exists a difference between years since PhD and years of service of different rank professors. We will start with ploting the data.

ggplot(data = sal1, aes(x = yrs.since.phd, y = yrs.service)) + geom_point() + 
  facet_grid(~rank) + xlab('Years since PhD') + ylab('Years of service') + 
  ggtitle('Years since PhD vs years of service for different Professors')

To check last exercise, we will create a new variable called diff. This will keep the value of difference between years since PhD - years of service. Positive value will correspond to a professor that started his/her learning adventure after completing PhD. Negative value will mean that he/she was teaching before that.

saldiff <- sal1 %>% mutate(diff = yrs.since.phd - yrs.service)

a <- saldiff %>%
  group_by(rank) %>%
  get_summary_stats(diff, type = "mean")

Summary of data can show us that there occurs to be some difference between ranks. Professors tend to have bigger mean difference between years since PhD and years of service. We must take into account however, that these people are usually older and tend to have worked longer. However, we cannot check that using our dataset.

ggplot(data = saldiff, aes(y = diff, x = rank, color = rank)) + geom_boxplot() + 
  theme(legend.position = 'none') + xlab('') + ylab('yrs since PhD - yrs of service') + 
  ggtitle('Difference between years since PhD and years of service for different Professors')

Assumptions

In this situation we would like to perform 1-way ANOVA. However, first we must take a look at assumptions.

• samples are independent
  
• samples come from populations that are normally distributed - normality check
  
• homogeneity of variances in groups - check of homoscedasticity
• not significant outliers

1. Outliers check

outliers_data <- saldiff %>%
  group_by(sex, rank, discipline) %>%
  identify_outliers(salary)

x <- outliers_data %>% filter(is.extreme == 'TRUE')
x <- x[, 4]
x <- as.numeric(unlist(x))

saldiff2 <- saldiff[-which(saldiff$v %in% x),]

As we found some significant outliers, we will delete them from our dataset.

2. Normality check

a <- saldiff2 %>%
  group_by(rank) %>%
  shapiro_test(diff)

qqnorm(saldiff2$diff, main="QQ plot", pch=19)
qqline(saldiff2$diff)

As we can observe, we are not dealing with normal data. We will check homoscedasticity, however even then we will use non - parametric test.

3. Homoscedasticity

a <- leveneTest(data = saldiff2, diff ~ rank)

Our data do not contain homogeneous variance. Hence, we must perform non parametric test - Kruskal-Wallis.

Kruskal-Wallis Test

kruskal.test(data = saldiff2, diff ~ rank)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  diff by rank
## Kruskal-Wallis chi-squared = 8.4193, df = 2, p-value = 0.01485

We obtained a very nice value of p 0.01 -> meaning with our significance level 0.05 we can reject null hypothesis. Hence, there is a difference in differences between years since PhD and years of service for different rank professors.