---
title: Drawing Statistical Conclusions
author: J. Taylor
date: Spring 2024
engine: knitr
format:
  revealjs:
    html-math-method: mathjax
    fig-width: 4
    header-logo: https://web.stanford.edu/class/stats191/_static/logo.png
    header-logo-link: 'https://stats191.stanford.edu'
    scrollable: true
    smaller: true
    css: ../revealstyle.css
webr:
  packages:
    - ggplot2
  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.15.0
  kernelspec:
    display_name: R
    language: R
    name: ir
---


## Case Study A: Motivation for creative writers 

- Creative writing students randomly assigned to *intrinsic* vs. *extrinsic* priming questionnaires.



```{webr-r}
#| context: setup
require(ggplot2)
creativity = read.csv('https://raw.githubusercontent.com/StanfordStatistics/stats191-data/main/Sleuth3/creativity.csv', header=TRUE)
salaries = read.csv('https://raw.githubusercontent.com/StanfordStatistics/stats191-data/main/Sleuth3/salaries.csv', header=TRUE)
set.seed(0)
```

```{webr-r}
head(creativity)
```

## Summarizing the groups

### Extrinsic Group
```{webr-r}
extrinsic = creativity$Score[creativity$Treatment == 'Extrinsic']
summary(extrinsic)
```

### Intrinsic Group

```{webr-r}
intrinsic = creativity$Score[creativity$Treatment == 'Intrinsic']
summary(intrinsic)
```

## Histogram of `Score` stratified by `Sex`

```{webr-r}
#| context: output
# this plot is for visualization only
# students not expected to reproduce
fig <- (ggplot(creativity, aes(x=Score, fill=Treatment)) +
        geom_histogram(aes(y=after_stat(density)),
	color="#e9ecef",
	alpha=0.6,
	position='identity',
	bins=10) +
        labs(fill=""))
fig
```

## Case Study B: Difference in salaries between male and female employees

- Salaries from Harris Trust and Bank over years 1969-1977

```{webr-r}
head(salaries)
```

### Females

```{webr-r}
female = salaries$Salary[salaries$Sex == 'Female']
summary(female)
```

### Males

```{webr-r}
male = salaries$Salary[salaries$Sex == 'Male']
summary(male)
```

## Histogram of `Salary` stratified by `Sex`

```{webr-r}
#| context: output

# this plot is for visualization only
# students not expected to reproduce

fig <- (ggplot(salaries, aes(x=Salary, fill=Sex)) +
        geom_histogram(aes(y=after_stat(density)),
                       color="#e9ecef",
                       alpha=0.6,
                       position='identity',
                       bins=10) +
        labs(fill=""))
fig
```

## Boxplot of `Salary` stratified by `Sex`

```{webr-r}
boxplot(Salary ~ Sex, data=salaries, col='orange')
```


## Key differences between the studies 

- Creative writing study was a *randomized experiment*.

- Salary dataset was an *observational study*.

### Implications

- Differences in strength of conclusions: randomized experiments like `creativity` can admit *causal* conclusions

- Generalizability: what population are the data from?

    1. If we consider Harris a typical bank, then `salaries` represents a sample of starting salaries.

## Making statistical inferences

