---
title: "Quack Sampling"
author: "Sophie Regan"
date: "2025-11-03"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
```

################### Warm up #########################
```{r}
# 1. Write a for loop that makes a string indicating the output file paths for
# each of the following participants and saves it them in a vector called filepaths:
subj <- c("subj01", "subj02", "subj03", "subj04", "subj05", "subj06", "subj07", "subj08")

# The file path for a participants should also include the condition that they were in:
cond1 <- c("subj01", "subj04", "subj07", "subj08")
cond2 <- c("subj02", "subj03", "subj05", "subj06")

# For example, the output for participant 1 would be: "../data/cond1/data_subj01"
# and for participant 2 it would be: "../data/cond2/data_subj02"



# 2. Look at the function rnorm() in the Help window (you can use ?rnorm()) or google it. What does rnorm do?
```



```{r}
#### We can simulate the data collection process in R ####

# Define our random variable with a function



# Plot our observations


# The exact values change but the definition of heights is consistent. 

# In practice we don't want the values to change each time. To keep our results consistent, we want to simulate random data collection ONCE.

set.seed(342)


# This lets us reproduce the same sequence of random numbers every time we run the analysis. 
# Try changing the seed and see what happens.
```

```{r}
#### Estimating means of a population ####

# Imagine that you live on a planet called Quackland and are baby researchers
# studying gyrification of different regions of the brain. We want to know the
# average gyrification of a certain region of interest (r1) in 1-year-old
# children.

# How would we go about estimating the average gyrification of r1 in
# 1-year-olds?








# Now imagine that Willa and Elena are magical wizards who happen to know the
# actual gyrification of r1 in all 40,000 1-year-olds who live on Quackland.
# Lucky for you, we posted it on OSF just for us to use in this class!

# Load in our #openscience data!
# We are interested in r1
# Data adapted by Elena from a CSV file used by Daniel Catterson in Psych 101
baby <- read.csv("../data/brainwavebabydata.csv")



# Let's randomly sample from our population (as we described above) and estimate
# our mean.




# What if we sample just 10 babies? Plot a histogram. What is the mean
# gyrification? The sd?



# As all knowing wizards, Elena and Willa have granted you the power to check
# the TRUE MEAN gyrification of the population. Check that value. How does it
# compare to the estimated mean from the 10-baby sample?



# Now, what if we took a bunch of samples (5) of 10 babies and took the mean of
# the means (e.g., 5 follow-up studies). How does this value compare to the TRUE
# population mean? Plot a histogram of the means. What is the mean gyrification
# of the means? The sd of the means?




# What if we could take 100 samples of 10 babies? How does this value compare to
# the TRUE population mean? Plot a histogram of the means. What is the mean
# gyrification of the means? The sd of the means?
```

```{r}
## Practice: Now try this all again, but now drawing samples of 100 babies. ##

# Draw one sample of 100 babies. How does this compare to the mean r1 depth from the 1 sample of 10 babies? To the true population mean?


# Draw 100 babies 5 times and compare the means


# Draw 100 babies 100 times and compare the means
```

```{r}
#### The Central Limit Theorem! ####



#### How reliable is our measure? ####

# More realistically, we we can only collect one sample, or maybe up to a few.
# Take our penguins for example.
penguins <- read.csv("../data/penguins_clean.csv") # Only 2008 sample!


# Suppose we want to know how reliable is our measure of interest (e.g., mean
# body weight, relation between body weight and bill length, etc.). 

#     What is the *ideal* way to test this?


#     What is a possible work-around?


```


```{r}
# We could take OUR SAMPLE and pretend it is the POPULATION and then sample *from
# it*!!!


# How reliable is the mean body mass?
mean(penguins$body_mass_g)
sd(penguins$body_mass_g)

# In other words, how much would it vary across samples?

# Ideally we would collect many many samples of penguins to do this. But... We can't. So the next best thing is to use our own data as if they were the whole population, and *sample from it*!!!

# Key here is: SAMPLE WITH REPLACEMENT!


# Let's sample 1000 times from our data and take the mean.




# Plot our means



# What is the mean and sd of our means?




# This is called "Bootstrapping"!!! Sampling with replacement + computing some
# metric of interest.



