Testing for difference in means with bimodally distributed data? I have two bimodal distributions of data with two peaks (one around 0 and the other around 1). I have provided an example of one of the distributions.

Although their means and variances are different, they both have peaks around the same x-axis markers. How do I test for difference in means using statistically rigorous approaches? Can I still use a t-test even though the data is not normally distributed?
 A: This looks like a textbook example of using beta regression, which can model values between 0 and 1. It is often used for modeling rates and proportions. The distribution can take on various different shapes, and your distribution looks a lot like a beta distribution with shape parameters of .1 and .2. In R:
set.seed(1839); hist(rbeta(1000, .1, .2)) 


I can simulate two samples from different beta distributions, one with a higher mean, and then do a significance test using the betareg function from the betareg package in R:
set.seed(1839)
y <- c(rbeta(300, .1, .3), rbeta(300, .2, .1))
y <- ifelse(y == 1, .999999, y) # because values must be BETWEEN 0 and 1
x <- rep(c(1, 0), each = 300)
mod <- betareg::betareg(y ~ x)
summary(mod)

This returns:
Call:
betareg::betareg(formula = y ~ x)

Standardized weighted residuals 2:
    Min      1Q  Median      3Q     Max 
-3.8917 -0.5212  0.0366  0.5800  3.0307 

Coefficients (mean model with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.55237    0.08089   6.828 8.59e-12 ***
x           -1.48921    0.12093 -12.315  < 2e-16 ***

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)    
(phi)  0.33130    0.01566   21.16   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood:  2806 on 3 Df
Pseudo R-squared: 0.2134
Number of iterations: 19 (BFGS) + 1 (Fisher scoring) 

I have written explainers on beta regression, with some citations you can check out, both here (using betareg) and here (using gamlss)

But what if you were to assume normality and do a t-test? Usually, people worry about Type I error with violating assumptions. I simulated data where there was no difference at all between conditions 1000 times with a sample size of 1000:
set.seed(1839)
iter <- 1000
n <- 1000
results <- lapply(seq_len(iter), function(zzz) {
  y <- rbeta(n, .2, .3)
  y <- ifelse(y == 1, .999999, y)
  x <- rbinom(n, 1, .5)
  beta_reg <- summary(betareg::betareg(y ~ x))$coef$mean[2, 4] < .05
  t_test <- t.test(y ~ x)$p.value < .05
  c(beta_reg = beta_reg, t_test = t_test)
})
colMeans(do.call(rbind, results))

Type I error rates show that it is about .05 for both. So, you're probably fine doing a t-test. But (a) you should ask yourself if the two poles are worth modeling themselves—you could trichotomize the data into low pole/between/high pole if those poles are theoretically or practically meaningful—and (b) beta regression can provide a sanity check; do both a t-test and beta regression. If they agree, then report that.
