anova.Rmd

---
title: "ANOVA"
author: "May"
date: "2/20/2021"
output: 
  html_document:
    self_contained: false
    toc: true
    toc_float: true
    toc_depth: 4
    keep_tex: yes
    theme: united
    highlight: tango
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, include=F}
colorize <- function(x, color) {
  if (knitr::is_latex_output()) {
    sprintf("\\textcolor{%s}{%s}", color, x)
  } else if (knitr::is_html_output()) {
    sprintf("<span style='color: %s;'>%s</span>", color, 
      x)
  } else x
}
```

```{r, warning=F, message=F}
library(tidyverse)
library(ggplot2)
```


# Diet Data Set

This data set contains information on 76 people using one of three diets.

**Source** 
[link](https://www.sheffield.ac.uk/mash/statistics/datasets) 


**Format** 

- `gender`:   1 = male, 0 = female
- `Age`: Age (years)
- `Height`: Height (cm)
- `preweight`: Weight before the diet (kg)
-  `Diet`:  Diet Type
- `weight6weeks` - Weight after 6 weeks (kg)


```{r}
diet = read.csv("datasets/Diet_R.csv")
diet$weight_loss = diet$pre.weight - diet$weight6weeks
diet$Diet = factor(diet$Diet)
```

```{r}
str(diet)
```

```{r Histogram of Weight Lost}
diet %>% ggplot(aes(x=weight_loss)) +
  geom_histogram(bins = 25, color = 'white', fill = 'darkblue') +
  labs(x= "Weight Lost", title = "Histogram of Weight Lost")
```

```{r Histogram of Weight Lost by Diet , warning=F}
diet %>% ggplot(aes(x=weight_loss, fill = Diet)) +
  geom_histogram(alpha= .5, color = 'white', position =  'identity') +
  labs(x= "Weight Lost", title = "Histogram of Weight Lost by Diet")
```

```{r}
diet %>% group_by(Diet) %>% summarise(mean= mean(weight_loss),
                                      sd = sd(weight_loss))
```



## One-way ANOVA


```{r Box plot weight loss}
ggplot(data = diet, aes(x=Diet, y= weight_loss, color=Diet)) +
  geom_boxplot(outlier.colour="red", outlier.shape=8,outlier.size=3) +
  labs(x="Diet Type", y= "Weight Loss") +
  theme(legend.position = "")
```

### Analysis of Variance 

```{r}
diet.aov = aov(weight_loss~Diet,data=diet)
summary(diet.aov)
```

### Checking Assumptions

#### Levin's test
```{r, warning=F}
library(lawstat)
levene.test(diet$weight_loss, diet$Diet, location = "mean")



# manually

# adding column of the group mean of weight loss 
diet = diet %>% group_by(Diet) %>% mutate(group_mean = mean(weight_loss))
# adding Z column: 
diet$z = abs(diet$weight_loss - diet$group_mean)
  
aov_equal_var = aov(z~Diet,data=diet)
summary(aov_equal_var)
```


#### Diagnostic plots

```{r Checking Assumptions r}
plot(diet.aov)
```

We can also create these plot manually:
```{r Checking Assumptions }
# residudal vs fitted plot:
diet.resid = resid(diet.aov)
diet.predict = predict(diet.aov)
plot(diet.resid~diet.predict)
abline(0,0, col="red")

## qqplot:
qqnorm(diet.resid)
qqline(diet.resid)

```


### Post-hoc analysis 

#### Estimated Expectations
```{r}
tapply(X = diet$weight_loss, INDEX =  diet$Diet, FUN = mean )
```

#### Finding which group is different from all other groups 

`r colorize("????? ", "red")` 




#### Pairwise Comparisons

##### Regular (No multiple correction)
```{r, results='hide' , message=F, warning=F}
library(multcomp)
```


```{r }
contr = rbind(
  "2 - 1" = c(-1,1,0),
  "3 - 1" = c(-1,0,1),
  "3 - 2" = c(0,-1,1)
)

### set up general linear hypothesis
pairwise_res = glht(diet.aov,  linfct = mcp(Diet = contr))
summary(pairwise_res,test = adjusted(type="none"))

