# HEWL S-SAD Merging Statistics

Merging statistics are a useful means to assess data quality in crystallography. However, each statistic has inherent shortcomings. For example, R-merge will appear inflated if the multiplicity is high, and the Pearson correlation coefficients used for $CC_{1/2}$ are very sensitive to outliers. 

Most scaling and merging programs output multiple merging statistics to get around these shortcomings. However, one can imagine that it could also be useful to customize certain parameters, such as how many resolution bins are used. Or, perhaps a better statistic will be developed that is worth implementing. 

In this notebook, ``reciprocalspaceship`` is used to compute half-dataset correlation coefficients ($CC_{1/2}$ and $CC_{anom}$) for a dataset collected from a tetragonal hen egg-white lysozyme (HEWL) crystal at 6550 eV. 
These data are unmerged, but were scaled in AIMLESS. 
They contain sufficient sulfur anomalous signal to determine a solution by the SAD method. 
This illustrates the use of ``rs`` to implement a merging routine, create a custom analysis, and could be useful as a template for other exploratory crystallographic data analyses. 
As an example, we will compare half-dataset correlations computed using both Pearson and Spearman correlation coefficients.


```python
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_context("notebook", font_scale=1.3)
import numpy as np
```


```python
import reciprocalspaceship as rs
```


```python
print(rs.__version__)
```

    0.9.4


---
### Load scaled, unmerged data  

This data has been scaled in AIMLESS. The data includes the image number and the scaled **I** and **SIGI** values. 


```python
hewl = rs.read_mtz("data/HEWL_unmerged.mtz")
```


```python
hewl.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th></th>
      <th></th>
      <th>BATCH</th>
      <th>I</th>
      <th>SIGI</th>
      <th>PARTIAL</th>
    </tr>
    <tr>
      <th>H</th>
      <th>K</th>
      <th>L</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th rowspan="5" valign="top">0</th>
      <th rowspan="5" valign="top">0</th>
      <th>4</th>
      <td>137</td>
      <td>696.5212</td>
      <td>87.83294</td>
      <td>False</td>
    </tr>
    <tr>
      <th>4</th>
      <td>520</td>
      <td>710.6812</td>
      <td>88.107025</td>
      <td>False</td>
    </tr>
    <tr>
      <th>4</th>
      <td>856</td>
      <td>672.05634</td>
      <td>87.75671</td>
      <td>False</td>
    </tr>
    <tr>
      <th>4</th>
      <td>1239</td>
      <td>642.47485</td>
      <td>87.90302</td>
      <td>False</td>
    </tr>
    <tr>
      <th>4</th>
      <td>2160</td>
      <td>655.71783</td>
      <td>87.74394</td>
      <td>False</td>
    </tr>
  </tbody>
</table>
</div>




```python
print(f"Number of observed reflections: {len(hewl)}")
```

    Number of observed reflections: 816804


---
### Merging with Inverse-Variance Weights

Since the input data are unmerged, we will implement the inverse-variance weighting scheme used by AIMLESS to merge the observations. 
The weighted average is a better estimator of the true mean than the raw average, and this weighting scheme corresponds to the maximum likelihood estimator of the true mean if we assume that the observations are normally-distributed about the true mean. 

The merged intensity for each reflection, $I_h$, can be determined from the observed intensities, $I_{h,i}$, and error estimates, $\sigma_{h,i}$, as follows:

\begin{equation}
I_h = \frac{\sum_{i}w_{h,i} I_{h,i}}{\sum_{i} w_{h,i}}
\end{equation}

where the weight for each observation, $w_{h,i}$ is given by:

\begin{equation}
w_{h,i} = \frac{1}{(\sigma_{h,i})^2}
\end{equation}

The updated estimate of the uncertainty, $\sigma_{h}$, is given by:

\begin{equation}
\sigma_{h} = \sqrt{\frac{1}{\sum_{i} w_{h,i}}}
\end{equation}

Let's start by implementing the above equations in a function that will compute the merged $I_h$ and $\sigma_h$. We will use the [Pandas groupby](https://pandas.pydata.org/pandas-docs/stable/user_guide/groupby.html) methods to apply this function across the unique Miller indices in the ``DataSet``. If we group Friedel pairs together, we will refer to the quantity as **IMEAN**, and if we keep the Friedel pairs separate we will refer to the quantities as **I(+)**  and **I(-)**.


```python
def merge(dataset, anomalous=False):
    """
    Merge dataset using inverse-variance weights.
    
    Parameters
    ----------
    dataset : rs.DataSet
        DataSet to be merged containing scaled I and SIGI
    anomalous : bool
        If True, I(+) and I(-) will be reported. If False,
        IMEAN will be reported
        
    Returns
    -------
    rs.DataSet
        Merged DataSet object
    """
    ds = dataset.hkl_to_asu(anomalous=anomalous)
    ds["w"]  = ds['SIGI']**-2
    ds["wI"] = ds["I"] * ds["w"]
    g = ds.groupby(["H", "K", "L"])
    
    result = g[["w", "wI"]].sum()
    result["I"] = result["wI"] / result["w"]
    result["SIGI"] = np.sqrt(1 / result["w"])
    result = result.loc[:, ["I", "SIGI"]]
    result.merged = True
    
    if anomalous:
        result = result.unstack_anomalous()
    
    return result
