Lab 9 - cluster_validation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d0fa5561",
   "metadata": {},
   "source": [
    "# Lab 9 - Cluster Validation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63fe6ef3",
   "metadata": {},
   "source": [
    "In this last notebook of the Winter 2025 series, we consider cluster validation measures. These can be classified as *internal* validation measures, such as the measures of fit we already considered, as well as *external* validation measures, which provide a way to compare cluster solutions to a given (known, or assumed known) reference.\n",
    "\n",
    "In addition to measures of fit, we also consider indicators of the balance of cluster solutions, i.e., the evenness of the number of observations in each cluster. Such measures include *entropy* and *Simpson's index*. Other measures, introduced in Anselin (2024) are based on the *spatial* properties of the clusters. These include the *join count ratio*, an indicator of how many *neighbors* of each observation in a cluster are also members of the cluster. For a spatially compact cluster solution, this measure should equal one (except for boundary effects). For non spatially constrained clusters, it indicates how closely it approximates a spatial solution.\n",
    "\n",
    "For spatially constrained cluster solutions, compactness is a key characteristic. This can be quantified by means of the *isoperimeter quotient (IPQ)*, the ratio of the area of a cluster shape to that of a circle with equal perimeter. A final measure of compactness introduced in Anselin (2024) is the *diameter* of the unweighted graph representation of the spatial weights matrix. To obtain a relative measure, the diameter is rescaled by the number of observations in the cluster. The latter measures are only applicable to spatially constrained clusters.\n",
    "\n",
    "In addition, we also consider two classic indicators of external validity, i.e., the *Adjusted Rand Index (ARI)*, based on counting pairs, and the *Normalized Information Distance (NID)*, derived from measures of entropy.\n",
    "\n",
    "The material is part of the Spatial Cluster Analysis course taught at the University of Chicago in the Winter Quarter of 2025.\n",
    "\n",
    "Prepared by: Luc Anselin (anselin@uchicago.edu) and Pedro Amaral (pedroamaral@uchicago.edu)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06c7b3cb",
   "metadata": {},
   "source": [
    "## Preliminaries\n",
    "The empirical illustration presented here is based on the material in Chapter 12 of the Spatial Cluster book. However, it is adapted to reflect the max-p solution of p=13 obtained with `pygeoda`, which is different from the p=12 listed in the Spatial Clustering book."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d76e4c1e",
   "metadata": {},
   "source": [
    "### Required packages"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5be87c72",
   "metadata": {},
   "source": [
    "The conda enviroment used for this exercise was created from a yml file with the same specification as in the previous notebooks:"
   ]
  },
  {
   "cell_type": "raw",
   "id": "86259abc",
   "metadata": {},
   "source": [
    "name: sp_clusters\n",
    "channels:\n",
    "  - conda-forge\n",
    "  - defaults\n",
    "dependencies:\n",
    "  - python=3.12\n",
    "  - numpy \n",
    "  - nb_conda\n",
    "  - notebook\n",
    "  - geopandas\n",
    "  - scikit-learn\n",
    "  - matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb97d2ab",
   "metadata": {},
   "source": [
    "In addition to the usual `numpy`, `pandas` and `geopandas`, we also import several specialized packages from scikit-learn and `pygeoda` to carry out the cluster analysis. Specifically, to carry out variable standardization we import `StandardScaler` from `sklearn.preprocessing`. The specific clustering methods are `AgglomerativeClustering` and `KMeans` from `sklearn.cluster`. The other clustering solutions are obtained with `pygeoda`. The external validation measures are contained in `sklearn.metrics.adjusted_rand_score` and `sklearn.metrics.adjusted_mutual_info_score`.\n",
    "\n",
    "The new internal validation measures are based on `pygeoda.spatial_validation`. Several helper functions contained in the `spatial_cluster_course` module extract the relevant information and present it as a pandas data frame: `cluster_fragmentation`, `cluster_joincount`, `cluster_compactness` and `cluster_diameter`.\n",
    "\n",
    "\n",
    "For the spatially constrained clustering methods, make sure the latest version of `pygeoda` is installed by using `pip install -U pygeoda`.\n",
    "\n",
    "In addition to the helper functions mentioned above, we also import:\n",
    "\n",
    "- `cluster_stats`\n",
    "- `cluster_fit`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e49b893f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from sklearn.cluster import AgglomerativeClustering, KMeans\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.metrics import adjusted_rand_score, adjusted_mutual_info_score\n",
    "\n",
    "from spatial_cluster_course import cluster_stats, cluster_fit, cluster_fragmentation\n",
    "from spatial_cluster_course import cluster_joincount, cluster_compactness, cluster_diameter\n",
    "\n",
    "import pygeoda"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "737ad673",
   "metadata": {},
   "source": [
    "### Load data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51281568",
   "metadata": {},
   "source": [
    "For this final exercise, we will again use a data set on Zika and Microcephaly infections and socio-economic profiles for 2013-2016 in municipalities in the State of Ceará, Brazil. This is also a GeoDa sample data set. Detailed source and information available at https://geodacenter.github.io/data-and-lab/Ceara-Zika/\n",
    "\n",
    "The following files will be used:\n",
    "- **ceara.shp,shx,dbf,prj**: shape file (four files) for 184 municipalities\n",
    "\n",
    "We follow the usual practice of setting a path (if needed), reading the data from the shape file and a quick check of its contents (`head`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2bd865b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(184, 36)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>code7</th>\n",
       "      <th>mun_name</th>\n",
       "      <th>state_init</th>\n",
       "      <th>area_km2</th>\n",
       "      <th>state_code</th>\n",
       "      <th>micro_code</th>\n",
       "      <th>micro_name</th>\n",
       "      <th>inc_mic_4q</th>\n",
       "      <th>inc_zik_3q</th>\n",
       "      <th>inc_zik_2q</th>\n",
       "      <th>...</th>\n",
       "      <th>gdp</th>\n",
       "      <th>pop</th>\n",
       "      <th>gdpcap</th>\n",
       "      <th>popdens</th>\n",
       "      <th>zik_1q</th>\n",
       "      <th>ziq_2q</th>\n",
       "      <th>ziq_3q</th>\n",
       "      <th>zika_d</th>\n",
       "      <th>mic_d</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2300101.0</td>\n",
       "      <td>Abaiara</td>\n",
       "      <td>CE</td>\n",
       "      <td>180.833</td>\n",
       "      <td>23</td>\n",
       "      <td>23019</td>\n",
       "      <td>19Âª RegiÃ£o Brejo Santo</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.00</td>\n",
       "      <td>...</td>\n",
       "      <td>35974.0</td>\n",
       "      <td>10496.0</td>\n",
       "      <td>3.427</td>\n",
       "      <td>58.043</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>POLYGON ((5433729.65 9186242.97, 5433688.546 9...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2300150.0</td>\n",
       "      <td>Acarape</td>\n",
       "      <td>CE</td>\n",
       "      <td>130.002</td>\n",
       "      <td>23</td>\n",
       "      <td>23003</td>\n",
       "      <td>3Âª RegiÃ£o MaracanaÃº</td>\n",
       "      <td>6.380399</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.00</td>\n",
       "      <td>...</td>\n",
       "      <td>68314.0</td>\n",
       "      <td>15338.0</td>\n",
       "      <td>4.454</td>\n",
       "      <td>117.983</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>POLYGON ((5476916.288 9533405.667, 5476798.561...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2300200.0</td>\n",
       "      <td>AcaraÃº</td>\n",
       "      <td>CE</td>\n",
       "      <td>842.471</td>\n",
       "      <td>23</td>\n",
       "      <td>23012</td>\n",
       "      <td>12Âª RegiÃ£o AcaraÃº</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.63</td>\n",
       "      <td>...</td>\n",
       "      <td>309490.0</td>\n",
       "      <td>57551.0</td>\n",
       "      <td>5.378</td>\n",
       "      <td>68.312</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>POLYGON ((5294389.783 9689469.144, 5294494.499...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3 rows × 36 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       code7 mun_name state_init  area_km2  state_code  micro_code  \\\n",
       "0  2300101.0  Abaiara         CE   180.833          23       23019   \n",
       "1  2300150.0  Acarape         CE   130.002          23       23003   \n",
       "2  2300200.0  AcaraÃº         CE   842.471          23       23012   \n",
       "\n",
       "                 micro_name  inc_mic_4q  inc_zik_3q  inc_zik_2q  ...  \\\n",
       "0  19Âª RegiÃ£o Brejo Santo    0.000000         0.0        0.00  ...   \n",
       "1    3Âª RegiÃ£o MaracanaÃº    6.380399         0.0        0.00  ...   \n",
       "2      12Âª RegiÃ£o AcaraÃº    0.000000         0.0        1.63  ...   \n",
       "\n",
       "        gdp      pop  gdpcap  popdens  zik_1q  ziq_2q  ziq_3q  zika_d  mic_d  \\\n",
       "0   35974.0  10496.0   3.427   58.043     0.0     0.0     0.0     0.0    0.0   \n",
       "1   68314.0  15338.0   4.454  117.983     0.0     0.0     0.0     0.0    1.0   \n",
       "2  309490.0  57551.0   5.378   68.312     0.0     1.0     0.0     1.0    0.0   \n",
       "\n",
       "                                            geometry  \n",
       "0  POLYGON ((5433729.65 9186242.97, 5433688.546 9...  \n",
       "1  POLYGON ((5476916.288 9533405.667, 5476798.561...  \n",
       "2  POLYGON ((5294389.783 9689469.144, 5294494.499...  \n",
       "\n",
       "[3 rows x 36 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Setting working folder:\n",
    "#path = \"/your/path/to/data/\"\n",
    "path = \"\"\n",
    "\n",
    "# Load the Ceará data:\n",
    "dfs = gpd.read_file(path+\"ceara/ceara.shp\")\n",
    "print(dfs.shape)\n",
    "dfs.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4024a438",
   "metadata": {},
   "source": [
    "#### Selecting variables and checking their correlation"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "c63b04c2-f618-47f4-94fa-a87e77adce0e",
   "metadata": {},
   "source": [
    "Following Chapter 12 of Anselin (2024) (https://lanselin.github.io/introbook_vol2/CHclustervalidation.html), we select the following variables from the Ceará sample data set.\n",
    "\n",
    "List of variables:\n",
    "| Column Name | Description                                      |\n",
    "|-------------|--------------------------------------------------|\n",
    "| mobility    | Mobility index                                  |\n",
    "| environ     | Environment index                               |\n",
    "| housing     | Housing index                                  |\n",
    "| sanitation  | Sanitation index                               |\n",
    "| infra       | Infrastructure index                           |\n",
    "| gdpcap      | GDP per capita                                  | \n",
    "\n",
    "We carry out the by now familiar manipulation to create the required input data for the cluster routines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8d74ec17",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Weights Meta-data:\n",
      " number of observations:                  184\n",
      "           is symmetric:                 True\n",
      "               sparsity:  0.02953686200378072\n",
      "        # min neighbors:                    1\n",
      "        # max neighbors:                   13\n",
      "       # mean neighbors:    5.434782608695652\n",
      "     # median neighbors:                  5.0\n",
      "           has isolates:                False\n",
      "\n"
     ]
    }
   ],
   "source": [
    "varlist = ['mobility', 'environ', 'housing', 'sanitation', 'infra', 'gdpcap']\n",
    "ceara_g = pygeoda.open(dfs)\n",
    "queen_w = pygeoda.queen_weights(ceara_g)\n",
    "print(queen_w)\n",
    "data = dfs[varlist]\n",
    "data_g = ceara_g[varlist]\n",
    "n_clusters = 13"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25b9fd9d",
   "metadata": {},
   "source": [
    "## Cluster Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f1a68c2",
   "metadata": {},
   "source": [
    "Before considering the validation measures, we compute the cluster solutions for hierarchical clustering (using `sklearn.AgglomerativeClustering`), K-Means (using `sklearn.KMeans`), and the spatially constrained clustering methods using `pygeoda`. For details on the arguments and helper functions, see the previous notebooks.\n",
    "\n",
    "For each cluster, we extract the `labels` and the `fit` using the respective helper functions.\n",
    "\n",
    "To illustrate the internal validity measures, we will focus on Ward's agglomerative clustering as an example of a standard method and on AZP with SCHC initial solution as an example of a spatially constrained cluster solution. For this, we also generate the cluster cardinalities. This is not pursued for the other cluster solutions, but can be readily implemented."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e3f6b9f",
   "metadata": {},
   "source": [
    "### Hierarchical clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8054882a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Labels  Cardinality\n",
      "      0           16\n",
      "      1           26\n",
      "      2           29\n",
      "      3            9\n",
      "      4           14\n",
      "      5            3\n",
      "      6           26\n",
      "      7           11\n",
      "      8           20\n",
      "      9            4\n",
      "     10           15\n",
      "     11           10\n",
      "     12            1\n"
     ]
    }
   ],
   "source": [
    "method = 'ward'\n",
    "X = StandardScaler().fit_transform(data)\n",
    "\n",
    "agg_clusters = AgglomerativeClustering(n_clusters=n_clusters, linkage=method, compute_distances=True)\n",
    "agg_clusters.fit(X)\n",
    "agg_labels = tuple(int(label) for label in agg_clusters.labels_)\n",
    "\n",
    "agg_clusters_fit = cluster_fit(data=data,clustlabels=agg_clusters.labels_,\n",
    "                 n_clusters=n_clusters, printopt=False)\n",
    "agg_stats = cluster_stats(agg_labels)\n",
    "print(agg_stats.to_string(index=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6235757f",
   "metadata": {},
   "source": [
    "### K-means clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "66e42867",
   "metadata": {},
   "outputs": [],
   "source": [
    "kmeans_clusters = KMeans(n_clusters=n_clusters, n_init=150, random_state=1234567).fit(X) \n",
    "kmeans_labels = tuple(int(label) for label in kmeans_clusters.labels_)\n",
    "kmeans_clusters_fit = cluster_fit(data=data,clustlabels=kmeans_clusters.labels_,\n",
    "                 n_clusters=n_clusters, printopt=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc24e11a",
   "metadata": {},
   "source": [
    "### SCHC with Ward’s linkage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b45e8390",
   "metadata": {},
   "outputs": [],
   "source": [
    "schc_clusters = pygeoda.schc(n_clusters, queen_w, data_g, \"ward\")\n",
    "schc_labels = schc_clusters['Clusters']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c147255",
   "metadata": {},
   "source": [
    "### SKATER"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "44102708",
   "metadata": {},
   "outputs": [],
   "source": [
    "skater_clusters = pygeoda.skater(n_clusters, queen_w, data_g)\n",
    "skater_labels = skater_clusters['Clusters']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84bdd777",
   "metadata": {},
   "source": [
    "### REDCAP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9831e9c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "redcap_clusters = pygeoda.redcap(n_clusters, queen_w, data_g, method='fullorder-wardlinkage')\n",
    "redcap_labels = redcap_clusters['Clusters']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e76a8e7",
   "metadata": {},
   "source": [
    "### AZP with simulated annealing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d7c0650a",
   "metadata": {},
   "outputs": [],
   "source": [
    "azp_sa_clusters = pygeoda.azp_sa(n_clusters, queen_w, data_g, cooling_rate=0.8, sa_maxit=5)\n",
    "azp_sa_labels = azp_sa_clusters['Clusters']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0649f6db",
   "metadata": {},
   "source": [
    "### AZP with SCHC as initial solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b98dd812",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Labels  Cardinality\n",
      "      1           89\n",
      "      2           43\n",
      "      3           15\n",
      "      4           14\n",
      "      5            6\n",
      "      6            4\n",
      "      7            4\n",
      "      8            3\n",
      "      9            2\n",
      "     10            1\n",
      "     11            1\n",
      "     12            1\n",
      "     13            1\n"
     ]
    }
   ],
   "source": [
    "azp_schc_clusters = pygeoda.azp_sa(n_clusters, queen_w, data_g, cooling_rate=0.8, sa_maxit=5,\n",
    "                                init_regions=schc_labels)\n",
    "azp_schc_labels = azp_schc_clusters['Clusters']\n",
    "azp_schc_stats = cluster_stats(azp_schc_labels)\n",
    "print(azp_schc_stats.to_string(index=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35286ad7",
   "metadata": {},
   "source": [
    "### Max-p Regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1a2c31cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "maxp_sa_clusters = pygeoda.maxp_sa(queen_w, data_g, \n",
    "                                      bound_variable=dfs['pop'], \n",
    "                                      min_bound=dfs['pop'].sum()*0.05,\n",
    "                                      iterations=9999,\n",
    "                                      cooling_rate=0.9,\n",
    "                                      sa_maxit=5)\n",
    "maxp_sa_labels = maxp_sa_clusters['Clusters']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c632b9d",
   "metadata": {},
   "source": [
    "## Internal Validation Measures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d277161f",
   "metadata": {},
   "source": [
    "As mentioned, in addition to the classic measures of fit, we also consider *fragmentation*, the *join count ratio*, and, for spatially constrained cluster solutions, the *compactness* and *diameter*.\n",
    "\n",
    "These measures are provided as attributes in the solution object created by `pygeoda.spatial_validation`. This requires the `pygeoda` data set, the cluster labels and the spatial weights as arguments.\n",
    "\n",
    "We illustrate this for Ward's agglomerative clustering, with `agg_labels` as the cluster labels, and for AZP-SCHC, with `azp_schc_labels` as the cluster labels. For both, the data set is `ceara_g` and the spatial weights are contained in `queen_w`.\n",
    "\n",
    "We store the results in, respectively, `agg_validation` and `azp_schc_validation`. These objects will then be used as arguments to the helper functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ba5c9ff1",
   "metadata": {},
   "outputs": [],
   "source": [
    "agg_validation = pygeoda.spatial_validation(ceara_g, agg_labels, queen_w)\n",
    "azp_schc_validation = pygeoda.spatial_validation(ceara_g, azp_schc_labels, queen_w)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cbe6ca9",
   "metadata": {},
   "source": [
    "### Fragmentation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bbeab36",
   "metadata": {},
   "source": [
    "The fragmentation measures are computed from the makeup of the cluster components. An ideally balanced cluster is when each component has the same number of observations. This is quantified by means of entropy and its standardized counterpart, as well as by Simpson's index and its standardized counterpart. For entropy, larger values suggest a greater balance, whereas for Simpson's index, it is the other way around.\n",
    "\n",
    "The `pygeoda.spatial_validation` return object includes the fragmentation information in two attributes: `fragmentation` and `cluster_fragmentation`. The first contains the overall measures, as well as the number of clusters, in `fragmentation.n`, `fragmentation.entropy`, `fragmentation.std_entropy`, `fragmentation.simpson` and `fragmentation.std_simpson`. The second has the same information, organized as a list by cluster. It shows the within-cluster fragmentation for clusters that are not spatially constrained.\n",
    "\n",
    "The `cluster_fragmentation` helper function, takes a data frame with the labels and cardinalities (created by `cluster_stats`), and the `cluster_fragmentation`, `fragmentation` and `spatially_constrained` attributes from the validation object. The `spatially_constrained` flag is used to limit the fragmentation output to the totals only for spatially constrained clusters.\n",
    "\n",
    "This is illustrated for Ward's agglomerative clustering and AZP-SCHC. Note that for the latter, only the totals are given, since it is a spatially constrained solution.\n",
    "\n",
    "The other cluster solutions can be analyzed in the same way."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44825e1c",
   "metadata": {},
   "source": [
    "#### Agglomerative clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "690e53ca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fragmentation\n",
      "Label   N Sub  Entropy  Entropy*  Simpson  Simpson*\n",
      "    0  16   9 1.751176  0.796994 0.250892  2.258026\n",
      "    1  26  14 2.397937  0.908634 0.115385  1.615385\n",
      "    2  29  14 2.425806  0.919194 0.109467  1.532544\n",
      "    3   9  11 2.250260  0.938431 0.120000  1.320000\n",
      "    4  14   9 1.923066  0.875225 0.187500  1.687500\n",
      "    5   3   8 1.933810  0.929966 0.164444  1.315556\n",
      "    6  26  10 2.168223  0.941647 0.132653  1.326531\n",
      "    7  11   9 2.145842  0.976615 0.123967  1.115702\n",
      "    8  20   6 1.609438  0.898244 0.240000  1.440000\n",
      "    9   4   7 1.831020  0.940958 0.185185  1.296296\n",
      "   10  15   0 0.000000  0.000000 0.000000  0.000000\n",
      "   11  10   3 1.098612  1.000000 0.333333  1.000000\n",
      "   12   1   0 0.000000  0.000000 0.000000  0.000000\n",
      "  All 184     2.351159  0.916649 0.106274  1.381557\n"
     ]
    }
   ],
   "source": [
    "agg_frag = cluster_fragmentation(agg_stats,agg_validation.cluster_fragmentation,\n",
    "                   agg_validation.fragmentation,agg_validation.spatially_constrained)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68fbbb8b",
   "metadata": {},
   "source": [
    "#### AZP-SCHC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "30ac358b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fragmentation\n",
      "Label   N Sub  Entropy  Entropy*  Simpson  Simpson*\n",
      "  All 184     1.599116  0.623449 0.303521   3.94577\n"
     ]
    }
   ],
   "source": [
    "azp_schc_frag = cluster_fragmentation(azp_schc_stats,azp_schc_validation.cluster_fragmentation,\n",
    "                   azp_schc_validation.fragmentation,azp_schc_validation.spatially_constrained)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d81ad2b6",
   "metadata": {},
   "source": [
    "### Join Count Ratio"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7672c2b2",
   "metadata": {},
   "source": [
    "The join count ratio is a *spatial* measure of the degree of internal connectedness in a cluster solution. It is computed for each cluster separately as well as for the cluster solution as a whole. It is a count of how many neighbors of observations in a cluster are also members of that cluster.\n",
    "\n",
    "The relevant measures are included in the `joincount_ratio` and `all_joincount_ratio` attributes of the `pygeoda.spatial_validation` solution object. The former is a list with k entries, which are themselves objects, containing attributes `n`, `neighbors`, `join_count` and `ratio`. The latter is the same for the overall cluster solution.\n",
    "\n",
    "The result is provided by the `cluster_joincount` helper function. It takes the data frame with cluster cardinalities and the `joincount_ratio` and `all_joincount_ratio` attributes from the validation solution.\n",
    "\n",
    "We use the same examples in the illustration below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94c909fa",
   "metadata": {},
   "source": [
    "#### Agglomerative clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "bcfc5dee",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Join Count Ratio\n",
      "Label   N  Neighbors  Join Count  Ratio\n",
      "    0  16         96          20  0.208\n",
      "    1  26        134          26  0.194\n",
      "    2  29        158          64  0.405\n",
      "    3   9         47           6  0.128\n",
      "    4  14         80           8  0.100\n",
      "    5   3         15           0  0.000\n",
      "    6  26        158          30  0.190\n",
      "    7  11         43           4  0.093\n",
      "    8  20        116          18  0.155\n",
      "    9   4         21          10  0.476\n",
      "   10  15         83          18  0.217\n",
      "   11  10         46           8  0.174\n",
      "   12   1          3           0  0.000\n",
      "  All 184       1000         212  0.212\n"
     ]
    }
   ],
   "source": [
    "agg_jc = cluster_joincount(agg_stats,agg_validation.joincount_ratio,\n",
    "                   agg_validation.all_joincount_ratio)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2e36fee",
   "metadata": {},
   "source": [
    "#### AZP-SCHC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "f4656f10",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Join Count Ratio\n",
      "Label   N  Neighbors  Join Count  Ratio\n",
      "    0  89        474         342  0.722\n",
      "    1  43        225         128  0.569\n",
      "    2  15         96          42  0.438\n",
      "    3  14         78          44  0.564\n",
      "    4   6         37          10  0.270\n",
      "    5   4         17           6  0.353\n",
      "    6   4         21          10  0.476\n",
      "    7   3         22           4  0.182\n",
      "    8   2         12           2  0.167\n",
      "    9   1          3           0  0.000\n",
      "   10   1          3           0  0.000\n",
      "   11   1          7           0  0.000\n",
      "   12   1          5           0  0.000\n",
      "  All 184       1000         588  0.588\n"
     ]
    }
   ],
   "source": [
    "azp_schc_jc = cluster_joincount(azp_schc_stats,azp_schc_validation.joincount_ratio,\n",
    "                   azp_schc_validation.all_joincount_ratio)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3bf6aba",
   "metadata": {},
   "source": [
    "### Compactness"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "586d257b",
   "metadata": {},
   "source": [
    "Compactness is a criterion that is only applicable to spatially constrained cluster solutions. It measures the ratio of the perimeter of the cluster to that of a circle with the same area. The `compactness` attribute of the `spatial_validation` object is a list with k items, each an object with attributes `area`, `perimeter` and `isoperimeter_quotient`. The closer the IPQ is to one, the more compact is the cluster shape.\n",
    "\n",
    "The helper function `cluster_compactness` extracts this information. Its attributes are the data frame with cluster cardinalities (from `cluster_stats`), the `compactness` attribute, and the `spatially_constrained` attribute. Compactness is not relevant for clusters that are not spatially constrained and therefore the helper function will yield an error message when this happens.\n",
    "\n",
    "We continue with the same two examples. Note that for Ward's agglomerative cluster, an error message is generated."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea631e33",
   "metadata": {},
   "source": [
    "#### Agglomerative clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "d5c2a0b3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error: Compactness is only applicable to spatially constrained clusters\n"
     ]
    }
   ],
   "source": [
    "agg_compactness = cluster_compactness(agg_stats,agg_validation.compactness,\n",
    "                                      agg_validation.spatially_constrained)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bd748b4",
   "metadata": {},
   "source": [
    "#### AZP-SCHC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "dd82ab5b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compactness\n",
      " Label  N         Area    Perimeter      IPQ\n",
      "     0 89 8.672968e+10 1.536222e+07 0.004618\n",
      "     1 43 2.797625e+10 6.198588e+06 0.009150\n",
      "     2 15 1.053291e+10 2.220946e+06 0.026834\n",
      "     3 14 1.138704e+10 2.122643e+06 0.031759\n",
      "     4  6 4.038275e+09 9.257796e+05 0.059209\n",
      "     5  4 3.085985e+09 5.499120e+05 0.128238\n",
      "     6  4 1.779537e+09 4.240068e+05 0.124386\n",
      "     7  3 9.969357e+08 3.168568e+05 0.124782\n",
      "     8  2 1.474872e+09 3.229358e+05 0.177718\n",
      "     9  1 6.171934e+08 1.199243e+05 0.539283\n",
      "    10  1 7.908164e+07 5.056877e+04 0.388616\n",
      "    11  1 8.449201e+08 1.739887e+05 0.350738\n",
      "    12  1 4.124364e+08 1.255397e+05 0.328855\n"
     ]
    }
   ],
   "source": [
    "azp_schc_compactness = cluster_compactness(azp_schc_stats,azp_schc_validation.compactness,\n",
    "                                      azp_schc_validation.spatially_constrained)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd7936d6",
   "metadata": {},
   "source": [
    "### Diameter"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b23b927",
   "metadata": {},
   "source": [
    "The diameter of a spatially constrained cluster is an alternative measure of compactness, based on the network structure reflected in the spatial weights. The diameter of a cluster is the number of steps in the spatial weights graph that corresponds with the longest shortest path between any pair of observations (Newman 2018). Since this number will increase with cluster size, it is also standardized by dividing by the number of cluster members. Note that when a cluster is a singleton, the diameter will be zero.\n",
    "\n",
    "The `diameter` attribute of the `pygeoda.spatial_validation` object is a list with k items, one for each cluster, as an object with attributes `steps` and `ratio`. The helper function `cluster_diameter` extracts this information as a data frame. It takes as arguments the cluster cardinalities (from `cluster_stats`), and the `diameter` and `spatially_constrained` attributes from the `pygeoda.spatial_validation` object.\n",
    "\n",
    "As in the case of compactness, an error message is generated for clusters that are not spatially constrained.\n",
    "\n",
    "Again, we use the same two examples."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "635999dd",
   "metadata": {},
   "source": [
    "#### Agglomerative clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "4885a986",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error: Diameter is only applicable to spatially constrained clusters\n"
     ]
    }
   ],
   "source": [
    "agg_diam = cluster_diameter(agg_stats,agg_validation.diameter,\n",
    "                            agg_validation.spatially_constrained)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9af797b4",
   "metadata": {},
   "source": [
    "#### AZP-SCHC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "27b3b1cd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diameter\n",
      " Label  N  Steps    Ratio\n",
      "     0 89     22 0.247191\n",
      "     1 43     17 0.395349\n",
      "     2 15      6 0.400000\n",
      "     3 14      7 0.500000\n",
      "     4  6      3 0.500000\n",
      "     5  4      3 0.750000\n",
      "     6  4      2 0.500000\n",
      "     7  3      2 0.666667\n",
      "     8  2      1 0.500000\n",
      "     9  1      0 0.000000\n",
      "    10  1      0 0.000000\n",
      "    11  1      0 0.000000\n",
      "    12  1      0 0.000000\n"
     ]
    }
   ],
   "source": [
    "azp_schc_diam = cluster_diameter(azp_schc_stats,azp_schc_validation.diameter,\n",
    "                            azp_schc_validation.spatially_constrained)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5a1ee80",
   "metadata": {},
   "source": [
    "### Overall Comparison of Internal Validation Measures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae11b981",
   "metadata": {},
   "source": [
    "We conclude our discussion of internal validation measures with an overview of the main non-spatial measures for all the cluster solutions considered above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "6586219c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      Method  Spat. Const.      WSS  BSS/TSS  Join Count  Entropy  Simpson\n",
      "Hierarchical         False 349.0000   0.6800       0.212   2.3512   0.1063\n",
      "     K-Means         False 335.7500   0.7000       0.210   2.3712   0.1015\n",
      "        SCHC          True 568.0173   0.4827       0.668   1.6394   0.2730\n",
      "      SKATER          True 604.2209   0.4497       0.784   1.3661   0.3901\n",
      "      REDCAP          True 562.7003   0.4875       0.660   1.6109   0.2769\n",
      "         AZP          True 617.0010   0.4381       0.526   1.9521   0.2016\n",
      " AZP_Initial          True 538.3477   0.5097       0.588   1.5991   0.3035\n",
      "       Max-p          True 745.2474   0.3213       0.496   2.4175   0.0959\n"
     ]
    }
   ],
   "source": [
    "clusters = [\n",
    "    agg_clusters_fit, kmeans_clusters_fit, schc_clusters, skater_clusters,\n",
    "    redcap_clusters, azp_sa_clusters, azp_schc_clusters, maxp_sa_clusters\n",
    "]\n",
    "labels = [\n",
    "    agg_labels, kmeans_labels, schc_labels, skater_labels,\n",
    "    redcap_labels, azp_sa_labels, azp_schc_labels, maxp_sa_labels\n",
    "]\n",
    "label_names = [\n",
    "    'Hierarchical', 'K-Means', 'SCHC', 'SKATER',\n",
    "    'REDCAP', 'AZP', 'AZP_Initial', 'Max-p'\n",
    "]\n",
    "\n",
    "results = []\n",
    "\n",
    "# Run pygeoda.spatial_validation for each label set\n",
    "for cluster, label, name in zip(clusters, labels, label_names):\n",
    "    result = pygeoda.spatial_validation(ceara_g, label, queen_w)\n",
    "    \n",
    "    try:\n",
    "        wss = np.round(cluster['Total within-cluster sum of squares'],4)\n",
    "        bss_tss = np.round(cluster['The ratio of between to total sum of squares'],4)\n",
    "    except:\n",
    "        try:\n",
    "            wss = np.round(cluster[\"WSS\"],2)\n",
    "            bss_tss = np.round(cluster[\"Ratio\"],2)\n",
    "        except:\n",
    "            wss = None\n",
    "            bss_tss = None\n",
    "\n",
    "    spatially_constrained = result.spatially_constrained\n",
    "    all_join_count_ratio = np.round(result.all_joincount_ratio.ratio,4)\n",
    "    entropy = np.round(result.fragmentation.entropy,4)\n",
    "    simpson = np.round(result.fragmentation.simpson,4)\n",
    "    \n",
    "    results.append({\n",
    "        'Method': name,\n",
    "        'Spat. Const.': spatially_constrained,\n",
    "        'WSS': wss,\n",
    "        'BSS/TSS': bss_tss,\n",
    "        'Join Count': all_join_count_ratio,\n",
    "        'Entropy': entropy,\n",
    "        'Simpson': simpson\n",
    "    })\n",
    "\n",
    "validation = pd.DataFrame(results)\n",
    "print(validation.to_string(index=False))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b93f7859",
   "metadata": {},
   "source": [
    "## External Validation Measures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5924af07",
   "metadata": {},
   "source": [
    "External validation measures are designed to compare a cluster solution to a *truth*, but they can also be employed to compare several cluster solutions to each other. The validation indices reveal how close the cluster solutions are. We consider two measures, the *adjusted rand index* and the *normalized information distance*."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "959e9826",
   "metadata": {},
   "source": [
    "### Adjusted Rand Index (ARI)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b44796ec",
   "metadata": {},
   "source": [
    "The adjusted rand index is based on counting how many pairs of observations are in the same grouping for two cluster solutions. It can be computed by `sklearn.metrics.adjusted_rand_score`. The two arguments are numpy arrays of the labels of the *reference* (first argument) and the labels of the cluster to be compared.\n",
    "\n",
    "Note that we need to convert our `labels` solution to a numpy array to make this work.\n",
    "\n",
    "For example, the ARI between Ward's agglomerative solution and K-Means is found by passing numpy arrays for `agg_labels` and `kmeans_labels`. The result of 0.414 suggests only low correspondence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "107f3fa3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.414\n"
     ]
    }
   ],
   "source": [
    "ari = adjusted_rand_score(np.array(agg_labels),np.array(kmeans_labels))\n",
    "print(np.round(ari,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3637c7b3",
   "metadata": {},
   "source": [
    "We can now compute all pairwise indices with a simple loop. We first recreate the `labels` and `label_names` lists from above (so this can be run without the internal validation measures)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "84e3f9a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = [\n",
    "    agg_labels, kmeans_labels, schc_labels, skater_labels,\n",
    "    redcap_labels, azp_sa_labels, azp_schc_labels, maxp_sa_labels\n",
    "]\n",
    "label_names = [\n",
    "    'Hierarchical', 'K-Means', 'SCHC', 'SKATER',\n",
    "    'REDCAP', 'AZP', 'AZP_Initial', 'Max-p'\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b270393",
   "metadata": {},
   "source": [
    "In a simple loop, we compute all pairwise indices and populate a matrix. This is then turned into a data frame and printed. Note that the matrix is symmetric and the diagonal values of 1.0 can be ignored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "791ce3bf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              Hierarchical  K-Means   SCHC  SKATER  REDCAP    AZP  \\\n",
      "Hierarchical         1.000    0.414  0.131   0.079   0.145  0.160   \n",
      "K-Means              0.414    1.000  0.151   0.073   0.155  0.166   \n",
      "SCHC                 0.131    0.151  1.000   0.424   0.918  0.504   \n",
      "SKATER               0.079    0.073  0.424   1.000   0.384  0.258   \n",
      "REDCAP               0.145    0.155  0.918   0.384   1.000  0.516   \n",
      "AZP                  0.160    0.166  0.504   0.258   0.516  1.000   \n",
      "AZP_Initial          0.155    0.147  0.735   0.452   0.727  0.577   \n",
      "Max-p                0.089    0.100  0.194   0.173   0.185  0.218   \n",
      "\n",
      "              AZP_Initial  Max-p  \n",
      "Hierarchical        0.155  0.089  \n",
      "K-Means             0.147  0.100  \n",
      "SCHC                0.735  0.194  \n",
      "SKATER              0.452  0.173  \n",
      "REDCAP              0.727  0.185  \n",
      "AZP                 0.577  0.218  \n",
      "AZP_Initial         1.000  0.163  \n",
      "Max-p               0.163  1.000  \n"
     ]
    }
   ],
   "source": [
    "h = len(labels)\n",
    "allari = np.zeros((h,h))\n",
    "for i in range(h):\n",
    "    labi = np.array(labels[i])\n",
    "    for j in range(h):\n",
    "        labj = np.array(labels[j])\n",
    "        allari[i,j] = adjusted_rand_score(labi,labj)\n",
    "dfari = pd.DataFrame(allari, columns = label_names, index = label_names)\n",
    "print(np.round(dfari,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccbadb09",
   "metadata": {},
   "source": [
    "As in Chapter 12, the matrix reveals a much closer correspondence among the non-spatial solutions and the spatial solutions respectively. The greatest correspondence is between SCHC and Redcap, with an ARI of 0.918."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "571fc29a",
   "metadata": {},
   "source": [
    "### Normalized Information Distance (NID)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6526c16c",
   "metadata": {},
   "source": [
    "The second external validation measure is based on information-theoretic considerations, such as entropy. In Chapter 12, the normalized information distance is introduced (NID). A close counterpart can be computed by means of `sklearn.metrics.adjusted_mutual_info_score`. The arguments are the same as for ARI. However, in contrast to NID as presented in Chapter 12, a higher value for the adjusted mutual information score indicates closer similarity.\n",
    "\n",
    "We first illustrate this for Ward's agglomerative clustering and K-Means, passing numpy arrays of `agg_labels` and `kmeans_labels`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "9a61c8f1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.601\n"
     ]
    }
   ],
   "source": [
    "nid = adjusted_mutual_info_score(np.array(agg_labels),np.array(kmeans_labels))\n",
    "print(np.round(nid,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd63138f",
   "metadata": {},
   "source": [
    "Finally, we run the same loop as for ARI to compute all pairwise NID scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "e473b192",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              Hierarchical  K-Means   SCHC  SKATER  REDCAP    AZP  \\\n",
      "Hierarchical         1.000    0.601  0.257   0.219   0.267  0.306   \n",
      "K-Means              0.601    1.000  0.299   0.220   0.297  0.299   \n",
      "SCHC                 0.257    0.299  1.000   0.548   0.897  0.581   \n",
      "SKATER               0.219    0.220  0.548   1.000   0.496  0.425   \n",
      "REDCAP               0.267    0.297  0.897   0.496   1.000  0.595   \n",
      "AZP                  0.306    0.299  0.581   0.425   0.595  1.000   \n",
      "AZP_Initial          0.317    0.320  0.727   0.545   0.727  0.646   \n",
      "Max-p                0.177    0.207  0.388   0.423   0.372  0.445   \n",
      "\n",
      "              AZP_Initial  Max-p  \n",
      "Hierarchical        0.317  0.177  \n",
      "K-Means             0.320  0.207  \n",
      "SCHC                0.727  0.388  \n",
      "SKATER              0.545  0.423  \n",
      "REDCAP              0.727  0.372  \n",
      "AZP                 0.646  0.445  \n",
      "AZP_Initial         1.000  0.367  \n",
      "Max-p               0.367  1.000  \n"
     ]
    }
   ],
   "source": [
    "allnid = np.zeros((h,h))\n",
    "for i in range(h):\n",
    "    labi = np.array(labels[i])\n",
    "    for j in range(h):\n",
    "        labj = np.array(labels[j])\n",
    "        allnid[i,j] = adjusted_mutual_info_score(labi,labj)\n",
    "dfnid = pd.DataFrame(allnid, columns = label_names, index = label_names)\n",
    "print(np.round(dfnid,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f117d025",
   "metadata": {},
   "source": [
    "As for ARI, we find the closest correspondence between SCHC and Redcap."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "myenv12",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}