### t.test
pairwise.t.test(diet$weight_loss, 
                diet$Diet,
                p.adjust.method = "none")



# calpha = univariate_calpha() - without correction 
pairwise_ci = confint(pairwise_res, calpha = univariate_calpha())
pairwise_ci
```

```{r Pairwise Comparisons}
plot(pairwise_ci, xlab = "Weight loss diff")
```


##### Bonferroni Correction
```{r Bonferroni Correction }

summary(pairwise_res, test = adjusted("bonferroni"))

### t.test
pairwise.t.test(diet$weight_loss, 
                diet$Diet,
                p.adjust.method = "bonferroni")



alpha = 0.05 ; m = 3 
pairwise_ci_bonferroni  =  confint(pairwise_res, level = 1- (alpha/m) , calpha = univariate_calpha())
pairwise_ci_bonferroni
plot(pairwise_ci_bonferroni)
```




##### Tukey's Method
```{r Tukeys Method }
pairwise_res_tukey = TukeyHSD(diet.aov)
pairwise_res_tukey
```

```{r}
plot(pairwise_res_tukey)
```




#### Multiple Comparisons

##### Prespecified Comparison

`r colorize("Any Other Idea???", "red")`    

$$
H_0: \frac{\alpha_1+\alpha_2}{2}-\alpha_3=0 \leftrightarrow  \frac{\mu_1 + \mu_2}{2} - \mu_3 = 0 
$$

```{r}
contr1 = (c(1/2,1/2,-1))
glht_test = glht(diet.aov,  linfct = mcp(Diet = contr1))
summary(glht_test)

# other option:
lvls = data.frame(diet_type = unique(diet$Diet))
X = model.matrix(~ diet_type  , data = lvls)
summary( glht(diet.aov,  linfct = contr1 %*% X ))

```



##### Scheffe Method

`r colorize("Any Other Idea???", "red")`  

```{r, warning=F}
library(DescTools)
ScheffeTest(diet.aov , contrasts = contr1, 
            conf.level = 0.95) 
```


## Two-way ANOVA

**Effect of diet and gender on weight lost** 

Gender:  
1- male, 0 - female   


```{r diet Boxplot of weight lost by gender }
diet$gender = factor(diet$gender)
ggplot(data = diet, aes(x=Diet, y= weight_loss, color=gender)) +
  geom_boxplot(outlier.colour="red", outlier.shape=8,outlier.size=2) +
  labs(x="Diet Type", y= "Weight Loss") +
  scale_color_discrete(name = "Gender", labels = c("Female", "Male"))
```



```{r diet two-way checking assumptions }
diet.aov2 = aov( weight_loss~Diet + gender ,data=diet)
summary(diet.aov2) 

plot(diet.aov2)
```




**F tests:**
```{r}
diet.aov = aov(weight_loss~Diet,data=diet[!is.na(diet$gender),] )
results = anova( diet.aov, diet.aov2)
print(results)



# note it is equivalent to aov() even if the data is unbalanced
results2 = summary(diet.aov2)
results$`Pr(>F)`[2]
results2[[1]]$`Pr(>F)`[2] 

```




## Intercations

```{r}
diet.aov2_intr = aov( weight_loss~Diet * gender ,data=diet)
summary(diet.aov2_intr) 
```

**Interaction Plot**
```{r interatcion plot diet }
interaction.plot(x.factor = diet$Diet,
                 trace.factor = diet$gender, 
                 response = diet$weight_loss,
                 fun = mean, 
                 type = "b", legend = TRUE, 
                 xlab = "Diet", ylab="Weight Loss",
                 trace.label = "Gender",
                 pch=c(1,8), col = c("#FF6347", "#3CB371"))
```

```{r}
plot(diet.aov2_intr)
```


**Levin's test**

```{r, warning=F}
diet = diet[!is.na(diet$gender),]
groups_df = cbind ( expand.grid(Diet = unique(diet$Diet),gender =  unique(diet$gender)), group = factor(1:6))
diet2 = inner_join(diet ,groups_df, by = c("Diet", "gender") )

# option 1 
lawstat::levene.test(y = diet2$weight_loss, group = diet2$group , location = "mean" )

