Activity: Power and Type I & II Error
Power and error rates
Worksheet: Power and Type I & II Error
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.
Part 1 · Setup — Load and Explore
▶ Run this in your Script:
library(patchwork)
library(car) # For diagnostic tests
library(tidyverse) # For data manipulation and visualization
library(readxl)
pine_switch_df <- read_excel("data/class_pine needle length switched.xlsx")
head(pine_switch_df)
ps_df <- pine_switch_df %>%
group_by(group, tree_no, tree_char, side) %>%
summarise(length_mm = mean(length_mm, na.rm = TRUE))
ps_shady_df <- ps_df %>% filter(side == "shady")
ps_sunny_df <- ps_df %>% filter(side == "sunny")
head(ps_df)✏️ Your turn: Why do we average by group, tree_no, tree_char, and side before doing anything else? ________________________
Part 2 · Hypothesis Testing Recap
Before we test anything, review the framework from lecture: H₀, Hₐ, and α are all decided before you look at the data.
✏️ Your turn: For the question “Do needle lengths differ between the sunny and shady sides?”, write:
H0 = ________________________________________________________
Ha = ________________________________________________________
alpha = ________✏️ Your turn: In one sentence, what’s the difference between a Type I error and a Type II error?
________________________________________________________🚀 If you finish early: With df = 14 and a two-tailed α = 0.05, find the critical t-value with qt(0.975, df = 14). Then use pt() to find the p-value for an observed t of 2.8 on the same df: 2 * (1 - pt(2.8, df = 14)). Is 2.8 inside or outside the rejection region?
# Write your code here:Part 3 · Exploratory Data Analysis and Effect Size
▶ Run this:
pine_summary <- ps_df %>%
group_by(side) %>%
summarize(
mean_length = mean(length_mm, na.rm = TRUE),
sd_length = sd(length_mm, na.rm = TRUE),
n = sum(!is.na(length_mm)),
se_length = sd_length / (n^0.5),
t_critical = qt(0.975, df = n - 1), # 95% CI uses 0.975 (two-tailed)
ci_lower = mean_length - t_critical * se_length,
ci_upper = mean_length + t_critical * se_length
)
print(pine_summary)▶ Run this — visualize the distribution:
ggplot(ps_df, aes(x = length_mm, fill = side)) +
geom_histogram(binwidth = 2) +
labs(title = "Distribution of Pine Needle Lengths",
x = "Length (mm)", y = "Frequency") +
theme_minimal() +
facet_wrap("side", ncol = 1)▶ Run this — the effect size (difference in means):
pine_summary %>%
summarize(difference = mean_length[side == "sunny"] - mean_length[side == "shady"])✏️ Your turn: Looking at the histogram, do the two distributions look like they overlap a lot, or are they clearly separated? Does that match the size of the difference you just calculated? ________________________
🚀 If you finish early: Recompute pine_summary, but group by group (field team) instead of side. Do all teams show a similar mean and spread, or does one team look different?
# Write your code here:Part 4 · The Two-Sample t-Test
H₀: μ₁ = μ₂ (the mean needle lengths are equal)
Hₐ: μ₁ ≠ μ₂ (the mean needle lengths are different)
Deciding between the standard t-test (equal variances) and Welch’s t-test (unequal variances) depends on a Levene’s test result — review the T-Tests worksheet if you need a refresher.
▶ Run this:
t_test_result <- t.test(length_mm ~ side, data = ps_df, var.equal = TRUE)
print("Standard two-sample t-test:")
print(t_test_result)✏️ Your turn: What is the t-statistic, degrees of freedom, and p-value? Do you reject or fail to reject H₀ at α = 0.05?
t = ________ df = ________ p = ________
Decision: ________________________🚀 If you finish early: Rerun the test with var.equal = FALSE (Welch’s t-test) and compare the p-value to the standard version above. How different are they?
# Write your code here:Part 5 · Statistical Power
Statistical power is the probability of detecting a true effect — rejecting H₀ when it is actually false. A power analysis is typically done to determine the required sample size before a study, to evaluate whether a completed study’s sample size was adequate, or to find the minimum detectable effect size for a given sample. 80% power is generally considered acceptable.
\[s_p = \sqrt{\frac{(n_1 - 1)s_1^2 + (n_2 - 1)s_2^2}{n_1 + n_2 - 2}}\]
▶ Run this — power to detect a 1mm difference:
side_diff <- 1
sunny_n <- nrow(ps_sunny_df)
shady_n <- nrow(ps_shady_df)
sun_sd_pooled <- sqrt((var(ps_sunny_df$length_mm) * (sunny_n - 1) +
var(ps_shady_df$length_mm) * (shady_n - 1)) /
(sunny_n + shady_n - 2))
sun_effect_size <- side_diff / sun_sd_pooled
sun_df <- sunny_n + shady_n - 2
sun_alpha <- 0.05
sun_power <- power.t.test(n = min(sunny_n, shady_n),
delta = side_diff, # Raw difference, not effect size
sd = sun_sd_pooled, # Use pooled SD
sig.level = sun_alpha,
type = "two.sample",
alternative = "two.sided")
print("Sample sizes:")
print(paste("Sunny:", sunny_n, "Shady:", shady_n))
print(paste("Pooled SD:", round(sun_sd_pooled, 3)))
print(paste("Effect size (Cohen's d):", round(sun_effect_size, 3)))
print("Power analysis results:")
print(sun_power)✏️ Your turn: Is the power above 80%? If our study is underpowered, what are two things we could change to increase power? ________________________
🚀 If you finish early: power.t.test() can also solve for sample size instead of power — leave n out and add power = 0.80 instead. How many needles per side would you need to reliably (80% power) detect a 1mm difference?
power.t.test(delta = side_diff,
sd = sun_sd_pooled,
sig.level = sun_alpha,
power = 0.80,
type = "two.sample",
alternative = "two.sided")Part 6 · A Publication-Quality Plot
Typically we present a plot with the mean and standard error to represent the data.
▶ Run this:
pine_mean_se <- ps_df %>%
ggplot(aes(side, length_mm, color = side)) +
stat_summary(fun = "mean", na.rm = TRUE, geom = "point", size = 3) +
stat_summary(fun.data = "mean_se", width = 0.2, geom = "errorbar")
pine_mean_se▶ Run this — a custom theme you can reuse on any plot:
theme_class <- function(base_size = 14, base_family = "Sans") {
theme(
panel.background = element_rect(fill = "transparent", colour = "transparent"),
plot.background = element_rect(fill = "transparent", colour = NA),
panel.grid.major = element_line(linetype = "blank"),
panel.grid.minor = element_line(linetype = "blank"),
axis.text = element_text(colour = "black"),
axis.title.x = element_text(size = 18, face = "bold"),
axis.title.y = element_text(size = 18, face = "bold"),
axis.text.x = element_text(size = 16, face = "bold", angle = 0, vjust = .5, hjust = .5),
axis.text.y = element_text(size = 16, face = "bold"),
axis.ticks = element_line(colour = "black"),
axis.line.x = element_line(color = "black", linewidth = 0.5, linetype = "solid"),
axis.line.y = element_line(color = "black", linewidth = 0.5, linetype = "solid"),
legend.text = element_text(colour = "black", size = 16, face = "bold"),
legend.title = element_text(colour = "black", size = 18, face = "bold"),
legend.position = "right",
legend.key = element_rect(fill = "transparent", colour = "transparent"),
legend.background = element_rect(fill = "transparent", colour = "transparent"),
panel.border = element_blank(),
plot.title = element_text(hjust = 0, vjust = 2.12),
plot.caption = element_text(hjust = 0, vjust = 1.12)
)
}⚠️ Watch out! Run the whole
theme_class <- function(...) {...}block at once (select it all and run) — running just the inside line by line will fail, since those lines only make sense as part of the function.
▶ Run this — apply your theme to the plot (note +, not <-, so the original pine_mean_se object is untouched):
pine_mean_se +
theme_class()▶ Run this — put it all together with proper labels, a fixed y-axis range, and custom colors:
pine_mean_se <- ps_df %>%
ggplot(aes(side, length_mm, color = side)) +
stat_summary(fun = "mean", na.rm = TRUE, geom = "point", size = 4) +
stat_summary(fun.data = "mean_se", geom = "errorbar", width = 0.1, linewidth = 0.35) +
labs(x = "Side of Tree", y = "Mean Length (mm ± 1 SE)") +
coord_cartesian(ylim = c(15, 20)) +
theme_class() +
scale_color_manual(
name = "Side of tree",
labels = c("shady" = "Shady Side", "sunny" = "Sunny Side"),
values = c("shady" = "darkgreen", "sunny" = "coral")
)
pine_mean_se✏️ Your turn: Change the two hex-ish color names in values = to two colors of your own choosing (e.g., "darkblue", "darkorange"). Re-run and confirm the plot updates.
# Write your code here:🚀 If you finish early: pine_mean_se currently uses color =, which only affects points and error bars — not a filled shape. Build a boxplot version using aes(side, length_mm, fill = side) + geom_boxplot() + scale_fill_manual() with the same labels/values pattern.
# Write your code here:Summary and Conclusions
In this activity, we’ve:
- Formulated hypotheses about pine needle length
- Reviewed Type I and Type II error and statistical power
- Conducted a two-sample t-test
- Calculated statistical power for our observed effect size
- Built a publication-quality mean ± SE plot with a custom theme and custom colors
Key takeaways:
- Always check assumptions before conducting tests
- Visualize your data to understand patterns
- Report results comprehensively, including effect size and power when relevant
- Consider alternatives when assumptions are violated
Reflection Questions
- How does sample size affect our confidence in estimating the population mean?
- Why is the t-distribution more appropriate than the normal distribution when working with small samples?
- When comparing two populations, what can we learn from confidence intervals that a t-test alone doesn’t tell us?
- How would you explain the concept of statistical power to someone who has never taken a statistics course?
- What do we do if assumptions fail?
Review and checkpoint
At this point you can:
📤 What to turn in before next class
Upload both of these to the course management system:
- Your code — the
scripts/folder (or just07_power_and_error.R) - This worksheet, with your written answers
Getting unstuck
When code breaks — and it will, that is normal:
- Read the error message out loud. R usually names the line and the problem.
- Check the usual suspects: did you run
library(tidyverse),library(car), andlibrary(readxl)? Spelling? A missing)or%>%at the start of a line? ?function_nameopens the built-in help page.- 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.