```

Using `anomalous=False`, this function can be used to compute **IMEAN** and **SIGIMEAN** by including both Friedel pairs:


```python
result1 = merge(hewl, anomalous=False)
```


```python
result1.sample(5)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th></th>
      <th></th>
      <th>I</th>
      <th>SIGI</th>
    </tr>
    <tr>
      <th>H</th>
      <th>K</th>
      <th>L</th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>34</th>
      <th>14</th>
      <th>9</th>
      <td>238.14276</td>
      <td>5.412156</td>
    </tr>
    <tr>
      <th>5</th>
      <th>5</th>
      <th>10</th>
      <td>2619.72</td>
      <td>43.076794</td>
    </tr>
    <tr>
      <th>43</th>
      <th>8</th>
      <th>1</th>
      <td>45.92524</td>
      <td>2.000772</td>
    </tr>
    <tr>
      <th>22</th>
      <th>2</th>
      <th>7</th>
      <td>92.24401</td>
      <td>1.229599</td>
    </tr>
    <tr>
      <th>42</th>
      <th>10</th>
      <th>0</th>
      <td>53.67026</td>
      <td>3.489320</td>
    </tr>
  </tbody>
</table>
</div>



Using `anomalous=True`, this function can be used to compute **I(+)**, **SIGI(+)**, **I(-)**, and **SIGI(-)** by separating Friedel pairs:


```python
result2 = merge(hewl, anomalous=True)
```