# option 2 
DescTools::LeveneTest(weight_loss ~ Diet * factor(gender)  ,data=diet )

```





# Interaction Plots Examples

**Invented example:**  
Ranking of movies from different genres by women and men.

## No interaction 

```{r rating interaction plot 1, warning=F}
set.seed(22)
group <- gl(n = 2, k = 20, labels = c("Female","Male"))
trt <- rep(rep(c("Adventure","Comedy"), each=10),2)


resp <- c(
  rnorm(n = 20, mean = rep(c(3.60,3.65), each = 10), sd = .01),
  rnorm(n = 20, mean = rep(c(3.62,3.67), each = 10), sd = .01)
)


df1 <- data.frame(group, trt, resp)
interaction.plot(x.factor = df1$trt, 
                 trace.factor = df1$group,
                 response = df1$resp,
                 type = "b", legend = TRUE,
                 xlab = "Genre", ylab="Movie Rating",
                 trace.label = "Gender",
                 pch=c(1,8), col = c("#FF6347", "#3CB371") )




```

## with interaction

```{r rating interaction plot 2 }
set.seed(22)
group <- gl(n = 2, k = 20, labels = c("Female","Male"))
trt <- rep(rep(c("Sci-fi","Romance"), each=10),2)

resp <- c(
  rnorm(n = 20, mean = rep(c(3.60,3.65), each = 10),  sd = .01),
  rnorm(n = 20, mean = rep(c(3.65,3.60), each = 10),  sd = .01)
)

df2 <- data.frame(group, trt, resp)
interaction.plot(x.factor = df2$trt, 
                 trace.factor = df2$group,
                 response = df2$resp,
                 type = "b", legend = TRUE,
                 xlab = "Genre", ylab="Movie Rating",
                 trace.label = "Gender",
                 pch=c(1,8), col = c("#FF6347", "#3CB371") )


```

```{r rating interaction plot 3}
set.seed(22)
trt <- rep(rep(c("Romantic Comedy", "Musicals"), each=10),2)

resp <- c(
  rnorm(n = 20, mean = rep(c(3.25,3.17), each = 10),  sd = .01),
  rnorm(n = 20, mean = rep(c(3.10,3.13), each = 10),  sd = .01)
)

df3 <- data.frame(group, trt, resp)
interaction.plot(x.factor = df3$trt, 
                 trace.factor = df3$group,
                 response = df3$resp,
                 type = "b", legend = TRUE,
                 xlab = "Genre", ylab="Movie Rating",
                 trace.label = "Gender",
                 pch=c(1,8), col = c("#FF6347", "#3CB371") )


```


# Ricci: Firefighter Promotion Exam Scores Dataset


**Details**  

The city of New Haven, Connecticut administered exams (both written and oral) in November and December of 2003 to firefighters hoping to qualify for promotion to either Lieutenant or Captain in the city fire department. A final score consisting of a 60% weight for the written exam and a 40% weight for the oral exam was computed for each person who took the exam. For each person who took the exams, there are measurements on their race (black, white, or Hispanic), which position they were trying for (Lieutenant, Captain), scores on the oral and written exams, and the combined score. These data were used as part of a court case (Ricci v.DeStefano) dealing with racial discrimination

**Source** 

An article on using these data: Miao, W. (2011) "Did the Results of Promotion Exams Have a Disparate Impact on Minorities? Using Statistical Evidence in Ricci v. DeStefano," JSE 19:1 at 
[link](http://jse.amstat.org/v18n3/miao.pdf) 


**Format** 

- `Race`:  Race of firefighter (B=black, H=Hispanic, or W=white)
- `Position`: Promotion desired (Captain or Lieutenant)
- `Oral`: Oral exam score
- `Written`: Written exam score
-  `Combine`:  Combined score (written exam gets 60% weight)

```{r}
library(Stat2Data)

data(Ricci)

str(Ricci)
```



```{r ricci boxplot of score by race}
ggplot(data = Ricci, 
       aes(x=Race, y= Combine, color=Position)) +
  geom_boxplot(outlier.colour="red", outlier.shape=8,outlier.size=3) +
  labs(x="Race", y= "Combined score") 