![](http://www.stanford.edu/class/stats191/Chapter1/figs/conclusions.png){width=600 fig-align="center"}

- What sort of conclusions are we entitled to make?

## Modelling uncertainty in `creativity` 

- The observed difference (i.e. `treatment_effect`) is not 0. Is there a real difference?

- We need a *(statistical) model* to draw statistical inferences!

## Mental model: the world before randomization

![](http://www.stanford.edu/class/stats191/Chapter1/figs/potential_outcomes.png){width=600 fig-align="center"}

- Potential outcomes before randomization

## Mental model: the null hypothesis

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

- $H_0$: <font color="green">green</font> outcome identical to <font color="red">red</font>


## Mental model: the world after randomization

![](http://www.stanford.edu/class/stats191/Chapter1/figs/potential_outcomes_randomized.png){width=800 fig-align="center"}

- Observed outcomes after randomization


## Computing difference via `t.test`

```{webr-r}
estimates = t.test(Score ~ Treatment, data=creativity)$estimate
estimates
```

### Effect

```{webr-r}
treatment_effect = treatment_effect=estimates[2] - estimates[1]
treatment_effect
```


## Null hypothesis: no difference in `Score` between the groups

```{webr-r}
null_treatment = sample(creativity$Treatment, 47, replace=FALSE)
null_data = data.frame(Score=creativity$Score,
                       Treatment=null_treatment)
null_estimates = t.test(Score ~ Treatment, data=null_data)$estimate
null_estimates
```

```{webr-r}
null_treatment_effect = null_estimates[2] - null_estimates[1]
null_treatment_effect
```


## Repeated 10000 times

```{webr-r}
#| context: output
#| fig-height: 4

# this code/plot is for visualization only
# students not expected to reproduce

estimates = t.test(Score ~ Treatment, data=creativity)$estimate
treatment_effect = treatment_effect=estimates[2] - estimates[1]

null_treatment_effect = rep(NA, 10000)
for (i in 1:10000) {
    null_treatment = sample(creativity$Treatment, 47, replace=FALSE)
    null_data = data.frame(Score=creativity$Score,
                           Treatment=null_treatment)
    null_estimates = t.test(Score ~ Treatment, data=null_data)$estimate
    null_treatment_effect[i] = null_estimates[2] - null_estimates[1]
}

treatment_effect = data.frame(treatment_effect=treatment_effect)
fig <- (ggplot(data.frame(null_treatment_effect),
               aes(x=null_treatment_effect)) +
    geom_histogram(aes(y=after_stat(density)),
                   color="#e9ecef",
                   alpha=0.6, bins=30) +
    geom_vline(aes(xintercept=treatment_effect), treatment_effect,
               color='red', linewidth=2) +
    geom_vline(aes(xintercept=-treatment_effect), treatment_effect,
               color='red', alpha=0.5, linewidth=1) +
    geom_density() +  
    labs(fill=""))
fig
treatment_effect = treatment_effect[1,]
```

```{webr-r}
length(null_treatment_effect)
p_value = mean(abs(null_treatment_effect) > treatment_effect)
p_value
```


## Modelling uncertainty in `salaries` 

- The difference is not 0. Is the difference real?

- We need a *model* to draw statistical inferences!


## Mental model: `Male` and `Female` salaries

![](http://www.stanford.edu/class/stats191/Chapter1/figs/two_sample_salaries.png){width=600 fig-align="center"}

- There are two populations of salaries

## Mental model: `Male` and `Female` salaries

![](http://www.stanford.edu/class/stats191/Chapter1/figs/two_sample_salaries.png){width=600 fig-align="center"}

- $H_0$: distribution of <font color="orange">orange</font> box identical to
<font color="purple">purple</font>

## Computing difference via `t.test`

```{webr-r}
sex_estimates = t.test(Salary ~ Sex, data=salaries)$estimate
sex_effect = sex_estimates[2] - sex_estimates[1]
sex_effect
```

## Repeated 10000 times 

```{webr-r}
#| context: output
#| fig-height: 4

# this code/plot is for visualization only
# students not expected to reproduce

sex_estimates = t.test(Salary ~ Sex, data=salaries)$estimate
sex_effect = sex_estimates[2] - sex_estimates[1]

null_sex_effect = rep(NA, 10000)
for (i in 1:10000) {
    null_sex = sample(salaries$Sex, length(salaries$Sex), replace=FALSE)
    null_data = data.frame(Salary=salaries$Salary,
                           Sex=null_sex)
    null_estimates = t.test(Salary ~ Sex, data=null_data)$estimate
    null_sex_effect[i] = null_estimates[2] - null_estimates[1]
}

sex_effect = data.frame(sex_effect=sex_effect)
fig <- (ggplot(data.frame(null_sex_effect),
               aes(x=null_sex_effect)) +
        geom_histogram(aes(y=after_stat(density)),
                       color="#e9ecef",
	               alpha=0.6, bins=30) +
    geom_vline(aes(xintercept=sex_effect), sex_effect,
               color='red', linewidth=2) +
    geom_vline(aes(xintercept=-sex_effect), sex_effect,
               color='red', alpha=0.5, linewidth=1) +
    geom_density() +  
    labs(fill=""))
fig
sex_effect = sex_effect[1,]
```

```{webr-r}
length(null_sex_effect)
p_value = mean(abs(null_sex_effect) > sex_effect)
p_value
```

## Other issues 

- We used the same method even for these different studies... does this make sense?

- Terminology:

    1. *Parameter*: a property of the probability model (often written $\theta$)

    2. *Estimate*: a function of the sample data (often written $\hat{\theta}$)

    3. Goal of *statistical inference* is to learn about the parameter $\theta$ from the estimate $\hat{\theta}$

## Other issues

- Experimental design:

    1. *Randomization*: individuals were randomly assigned `Treatment` in `creative` study

    2. *Simple random sample*: a way of sampling $n$ from a population such that every $n$ points are equally likely.

    3. Other sampling mechanisms: *systematic sampling*, *cluster sampling*.