```python
result2.sample(5)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th></th>
      <th></th>
      <th>I(+)</th>
      <th>SIGI(+)</th>
      <th>I(-)</th>
      <th>SIGI(-)</th>
    </tr>
    <tr>
      <th>H</th>
      <th>K</th>
      <th>L</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>11</th>
      <th>7</th>
      <th>16</th>
      <td>188.28104</td>
      <td>6.185037</td>
      <td>173.23112</td>
      <td>6.336241</td>
    </tr>
    <tr>
      <th>24</th>
      <th>24</th>
      <th>5</th>
      <td>17.531303</td>
      <td>1.161394</td>
      <td>17.531303</td>
      <td>1.161394</td>
    </tr>
    <tr>
      <th>36</th>
      <th>11</th>
      <th>2</th>
      <td>39.90348</td>
      <td>1.855634</td>
      <td>46.158234</td>
      <td>1.342185</td>
    </tr>
    <tr>
      <th>27</th>
      <th>11</th>
      <th>11</th>
      <td>310.0335</td>
      <td>8.241861</td>
      <td>327.22</td>
      <td>7.975232</td>
    </tr>
    <tr>
      <th>13</th>
      <th>13</th>
      <th>12</th>
      <td>124.78872</td>
      <td>2.328268</td>
      <td>124.78872</td>
      <td>2.328268</td>
    </tr>
  </tbody>
</table>
</div>



A variant of the above function is implemented in `rs.algorithms`, and we will use that implementation in the next section for computing merging statistics. This function computes **IMEAN**, **I(+)**, **I(-)**, and associated uncertainties for each unique Miller index.


```python
result3 = rs.algorithms.merge(hewl)
```


```python
result3.sample(5)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th></th>
      <th></th>
      <th>IMEAN</th>
      <th>SIGIMEAN</th>
      <th>I(+)</th>
      <th>SIGI(+)</th>
      <th>I(-)</th>
      <th>SIGI(-)</th>
      <th>N(+)</th>
      <th>N(-)</th>
    </tr>
    <tr>
      <th>H</th>
      <th>K</th>
      <th>L</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>11</th>
      <th>5</th>
      <th>9</th>
      <td>445.99692</td>
      <td>5.5288916</td>
      <td>446.19846</td>
      <td>7.563336</td>
      <td>445.76562</td>
      <td>8.102551</td>
      <td>64</td>
      <td>56</td>
    </tr>
    <tr>
      <th>30</th>
      <th>8</th>
      <th>3</th>
      <td>81.085785</td>
      <td>1.0904284</td>
      <td>72.2529</td>
      <td>1.5436918</td>
      <td>89.8823</td>
      <td>1.5405104</td>
      <td>60</td>
      <td>64</td>
    </tr>
    <tr>
      <th>15</th>
      <th>1</th>
      <th>5</th>
      <td>2664.3076</td>
      <td>32.183956</td>
      <td>2670.4834</td>
      <td>47.109356</td>
      <td>2658.9023</td>
      <td>44.07226</td>
      <td>56</td>
      <td>64</td>
    </tr>
    <tr>
      <th>8</th>
      <th>5</th>
      <th>0</th>
      <td>917.98047</td>
      <td>15.161337</td>
      <td>917.98047</td>
      <td>15.161337</td>
      <td>917.98047</td>
      <td>15.161337</td>
      <td>64</td>
      <td>64</td>
    </tr>
    <tr>
      <th>34</th>
      <th>27</th>
      <th>0</th>
      <td>50.074245</td>
      <td>4.3588166</td>
      <td>50.074245</td>
      <td>4.3588166</td>
      <td>50.074245</td>
      <td>4.3588166</td>
      <td>4</td>
      <td>4</td>
    </tr>
  </tbody>
</table>
</div>



---
### Merging with 2-fold Cross-Validation  

To compute correlation coefficients we will repeatedly split our data into half-datasets. We will do this by randomly splitting the data using the image number. These half-datasets will be merged independently and used to determine uncertainties in the correlation coefficients. We will first write a method to randomly split our data, and then we will then write a method that automates the sampling and merging of multiple half-datasets. 


```python
def sample_halfdatasets(data):
    """Randomly split DataSet into two equal halves by BATCH"""
    batch = data.BATCH.unique().to_numpy(dtype=int)
    np.random.shuffle(batch)
    halfbatch1, halfbatch2 = np.array_split(batch, 2)
    half1 = data.loc[data.BATCH.isin(halfbatch1)]
    half2 = data.loc[data.BATCH.isin(halfbatch2)]
    return half1, half2
```


```python
def merge_dataset(dataset, nsamples):
    """
    Merge DataSet using inverse-variance weighting scheme. This represents the 
    maximum-likelihood estimator of the mean of the observed intensities assuming 
    they are independent and normally distributed with the same mean. 
    
    Sample means across half-datasets can be used to compute the merging statistics CC1/2 and CCanom.
    """
    dataset = dataset.copy()
    samples = []
    for n in range(nsamples):
        half1, half2 = sample_halfdatasets(dataset)
        mergedhalf1 = rs.algorithms.merge(half1)
        mergedhalf2 = rs.algorithms.merge(half2)
        result = mergedhalf1.merge(mergedhalf2, left_index=True, right_index=True, suffixes=(1, 2))
        result["sample"] = n
        samples.append(result)
    return rs.concat(samples).sort_index()
```

