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 $2sigma$) 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_srightarrow mu^+mu^-) = 3.81(18) times 10^{-9}$ and ${text{Br}}(B_drightarrow 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.