Activity: Study Design & Sampling

Design and sampling

Hands-on activity: natural experiments, pseudoreplication, and a priori/post hoc power analysis on a simulated waterfall fish dataset.

Worksheet: Study Design and Power Analysis

How to use this worksheet

Work through each part in order, at your own pace. Type every line of code yourself into a plain R script — do not copy-paste. Blocks marked ▶ Run this are code you should type and execute. Blocks marked ✏️ Your turn ask you to write, modify, or answer something. Boxes marked 🚀 If you finish early are optional bonus material that goes a bit further than what we covered in lecture.

This worksheet uses a simulated fish dataset — we build the data ourselves in R, rather than loading a file — so you can see exactly what assumptions go into it. We’ll study fish mass above and below waterfalls across multiple streams, and use it to practice recognizing proper replication vs. pseudoreplication and running power analysis before and after data collection.

---

Part 1 · Setup and the Research Question

▶ Run this in your Script:

library(tidyverse)
library(pwr)

set.seed(42)

Main Question: Do waterfalls affect fish mass in stream ecosystems?

Hypothesis: Fish below waterfalls are larger, due to better feeding opportunities from nutrients and prey coming down from the lake.

✏️ Your turn: Write this as a formal null and alternative hypothesis:

H0 = ________________________________________________________
Ha = ________________________________________________________

---

Part 2 · Study Design — a Natural Experiment

We’ll sample fish above and below waterfalls in several streams. This is a natural experiment, since we can’t manipulate whether a waterfall is present.

▶ Run this — simulate the data:

# ADJUST THESE VARIABLES TO EXPLORE DIFFERENT SCENARIOS:
n_streams <- 5         # Number of streams to sample
n_fish_per_site <- 6   # Number of fish per site (above/below)
waterfall_effect <- 15 # Effect size: how much larger fish are below (in grams)

stream_df <- tibble(
  stream_id = 1:n_streams,
  stream_mean = rnorm(n_streams, mean = 85, sd = 5)  # Random baseline for each stream
)

f_df <- tibble(
  stream_id = rep(1:n_streams, each = n_fish_per_site * 2),
  position = rep(c("above", "below"), each = n_fish_per_site, times = n_streams),
  fish_id = 1:(n_streams * n_fish_per_site * 2)
) %>%
  left_join(stream_df, by = "stream_id") %>%
  mutate(
    expected_mass = ifelse(position == "above", stream_mean, stream_mean + waterfall_effect),
    mass_g = rnorm(n(), mean = expected_mass, sd = 12)
  ) %>%
  select(stream_id, position, fish_id, mass_g)

head(f_df, 12)

▶ Run this — visualize the pooled data, then by stream:

basic_plot <- f_df %>%
  ggplot(aes(x = position, y = mass_g, fill = position)) +
  geom_boxplot() +
  geom_jitter(width = 0.2, alpha = 0.6) +
  labs(x = "Position Relative to Waterfall", y = "Fish Mass (g)") +
  theme_minimal()
basic_plot

stream_plot <- f_df %>%
  ggplot(aes(x = position, y = mass_g, group = stream_id, color = as.factor(stream_id))) +
  stat_summary(fun = mean, geom = "point") +
  stat_summary(fun = mean, geom = "line") +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1) +
  facet_wrap(~stream_id) +
  labs(x = "Position", y = "Fish Mass (g)", color = "Stream") +
  theme_minimal()
stream_plot

✏️ Your turn: Looking at stream_plot, does every stream show fish getting larger below the waterfall, or does the pattern vary stream to stream? ________________________

Tip

🚀 If you finish early: Change waterfall_effect to 0 and rerun the whole simulation. Does basic_plot still look like there’s a difference? This is what “no true effect” looks like in this simulated world.

# Write your code here:

---

Part 3 · Replication Issues

Question: In our fish study, what is the true experimental unit?

1. Individual fish
2. Stream locations (above/below pairs)
3. Individual streams
4. All fish combined

✏️ Your turn: Circle your answer above, then explain why in one sentence.

________________________________________________________

