Statistical model of ligand substitution
I think the easy way out is to invoke $S_\mathrm m = R \ln \Omega$. If we assume that for a generic complex $\ce{MA_{n}B_{$N-n$}}$,

$$\Omega = {N \choose n} = \frac{N!}{n!(N-n)!} \quad \left[ = {N \choose N-n} \right]$$

and that for the individual molecules $\ce{A}$ and $\ce{B}$, $\Omega = 1$, then the equilibrium constant $K$ for 

$$\ce{MA_{n}B_{$N-n$} + A &lt;=&gt; MA_{n + 1}B_{$N-n-1$} + B}$$

is given by

$$\begin{align}
K &amp;= \exp\left(\frac{-\Delta_\mathrm r G}{RT}\right) \\
&amp;= \exp\left(\frac{\Delta_\mathrm r S}{R}\right) \\
&amp;= \exp\left(\frac{S_\mathrm{m}(\ce{MA_{n + 1}B_{$N-n-1$}}) + S_\mathrm{m}(\ce{B}) - S_\mathrm{m}(\ce{MA_{n}B_{$N-n$}}) - S_\mathrm{m}(\ce{A})}{R}\right) \\
&amp;= \exp[\ln\Omega(\ce{MA_{n + 1}B_{$N-n-1$}}) + \ln\Omega(\ce{B}) - \ln\Omega(\ce{MA_{n}B_{$N-n$}}) - \ln\Omega(\ce{A})] \\
&amp;= \exp\left[\ln\left(\frac{\Omega(\ce{MA_{n + 1}B_{$N-n-1$}})\Omega(\ce{B})}{\Omega(\ce{MA_{n}B_{$N-n$}})\Omega(\ce{A})}\right)\right] \\
&amp;= \frac{\Omega(\ce{MA_{n + 1}B_{$N-n-1$}})\Omega(\ce{B})}{\Omega(\ce{MA_{n}B_{$N-n$}})\Omega(\ce{A})} \\
&amp;= \frac{N!}{(n+1)!(N-n-1)!} \cdot \frac{n!(N-n)!}{N!} \\
&amp;= \frac{N-n}{n+1}
\end{align}$$

The reason for ignoring $\Delta_\mathrm r H$ is because we are only interested in statistical effects, i.e. entropy, and we don't care about the actual stability of the complex or the strength of the M–L bonds. However, the exact justification for assuming this form for $\Omega$ still eludes me. It makes intuitive sense (that there are $N!/(n!(N-n)!)$ ways to arrange $n$ different ligands in $N$ different coordination sites), but I can't convince myself (and don't want to attempt to convince you) that it's entirely rigorous. In particular, I feel like symmetry should play a role here; maybe it is simply that the effects of any symmetry eventually cancel out.