Abstract Recently, many different pulsar timing array (PTA) collaborations have reported strong evidence for a common stochastic process in their data sets. The reported amplitudes are in tension with previously computed upper limits. In this paper, we investigate how using a subset of a set of pulsars biases Bayesian upper limit recovery. We generate 500 simulated PTA data sets, based on the NANOGrav 11 yr data set with an injected stochastic gravitational-wave background (GWB). We then compute the upper limits by sampling the individual pulsar likelihoods, and combine them through a factorized version of the PTA likelihood to obtain upper limits on the GWB amplitude, using different numbers of pulsars. We find that it is possible to recover an upper limit (95% credible interval) below the injected value, and that it is significantly more likely for this to occur when using a subset of pulsars to compute the upper limit. When picking pulsars to induce the maximum possible bias, we find that the 95% Bayesian upper limit recovered is below the injected value in 10.6% of the realizations (53 of 500). Further, we find that if we choose a subset of pulsars in order to obtain a lower upper limit than when using the full set of pulsars, the distribution of the upper limits obtained from these 500 realizations is shifted to lower-amplitude values.