  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-year data set with an injected stochastic gravitational wave
background (GWB). We then compute upper limits by sampling 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) \textit{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\% (53 of 500) of realizations. 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 upper limits obtained from
these 500 realizations is shifted to lower amplitude values.
