# Fitting a (maybe Poisson) regression to my data in R

I'm trying to understand the probability distribution for my data. Each point represents an individual patient, the 0 or 1 label indicates whether they have a particular disease, and the number of detections indicates how many times each activity was detected in a patient. I only have 35 data points.

Basically, I want

$$P(\textrm{number of detections} | \textrm{label})$$.

There isn't enough data here to use something like fitdistr in R, I think, so I'm looking to build a glm in R with the poisson family, i.e.

mdl <- glm(detections ~ label, data = ata, family = poisson(link='log'))

but this isn't giving a great model (I'm guessing because the data is overdispersed).

Any suggestions for what to try next? A negative binomial distribution? Quasi-poisson? And suggestions for goodness-of-fit?