▶ Run this — the WRONG analysis, treating every fish as independent:

# This ignores that fish within a stream may be similar
pseudo_test <- t.test(mass_g ~ position, data = f_df)
pseudo_test

▶ Run this — the CORRECT analysis, averaging by stream first:

stream_means <- f_df %>%
  group_by(stream_id, position) %>%
  summarize(mean_mass = mean(mass_g), .groups = "drop")

stream_wide <- stream_means %>%
  pivot_wider(names_from = position, values_from = mean_mass)

proper_test <- t.test(stream_wide$below, stream_wide$above, paired = TRUE)
proper_test

The proper analysis uses streams as replicates (n = 5, if you kept the default n_streams), not individual fish (n = 36 per group).

✏️ Your turn: Compare the p-values from pseudo_test and proper_test. Which one is smaller? Why does pretending fish are independent make the result look more significant than it should? ________________________

Tip

🚀 If you finish early: Rerun Part 2’s simulation with n_fish_per_site <- 50 instead of 6 (keep everything else the same). Does pseudo_test’s p-value get even smaller, even though the true number of streams (and true replicates) hasn’t changed?

# Write your code here:

---

Part 4 · Power Analysis — Planning Phase

Cohen’s d measures the difference between two group means in standard deviations. It helps understand the magnitude of a difference beyond statistical significance. A Cohen’s d of 0.2 is a small effect, 0.5 a moderate effect, 0.8 a large effect.

▶ Run this — set expected values from pilot data, and compute the effect size:

above_mean <- 85
below_mean <- 110
pooled_sd <- 20

effect_val <- abs(below_mean - above_mean) / pooled_sd
effect_val

▶ Run this — required sample size for 80% power:

power_result <- pwr.t.test(
  d = effect_val,
  sig.level = 0.05,
  power = 0.8,
  type = "paired"
)
power_result

✏️ Your turn: According to power_result, how many streams do we need for 80% power? ________________________

▶ Run this — visualize the power curve:

power_curve_df <- tibble(streams = 3:15) %>%
  rowwise() %>%
  mutate(power = pwr.t.test(n = streams,
                            d = effect_val,
                            sig.level = 0.05,
                            type = "paired")$power) %>%
  ungroup()

curve_plot <- ggplot(power_curve_df, aes(x = streams, y = power)) +
  geom_line(linewidth = 1.2, color = "blue") +
  geom_hline(yintercept = 0.8, linetype = "dashed", color = "red") +
  geom_vline(xintercept = ceiling(power_result$n), linetype = "dashed", color = "red") +
  labs(x = "Number of Streams", y = "Statistical Power") +
  theme_minimal()
curve_plot

✏️ Your turn — experiment with different scenarios. Change the values below and rerun:

above_mean <- 85
below_mean <- 100    # Try changing this
pooled_sd <- 25       # Try changing this

effect_val <- abs(below_mean - above_mean) / pooled_sd

power_result <- pwr.t.test(
  d = effect_val,
  sig.level = 0.05,
  power = 0.8,
  type = "paired"
)
power_result

1. What happens to the required sample size if the effect is smaller (below_mean = 95)?
2. What if variation is higher (pooled_sd = 30)?
3. What if we want 90% power instead of 80%?

________________________________________________________

Tip

🚀 If you finish early: Re-derive power_result using type = "two.sample" instead of "paired", keeping the same effect size. Does the paired or unpaired design need fewer streams for 80% power? Does that match what lecture said about paired designs and power?

# Write your code here:

---

Part 5 · Post-Hoc Power Analysis

Now let’s analyze the power we actually had with our streams from Part 2.

▶ Run this:

observed_above <- mean(stream_wide$above)
observed_below <- mean(stream_wide$below)
observed_diff <- observed_below - observed_above
observed_sd <- sd(stream_wide$below - stream_wide$above)

observed_effect <- observed_diff / observed_sd
observed_effect

actual_power <- pwr.t.test(
  n = 6,
  d = observed_effect,
  sig.level = 0.05,
  type = "paired"
)
actual_power

