[![](../rcode_button.png){fig-alt="R Code Button" fig-align="left" width=20%}](lab07.Rmd)

The goal of this lab is to introduce the RMarkdown file format.
This format is convenient for including R code and output into an html 
or pdf document. 
By default the Rmd file you downloaded produces an html output, because that
is the format you will find the lab in on the external course website. 
We can change the output type to `pdf_document` if we want to create a pdf. 

```{r metadata, include=FALSE}
# Author: Jarad Niemi
# Date:   2024-10-17
# Purpose: Introduce RMarkdown file format using Module 2 content
#-------------------------------------------------------------------------------

```

We will need the following packages for this file. 

```{r packages}
library("tidyverse")
library("Sleuth3")   # data
library("rmarkdown") # compiling RMarkdown files
```

If you want to compile to a pdf, you will need a 
LaTeX installation. 
On a Mac, I use [MacTex](https://www.tug.org/mactex/).
On a PC, I have used [MikTex](https://miktex.org/).
On either, you should be able to use the R package
[tinytex](https://yihui.org/tinytex/) which is available on CRAN
and following the [instructions](https://yihui.org/tinytex/).


## Statistical analyses

### One binomial population

A binomial random variable is one where you have the **count of the number of
successes** out of some number of attempts where each attempt is 
**independent** and has the same **probability of success**.

Many datasets can be converted to a binomial random variable. 
For example, consider the `case0401` data set in the 
[Sleuth3]() R package (make sure you have installed and loaded it). 
You can find information about this data set using `?case0401`.

```{r case0401}
case0401
```

Let $y$ the number of non-incident Warm launches and $n$ be the total number
of warm launches. 
Assume $Y \sim Bin(n,\theta)$. 

```{r one-binomial-data}
n <- sum(case0401$Launch == "Warm")
y <- sum(case0401$Incidents[case0401$Launch == "Warm"] == 0)
```

We can use some built-in functions to calculate quantities for us. 

```{r one-binomial-frequentist}
## Frequentist analysis
binom.test(y, n)
prop.test( y, n) # based on CLT (not reasonable here)
```

On homework, rather than printing out the entire analysis, 
you should print out relevant quantities for the questions asked.
For example, if question 1.5 requested an approximate 95% confidence interval,
you could use the following

```{r one-binomial-confidence-interval}
binom.test(y, n)$conf.int # Q1.5
```

We can also do more manual calculations, like those for a Bayesian analysis. 

```{r one-binomial-bayesian}
# Bayesian analysis
(1 + y) / (2 + n)
qbeta(c(.025, .975), 1 + y, 1 + n - y) # Credible interval
1 - pbeta(0.9, 1 + y, 1 + n - y)       # P(theta > 0.9 | y)
```

Make sure your code and text help the reader understand what you are doing. 

### One normal population

We use a normal model when our data are continuous. 
Consider the data in `case0101`. 

Let $y_i$ be the creativity score for the $i$th individual in the Extrinsic
group. 
Assume $Y_i \stackrel{ind}{\sim} N(\mu,\sigma^2)$. 

```{r case0101-summary}
extrinsic <- case0101 |> 
  filter(Treatment == "Extrinsic") 

summary(extrinsic)
```

Frequentist and Bayesian analyses are both accomplished through the `t.test`
function.

```{r one-normal-t-test}
t.test(extrinsic$Score)
```

The default hypothesized value for the mean is $0$. 
Since the possible range for creativity scores is likely 0 to 30,
this hypothesized value is completely unreasonable. 
Therefore, the p-value in this output is meaningless. 


### Two binomial populations

Analyses get more interesting when we compare two different populations. 
Using the shuttle launch data (`case0401`), 
we might be interested in comparing the probability of no O-ring issues depending
on temperature which is recorded as Warm or Cool.
Thus, we have two populations (Warm and Cool) and within each we have the 
number of successes out of some number of attempts. 

Let $Y_g$ be the number of non-incident launches out of $n_g$ attempts in 
group $g$. Define $g=1$ to be Cool and $g=2$ to be Warm. 
Assume $Y_g \stackrel{ind}{\sim} Bin(n_g, \theta_g)$. 
We are interested in the quantity $\theta_1 - \theta_2$ (or the reverse). 

```{r two-binomial-data}
# Compute summary statistics from data
d <- case0401 |>
  group_by(Launch) |>
  summarize(
    n = n(),
    y = sum(Incidents == 0),
    p = y/n
  )

d
```

Unfortunately, frequentist analyses generally require data asymptotics,
i.e. large sample sizes. 

```{r two-binomial-frequentist}
prop.test(d$y, d$n)
```
Fortunately, `prop.test` warns us (with a cryptic message) that the sample size
is not large enough. 

We can perform a Bayesian analysis with any sample size,
but our prior will have more impact on the analysis when the sample size is 
small.

```{r two-binomial-bayesian}
# Simulate from the posterior for each probability of success
nreps      <- 1e5                                           # make this large
prob_cool <- rbeta(nreps, 1 + d$y[1], 1 + d$n[1] - d$y[1]) 
prob_warm <- rbeta(nreps, 1 + d$y[2], 1 + d$n[2] - d$y[2])
diff       <- prob_cool - prob_warm

# Calculate quantities of interest
mean(diff < 0)                          # P(prob_warm > prob_cool | y)
quantile(-diff, probs = c(.025, 0.975)) # 95% credible interval for prob_warm - prob_cool
```

This interval is pretty wide for the difference between two probabilities,
but with only 24 binomial observations, what did you expect?


### Two normal populations

When we have continuous data that arises from two populations, 
we can often use two normal distributions to compare them. 

Consider the `case0101` data with creativity scores that may depend on treatment
an individual received. 
We can compute summary statistics. 

```{r two-normal-data}
# Summary statistics
case0101 |>
  group_by(Treatment) |>
  summarize(
    n    = n(),
    mean = mean(Score),
    sd   = sd(Score)
  )
```
We can also plot the data

```{r two-normal-plot}
ggplot(case0101,
       aes(x = Treatment, 
           y = Score,
           color = Treatment,
           shape = Treatment)) +
  geom_jitter(width = 0.1)       # jitter the points makes sure points don't overlap
```

A common assumption is whether the two populations
have the same population variance.
We can check this by looking at the sample standard deviations as well as the
plots above. 
here the variability looks pretty similar between the two groups. 

We can perform both frequentist and Bayesian analyses using `t.test`.

```{r two-normal-analyses}
(tt <- t.test(Score ~ Treatment, data = case0101, var.equal = TRUE))
```

Here hypothesized difference between the means of the two groups is zero and 
thus the p-value is relevant. 
Before we are rash in our conclusions, 
we should check and make sure independence and normality are reasonable. 
 




## Inference

For the above two populations, 
consider whether the data are a random sample and whether the treatment has
been randomly applied. 



### Random sample

The shuttle launch data is probably the entire population over some time 
period. 
With that being said, we can consider that it is a sample of shuttle launches
that could have taken place. 
But it is unlikely that this sample is random. 
Ultimately, for these data the statistics provide us a reasonable summary of 
the data, but cannot be used to infer something about the larger population. 

Psychology experiments like `case0101` are typically performed with an 
convenience sample of undergraduate students at the institution where the 
researchers are. 
Thus, the sample is not random and certainly not representative of the 
population at large.
Thus, we cannot infer results about some larger population. 

### Randomized treatment

The shuttle launch data collection does not have a randomized treatment. 
That is, we did not randomly assign the temperature to a given day since
that is impossible. 
Thus we may only say there is an association (a strong one here) between 
temperature and probability of no O-ring incidents on a given launch. 
Even though the association is strong, that does not (on its own) imply there
is a causal relationship between temperature and probability of no incident.

```{r my-round, echo=FALSE}
my_round <- function(x, nd = 2) {
  format(round(x, nd), nsmall = nd)
}
```

As indicated in the helpfile for `case0101`, 
individuals were randomized to the two treatment groups. 
Since treatment was randomized, we can make a causal claim about the effect 
of motivation type (extrinsic vs intrinsic) on creativity score. 
In this case, intrinsic motivation score appears to increase creativity scores
on average by (`r my_round(-tt$conf.int[2],1)`,`r my_round(-tt$conf.int[1],1)`)
on a 30 point scale.
