library(admix)
The clustering of populations following admixture models is, for now, based on the K-sample test theory. Consider \(K\) samples. For \(i=1,...,K\), sample \(X^{(i)} = (X_1^{(i)}, ..., X_{n_i}^{(i)})\) follows \[L_i(x) = p_i F_i(x) + (1-p_i) G_i, \qquad x \in \mathbb{R}.\]
We still use IBM approach to perform pairwise hypothesis testing. The idea is to adapt the K-sample test procedure to obtain a data-driven method that cluster the \(K\) populations into \(N\) subgroups, characterized by a common unknown mixture component. The advantages of such an approach is twofold:
This clustering technique thus allows to cluster unobserved subpopulations instead of individuals. We call this algorithm the K-sample 2-component mixture clustering (K2MC).
We now detail the steps of the algorithm.
If $H_0$ is not rejected then $S_1 = \{x,y\}$,\\
Else $S_1 = \{x\}$, $S_{c+1} = \{y\}$ and then $c=c+1$.
Select $u={\rm argmin}\{d(i,j); i\in S_c, j\in S\setminus \bigcup_{k=1}^c S_k\}$;
Test $H_0$ the simultaneous equality of all the $f_j$, $j\in S_c$ :\\
If $H_0$ not rejected, then put $S_c=S_c\bigcup \{u\}$;\\
Else $S_{c+1} = \{u\}$ and $c = c+1$.
We present a case study with 5 populations to cluster on \(\mathbb{R}^+\), with Gamma-Exponential and Gamma-Gamma mixtures.
## Simulate data (chosen parameters indicate 3 clusters (populations (1,3), (2,5) and 4)!):
<- list(f1 = "gamma", g1 = "exp",
list.comp f2 = "gamma", g2 = "exp",
f3 = "gamma", g3 = "gamma",
f4 = "exp", g4 = "exp",
f5 = "gamma", g5 = "exp")
<- list(f1 = list(shape = 16, rate = 4), g1 = list(rate = 1/3.5),
list.param f2 = list(shape = 14, rate = 2), g2 = list(rate = 1/5),
f3 = list(shape = 16, rate = 4), g3 = list(shape = 12, rate = 2),
f4 = list(rate = 1/2), g4 = list(rate = 1/7),
f5 = list(shape = 14, rate = 2), g5 = list(rate = 1/6))
<- rsimmix(n=6000, unknownComp_weight=0.8, comp.dist = list(list.comp$f1,list.comp$g1),
A.sim comp.param = list(list.param$f1, list.param$g1))$mixt.data
<- rsimmix(n=6000, unknownComp_weight=0.5, comp.dist = list(list.comp$f2,list.comp$g2),
B.sim comp.param = list(list.param$f2, list.param$g2))$mixt.data
<- rsimmix(n=6000, unknownComp_weight=0.65, comp.dist = list(list.comp$f3,list.comp$g3),
C.sim comp.param = list(list.param$f3, list.param$g3))$mixt.data
<- rsimmix(n=6000, unknownComp_weight=0.75, comp.dist = list(list.comp$f4,list.comp$g4),
D.sim comp.param = list(list.param$f4, list.param$g4))$mixt.data
<- rsimmix(n=6000, unknownComp_weight=0.55, comp.dist = list(list.comp$f5,list.comp$g5),
E.sim comp.param = list(list.param$f5, list.param$g5))$mixt.data
## Look for the clusters:
<- list(f1 = NULL, g1 = "exp",
list.comp f2 = NULL, g2 = "exp",
f3 = NULL, g3 = "gamma",
f4 = NULL, g4 = "exp",
f5 = NULL, g5 = "exp")
<- list(f1 = NULL, g1 = list(rate = 1/3.5),
list.param f2 = NULL, g2 = list(rate = 1/5),
f3 = NULL, g3 = list(shape = 12, rate = 2),
f4 = NULL, g4 = list(rate = 1/7),
f5 = NULL, g5 = list(rate = 1/6))
<- admix_clustering(samples = list(A.sim,B.sim,C.sim,D.sim,E.sim), n_sim_tab = 10,
clusters comp.dist=list.comp, comp.param=list.param, parallel=FALSE, n_cpu=2)
clusters#> $n_clust
#> [1] 3
#>
#> $clusters
#> [,1] [,2] [,3] [,4] [,5]
#> Id_sample 1 2 3 4 5
#> Id_cluster 3 1 3 2 1
#>
#> $discrepancy_matrix
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.0000000 106.2503840 0.4559752 17.198006 212.5783544
#> [2,] 106.2503840 0.0000000 68.5687703 5.998583 0.2265777
#> [3,] 0.4559752 68.5687703 0.0000000 227.404200 74.4947230
#> [4,] 17.1980060 5.9985831 227.4041998 0.000000 9.5940432
#> [5,] 212.5783544 0.2265777 74.4947230 9.594043 0.0000000