  We consider random singlet phases of spin-$\frac{1}{2}$, random,
antiferromagnetic spin chains, in which the universal leading-order divergence
$\frac{\ln 2}{3}\ln\ell$ of the average entanglement entropy of a block of
$\ell$ spins, as well as the closely related leading term $\frac{2}{3}l^{-2}$
in the distribution of singlet lengths are well known by the strong-disorder
renormalization group (SDRG) method. Here, we address the question of how large
the subleading terms of the above quantities are. By an analytical calculation
performed along a special SDRG trajectory of the random XX chain, we identify a
series of integer powers of $1/l$ in the singlet-length distribution with the
subleading term $\frac{4}{3}l^{-3}$. Our numerical SDRG analysis shows that,
for the XX fixed point, the subleading term is generally $O(l^{-3})$ with a
non-universal coefficient and also reveals terms with half-integer powers:
$l^{-7/2}$ and $l^{-5/2}$ for the XX and XXX fixed points, respectively. We
also present how the singlet lengths originating in the SDRG approach can be
interpreted and calculated in the XX chain from the one-particle states of the
equivalent free-fermion model. These results imply that the subleading term
next to the logarithmic one in the entanglement entropy is $O(\ell^{-1})$ for
the XX fixed point and $O(\ell^{-1/2})$ for the XXX fixed point with
non-universal coefficients. For the XX model, where a comparison with exact
diagonalization is possible, the order of the subleading term is confirmed but
we find that the SDRG fails to provide the correct non-universal coefficient.
