The nuclear-electronic orbital explicitly correlated Hartree-Fock (NEO-XCHF) approach is extended and applied to the positronic systems PsH, LiPs, and e(+)LiH. In this implementation, all electrons and positrons are treated quantum mechanically, and all nuclei are treated classically. This approach utilizes molecular orbital techniques with Gaussian basis sets for the electrons and positrons and includes electron-positron correlation with explicitly correlated Gaussian-type geminal functions. An efficient strategy is developed to reduce the number of variational parameters in the NEO-XCHF calculations. The annihilation rates, electron and positron densities, and electron-positron contact densities are compared to available results from higher-level calculations. Our analysis illustrates that the NEO-XCHF method produces qualitative to semi-quantitative results for these properties at a relatively low computational cost by treating only the essential electron-positron correlation explicitly. The NEO-HF method, which does not include explicit correlation and therefore is extremely efficient, is found to provide qualitatively accurate electron-positron contact densities for the e(+)LiH system but not for the LiPs system. Thus, the utility of the NEO-HF method for determining where annihilation occurs is system dependent and not generally reliable. The NEO-XCHF method, however, provides a computationally practical and reliable approach for determining where annihilation will occur in positronic systems.