Plots

First, we want to review the data to view possible patterns. We have three independent variables and one response = pain score value.

x <- sample(1:nrow(headache), 10)
x <- sort(x)

kbl(headache[x, ], caption = "Random 10 rows from headache dataset", col.names = c("ID", "Gender", "Risk", "Treatment", "Pain score")) %>%
  kable_minimal() %>%
  column_spec(5, background = spec_color(headache$pain_score[1:10], end = 0.8))

• id => individual identifier of the participant, quantitative variable.
• gender => qualitative nominal variable ‘male’ or ‘female’
• risk => qualitative nominal variable ‘low’ or ‘high’ (describes how likely a person is to have migraine)
• treatment => qualitative nominal variable ‘X’, ‘Y’, ‘Z’ describes three drugs given to the participants
• pain_score => quantitative variable describing the pain score given by the participant.

Now we want to see which of three independent variables affect pain score mostly.

plot.design(pain_score ~ gender * risk * treatment, data = headache, col = "brown", ylab = "Mean of pain score", main = "Impact of factors on pain score") 

Using plot.desing function we can observe that the biggest influence on pain_score value has risk level = whether the person is likely or less likely to have migraine. However, the impact of the treatment is much smaller - it is actually similar to the gender factor. We will look into that later with ANOVA analysis.

Another way to look into these three factors is by looking at the connection of them - for example whether gender and risk level together impact the pain score given by the participants.

interaction.plot(headache$gender, headache$risk, headache$pain_score, col = c("blue", "purple"), 
                 lty = 1, lwd = 2, ylim = c(70, 85), trace.label = "Risk",
                 xlab = "Gender", ylab = "Mean of pain score")

interaction.plot(headache$gender, headache$treatment, headache$pain_score, col = c("blue", "pink", "purple"), 
                 lty = 1, lwd = 2, ylim = c(70, 85), trace.label = "Treatment",
                 xlab = "Gender", ylab = "Mean of pain score")

interaction.plot(headache$risk, headache$treatment, headache$pain_score, col = c("blue", "pink", "purple"), 
                 lty = 1, lwd = 2, ylim = c(70, 85), trace.label = "Treatment",
                 xlab = "Risk level", ylab = "Mean of pain score")

Here in the three plots we can observe that:
1. Gender and risk level of having migraine. are definitely correlated. The mean pain score for men was higher that for women in both cases, when participants had high and also low risk of having migraine.
2. Gender and treatment performed showed significant difference in case of only one treatment - drug X when value of mean pain score is much higher for males than foe females. When using drug Z and Y, the difference between males and females pains scores is not as big.
3. Risk level of having migraine and treatment performed also showed some collocation. All three drugs showed smaller pain score for people not as likely to have migraine.

By connecting all data on one boxplot we can observe all at one place. Hence, the distinction between low and high risk level definitely affects pain score. The difference between gender is not as visible - however, the use of drug X with high risk shows a significant difference between genders (as we saw previously on interaction plots).

ggplot(aes(x = pain_score, y = risk, color = gender), data = headache) + geom_boxplot() + facet_grid(~treatment) + xlab("Pain score") + ylab("Risk level") + ggtitle("Pain score distribution grouped by gender, risk and treatment")

Here we will add new variable that sums up the independent variables in following pattern:
male + low risk + treatment X -> mlX

Then we are able to create a boxplot of all distributions separately.

h1 <- headache %>% mutate('Total' = ifelse(gender == 'male' & risk == 'low' & treatment == 'X', 'mlX', 
                                    ifelse(gender == 'male' & risk == 'low' & treatment == 'Y', 'mlY', 
                                    ifelse(gender == 'male' & risk == 'low' & treatment == 'Z', 'mlZ', 
                                    ifelse(gender == 'male' & risk == 'high' & treatment == 'X', 'mhX', 
                                    ifelse(gender == 'male' & risk == 'high' & treatment == 'Y', 'mhY', 
                                    ifelse(gender == 'male' & risk == 'high' & treatment == 'Z', 'mhZ', 
                                    ifelse(gender == 'female' & risk == 'low' & treatment == 'X', 'flX', 
                                    ifelse(gender == 'female' & risk == 'low' & treatment == 'Y', 'flY', 
                                    ifelse(gender == 'female' & risk == 'low' & treatment == 'Z', 'flZ', 
                                    ifelse(gender == 'female' & risk == 'high' & treatment == 'X', 'fhX', 
                                    ifelse(gender == 'female' & risk == 'high' & treatment == 'Y', 'fhY', 
                                    ifelse(gender == 'female' & risk == 'high' & treatment == 'Z', 'fhZ', 0)))))))))))))