---
### Merge HEWL data

We will now merge the HEWL data, repeatedly sampling across half-datasets in order to assess the distribution of correlation coefficients.


```python
# This cell takes a few minutes with nsamples=15
merged = merge_dataset(hewl, 15)
```


```python
merged
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th></th>
      <th></th>
      <th>IMEAN1</th>
      <th>SIGIMEAN1</th>
      <th>I(+)1</th>
      <th>SIGI(+)1</th>
      <th>I(-)1</th>
      <th>SIGI(-)1</th>
      <th>N(+)1</th>
      <th>N(-)1</th>
      <th>IMEAN2</th>
      <th>SIGIMEAN2</th>
      <th>I(+)2</th>
      <th>SIGI(+)2</th>
      <th>I(-)2</th>
      <th>SIGI(-)2</th>
      <th>N(+)2</th>
      <th>N(-)2</th>
      <th>sample</th>
    </tr>
    <tr>
      <th>H</th>
      <th>K</th>
      <th>L</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th rowspan="5" valign="top">0</th>
      <th rowspan="5" valign="top">0</th>
      <th>4</th>
      <td>662.71674</td>
      <td>27.763823</td>
      <td>662.71674</td>
      <td>27.76382</td>
      <td>662.71674</td>
      <td>27.76382</td>
      <td>10</td>
      <td>10</td>
      <td>658.93616</td>
      <td>35.859802</td>
      <td>658.93616</td>
      <td>35.859802</td>
      <td>658.93616</td>
      <td>35.859802</td>
      <td>6</td>
      <td>6</td>
      <td>0</td>
    </tr>
    <tr>
      <th>4</th>
      <td>673.3035</td>
      <td>39.26171</td>
      <td>673.3035</td>
      <td>39.26171</td>
      <td>673.3035</td>
      <td>39.26171</td>
      <td>5</td>
      <td>5</td>
      <td>655.83997</td>
      <td>26.47925</td>
      <td>655.83997</td>
      <td>26.47925</td>
      <td>655.83997</td>
      <td>26.47925</td>
      <td>11</td>
      <td>11</td>
      <td>1</td>
    </tr>
    <tr>
      <th>4</th>
      <td>652.56683</td>
      <td>29.268377</td>
      <td>652.56683</td>
      <td>29.268377</td>
      <td>652.56683</td>
      <td>29.268377</td>
      <td>9</td>
      <td>9</td>
      <td>672.53235</td>
      <td>33.193497</td>
      <td>672.53235</td>
      <td>33.193497</td>
      <td>672.53235</td>
      <td>33.193497</td>
      <td>7</td>
      <td>7</td>
      <td>2</td>
    </tr>
    <tr>
      <th>4</th>
      <td>651.18304</td>
      <td>29.259443</td>
      <td>651.18304</td>
      <td>29.259443</td>
      <td>651.18304</td>
      <td>29.259443</td>
      <td>9</td>
      <td>9</td>
      <td>674.3303</td>
      <td>33.20654</td>
      <td>674.3303</td>
      <td>33.20654</td>
      <td>674.3303</td>
      <td>33.20654</td>
      <td>7</td>
      <td>7</td>
      <td>3</td>
    </tr>
    <tr>
      <th>4</th>
      <td>664.65784</td>
      <td>29.259375</td>
      <td>664.65784</td>
      <td>29.259375</td>
      <td>664.65784</td>
      <td>29.259375</td>
      <td>9</td>
      <td>9</td>
      <td>656.9748</td>
      <td>33.206642</td>
      <td>656.9748</td>
      <td>33.206642</td>
      <td>656.9748</td>
      <td>33.206642</td>
      <td>7</td>
      <td>7</td>
      <td>4</td>
    </tr>
    <tr>
      <th>...</th>
      <th>...</th>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th rowspan="5" valign="top">45</th>
      <th rowspan="5" valign="top">10</th>
      <th>2</th>
      <td>7.669444</td>
      <td>6.6297607</td>
      <td>7.669444</td>
      <td>6.6297607</td>
      <td>NaN</td>
      <td>NaN</td>
      <td>1</td>
      <td>0</td>
      <td>22.89468</td>
      <td>4.1084785</td>
      <td>22.89468</td>
      <td>4.1084785</td>
      <td>NaN</td>
      <td>NaN</td>
      <td>3</td>
      <td>0</td>
      <td>10</td>
    </tr>
    <tr>
      <th>2</th>
      <td>14.465271</td>
      <td>4.4560146</td>
      <td>14.465271</td>
      <td>4.4560146</td>
      <td>NaN</td>
      <td>NaN</td>
      <td>2</td>
      <td>0</td>
      <td>25.364758</td>
      <td>5.6226034</td>
      <td>25.364758</td>
      <td>5.6226034</td>
      <td>NaN</td>
      <td>NaN</td>
      <td>2</td>
      <td>0</td>
      <td>11</td>
    </tr>
    <tr>
      <th>2</th>
      <td>18.59031</td>
      <td>8.728156</td>
      <td>18.59031</td>
      <td>8.728156</td>
      <td>NaN</td>
      <td>NaN</td>
      <td>1</td>
      <td>0</td>
      <td>18.685297</td>
      <td>3.8105915</td>
      <td>18.685297</td>
      <td>3.8105915</td>
      <td>NaN</td>
      <td>NaN</td>
      <td>3</td>
      <td>0</td>
      <td>12</td>
    </tr>
    <tr>
      <th>2</th>
      <td>14.465271</td>
      <td>4.4560146</td>
      <td>14.465271</td>
      <td>4.4560146</td>
      <td>NaN</td>
      <td>NaN</td>
      <td>2</td>
      <td>0</td>
      <td>25.364758</td>
      <td>5.6226034</td>
      <td>25.364758</td>
      <td>5.6226034</td>
      <td>NaN</td>
      <td>NaN</td>
      <td>2</td>
      <td>0</td>
      <td>13</td>
    </tr>
    <tr>
      <th>2</th>
      <td>15.318142</td>
      <td>3.9687188</td>
      <td>15.318142</td>
      <td>3.9687188</td>
      <td>NaN</td>
      <td>NaN</td>
      <td>3</td>
      <td>0</td>
      <td>30.170218</td>
      <td>7.3511124</td>
      <td>30.170218</td>
      <td>7.3511124</td>
      <td>NaN</td>
      <td>NaN</td>
      <td>1</td>
      <td>0</td>
      <td>14</td>
    </tr>
  </tbody>
</table>
<p>187328 rows × 17 columns</p>
</div>



---
### Compute $CC_{1/2}$ and $CC_{anom}$  

We will first assign each reflection to a resolution bin and then we will compute the correlation coefficients


```python
merged, labels = merged.assign_resolution_bins(bins=15)
```


```python
groupby1 = merged.groupby(["sample", "bin"])[["IMEAN1", "IMEAN2"]]
pearson1 = groupby1.corr(method="pearson").unstack().loc[:, ("IMEAN1", "IMEAN2")]
pearson1.name = "Pearson"
spearman1 = groupby1.corr(method="spearman").unstack().loc[:, ("IMEAN1", "IMEAN2")]
spearman1.name = "Spearman"
results1 = rs.concat([pearson1, spearman1], axis=1)
results1 = results1.groupby("bin").agg(["mean", "std"])
```


```python
results1.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead tr th {
        text-align: left;
    }

    .dataframe thead tr:last-of-type th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr>
      <th></th>
      <th colspan="2" halign="left">Pearson</th>
      <th colspan="2" halign="left">Spearman</th>
    </tr>
    <tr>
      <th></th>
      <th>mean</th>
      <th>std</th>
      <th>mean</th>
      <th>std</th>
    </tr>
    <tr>
      <th>bin</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.997269</td>
      <td>0.000416</td>
      <td>0.998856</td>
      <td>0.000104</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0.998025</td>
      <td>0.000334</td>
      <td>0.999361</td>
      <td>0.000059</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0.998973</td>
      <td>0.000313</td>
      <td>0.999649</td>
      <td>0.000054</td>
    </tr>
    <tr>
      <th>3</th>
      <td>0.999317</td>
      <td>0.000113</td>
      <td>0.999661</td>
      <td>0.000027</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0.999386</td>
      <td>0.000075</td>
      <td>0.999680</td>
      <td>0.000015</td>
    </tr>
  </tbody>
