11 Mar

############################################################
# MULTIPLE TESTING CORRECTION DEMONSTRATION
# Purpose: Demonstrate multiple testing and corrections in R
############################################################

############################################################
# 1. SETUP
############################################################

[Link](123)

cat("Multiple Testing Demonstration\n")

cat("------------------------------\n")

############################################################
# 2. FALSE POSITIVES UNDER THE NULL
############################################################

# Number of hypothesis tests

n_tests <- 1000

# Simulate p-values when the null hypothesis is true

pvals <- runif(n_tests)

# Count how many are significant

false_positives <- sum(pvals < 0.05)

cat("Number of tests:", n_tests, "\n")

cat("Significant at p < 0.05:", false_positives, "\n\n")

############################################################
# 3. VISUALIZE P-VALUE DISTRIBUTION
############################################################

hist(pvals,
breaks = 20,
main = "P-value distribution under the null",
xlab = "p-value",
col = "lightblue")
############################################################
# 4. HOW FALSE POSITIVES SCALE
############################################################

test_sizes <- c(100, 1000, 10000)

for(n in test_sizes)
{
pvals <- runif(n)
fp <- sum(pvals < 0.05)

cat("Tests:", n,
" False positives:", fp,
" Expected:", n*0.05, "\n")
}

############################################################
# 5. FAMILY-WISE ERROR RATE (FWER)
############################################################

# Probability of at least one false positive

alpha <- 0.0005

n <- 100

fwer <- 1 - (1 - alpha)^n

cat("\nProbability of >=1 false positive with",

n, "tests:", fwer, "\n\n")

############################################################
# 6. BONFERRONI CORRECTION
############################################################

n_tests <- 1000

pvals <- runif(n_tests)

pvals<-c(0.0000000000000001,pvals)

# Bonferroni adjusted p-values

p_bonf <- [Link](pvals, method = "bonferroni")

significant_bonf <- sum(p_bonf < 0.05)

cat("Significant after Bonferroni:", significant_bonf, "\n")

par(mfrow=c(1,2))
hist(pvals,
breaks = 20,
main = "P-value distribution under the null",
xlab = "p-value",
col = "lightblue")

hist(p_bonf,
breaks = 200,
main = "P-value distribution after Bonferroni",
xlab = "p-value",
col = "lightblue")

c(p_bonf[p_bonf<1],0.0000000000000001)
############################################################
# 7. HOLM CORRECTION
############################################################

p_holm <- [Link](pvals, method = "holm")

significant_holm <- sum(p_holm < 0.05)

cat("Significant after Holm:", significant_holm, "\n")

c(p_holm[p_holm<1],p_bonf[p_bonf<1],0.0000000000000001)

############################################################
# 8. FALSE DISCOVERY RATE (BENJAMINI-HOCHBERG)
############################################################

p_bh <- [Link](pvals, method = "BH")

significant_bh <- sum(p_bh < 0.05)

cat("Significant after BH (FDR):", significant_bh, "\n")

############################################################
# 9. COMPARISON OF METHODS
############################################################

methods <- c("bonferroni", "holm", "BH")

for(m in methods)
{
adj <- [Link](pvals, method = m)
sig <- sum(adj < 0.05)

cat("Method:", m,
" Significant:", sig, "\n")
}

############################################################
# 10. PLOT RAW VS ADJUSTED P-VALUES
############################################################

adj <- [Link](pvals, method = "BH")

plot(pvals,
adj,
pch = 19,
cex = 0.5,
xlab = "Raw p-values",
ylab = "BH adjusted p-values",
main = "Raw vs Adjusted p-values")

abline(0,1,col="red")

############################################################
# 11. SIMULATION WITH TRUE SIGNAL
############################################################

[Link](123)

n_tests <- 1000

n_signal <- 100
n_null <- n_tests - n_signal
# null p-values
p_null <- runif(n_null)

# signal p-values (skewed toward zero)

p_signal <- rbeta(n_signal, 0.3, 1)

pvals <- c(p_null, p_signal)

############################################################
# 12. VISUALIZE SIGNAL VS NULL
############################################################

hist(pvals,
breaks = 30,
col = "lightgreen",
main = "P-values with True Signal",
xlab = "p-value")

############################################################
# 13. APPLY MULTIPLE CORRECTIONS
############################################################

methods <- c("bonferroni", "holm", "BH")

results <- [Link]()

for(m in methods)
{
adj <- [Link](pvals, method = m)

discoveries <- sum(adj < 0.05)

results <- rbind(results,

[Link](method=m,
discoveries=discoveries))
}

print(results)

############################################################
# 14. ESTIMATE FALSE DISCOVERIES
############################################################

# Suppose first 900 tests are null

null_indices <- 1:n_null

