R Fitting data from randomization test to simple linear regression In this dataset, there are variable groups referring to math and physics and variable improvement_score
A hypothesis test on two groups can be conducted to compare mean score improvements after online training with physics and mathematics.
test_stat<-as.numeric(my_data %>%
  group_by(group) %>%
  summarise(means = mean(improvement_score), .groups='drop') %>% 
    # printed, but doesn't change behaviour
  summarise(value = diff(means)))

# conduct randomization test
simulated_values <- rep(NA, 1000)

for (i in 1:1000) {
  sim <- study_dat %>% mutate(group = sample(group))
  sim_value <- sim %>%
    group_by(group) %>%
    summarise(means = mean(improvement_score), .groups='drop') %>%
    summarise(value = diff(means))
  simulated_values[i] <- as.numeric(sim_value)
}

sim <- tibble(mean_diff = simulated_values)

ggplot(sim, aes(x=mean_diff)) +
  geom_histogram(col="black",fill="gray", binwidth=0.2) +
  labs(x = "Simulated mean improvements differences (assuming H0)")


sim <- tibble(mean_diff = simulated_values)
sim %>%
  filter(mean_diff >= abs(test_stat) |
           mean_diff <= -1*abs(test_stat)) %>%
  summarise(p_value = n() / 1000)

The estimated p-value based on this randomization test is 0 so there is very strong evidence against the hypothesis that the mean Score improvement is the same for those training online physics and math.
I can't understand how to use a simple linear regression model to test for a difference in mean score improvements for those training online physics and math.
Any guesses how the code should look like?

The t-test can be framed as a comparison of regression models.
Let's make a simple model L0 that always predicts the pooled mean taken from both math and physics. Let's make a more complex model that predicts the mean of the math scores if the test was a math test and the mean of the physics scores if the test was a physics test. In other words, we do a regression on an indicator variable for the type of test.
We then compare the two models and get a p-value for the model comparison using something called a Wald test.
An R simulation will show these approaches to be equivalent.
install.packages('lmtest')
set.seed(2021)
N <- 10
x <- c(rep(0, N), rep(1, N + 7)) # unequal group sizes
y <- x + rnorm(length(x))
L1 <- lm(y ~ x)
L0 <- lm(y ~ 1)
lmtest::waldtest(L0, L1)$`Pr(>F)`[2] # = 0.115670428569383
t.test(y[x == 0], y[x == 1], var.equal = T)$p.value # = 0.115670428569383

This is what the summary of L1 does to calculate the p-value for the variable.
summary(L1)$coef[2, 4] # = 0.115670428569383

