library(tidyverse)
theme_set(theme_classic())
# simulate two populations A and B
sigma <- 1.5
mu_df <- data.frame(Group=c("A", "B"),
Lab=c("mu[A]", "mu[B]"),
mu=c(12, 17)) |>
mutate(density=dnorm(mu, mu, sigma))
pop_df <- data.frame(value=seq(min(mu_df$mu) - 3*sigma,
max(mu_df$mu) + 3*sigma,
length.out=1e3)) |>
mutate(A=dnorm(value, mu_df$mu[1], sigma),
B=dnorm(value, mu_df$mu[2], sigma)) |>
pivot_longer(2:3, names_to="Group", values_to="density")
x_lim <- xlim(min(pop_df$value), max(pop_df$value))
y_lim <- ylim(-max(pop_df$density)*0.1, max(pop_df$density)*1.1)
# simulate sampling
N <- 10
set.seed(1)
samp_df <- data.frame(Group=rep(c("A", "B"), each=N),
y=c(rnorm(N, mu_df$mu[1], sigma),
rnorm(N, mu_df$mu[2], sigma))) |>
mutate(ypos=seq(max(mu_df$density)*0.05, max(mu_df$density)*0.6, length.out=n())) |>
group_by(Group) |>
mutate(y_bar=mean(y),
y_bar_pos=mean(ypos)) |>
ungroup() |>
mutate(y_bar_lab=max(ypos)*1.1,
y_bar_bar=mean(y),
Lab=paste0("bar(y)[", Group, "]"))
# plot
ggplot(samp_df, aes(y, colour=Group)) +
# population in background
geom_hline(yintercept=0, linewidth=0.2, colour="grey80") +
geom_line(data=pop_df, aes(value, density), linewidth=1, alpha=0.25) +
geom_text(data=mu_df,
aes(x=mu, y=-density*0.05, label=Lab),
vjust=1, parse=T, size=6, alpha=0.5) +
geom_segment(data=mu_df,
aes(x=mu, xend=mu, y=-density*0.05, yend=0),
linewidth=2, linetype=1, alpha=0.5) +
# observations
geom_point(aes(y, ypos), size=2.5, shape=1, stroke=1) +
# group means
geom_segment(aes(x=y_bar, xend=y_bar, yend=0, y=y_bar_lab), linewidth=0.75, linetype=2) +
geom_text(aes(x=y_bar, y=y_bar_lab, label=Lab), nudge_y=0.02, parse=T, size=6) +
# grand mean
geom_segment(aes(x=y_bar_bar, xend=y_bar_bar, yend=0, y=y_bar_lab), colour="black",
linewidth=0.75, linetype=2) +
annotate("text", x=samp_df$y_bar_bar[1], y=samp_df$y_bar_lab[1]+0.02, label="bar(bar(y))",
parse=T, size=6, colour="black") +
# within-group variation
geom_segment(aes(x=y_bar, xend=y, y=ypos, yend=ypos),
arrow=arrow(ends="both", length=unit(0.05, "inches")), colour="black") +
# among-group variation
geom_segment(aes(x=y_bar, xend=y_bar_bar, y=y_bar_lab*0.9, yend=y_bar_lab*0.9),
arrow=arrow(ends="both", length=unit(0.1, "inches")),
colour="black", linewidth=1) +
scale_colour_manual(values=c("steelblue2", "firebrick"), guide="none") +
x_lim + y_lim +
labs(x="Response variable", y="")
