---
title: "Chapter 13: Factorial designs"
author: "Andy Field"
date: "2025-05-06"
date-modified: today
format: 
  html:
    theme: journal
    embed-resources: true
    toc: true
    code-overflow: wrap
execute: 
    message: false
    warning: false
---

```{r}
#| label: setup

library(easystats)
library(rstanarm)
library(tidyverse)


# load data
goggles_tib <- here::here("data/goggles.csv") |>
  read_csv() |>
  mutate(
    facetype = as_factor(facetype),
    alcohol = as_factor(alcohol)
  )

# or load from the discovr package
# goggles_tib <- discovr::goggles

```


We want the levels of `facetype` to be ordered unattractive then attractive and the levels of `alcohol` should be placebo, low dose, high dose. We can make sure this is the case using the code below.

```{r}
#| label: chk_levels

# ensure that the factors are ordered as we want them
goggles_tib <- goggles_tib |>
    mutate(
    facetype = fct_relevel(facetype, "Unattractive"),
    alcohol = fct_relevel(alcohol, "Placebo", "Low dose", "High dose")
  )

```


## Self-test: Fitting a linear model

```{r}
#| label: st_lm
#| 

gogg_lm <- goggles_tib |>
  filter(alcohol != "Low dose") |>
  lm(attractiveness ~ facetype*alcohol, data = _)

model_parameters(gogg_lm) |> 
  display(digits = 3)
```


## Linear models for factorial designs

### Plot the data (self-test)

```{r}
#| label: st_plot

ggplot(goggles_tib, aes(x = alcohol, y = attractiveness, colour = facetype, shape = facetype)) +
  stat_summary(fun.data = "mean_cl_normal", geom = "pointrange", position = position_dodge(width = 0.2)) +
  coord_cartesian(ylim = c(0,8)) +
  scale_y_continuous(breaks = 0:8) +
  scale_colour_viridis_d(begin = 0.3, end = 0.85) +
  labs(x = "Alcohol consumption", y = "Attractiveness (0-10)", colour = "Facial stimulus", shape = "Facial stimulus") +
  theme_minimal()
```

### Summary statistics (self-test)

```{r}
#| label: st_describe

goggles_tib |> 
  group_by(alcohol, facetype) |> 
  describe_distribution() |> 
  data_remove(Variable) |> # optional line: removes a column for the purpose of my screenshot! 
  display()
```


## Fitting the model using `afex::aov_4()`

```{r}
#| label: afex

# Raw analysis:
goggles_afx <- afex::aov_4(attractiveness ~ facetype*alcohol + (1|id), data = goggles_tib)
# extract and display the results
model_parameters(goggles_afx, es_type = "omega", ci = 0.95) |> 
  display(use_symbols = TRUE)

# Correcting *p*-values:
model_parameters(goggles_afx, es_type = "omega", ci = 0.95, p_adjust = "bonferroni") |> 
  display(use_symbols = TRUE)

# Plotting the data

afex::afex_plot(goggles_afx, "alcohol", "facetype", legend_title = "Face type") +
  labs(x = "Alcohol consumption", y = "Attractiveness rating (0-10)") +
  theme_minimal()
```


## Fitting the model using `lm()`


```{r}
#| label: lm
#Set contrasts for `facetype`:
  
att_vs_unatt <- c(-1/2, 1/2)
contrasts(goggles_tib$facetype) <- att_vs_unatt

#Set contrasts for `alcohol`:
alcohol_vs_none <- c(-2/3, 1/3, 1/3)
high_vs_low <- c(0, -1/2, 1/2)
contrasts(goggles_tib$alcohol) <- cbind(alcohol_vs_none, high_vs_low)

# Check the contrasts:
contrasts(goggles_tib$facetype)
contrasts(goggles_tib$alcohol)

# Built up the model
incpt_lm <- lm(attractiveness ~ 1, data = goggles_tib) # intercept only
face_lm <- lm(attractiveness ~ facetype, data = goggles_tib) # adds facetype 
alcohol_lm <- lm(attractiveness ~ facetype + alcohol, data = goggles_tib) # adds alcohol
goggles_lm <- lm(attractiveness ~ facetype*alcohol, data = goggles_tib) # add interaction

# compare models
test_wald(incpt_lm, face_lm, alcohol_lm, goggles_lm) |> 
  display()


# should you wish to mimic afex:
car::Anova(goggles_lm, type = 3) |> 
  model_parameters(es_type = "omega", ci = 0.95, p_adjust = "bonferroni") |> 
  display()
```

## Interpreting effects

### Violin plot of facetype (self-test)

Plot the main effect of type of face

```{r}
#| label: st_facetype

ggplot(goggles_tib, aes(x = facetype, y= attractiveness)) +
  geom_violin(colour = "#999933", fill = "#DDCC77", alpha = 0.3) +
  stat_summary(fun.data = "mean_cl_normal", colour = "#999933") +
  scale_y_continuous(breaks = 0:8) +
  labs(y = "Attractiveness (out of 10)", x = "Type of face") +
  theme_minimal()
```

Getting the estimated marginal means for the model created using `afex::aov_4()`