```






## One-way ANOVA

### Option 1 

```{r}
aov_captain = aov(Combine  ~ Race , data = Ricci %>% filter(Position == "Captain"))
summary(aov_captain)
```

#### Tuke's 95% CI's

```{r  }
ricci_capt_tukey = TukeyHSD(aov_captain)
ricci_capt_tukey
```

```{r}
plot(ricci_capt_tukey)
```



### Option 2

```{r}
aov_lieutenant = aov(Combine  ~ Race , data = Ricci %>% filter(Position == "Lieutenant"))
summary(aov_lieutenant)
```

```{r }
ricci_lie_tukey = TukeyHSD(aov_lieutenant)
ricci_lie_tukey
```

```{r}
plot(ricci_lie_tukey)
```


## Two-way ANOVA

```{r}
summary(aov(Combine  ~ Race*Position , data = Ricci))
```


```{r ricci two way anova interaction plot}
# interaction plot 
interaction.plot(x.factor = Ricci$Race,
                 trace.factor = Ricci$Position, 
                 response = Ricci$Combine,
                 fun = mean, 
                 type = "b", legend = TRUE, 
                 xlab = "Race", ylab="Combine Score",
                 trace.label = "Position",
                 pch=c(1,8), col = c("#FF6347", "#3CB371"))
```


```{r ricci two way anova interaction plot2}
interaction.plot(x.factor = Ricci$Position,
                 trace.factor = Ricci$Race, 
                 response = Ricci$Combine,
                 fun = mean, 
                 type = "b", legend = TRUE, 
                 xlab = "Position", ylab="Combine Score",
                 trace.label = "Race",
                 pch=c(1,8, 3), col = c("#FF6347", "#3CB371","#2980B9" ))
```



# One-factor-at-a-time 

Example:

```{r One-factor-at-a-time}
res_s = c(13, 24, 18, 16 ,14, 27, 30, 14, 14)
df = cbind(expand.grid(Treatment_1 = c("none", "low", "high"),
            Treatment_2 = c("none", "low", "high")), y_exp =res_s)
ggplot(data = df, aes(x=Treatment_1, y = Treatment_2, size = y_exp, color = factor(y_exp))) +
  geom_point() + scale_size(guide = "none") +  guides(color=guide_legend(title="E[Y]")) +
  labs(x= "Treatment 1", y= "Treatment 2") + theme_bw() 


```



# Mixed Effect Model


## Positive Negative Feedback  - 

**The effects of positive and negative verbal feedback on repeated force production**

#### Arranging the data in a long format


```{r}
mvc = read.csv("datasets/pos_neg_feedback.csv")
emg_bi = read.csv("datasets/BICEPS.csv")
emg_tr = read.csv("datasets/TRICEPS.csv")

emg_bi = emg_bi %>% rename_all(janitor::make_clean_names) %>%
  rename_all(~str_remove(., "_")) %>% rename( "gender" = "a")

emg_tr = emg_tr %>% rename_all(janitor::make_clean_names) %>%
  rename_all(~str_remove(., "_"))  %>% rename( "gender" = "a")
  
```

```{r, eval=F}
to_long = function(data, activity) {
  neg_long = data %>% mutate("id"= row_number()) %>%
  select( neg1:neg12, id, gender) %>% 
  pivot_longer(cols = neg1:neg12, values_to = "performance", names_to = "time") %>%
  mutate_at("time", ~(str_extract(., pattern = "[0-9]+"))) %>%
  mutate(feedback = "negative", activity = activity)

  pos_long =  data %>% mutate("id"= row_number()) %>%
  select( pos1:pos12, id, gender) %>% 
  pivot_longer(cols = pos1:pos12, values_to = "performance", names_to = "time") %>%
  mutate_at("time", ~(str_extract(., pattern = "[0-9]+"))) %>%
  mutate(feedback = "positive",  activity = activity)

  con_long =  data %>% mutate("id"= row_number()) %>%
  select( con1:con12, id, gender) %>% 
  pivot_longer(cols = con1:con12, values_to = "performance", names_to = "time") %>%
  mutate_at("time", ~(str_extract(., pattern = "[0-9]+"))) %>%
  mutate(feedback = "no feedback",  activity = activity)

  return( neg_long %>% bind_rows(pos_long) %>% bind_rows(con_long) )
}

