---
title: Alternatives to t-tests
author: J. Taylor
date: Spring 2024
engine: knitr
format:
  revealjs:
    smaller: true
    scrollable: true
    html-math-method: mathjax
    header-logo: https://web.stanford.edu/class/stats191/_static/logo.png
    header-logo-link: "https://stats191.stanford.edu"
    css: ../revealstyle.css
webr:
  show-startup-message: false
filters:
  - webr
  - reveal-header
jupyter:
  jupytext:
    formats: 'qmd:quarto,ipynb'
    text_representation:
      extension: .qmd
      format_name: quarto
      format_version: '1.0'
      jupytext_version: 1.16.1
  kernelspec:
    display_name: R
    language: R
    name: ir
---

## Outline 

- Case studies:

    1. Shuttle O-ring incidents

    2. Cognitive load in math problems

- Robustness and resistance of two-sample $t$-tests

- Transformations


```{webr-r}
#| context: setup
require(ggplot2)
set.seed(0)
```

## Case study A: O-rings failing to seal fuel

- Is there a difference between <font color="Blue">Cool(<65F)</font> and <font color="Red">Warm (>65F)</font>?

- Temperature on launch day: 29F.

```{webr-r}
orings = read.csv('https://raw.githubusercontent.com/StanfordStatistics/stats191-data/main/Sleuth3/orings.csv', header=TRUE)
orings$Incidents[orings$Launch == 'Cool']
orings$Incidents[orings$Launch == 'Warm']
```


## What test to use?

- Samples are really small

- No sense in which they could be normally distributed...

### Permutation test

- Decide a test statistic $T$ comparing <font color="Blue">Cool(<65F)</font> to <font color="Red">Warm (>65F)</font>. Could be the two-sample $t$ test statistic

- Under $H_0:$ <font color="Blue">Cool(<65F)</font> has the same distribution <font color="Red">Warm (>65F)</font> suppose we:

     1. Shuffle the response in the two samples by a random ordering $\sigma$ (called a permutation)

     2. Recompute the statistic yielding $T(\sigma)$

## Permutation test in practice

### Build up a reference distribution under $H_0$

```{webr-r}
B = 10000 # n_permutations
null_stats = rep(NA, B)
N = nrow(orings)
for (i in 1:B) {
    idx = sample(N, N, replace=FALSE)
    orings_star = data.frame(Incidents=orings$Incidents[idx],
                             Launch=orings$Launch)
    null_stats[i] = t.test(Incidents ~ Launch, data=orings_star)$stat
}
```


## Histogram of $t$-statistic under $H_0$

- Doesn't look $t$-shaped at all... but test is valid. <font color="red">Why?</font>

- <font color="red">Under $H_0$: shuffling the response doesn't change the distribution!</font>

```{webr-r}
hist(null_stats)
```

---

### Compute a $p$-value

```{webr-r}
observed = t.test(Incidents ~ Launch, data=orings)$stat
observed
mean(null_stats > observed)
```

## Rank Sum Test

- One can use *any* test statistic for permutation test...

- Only tests null $H_0$: distributions are the same.

- What should we use for $H_a$?


## Case study B: difference in cognitive load

- Time to solve math problems under two different presentations

- Small sample size, not normal?

```{webr-r}
cognitive_load = read.csv('https://raw.githubusercontent.com/StanfordStatistics/stats191-data/main/Sleuth3/cognitive_load.csv', header=TRUE)
cognitive_load
```

## Rank sum test

- Rank the outcome

- $T$ = Sum the ranks in one of the groups (14 of 28 are `Modified`)

- Under $H_0:$ distributions are identical, $T$ has distribution the sum of 14 ranks in a random permutation...

```{webr-r}
cognitive_load$R = rank(cognitive_load$Time)
rank_sum = with(cognitive_load, sum(R[Treatment == 'Modified']))
rank_sum
W = 14^2 + 14*15/2 - rank_sum
c(W, 14*14-W)
```

## Rank sum test

- Using a builtin test

```{webr-r}
wilcox.test(Time ~ Treatment,
            data=cognitive_load,
            alternative='greater')
```

## Rank sum test

### Null distribution (ignoring ties)

```{webr-r}
null_sample = runif(28)
null_treat = c(rep('M', 14), rep('U', 14))
null_df = data.frame(S=null_sample, T=null_treat)
null_df$R = rank(null_df$S)
with(null_df, sum(R[T=='M']))
```

- `R` uses a distributional approximation...

---

### Confidence interval

```{webr-r}
wilcox.test(Time ~ Treatment,
            data=cognitive_load,
            alternative='two.sided',
            conf.int=TRUE)
```

- What parameter is this estimating?


## Shift alternative

![](http://www.stanford.edu/class/stats191/Chapter4/figs/shift_alternative.png){height=400 fig-align="center"}

- <font color="green">$H_a:$</font> the difference between groups is a simple shift

- Other interpretations: estimating median of difference...


## Paired data

- There are plenty of other non-parametric tests out there

### Signed rank test for paired data (i.e. `schizophrenia`)

- <font color="red">$H_0:$</font> distributions are
the same in each group..

- $\implies$ **differences are symmetrically distributed around 0!**

## Symmetric null

![](http://www.stanford.edu/class/stats191/Chapter4/figs/sign_null.png){height=400 fig-align="center"}

- <font color="red">$H_0:$</font> distribution is symmetric around 0.

## Symmetric alternative

![](http://www.stanford.edu/class/stats191/Chapter4/figs/sign_alternative.png){height=400 fig-align="center"}

- <font color="green">$H_a:$</font> distribution is symmetric around $\delta$.

---

### Example of signed rank `schizophrenia`

```{webr-r}
schizophrenia = read.csv('https://raw.githubusercontent.com/StanfordStatistics/stats191-data/main/Sleuth3/schizophrenia.csv', header=TRUE)
with(schizophrenia, wilcox.test(Affected, Unaffected), paired=TRUE)
```

### Compare to `t.test`

```{webr-r}
with(schizophrenia, t.test(Affected - Unaffected))
```