  We present a study of unprecedented statistical power regarding the
halo-to-halo variance of dark matter substructure. Using a combination of
N-body simulations and a semi-analytical model, we investigate the variance in
subhalo mass fractions and subhalo occupation numbers, with an emphasis on how
these statistics scale with halo formation time. We demonstrate that the
subhalo mass fraction, f_sub, is mainly a function of halo formation time, with
earlier forming haloes having less substructure. At fixed formation redshift,
the average f_sub is virtually independent of halo mass, and the mass
dependence of f_sub is therefore mainly a manifestation of more massive haloes
assembling later. We compare observational constraints on f_sub from
gravitational lensing to our model predictions and simulation results. Although
the inferred f_sub are substantially higher than the median LCDM predictions,
they fall within the 95th percentile due to halo-to-halo variance. We show that
while the halo occupation distribution of subhaloes, P(N|M), is
super-Poissonian for large <N>, a well established result, it becomes
sub-Poissonian for <N> < 2. Ignoring the non-Poissonity results in systematic
errors of the clustering of galaxies of a few percent, and with a complicated
scale- and luminosity-dependence. Earlier-formed haloes have P(N|M) closer to a
Poisson distribution, suggesting that the dynamical evolution of subhaloes
drives the statistics towards Poissonian. Contrary to a recent claim, the
non-Poissonity of subhalo occupation statistics does not vanish by selecting
haloes with fixed mass and fixed formation redshift. Finally, we use subhalo
occupation statistics to put loose constraints on the mass and formation
redshift of the Milky Way halo. Using observational constraints on the V_max of
the most massive satellites, we infer that 0.25<M_vir/10^12M_sun/h<1.4 and
0.1<z_f<1.4 at 90% confidence.