mvc_long = to_long(mvc, "mvc")
emg_bi_long = to_long(emg_bi, "emg_bi")
emg_tr_long = to_long(emg_tr, "emg_tr")


feedback_long_df = mvc_long %>% bind_rows(emg_bi_long) %>% bind_rows(emg_tr_long)
# write_csv(feedback_long_df, "datasets/feedback_long_df.csv")
```




let's filter only mvc performance:


```{r}
feedback_long_df = read.csv("datasets/feedback_long_df.csv")
feedback_long_df_mvc  = feedback_long_df %>% filter(activity == "mvc")

str(feedback_long_df_mvc)
```


```{r}
feedback_long_df_mvc = feedback_long_df_mvc %>% mutate_at(vars(id, gender, feedback), factor) %>% mutate_at("time",  as.integer)
```

let's take a look:

```{r}
feedback_long_df_mvc %>% group_by(id, feedback) %>%
  summarise(mean_per = mean(performance)) %>%
  ggplot(aes(x=feedback , y=mean_per, color = feedback)) +
  geom_boxplot() + theme(legend.position = "") +
  labs(title = "Boxplots of Average Relative Performance by Feedback" , y = "Average Relative Performance")
```

```{r}
feedback_long_df_mvc %>% 
  ggplot(aes(x=feedback , y=performance, color = feedback)) +
  geom_boxplot() + theme(legend.position = "") +
  labs(title = "Boxplots of Relative Performance by Feedback" , y = "Relative Performance")
```




```{r}
feedback_long_df_mvc %>% group_by(id, feedback, gender) %>%
  summarise(mean_per = mean(performance)) %>%
  ggplot(aes(x=feedback , y=mean_per, color = gender)) +
  geom_boxplot() +
  labs(title = "Boxplots of Average Relative Performance by Feedback" , y = "Average Relative Performance")
```

```{r}
feedback_long_df_mvc %>% 
  ggplot(aes(x=feedback , y=performance, color = gender)) +
  geom_boxplot() +
  labs(title = "Boxplots of  Relative Performance by Feedback" , y = " Relative Performance")
```




```{r}
feedback_long_df_mvc %>% group_by(feedback, time) %>% summarise(mean_perf = mean(performance)) %>%
  ggplot(aes(x=factor(time) , y=mean_perf, color = feedback, group = feedback)) +
  geom_point() + 
  geom_line(alpha=.5, size = 1.5) + 
  labs(title = "Average Relative Performance by Feedback & Time" ,
       x = "Repetitions" , y = "Average Relative Performance")
```

```{r}
feedback_long_df_mvc %>% group_by(feedback, time, gender) %>% summarise(mean_perf = mean(performance)) %>%
  ggplot(aes(x=factor(time) , y=mean_perf, color = feedback, group = feedback)) +
  geom_point() + 
  geom_line(alpha=.5, size = 1.5) + 
  facet_wrap(~gender) + theme(legend.position = "bottom") +
  labs(title = "Average Relative Performance by Feedback & Time" , subtitle = "Facet by gender",
       x = "Repetitions" , y = "Average Relative Performance")


```


### Analysis


#### Fixed Effect ANOVA

```{r}
fe.fit.1 = aov(performance ~ feedback, data = feedback_long_df_mvc)
anova(fe.fit.1)

```

```{r}
fe.fit.2 = aov(performance ~ feedback + gender , data = feedback_long_df_mvc)
anova(fe.fit.2)
```

**Interaction**

```{r}
ggplot(feedback_long_df_mvc, aes(x = gender, y = performance, group = feedback, col = feedback)) + 
   stat_summary(fun.y = mean, geom = "line", size = 1.5) +
  labs(x= "Gender", y="Relative Performance")
```

```{r}
ggplot(feedback_long_df_mvc, aes(x = feedback, y = performance, group = gender, col = gender)) + 
   stat_summary(fun.y = mean, geom = "line", size = 1.5) +
  labs(x= "Feedback", y="Relative Performance")