```{r}
#| label: facetype_emm

# Predicted means for facetytype (afex)
estimate_means(goggles_afx, "facetype") |> 
  display()

# Predicted means for facetytype (lm)
estimate_means(goggles_lm, "facetype") |> 
  display()
```

### Violin plot of attractiveness (self-test)

Plot the main effect of alcohol consumption

```{r}
#| label: st_alcohol

ggplot(goggles_tib, aes(x = alcohol, y= attractiveness)) +
   geom_violin(colour = "#117733", fill = "#44AA99", alpha = 0.3) +
  stat_summary(fun.data = "mean_cl_normal", colour = "#117733") +
  scale_y_continuous(breaks = 0:8) +
  labs(y = "Attractiveness (out of 10)", x = "Dose of alcohol") +
  theme_minimal()
```


Get estimated marginal means

```{r}
#| label: alcohol_emm

# Get estimated marginal means from the `afex' model
estimate_means(goggles_afx, "alcohol") |> 
  display()

# Get estimated marginal means from the `lm' model
estimate_means(goggles_lm, "alcohol") |> 
  display()

```

### The interaction effect

Plot the interaction using `afex`

```{r}
#| label: int_afx

afex::afex_plot(object = goggles_afx,
                x = "alcohol",
                trace = "facetype",
                legend_title = "Face type") + 
  labs(x = "Alcohol consumption", y = "Attractiveness rating (0-10)") +
  theme_minimal()
```

Get estimated marginal means 

```{r}
#| label: emm_int

# Get estimated marginal means from the `afex` model
estimate_means(goggles_afx, c("alcohol", "facetype")) |> 
  display()

# Get estimated marginal means from the `lm` model
estimate_means(goggles_lm, c("alcohol", "facetype")) |> 
  display()
```

### Contrasts

View the contrasts set for the model fitted with `lm()`

```{r}
#| label: contrasts
#| 
model_parameters(goggles_lm) |>
  display()
```

### Simple effects

To look at the effect of facetype at each level of alcohol, we'd execute:

```{r}
#| label: se_alcohol

# afex model
estimate_contrasts(goggles_afx, contrast = "facetype", by = "alcohol", comparison = "joint", p_adjust = "bonferroni") |> 
  display()

# lm model
estimate_contrasts(goggles_lm, contrast = "facetype", by = "alcohol", comparison = "joint", p_adjust = "bonferroni") |> 
  display()
```


To look at the effect of alcohol at each level of facetype, we'd execute:

```{r}
#| label: se_facetype

# afex model
estimate_contrasts(goggles_afx, contrast = "alcohol", by = "facetype", comparison = "joint", p_adjust = "bonferroni") |> 
  display()

# lm model
estimate_contrasts(goggles_lm, contrast = "alcohol", by = "facetype", comparison = "joint", p_adjust = "bonferroni") |> 
  display()

```

### *Post hoc* tests

```{r}
#| label: post_hoc

# afex model
estimate_contrasts(goggles_afx, contrast = "alcohol", by = "facetype", p_adjust = "bonferroni") |> 
  display()

# lm model
estimate_contrasts(goggles_lm, contrast = "alcohol", by = "facetype", p_adjust = "bonferroni") |> 
  display()

```



## Self test: Diagnostic plots 

Plot residuals from the `goggles_lm` model. Remember to execute `library(ggfortify)` before running this code.

```{r}
#| label: resid_plot


# afex model
check_model(goggles_afx)
# lm model
check_model(goggles_lm)

```

## Robust models

### Robust parameter estimates

Fit a model using `lmRob()`

```{r}
#| label: lmrob_gogg

goggles_rob <- robust::lmRob(attractiveness ~ facetype*alcohol, data = goggles_tib)
summary(goggles_rob)
```

### Robust standard errors

HC4 standard errors

```{r}
#| label: hc4

model_parameters(goggles_lm, vcov = "HC4") |> 
  display()
```

### Bootstrap standard errors

```{r}
#| label: boot

model_parameters(goggles_lm, bootstrap = TRUE) |> 
  display()
```

### Robust overall tests

Fit a model based on 20% trimmed means.

```{r}
#| label: wrs2

WRS2::t2way(attractiveness ~ alcohol*facetype, data = goggles_tib)
goggles_mcp2atm <- WRS2::mcp2atm(attractiveness ~ alcohol*facetype, data = goggles_tib)
goggles_mcp2atm
goggles_mcp2atm$contrasts
```

Fit a model based on an M-estimator.

```{r}
#| label: wrs2_more


WRS2::pbad2way(attractiveness ~ alcohol*facetype, data = goggles_tib, nboot = 1000)
goggles_mcp2a <- WRS2::mcp2a(attractiveness ~ alcohol*facetype, data = goggles_tib, nboot = 1000)
goggles_mcp2a
goggles_mcp2a$contrasts
```

## Bayesian model and Bayes factors

```{r}
# Bayesian model with default priors and bayes factors
goggles_rtb <- stan_glm(attractiveness ~ alcohol*facetype, data = goggles_tib)

model_parameters(goggles_rtb, test = c("pd", "bf"), null = 0) |> 
  display()
```