✏️ Your turn: What was our actual power? Is it above or below the 80% “acceptable” threshold from lecture? ________________________

Tip

🚀 If you finish early: actual_power used n = 6 hardcoded. Replace it with n = n_streams so it automatically matches whatever you set in Part 2. Rerun Part 2 with n_streams <- 10 and confirm actual_power updates accordingly.

# Write your code here:

---

Part 6 · Alternative Analysis Approaches

▶ Run this — what if we ignored the stream pairing?

unpaired_test <- t.test(mean_mass ~ position, data = stream_means)
unpaired_test

▶ Run this — summary statistics at both levels:

# By individual fish (the WRONG level of replication)
summary_df <- f_df %>%
  group_by(position) %>%
  summarize(
    n_fish = n(),
    mean_mass = mean(mass_g),
    sd_mass = sd(mass_g),
    se_mass = sd_mass / sqrt(n_fish)
  )
summary_df

# By stream (the proper level of replication)
stream_summary <- stream_means %>%
  group_by(position) %>%
  summarize(
    n_streams = n(),
    mean_mass = mean(mean_mass),
    sd_mass = sd(mean_mass),
    se_mass = sd_mass / sqrt(n_streams)
  )
stream_summary

✏️ Your turn: Compare the se_mass values in summary_df (fish-level) vs. stream_summary (stream-level). Which is smaller, and why does that make fish-level analysis look artificially more precise? ________________________

Tip

🚀 If you finish early: Compare the p-value of unpaired_test (Part 6) to proper_test (Part 3, which was paired). Which is more powerful for this same data — the paired or unpaired analysis at the stream level?

# Write your code here:

---

Part 7 · Design Your Own Study

Scenario: You want to study if fish size differs between fast and slow water areas.

Design questions:

1. What is your experimental unit?

________________________________________________________

2. How many replicates do you need? Use these values: fast water mean = 75g, slow water mean = 85g, standard deviation = 15g, desired power = 80%.

fast_mean <- 75
slow_mean <- 85
sd_val <- 15

effect_size <- abs(slow_mean - fast_mean) / sd_val

sample_result <- pwr.t.test(
  d = effect_size,
  sig.level = 0.05,
  power = 0.8,
  type = "two.sample"  # or "paired" if appropriate
)
sample_result

3. What could cause pseudoreplication in this design?

________________________________________________________

4. How would you avoid it?

________________________________________________________

Tip

🚀 If you finish early: If fast and slow water sites can be paired within the same stream (like above/below waterfall was), redo the calculation with type = "paired". How many fewer replicates do you need?

# Write your code here:

---

Summary

Key points:

• Experimental unit — the thing that receives the treatment independently
• Replication — must be at the level of the experimental unit
• Power analysis — plan sample size before collecting data
• Paired vs. unpaired — paired tests are more powerful when appropriate
• Effect size — larger effects need fewer samples to detect

Common mistakes:

• Treating subsamples as independent replicates
• Analyzing data at the wrong level
• Collecting data before planning the analysis
• Ignoring natural pairing in the design

Review and checkpoint

At this point you can:

• Identify the true experimental unit in a nested/repeated-measures design
• Recognize pseudoreplication and fix it by aggregating to the correct level
• Run an a priori power analysis to plan sample size with pwr.t.test()
• Run a post hoc power analysis on completed data
• Explain why a paired design is often more powerful than an unpaired one

Note

📤 What to turn in before next class

Upload both of these to the course management system:

1. Your code — the scripts/ folder (or just 09_study_design_sampling.R)
2. This worksheet, with your written answers

---

Getting unstuck

When code breaks — and it will, that is normal:

1. Read the error message out loud. R usually names the line and the problem.
2. Check the usual suspects: did you run library(tidyverse) and library(pwr)? Spelling? A missing ) or %>% at the start of a line?
3. ?function_name opens the built-in help page.
4. Bring the exact error (copy-paste it) to class or office hours.

💡 Key idea: Every working scientist googles error messages daily. Getting stuck is not failing — it is the job.

---