Revise docs around nan values

Author: apragsdale

docs/stats.md

@@ -836,26 +836,11 @@ or tuple will produce a single two-dimensional matrix, while list of these
 will produce a three-dimensional array, with the outer dimension of length
 equal to the length of the list.
 
-For concreteness, we would expect the following dimensions with the specified
-`sample_sets` and `indexes` arguments.
+For example, to compute the {math}`r^2` LD matrix over a subset of samples in
+the tree sequence (such as sample nodes 0 through 7), we would specify the
+samples as follows:
 
 ```
-# one-way
-ts.ld_matrix(sample_sets=None) -> 2 dimensions
-ts.ld_matrix(sample_sets=[0, 1, 2, 3]) -> 2 dimensions
-ts.ld_matrix(sample_sets=[[0, 1, 2, 3]]) -> 3 dimensions
-# two-way
-ts.ld_matrix(sample_sets=[[0, 1, 2, 3], [4, 5, 6, 7]], indexes=(0, 1)) -> 2 dimensions
-ts.ld_matrix(sample_sets=[[0, 1, 2, 3], [4, 5, 6, 7]], indexes=[(0, 1)]) -> 3 dimensions
-```
-
-By computing the LD statistic for a subset of samples in the tree sequence, it
-is possible that statistics are computed over pairs of sites that are not
-variable within the sample set. This can result in zero- or nan-valued statistics
-as outputs. To avoid this, the tree sequence may first be simplified to the
-sample subset, which will remove sites that are invariant within that subset:
-
-```{code-cell} ipython3
 ts = msprime.sim_ancestry(
     20,
     population_size=10000,
@@ -864,32 +849,38 @@ ts = msprime.sim_ancestry(
     random_seed=12)
 ts = msprime.sim_mutations(ts, rate=2e-8, random_seed=12)
 
-ld_all = ts.ld_matrix(mode="site")
-print(ld_all)
+ld = ts.ld_matrix(mode=“site”, sample_sets=range(8))
+print(ld)
 ```
 
-Considering a subset of samples:
+For concreteness, we would expect the following dimensions with the specified
+`sample_sets` and `indexes` arguments.
 
-```{code-cell} ipython3
-ld_sub = ts.ld_matrix(mode="site", sample_sets=range(8))
-print(ld_sub)
 ```
-
-Simplifying and computing the LD matrix avoids computing {math}`r^2` for
-invariant sites in the subsample:
-
-```{code-cell} ipython3
-ts_simp = ts.simplify(samples=range(8), filter_sites=True)
-ld_sub = ts_simp.ld_matrix(mode="site", sample_sets=range(8))
-print(ld_sub)
+# one-way
+ts.ld_matrix(sample_sets=None) -> 2 dimensions
+ts.ld_matrix(sample_sets=[0, 1, 2, 3]) -> 2 dimensions
+ts.ld_matrix(sample_sets=[[0, 1, 2, 3]]) -> 3 dimensions
+# two-way
+ts.ld_matrix(sample_sets=[[0, 1, 2, 3], [4, 5, 6, 7]], indexes=(0, 1)) -> 2 dimensions
+ts.ld_matrix(sample_sets=[[0, 1, 2, 3], [4, 5, 6, 7]], indexes=[(0, 1)]) -> 3 dimensions
 ```
 
-Any remaining sites that produce {math}`r^2` with a nan values are those with
-a mutation that occurs above the root of the tree, so that all samples in the
-subset carry the derived allele (so that its allele frequency is 1). Also note
-that under certain mutation models, back mutations may create sites that are
-invariant in a sample, and this will also result in nan values observed in the
-LD matrix.
+#### Why are there nan values in the LD matrix?
+
+For some statistics, it is possible to observe nan entries in the LD matrix,
+which can be surprising or numerically impact downstream analyses. A nan entry
+may occur when computing a ratio statistic (such as {math}`r` or {math}`r^2`)
+with a denominator of zero, indicating that one or both sites in the pair are
+monomorphic. This can happen for a number of reasons:
+
+- The mutation models allows for reversible mutations, so a back mutation at
+  a site has resulted in a single allele despite multiple mutations in the
+  history of the sample.
+- LD is computed for a subsample of individuals, and some sites are not
+  variable among the sample nodes in the subsample.
+- A mutation exists above the root of the local tree, so that all samples carry
+  the mutation, and one or more sites are not variable.
 
 (sec_stats_two_locus_sample_one_way_stats)=