```

```{r}
fe.fit.3 = aov(performance ~ feedback * gender , data = feedback_long_df_mvc)
anova(fe.fit.3)
```




```{r, warning=F, message=F , results='hide'}
library("lmerTest")
```

#### One-way Random Effect

let's filter only control group: feedback == "no feedback" 

```{r Feedback One way Random Effect}
nofeed_df =  feedback_long_df_mvc[(feedback_long_df_mvc$feedback=="no feedback"),]

feed_fit.1 =lmer(performance ~  (1 | id) , data = nofeed_df )
```

```{r}
summary(feed_fit.1)
```

```{r}
# profile likelihood condifence intervals
confint(feed_fit.1, oldNames=FALSE)
```

**Checking Assumptions**

```{r Feedback One way Random Effect  checking assumptions}

feed_fit.1.resid = resid(feed_fit.1)
feed_fit.1.predict = predict(feed_fit.1)
plot(feed_fit.1.predict,feed_fit.1.resid, xlab = "Fitted Values", ylab = "Residuals")
abline(0,0, col = "red")
```

```{r Feedback One way Random Effect  checking assumptions2}
sresid = feed_fit.1.resid / sd(feed_fit.1.resid)
qqnorm(sresid) 
abline(a = 0 , b = 1) # add line: y=x 

# we can also use: 
# qqline(sresid)
```

We can also get “estimates” (conditional means) of the random effects $a_i$ with the function `ranef`:
```{r Feedback One way Random Effect  checking assumptions3}
sresid_ref = ranef(feed_fit.1)$id[, "(Intercept)"] / sd(ranef(feed_fit.1)$id[, "(Intercept)"])
qqnorm(sresid_ref, main = "Random effects")
qqline(sresid_ref)
```


#### Mixed Model Random Effect 

##### **Option 1**  
We will continue with the data for the control group only and add a fixed effect to gender:
```{r}
feed_fit.2 = lmer(performance ~  gender + (1 | id) , data = nofeed_df )
anova(feed_fit.2)
summary(feed_fit.2)
```


```{r}
confint(feed_fit.2, oldNames=FALSE)
```

**Checking Assumptions**

```{r Feedback  Mixed Model Random Effect  checking assumptions }
resid_plot = function(model) {
  model.resid = resid(model)
  model.predict = predict(model)
  plot(model.predict,model.resid, xlab = "Fitted Values", ylab = "Residuals")
  abline(0,0, col = "red")
}

qq_plot = function(model) {
  s.model.resid = resid(model) / sd(resid(model))
  qqnorm(s.model.resid)
  qqline(s.model.resid)
}

qq_plot_ref = function(model) {
  sresid = (ranef(model)$id[, "(Intercept)"]) / sd((ranef(model)$id[, "(Intercept)"]))
  qqnorm(sresid, main = "Random effects")
  qqline(sresid)
}

```


```{r Feedback  Mixed Model Random Effect  checking assumptions 1}
resid_plot(feed_fit.2)
qq_plot(feed_fit.2)
qq_plot_ref(feed_fit.2)
```




##### **Option 2**

```{r}
feed_fit.3 = lmer(performance ~ feedback +  (1 | id) , data = feedback_long_df_mvc)
anova(feed_fit.3)
summary(feed_fit.3)
```


```{r}
confint(feed_fit.3, oldNames=FALSE)
```

**Checking Assumptions**

```{r Feedback  Mixed Model Random Effect  checking assumptions 2 }
resid_plot(feed_fit.3)
qq_plot(feed_fit.3)
qq_plot_ref(feed_fit.3)
```



##### **Option 3**



```{r}
feed_fit.4 = lmer(performance ~ feedback + gender + feedback * gender  + (1 | id) , data = feedback_long_df_mvc)
anova(feed_fit.4)
summary(feed_fit.4)
```

```{r}
feedback_long_df
```

**Checking Assumptions**

```{r Feedback Mixed Model Random Effect checking assumptions 3 }
resid_plot(feed_fit.4)
qq_plot(feed_fit.4)
qq_plot_ref(feed_fit.4)
```