ggplot(h1, aes(x = pain_score, y = Total, fill = Total)) + geom_boxplot() + scale_fill_brewer(name = "", palette="Set3") + xlab("Pain score") + ylab("") + ggtitle("Pain score distribution grouped by gender, risk and treatment") + theme(legend.position = 'none')

Analysis

x <- summary(headache[c('gender', 'risk', 'treatment', 'pain_score')])
kbl(x)

By looking at the summary of the data we can view how many observations fall into each category. The data are equally distributed, meaning there are the same numbers of observations grouped by gender, risk and treatment in each 2 (or 3 in treatment) subgroups.
The mean, median, percentiles, minimum and maximum values of the pain score are also given.

Now too look at the descriptive statistics individually grouped into categories:

kbl(headache %>%
  group_by(gender, risk, treatment) %>%
  get_summary_stats(pain_score, type = "full"))

Outliers

We will try by checking if there are any significant outliers it the data.

outliers_data <- headache %>%
  group_by(gender, risk, treatment) %>%
  identify_outliers(pain_score)

kbl(outliers_data)

We can observe that there are 4 values that strongly influence the data. As we could have seen on the graph above this points are easily visible.

However, before we delete those values we will check normality assumption. If it is fulfilled, there is no point in removing observations. We must take into consideration that we are dealing with really small dataset and removing observations (as done with headache1) can lead us to completely getting rid of all data in a subgroup.

Now trying a deletion of these outliers:

x <- outliers_data[, 4]
x <- unlist(x) 

headache1 <- headache %>% filter(id != x)

ggplot(aes(x = pain_score, y = risk, color = gender), data = headache1) + geom_boxplot() + facet_grid(~treatment) 

As we see, this still leaves us with another outliers.

Normality

First, lets look at the normality assumption. We will use Shapiro test. The null hypothesis in the test is H₀ : Sample data come from normal population.

x <- mshapiro_test(headache$pain_score)
kbl(x)

x <- headache %>%
  group_by(gender) %>%
  shapiro_test(pain_score)
kbl(x)

x <- headache %>%
  group_by(risk) %>%
  shapiro_test(pain_score)
kbl(x)

x <- headache %>%
  group_by(treatment) %>%
  shapiro_test(pain_score)
kbl(x)

We can see that p value is higher than our significance level. Hence, we do not reject the null hypothesis and consider our data normally distributed. Below we create a plot to view residuals and how they follow normal distribution.

ggqqplot(residuals(lm(pain_score ~ gender*risk*treatment, data = headache))) + theme_minimal() + ggtitle("Residuals for checking normality assumption")

As our outliers do not disturb with normal distribution of the dataset, we will continue using headache without deleted outliers.

Homogeneity of variance

To check the homogeneity of variances we can use two tests: Bartlett’s Test and Levene’s Test. Levene’s test is more robust and not as sensitive to outliers. Hence, it is better to use in this case.

x <- leveneTest(pain_score ~ risk * treatment * gender, data = headache)
kbl(x)

From the Levene’s Test we can see that p value (almost 1) is much higher than our alpha level 0.05. Hence, the assumption of homogeneity of variance is fulfilled.

Two-way interactions

headacheX <- headache %>% filter(treatment == 'X')
headacheZ <- headache %>% filter(treatment == 'Z')
headacheY <- headache %>% filter(treatment == 'Y')

modelX <- aov(data = headacheX, pain_score ~ gender * risk)
kbl(anova(modelX))

modelZ <- aov(data = headacheZ, pain_score ~ gender * risk)
kbl(anova(modelZ))

modelY <- aov(data = headacheY, pain_score ~ gender * risk)
kbl(anova(modelY))

We can now observe that significant two-way interaction was observed only for modelX being the use of drug X.

Another way to separate the three-way interaction is by dividing it into gender.

headachemen <- headache %>% filter(gender == 'male')
headachewomen <- headache %>% filter(gender == 'female')

modelmen <- aov(data = headachemen, pain_score ~ treatment * risk)
kbl(anova(modelmen))

modelwomen <- aov(data = headachewomen, pain_score ~ treatment * risk)
kbl(anova(modelwomen))

Here the only significant two-way interaction is with male participants.

Last division is by risk level.

headachelow <- headache %>% filter(risk == 'low')
headachehigh <- headache %>% filter(risk == 'high')

modellow <- aov(data = headachelow, pain_score ~ gender * treatment)
kbl(anova(modellow))

modelhigh <- aov(data = headachehigh, pain_score ~ gender * treatment)
kbl(anova(modelhigh))

Hence, the only significant two-way interaction is when patients have high risk of migraine.

Main effects

Now, to find simple main effect we will perform those significant two-way interactions from the previous step checking the second variable.

1. Significant two-way interaction gender:risk (treatment X)

