---
title: "The Modelisation Approach"
output: 
  rmarkdown::html_vignette:
    toc: true
    fig_width: 10.08
    fig_height: 6
tags: [r, estimate, estimate link, predictions]
vignette: >
  %\VignetteIndexEntry{The Modelisation Approach}
  \usepackage[utf8]{inputenc}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
bibliography: bibliography.bib
---

```{r warning=FALSE, message=FALSE, include=FALSE}
library(knitr)
options(knitr.kable.NA = "")
knitr::opts_chunk$set(
  comment = ">",
  warning = FALSE,
  message = FALSE,
  dpi = 300
)
options(digits = 2)

if (!requireNamespace("ggplot2", quietly = TRUE) ||
  !requireNamespace("see", quietly = TRUE) ||
  !requireNamespace("gganimate", quietly = TRUE) ||
  !requireNamespace("rstanarm", quietly = TRUE) ||
  !requireNamespace("datawizard", quietly = TRUE)) {
  knitr::opts_chunk$set(eval = FALSE)
}

set.seed(333)
```

For a long time, applied statistics, especially in psychology, revolved around
the idea of assessing differences between groups, using ***t*-tests** or
**ANOVAs**, or by investigating the existence of a significant linear
relationship between two continuous variables using **correlations**.

Every psychology student heard that *"ANOVAs, t-tests and correlations are all
linear models"* (aka, linear regressions). And thus, one can only wonder why
these different tests have not been yet replaced by a unified model.

```{r echo=FALSE, fig.align='center', out.width="50%"}
knitr::include_graphics("https://github.com/easystats/easystats/raw/main/hextools/modelbased/allregressions.png")
```

The tests mentioned above have major advantages. For instance, they are very easy to
compute (one can even figure them out by hand). This partially explains their
historical popularity, as these tests became the norm at an age when computers
and fancy programs weren't available. Moreover, they are easy to visualize (or
so it is *believed*, as it's not always true - but we'll come to that later)
and, most importantly, they have straightforward (or at least, *conventionally
accepted*) interpretations. After all, most of the time, statistics are used
with the aim to draw some conclusions from the data.

But here are some problems with this approach:

1. We don't know whether the model lying underlying the test (the *t*-test, the
   ANOVAs or the correlation) is any good to "represent" the data, preventing in
   turn to adjust the confidence in the results. For instance, we know that reaction
   times (RTs) have a non-normal distribution, and yet we often see the use of 
   unsuited methods, like outliers-removal, scale-transformations, and computation
   of empirical scores like "mean" or "SD" that do not represent the data well.

2. The visualisation are usually based on the data, rather than on the model
   (which is problematic in the case of complex models).

3. New indices, such as *group means*, *contrasts*, etc., are also derived from
   the data, instead of using the model.

Thus, we will describe below another approach to statistics, centered around
statistical modelisation.

# The Empirical Approach (Classic)

```{r}
library(datawizard)
library(emmeans)
library(parameters)
library(modelbased)
library(ggplot2)
```

## Data Simulation


First we will run the function below to simulate some data. There's no need to
understand the hows and whys, as we will explain everything in due time.

<details> <summary>Click here to see the code to generate the data.</summary>

```{r}
generate_data <- function(effect = 5, noise = 0.5) {
  data <- data.frame()
  n <- 100
  for (i in 1:length(effect)) {
    participant <- data.frame(Experimental_Variable = c(seq(-3, 3, length = n / 2), seq(-3, 3, length = n / 2)))
    participant$RT <- c(participant$Experimental_Variable[1:(n / 2)]**2 - effect[i], (participant$Experimental_Variable[(n / 2 + 1):n] + effect[i])) + rnorm(n, 0, abs(noise[i]))
    participant$Condition <- rep(c("A", "B"), each = n / 2)
    participant$Participant <- paste0("S", i)
    data <- rbind(data, participant)
  }
  data$RT <- (100 + data$RT) * 10
  data
}

data <- generate_data(effect = rnorm(30, 2, 2), noise = rnorm(30, 0, 0.4))
#
# library(rtdists)
#
# data <- data.frame(
#   Participant = paste0("S", speed_acc$id),
#   Item = as.numeric(speed_acc$stim),
#   Condition = speed_acc$condition,
#   Correct = ifelse(as.character(speed_acc$stim_cat) == as.character(speed_acc$response), 1, 0),
#   RT = speed_acc$rt * 1000
# )
```

</details>


## Preprocessing 


## ANOVAs 

In ANOVAs, it's all about groups. Even though people also add in continuous
variables (creating these ANCOVAs, MANOVAs and other monstrosities), it's not
really "in their spirit": ANOVAs were made to compare groups.

So we will take, for each participant, its 20...

**TO DO: COMPLETE THIS VIGNETTE**.


```{r}
data_anova <- data
data_anova$Category <- recode_into(
  Experimental_Variable < -1.5 ~ "Low",
  Experimental_Variable > 1.5 ~ "High",
  default = "Middle",
  data = data_anova
)
data_anova$Category <- factor(data_anova$Category, levels = c("Low", "Middle", "High"))

data_anova <- data_anova |>
  data_group(c("Participant", "Condition", "Category")) |>
  data_summary(RT = mean(RT))

results <- aov(RT ~ Condition * Category + Error(Participant), data = data_anova)
parameters(results)
```

What can we conclude from that? Absolutely **nothing**! We need to investigate in more details, for instance by running post-hoc comparison tests. This uninformativeness is one reason why ANOVA should be banned from psychological science.

## Post-hoc comparison tests

```{r}
posthoc <- model_emmeans(results, by = c("Condition", "Category")) |>
  pairs()
parameters(posthoc)
```

```{r}
data_anova |>
  ggplot(aes(x = Category, y = RT, fill = Condition)) +
  geom_boxplot() +
  see::theme_modern()
```


# The Modelisation Approach

A model is made of parameters, which have 'real' meaning, as opposed to indices
of significance (which are abstract).

## Summary 

## 1. Draw what you want to visualize

We can use `geom_smooth()`, which can fit non-linear relationships in an empirical way, to give us an idea of the shape of the relationships.

```{r}
data |>
  data_group(c("Participant", "Condition")) |>
  ggplot(aes(x = Experimental_Variable, y = RT, color = Condition)) +
  geom_jitter(alpha = 0.4) +
  geom_smooth(method = "loess", se = FALSE) +
  see::theme_modern()
```



## 2. Make models for it



## 3. Select the best Model

## 4. Visualize the best model

## 5. Investigate its Parameters

# References

<!-- ESTIMATION OF INTERINDIVIDUAL VARIABILITY see Williams, D. R., Mulder, J., Rouder, J. N., & Rast, P. (2020). Beneath the surface: Unearthing within-person variability and mean relations with Bayesian mixed models. Psychological Methods. -->
