Iterative list-mode PET reconstruction requires the pre-calculation of a sensitivity image, in theory by back-projecting of all possible detector lines of response (LORs). For systems with a large number of LORs, or because of subject motion compensation, computing the exact sensitivity image is computationally expensive, so sub-sampling the LORs is attractive to reduce the work. The resulting inaccuracies in the sensitivity image propagate to the reconstructed emission image, which motivates looking at different sub-sampling schemes to reduce the errors. We extended the MOLAR reconstruction, originally implemented using uniform random (Uniform) subsampling of LORs, by adding a quasi-random (Halton) sequence based sampler and also a previously described Monte Carlo-based (MonteCarlo) non-uniform sampler. We also evaluated sampling using both techniques simultaneously (Halton+MonteCarlo). A simple post-smoothing of the sensitivity image can be used to improve variability in exchange for a loss in sensitivity spatial resolution. We included this method (Uniform Smoothed) in our comparison. Applying the five sampling methods (Uniform, Uniform-Smoothed, Halton, MonteCarlo, Halton+MonteCarlo), we obtain region of interest coefficients of variation (CoV) for replicated reconstructions of a human brain. Relative to Uniform, we observe two to three-fold improvement in CoV using Halton with a 4% increase in time to compute the sensitivity image, and slightly better improvements for Halton+MonteCarlo with a 163% increased computational investment. Halton and Halton+MonteCarlo produced lower CoVs than Uniform-Smoothed, without smoothing-induced bias in the reconstructed images.