babies
walk program 1 9 Special 2 9 Special 3 10 Special 4 10 Special 5 13 Special 6 11 Regular 7 10 Regular 8 10 Regular 9 12 Regular 10 10 Regular 11 11 Weekly report 12 12 Weekly report 13 9 Weekly report 14 12 Weekly report 15 13 Weekly report 16 13 Single report 17 11 Single report 18 12 Single report 19 13 Single report 20 11 Single report
library(ggplot2) library(ggbeeswarm) babies |> ggplot(aes(x = program, y = walk)) + geom_beeswarm(cex = 3, size = 3)
library(dplyr) group_means <- babies |> group_by(program) |> summarize(avg = mean(walk)) group_means
# A tibble: 4 × 2 program avg <fct> <dbl> 1 Special 10.2 2 Regular 10.6 3 Weekly report 11.4 4 Single report 12
group_means$avg[1] - group_means$avg[2]
[1] -0.4
group_means$avg[1] - group_means$avg[3]
[1] -1.2
group_means$avg[1] - group_means$avg[4]
[1] -1.8
group_means$avg[2] - group_means$avg[3]
[1] -0.8
group_means$avg[2] - group_means$avg[4]
[1] -1.4
group_means$avg[3] - group_means$avg[4]
[1] -0.6
overall_mean <- mean(babies$walk) a_hat1 <- group_means$avg[1] - overall_mean a_hat2 <- group_means$avg[2] - overall_mean a_hat3 <- group_means$avg[3] - overall_mean a_hat4 <- group_means$avg[4] - overall_mean c(a_hat1, a_hat2, a_hat3, a_hat4)
[1] -0.85 -0.45 0.35 0.95
The base R engine for randomization: sample().
sample()
a <- c(3, 5, 7, 9) sample(x = a, size = 2, replace = FALSE)
[1] 3 9
sample(a, 4)
[1] 5 3 9 7
sample(a)
[1] 5 9 3 7
shuffle4 <- function(y, z) { Z <- sample(z) # source of randomness babies_sim <- data.frame(walk = y, program = Z) group_means <- babies_sim |> group_by(program) |> summarize(avg = mean(walk)) overall_mean <- mean(babies_sim$walk) stats <- c(group_means$avg[1] - overall_mean, group_means$avg[2] - overall_mean, group_means$avg[3] - overall_mean, group_means$avg[4] - overall_mean) stats }
null_stats <- replicate(5, shuffle4(y = babies$walk, z = babies$program)) null_stats
[,1] [,2] [,3] [,4] [,5] [1,] 0.15 -0.05 0.35 0.15 0.35 [2,] -0.85 0.15 -0.45 0.15 -0.45 [3,] 0.75 -0.05 0.75 -0.25 -0.05 [4,] -0.05 -0.05 -0.65 -0.05 0.15
To access simulations from \(\hat{\alpha}_1\):
null_stats[1, ] # extract first row as a vector.
[1] 0.15 -0.05 0.35 0.15 0.35
null_stats <- replicate(500, shuffle4(y = babies$walk, z = babies$program)) null_1 <- null_stats[1, ] null_2 <- null_stats[2, ] null_3 <- null_stats[3, ] null_4 <- null_stats[4, ]
null1 <- data.frame(stat = null_1) ggplot(null1, aes(x = stat)) + geom_histogram()
null1 <- data.frame(stat = null_1) ggplot(null1, aes(x = stat)) + geom_histogram() + geom_vline(xintercept = a_hat1, color = "tomato", lwd = 2)
null1 |> summarize(pval = mean(stat < a_hat1) * 2)
pval 1 0.088
What is our decision regarding \(H_0\)? . . .
At threshold of .05, we fail to reject the null.
null3 <- data.frame(stat = null_3) ggplot(null3, aes(x = stat)) + geom_histogram() + geom_vline(xintercept = a_hat3, color = "tomato", lwd = 2)
babies |> group_by(program) |> summarize(var = var(walk))
# A tibble: 4 × 2 program var <fct> <dbl> 1 Special 2.7 2 Regular 0.8 3 Weekly report 2.3 4 Single report 1
babies |> group_by(program) |> summarize(var = var(walk)) |> summarize(avg_var = mean(var))
# A tibble: 1 × 1 avg_var <dbl> 1 1.7
avg_var <- babies |> group_by(program) |> summarize(var = var(walk)) |> summarize(avg_var = mean(var)) |> pull(avg_var) avg_var
[1] 1.7
babies |> group_by(program) |> summarize(avg = mean(walk))
babies |> group_by(program) |> summarize(avg = mean(walk)) |> summarize(var_avg = var(avg))
# A tibble: 1 × 1 var_avg <dbl> 1 0.650
var_avg <- babies |> group_by(program) |> summarize(avg = mean(walk)) |> summarize(var_avg = var(avg)) |> pull(var_avg) var_avg
[1] 0.65
f_stat <- var_avg / avg_var
shuffle4_f <- function(y, z) { Z <- sample(z) # source of randomness babies_sim <- data.frame(walk = y, program = Z) avg_var <- babies_sim |> group_by(program) |> summarize(var = var(walk)) |> summarize(avg_var = mean(var)) |> pull(avg_var) var_avg <- babies_sim |> group_by(program) |> summarize(avg = mean(walk)) |> summarize(var_avg = var(avg)) |> pull(var_avg) var_avg / avg_var }
null_stats <- replicate(5, shuffle4_f(y = babies$walk, z = babies$program)) null_stats
[1] 0.12048193 0.06590038 0.34095238 0.22982456 0.44923077
null_stats <- replicate(500, shuffle4_f(y = babies$walk, z = babies$program)) null <- data.frame(stats = null_stats) ggplot(null, aes(x = stats)) + geom_histogram() + geom_vline(xintercept = f_stat, color = "tomato", lwd = 2)
null |> summarize(pval = mean(stats > f_stat))
pval 1 0.154
The workhorse function: aov().
aov()
babies_anova <- aov(walk ~ factor(program), data = babies) anova_table <- summary(babies_anova) anova_table
Df Sum Sq Mean Sq F value Pr(>F) factor(program) 3 9.75 3.25 1.912 0.168 Residuals 16 27.20 1.70
obs_f_stat <- anova_table[[1]]$`F value`[1] obs_f_stat
[1] 1.911765
ggplot(null_stats, aes(x = f_vec)) + geom_histogram(aes(y = ..density..)) + stat_function(fun = df, args = list(df1 = 3, df2 = 16), color = "goldenrod", size = 2) + geom_vline(xintercept = obs_f_stat, color = "tomato", lwd = 2)
