2024-04-03_UseAndAbuseANOVA.qmd

---
title-block-banner: true
title: "Avoiding abuse and misuse of T-test and ANOVA: Regression for categorical responses"
subtitle: "An approach using Bayesian regression with brms"
date: today
date-format: full
author: 
  - name: "Daniel Manrique-Castano"
    orcid: 0000-0002-1912-1764
    degrees:
      - PhD
    affiliation: 
      - name: Univerisity Laval 
        department: Psychiatry and Neuroscience
        group: Laboratory of neurovascular interactions 
note: "GitHub: https://daniel-manrique.github.io/"
keywords: 
  - Ordinal regression distribution
  - Bayesian modeling
  - brms 
   
license: "CC BY"

format:
   pdf: 
    toc: true
    number-sections: true
    colorlinks: true
   html:
    code-fold: true
    embed-resources: true
    toc: true
    toc-depth: 2
    toc-location: left
    number-sections: true
    theme: spacelab

knitr:
  opts_chunk: 
    warning: false
    message: false
    
csl: science.csl
bibliography: References.bib
---

In neuroscience and other biomedical sciences, it is common to use behavioral tests to assess responses to experimental conditions or treatments. We can assess many aspects, from basic motor and exploratory behaviors to memory and learning. Many of these variables are continuous (numerical) responses, i.e. they can take (finite or infinite) values in a given range. Time in the open field, animal weight, or the number of cells in a brain region can generally be thought of as continuous numerical variables.

However, there are other types of variables that we researchers record in our experiments. A very common one is ordered categorical variables, also called ordinal variables. These are categorical variables that have a natural order, analogous to the well-known surveys in which we answer whether we agree or disagree with a statement, with 0 being strongly disagree and 5 being strongly agree. To facilitate the recording of these variables in printed or digital datasheets, we codify them (by convention) as numbers. This is the case of the 5-point Bederson score used in the context of cerebral ischemia research in rodent models [@bieber2019], which is coded as follows:

-   0 = no observable deficit

-   1 = forelimb flexion

-   2 = forelimb flexion and decreased resistance to lateral push

-   3 = circling

-   4 = circling and spinning around the cranial-caudal axis

-   5 = no spontaneous movement.

Note that the numbers are just simple conventions. We could also use a,b,c,d,e; excellent, good, not so good, bad, very bad, almost dead; etc. I really do not think it is presumptuous to stress the self-evident nature of the matter. Surprisingly, however, [@fig-sample; @onufriev2021; @liu2022] represents a particular and widespread malpractice that scientists in this field have been practicing for years: They use t-tests and ANOVAS to analyze such ordered categorical variables.

