Python或R函数根据所有测试变量的相似度对多变量化学数据进行分组

Question

我对编程还是个新手，但我使用GIS软件来绘制样本的空间分布（样本数量为2000），这些样本有多个变量（总共有16个元素）进行测试。我想用Python、R或者其他任何语言，来根据所有测试元素的相对相似度（%）来统计出化学相关的分组。

举个例子，使用值的正负10%：

          Al    Na    K     Si
Sample1   1     2     3     4
Sample2   1.5   2.5   3.5   4.5 
Sample3   1.05  2.1   3.1   1
Sample4   0.95  1.9   3.1   3.8

期望的结果：

• 样本1会和样本4在同一组（因为它们的四个值都在10%的范围内）。
• 样本2会单独成组。
• 样本3也会单独成组，因为虽然前面三个元素和样本1的值在10%范围内，但最后一个元素（Si）不在这个范围内，所以它不被分到同一组。

我看过一些讨论，但似乎没有一个能把所有变量都考虑在内，只是根据一个变量进行分组。

Answer 1

利用递归的方式：

R

get_groups <- function(x, a = list()){
 r <-  colSums(abs(t(x) - (y<-unlist(x[1,])))/y > 0.1) == 0
 if(all(r)) c(a, list(x))
 else get_groups(x[!r, , drop = FALSE], c(a, list(x[r,, drop = FALSE])))
}

res <- get_groups(df)
res
[[1]]
          Al  Na   K  Si
Sample1 1.00 2.0 3.0 4.0
Sample4 0.95 1.9 3.1 3.8

[[2]]
         Al  Na   K  Si
Sample2 1.5 2.5 3.5 4.5

[[3]]
          Al  Na   K Si
Sample3 1.05 2.1 3.1  1

---

Python

def get_groups(x, a = None):
  if a is None: a = []
  r = (x - (y := x.iloc[0])).abs().div(y).gt(0.1).sum(1) == 0
  if r.all():  return a + [x]
  return get_groups(x[~r], a + [x[r]])

res = get_groups(df)

res[0]
           Al   Na    K   Si
Sample1  1.00  2.0  3.0  4.0
Sample4  0.95  1.9  3.1  3.8

res[1]
          Al   Na    K   Si
Sample2  1.5  2.5  3.5  4.5

res[2]
           Al   Na    K   Si
Sample3  1.05  2.1  3.1  1.0

Answer 2

因为你想要逐个比较数值，所以像DBSCAN或K-means这样的标准聚类方法不适用，因为它们是通过计算点之间的整体距离来工作的。

为了解决这个问题，你可以通过逐个比较每个项目与其他所有项目，找到每个数值之间的最大差异（以百分比表示）。函数max_pct_difference就是用来计算最大百分比差异的。这些X与Y之间的差异会形成你感兴趣的“距离”矩阵的下半部分。

通过提前计算一个距离矩阵，你可以在聚类算法（比如DBSCAN）中使用它。这样可以把相似的项目聚在一起，只要它们之间的距离小于一个参数值：epsilon。不过你要注意，并不是所有在同一组里的项目之间的距离都小于10%。这是因为可能存在这样的情况：A与B的距离在10%以内，B与C的距离也在10%以内，但A与C之间的距离却不在10%以内：

A   1.0  1.0  1.0  2.1
B   1.0  1.0  1.0  2.3
C   1.0  1.0  1.0  2.4

这里有一些代码可以帮助你入门：

import pandas as pd
from itertools import combinations_with_replacement

df = pd.DataFrame(
 {'Al': {'Sample1': 1.0, 'Sample2': 1.5, 'Sample3': 1.05, 'Sample4': 0.95},
  'Na': {'Sample1': 2.0, 'Sample2': 2.5, 'Sample3': 2.1, 'Sample4': 1.9},
  'K': {'Sample1': 3.0, 'Sample2': 3.5, 'Sample3': 3.1, 'Sample4': 3.1},
  'Si': {'Sample1': 4.0, 'Sample2': 4.5, 'Sample3': 1.0, 'Sample4': 3.8}}
)

def max_pct_difference(s1: tuple, s2: tuple):
    max_diff = 0.0
    for x, y in zip(s1[1:], s2[1:]):
        # handle division by zero
        if x == 0 and y == 0:
            continue
        diff = abs(x - y) / max(x, y)
        max_diff = max(diff, max_diff)
    return diff


pairs = []
for s1, s2 in combinations_with_replacement(df.itertuples(), 2):
    max_diff = max_difference(s1, s2)
    # append s1, s2 and then s2, s1 to make the matrix symmetric
    pairs.append((s1.Index, s2.Index, max_diff))
    if s1.Index != s2.Index:
        pairs.append((s2.Index, s1.Index, max_diff))


# this creates the distance matrix
dist_df = pd.DataFrame(pairs, columns=['X', 'Y', 'd']).set_index(['X','Y']).unstack(level=1)

在创建了距离矩阵后，你可以使用DBSCAN来找到在距离阈值epsilon内的聚类/组。对于参数，你可以设置eps=0.1，并将每个聚类的最小点数设置为2，这样可以保留所有至少有2个相近点的组。对于度量标准，我们将使用"precomputed"，因为我们已经以自定义方式计算了距离。

from sklearn.cluster import DBSCAN

dbs = DBSCAN(eps=0.1, min_samples=2, metric='precomputed')

# get the group label for each item in the distance data frame
groups = dbs.fit_predict(dist_df)

# assign the group label to the corresponding sample
df['Group'] = groups

print(df)
# prints:
#            Al   Na    K   Si  Group
# Sample1  1.00  2.0  3.0  4.0      0
# Sample2  1.50  2.5  3.5  4.5     -1
# Sample3  1.05  2.1  3.1  1.0     -1
# Sample4  0.95  1.9  3.1  3.8      0

在这个例子中只有一个组，组0，它包含样本1和4。标记为组-1的行没有被分到任何组。

你可以通过调整eps或者改变距离的计算方式来控制组的归属。