The Lattice Boltzmann Method (LBM) is a numerical method based on computational statistical mechanics that is well-suited for approximating complex flow behaviors such as non-Newtonian, free surface, and multiphase multicomponent flow. LBM is typically applied to simulate flow through a series of time steps, each consisting of streaming particle distributions to neighboring nodes and collisions of particle distributions at each node through a collision operator. The collision operator is of interest because it, along with the equilibrium distribution function, determines the physics that are simulated (e.g. constitutive laws, interfacial dynamics, etc.) and it has implications on numerical stability and computational efficiency. This work examines various collision operators and methods for stability enhancement for their suitability for simulating non-Newtonian fluid flows in terms of their accuracy, numerical stability and computational efficiency. The investigation was carried out as a numerical study looking for qualitative, yet practical, results, including testing the BGK and MRT collision operators, with and without entropic filtering, as applied to Bingham plastics and power-law fluids. Two different benchmark problems were chosen for the flows: Poiseuille flow, and lid-driven square cavity flow. The results of the numerical study showed that the MRT collision operator can have an advantage in terms of stability and accuracy for a variety of non-Newtonian flow behaviors, but at an increased computing cost that was, in some cases, as much as five times greater than the BGK collision operator. It was also shown that, although it introduces error in the constitutive response of the fluid (and therefore, may not accurately capture the physics of the flow), artificial dissipation can be an effective technique for stabilizing the numerics of non-Newtonian, lid-driven cavity flow simulations, and is particularly effective for stabilizing shear-thinning fluids.