</table>
</div>




```python
plt.figure(figsize=(8, 4))
plt.errorbar(results1.index, results1[("Pearson", "mean")], 
             yerr=results1[("Pearson", "std")], 
             color="#1b9e77", 
             label=r"$CC_{1/2}$ (Pearson)")
plt.errorbar(results1.index, results1[("Spearman", "mean")], 
             yerr=results1[("Spearman", "std")], 
             color="#d95f02", 
             label=r"$CC_{1/2}$ (Spearman)")
plt.xticks(results1.index, labels, rotation=45, ha='right', rotation_mode='anchor')
plt.ylabel("Correlation Coefficient")
plt.xlabel(r"Resolution Bin ($\AA$)")
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid(axis="y", linestyle='--')
plt.show()
```

It is important to note the scale on the y-axis -- this dataset is edge-limited, and as such the $CC_{1/2}$ is very high across all resolution bins. The Spearman CC appears higher across all resolution bins except at high resolution, and overall has a lower standard deviation among samples.
This is consistent with our expectation that Spearman CCs are a more robust estimator of correlation than Pearson CCs.

Let's now repeat this for the anomalous data, computing $CC_{anom}$:


```python
merged["ANOM1"] = merged["I(+)1"] - merged["I(-)1"]
merged["ANOM2"] = merged["I(+)2"] - merged["I(-)2"]

# Similar to CChalf, but we will only look at acentric reflections
groupby2 = merged.acentrics.groupby(["sample", "bin"])[["ANOM1", "ANOM2"]]
pearson2 = groupby2.corr(method="pearson").unstack().loc[:, ("ANOM1", "ANOM2")]
pearson2.name = "Pearson"
spearman2 = groupby2.corr(method="spearman").unstack().loc[:, ("ANOM1", "ANOM2")]
spearman2.name = "Spearman"
results2 = rs.concat([pearson2, spearman2], axis=1)
results2 = results2.groupby("bin").agg(["mean", "std"])
```