![Left: Neurological score by Onufriev et. al (2021) (CC-BY). Right: Neurological score by Liu et al. (2022) (CC-BY).](Plots/2024-04-03_UseAndAbuseANOVA/Test%20example.png){#fig-sample fig-align="center"}

I still cannot find a logical explanation why dozens of authors, reviewers, and editors feel comfortable with this scenario. Not only is it not logical, but more importantly, it leads to biased effect sizes, low detection rates, and type I error rates, among other things [@bürkner2019].

When I have played the role of reviewer 2 and emphasized this point to the authors, asking them to explain why they evaluate a categorical response with a statistical test designed to deal with continuous numerical variables, what I get is a long list of published articles that follow this irrational practice. So I finally found an answer to why they do it: what Gerd Gigerenzer (2004) calls "the ritual of mindless statistics" [@gigerenzer2004]. In fact, most of us scientists have little idea what we are doing with our data, and are simply repeating common malpractices that have been passed down from generation to generation.

In this article, we will then look at a more viable alternative for analyzing ordered categorical variables using `R`, `brms` [@brms-2]and elements from the `tidyverse` [@tidyverse].

# Facing the ritual of mindless statistics

To avoid the ritual of mindless statistics, we'll recreate the data set from Liu et. al (2022) by eyeballing the data points and organizing them in a table:

```{r}
#| label: CrateData
#| include: true
#| warning: false
#| message: false
#| column: margin

# We Define the observations
observations <- list(
  Sham = c(rep(0, 7)),
  tMCAO = c(rep(2, 3), rep(3, 2), rep(4, 2)),
  tMCAO_C = c(rep(1, 3), rep(2, 3), rep(3, 1))
)

# We create an empty data frame to populate
df <- data.frame(Group = character(), Response = integer())

# We populate the data frame
for (group in names(observations)) {
  for (response in unique(observations[[group]])) {
    df <- rbind(df, data.frame(Group = rep(group, sum(observations[[group]] == response)), 
                               Response = rep(response, sum(observations[[group]] == response))))
  }
}

head(df)
```

If you look at the R code at this point, you will notice that the variable `Response' is identified as a number ranging from 0 to 4. Of course, we are fully aware that this response is not numeric, but an ordered (ordinal) categorical variable. So we do the conversion explicitly:

```{r}
#| label: ConvertData
#| include: true
#| warning: false
#| message: false
#| column: margin

df$Response <- factor(df$Response, levels = c("0", "1", "2", "3", "4"), ordered = TRUE)

str(df)

```

Now we can verify that it is recognized as an ordered categorical variable. With this in hand, we can easily move on to visualization and modeling.

## Exploratory data visualization

First, let's load the necessary libraries and create a visual theme for our plots.

```{r}
#| label: LoadPack
#| include: true
#| warning: false
#| message: false

library(ggplot2)
library(brms)
library(ggdist)
library(easystats)
library(dplyr)
library(tibble)

Plot_theme <- theme_classic() +
  theme(
      plot.title = element_text(size=18, hjust = 0.5, face="bold"),
      plot.subtitle = element_text(size = 10, color = "black"),
      plot.caption = element_text(size = 12, color = "black"),
      axis.line = element_line(colour = "black", size = 1.5, linetype = "solid"),
      axis.ticks.length=unit(7,"pt"),
     
      axis.title.x = element_text(colour = "black", size = 16),
      axis.text.x = element_text(colour = "black", size = 16, angle = 0, hjust = 0.5),
      axis.ticks.x = element_line(colour = "black", size = 1),
      
      axis.title.y = element_text(colour = "black", size = 16),
      axis.text.y = element_text(colour = "black", size = 16),
      axis.ticks.y = element_line(colour = "black", size = 1),
      
      legend.position="right",
      legend.direction="vertical",
      legend.title = element_text(colour="black", face="bold", size=12),
      legend.text = element_text(colour="black", size=10),
      
      plot.margin = margin(t = 10,  # Top margin
                             r = 2,  # Right margin
                             b = 10,  # Bottom margin
                             l = 10) # Left margin
      ) 
```

A simple way to visualize categorical data is with a bar graph.

```{r}
#| label: fig-Fig2
#| include: true
#| warning: false
#| message: false
#| fig-cap: Response colored by group.
#| fig-height: 5
#| fig-width: 6

ggplot(df, aes(x = factor(Response), fill = Group)) +
  geom_bar(position = "dodge") +
  labs(x = "Response", y = "Count", title = "Response by Group") +
  theme_minimal() +
  scale_fill_brewer(palette = "Set1") +
  Plot_theme +
theme(legend.position = "top", legend.direction = "horizontal")
```

@Fig-Fig2 shows the frequency per group. This is more in line with the logic of a categorical variable than box plots, which are more relevant for continuous numerical variables. Now let's run a regression to unravel the mysteries (if any) of this data set.

# Fitting an Ordinal Regression with brms

We'll use `brms` to fit a cumulative model. This model assumes that the neurological score Y is derived from the categorization of a (presumably) latent (but not observable or measured) continuous variable Y˜ [@bürkner2019]. As usual in most `brms` tutorials, I must apologize for skipping the "priors" issue. Let's assume an "I have no idea" mentality and let the default (flat) `brms` priors do the dirty work.

We fit the `cumulative` model by following the usual formula syntax and adding `cumulative("probit")` as a family (assuming the latent variable and the corresponding error term are normally distributed. We have only one predictor variable, the experimental group to which each animal belongs.

```{r}
#| label: OrdinalFit
#| include: true
#| warning: false
#| message: false
#| results: false

Ordinal_Fit <- brm(Response ~ Group, 
           data = df, 
           family = cumulative("probit"),
           # seed for reproducibility purposes
           seed = 8807,
           control = list(adapt_delta = 0.99),
           # this is to save the model in my laptop
           file    = "Models/2024-04-03_UseAndAbuseANOVA/Ordinal_Fit.rds",
           file_refit = "never")

# Add loo for model comparison
Ordinal_Fit <- 
  add_criterion(Ordinal_Fit, c("loo", "waic", "bayes_R2"))
```

Before we look at the results, let's do some model diagnostics to compare the observations and the model predictions.

## Model diagnostics

```{r}
#| label: fig-OrdinalDiag1
#| include: true
#| warning: false
#| message: false
#| fig-cap: Model diagnostics for the ordinal regression
#| fig-height: 4
#| fig-width: 5

set.seed(8807)

pp_check(Ordinal_Fit, ndraws = 100) +
  labs(title = "Ordinal regression") +
  theme_classic()
```

From @fig-OrdinalDiag1, I speculate that the uneven distribution (variance) of scores across groups may cause such deviations. Later, we'll see if predicting the variance yields better estimates. For now, we are good to go.

## Checking the results for the ordinal regression

Let's take a look at the posterior distribution using the `describe_posterior` function from the `bayestestR` package [@bayestestR], as an alternative to the typical `summary` function.

```{r}
#| label: Student_Posterior
#| include: true
#| warning: false
#| message: false

describe_posterior(Ordinal_Fit,
                   centrality = "mean",
                   dispersion = TRUE,
                   ci_method = "HDI",
                   test = "rope",
                   )
```

The thresholds (score thresholds) are labeled as "intercepts" in this model, which apply to our baseline "sham" condition. The coefficients for GrouptMCAO and GrouptMCAO_C indicate the difference from sham animals on the latent Y˜ scale. Thus, we see that GrouptMCAO has 8.6 standard deviations higher scores on the latent Y˜ scale.

I want to point out a crucial (and not trivial) aspect here as a reminder for a future post (stay tuned). In my opinion, it is irrelevant to make comparisons between a group that does not have a distribution (which is 0 in all cases), such as the sham group, and a group that does have a distribution (the experimental groups). If you think about it carefully, the purpose of this procedure is null. But let us follow the thread of modern scientific culture and not think too much about it.

Of course, appreciating the difference between the groups becomes more feasible if we look at it visually with the `conditional effects' function of `brms'.

```{r}
#| label: fig-OrdinalEff
#| include: true
#| warning: false
#| message: false
#| fig-cap: Conditional effects for the Ordinal model
#| fig-height: 5
#| fig-width: 7


Ordinal_CondEffects <- 
  conditional_effects(Ordinal_Fit, "Group", categorical = TRUE)

Ordinal_CondEffects <- plot(Ordinal_CondEffects, 
       plot = FALSE)[[1]]

Ordinal_CondEffects + 
  Plot_theme +
  theme(legend.position = "bottom", legend.direction = "horizontal")
```

Curiously (or perhaps not so much from the MCMC simulation framework), the model estimates that the dummy group can have values of 1, although we have not seen this in the data. The reason is no doubt that the model estimates a common variance for all groups. We will see if this changes if our model also includes the variance as a response.

Otherwise, the estimates we obtain for the tMCAO and tMCO_C groups are much closer to the data. This allows us to make more precise statements instead of running an ANOVA (incorrectly for ordinal variables) and saying that there is "a significant difference" between one group and the other, which is the same as saying nothing. For example, the model tells us that scores of 2-4 for the tMCAO group have a very similar probability (about 25%). The case is different for the tMCAO_C group, where the probability of 2 in the neurological score is higher (with considerable uncertainty) than for the other scores. If I were confronted with this data set and this model, I would claim that the probability that the tMCAO_C group reflects less neurological damage (based on scores 1 and 2) is higher than for the tMCAO group.

Can we get precise numbers that quantify the differences between the scores of the different groups and their uncertainty? Of course we can! using the `emmeans` package [@emmeans]. But that will be the subject of another post (stay tuned).

# Including the variance as a response variable

For this type of cumulative model, there is no `sigma` parameter. Instead, to account for unequal variances, we need to use a parameter called `disc`. For more on this aspect, see the Ordinal Regression tutorial by Paul-Christian Bürkner, the creator of `brms` [@bürkner2019].

```{r}
#| label: OrdinalFit_Sigma
#| include: true
#| warning: false
#| message: false
#| results: false

Ordinal_Mdl2 <- bf (Response ~ Group) +
                     lf(disc ~  0 + Group, cmc = FALSE)

Ordinal_Fit2 <- brm( 
           formula = Ordinal_Mdl2,
           data = df, 
           family = cumulative("probit"),
           # seed for reproducibility purposes
           seed = 8807,
           control = list(adapt_delta = 0.99),
           # this is to save the model in my laptop
           file    = "Models/2024-04-03_UseAndAbuseANOVA/Ordinal_Fit2.rds",
           file_refit = "never")

# Add loo for model comparison
Ordinal_Fit2 <- 
  add_criterion(Ordinal_Fit2, c("loo", "waic", "bayes_R2"))
```

## Model diagnostics

We perform the model diagnostics as done previously: 

```{r}
#| label: fig-OrdinalDiag2
#| include: true
#| warning: false
#| message: false
#| fig-cap: Model diagnostics for our model predicting the variance
#| fig-height: 4
#| fig-width: 5

set.seed(8807)

pp_check(Ordinal_Fit2, ndraws = 100) +
  labs(title = "Student-t") +
  theme_classic()
```

@fig-OrdinalDiag2 shows that, in fact, my expectations were not met. Including the variance as a response in the model does not improve the fit of the predictions to the data. The trend is maintained, but the predictions vary significantly. Nevertheless, this is a candidate for a generative model.

## Checking the results for our new model

We visualize the posterior distribution for this model:

```{r}
#| label: Ordinal_Posterior2
#| include: true
#| warning: false
#| message: false

describe_posterior(Ordinal_Fit2,
                   centrality = "mean",
                   dispersion = TRUE,
                   ci_method = "HDI",
                   test = "rope",
                   )
```

We can see a meaningful difference in the coefficients compared to the first model. The coefficient for "GrouptMCAO" increases from 8.6 to 15.9 and that for "GrouptMCAO_C" from 7 to 10.8. Undoubtedly, this model gives us a different picture. Otherwise, the variance term is presented under the names "disc_GrouptMCAO" and "disc_GrouptMCAO_C". We can see that both variances are very different from our "sham" baseline.

Let's plot the results:

```{r}
#| label: fig-OrdinalCondEffects2
#| include: true
#| warning: false
#| message: false
#| fig-cap: Conditional effects for our model including the variance as a response
#| fig-height: 5
#| fig-width: 5
#| layout-nrow: 1

Ordinal_CondEffects2 <- 
  conditional_effects(Ordinal_Fit2, categorical = TRUE)

Ordinal_CondEffects2 <- plot(Ordinal_CondEffects2, 
       plot = FALSE)[[1]]

Ordinal_CondEffects2 +
  Plot_theme +
  theme(legend.position = "bottom", legend.direction = "horizontal")

```

Contrary to what I expected, the model still predicts that sham animals have a small probability of scoring 1. What we have confirmed here is that this prediction is not based on the (false) assumption that all groups have the same variance. Nevertheless, it is still a logical prediction (not irracional) within this framework (ordinal regression) based on thresholds. If we go to Richard McElreath's Statistical Rethinking [@mcelreath2020k] we find the same situation with monkeys pulling levers. Fitting a more constrained model will require the use of informative priors. Leave that for a future post. I know I made three promises here, but I will keep them.

In this model, the probabilities for the different outcomes in the tMCAO group shift slightly. However, given the high uncertainty, I will not change my conclusion about the performance of this group based on this model. For the tMCAO_C group, on the other hand, the predictions do not shift in a way that is readily apparent to the eye. Let's conclude this blog post by comparing the two models.

# Model comparison

We do the model comparison using the the `loo` package [@loo; @vehtari2016] for leave-one-out cross validation. For an alternative approach using the WAIC criteria [@gelman2013] I suggest you read [this](https://medium.com/towards-data-science/do-not-over-think-about-outliers-use-a-student-t-distribution-instead-b6c584b91d5c) post also published by Towards Data Science.

```{r}
#| label: Models_Compare
#| include: true
#| warning: false
#| message: false
#| results: false

loo(Ordinal_Fit, Ordinal_Fit2)
```

Under this scheme, the models have very similar performance. In fact, the first model is slightly better for out-of-sample predictions. Accounting for variance did not help much in this particular case, where perhaps relying on informative priors can unlock the next step of scientific inference.

I would appreciate your comments or feedback letting me know if this journey was useful to you. If you want more quality content on data science and other topics, you might consider becoming a [medium member](https://medium.com/membership).

In the future, you can find an updated version of this post on my [GitHub site](https://github.com/daniel-manrique/MediumBlog/blob/main/2024-04-03_UseAndAbuseANOVA.qmd). 

* All images, unless otherwise stated,  were generated using the displayed R code. 

# References

::: {#refs}
:::

```{r}
sessionInfo()
```