  The search for flavour-changing neutral current effects in $B$-meson decays
is a powerful probe of physics beyond the Standard Model. Deviations from SM
behaviour are often quantified by extracting the preferred values of the Wilson
coefficients of an operator product expansion. We use the FlavBit module of the
GAMBIT package to perform a simultaneous global fit of the Wilson coefficients
$C_7$, $C_9$, and $C_{10}$ using a combination of all current data on $b \to s
\mu^+ \mu^-$ transitions. We further extend previous analyses by accounting for
the correlated theoretical uncertainties at each point in the Wilson
coefficient parameter space, rather than deriving the uncertainties from a
Standard Model calculation. We find that the best fit deviates from the SM
value with a significance of 6.6$\sigma$. The largest deviation is associated
with a vector coupling of muons to $b$ and $s$ quarks.