```python
results2.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead tr th {
        text-align: left;
    }

    .dataframe thead tr:last-of-type th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr>
      <th></th>
      <th colspan="2" halign="left">Pearson</th>
      <th colspan="2" halign="left">Spearman</th>
    </tr>
    <tr>
      <th></th>
      <th>mean</th>
      <th>std</th>
      <th>mean</th>
      <th>std</th>
    </tr>
    <tr>
      <th>bin</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.405502</td>
      <td>0.048131</td>
      <td>0.492206</td>
      <td>0.029448</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0.255115</td>
      <td>0.047082</td>
      <td>0.321767</td>
      <td>0.030535</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0.343341</td>
      <td>0.052632</td>
      <td>0.432287</td>
      <td>0.020717</td>
    </tr>
    <tr>
      <th>3</th>
      <td>0.470596</td>
      <td>0.051916</td>
      <td>0.551680</td>
      <td>0.025049</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0.586280</td>
      <td>0.043609</td>
      <td>0.615225</td>
      <td>0.024101</td>
    </tr>
  </tbody>
</table>
</div>




```python
plt.figure(figsize=(8, 4))
plt.errorbar(results2.index, results2[("Pearson", "mean")], 
             yerr=results2[("Pearson", "std")], 
             color="#1b9e77", 
             label=r"$CC_{anom}$ (Pearson)")