headacheXlow <- headacheX %>% filter(risk == 'low')
headacheXhigh <- headacheX %>% filter(risk == 'high')

modelXlow <- aov(data = headacheXlow, pain_score ~ gender)
kbl(anova(modelXlow))

modelXhigh <- aov(data = headacheXhigh, pain_score ~ gender)
kbl(anova(modelXhigh))

Significant main effect: gender for treatment X + high risk level.

headacheXmen <- headacheX %>% filter(gender == 'male')
headacheXfem <- headacheX %>% filter(gender == 'female')

modelXmen <- aov(data = headacheXmen, pain_score ~ risk)
kbl(anova(modelXmen))

modelXfem <- aov(data = headacheXfem, pain_score ~ risk)
kbl(anova(modelXfem))

Significant main effect: risk for treatment X + gender male.

2. Significant two-way interaction risk:treatment (male)

headachemenX <- headachemen %>% filter(treatment == 'X')
headachemenY <- headachemen %>% filter(treatment == 'Y')
headachemenZ <- headachemen %>% filter(treatment == 'Z')

modelmenX <- aov(data = headachemenX, pain_score ~ risk)
kbl(anova(modelmenX))

modelmenY <- aov(data = headachemenY, pain_score ~ risk)
kbl(anova(modelmenY))

modelmenZ <- aov(data = headachemenZ, pain_score ~ risk)
kbl(anova(modelmenZ))

Significant main effect: risk for treatment X + gender male.
Significant main effect: risk for treatment Y + gender male.

headachemenlow <- headachemen %>% filter(risk == 'low')
headachemenhigh <- headachemen %>% filter(risk == 'high')

modelmenlow <- aov(data = headachemenlow, pain_score ~ treatment)
kbl(anova(modelmenlow))

modelmenhigh <- aov(data = headachemenhigh, pain_score ~ treatment)
kbl(anova(modelmenhigh))

Significant main effect: treatment for high risk + gender male.

3. Significant two-way interaction gender:treatment (high risk)

headachehighmen <- headachehigh %>% filter(gender == 'male')
headachehighfem <- headachehigh %>% filter(gender == 'female')

modelhighmen <- aov(data = headachehighmen, pain_score ~ treatment)
kbl(anova(modelhighmen))

modelhighfem <- aov(data = headachehighfem, pain_score ~ treatment)
kbl(anova(modelhighfem))

Significant main effect: treatment for high risk + gender male.

headachehighX <- headachehigh %>% filter(treatment == 'X')
headachehighY <- headachehigh %>% filter(treatment == 'Y')
headachehighZ <- headachehigh %>% filter(treatment == 'Z')

modelhighX <- aov(data = headachehighX, pain_score ~ gender)
kbl(anova(modelhighX))

modelhighY <- aov(data = headachehighY, pain_score ~ gender)
kbl(anova(modelhighY))

modelhighZ <- aov(data = headachehighZ, pain_score ~ gender)
kbl(anova(modelhighZ))

Significant main effect: gender for high risk + treatment X.

SUMMARY: From all these ANOVA tests we can conclude that our main effects are:

• gender for high risk + treatment X
  
• treatment for high risk + gender male
  
• risk for treatment X + gender male
  
• risk for treatment Y + gender male

Pairwise comparisons

Lastly, we can perform TukeyHSD post-hoc test to view the summary in a different form.

pairwise_t_test(data = h1, pain_score ~ Total, p.adjust.method = "bonferroni")

## # A tibble: 66 x 9
##    .y.        group1 group2    n1    n2          p p.signif   p.adj p.adj.signif
##  * <chr>      <chr>  <chr>  <int> <int>      <dbl> <chr>      <dbl> <chr>       
##  1 pain_score fhX    fhY        6     6 0.367      ns       1   e+0 ns          
##  2 pain_score fhX    fhZ        6     6 0.396      ns       1   e+0 ns          
##  3 pain_score fhY    fhZ        6     6 0.956      ns       1   e+0 ns          
##  4 pain_score fhX    flX        6     6 0.0686     ns       1   e+0 ns          
##  5 pain_score fhY    flX        6     6 0.00757    **       5   e-1 ns          
##  6 pain_score fhZ    flX        6     6 0.00879    **       5.8 e-1 ns          
##  7 pain_score fhX    flY        6     6 0.000112   ***      7.38e-3 **          
##  8 pain_score fhY    flY        6     6 0.00000448 ****     2.95e-4 ***         
##  9 pain_score fhZ    flY        6     6 0.00000548 ****     3.62e-4 ***         
## 10 pain_score flX    flY        6     6 0.0261     *        1   e+0 ns          
## # ... with 56 more rows

Looking at the pairwise comparison we can observe that there are some significant correlations marked with * stars. Rest of them, marked with ‘ns’ means not significant values.