In this paper, motivated by the observation that the Standard Model predictions are now above the experimental data for the mass difference ∆M s(d) , we perform a detailed study of B s(d) −B s(d) mixing and B s → µ + µ − decay in the Z 3 -invariant NMSSM with non-minimal flavour violation, using the recently developed procedure based on the Flavour Expansion Theorem, with which one can perform a purely algebraic mass-insertion expansion of an amplitude written in the mass eigenstate basis without performing any diagrammatic calculations in the interaction/flavour basis. Specifically, we consider the finite orders of mass insertions for neutralinos but the general orders for squarks and charginos, under two sets of assumptions for the squark flavour structures (i.e., while the flavour-conserving off-diagonal element δ LR 33 is kept in both of these two sectors, only the flavour-violating off-diagonal elements δ LL 23 and δ RR i3 (i = 1, 2) are kept in the LL and RR sectors, respectively). Our analytic results are then expressed directly in terms of the initial Lagrangian parameters in the interaction/flavour basis, making it easy to impose the experimental bounds on them. It is found numerically that the NMSSM effects with the above two assumptions for the squark flavour structures can accommodate the observed deviation for ∆M s(d) , while complying with the experimental constraints from the branching ratios of B s → µ + µ − and B → X s γ decays.