We present an explicit detailed theoretical and observational investigation of an anisotropic massive Brans-Dicke (BD) gravity extension of the standard $Lambda$CDM model, wherein the extension is characterized by two additional degrees of freedom; the BD parameter, $omega$, and the present day density parameter corresponding to the shear scalar, $Omega_{sigma^2,0}$. The BD parameter, determining the deviation from general relativity (GR), by alone characterizes both the dynamics of the effective dark energy (DE) and the redshift dependence of the shear scalar. These two affect each other depending on $omega$, namely, the shear scalar contributes to the dynamics of the effective DE, and its anisotropic stress --which does not exist in scalar field models of DE within GR-- controls the dynamics of the shear scalar deviating from the usual $propto(1+z)^6$ form in GR. We mainly confine the current work to non-negative $omega$ values as it is the right sign --theoretically and observationally-- for investigating the model as a correction to the $Lambda$CDM. By considering the current cosmological observations, we find that $omegagtrsim 250$, $Omega_{sigma^2,0}lesssim 10^{-23}$ and the contribution of the anisotropy of the effective DE to this value is insignificant. We conclude that the simplest anisotropic massive BD gravity extension of the standard $Lambda$CDM model exhibits no significant deviations from it all the way to the Big Bang Nucleosynthesis. We also point out the interesting features of the model in the case of negative $omega$ values; for instance, the constraints on $Omega_{sigma^2,0}$ could be relaxed considerably, the values of $omegasim-1$ (relevant to string theories) predict dramatically different dynamics for the expansion anisotropy.