Count data that has two peaks. How would I model this? We did an experiment where people came in the lab and engaged in a helping task. They were told they could help with as many puzzles as they wanted to and a peer would finish them. The DV was the number of puzzles that were completed out of 20. We also had an experimental manipulation (IV was a categorical variable with two conditions). We expected people in the control group would do approximately half of the puzzles and those in the experimental group would complete more than half of the puzzles.  
I want to tackle this using a Bayesian approach, which is very new to me. My main roadblock is understanding the appropriate likelihood to use. The DV is a count variable that has a bimodal distribution around 10 and 20.
The data are here: 
puzzles <- c(3,9,7,9,20,13,20,12,10,10,12,10,13,20,10,15,10,20,9,10,20,3,7,10,12,20,7,20,20,16,16,12,20,6)
condition <- c('C','C','C','C','C','E','E','E','C','E','E','C','C','E','E','E','C','E','C','E','C','E','E','C','E','C','C','E','C','E','E','C','E','C')
df <- data.frame(puzzles,condition)

I found this post, which seems relevant, but I'm not quite sure it is the exact situation. Also, I think it is important that my data has multiple peaks to it. 
If it is relevant, I am using the brms package in R. 
Any help would be greatly appreciated!
 A: What you have is a bit similar to the count model in the vignette, under zero-inflated model. At least in the example above it seems overdispersed, you can check with a glm quasipoisson:
summary(glm(puzzles ~ condition ,family=quasipoisson,data=df))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.4398     0.1063  22.948   <2e-16 ***
conditionE    0.1908     0.1437   1.328    0.194    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 2.204227)

    Null deviance: 77.242  on 33  degrees of freedom
Residual deviance: 73.335  on 32  degrees of freedom
AIC: NA

The output above, there's a line saying "Dispersion parameter for quasipoisson family taken to be 2.204227". Normally dispersion parameter > 1 indicates over-dispersion,so you can model it in brm with a negative binomial:
fit <- brm(puzzles ~ condition ,family=negbinomial,data=df)
summary(fit)

 Family: negbinomial 
  Links: mu = log; shape = identity 
Formula: puzzles ~ condition 
   Data: df (Number of observations: 34) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept      2.43      0.10     2.23     2.64 1.00     3346     2582
conditionE     0.20      0.14    -0.09     0.47 1.00     3039     2736

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape    12.04      8.43     4.31    32.08 1.00     2655     1842

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Results in this case, seems to suggest the mean in groupE might be slightly higher than C. You can visualize the conditional means:
plot(conditional_effects(fit),points=TRUE)


The nature of your data does not require that many levels of modeling, so most likely you get conclusions similar to that of using a glm.. Or maybe you can elaborate what you would like to use brm for.
