  We calculate the bag parameters for neutral $B$-meson mixing in and beyond
the Standard Model, in full four-flavour lattice QCD for the first time. We
work on gluon field configurations that include the effect of $u$, $d$, $s$ and
$c$ sea quarks with the Highly Improved Staggered Quark (HISQ) action at three
values of the lattice spacing and with three $u/d$ quark masses going down to
the physical value. The valence $b$ quarks use the improved NRQCD action and
the valence light quarks, the HISQ action. Our analysis was blinded. Our
results for the bag parameters for all five operators are the most accurate to
date. For the Standard Model operator between $B_s$ and $B_d$ mesons we find:
$\hat{B}_{B_s}=1.232(53)$, $\hat{B}_{B_d}=1.222(61)$. Combining our results
with lattice QCD calculations of the decay constants using HISQ quarks from the
Fermilab/MILC collaboration and with experimental values for $B_s$ and $B_d$
oscillation frequencies allows determination of the CKM elements $V_{ts}$ and
$V_{td}$. We find $V_{ts} = 0.04189(93)$, $V_{td} = 0.00867(23)$ and
$V_{ts}/V_{td} = 0.2071(27)$. Our results agree well (within $2\sigma$) with
values determined from CKM unitarity constraints based on tree-level processes
(only). Using a ratio to $\Delta M$ in which CKM elements cancel in the
Standard Model, we determine the branching fractions
${\text{Br}}(B_s\rightarrow \mu^+\mu^-) = 3.81(18) \times 10^{-9}$ and
${\text{Br}}(B_d\rightarrow \mu^+\mu^-) = 1.031(54) \times 10^{-10}$. We also
give results for matrix elements of the operators $R_0$, $R_1$ and
$\tilde{R}_1$ that contribute to neutral $B$-meson width differences.