for(m in methods)
{
adj <- [Link](pvals, method = m)

sig <- which(adj < 0.05)

false_discoveries <- sum(sig %in% null_indices)

cat("\nMethod:", m, "\n")
cat("Discoveries:", length(sig), "\n")
cat("False discoveries:", false_discoveries, "\n")

if(length(sig) > 0)
{
cat("False discovery rate:",
false_discoveries/length(sig), "\n")
}
}

############################################################
# 15. GENOMICS SCALE EXAMPLE
############################################################

[Link](123)

n_genes <- 20000

n_true <- 500

p_null <- runif(n_genes - n_true)

p_signal <- rbeta(n_true, 0.2, 1)

pvals <- c(p_null, p_signal)

p_adj <- [Link](pvals, method = "BH")

discoveries <- sum(p_adj < 0.05)

cat("\nGENOME SCALE EXAMPLE\n")
cat("Genes tested:", n_genes, "\n")
cat("True signals:", n_true, "\n")
cat("Significant genes (FDR<0.05):", discoveries, "\n")

############################################################
# 16. SUMMARY TABLE
############################################################

methods <- c("bonferroni","holm","BH")

summary_table <- [Link]()

for(m in methods)
{
adj <- [Link](pvals, method=m)

discoveries <- sum(adj < 0.05)

summary_table <- rbind(summary_table,

[Link](Method=m,
Discoveries=discoveries))
}

print(summary_table)

############################################################
# END OF SCRIPT
############################################################

############################################################
# INTERACTIVE MULTIPLE TESTING SIMULATION
############################################################

