$\newcommand\la\lambda\newcommand\La\Lambda\newcommand{\Ga}{\Gamma}\newcommand{\R}{\mathbb R}$Some preliminaries: To make it explicit that the distribution of $Y$ depends on the parameter $r$, let $Y_r:=Y$, so that
\begin{equation*}
\Pr(Y_r=j) = \binom{-r}j (-q)^j p^r \tag{1}\label{1}
\end{equation*}
for $j=0,1,\dots$. We know that the number $N$ of purchasers has the Poisson distribution with some parameter $\mu$.
We have
\begin{equation*}
Y_r=\sum_{i=1}^N X_i,
\end{equation*}
where $X_i$ is the number of Items bought by the $i$th purchaser. Suppose that the purchasers act independently and are statistically interchangeable, so that the $X_i$'s are iid. We also suppose that the number $N$ of purchasers is independent of the numbers $X_i$ of Items bought by purchasers. It follows that the distribution of $Y_r$ is a compound Poisson distribution:
\begin{equation*}
P_{Y_r}=\sum_{k=0}^\infty \Pr(N=k)\,P_X^{*k}
=\sum_{k=0}^\infty \frac{\mu^k}{k!}\,e^{-\mu}\,P_X^{*k}, \tag{2}\label{2}
\end{equation*}
where $X:=X_1$, $P_Z$ denotes the distribution of a random variable (r.v.) $Z$, and $Q^{*k}:=\underbrace{Q*\cdots*Q}_{k}$.
We have to plausibly/"intuitively" recover the distribution $P_X$ of $X$ given \eqref{2}. In view of the infinite divisibility of the negative binomial distribution of $Y_r$ and by \eqref{2}, we have
\begin{equation*}
\lim_{n\to\infty}P_{Y_{r/n}}^{*n}=P_{Y_r}
=\lim_{n\to\infty}\big(P_0+\tfrac\mu n\,(P_X-P_0)\big)^{*n}.
\end{equation*}
So, it is plausible that for large $n$
\begin{equation*}
(1-\tfrac\mu n)P_0+\tfrac\mu n\,P_X=P_0+\frac\mu n\,(P_X-P_0)
\approx P_{Y_{r/n}}.
\end{equation*}
Therefore and because $X\ge1$, we see that $P_X$ is the limit as $n\to\infty$ of the conditional distribution of $Y_{r/n}$ given $Y_{r/n}\ne0$. So, for integers $j\ge1$, in view of \eqref{1},
\begin{equation*}
\Pr(X=j)=\lim_{s\downarrow0}\frac{\Pr(Y_s=j)}{1-\Pr(Y_s=0)} \\
=\lim_{s\downarrow0}\frac{(j-1+s)(j-2+s)\cdots(1+s)s\; q^j p^s/j!}{1-p^s}
=\frac{q^j/j}{-\ln p},
\end{equation*}
as desired.
$$\mbox{********************************}$$
Answer to the changed question
In the initial version of the question, the negative binomial distribution of $Y_r$ appeared to be assumed, and the question appeared to be asking for a more intuitive derivation of the logarithmic series distribution for $X$ based on that "negative binomial" assumption.
After the above answer was posted, the OP changed the question, which now calls for a derivation of the logarithmic series distribution for $X$ without assuming the negative binomial distribution of $Y_r$.
Let us try to answer the changed question as well. Here I will be following lines of Fisher's paper (pp. 54--58).
Let $Z$ denote the random number of books purchased by a potential buyer (not necessarily an actual purchaser). Assume that $Z$ depends on some "hidden" positive r.v. $\La$ so that the conditional distribution of $Z$ given $\La$ is the Poisson distribution with mean $\La$:
\begin{equation*}
\Pr(Z=z\mid\La)=\frac{\La^z}{z!}\,e^{-\La}\,1(z=0,1,\dots). \tag{10}\label{10}
\end{equation*}
In turn, assume that $\La$ has the gamma distribution with positive parameters $a$ and $b$, so that
\begin{equation*}
\Pr(\La\in d\la)=\frac1{\Ga(a)b^a}\,\la^{a-1}e^{-\la/b}\,1(\la>0)\,d\la. \tag{20}\label{20}
\end{equation*}
Then
\begin{align}
\Pr(Z=z)&=\int_\R\Pr(\La\in d\la)\Pr(Z=z\mid\La) \notag \\
&=\frac{\Ga(a+z)}{\Ga(a)}\frac1{z!}\,\frac{1}{b^a (1+1/b)^{a+z}}\,1(z=0,1,\dots) \notag \\
&=\frac{a(a+1)\cdots(a+z-1)}{z!}\,p^a(1-p)^z\,1(z=0,1,\dots), \tag{30}\label{30}
\end{align}
where
\begin{equation*}
p:=\frac1{1+b}.
\end{equation*}
The "Poisson" assumption \eqref{10} seems natural.
Unfortunately, just as in the mentioned paper by Fisher, I cannot offer a reason to model the mixing distribution of $\La$ as a gamma distribution, except that then it is easy to get \eqref{30}. One may note here, though, that the gamma distribution is the conjugate prior distribution with respect to the Poisson family of distribution -- a fact prized by Bayesians.
Anyhow, once this modeling is accepted, we next note that the expected number of purchases by a potential buyer is $EZ=EE(Z\mid\La)=E\La=ab$, which may be reasonably assumed small. If $b>0$ is fixed (that is, if $p=\frac1{1+b}\in(0,1)$ is fixed), then it follows that $a>0$ is small.
So, the distribution of $Z$ will be close to the Dirac distribution supported on the singleton set $\{0\}$. However, since we are interested in those potential buyers who are actually purchasers, we consider the conditional distribution of $Z$ given $Z>0$. In view of \eqref{30}, the probability mass function (pmf) -- say $f_{a,p}$ -- of the just mentioned conditional distribution is given by the formula
\begin{align*}
& f_{a,p}(z)=\frac{\Pr(Z=z)\,1(z>0)}{1-\Pr(Z=0)} \\[6pt]
= {} & \frac{a(a+1)\cdots(a+z-1)}{z!}\,\frac{p^a(1-p)^z}{1-p^a} \, 1(z=1,2,\dots).
\end{align*}
Since $a$ is small, it is natural to consider the limit
\begin{equation*}
f_{0+,p}(z):=\lim_{a\downarrow0}f_{a,p}(z)
=\frac1{-\ln p}\frac{(1-p)^z}z\,1(z=1,2,\dots).
\end{equation*}
So, $f_{0+,p}$ is the pmf of the logarithmic series distribution, as desired.