plt.errorbar(results2.index, results2[("Spearman", "mean")],
             yerr=results2[("Spearman", "std")], 
             color="#d95f02", 
             label=r"$CC_{anom}$ (Spearman)")
plt.xticks(results2.index, labels, rotation=45, ha='right', rotation_mode='anchor')
plt.ylabel("Correlation Coefficient")
plt.xlabel(r"Resolution Bin ($\AA$)")
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid(axis="y", linestyle='--')
plt.show()
```

There is significant anomalous signal across all but the highest resolution bins. The Spearman CCs have smaller error bars than the corresponding Pearson CCs and the Spearman CCs are also higher in most bins. These differences highlight the influence of outlier measurements on the different correlation coefficients.

---
### Summary  

We have used `reciprocalspaceship` to merge a dataset using inverse-variance weights. As part of this analysis, we performed repeated 2-fold cross-validation to compute Pearson and Spearman correlation coefficients and associated uncertainties. This relatively simple procedure to obtain uncertainty estimates for correlation coefficients is seldom done when analyzing merging quality. However, we can see that the standard deviation for computed correlation coefficients can be nearly $\pm0.1$ for quantities such as $CC_{anom}$. This is worth keeping in mind when analyzing SAD experiments because this dataset is very high quality (edge-limited) when one considers the $CC_{1/2}$. Putting these two correlation coefficients on the same axes emphasizes this point:


```python
plt.figure(figsize=(9, 6))
plt.errorbar(results1.index, results1[("Pearson", "mean")], 
             yerr=results1[("Pearson", "std")], 
             color='#1b9e77', 
             label=r"$CC_{1/2}$ (Pearson)")
plt.errorbar(results1.index, results1[("Spearman", "mean")], 
             yerr=results1[("Spearman", "std")], 
             color='#d95f02', 
             label=r"$CC_{1/2}$ (Spearman)")
plt.errorbar(results2.index, results2[("Pearson", "mean")], 
             yerr=results2[("Pearson", "std")], 
             color='#1b9e77', 
             linestyle="--", 
             label=r"$CC_{anom}$ (Pearson)")
plt.errorbar(results2.index, results2[("Spearman", "mean")], 
             yerr=results2[("Spearman", "std")], 
             color='#d95f02', 
             linestyle="--", 
             label=r"$CC_{anom}$ (Spearman)")
plt.xticks(results1.index, labels, rotation=45, ha='right', rotation_mode='anchor')
plt.ylabel(r"Correlation Coefficient")
plt.xlabel(r"Resolution Bin ($\AA$)")
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid(axis="y", linestyle='--')
plt.tight_layout()
```

Even a simple change to the procedure for computing merging statistics, such as substituting a Spearman correlation coefficient for a Pearson one, can alter the apparent quality of a dataset. By lowering the barrier to implementing new analyses, we hope that `reciprocalspaceship` can encourage the development of more robust indicators of crystallographic data quality.

d92bfadf94d77eeb25587a9a07ec27f35a073750
ipynb
Jupyter Notebook
docs/examples/2_mergingstats.ipynb
marinegor/reciprocalspaceship
6a5ef478cea85183d178c35963016bdbe2c582bf
MIT
docs/examples/2_mergingstats.ipynb
marinegor/reciprocalspaceship
6a5ef478cea85183d178c35963016bdbe2c582bf
MIT
docs/examples/2_mergingstats.ipynb
marinegor/reciprocalspaceship
6a5ef478cea85183d178c35963016bdbe2c582bf
MIT
Qwen/Qwen-72B
1. YES
2. YES