simulate_tests <- function(n_tests = 10000,

prop_signal = 0.05,
effect_strength = 0.3,
alpha = 0.05)
{
n_signal <- round(n_tests * prop_signal)
n_null <- n_tests - n_signal

# Null p-values
p_null <- runif(n_null)

# Signal p-values
p_signal <- rbeta(n_signal, effect_strength, 1)

pvals <- c(p_null, p_signal)

methods <- c("none","bonferroni","holm","BH")

results <- [Link]()

for(m in methods)
{
if(m == "none")
{
adj <- pvals
} else {
adj <- [Link](pvals, method=m)
}

sig <- which(adj < alpha)

discoveries <- length(sig)

false_discoveries <- sum(sig <= n_null)

if(discoveries > 0)
{
fdr <- false_discoveries / discoveries
} else {
fdr <- 0
}

results <- rbind(results,

[Link](Method=m,
Discoveries=discoveries,
False_Discoveries=false_discoveries,
FDR=round(fdr,3)))
}

print(results)
 # Plot p-value distribution
hist(pvals,
breaks=40,
col="lightblue",
main=paste("P-value Distribution (",n_tests,"tests)",sep=""),
xlab="p-value")

############################################################
# RUN THE SIMULATION
############################################################

simulate_tests(
n_tests = 20000,
prop_signal = 0.05,
effect_strength = 0.3
)

13 Mar

#A=B=C:

x = rnorm(10, mean = 5, sd = 10)

y = rnorm(10, mean = 5, sd = 10)
z = rnorm(10, mean = 5, sd = 10)

abcBox <- [Link](

values = c(x,y,z),
variable = c(
rep("A", length(x)),
rep("B", length(y)),
rep("C", length(z))
)
)

anova_result <- aov(values ~ variable, data = abcBox)

summary(anova_result)

#A>B=C
x = rnorm(10, mean = 10, sd = 10)
y = rnorm(10, mean = 5, sd = 10)
z = rnorm(10, mean = 5, sd = 10)
abcBox <- [Link](
values = c(x,y,z),
variable = c(
rep("A", length(x)),
rep("B", length(y)),
rep("C", length(z))
)
)

anova_result <- aov(values ~ variable, data = abcBox)

summary(anova_result)
#A>B>C
x = rnorm(10, mean = 15, sd = 10)
y = rnorm(10, mean = 10, sd = 10)
z = rnorm(10, mean = 5, sd = 10)
abcBox <- [Link](
values = c(x,y,z),
variable = c(
rep("A", length(x)),
rep("B", length(y)),
rep("C", length(z))
)
)

anova_result <- aov(values ~ variable, data = abcBox)

summary(anova_result)

#A>B<C
x = rnorm(10, mean = 12, sd = 10)
y = rnorm(10, mean = 5, sd = 10)
z = rnorm(10, mean = 11, sd = 10)
abcBox <- [Link](
values = c(x,y,z),
variable = c(
rep("A", length(x)),
rep("B", length(y)),
rep("C", length(z))
)
)

anova_result <- aov(values ~ variable, data = abcBox)

summary(anova_result)
#A=B>C
x = rnorm(10, mean = 10, sd = 10)
y = rnorm(10, mean = 10, sd = 10)
z = rnorm(10, mean = 5, sd = 10)
abcBox <- [Link](
values = c(x,y,z),
variable = c(
rep("A", length(x)),
rep("B", length(y)),
rep("C", length(z))
)
)

anova_result <- aov(values ~ variable, data = abcBox)

summary(anova_result)

#A=C>B
x = rnorm(10, mean = 10, sd = 10)
y = rnorm(10, mean = 5, sd = 10)
z = rnorm(10, mean = 10, sd = 10)
abcBox <- [Link](
values = c(x,y,z),
variable = c(
rep("A", length(x)),
rep("B", length(y)),
rep("C", length(z))
)
)

anova_result <- aov(values ~ variable, data = abcBox)

summary(anova_result)

par(mfrow = c(2,3))

# 1) A = B = C
x = rnorm(10,5,10)
y = rnorm(10,5,10)
z = rnorm(10,5,10)

abcBox = [Link](
values=c(x,y,z),
variable=c(rep("A",10),rep("B",10),rep("C",10))
)
boxplot(values~variable,data=abcBox,main="A=B=C",
col = "pink")

# 2) A > B = C
x = rnorm(10,10,10)
y = rnorm(10,5,10)
z = rnorm(10,5,10)

abcBox = [Link](
values=c(x,y,z),
variable=c(rep("A",10),rep("B",10),rep("C",10))
)

boxplot(values~variable,data=abcBox,main="A>B=C",
col = "blue")

# 3) A > B > C
x = rnorm(10,15,10)
y = rnorm(10,10,10)
z = rnorm(10,5,10)

abcBox = [Link](
values=c(x,y,z),
variable=c(rep("A",10),rep("B",10),rep("C",10))
)

boxplot(values~variable,data=abcBox,main="A>B>C",
col = "yellow")

# 4) A > B < C
x = rnorm(10,12,10)
y = rnorm(10,5,10)
z = rnorm(10,11,10)

abcBox = [Link](
values=c(x,y,z),
variable=c(rep("A",10),rep("B",10),rep("C",10))
)

boxplot(values~variable,data=abcBox,main="A>B<C",
 col = "purple")

# 5) A = B > C
x = rnorm(10,10,10)
y = rnorm(10,10,10)
z = rnorm(10,5,10)

abcBox = [Link](
values=c(x,y,z),
variable=c(rep("A",10),rep("B",10),rep("C",10))
)

boxplot(values~variable,data=abcBox,main="A=B>C",
col = "green")

# 6) A = C > B
x = rnorm(10,10,10)
y = rnorm(10,5,10)
z = rnorm(10,10,10)

abcBox = [Link](
values=c(x,y,z),
variable=c(rep("A",10),rep("B",10),rep("C",10))
)

boxplot(values~variable,data=abcBox,main="A=C>B",
col = "red")

One way anova

16 Mar

##16th March 2026*#

data(ToothGrowth)

# Convert to factors
ToothGrowth$dose <- [Link](ToothGrowth$dose)
ToothGrowth$supp <- [Link](ToothGrowth$supp)

# Create extra factor for 3-way ANOVA

[Link](1)
ToothGrowth$group <- factor(sample(c("G1","G2"), nrow(ToothGrowth), replace=TRUE))

head(ToothGrowth)

#Additive Model (no interaction)

model_2_add <- aov(len ~ dose + supp, data=ToothGrowth)
summary(model_2_add)

#Multiplicative Model
model_2_mul <- aov(len ~ dose * supp, data=ToothGrowth)
summary(model_2_mul)

#Compare additive vs multiplicative

anova(model_2_add, model_2_mul)

#Pairwise t-test:
[Link](ToothGrowth$len,
interaction(ToothGrowth$dose, ToothGrowth$supp),
[Link]="bonferroni")

#Tukey Post-Hoc test

TukeyHSD(model_2_mul)

#THREE WAY ANOVA:

#Additive:
model_3_add <- aov(len ~ dose + supp + group, data=ToothGrowth)
summary(model_3_add)

#Multiplicative (with interactions)

model_3_mul <- aov(len ~ dose * supp * group, data=ToothGrowth)
summary(model_3_mul)

#Comparing:
anova(model_3_add, model_3_mul)

#Multiway ANOVA:
model_multi <- aov(len ~ dose * supp * group, data=ToothGrowth)
summary(model_multi)

#Tukey for multiway

TukeyHSD(model_multi)

summary(model_2_mul)
summary(model_3_mul)

anova(model_2_add, model_2_mul)
summary(model_2_mul)
#Two-way ANOVA was fitted using additive and multiplicative models.
#The interaction term dose:supp was significant (p < 0.05),
#therefore the multiplicative model was preferred over the additive model.

summary(model_3_mul)

##Three-way ANOVA was fitted using dose, supp and group.

#The factor group and all interactions involving group were not significant.
#Only the interaction between dose and supp was significant (p < 0.05).
#Therefore the full multiplicative three-way model was not required.
#The best model was the two-way multiplicative model with dose and supp.

par(mfrow=c(2,2))
plot(model_2_mul)
par(mfrow=c(1,1))
res <- residuals(model_2_mul)
fit <- fitted(model_2_mul)
[Link](res)
qqnorm(res)
qqline(res)

#Homogeneity of variance (Bartlett test)

[Link](len ~ interaction(dose,supp), data=ToothGrowth)

#Levene test (better than Bartlett)

library(car)

leveneTest(len ~ dose*supp, data=ToothGrowth)

#Final ANOVA validation check

anova(model_2_mul)
summary(model_2_mul)
TukeyHSD(model_2_mul)

#Diagnostic plots were examined for the multiplicative two-way ANOVA model.
#Residual plots showed no serious deviation from assumptions.
#Shapiro-Wilk test indicated normality of residuals.
#Bartlett / Levene test showed homogeneity of variance.
#Therefore ANOVA assumptions were satisfied.
#The final selected model was the two-way multiplicative ANOVA with interaction between dose
and supp.

#Interaction Plot:
[Link](ToothGrowth$dose,
ToothGrowth$supp,
ToothGrowth$len)
#Means plot
plot(tapply(ToothGrowth$len,
interaction(ToothGrowth$dose, ToothGrowth$supp),
mean),
type="b",
xlab="Dose-Supp",
ylab="Mean length")

#Boxplots by factors:
boxplot(len ~ dose * supp,
data = ToothGrowth,
col="lightblue")

#Residual vs fitted (separate)

plot(fitted(model_2_mul),
residuals(model_2_mul),
xlab="Fitted",
ylab="Residuals")

abline(h=0)

18 Mar

##ADDITIONAL------p value simulations

####################################
simulate_pvalues <- function(effect_size = 10, n_sim = 1000, add_outliers = FALSE) {

p_anova <- numeric(n_sim)

p_kw <- numeric(n_sim)

for (i in 1:n_sim) {

# Generate data
a <- rnorm(30, mean = 50, sd = 10)
b <- rnorm(30, mean = 50 + effect_size, sd = 10)
c <- rnorm(30, mean = 50 + 2*effect_size, sd = 10)
 # Add outliers if required
if (add_outliers) {
a <- c(a, 300)
b <- c(b, -200)
c <- c(c, 400)
}

data <- [Link](

values = c(a, b, c),
group = factor(rep(c("A","B","C"),
times = c(length(a), length(b), length(c))))
)

# ANOVA
p_anova[i] <- summary(aov(values ~ group, data = data))[[1]][["Pr(>F)"]][1]

# Kruskal-Wallis
p_kw[i] <- [Link](values ~ group, data = data)$[Link]
}

return(list(anova = p_anova, kw = p_kw))

}
[Link](123)

# Effect size = 10
res_small <- simulate_pvalues(effect_size = 10)

# Effect size = 50
res_large <- simulate_pvalues(effect_size = 50)

# With outliers
res_outliers <- simulate_pvalues(effect_size = 10, add_outliers = TRUE)
par(mfrow = c(3,2))

# Small effect size

hist(res_small$anova, main="ANOVA (Effect=10)", col="lightblue")
hist(res_small$kw, main="Kruskal (Effect=10)", col="lightgreen")

# Large effect size

hist(res_large$anova, main="ANOVA (Effect=50)", col="lightblue")
hist(res_large$kw, main="Kruskal (Effect=50)", col="lightgreen")

# With outliers
hist(res_outliers$anova, main="ANOVA (Outliers)", col="lightblue")
hist(res_outliers$kw, main="Kruskal (Outliers)", col="lightgreen")
mean(res_small$anova < 0.05)
mean(res_small$kw < 0.05)

mean(res_large$anova < 0.05)

mean(res_large$kw < 0.05)

mean(res_outliers$anova < 0.05)

mean(res_outliers$kw < 0.05)

20 Mar

library(ggplot2)
library(dplyr)

[Link](42)

# =========================================
# PART A: STATISTICAL POWER COMPARISON
# =========================================

simulate_power <- function(n = 20, effect_sizes = seq(0, 1.5, by = 0.15), n_sim = 500) {

results <- [Link]()

for (d in effect_sizes) {

t_sig <- 0
w_sig <- 0

for (i in 1:n_sim) {
a <- rnorm(n, mean = 0, sd = 1)
b <- rnorm(n, mean = d, sd = 1)

t_p <- [Link](a, b)$[Link]

w_p <- [Link](a, b)$[Link]

if (t_p < 0.05) t_sig <- t_sig + 1

if (w_p < 0.05) w_sig <- w_sig + 1
}

results <- rbind(results, [Link](

effect_size = d,
 test = "t-test",
power = t_sig / n_sim
))

results <- rbind(results, [Link](

effect_size = d,
test = "Wilcoxon",
power = w_sig / n_sim
))
}

return(results)
}

power_data <- simulate_power()

ggplot(power_data, aes(x = effect_size, y = power, color = test)) +

geom_line(size = 1.2) +
geom_point() +
labs(title = "Higher Statistical Power of Parametric vs Non-parametric Tests",
x = "Effect Size (Mean Difference)",
y = "Power") +
theme_minimal()

# =========================================
# PART B: OUTLIER SENSITIVITY
# =========================================

generate_outliers <- function(data, type) {

if (type == "extreme_high") {
data <- c(data, 100, 120, 150)
}

if (type == "extreme_low") {
data <- c(data, -100, -120, -150)
}

if (type == "both_sides") {
data <- c(data, -120, -100, 100, 120)
}

if (type == "skewed_cluster") {
 data <- c(data, rnorm(5, mean = 80, sd = 2))
}

if (type == "mild_outliers") {
data <- c(data, 65, 70, 75)
}

return(data)
}

outlier_types <- c("extreme_high", "extreme_low", "both_sides",

"skewed_cluster", "mild_outliers")

plot_list <- list()

for (otype in outlier_types) {

# BASE DATA (your structure)

a <- rnorm(20, mean = 50, sd = 10)
b <- rnorm(20, mean = 55, sd = 10)

# ADD OUTLIERS
a_out <- generate_outliers(a, otype)
b_out <- generate_outliers(b, otype)

# TESTS
t_p <- [Link](a_out, b_out)$[Link]
w_p <- [Link](a_out, b_out)$[Link]

df <- [Link](
values = c(a_out, b_out),
group = rep(c("A", "B"), c(length(a_out), length(b_out)))
)

p <- ggplot(df, aes(x = group, y = values)) +

geom_boxplot([Link] = NA) +
geom_jitter(width = 0.2, alpha = 0.5) +
ggtitle(paste0(otype,
"\n t-test p=", round(t_p, 4),
" | Wilcoxon p=", round(w_p, 4))) +
theme_minimal()

plot_list[[otype]] <- p
}

# Display plots
for (p in plot_list) print(p)

plot_list$extreme_high
plot_list$extreme_low
plot_list$both_sides
plot_list$skewed_cluster
plot_list$mild_outliers

# =========================================
# BONUS: DISTRIBUTION DISTORTION
# =========================================

for (otype in outlier_types) {

base <- rnorm(100, mean = 50, sd = 10)

altered <- generate_outliers(base, otype)

df <- [Link](
values = c(base, altered),
type = rep(c("Original", otype), c(length(base), length(altered)))
)

p <- ggplot(df, aes(x = values, fill = type)) +

geom_density(alpha = 0.4) +
ggtitle(paste("Distribution distortion:", otype)) +
theme_minimal()

print(p)
}

30 Mar

############################################################
# POLISHED TEACHING SCRIPT:
# Slope, spread, correlation coefficient, Spearman correlation,
# p-values, sample size, and nonlinear relationships
#
# Goal:
# Visually demonstrate how changes in:
# - slope
# - spread around a line
# - sample size
# - relationship shape
# affect:
# - Pearson correlation (r)
# - Spearman correlation (rho)
# - p-values
#
# Output:
# Figure 1: Conceptual overview panels
# Figure 2: Sample-size sweep
# Figure 3: Noise sweep
# Figure 4: Repeated simulation summary
#
# Author: ChatGPT
############################################################

##############################
# 1. Packages
##############################

# Install if needed:
# [Link](c("ggplot2", "dplyr", "patchwork", "tibble"))

library(ggplot2)
library(dplyr)
library(patchwork)
library(tibble)

##############################
# 2. Global settings
##############################

[Link](1234)

theme_teaching <- function(base_size = 12) {

theme_bw(base_size = base_size) +
theme(
[Link] = element_blank(),
[Link] = element_line(linewidth = 0.25),
[Link] = element_rect(fill = "grey95", colour = "black"),
[Link] = element_text(face = "bold"),
[Link] = element_text(face = "bold"),
[Link] = element_text(size = base_size * 0.95),
[Link] = element_text(face = "bold"),
 [Link] = element_text(face = "bold", size = base_size * 1.25)
)
}

##############################
# 3. Helper functions
##############################

simulate_linear <- function(n = 50,

slope = 1,
intercept = 0,
noise_sd = 1,
x_min = 0,
x_max = 10) {
x <- seq(x_min, x_max, [Link] = n)
y <- intercept + slope * x + rnorm(n, mean = 0, sd = noise_sd)
tibble(x = x, y = y)
}

simulate_nonlinear_quadratic <- function(n = 80,

noise_sd = 3,
x_min = -5,
x_max = 5) {
x <- seq(x_min, x_max, [Link] = n)
y <- x^2 + rnorm(n, mean = 0, sd = noise_sd)
tibble(x = x, y = y)
}

simulate_nonlinear_monotonic <- function(n = 80,

noise_sd = 0.25,
x_min = 0.5,
x_max = 10) {
x <- seq(x_min, x_max, [Link] = n)
y <- log(x) + rnorm(n, mean = 0, sd = noise_sd)
tibble(x = x, y = y)
}

get_stats <- function(df) {

pearson_test <- [Link](df$x, df$y, method = "pearson")
spearman_test <- [Link](df$x, df$y, method = "spearman", exact = FALSE)
lm_fit <- lm(y ~ x, data = df)

tibble(
slope_est = coef(lm_fit)[2],
 intercept_est = coef(lm_fit)[1],
pearson_r = unname(pearson_test$estimate),
pearson_p = pearson_test$[Link],
spearman_rho = unname(spearman_test$estimate),
spearman_p = spearman_test$[Link]
)
}

format_p <- function(p) {

if ([Link](p)) return("p = NA")
if (p < 2.2e-16) {
"p < 2.2e-16"
} else if (p < 0.001) {
paste0("p = ", format(p, scientific = TRUE, digits = 2))
} else {
paste0("p = ", signif(p, 3))
}
}

make_annotation <- function(stats_row,

include_slope = TRUE,
include_spearman = TRUE) {
pieces <- c()
if (include_slope) {
pieces <- c(pieces, paste0("slope = ", round(stats_row$slope_est, 2)))
}
pieces <- c(
pieces,
paste0("Pearson r = ", round(stats_row$pearson_r, 2)),
format_p(stats_row$pearson_p)
)
if (include_spearman) {
pieces <- c(
pieces,
paste0("Spearman rho = ", round(stats_row$spearman_rho, 2))
)
}
paste(pieces, collapse = "\n")
}

panel_plot <- function(df,

title = "",
subtitle = NULL,
show_lm = TRUE,
 show_loess = FALSE,
annotate_text = NULL,
tag = NULL) {

p <- ggplot(df, aes(x = x, y = y)) +

geom_point(size = 2, alpha = 0.8)

if (show_lm) {
p <- p + geom_smooth(method = "lm", se = FALSE, linewidth = 0.9)
}

if (show_loess) {
p <- p + geom_smooth(method = "loess", se = FALSE, linewidth = 0.9, linetype = 2)
}

p <- p +
labs(
title = title,
subtitle = subtitle,
x = "X",
y = "Y",
tag = tag
)+
theme_teaching()

if () {
x_pos <- min(df$x) + 0.03 * diff(range(df$x))
y_pos <- max(df$y) - 0.04 * diff(range(df$y))
p <- p + annotate(
"label",
x = x_pos,
y = y_pos,
hjust = 0,
vjust = 1,
size = 3.5,
label = annotate_text
)
}

p
}

##############################
# 4. Figure 1:
# Conceptual overview panels
##############################

# A. Same high correlation, shallow slope

df_A <- simulate_linear(n = 70, slope = 0.4, noise_sd = 0.30)
st_A <- get_stats(df_A)
ann_A <- make_annotation(st_A)

# B. Similar correlation, steeper slope

# Scaling both slope and noise keeps the pattern similarly tight
df_B <- simulate_linear(n = 70, slope = 1.8, noise_sd = 1.35)
st_B <- get_stats(df_B)
ann_B <- make_annotation(st_B)

# C. Same slope as B, but much more spread, so lower correlation

df_C <- simulate_linear(n = 70, slope = 1.8, noise_sd = 4.5)
st_C <- get_stats(df_C)
ann_C <- make_annotation(st_C)

# D. Nonlinear quadratic relationship:

# strong structure, but Pearson can be near 0
df_D <- simulate_nonlinear_quadratic(n = 90, noise_sd = 2.8)
st_D <- get_stats(df_D)
ann_D <- make_annotation(st_D)

pA <- panel_plot(
df_A,
title = "Shallow slope, tight spread",
subtitle = "High linear association despite shallow slope",
show_lm = TRUE,
annotate_text = ann_A,
tag = "A"
)

pB <- panel_plot(
df_B,
title = "Steeper slope, similarly tight spread",
subtitle = "Slope changes, but correlation can remain high",
show_lm = TRUE,
annotate_text = ann_B,
tag = "B"
)

pC <- panel_plot(
 df_C,
title = "Same slope, much wider spread",
subtitle = "More scatter around the line reduces correlation",
show_lm = TRUE,
annotate_text = ann_C,
tag = "C"
)

pD <- panel_plot(
df_D,
title = "Nonlinear relationship",
subtitle = "Pearson may be small even when the relationship is strong",
show_lm = TRUE,
show_loess = TRUE,
annotate_text = ann_D,
tag = "D"
)

fig1 <- (pA | pB) / (pC | pD) +

plot_annotation(
title = "Figure 1. Slope, spread, and the meaning of correlation",
subtitle = paste(
"Panels A and B show that slope alone does not determine correlation.",
"Panel C shows that greater spread around the line lowers correlation.",
"Panel D shows that Pearson correlation can miss nonlinear structure."
)
)

print(fig1)

ggsave("Figure1_conceptual_correlation_panels.png", fig1, width = 13, height = 10, dpi = 300)

ggsave("Figure1_conceptual_correlation_panels.pdf", fig1, width = 13, height = 10)

##############################
# 5. Figure 2:
# Sample-size sweep
##############################

sample_sizes <- c(8, 12, 20, 40, 80, 160)

sample_panels <- lapply(sample_sizes, function(nn) {

df <- simulate_linear(n = nn, slope = 0.8, noise_sd = 2.4)
st <- get_stats(df)
 df$panel <- paste0("n = ", nn)
df$annotation <- make_annotation(st)
df
})

sample_df <- bind_rows(sample_panels)

sample_labels <- sample_df %>%

group_by(panel) %>%
summarise(
x = min(x) + 0.03 * diff(range(x)),
y = max(y) - 0.04 * diff(range(y)),
annotation = first(annotation),
.groups = "drop"
)

fig2 <- ggplot(sample_df, aes(x = x, y = y)) +

geom_point(size = 1.8, alpha = 0.8) +
geom_smooth(method = "lm", se = FALSE, linewidth = 0.85) +
facet_wrap(~ panel, ncol = 3, scales = "fixed") +
geom_label(
data = sample_labels,
aes(x = x, y = y, label = annotation),
[Link] = FALSE,
hjust = 0,
vjust = 1,
size = 3.2,
[Link] = 0.25
)+
labs(
title = "Figure 2. Effect of sample size on correlation and p-value",
subtitle = paste(
"Slope and noise are similar across panels.",
"Observed p-values often become smaller as sample size increases,",
"even when the underlying pattern is conceptually similar."
),
x = "X",
y = "Y"
)+
theme_teaching()

print(fig2)

ggsave("Figure2_sample_size_effect.png", fig2, width = 12, height = 8.5, dpi = 300)

ggsave("Figure2_sample_size_effect.pdf", fig2, width = 12, height = 8.5)

##############################
# 6. Figure 3:
# Noise sweep
##############################

noise_levels <- c(0.3, 0.7, 1.5, 2.5, 4.0, 6.0)

noise_panels <- lapply(noise_levels, function(ns) {

df <- simulate_linear(n = 50, slope = 1.2, noise_sd = ns)
st <- get_stats(df)

df$panel <- paste0("noise SD = ", ns)

df$annotation <- make_annotation(st)
df
})

noise_df <- bind_rows(noise_panels)

noise_labels <- noise_df %>%

group_by(panel) %>%
summarise(
x = min(x) + 0.03 * diff(range(x)),
y = max(y) - 0.04 * diff(range(y)),
annotation = first(annotation),
.groups = "drop"
)

fig3 <- ggplot(noise_df, aes(x = x, y = y)) +

geom_point(size = 1.8, alpha = 0.8) +
geom_smooth(method = "lm", se = FALSE, linewidth = 0.85) +
facet_wrap(~ panel, ncol = 3, scales = "fixed") +
geom_label(
data = noise_labels,
aes(x = x, y = y, label = annotation),
[Link] = FALSE,
hjust = 0,
vjust = 1,
size = 3.2,
[Link] = 0.25
)+
labs(
title = "Figure 3. Effect of spread around the line",
 subtitle = paste(
"Slope is held constant. Increasing noise weakens the alignment of points",
"with the line, which usually lowers Pearson r and often also affects significance."
),
x = "X",
y = "Y"
)+
theme_teaching()

print(fig3)

ggsave("Figure3_noise_effect.png", fig3, width = 12, height = 8.5, dpi = 300)

ggsave("Figure3_noise_effect.pdf", fig3, width = 12, height = 8.5)

##############################
# 7. Figure 4:
# Repeated simulation summary
##############################

repeat_sim <- function(n = 20, slope = 0.6, noise_sd = 3, reps = 500) {

out <- vector("list", reps)

for (i in seq_len(reps)) {
df <- simulate_linear(n = n, slope = slope, noise_sd = noise_sd)
st <- get_stats(df)

out[[i]] <- tibble(

rep = i,
pearson_r = st$pearson_r,
pearson_p = st$pearson_p,
spearman_rho = st$spearman_rho,
spearman_p = st$spearman_p,
significant = st$pearson_p < 0.05
)
}

bind_rows(out)
}

rep_df <- repeat_sim(n = 20, slope = 0.6, noise_sd = 3, reps = 500)

p_r_hist <- ggplot(rep_df, aes(x = pearson_r)) +

geom_histogram(bins = 30) +
labs(
 title = "Observed Pearson correlation across repeated simulations",
subtitle = "The same underlying process can yield a range of observed r values",
x = "Observed Pearson r",
y = "Count"
)+
theme_teaching()

p_p_hist <- ggplot(rep_df, aes(x = pearson_p)) +

geom_histogram(bins = 30) +
geom_vline(xintercept = 0.05, linetype = 2, linewidth = 0.9) +
labs(
title = "Observed p-values across repeated simulations",
subtitle = "Significance can vary from sample to sample",
x = "Observed p-value",
y = "Count"
)+
theme_teaching()

fig4 <- p_r_hist | p_p_hist

print(fig4)

ggsave("Figure4_repeated_simulations.png", fig4, width = 12, height = 5.5, dpi = 300)

ggsave("Figure4_repeated_simulations.pdf", fig4, width = 12, height = 5.5)

##############################
# 8. Optional extra:
# Monotonic nonlinear example
# Pearson may still work somewhat, Spearman often tracks monotonicity well
##############################

df_extra <- simulate_nonlinear_monotonic(n = 80, noise_sd = 0.18)

st_extra <- get_stats(df_extra)
ann_extra <- make_annotation(st_extra)

fig_extra <- panel_plot(

df_extra,
title = "Optional example: monotonic but nonlinear relationship",
subtitle = "Spearman often captures monotonic structure better than Pearson",
show_lm = TRUE,
show_loess = TRUE,
annotate_text = ann_extra,
tag = "E"
)
print(fig_extra)

ggsave("Figure_extra_monotonic_nonlinear.png", fig_extra, width = 6.5, height = 5.5, dpi = 300)

ggsave("Figure_extra_monotonic_nonlinear.pdf", fig_extra, width = 6.5, height = 5.5)

##############################
# 9. Summary tables
##############################

summarise_case <- function(name, df, scenario) {

st <- get_stats(df)
tibble(
figure = name,
scenario = scenario,
n = nrow(df),
fitted_slope = round(st$slope_est, 3),
pearson_r = round(st$pearson_r, 3),
pearson_p = signif(st$pearson_p, 3),
spearman_rho = round(st$spearman_rho, 3),
spearman_p = signif(st$spearman_p, 3)
)
}

summary_table <- bind_rows(

summarise_case("Figure 1A", df_A, "Shallow slope, tight spread"),
summarise_case("Figure 1B", df_B, "Steep slope, tight spread"),
summarise_case("Figure 1C", df_C, "Steep slope, wide spread"),
summarise_case("Figure 1D", df_D, "Quadratic nonlinear"),
summarise_case("Extra", df_extra, "Monotonic nonlinear")
)

cat("\n================ SUMMARY TABLE ================\n")

print(summary_table)

sample_summary <- sample_df %>%

group_by(panel) %>%
summarise(
n = n(),
pearson_r = round(cor(x, y, method = "pearson"), 3),
spearman_rho = round(cor(x, y, method = "spearman"), 3),
pearson_p = signif([Link](x, y, method = "pearson")$[Link], 3),
.groups = "drop"
)
cat("\n=========== SAMPLE SIZE SUMMARY ===========\n")
print(sample_summary)

noise_summary <- noise_df %>%

group_by(panel) %>%
summarise(
n = n(),
pearson_r = round(cor(x, y, method = "pearson"), 3),
spearman_rho = round(cor(x, y, method = "spearman"), 3),
pearson_p = signif([Link](x, y, method = "pearson")$[Link], 3),
.groups = "drop"
)

cat("\n============= NOISE SUMMARY ===============\n")

print(noise_summary)

cat("\n======= REPEATED SIMULATION SUMMARY =======\n")

cat("Mean observed Pearson r: ", round(mean(rep_df$pearson_r), 3), "\n")
cat("SD observed Pearson r: ", round(sd(rep_df$pearson_r), 3), "\n")
cat("Mean observed p-value: ", round(mean(rep_df$pearson_p), 3), "\n")
cat("Proportion p < 0.05: ", round(mean(rep_df$significant), 3), "\n")

##############################
# 10. Interpretation notes
##############################

cat("
INTERPRETATION NOTES

1. Correlation is not the same thing as slope.

A shallow slope can still have very high correlation if the points lie tightly around a line.

2. Increasing spread around the line usually lowers the magnitude of Pearson correlation.

3. P-values depend on both:

- the strength of the observed association
- the sample size

4. A moderate correlation may be non-significant when n is small,

but highly significant when n is large.

5. Pearson correlation measures linear association.

It can miss strong nonlinear structure.
6. Spearman correlation is often useful when the relationship is monotonic
but not strictly linear.
\n")