Solution Exercise 2a

Data Simulation with Monte Carlo Methods

Marko Bachl

University of Hohenheim

Group exercise 2

Exercise 2: Welch’s t-test for count data

Part 1: False discovery rates

It is generally recommended to model count outcomes (e.g., number of likes of social media posts) with Poisson models (or negative binomial models, but that’s for another simulation).
However, if we only want to compare two groups on a count outcome, Welch’s t-test might be an easier option and should also do well:
- The mean is a consistent estimator of the rate parameter (CLT; see exercise 1a).
- Welch’s t-test is robust if the group standard deviations are unequal, which, by definition, would be the case if the counted event occurred at different rates in the two groups.

Exercise 2a: False discovery rates

We can test this assumption by simulation, similar to the comparison of Student’s and Welch’s t-test. More specifically, we first want to compare the false discovery rates of
- Student’s t-test,
- Welch’s t-test, and
- Poisson regression with a dummy predictor
if the true data generating process is drawing from a Poisson distribution.
Extend and adapt the previous simulation.

Simulation function

# No standard deviation ratio any more because SD fixed at sqrt(lambda)
# M_diff not necessary at this stage, but useful later
sim_ttest_glm = function(n = 200, GR = 1, lambda1 = 1, M_diff = 0) {
  n1 = round(n / (GR + 1))
  n2 = round(n1 * GR)
  lambda2 = lambda1 + M_diff
  g1 = rpois(n = n1, lambda = lambda1)
  g2 = rpois(n = n2, lambda = lambda2)
  Welch = t.test(g1, g2)$p.value
  Student = t.test(g1, g2, var.equal = TRUE)$p.value
  GLM = glm(outcome ~ group,
            data = data.frame(outcome = c(g1, g2),
                              group = c(rep("g1", n1), rep("g2", n2))),
            family = poisson)
  res = lst(Welch, Student, GLM = coef(summary(GLM))[2, 4])
  return(res)
}
sim_ttest_glm()

$Welch
[1] 0.7795355

$Student
[1] 0.7795354

$GLM
[1] 0.781524

Conditions

# It makes sense to vary group size ratio and lambda
conditions = expand_grid(GR = c(0.5, 1, 2), 
                         n = 600,
                         lambda1 = c(1, 10, 20),
                         M_diff = 0) %>% 
  rowid_to_column(var = "condition")
conditions

# A tibble: 9 × 5
  condition    GR     n lambda1 M_diff
      <int> <dbl> <dbl>   <dbl>  <dbl>
1         1   0.5   600       1      0
2         2   0.5   600      10      0
3         3   0.5   600      20      0
4         4   1     600       1      0
5         5   1     600      10      0
6         6   1     600      20      0
7         7   2     600       1      0
8         8   2     600      10      0
9         9   2     600      20      0

Run simulation

# use smaller i for now because additional glm() needs time
set.seed(35)
i = 1000
tic()
sims = map_dfr(1:i, ~ conditions) %>%
  rowid_to_column(var = "sim") %>%
  rowwise() %>%
  mutate(p.value = list(sim_ttest_glm(n = n, M_diff = M_diff,
                                 GR = GR, lambda1 = lambda1))) %>% 
  unnest_longer(p.value, indices_to = "method")
toc()

29.318 sec elapsed

Results: Histogram of p-values

sims %>% 
  ggplot(aes(p.value)) +
  geom_histogram(binwidth = 0.1, boundary = 0) + 
  facet_grid(method ~ lambda1 + GR, labeller = label_both) + 
  scale_x_continuous(breaks = c(1, 3)/4)

Results: Histogram of p-values

Results: QQ-Plot of p-values

sims %>% 
  ggplot(aes(sample = p.value, color = method)) +
  geom_qq(distribution = stats::qunif, geom = "line") + 
  geom_abline(slope = 1) + 
  facet_grid(GR ~ lambda1, labeller = label_both)

Results: QQ-Plot of p-values

Results: Rates of p < .05

Code

sims %>% 
  group_by(GR, lambda1, method) %>% 
  summarise(P_p05 = mean(p.value < 0.05)) %>% 
  ggplot(aes(factor(GR), P_p05, color = method, group = method)) + 
  geom_point() + geom_line() + 
  geom_hline(yintercept = 0.05, linetype = 2) + 
  facet_wrap(vars(lambda1), labeller = label_both) +
  labs(x = "Group size ratio", y = "Proportion p < 0.05")

Summary

All tests performed very similarly in terms of false discovery rates.
Not surprising: Previous simulation study showed that Student’s t-test breaks if group sizes and group standard deviations are unequal. Group standard deviations were equal by definition of the data-generating process.
t-tests and Poisson regression did also very similarly.