#### Is the difference between two groups meaningful and reliable? ####

mean(penguins[penguins$sex == "female",]$body_mass_g)
mean(penguins[penguins$sex == "male",]$body_mass_g)


# What are some ways that we could test whether this difference is reliable?

```


```{r}
######################### How accurate is our measure? #########################

# In the real world, we can't know the true population mean (or any statistic).
# We also can't draw an infinite number of samples from the population in order
# to estimate the population mean!

# The best we can do is take one sample, may a few.

# Take our penguins data for example. It was very costly to collect data on so
# many penguins!

# What is the mean body mass of all the penguins?
penguins <- read.csv("../data/penguins_clean.csv") # (Note: only 2008 sample)




# Suppose we want to know how accurate our measure of interest is (e.g., mean
# body weight, relation between body weight and bill length, etc.). What is the
# range of values the mean can take on?

# What is the *ideal* way to test this since we can't measure the whole
# population?
#     The ideal way to test this would be to draw many samples from the
#     population, get the mean body_mass, and estimate the population mean from
#     the mean of the sample means!


# What is a possible work-around?
#     Since it is not realistic to draw so many real samples from the
#     population, we can use our current sample, pretend that it is the
#     population, and sample from it! The key here is that we resample from our
#     sample WITH REPLACEMENT! We assume that if we sampled radomly from the
#     population then our sample is representative of the population! Including
#     the frequencies with which different values of body mass come up. This
#     allows us to let our sample act as an approximation of our population in
#     order to estimate how much variation there is in our measure of interest
#     from sample to sample.



# Key insight: We can take OUR SAMPLE and pretend it is the POPULATION and then
# sample *from it*!!!


# How accurate is the mean body mass? In other words, how much would it vary across samples?
mean(penguins$body_mass_g)
sd(penguins$body_mass_g)

# Let's use resampling to estimate how much the mean body mass varies between samples

# Goal: Resample 1000 times *from our data* and take the mean every time.

# What do we need to do? (Pseudocode)
# 1) 
# 2) 
# 3)

# Code




# Plot our means


# What is the mean and sd of our means?




# What does the standard deviation tell you?




# When we report any statistic, we need to always report error!
# For example, our mean body mass is  +/-  grams. 
# Or better yet, report the 95% confidence interval:
# Lower bound:


# Upper bound: 


print(paste("Based on our sample, we estimate that the mean body mass of the population of penguins is between", round(lb, 3), "and", round(ub, 3), "grams."))

# OR
print(paste0("We found that penguin body mass is ", round(mean(penguins$body_mass_g), 4), " grams (95% CI: [", round(lb, 3), ", ", round(ub, 3), "])."))

# 95% CI means that we estimate that 95% of random samples we draw from our
# population will have mean body masses between the lower bound and upper bound
# values.



# This resampling technique is called "Bootstrapping"!!! Sampling with
# replacement from our collected data to find the amount of sampling error for
# some statistic of interest.

```

```{r}

############### Is the difference between two groups meaningful? ###############
# Find the body mass means for each group
(f.mean <- with(penguins, mean(penguins[sex == "female",]$body_mass_g)))
(m.mean <- with(penguins, mean(penguins[sex == "male",]$body_mass_g)))

# Do it with tidyverse, too!



# Plot the body mass for each group to visualize this difference



# Let's quantify this difference in mean body mass between the group

# We want to know: is this difference different than 0? In other words, is there
# really a group difference? Or is this difference that we found just due to
# chance?


# What are some ways that we could test whether this difference is reliable?
# 






# Here's another idea:
# We randomly shuffle the labels of sex and see if the two groups are different or not. If the label is MEANINGFUL then there will be a difference in the original data but there will not be a difference in the shuffled data!





# Do this many times and see what happens!




# Plot all the differences in means



# What percent of our shuffled samples had differences in means more extreme than our difference in means?





# This is called a permutation test!
# Shuffle our labels WITHOUT REPLACEMENT! Calculate our test statistic, and do
# it many times. Test whether the test statistic that we got is more extreme
# than the distribution of test statistics from data sets with randomly shuffled
# labels.
```

