How to correctly simulate the idea of confidence interval in R

Question

I want to simulate/show the idea of confidence intervals but am not entirely sure if my approach is correct.

If I‘m not mistaken, a confidence interval can be understood as: if we would draw many samples, a certain proportion of these samples (e.g. 95% of them) would produce confidence intervals that encompass the true population parameter.

I want to show this idea dependent on two parameters: the number of samples I draw and the sample size.

My approach would be:

1. Create a vector of number of samples and sample sizes
2. Create some random population data
3. Calculate the true population mean
4. Take one of the sample sizes and loop through all number of samples. For all of these drawn samples I will calculate their mean, the standard error of that mean and then I will check if the true population mean lies within +/-1.96 times of that SE. Finally, I calculate the proportion for which this check is TRUE.
5. I repeat #4 for all other sample sizes.

And my code for that would be:

library(tidyverse)

m0 <- 2.5
std.d0 <- 1
pop <- rnorm(1e6, m0, std.d0)
m <- mean(pop)
std.d <- sum((m - pop)^2) / 1e6

n <- as.integer(c(5, 10, 15, 20, 30, 50, 100, 200, 500))
times <- as.integer(c(3, 5, 10, 30, 100, 200))

sampling <- function(pop, n, mean, sd, times)
{
  map(.x = times,
      .f = ~replicate(.x, sample(pop, size = n, replace = FALSE)) %>%
        as.data.frame() %>%
        pivot_longer(cols = everything()) %>%
        group_by(name) %>%
        summarize(draw_mean = mean(value),
                  draw_se = sd(value) / sqrt(n),
                  hit = between(!!m, draw_mean - 1.96 * draw_se, draw_mean + 1.96 * draw_se), .groups = 'drop') %>%
        summarize(perc_hit = mean(hit),
                  n = n,
                  rep = .x)) %>%
     bind_rows()
}

results <- map(.x = n,
    .f = ~sampling(pop = pop, n = .x, mean = m, sd = std.d, times = times)) %>%
  bind_rows() %>%
  print(n = Inf)

results %>%
  mutate(rep = as.factor(rep),
         n = as.factor(n)) %>%
  ggplot() +
  geom_point(aes(x = rep, y = perc_hit, color = n), position = position_jitter(width = 0.1))