Exploring_Linguistic_Diversity_Using_R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Exploring Linguistic Diversity Using R\n",
    "\n",
    "## Amir Barghi\n",
    "\n",
    "### Mathematics and Statistics \n",
    "\n",
    "-----\n",
    "\n",
    "## Social Science Seminar\n",
    "\n",
    "### Thursday, October 20, 2022"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Greenberg (1956) defined linguistic diversity as the probability that two randomly chosen individuals (with replacement) from a population have different first languages. In other words, if there are $n$ languages in a population and $p_i$ denotes the probability that an individual chosen at random from this population has the $i$th language as their first language, for $i \\in \\{1,\\ldots,n\\}$, then Greenberg's index of linguistic diversity (LDI) is calculated as\n",
    "$$ H_2 = 1 - \\sum_{i=1} ^n p_i^2.$$\n",
    "\n",
    "----\n",
    "- M. J. H. Greenberg, The measurement of linguistic diversity, *Language*, vol. 32, pp. 109--115, Jan.-Mar. 1956."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In ecology, $H_2$ is known as Gini-Simpson index (named after Italian Statistician, <a href=\"https://www.britannica.com/biography/Corrado-Gini\">Corrado Gini</a>, and after Edward Hugh Simpson (1949)), and as discussed by Jost (2006), it is measure of entropy (or uncertainty) and not diversity. According to Jost (2006), a diversity index should measure \"the effective number of elements in a system.\" \n",
    "In ecology, this is **the number of equiprobable species in a population having the same entropic index**. \n",
    "\n",
    "-----\n",
    "- E. Simpson, Measurement of diversity, *Nature*, vol 163, p. 688, April 1949.\n",
    "\n",
    "- L. Jost, Entropy and diversity, *Oikos*, vol. 113, pp. 363--375, Jan. 2006."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Jost (2006) shows that, in order to use the Gini-Simpson index (or equivalently, LDI) as a measure of diversity, one needs to consider the following transformation of $H_2$:\n",
    "$$ D_2 = \\frac{1}{1-H_2}  = \\frac{1}{\\sum_{i=1} ^n p_i ^2}.$$\n",
    "This is also known as the *effective number* of a population which was introduced by MacArthur (1965). \n",
    "\n",
    "-----\n",
    "\n",
    "\n",
    "- L. Jost, Entropy and diversity, *Oikos*, vol. 113, pp. 363--375, Jan. 2006.\n",
    "\n",
    "- R. H. MacArthur, Patterns of species diversity, *Biol. Rev. Camb. Philos. Soc.*, vol. 40, pp. 510--533, Nov. 1965."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Preliminaries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Here we discuss different entropic indices and how to transform them into diversity measures, i.e., effective numbers. Throughout this exploration, we will use diversity and effective number interchangeably. We also assume that there are $n$ languages in a population or a country, and $p_i$ denotes the probability that an individual chosen at random from this population has the $i$th language as their first language, for $i \\in \\{1,\\ldots,n\\}$. We denote an entropy and a diversity by $H_q$ and $D_q$, respectively, where $q$ denotes the order of the entropy and of the diversity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Richness\n",
    "\n",
    "Linguistic richness is the simplest measure of diversity and it counts the number of languages in a country. In mathematical terms, the entropic index is defined as \n",
    "$$ H_0 = \\sum_{i = 1} ^n p_i ^0,$$\n",
    "and it is not difficult to see that its effective number is \n",
    "$$D_0 = \\sum_{i = 1} ^n p_i ^0.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Shannon\n",
    "\n",
    "Shannon entropy, the very first measure of entropy in the context of information theory, was introduced by Shannon (1948). It is defined as \n",
    "$$H_1 = \\sum_{i=1} ^n p_i \\ln(p_i),$$\n",
    "and it measures the expected value of uncertainty in a system. Hill (1973) shows that the effective number of this entropy is calculated as \n",
    "$$ D_1 = \\exp(H_1) = \\exp\\left( - \\sum_{i=1} ^n p_i \\ln(p_i) \\right).$$\n",
    "\n",
    "----\n",
    "- C. E. Shannon, A mathematical theory of communication, *Bell System Tech. J.*, vol. 27, pp. 379--423, 623--656, July, October, 1948.\n",
    "\n",
    "- M. Hill, Diversity and evenness: a unifying notation and its consequences, *Ecology*, vol. 54, p. 427---432, Mar. 1973. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Greenberg/Gini-Simpson\n",
    "\n",
    "We saw earlier that Greenberg's index of linguistic diversity and Gini-Simpson index in ecology are the same, and we discussed how to transform them into effective numbers:\n",
    "$$ H_2 = 1 - \\sum_{i=1} ^n p_i^2$$\n",
    "and \n",
    "$$ D_2 = \\frac{1}{1-H_2}  = \\frac{1}{\\sum_{i=1} ^n p_i ^2}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## HCDT\n",
    "\n",
    "HCDT entropy was introduced by Tsallis (1988) in statistical mechanics and is defined as \n",
    "$$ H_q = \\frac{\\left(1 - \\sum_{i = 1} ^n p_i ^q \\right)}{q-1} ,$$\n",
    "for $q \\neq 1$.\n",
    "Jost (2006) implicitly shows that the effective number associated with this entropy is\n",
    "$$ D_q =  (1-(q-1)H_q)^{\\frac{1}{1-q}} = \\left(\\sum_{i = 1} ^n p_i ^q\\right)^{\\frac{1}{1-q}} ,$$\n",
    "which are known as a Hill number, introduced in ecology by Hill (1973). Note that richness and Greenberg diversity indices are special case of $D_q$ for $q = 0$ and $q=2$, respectively.\n",
    "\n",
    "----\n",
    "- C. Tsallis, Possible generalization of Boltzmann-Gibbs statistics, *J. Stat. Phys.*, vol. 52, pp. 479--487, Jul. 1988.\n",
    "\n",
    "- L. Jost, Entropy and diversity, *Oikos*, vol. 113, pp. 363--375, Jan. 2006.\n",
    "\n",
    "- M. Hill, Diversity and evenness: a unifying notation and its consequences, *Ecology*, vol. 54, p. 427---432, Mar. 1973. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Rényi Entropy\n",
    "\n",
    "Introduced by Rényi (1961) as a generalization of Shannon's entropy, Rényi entropy is defined as\n",
    "$$ H_q = \\frac{ - \\ln\\left( \\sum_{i = 1} ^n p_i ^q \\right)}{q-1},$$\n",
    "for $q \\neq 1$. In particular, Shannon entropy is the limit of $H_q$ as $q \\to 1$. Jost (2006) implicitly shows that the effective number for Renyi entropy is \n",
    "$$ D_q = \\exp(H) =  \\left(\\sum_{i = 1} ^n p_i ^q\\right)^{\\frac{1}{1-q}} . $$ \n",
    "\n",
    "----\n",
    "-  A. Rényi, On Measures of Entropy and Information, in *Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Contributions to the Theory of Statistics*, pp. 547--561, University of California Press, Berkeley, Calif., 1961.\n",
    "\n",
    "- L. Jost, Entropy and diversity, *Oikos*, vol. 113, pp. 363--375, Jan. 2006."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Hill Numbers as Measures of Linguistic Diversity\n",
    "\n",
    "Notice that the effective numbers computed from HCDT entropy and Renyi entropy are the same, and in this exploration, we would refer to them as the Hill numbers. We will consider two cases: $0 < q < 1$ and $1 < q$. As it is pointed out by Keylock (2005), when $0< q <1$, Hill numbers favor languages with smaller $p_i$'s (more rare languages), while for $1 < q$, they favor languages with larger $p_i$'s (more common languages). Also, as shown by Hill (1973), the effective number of Shannon entropy, i.e., $\\exp(H_1)$, is the limit of $D_q$ as $q \\to 1$. Moreover, recall that richness and Greenberg diversity indices are special case of $D_q$ for $q = 0$ and $q=2$, respectively. \n",
    "\n",
    "----\n",
    "- C. J. Keylock, Simpson diversity and the Shannon–Wiener index as special cases of a generalized entropy, *Oikos*, vol. 109, pp. 203--207, Apr. 2005. \n",
    "\n",
    "- M. Hill, Diversity and evenness: a unifying notation and its consequences, *Ecology*, vol. 54, p. 427---432, Mar. 1973. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Alpha, Beta, and Gamma\n",
    "\n",
    "Similar to linguistic effective numbers, alpha, beta, gamma diversities are computed based on alpha, beta, and gamma entropic indices, denoted by $H_{\\alpha}$, $H_{\\beta}$, and $H_{\\gamma}$, respectively. These diversities were briefly introduced and developed by Whitakker (1960, 1972). Given a group of countries or communities, an alpha entropy $H_{\\alpha}$ is the average of entropies within this group. This is also known as the conditional entropy as discussed by Jost (2006). A gamma entropy $H_{\\gamma}$ is the pooled entropy within this group, i.e., all the countries or communities combined, and a beta entropy $H_{\\beta}$ is the relative entropy between the countries or communities; see Jost (2006, 2008). All three entropies are then transformed into effective numbers using the appropriate transformations, resulting in $D_{\\alpha}$, $D_{\\gamma}$, and $D_{\\beta}$, respectively. Jost (2007) argues that\n",
    "$D_{\\alpha} D_{\\beta} = D_{\\gamma}$.\n",
    "\n",
    "----\n",
    "- R. H. Whittaker, Vegetation of the Siskiyou Mountains, Oregon and California, *Ecol. Monogr.*, vol. 30, pp. 279--338, Jul. 1960.\n",
    "\n",
    "- R. H. Whittaker, Evolution and Measurement of Species Diversity, *Taxon*, vol. 21, pp. 213--251, May 1972.\n",
    "\n",
    "- L. Jost, Entropy and diversity, *Oikos*, vol. 113, pp. 363--375, Jan. 2006.\n",
    "\n",
    "- L. Jost, Partitioning diversity into independent alpha and beta components, *Ecology*, vol. 88, pp. 2427–2439, Oct. 2007."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Weighted Alpha\n",
    "\n",
    "Let $N$ be number of samples or populations (in our exploration, countries) and let $S$ be the total number of species (in our case, languages). Let us denote the proportion of L1 speakers of the $j$th language in the $i$th country, where $1 \\leq j \\leq S$ and $1 \\leq i \\leq N$, by $p_{j,i}$. Suppose $w_i$ is the weight of the $i$th country. In case of unweighted alpha, beta, and gamma diversities, we assume that $w_i = 1/N$. However, for weighted diversities, we let $w_i$ to be the ratio between the $i$th country's population and the total population of the countries in this exploration. Jost (2007) shows that, for $q \\neq 1$,  \n",
    "$$ D_{w, \\alpha} ^q = \\left ( \\frac{\\sum_{i = 1} ^N w_i ^q \\sum_{j = 1}^S p_{j,i}^q}{\\sum_{i = 1} ^N w_i ^q} \\right)^{1/(1-q)}$$\n",
    "and, for $q = 1$,\n",
    "$$ D_{w, \\alpha} ^1 = \\exp \\left ( - \\sum_{i = 1} ^N w_i \\sum_{j = 1}^S p_{j,i} \\ln(p_{j,i}) \\right).$$\n",
    "\n",
    "----\n",
    "- L. Jost, Partitioning diversity into independent alpha and beta components, *Ecology*, vol. 88, pp. 2427–2439, Oct. 2007."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Unweighted Alpha\n",
    "\n",
    "When $w_i = 1/N$ for all $i$, we have\n",
    "$$ D_{u, \\alpha} ^q = \\left ( \\frac{\\sum_{i = 1}^N \\sum_{j = 1}^S p_{j,i}^q}{N} \\right)^{1/(1-q)}$$\n",
    "and, for $q = 1$,\n",
    "$$ D_{u, \\alpha} ^1 = \\exp \\left (- \\frac{\\sum_{i = 1}^N \\sum_{j = 1}^S p_{j,i} \\ln(p_{j,i})}{N} \\right).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Weighted Gamma\n",
    "\n",
    "For $q \\neq 1$,  \n",
    "$$ D_{w, \\gamma} ^q = \\left ( \\frac{\\sum_{j = 1}^S \\left( \\sum_{i = 1} ^N w_i p_{j,i} \\right)^q}{\\sum_{i = 1} ^N w_i ^q} \\right)^{1/(1-q)}$$\n",
    "and, for $q = 1$,\n",
    "$$ D_{w, \\gamma} ^1 = \\exp \\left ( - \\sum_{j = 1}^S \\left( \\sum_{i = 1} ^N w_i p_{j,i} \\right) \\ln\\left(\\sum_{i = 1} ^N w_i p_{j,i}\\right) \\right).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Unweighted Gamma\n",
    "\n",
    "When $w_i = 1/N$ for all $i$, we have\n",
    "$$ D_{u, \\gamma} ^q = \\left ( \\frac{ \\sum_{j = 1}^S \\left (\\sum_{i = 1}^N p_{j,i}\\right) ^q}{N} \\right)^{1/(1-q)}$$\n",
    "and, for $q = 1$,\n",
    "$$ D_{u, \\gamma} ^1 = \\exp \\left (- \\sum_{j = 1}^S \\left( \\frac{\\sum_{i = 1} ^N p_{j,i}}{N} \\right) \\ln\\left(\\frac{\\sum_{i = 1} ^N p_{j,i}}{N} \\right) \\right).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Weighted and Unweighted Beta\n",
    "\n",
    "For $q \\geq 0$, we have\n",
    "$$ D_{w, \\beta} ^q = \\frac{D_{w, \\gamma}^q}{D_{w, \\alpha}^q}$$\n",
    "and\n",
    "$$ D_{u, \\beta} ^q = \\frac{D_{u, \\gamma}^q}{D_{u, \\alpha}^q}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## MacArthur’s Homogeneity Measure and Relative Homogeneity\n",
    "\n",
    "Related to alpha, beta, gamma diversities are the notions of similarity and relative homogeneity. Let us assume that we are studying $N$ countries or communities. Since $D_{\\beta}$ is relative diversity between them, MacArthur's measure of homogeneity can be calculated as $1/D_{\\beta}$. According to Jost (2007), this \"answers the question 'What proportion of total diversity is found within the average community or sample?'\" As Jost points out, the range for $1/D_{\\beta}$ is the interval $[1/N,1]$, where $1/N$ represents heterogeneity and 1 homogeneity among the countries or communities. Consequently, Jost suggests \n",
    "$$ R = \\frac{1/D_{\\beta} - 1/N}{1 - 1/N}$$\n",
    "as a measure of relative homogeneity. \n",
    "\n",
    "----\n",
    "- R. H. MacArthur, Patterns of species diversity, *Biol. Rev. Camb. Philos. Soc.*, vol. 40, pp. 510--533, Nov. 1965.\n",
    "\n",
    "- L. Jost, Partitioning diversity into independent alpha and beta components, *Ecology*, vol. 88, pp. 2427--2439, Oct. 2007. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Sørensen-Dice and Jaccard Indices\n",
    "\n",
    "Both of these indices gauge similarity between two samples or populations. Sørensen-Dice index, introduced independently by Sørensen (1948) and Dice (1945), is computed as follows. Suppose $X$ and $Y$ are two samples or populations. Then, the Sørensen-Dice index for $X$ and $Y$ is\n",
    "$$ \\dfrac{2|X \\cap Y|}{|X| + |Y|} = \\dfrac{2|X \\cap Y|}{|X \\cup Y| + |X \\cap Y| }.$$\n",
    "\n",
    "\n",
    "-----\n",
    "- L. R. Dice, Measures of the amount of ecologic association between species, *Ecology*, vol. 26, pp. 297--302, Jul. 1945.\n",
    "\n",
    "- T. Sørensen, A method of establishing groups of equal amplitude in plant sociology based on similarity of species and its application to analyses of the vegetation on Danish commons, *Kongelige Danske Videnskabernes Selskab*, vol. 5, pp. 1--34, 1948"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Sørensen-Dice and Jaccard Indices\n",
    "\n",
    "Jaccard index, introduced in 1912, computes similarity between $X$ and $Y$ as\n",
    "$$ \\dfrac{|X \\cap Y|}{|X \\cup Y|} = \\dfrac{|X \\cap Y|}{|X| + |Y| - |X \\cap Y| }.$$\n",
    "\n",
    "\n",
    "-----\n",
    "- P. Jaccard, The distribution of the flora in the alpine zone, *New Phytol.*, vol. 11, pp. 37--50, Feb. 1912. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Setting Up Options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# configuring figure size\n",
    "options(repr.plot.width = 6, repr.plot.height = 6)\n",
    "\n",
    "# colorblind-friendly color palette\n",
    "cbPalette <- c('#999999', '#E69F00', '#56B4E9', '#009E73', \n",
    "               '#0072B2', '#D55E00', '#CC79A7', '#F0E442')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Loading Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "attributes": {
     "": "",
     "classes": [],
     "id": ""
    },
    "scrolled": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "library(tidyverse)\n",
    "\n",
    "library(latex2exp)\n",
    "\n",
    "library(maps) \n",
    "\n",
    "library(gganimate)\n",
    "\n",
    "library(gifski)\n",
    "\n",
    "require(viridis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Data Sets\n",
    "\n",
    "- <a href=\"https://www.ethnologue.com/ethnoblog/gary-simons/welcome-22nd-edition\">Ethnologue Global Dataset (22nd Edition)</a>\n",
    "\n",
    "    - `Table_of_Countries.tab`\n",
    "    - `Table_of_LICs.tab`\n",
    "- <a href=\"https://glottolog.org/\">Glottolog 4.6</a>\n",
    "    - `languoid.csv`\n",
    "    \n",
    "----\n",
    "- Eberhard, David M., Gary F. Simons, and Charles D. Fennig (eds.). 2019. Ethnologue: Languages of the World. Twenty-second edition. Dallas, Texas: SIL International.\n",
    "\n",
    "- Hammarström, Harald & Forkel, Robert & Haspelmath, Martin & Bank, Sebastian. 2022. Glottolog 4.6. Leipzig: Max Planck Institute for Evolutionary Anthropology. <a href=\"https://doi.org/10.5281/zenodo.6578297\">https://doi.org/10.5281/zenodo.6578297</a>, Accessed on 2022-08-19."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "# Loading and Cleaning Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Ethnologue Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "suppressMessages(countries_tbl <- read_tsv('Table_of_Countries.tab'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "#suppressMessages(languages_tbl <- read_csv('Table_of_Languages.csv')) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "names(countries_tbl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Changing the ISO_639 code for Man Nan Chinese from `nan` to `nzn`, because Python reads this as `NaN` (not a number, used for missing values as well). This might create issues with `NaN` (not a number, different than `NA` which is used for missing values) in R as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "#languages_tbl <- languages_tbl %>% \n",
    "#    mutate_at(vars(Uninverted_Name), \n",
    "#              list( ~ ifelse(ISO_639 == 'nan', 'Min Nan Chinese', .)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "suppressMessages(lic_tbl <- read_tsv('Table_of_LICs.tab'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "names(lic_tbl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "lic_tbl <- lic_tbl %>% \n",
    "    mutate_at(vars(All_Users, L1_Users, L2_Users), \n",
    "              list(~ ifelse(is.na(.), 0, .)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "lic_tbl <- lic_tbl %>% \n",
    "    mutate_at(vars(ISO_639), \n",
    "              list( ~ ifelse(Language_Name == 'Chinese, Min Nan', 'nzn', .)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "countries <- c('AE', 'AF', 'AM', 'AZ', 'BD', \n",
    "               'BH', 'BT', 'CY', 'GE', 'IL', \n",
    "               'IN', 'IQ', 'IR', 'JO', 'KG', \n",
    "               'KW', 'KZ', 'LB', 'LK', 'NP', \n",
    "               'OM', 'PK', 'PS', 'QA', 'SA', \n",
    "               'SY', 'TM', 'TR', 'TJ', 'UZ', 'YE')\n",
    "\n",
    "N <- length(countries)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "country_names_1 <- c('United Arab Emirates', 'Afghanistan', 'Armenia',\n",
    "                     'Azerbaijan','Bangladesh','Bahrain',\n",
    "                     'Bhutan','Cyprus','Georgia',\n",
    "                     'Israel','India','Iraq',\n",
    "                     'Iran','Jordan','Kyrgyzstan',\n",
    "                     'Kuwait','Kazakhstan','Lebanon',\n",
    "                     'Sri Lanka','Nepal','Oman',\n",
    "                     'Pakistan','Palestine','Qatar',\n",
    "                     'Saudi Arabia','Syria','Tajikistan',\n",
    "                     'Turkmenistan','Turkey','Uzbekistan','Yemen')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "main_tbl <- lic_tbl %>% filter(Country_Code %in% countries)\n",
    "\n",
    "main_tbl <- main_tbl %>% mutate(Region_New = case_when(\n",
    "    Country_Code %in% c('CY', 'IL', 'JO', 'LB', 'PS', 'SY') ~ 'Eastern Mediterranean',\n",
    "    Country_Code %in% c('AM', 'AZ', 'GE', 'IQ', 'IR', 'TR') ~ 'Western Asia',\n",
    "    Country_Code %in% c('AE', 'BH', 'KW', 'OM', 'QA', 'SA', 'YE') ~ 'the Arabian Peninsula',  \n",
    "    Country_Code %in% c('BD', 'BT', 'IN', 'LK', 'NP', 'PK') ~ 'Southern Asia',\n",
    "    Country_Code %in% c('AF', 'KG', 'KZ', 'TM', 'TJ', 'UZ') ~ 'Central Asia'))\n",
    "\n",
    "main_tbl <- main_tbl %>% \n",
    "    mutate_at(vars(Family), \n",
    "              list( ~ ifelse(. == 'Abkhaz-Adyghe', 'Abkhaz-Adyge', .)))\n",
    "\n",
    "main_tbl <- main_tbl %>% \n",
    "    mutate_at(vars(Family), \n",
    "              list( ~ ifelse(. == 'Austro-Asiatic', 'Austroasiatic', .)))\n",
    "\n",
    "main_tbl <- main_tbl %>% \n",
    "    mutate_at(vars(Family), \n",
    "              list( ~ ifelse(. == 'Mongolic', 'Mongolic-Khitan', .)))\n",
    "\n",
    "main_tbl <- main_tbl %>% \n",
    "    mutate_at(vars(Family), \n",
    "              list( ~ ifelse(. == 'Niger-Congo', 'Atlantic-Congo', .)))\n",
    "\n",
    "main_tbl <- main_tbl %>% \n",
    "    mutate_at(vars(Family), \n",
    "              list(~ ifelse(ISO_639 == 'bha', 'Indo-European', .)))\n",
    "\n",
    "main_tbl <- main_tbl %>% \n",
    "    mutate_at(vars(Family), \n",
    "              list(~ ifelse(ISO_639 %in% c('aml', 'asr', 'bfw', 'bix', 'biy', \n",
    "                                           'caq', 'cdz', 'crv', 'ekl', 'gaq', \n",
    "                                           'gbj', 'hoc', 'jun', 'juy', 'kfp', \n",
    "                                           'kfq', 'kha', 'khr', 'ksz', 'lyg', \n",
    "                                           'mef', 'mjx', 'mmj', 'ncb', 'nik', \n",
    "                                           'pbv', 'pcj', 'sat', 'srb', 'tef', \n",
    "                                           'trd', 'unr'), 'Austroasiatic', .)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "compact_df <- main_tbl %>% filter(L1_Users != 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# total population for each language in all the countries\n",
    "suppressMessages(totalpop_df <- compact_df %>% \n",
    "                     group_by(ISO_639) %>% \n",
    "                     summarise(Total_L1 = sum(L1_Users)) %>% \n",
    "                     ungroup())\n",
    "\n",
    "names(totalpop_df) <- c('ISO_639', 'Total_L1_Users')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "unique_lang <- unique(compact_df$ISO_639)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "families <- compact_df %>% select(Family) %>% distinct()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "families_ex <- c('Andamanese', 'Creole', 'Mixed language', 'Language isolate', 'Unclassified')\n",
    "\n",
    "families <- families %>% filter(!(Family %in% families_ex))\n",
    "\n",
    "fam_list <- do.call(rbind.data.frame, list('Jarawa-Onge')) \n",
    "\n",
    "names(fam_list) <- c('Family')\n",
    "\n",
    "families <- rbind(families, fam_list)\n",
    "\n",
    "length(families$Family)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Subsets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# can exclude\n",
    "countries_tbl %>% \n",
    "    filter(Country_Code %in% countries) %>%\n",
    "    select(Country_Code, Country_Name, Included, Diversity) %>% \n",
    "    glimpse()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# can exclude\n",
    "#languages_tbl %>% \n",
    "#    filter(Country_Code %in% countries) %>%\n",
    "#    select(ISO_639, Uninverted_Name, Country_Code, L1_Users, Family) %>%\n",
    "#    glimpse()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# must include\n",
    "lic_tbl %>% \n",
    "    filter(Country_Code %in% countries) %>%\n",
    "    select(ISO_639, Country_Code, L1_Users, Family) %>%\n",
    "    glimpse()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Glotolog Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "suppressMessages(languoid <- read_csv('languoid.csv'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "names(languoid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "languoid <- languoid %>% \n",
    "    mutate_at(vars(name), \n",
    "              list( ~ ifelse(iso639P3code == 'nan', 'Min Nan Chinese', .)))\n",
    "\n",
    "languoid <- languoid %>% \n",
    "    mutate_at(vars(iso639P3code), \n",
    "              list( ~ ifelse(iso639P3code == 'nan', 'nzn', .)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "languoid_compact <- languoid %>% \n",
    "    select(iso639P3code, \n",
    "           id, \n",
    "           name, \n",
    "           family_id, \n",
    "           parent_id, \n",
    "           level, \n",
    "           bookkeeping, \n",
    "           country_ids, \n",
    "           child_dialect_count)\n",
    "\n",
    "names(languoid_compact) <- c('ISO_639', \n",
    "                             'Glottocode', \n",
    "                             'Glotto_Name', \n",
    "                             'Glotto_Family', \n",
    "                             'Glotto_Parent', \n",
    "                             'Glotto_Level', \n",
    "                             'Bookkeeping', \n",
    "                             'Country_IDs', \n",
    "                             'Child_Dialect_Count')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Merging Data Sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# dataframe that contains the total population for each language in all the countries\n",
    "suppressMessages(compact_df <- compact_df %>% \n",
    "                     left_join(totalpop_df))\n",
    "\n",
    "suppressMessages(compact_df <- compact_df %>% \n",
    "                     inner_join(languoid_compact))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "glimpse(compact_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "names(compact_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Countries\n",
    "\n",
    "AE: United Arab Emirates; AF: Afghanistan; AM: Armenia; AZ: Azerbaijan; \n",
    "\n",
    "BD: Bangladesh; BH: Bahrain; BT: Bhutan; CY: Cyprus;\n",
    "\n",
    "GE: Georgia; IN: India; IR: Iran; IQ: Iraq;\n",
    "\n",
    "IL: Isreal; JO: Jordan; KZ: Kazakhstan; KW: Kuwait;\n",
    "\n",
    "KG: Kyrgyzstan; LB: Lebanon; LK: Sri Lanka; NP: Nepal;\n",
    "\n",
    "OM: Oman; PK: Pakistan; PS: Palestine; QA: Qatar;\n",
    "\n",
    "SA: Saudi Arabia; SY: Syria; TJ: Tajikistan; TM: Turkmenistan;\n",
    "\n",
    "TR: Turkey; UZ: Uzbekistan; YE: Yemen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Regions\n",
    "\n",
    "- The Arabian Peninsula: AE, BH, KW, OM, QA, SA, YE\n",
    "- Central Asia: AF, KG, KZ, TM, TJ, UZ\n",
    "- Eastern Mediterranean: CY, IL, JO, LB, PS, SY\n",
    "- Southern Asia: BD, BT, IN, LK, NP, PK\n",
    "- Western Asia: AM, AZ, GE, IQ, IR, TR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "countries_map <- map_data('world', region = country_names_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "countries_map <- countries_map %>% \n",
    "    mutate(Region = case_when(\n",
    "        region %in% c('Cyprus', 'Israel', 'Jordan', 'Lebanon', 'Palestine', 'Syria') ~ 'Eastern Mediterranean',\n",
    "        region %in% c('Armenia', 'Azerbaijan', 'Georgia', 'Iraq', 'Iran', 'Turkey') ~ 'Western Asia',\n",
    "        region %in% c('United Arab Emirates', 'Bahrain', 'Kuwait', 'Oman', 'Qatar', 'Saudi Arabia', 'Yemen') ~ 'Arabian Peninsula',  \n",
    "        region %in% c('Bangladesh', 'Bhutan', 'India', 'Sri Lanka', 'Nepal', 'Pakistan') ~ 'Southern Asia',\n",
    "        region %in% c('Afghanistan', 'Kyrgyzstan', 'Kazakhstan', 'Turkmenistan', 'Tajikistan', 'Uzbekistan') ~ 'Central Asia'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(countries_map, aes(x = long, y = lat, fill = Region)) +\n",
    "    geom_polygon(aes(group = group), color = 'white') +\n",
    "    scale_fill_manual(values = cbPalette) +\n",
    "    coord_map('polyconic') +\n",
    "    xlab('Longitude') +\n",
    "    ylab('Latitude')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Diversities\n",
    "\n",
    "## Richness, Shannon, Greenberg, and HCDT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# defining ranges\n",
    "range_1 <- seq(.001, 1, .01)\n",
    "\n",
    "range_2 <- seq(1.001, 5, .01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "suppressMessages(country_df <- compact_df %>% \n",
    "                     group_by(Country_Code) %>%\n",
    "                     mutate(Pop = sum(L1_Users), Prop = L1_Users / Pop) %>% \n",
    "                     ungroup())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# effective numbers for HCDT index, for 0 < q < 1\n",
    "main_1 <- NULL\n",
    "\n",
    "for (q in range_1) {\n",
    "    \n",
    "    suppressMessages(df <- country_df %>% \n",
    "                         group_by(Region_New, Country_Code, Country_Name) %>%\n",
    "                         summarise(Diversity = sum(Prop ** q) ** (1/ (1 - q)), q = q) %>%\n",
    "                         ungroup())\n",
    "                     \n",
    "    main_1 <- rbind(main_1, df)\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "for (region in unique(main_tbl$Region_New)) {\n",
    "    \n",
    "    print(main_1 %>% \n",
    "              filter(Region_New == region) %>%\n",
    "              ggplot(aes(x = q, y = log10(Diversity))) + \n",
    "              geom_line(aes(linetype = Country_Code, color = Country_Code)) + \n",
    "              scale_color_manual(values = cbPalette) +\n",
    "              labs(title = paste('Effective Numbers for Countries in', region),\n",
    "                   subtitle = 'Using HCDT Index',\n",
    "                   linetype = 'Country Code',\n",
    "                   color = 'Country Code') +\n",
    "              xlab('q') + \n",
    "              ylab(TeX('Effective Number on $\\\\log_{10}$-scale')) +\n",
    "              scale_y_continuous(limits = c(0, 3)) +\n",
    "              theme(legend.text = element_text(size = 12, face = 'bold')))\n",
    "        \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# effective numbers for HCDT index, for 1 < q < 5\n",
    "main_2 <- NULL\n",
    "\n",
    "for (q in range_2) {\n",
    "    \n",
    "    suppressMessages(df <- country_df %>% \n",
    "                         group_by(Region_New, Country_Code, Country_Name) %>% \n",
    "                         summarise(Diversity = sum(Prop ** q) ** (1/ (1 - q)), q = q) %>%\n",
    "                         ungroup())   \n",
    "    \n",
    "    main_2 <- rbind(main_2, df)\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "for (region in unique(main_tbl$Region_New)) {\n",
    "    \n",
    "    print(main_2 %>% \n",
    "              filter(Region_New == region) %>%\n",
    "              ggplot(aes(x = q, y = log10(Diversity))) + \n",
    "              geom_line(aes(linetype = Country_Code, color = Country_Code)) + \n",
    "              geom_vline(xintercept = 2, color = 'red') + \n",
    "              scale_color_manual(values = cbPalette) +\n",
    "              labs(title = paste('Effective Numbers for Countries in', region),\n",
    "                   subtitle = 'Using HCDT index',\n",
    "                   linetype = 'Country Code',\n",
    "                   color = 'Country Code') +\n",
    "              xlab('q') +\n",
    "              ylab(TeX('Effective Number on $\\\\log_{10}$-scale')) + \n",
    "              scale_y_continuous(limits = c(0, 3)) +\n",
    "              theme(legend.text = element_text(size = 12, face = 'bold')))\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# effective numbers for HCDT index, for 0 < q < 5\n",
    "main <- rbind(main_1, main_2)\n",
    "\n",
    "for (region in unique(main_tbl$Region_New)) {\n",
    "    \n",
    "    print(main %>% \n",
    "              filter(Region_New == region) %>%\n",
    "              ggplot(aes(x = q, y = log10(Diversity))) + \n",
    "              geom_line(aes(linetype = Country_Code, color = Country_Code)) + \n",
    "              geom_vline(xintercept = 1, color = 'yellow') + \n",
    "              geom_vline(xintercept = 2, color = 'red') + \n",
    "              scale_color_manual(values = cbPalette) +\n",
    "              labs(title = paste('Effective Numbers for Countries in', region),\n",
    "                   subtitle = 'Using HCDT index',\n",
    "                   linetype = 'Country Code',\n",
    "                   color = 'Country Code') +\n",
    "              xlab('q') +\n",
    "              ylab(TeX('Effective Number on $\\\\log_{10}$-scale')) + \n",
    "              scale_y_continuous(limits = c(0, 3)) +\n",
    "              theme(legend.text = element_text(size = 12, face = 'bold')))\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "main_1 %>% \n",
    "    ggplot(aes(x = q, y = log10(Diversity))) + \n",
    "    geom_line(color = 'blue') + \n",
    "    facet_wrap(vars(Country_Code)) + \n",
    "    scale_color_manual(values = cbPalette) +\n",
    "    labs(title = 'Effective Number Curves') +\n",
    "    xlab('q') +\n",
    "    ylab('Effective Number on Logarithmic-scale') +\n",
    "    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "main_2 %>% \n",
    "    ggplot(aes(x = q, y = log10(Diversity))) + \n",
    "    geom_line(color = 'blue') + \n",
    "    geom_vline(xintercept = 2, color = 'red') + \n",
    "    facet_wrap(vars(Country_Code)) + \n",
    "    scale_color_manual(values = cbPalette) +\n",
    "    labs(title = 'Effective Number Curves') +\n",
    "    xlab('q') +\n",
    "    ylab('Effective Number on Logarithmic-scale') +\n",
    "    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "main %>% \n",
    "    ggplot(aes(x = q, y = log10(Diversity))) + \n",
    "    geom_line(color = 'blue') + \n",
    "    geom_vline(xintercept = 1, color = 'yellow') + \n",
    "    geom_vline(xintercept = 2, color = 'red') + \n",
    "    facet_wrap(vars(Country_Code)) + \n",
    "    scale_color_manual(values = cbPalette) +\n",
    "    labs(title = 'Effective Number Curves') +\n",
    "    xlab('q') +\n",
    "    ylab('Effective Number on Logarithmic-scale') +\n",
    "    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "countries_tbl <- countries_tbl %>% \n",
    "    mutate(Greenberg = 1 / (1 - Diversity))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "suppressMessages(countries_agg_df <- country_df %>% \n",
    "                     group_by(Country_Code, Country_Name) %>%\n",
    "                     summarise(Richness = n(),\n",
    "                               Shannon = exp(-sum(Prop * log(Prop))),\n",
    "                               Greenberg = 1 / sum(Prop ** 2)) %>%\n",
    "                     ungroup())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_richness_ord <- countries_agg_df %>% \n",
    "    select(Country_Code, Richness) %>%\n",
    "    arrange(desc(Richness)) %>%\n",
    "    mutate(rank = row_number())\n",
    "\n",
    "df_richness_ord"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_shannon_ord <- countries_agg_df %>% \n",
    "    select(Country_Code, Shannon) %>%\n",
    "    arrange(desc(Shannon)) %>%\n",
    "    mutate(rank = row_number())\n",
    "\n",
    "df_shannon_ord"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_greenberg_ord <- countries_agg_df %>% \n",
    "    select(Country_Code, Greenberg) %>%\n",
    "    arrange(desc(Greenberg)) %>%\n",
    "    mutate(rank = row_number())\n",
    "\n",
    "df_greenberg_ord"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_richness_ord %>% \n",
    "    select(-Richness) %>% \n",
    "    left_join(., df_shannon_ord, by = 'Country_Code') %>%\n",
    "    select(-Shannon) %>%\n",
    "    left_join(., df_greenberg_ord, by = 'Country_Code') %>%\n",
    "    select(-Greenberg) %>%\n",
    "    mutate(Rank = 1/3*(rank + rank.x + rank.y)) %>%\n",
    "    select(Country_Code, Rank) %>%\n",
    "    arrange(Rank)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "tab1 <- left_join(countries_agg_df, countries_tbl, by = c('Country_Code', 'Country_Name')) %>% \n",
    "    select(Country_Code, Country_Name, Richness, Included, Greenberg.x, Greenberg.y) %>%\n",
    "    rename('Richness_Formulas' = 'Richness',\n",
    "           'Richness_SIL' = 'Included',\n",
    "           'Greenberg_Formulas' = 'Greenberg.x', \n",
    "           'Greenberg_SIL' = 'Greenberg.y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "country_names <- countries_agg_df %>% \n",
    "    select(Country_Name) %>%\n",
    "    unlist(use.names = FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "countries_map <- map_data('world', region = country_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "countries_map_df <- left_join(countries_map, \n",
    "                              countries_agg_df, \n",
    "                              by = c('region' = 'Country_Name'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(countries_map_df, aes(x = long, y = lat)) +\n",
    "    geom_polygon(aes(group = group, fill = log10(Richness)), color = 'white') +\n",
    "    scale_fill_viridis_c(limits = c(0, 2.7)) +\n",
    "    coord_map('polyconic') +\n",
    "    xlab('Longitude') +\n",
    "    ylab('Latitude')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(countries_map_df, aes(x = long, y = lat)) +\n",
    "    geom_polygon(aes(group = group, fill = log10(Shannon)), color = 'white') +\n",
    "    scale_fill_viridis_c(limits = c(0, 2.7)) +\n",
    "    coord_map('polyconic') +\n",
    "    xlab('Longitude') +\n",
    "    ylab('Latitude')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(countries_map_df, aes(x = long, y = lat)) +\n",
    "    geom_polygon(aes(group = group, fill = log10(Greenberg)), color = 'white') +\n",
    "    scale_fill_viridis_c(limits = c(0, 2.7)) +\n",
    "    coord_map('polyconic') +\n",
    "    xlab('Longitude') +\n",
    "    ylab('Latitude')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# https://conservancy.umn.edu/bitstream/handle/11299/220339/time-maps-tutorial-v2.html?sequence=3&isAllowed=y\n",
    "base_map <- ggplot(data = countries_map_df, aes(x = long, y = lat, group = group)) +\n",
    "    geom_polygon(fill = 'white', color = 'white') +\n",
    "    coord_map('polyconic') +\n",
    "    xlab('Longitude') +\n",
    "    ylab('Latitude')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "base_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "main_animate <- left_join(countries_map, \n",
    "                              main, \n",
    "                              by = c('region' = 'Country_Name')) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "range <- seq(.001, 5, .2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "filtered_main_animate <- main_animate %>%\n",
    "    filter(q %in% range)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "map_with_data <- base_map +\n",
    "    geom_polygon(data = filtered_main_animate, aes(x = long, y = lat, \n",
    "                                                   fill = log10(Diversity), group = group), color = 'white') +\n",
    "    scale_fill_viridis_c()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "map_with_animation <- map_with_data + \n",
    "    transition_time(filtered_main_animate$q) +\n",
    "    ggtitle('q: {frame_time}', subtitle = 'Frame {frame} of {nframes}')\n",
    "\n",
    "num_qs <- length(range)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "suppressMessages(animate(map_with_animation, nframes = num_qs, fps = 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Gamma Diversities "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Unweighted Gamma Diversities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "giml_df <- compact_df %>% \n",
    "    group_by(Country_Code) %>% \n",
    "    mutate(Country_Pop = sum(L1_Users)) %>% \n",
    "    ungroup()\n",
    "\n",
    "giml_u_m <- giml_df %>% \n",
    "    mutate(Total_Pop = sum(L1_Users), \n",
    "           Prop = L1_Users / Country_Pop, \n",
    "           Weight = 1 / N)\n",
    "\n",
    "giml_u_m <- giml_u_m %>% \n",
    "    select(-Country_Code) %>%\n",
    "    group_by(ISO_639) %>%\n",
    "    mutate(Weighted_Prop = sum(Prop * Weight)) %>% \n",
    "    ungroup() %>%\n",
    "    select(ISO_639, Weighted_Prop) %>%\n",
    "    unique()\n",
    "\n",
    "giml_u_m <- giml_u_m %>% \n",
    "    summarise(Gamma_Richness = n(),\n",
    "              Gamma_Shannon = exp(-sum(Weighted_Prop * log(Weighted_Prop))),\n",
    "              Gamma_Greenberg = 1 / sum(Weighted_Prop ** 2))\n",
    "\n",
    "giml_u_m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Weighted Gamma Diversities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "giml_w_m <- giml_df %>% \n",
    "    mutate(Total_Pop = sum(L1_Users), \n",
    "           Prop = L1_Users / Country_Pop, \n",
    "           Weight = Country_Pop / Total_Pop)\n",
    "\n",
    "giml_w_m <- giml_w_m %>%  \n",
    "    select(-Country_Code) %>%\n",
    "    group_by(ISO_639) %>%\n",
    "    mutate(Weighted_Prop = sum(Prop * Weight)) %>% \n",
    "    ungroup() %>%\n",
    "    select(ISO_639, Weighted_Prop) %>%\n",
    "    unique()\n",
    "\n",
    "giml_w_m <- giml_w_m %>% \n",
    "    summarise(Gamma_Richness = n(),\n",
    "              Gamma_Shannon = exp(-sum(Weighted_Prop * log(Weighted_Prop))),\n",
    "              Gamma_Greenberg = 1 / sum(Weighted_Prop ** 2))\n",
    "\n",
    "giml_w_m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Alpha Diversities "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Unweighted Alpha Diversities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "alep_df <- compact_df %>% mutate(Total_Pop = sum(L1_Users))\n",
    "\n",
    "alep_u_m <- alep_df %>% \n",
    "    group_by(Country_Code) %>% \n",
    "    mutate(Pop = sum(L1_Users), \n",
    "           Prop = L1_Users / Pop,\n",
    "           Total_Prop = L1_Users / Total_Pop,\n",
    "           Weight = 1 / N)\n",
    "\n",
    "suppressMessages(alep_u_m <- alep_u_m %>% \n",
    "                     summarise(Richness = n(), \n",
    "                               Shannon = -sum(Prop * log(Prop)),\n",
    "                               Greenberg = sum(Prop ** 2)) %>%\n",
    "                     ungroup() %>% \n",
    "                     unique())\n",
    "\n",
    "alep_u_m <- alep_u_m %>% \n",
    "    summarise(Alpha_Richness = mean(Richness), \n",
    "              Alpha_Shannon = exp(mean(Shannon)), \n",
    "              Alpha_Greenberg = 1 / mean(Greenberg))\n",
    "\n",
    "alep_u_m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Weighted Alpha Diversities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "alep_df <- compact_df %>% mutate(Total_Pop = sum(L1_Users))\n",
    "\n",
    "alep_w_m <- alep_df %>% \n",
    "    group_by(Country_Code) %>% \n",
    "    mutate(Country_Pop = sum(L1_Users), \n",
    "           Prop = L1_Users / Country_Pop,\n",
    "           Weight = Country_Pop / Total_Pop,\n",
    "           Richness = n(),\n",
    "           Shannon = -sum(Prop * log(Prop)),\n",
    "           Greenberg = sum(Prop ** 2)) %>%\n",
    "    ungroup() %>%\n",
    "    select(Country_Code, \n",
    "           Richness, \n",
    "           Shannon, \n",
    "           Greenberg, \n",
    "           Country_Pop, \n",
    "           Total_Pop, \n",
    "           Weight) %>%\n",
    "    unique()\n",
    "        \n",
    "suppressMessages(alep_w_m <- alep_w_m %>%\n",
    "                     summarise(Alpha_Richness = mean(Richness),\n",
    "                               Alpha_Shannon = exp(sum(Shannon * Weight)),\n",
    "                               Alpha_Greenberg = sum(Weight ** 2) / sum(Weight ** 2 * Greenberg )))\n",
    "\n",
    "alep_w_m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Beta Diversities "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Unweighted Beta Diversities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "bet_u_m <- giml_u_m / alep_u_m \n",
    "\n",
    "names(bet_u_m) <- c('Beta_Richness', 'Beta_Shannon', 'Beta_Greenberg')\n",
    "\n",
    "bet_u_m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Weighted Beta Diversities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "bet_w_m <- giml_w_m / alep_w_m\n",
    "\n",
    "names(bet_w_m) <- c('Beta_Richness', 'Beta_Shannon', 'Beta_Greenberg')\n",
    "\n",
    "bet_w_m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Visualizing Diversities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# defining two ranges: (0, 1) and (1, 10)\n",
    "range_1 <- seq(0.001, 1, .01)\n",
    "\n",
    "range_2 <- seq(1.001, 5, .01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Visualizing Unweighted Gamma Diversity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "qsum1_g_u_m <- NULL\n",
    "\n",
    "giml <- giml_df %>% \n",
    "    mutate(Total_Pop = sum(L1_Users), \n",
    "           Prop = L1_Users / Country_Pop, \n",
    "           Weight = 1 / N)\n",
    "\n",
    "for (q in range_1) {\n",
    "  \n",
    "    df2 <- NULL\n",
    "    \n",
    "    df2 <- giml %>%  \n",
    "        select(-Country_Code) %>%\n",
    "        group_by(ISO_639) %>%\n",
    "        mutate(Weighted_Prop = sum(Prop * Weight)) %>% \n",
    "        ungroup() %>%\n",
    "        select(ISO_639, Weighted_Prop) %>%\n",
    "        unique()\n",
    "    \n",
    "    df2 <- df2 %>% \n",
    "        summarise(Giml_Manual = sum(Weighted_Prop ** q) ** (1/ (1 - q)), \n",
    "                  q = q)\n",
    "    \n",
    "    qsum1_g_u_m <- rbind(qsum1_g_u_m, df2)\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print(qsum1_g_u_m %>% \n",
    "          ggplot(aes(x = q, y = Giml_Manual)) + \n",
    "          geom_line( color = 'blue') +\n",
    "          geom_hline(yintercept = giml_u_m$Gamma_Richness, color = 'green') +\n",
    "          geom_vline(xintercept = 0, color = 'green') + \n",
    "          geom_hline(yintercept = giml_u_m$Gamma_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          labs(title = 'Unweighted Gamma Diversity Computed Using Fromulas ') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\gamma$ Diversity'))) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "qsum2_g_u_m <- NULL\n",
    "\n",
    "giml <- giml_df %>% \n",
    "    mutate(Total_Pop = sum(L1_Users), \n",
    "           Prop = L1_Users / Country_Pop,\n",
    "           Weight = 1 / N)\n",
    "\n",
    "for (q in range_2) {\n",
    "  \n",
    "    df2 <- NULL\n",
    "    \n",
    "    df2 <- giml %>% \n",
    "        select(-Country_Code) %>%\n",
    "        group_by(ISO_639) %>%\n",
    "        mutate(Weighted_Prop = sum(Prop * Weight)) %>% \n",
    "        ungroup() %>%\n",
    "        select(ISO_639, Weighted_Prop) %>%\n",
    "        unique()\n",
    "    \n",
    "    df2 <- df2 %>% \n",
    "        summarise(Giml_Manual = sum(Weighted_Prop ** q) ** (1 / (1 - q)),\n",
    "                  q = q)\n",
    "    \n",
    "    qsum2_g_u_m <- rbind(qsum2_g_u_m, df2)\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print(qsum2_g_u_m %>% \n",
    "          ggplot(aes(x = q, y = Giml_Manual)) + \n",
    "          geom_line( color = 'blue') +\n",
    "          geom_hline(yintercept = giml_u_m$Gamma_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          geom_hline(yintercept = giml_u_m$Gamma_Greenberg, color = 'red') +\n",
    "          geom_vline(xintercept = 2, color = 'red') + \n",
    "          labs(title = 'Unweighted Gamma Diversity Computed Using Fromulas') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\gamma$ Diversity'))) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "qsum_g_u_m <- rbind(qsum1_g_u_m, qsum2_g_u_m)\n",
    "\n",
    "print(qsum_g_u_m %>% \n",
    "          ggplot(aes(x = q, y = Giml_Manual)) + \n",
    "          geom_line( color = 'blue') +\n",
    "          geom_hline(yintercept = giml_u_m$Gamma_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          geom_hline(yintercept = giml_u_m$Gamma_Greenberg, color = 'red') +\n",
    "          geom_vline(xintercept = 2, color = 'red') + \n",
    "          labs(title = 'Unweighted Gamma Diversity Computed Using Fromulas') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\gamma$ Diversity'))) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Visualizing Weighted Gamma Diversity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "qsum1_g_w_m <- NULL\n",
    "\n",
    "giml <- giml_df %>% \n",
    "    mutate(Total_Pop = sum(L1_Users), \n",
    "           Prop = L1_Users / Country_Pop, \n",
    "           Weight = Country_Pop / Total_Pop)\n",
    "\n",
    "for (q in range_1) {\n",
    "  \n",
    "    df2 <- NULL\n",
    "    \n",
    "    df2 <- giml %>%  \n",
    "        select(-Country_Code) %>%\n",
    "        group_by(ISO_639) %>%\n",
    "        mutate(Weighted_Prop = sum(Prop * Weight)) %>% \n",
    "        ungroup() %>%\n",
    "        select(ISO_639, Weighted_Prop) %>%\n",
    "        unique()\n",
    "    \n",
    "    df2 <- df2 %>% \n",
    "        summarise(Giml_Manual = sum(Weighted_Prop ** q) ** (1/ (1 - q)), \n",
    "                  q = q)\n",
    "    \n",
    "    qsum1_g_w_m <- rbind(qsum1_g_w_m, df2)\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print(qsum1_g_w_m %>% \n",
    "          ggplot(aes(x = q, y = Giml_Manual)) + \n",
    "          geom_line( color = 'blue') +\n",
    "          geom_hline(yintercept = giml_w_m$Gamma_Richness, color = 'green') +\n",
    "          geom_vline(xintercept = 0, color = 'green') + \n",
    "          geom_hline(yintercept = giml_w_m$Gamma_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          labs(title = 'Weighted Gamma Diversity Computed Using Fromulas ') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\gamma$ Diversity'))) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "qsum2_g_w_m <- NULL\n",
    "\n",
    "giml <- giml_df %>% \n",
    "    mutate(Total_Pop = sum(L1_Users), \n",
    "           Prop = L1_Users / Country_Pop, \n",
    "           Weight = Country_Pop / Total_Pop)\n",
    "\n",
    "for (q in range_2) {\n",
    "  \n",
    "    df2 <- NULL\n",
    "    \n",
    "    df2 <- giml %>% \n",
    "        select(-Country_Code) %>%\n",
    "        group_by(ISO_639) %>%\n",
    "        mutate(Weighted_Prop = sum(Prop * Weight)) %>% \n",
    "        ungroup() %>%\n",
    "        select(ISO_639, Weighted_Prop) %>%\n",
    "        unique()\n",
    "    \n",
    "    df2 <- df2 %>% \n",
    "        summarise(Giml_Manual = sum(Weighted_Prop ** q) ** (1 / (1 - q)),\n",
    "                  q = q)\n",
    "    \n",
    "    qsum2_g_w_m <- rbind(qsum2_g_w_m, df2)\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print(qsum2_g_w_m %>% \n",
    "          ggplot(aes(x = q, y = Giml_Manual)) + \n",
    "          geom_line(color = 'blue') +\n",
    "          geom_hline(yintercept = giml_w_m$Gamma_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          geom_hline(yintercept = giml_w_m$Gamma_Greenberg, color = 'red') +\n",
    "          geom_vline(xintercept = 2, color = 'red') + \n",
    "          labs(title = 'Weighted Gamma Diversity Computed Using Fromulas') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\gamma$ Diversity'))) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "qsum_g_w_m <- rbind(qsum1_g_w_m, qsum2_g_w_m)\n",
    "\n",
    "print(qsum_g_w_m %>% \n",
    "          ggplot(aes(x = q, y = Giml_Manual)) + \n",
    "          geom_line( color = 'blue') +\n",
    "          geom_hline(yintercept = giml_w_m$Gamma_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          geom_hline(yintercept = giml_w_m$Gamma_Greenberg, color = 'red') +\n",
    "          geom_vline(xintercept = 2, color = 'red') + \n",
    "          labs(title = 'Weighted Gamma Diversity Computed Using Fromulas') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\gamma$ Diversity'))) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "alep <- compact_df %>% mutate(Total_Pop = sum(L1_Users))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Visualizing Unweighted Alpha Diversity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "qsum1_a_u_m <- NULL\n",
    "\n",
    "for (q in range_1) {\n",
    "  \n",
    "    df2 <- NULL\n",
    "\n",
    "    df2 <- alep %>%\n",
    "        group_by(Country_Code) %>% \n",
    "        mutate(Country_Pop = sum(L1_Users), \n",
    "               Prop = L1_Users / Country_Pop,\n",
    "               Weight = 1 / N,\n",
    "               Smallqsum = sum((Prop * Weight) ** q),\n",
    "               q = q) %>%\n",
    "        ungroup() %>%\n",
    "        select(Country_Code, Smallqsum, Weight, q) %>%\n",
    "        unique()\n",
    "    \n",
    "    suppressMessages(df2 <- df2 %>% \n",
    "                         group_by(q) %>% \n",
    "                         summarise(Alep_Manual = (sum(Smallqsum) / sum(Weight ** q)) ** (1 / (1 - q)), \n",
    "                                   q = q))\n",
    "    \n",
    "    qsum1_a_u_m <- rbind(qsum1_a_u_m, df2)\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print(qsum1_a_u_m %>% \n",
    "          ggplot(aes(x = q, y = Alep_Manual)) + \n",
    "          geom_line( color = 'blue') +\n",
    "          geom_hline(yintercept = alep_u_m$Alpha_Richness, color = 'green') +\n",
    "          geom_vline(xintercept = 0, color = 'green') + \n",
    "          geom_hline(yintercept = alep_u_m$Alpha_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          labs(title = 'Unweighted Alpha Diversity Computed Using Fromulas ') +\n",
    "          xlab(TeX('$q$')) +\n",
    "         ylab(TeX('$\\\\alpha$ Diversity'))) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "qsum2_a_u_m <- NULL\n",
    "\n",
    "for (q in range_2) {\n",
    "    \n",
    "    df2 <- NULL\n",
    "    \n",
    "    df2 <- alep %>%\n",
    "        group_by(Country_Code) %>% \n",
    "        mutate(Country_Pop = sum(L1_Users), \n",
    "               Prop = L1_Users / Country_Pop,\n",
    "               Weight = 1 / N,\n",
    "               Smallqsum = sum((Prop * Weight) ** q),\n",
    "               q = q) %>%\n",
    "        ungroup() %>%\n",
    "        select(Country_Code, Smallqsum, Weight, q) %>%\n",
    "        unique()\n",
    "    \n",
    "    suppressMessages(df2 <- df2 %>% \n",
    "                         group_by(q) %>% \n",
    "                         summarise(Alep_Manual = (sum(Smallqsum) / sum(Weight ** q)) ** (1 / (1 - q)), \n",
    "                                   q = q))\n",
    "    \n",
    "    qsum2_a_u_m <- rbind(qsum2_a_u_m, df2)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print(qsum2_a_u_m %>% \n",
    "          ggplot(aes(x = q, y = Alep_Manual)) + \n",
    "          geom_line( color = 'blue') + \n",
    "          geom_hline(yintercept = alep_u_m$Alpha_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          geom_hline(yintercept = alep_u_m$Alpha_Greenberg, color = 'red') +\n",
    "          geom_vline(xintercept = 2, color = 'red') + \n",
    "          labs(title = 'Unweighted Alpha Diversity Computed Using Fromulas ') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\alpha$ Diversity'))) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "qsum_a_u_m <- rbind(qsum1_a_u_m, qsum2_a_u_m)\n",
    "\n",
    "print(qsum_a_u_m %>% \n",
    "          ggplot(aes(x = q, y = Alep_Manual)) + \n",
    "          geom_line( color = 'blue') +\n",
    "          geom_hline(yintercept = alep_u_m$Alpha_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          geom_hline(yintercept = alep_u_m$Alpha_Greenberg, color = 'red') +\n",
    "          geom_vline(xintercept = 2, color = 'red') + \n",
    "          labs(title = 'Unweighted Alpha Diversity Computed Using Fromulas') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\alpha$ Diversity'))) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Visualizing Weighted Alpha Diversity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "qsum1_a_w_m <- NULL\n",
    "\n",
    "for (q in range_1) {\n",
    "  \n",
    "    df2 <- NULL\n",
    "\n",
    "    df2 <- alep %>%\n",
    "        group_by(Country_Code) %>% \n",
    "        mutate(Country_Pop = sum(L1_Users), \n",
    "               Prop = L1_Users / Country_Pop,\n",
    "               Weight = Country_Pop / Total_Pop,\n",
    "               Smallqsum = sum((Prop * Weight) ** q),\n",
    "               q = q) %>%\n",
    "        ungroup() %>%\n",
    "        select(Country_Code, Smallqsum, Weight, q) %>%\n",
    "        unique()\n",
    "    \n",
    "    suppressMessages(df2 <- df2 %>% \n",
    "                         group_by(q) %>% \n",
    "                         summarise(Alep_Manual = (sum(Smallqsum) / sum(Weight ** q)) ** (1 / (1 - q)), \n",
    "                                   q = q))\n",
    "    \n",
    "    qsum1_a_w_m <- rbind(qsum1_a_w_m, df2)\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print(qsum1_a_w_m %>% \n",
    "          ggplot(aes(x = q, y = Alep_Manual)) + \n",
    "          geom_line( color = 'blue') + \n",
    "          geom_hline(yintercept = alep_w_m$Alpha_Richness, color = 'green') +\n",
    "          geom_vline(xintercept = 0, color = 'green') + \n",
    "          geom_hline(yintercept = alep_w_m$Alpha_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          labs(title = 'Weighted Alpha Diversity Computed Using Fromulas ') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\alpha$ Diversity'))) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "qsum2_a_w_m <- NULL\n",
    "\n",
    "for (q in range_2) {\n",
    "    \n",
    "    df2 <- NULL\n",
    "    \n",
    "    df2 <- alep %>%group_by(Country_Code) %>% \n",
    "        mutate(Country_Pop = sum(L1_Users), \n",
    "               Prop = L1_Users / Country_Pop,\n",
    "               Weight = Country_Pop / Total_Pop,\n",
    "               Smallqsum = sum((Prop * Weight) ** q),\n",
    "               q = q) %>%\n",
    "        ungroup() %>%\n",
    "        select(Country_Code, Smallqsum, Weight, q) %>%\n",
    "        unique()\n",
    "    \n",
    "    suppressMessages(df2 <- df2 %>%\n",
    "                         group_by(q) %>% \n",
    "                         summarise(Alep_Manual = (sum(Smallqsum) / sum(Weight ** q)) ** (1 / (1 - q)), \n",
    "                                   q = q))\n",
    "    \n",
    "    qsum2_a_w_m <- rbind(qsum2_a_w_m, df2)\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print(qsum2_a_w_m %>% \n",
    "          ggplot(aes(x = q, y = Alep_Manual)) + \n",
    "          geom_line( color = 'blue') + \n",
    "          geom_hline(yintercept = alep_w_m$Alpha_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          geom_hline(yintercept = alep_w_m$Alpha_Greenberg, color = 'red') +\n",
    "          geom_vline(xintercept = 2, color = 'red') + \n",
    "          labs(title = 'Weighted Alpha Diversity Computed Using Fromulas ') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\alpha$ Diversity')))  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "qsum_a_w_m <- rbind(qsum1_a_w_m, qsum2_a_w_m)\n",
    "\n",
    "print(qsum_a_w_m %>% \n",
    "          ggplot(aes(x = q, y = Alep_Manual)) + \n",
    "          geom_line( color = 'blue') +\n",
    "          geom_hline(yintercept = alep_w_m$Alpha_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          geom_hline(yintercept = alep_w_m$Alpha_Greenberg, color = 'red') +\n",
    "          geom_vline(xintercept = 2, color = 'red') + \n",
    "          labs(title = 'Weighted Alpha Diversity Computed Using Fromulas') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\alpha$ Diversity'))) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Visualizing Unweighted Beta Diversity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "qsum1_b_u_m <- inner_join(qsum1_g_u_m, qsum1_a_u_m, by = 'q')\n",
    "\n",
    "qsum1_b_u_m <- qsum1_b_u_m %>% mutate(Bet_Manual = Giml_Manual / Alep_Manual)\n",
    "\n",
    "print(qsum1_b_u_m %>%\n",
    "          ggplot(aes(x = q, y = Bet_Manual)) + \n",
    "          geom_line( color = 'blue') + \n",
    "          geom_hline(yintercept = bet_u_m$Beta_Richness, color = 'green') +\n",
    "          geom_vline(xintercept = 0, color = 'green') + \n",
    "          geom_hline(yintercept = bet_u_m$Beta_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          labs(title = 'Unweighted Beta Diversity Computed Using Fromulas') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\beta$ Diversity')))  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "qsum2_b_u_m <- inner_join(qsum2_g_u_m, qsum2_a_u_m, by = 'q')\n",
    "\n",
    "qsum2_b_u_m <- qsum2_b_u_m %>% mutate(Bet_Manual = Giml_Manual / Alep_Manual)\n",
    "\n",
    "print(qsum2_b_u_m %>% \n",
    "          ggplot(aes(x = q, y = Bet_Manual)) + \n",
    "          geom_line( color = 'blue') + \n",
    "          geom_hline(yintercept = bet_u_m$Beta_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          geom_hline(yintercept = bet_u_m$Beta_Greenberg, color = 'red') +\n",
    "          geom_vline(xintercept = 2, color = 'red') + \n",
    "          labs(title = 'Unweighted Beta Diversity Computed Using Fromulas') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\beta$ Diversity')))  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "qsum_b_u_m <- rbind(qsum1_b_u_m, qsum2_b_u_m)\n",
    "\n",
    "\n",
    "print(qsum_b_u_m %>% \n",
    "          ggplot(aes(x = q, y = Bet_Manual)) + \n",
    "          geom_line( color = 'blue') +\n",
    "          geom_hline(yintercept = bet_u_m$Beta_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          geom_hline(yintercept = bet_u_m$Beta_Greenberg, color = 'red') +\n",
    "          geom_vline(xintercept = 2, color = 'red') + \n",
    "          labs(title = 'Unweighted Beta Diversity Computed Using Fromulas') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\beta$ Diversity'))) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Visualizing Weighted Beta Diversity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "qsum1_b_w_m <- inner_join(qsum1_g_w_m, qsum1_a_w_m, by = 'q')\n",
    "\n",
    "qsum1_b_w_m <- qsum1_b_w_m %>% mutate(Bet_Manual = Giml_Manual / Alep_Manual)\n",
    "\n",
    "print(qsum1_b_w_m %>% \n",
    "          ggplot(aes(x = q, y = Bet_Manual)) + \n",
    "          geom_line( color = 'blue') + \n",
    "          geom_hline(yintercept = bet_w_m$Beta_Richness, color = 'green') +\n",
    "          geom_vline(xintercept = 0, color = 'green') + \n",
    "          geom_hline(yintercept = bet_w_m$Beta_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          labs(title = 'Weighted Beta Diversity Computed Using Fromulas') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\beta$ Diversity')))  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "qsum2_b_w_m <- inner_join(qsum2_g_w_m, qsum2_a_w_m, by = 'q')\n",
    "\n",
    "qsum2_b_w_m <- qsum2_b_w_m %>% mutate(Bet_Manual = Giml_Manual / Alep_Manual)\n",
    "\n",
    "print(qsum2_b_w_m %>% \n",
    "          ggplot(aes(x = q, y = Bet_Manual)) + \n",
    "          geom_line( color = 'blue') + \n",
    "          geom_hline(yintercept = bet_w_m$Beta_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          geom_hline(yintercept = bet_w_m$Beta_Greenberg, color = 'red') +\n",
    "          geom_vline(xintercept = 2, color = 'red') + \n",
    "          labs(title = 'Weighted Beta Diversity Computed Using Fromulas') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\beta$ Diversity')))  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "qsum_b_w_m <- rbind(qsum1_b_w_m, qsum2_b_w_m)\n",
    "\n",
    "print(qsum_b_w_m %>% \n",
    "          ggplot(aes(x = q, y = Bet_Manual)) + \n",
    "          geom_line( color = 'blue') +\n",
    "          geom_hline(yintercept = bet_w_m$Beta_Shannon, color = 'yellow') +\n",
    "          geom_vline(xintercept = 1, color = 'yellow') + \n",
    "          geom_hline(yintercept = bet_w_m$Beta_Greenberg, color = 'red') +\n",
    "          geom_vline(xintercept = 2, color = 'red') + \n",
    "          labs(title = 'Weighted Beta Diversity Computed Using Fromulas') +\n",
    "          xlab(TeX('$q$')) +\n",
    "          ylab(TeX('$\\\\beta$ Diversity'))) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# MacArthur’s Homogeneity Measure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Order 0 (Richness)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "alep_u_m$Alpha_Richness / giml_u_m$Gamma_Richness "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Order 1 (Shannon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "alep_u_m$Alpha_Shannon / giml_u_m$Gamma_Shannon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Order 2 (Greenberg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "alep_u_m$Alpha_Greenberg / giml_u_m$Gamma_Greenberg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Relative Homogeneity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Unweighted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "homog_u_m <- data.frame(list((1 / bet_u_m$Beta_Richness - 1 / N) / (1 - 1 / N),\n",
    "                             (1 / bet_u_m$Beta_Shannon - 1 / N) / (1 - 1 / N), \n",
    "                             (1 / bet_u_m$Beta_Greenberg - 1 / N) / (1 - 1 / N)))\n",
    "\n",
    "names(homog_u_m) <- c('Order 0 Homogeneity', \n",
    "                      'Order 1 Homogeneity', \n",
    "                      'Order 2 Homogeneity')\n",
    "\n",
    "homog_u_m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Weighted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "weights <- unique(giml$Weight)\n",
    "\n",
    "d_1_w_m <- exp(-sum(weights * log(weights))) \n",
    "\n",
    "(1 / bet_w_m$Beta_Shannon - 1 / d_1_w_m) / (1 - 1 / d_1_w_m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "# Turnover\n",
    "\n",
    "### Order 0 (Richness)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "(bet_u_m$Beta_Richness - 1) / (N - 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Order 1 (Shannon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "(bet_u_m$Beta_Shannon - 1) / (N - 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Order 2 (Greenberg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "(bet_u_m$Beta_Greenberg - 1) / (N - 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Sørensen-Dice and Jaccard Indices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Sørensen-Dice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "sorensen <- NULL\n",
    "\n",
    "for (region in unique(compact_df$Region_New)) {\n",
    "    \n",
    "    df <- compact_df %>% filter(Region_New == region)\n",
    "    \n",
    "    country_list <- df %>% select(Country_Code) %>% unique() %>% unlist(use.names = FALSE)\n",
    "    \n",
    "    comb <- combn(country_list, 2, simplify = FALSE)\n",
    "    \n",
    "    for (i in 1:length(comb)) {\n",
    "        \n",
    "        n <- comb[[i]][1]\n",
    "        \n",
    "        m <- comb[[i]][2]\n",
    "           \n",
    "        language_list_1 <- df %>% \n",
    "            filter(Country_Code == n) %>% \n",
    "            select(Glottocode) \n",
    "        \n",
    "        language_list_2 <- df %>% \n",
    "            filter(Country_Code == m) %>% \n",
    "            select(Glottocode)\n",
    "        \n",
    "        family_list_1 <- df %>% \n",
    "            filter(Country_Code == n) %>% \n",
    "            select(Glotto_Family) %>%\n",
    "            unique()\n",
    "        \n",
    "        family_list_2 <- df %>% \n",
    "            filter(Country_Code == m) %>% \n",
    "            select(Glotto_Family) %>%\n",
    "            unique()\n",
    "       \n",
    "        suppressMessages(inner_lang <- nrow(inner_join(language_list_1, language_list_2)))\n",
    "        \n",
    "        suppressMessages(outer_lang <- nrow(full_join(language_list_1, language_list_2)))\n",
    "        \n",
    "        suppressMessages(diff_lang_1 <- nrow(anti_join(language_list_1, language_list_2)))\n",
    "        \n",
    "        suppressMessages(diff_lang_2 <- nrow(anti_join(language_list_2, language_list_1)))\n",
    "        \n",
    "        suppressMessages(inner_fam <- nrow(inner_join(family_list_1, family_list_2)))\n",
    "        \n",
    "        suppressMessages(outer_fam <- nrow(full_join(family_list_1, family_list_2)))\n",
    "        \n",
    "        suppressMessages(diff_fam_1 <- nrow(anti_join(family_list_1, family_list_2)))\n",
    "        \n",
    "        suppressMessages(diff_fam_2 <- nrow(anti_join(family_list_2, family_list_1)))\n",
    "        \n",
    "        sorensen <- rbind(sorensen, data.frame(region, \n",
    "                                               n, \n",
    "                                               m, \n",
    "                                               inner_lang, \n",
    "                                               outer_lang, \n",
    "                                               diff_lang_1, \n",
    "                                               diff_lang_2, \n",
    "                                               2 * inner_lang / (diff_lang_1 + diff_lang_2 + 2 * inner_lang), \n",
    "                                               inner_fam, \n",
    "                                               outer_fam, \n",
    "                                               diff_fam_1, \n",
    "                                               diff_fam_2, \n",
    "                                               2 * inner_fam / (diff_fam_1 + diff_fam_2 + 2 * inner_fam)))\n",
    "    }\n",
    "}\n",
    "\n",
    "\n",
    "names(sorensen) <- c('Region_New', \n",
    "                     'Country_Code_1', \n",
    "                     'Country_Code_2', \n",
    "                     'Intersection_Language', \n",
    "                     'Union_Language', \n",
    "                     'Diff_Language_1', \n",
    "                     'Diff_Language_2', \n",
    "                     'Sorensen_Language',\n",
    "                     'Intersection_Family',\n",
    "                     'Union_Family', \n",
    "                     'Diff_Family_1', \n",
    "                     'Diff_Family_2', \n",
    "                     'Sorensen_Family')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "sorensen %>% \n",
    "    ggplot(aes(x = Sorensen_Language, fill = Region_New)) + \n",
    "    geom_histogram(binwidth = 0.1, show.legend = FALSE) + \n",
    "    scale_fill_manual(values = cbPalette) + \n",
    "    labs(title = 'Language-based Histogram for Sørensen-Dice Indices Based on Region') + \n",
    "    ylab('Count') +\n",
    "    xlab('Sørensen-Dice Index') + \n",
    "    guides(color = guide_legend(label.position = 'top')) + \n",
    "    facet_grid(rows = vars(Region_New)) + \n",
    "    theme(strip.text.y = element_text(size = 7.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "sorensen %>% \n",
    "    ggplot(aes(x = Sorensen_Language, color = Region_New)) + \n",
    "    geom_density(aes(linetype = Region_New)) + \n",
    "    scale_color_manual(values = cbPalette) + \n",
    "    labs(title = 'Language-based KDE Curves for Sørensen-Dice Indices Based on Region') + \n",
    "    ylab('Density') +\n",
    "    xlab('Sørensen-Dice Index') +     \n",
    "    guides(color = guide_legend(label.position = 'top', title = 'Region'), \n",
    "           linetype = guide_legend(label.position = 'top', title = 'Region'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "sorensen %>% \n",
    "    ggplot(aes(x = Sorensen_Family, fill = Region_New)) + \n",
    "    geom_histogram(binwidth = 0.2, show.legend = FALSE) + \n",
    "    scale_fill_manual(values = cbPalette) + \n",
    "    labs(title = 'Family-based Histogram for Sørensen-Dice Indices Based on Region') + \n",
    "    ylab('Count') +\n",
    "    xlab('Sørensen-Dice Index') + \n",
    "    guides(color = guide_legend(label.position = 'top')) + \n",
    "    facet_grid(rows = vars(Region_New)) + \n",
    "    theme(strip.text.y = element_text(size = 7.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "sorensen %>% \n",
    "    ggplot(aes(x = Sorensen_Family, color = Region_New)) + \n",
    "    geom_density(aes(linetype = Region_New)) + \n",
    "    scale_color_manual(values = cbPalette) + \n",
    "    labs(title = 'Family-based KDE Curves for Sørensen-Dice Indices Based on Region') +\n",
    "    ylab('Density') +\n",
    "    xlab('Sørensen-Dice Index') + \n",
    "    guides(color = guide_legend(label.position = 'top', title = 'Region'), \n",
    "           linetype = guide_legend(label.position = 'top', title = 'Region'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "sorensen %>% \n",
    "    arrange(desc(Sorensen_Language)) %>%\n",
    "    select(Region_New, Country_Code_1, Country_Code_2, Sorensen_Language) %>%\n",
    "    head(10) %>%\n",
    "    remove_rownames()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "sorensen %>% \n",
    "    arrange(desc(Sorensen_Language)) %>%\n",
    "    select(Region_New, Country_Code_1, Country_Code_2, Sorensen_Language) %>%\n",
    "    tail(10) %>%\n",
    "    remove_rownames()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "sorensen %>% \n",
    "    arrange(desc(Sorensen_Family)) %>%\n",
    "    select(Region_New, Country_Code_1, Country_Code_2, Sorensen_Family) %>%\n",
    "    head(10) %>%\n",
    "    remove_rownames()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "sorensen %>% \n",
    "    arrange(desc(Sorensen_Family)) %>%\n",
    "    select(Region_New, Country_Code_1, Country_Code_2, Sorensen_Family) %>%\n",
    "    tail(10) %>%\n",
    "    remove_rownames()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "suppressMessages(sorensen %>% group_by(Region_New) %>%\n",
    "    summarize(Mean_Sorensen_Language = mean(Sorensen_Language), \n",
    "              Mean_Sorensen_Family = mean(Sorensen_Family)) %>%\n",
    "    ggplot(aes(x = Mean_Sorensen_Language, y = reorder(Region_New, Mean_Sorensen_Language))) + \n",
    "    geom_segment(aes(x = 0, \n",
    "                     xend = Mean_Sorensen_Language, \n",
    "                     y = reorder(Region_New, Mean_Sorensen_Language), \n",
    "                     yend = reorder(Region_New, Mean_Sorensen_Language)), \n",
    "                     color = 'skyblue') +\n",
    "    geom_point(show.legend = FALSE, color = 'blue') +\n",
    "    labs(title = 'Language-based Lolliplot of Sørensen-Dice Indices in Each Region') +\n",
    "    ylab('Region') +\n",
    "    xlab('Mean Sørensen-Dice Index'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "suppressMessages(sorensen %>% group_by(Region_New) %>%\n",
    "    summarize(Mean_Sorensen_Language = mean(Sorensen_Language), \n",
    "              Mean_Sorensen_Family = mean(Sorensen_Family)) %>%\n",
    "    ggplot(aes(x = Mean_Sorensen_Family, y = reorder(Region_New, Mean_Sorensen_Family))) + \n",
    "    geom_segment(aes(x = 0, \n",
    "                     xend = Mean_Sorensen_Family, \n",
    "                     y = reorder(Region_New, Mean_Sorensen_Family), \n",
    "                     yend = reorder(Region_New, Mean_Sorensen_Family)), \n",
    "                     color = 'skyblue') +    \n",
    "    geom_point(show.legend = FALSE, color = 'blue') +\n",
    "    labs(title = 'Family-based Lolliplot of Sørensen-Dice Indices in Each Region') +\n",
    "    ylab('Region') +\n",
    "    xlab('Mean Sørensen-Dice Index'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "sorensen_2 <- NULL\n",
    "\n",
    "comb <- combn(unique(compact_df$Region_New), 2, simplify = FALSE)\n",
    "\n",
    "for (i in 1:length(comb)) {\n",
    "        \n",
    "    n <- comb[[i]][1]\n",
    "        \n",
    "    m <- comb[[i]][2]\n",
    "           \n",
    "    language_list_1 <- main_tbl %>% \n",
    "        filter(Region_New == n) %>% \n",
    "        select(ISO_639) \n",
    "        \n",
    "    language_list_2 <- main_tbl %>% \n",
    "        filter(Region_New == m) %>% \n",
    "        select(ISO_639)\n",
    "    \n",
    "    family_list_1 <- df %>% \n",
    "            filter(Country_Code == n) %>% \n",
    "            select(Glotto_Family) %>%\n",
    "            unique()\n",
    "        \n",
    "    family_list_2 <- df %>% \n",
    "            filter(Country_Code == m) %>% \n",
    "            select(Glotto_Family) %>%\n",
    "            unique()\n",
    "        \n",
    "    suppressMessages(inner <- nrow(inner_join(language_list_1, language_list_2)))\n",
    "        \n",
    "    suppressMessages(outer <- nrow(full_join(language_list_1, language_list_2)))\n",
    "        \n",
    "    suppressMessages(diff_1 <- nrow(anti_join(language_list_1, language_list_2)))\n",
    "        \n",
    "    suppressMessages(diff_2 <- nrow(anti_join(language_list_2, language_list_1)))\n",
    "        \n",
    "    sorensen_2 <- rbind(sorensen_2, data.frame(str_c(n, m, sep = '--'), \n",
    "                                               inner, \n",
    "                                               outer, \n",
    "                                               diff_1, \n",
    "                                               diff_2, \n",
    "                                               2 * inner / (diff_1 + diff_2 + 2 * inner)))\n",
    "\n",
    "}\n",
    "\n",
    "\n",
    "names(sorensen_2) <- c('Regions', \n",
    "                      'Intersection', \n",
    "                      'Union', \n",
    "                      'Diff_1', \n",
    "                      'Diff_2', \n",
    "                      'Sorensen')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "sorensen_2 %>% \n",
    "    arrange(desc(Sorensen)) %>%\n",
    "    select(Regions, Sorensen)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "sorensen_2 %>% \n",
    "    ggplot(aes(x = Sorensen, y = reorder(Regions, Sorensen))) + \n",
    "    geom_segment(aes(x = 0, \n",
    "                     xend = Sorensen, \n",
    "                     y = reorder(Regions, Sorensen), \n",
    "                     yend = reorder(Regions, Sorensen)), \n",
    "                     color = 'skyblue') +\n",
    "    geom_point(show.legend = FALSE, color = 'blue') + \n",
    "    labs(title = 'Lolliplot for Sørensen-Dice Indices between Regions') +\n",
    "    ylab('Pair of Regions') + \n",
    "    xlab('Sørensen-Dice Index')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Jaccard"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "jaccard <- NULL\n",
    "\n",
    "for (region in unique(compact_df$Region_New)) {\n",
    "    \n",
    "    df <- compact_df %>% filter(Region_New == region)\n",
    "    \n",
    "    country_list <- df %>% select(Country_Code) %>% unique() %>% unlist(use.names = FALSE)\n",
    "    \n",
    "    comb <- combn(country_list, 2, simplify = FALSE)\n",
    "    \n",
    "    for (i in 1:length(comb)) {\n",
    "        \n",
    "        n <- comb[[i]][1]\n",
    "        \n",
    "        m <- comb[[i]][2]\n",
    "           \n",
    "        language_list_1 <- df %>% \n",
    "            filter(Country_Code == n) %>% \n",
    "            select(Glottocode) \n",
    "        \n",
    "        language_list_2 <- df %>% \n",
    "            filter(Country_Code == m) %>% \n",
    "            select(Glottocode)\n",
    "        \n",
    "        family_list_1 <- df %>% \n",
    "            filter(Country_Code == n) %>% \n",
    "            select(Glotto_Family) %>%\n",
    "            unique()\n",
    "        \n",
    "        family_list_2 <- df %>% \n",
    "            filter(Country_Code == m) %>% \n",
    "            select(Glotto_Family) %>%\n",
    "            unique()\n",
    "       \n",
    "        suppressMessages(inner_lang <- nrow(inner_join(language_list_1, language_list_2)))\n",
    "        \n",
    "        suppressMessages(outer_lang <- nrow(full_join(language_list_1, language_list_2)))\n",
    "        \n",
    "        suppressMessages(diff_lang_1 <- nrow(anti_join(language_list_1, language_list_2)))\n",
    "        \n",
    "        suppressMessages(diff_lang_2 <- nrow(anti_join(language_list_2, language_list_1)))\n",
    "        \n",
    "        suppressMessages(inner_fam <- nrow(inner_join(family_list_1, family_list_2)))\n",
    "        \n",
    "        suppressMessages(outer_fam <- nrow(full_join(family_list_1, family_list_2)))\n",
    "        \n",
    "        suppressMessages(diff_fam_1 <- nrow(anti_join(family_list_1, family_list_2)))\n",
    "        \n",
    "        suppressMessages(diff_fam_2 <- nrow(anti_join(family_list_2, family_list_1)))\n",
    "        \n",
    "        jaccard <- rbind(jaccard, data.frame(region, \n",
    "                                           n, \n",
    "                                           m, \n",
    "                                           inner_lang, \n",
    "                                           outer_lang, \n",
    "                                           diff_lang_1, \n",
    "                                           diff_lang_2, \n",
    "                                           inner_lang / outer_lang, \n",
    "                                           inner_fam, \n",
    "                                           outer_fam, \n",
    "                                           diff_fam_1, \n",
    "                                           diff_fam_2, \n",
    "                                           inner_fam / outer_fam))\n",
    "    }\n",
    "}\n",
    "\n",
    "\n",
    "names(jaccard) <- c('Region_New', \n",
    "                   'Country_Code_1', \n",
    "                   'Country_Code_2', \n",
    "                   'Intersection_Language', \n",
    "                   'Union_Language', \n",
    "                   'Diff_Language_1', \n",
    "                   'Diff_Language_2', \n",
    "                   'Jaccard_Language',\n",
    "                   'Intersection_Family',\n",
    "                   'Union_Family', \n",
    "                   'Diff_Family_1', \n",
    "                   'Diff_Family_2', \n",
    "                   'Jaccard_Family')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "jaccard %>% \n",
    "    ggplot(aes(x = Jaccard_Language, fill = Region_New)) + \n",
    "    geom_histogram(binwidth = 0.1, show.legend = FALSE) + \n",
    "    scale_fill_manual(values = cbPalette) + \n",
    "    labs(title = 'Histogram for Language Jaccard Indices Based on Region') + \n",
    "    ylab('Count') +\n",
    "    xlab('Jaccard Index') +\n",
    "    guides(color = guide_legend(label.position = 'top')) + \n",
    "    facet_grid(rows = vars(Region_New)) + \n",
    "    theme(strip.text.y = element_text(size = 7.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "jaccard %>% \n",
    "    ggplot(aes(x = Jaccard_Language, color = Region_New)) + \n",
    "    geom_density(aes(linetype = Region_New)) + \n",
    "    scale_color_manual(values = cbPalette) + \n",
    "    labs(title = 'Language-based KDE Curves for Jaccard Indices Based on Region') + \n",
    "    ylab('Density') +\n",
    "    xlab('Jaccard Index') +\n",
    "    guides(color = guide_legend(label.position = 'top', title = 'Region'), \n",
    "           linetype = guide_legend(label.position = 'top', title = 'Region'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "jaccard %>% \n",
    "    ggplot(aes(x = Jaccard_Family, fill = Region_New)) + \n",
    "    geom_histogram(binwidth = 0.2, show.legend = FALSE) + \n",
    "    scale_fill_manual(values = cbPalette) + \n",
    "    labs(title = 'Family-based Histogram for Jaccard Indices Based on Region') + \n",
    "    ylab('Count') +\n",
    "    xlab('Jaccard Index') +\n",
    "    guides(color = guide_legend(label.position = 'top')) + \n",
    "    facet_grid(rows = vars(Region_New)) + \n",
    "    theme(strip.text.y = element_text(size = 7.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "jaccard %>% \n",
    "    ggplot(aes(x = Jaccard_Family, color = Region_New)) + \n",
    "    geom_density(aes(linetype = Region_New)) + \n",
    "    scale_color_manual(values = cbPalette) + \n",
    "    labs(title = 'Famiy-based KDE Curves for Jaccard Indices Based on Region') + \n",
    "    ylab('Density') +\n",
    "    xlab('Jaccard Index') +\n",
    "    guides(color = guide_legend(label.position = 'top', title = 'Region'), \n",
    "           linetype = guide_legend(label.position = 'top', title = 'Region'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "jaccard %>% \n",
    "    arrange(desc(Jaccard_Language)) %>%\n",
    "    select(Region_New, Country_Code_1, Country_Code_2, Jaccard_Language) %>%\n",
    "    head(10) %>%\n",
    "    remove_rownames()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "jaccard %>% \n",
    "    arrange(desc(Jaccard_Language)) %>%\n",
    "    select(Region_New, Country_Code_1, Country_Code_2, Jaccard_Language, Jaccard_Family) %>%\n",
    "    tail(10)  %>%\n",
    "    remove_rownames()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "jaccard %>% \n",
    "    arrange(desc(Jaccard_Family)) %>%\n",
    "    select(Region_New, Country_Code_1, Country_Code_2, Jaccard_Family) %>%\n",
    "    head(10) %>%\n",
    "    remove_rownames()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "jaccard %>% \n",
    "    arrange(desc(Jaccard_Family)) %>%\n",
    "    select(Region_New, Country_Code_1, Country_Code_2, Jaccard_Family) %>%\n",
    "    tail(10) %>%\n",
    "    remove_rownames()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "suppressMessages(jaccard %>% group_by(Region_New) %>%\n",
    "    summarize(Mean_Jaccard_Language = mean(Jaccard_Language), \n",
    "              Mean_Jaccard_Family = mean(Jaccard_Family)) %>%\n",
    "    ggplot(aes(x = Mean_Jaccard_Language, y = reorder(Region_New, Mean_Jaccard_Language))) + \n",
    "    geom_segment(aes(x = 0, \n",
    "                     xend = Mean_Jaccard_Language, \n",
    "                     y = reorder(Region_New, Mean_Jaccard_Language), \n",
    "                     yend = reorder(Region_New, Mean_Jaccard_Language)), \n",
    "                     color = 'skyblue') +\n",
    "    geom_point(show.legend = FALSE, color = 'blue') +\n",
    "    labs(title = 'Language-based Lolliplot of Jaccard Indices in Each Region') +\n",
    "    ylab('Region') +\n",
    "    xlab('Mean Jaccard Index'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "suppressMessages(jaccard %>% group_by(Region_New) %>%\n",
    "    summarize(Mean_Jaccard_Language = mean(Jaccard_Language), \n",
    "              Mean_Jaccard_Family = mean(Jaccard_Family)) %>%\n",
    "    ggplot(aes(x = Mean_Jaccard_Family, y = reorder(Region_New, Mean_Jaccard_Family))) + \n",
    "    geom_segment(aes(x = 0, \n",
    "                     xend = Mean_Jaccard_Family, \n",
    "                     y = reorder(Region_New, Mean_Jaccard_Family), \n",
    "                     yend = reorder(Region_New, Mean_Jaccard_Family)), \n",
    "                     color = 'skyblue') +    \n",
    "    geom_point(show.legend = FALSE, color = 'blue') +\n",
    "    labs(title = 'Family-based Lolliplot of Jaccard Indices in Each Region') +\n",
    "    ylab('Region') +\n",
    "    xlab('Mean Jaccard Index'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "jaccard_2 <- NULL\n",
    "\n",
    "comb <- combn(unique(compact_df$Region_New), 2, simplify = FALSE)\n",
    "\n",
    "for (i in 1:length(comb)) {\n",
    "        \n",
    "    n <- comb[[i]][1]\n",
    "        \n",
    "    m <- comb[[i]][2]\n",
    "           \n",
    "    language_list_1 <- main_tbl %>% \n",
    "        filter(Region_New == n) %>% \n",
    "        select(ISO_639) \n",
    "        \n",
    "    language_list_2 <- main_tbl %>% \n",
    "        filter(Region_New == m) %>% \n",
    "        select(ISO_639)\n",
    "    \n",
    "    family_list_1 <- df %>% \n",
    "            filter(Country_Code == n) %>% \n",
    "            select(Glotto_Family) %>%\n",
    "            unique()\n",
    "        \n",
    "    family_list_2 <- df %>% \n",
    "            filter(Country_Code == m) %>% \n",
    "            select(Glotto_Family) %>%\n",
    "            unique()\n",
    "        \n",
    "    suppressMessages(inner <- nrow(inner_join(language_list_1, language_list_2)))\n",
    "        \n",
    "    suppressMessages(outer <- nrow(full_join(language_list_1, language_list_2)))\n",
    "        \n",
    "    suppressMessages(diff_1 <- nrow(anti_join(language_list_1, language_list_2)))\n",
    "        \n",
    "    suppressMessages(diff_2 <- nrow(anti_join(language_list_2, language_list_1)))\n",
    "        \n",
    "    jaccard_2 <- rbind(jaccard_2, data.frame(str_c(n, m, sep = '--'), \n",
    "                                               inner, \n",
    "                                               outer, \n",
    "                                               diff_1, \n",
    "                                               diff_2, \n",
    "                                               inner / outer))\n",
    "\n",
    "}\n",
    "\n",
    "\n",
    "names(jaccard_2) <- c('Regions', \n",
    "                     'Intersection', \n",
    "                     'Union', \n",
    "                     'Diff_1', \n",
    "                     'Diff_2', \n",
    "                     'Jaccard')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "jaccard_2 %>% \n",
    "    arrange(desc(Jaccard)) %>%\n",
    "    select(Regions, Jaccard)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "jaccard_2 %>% \n",
    "    ggplot(aes(x = Jaccard, y = reorder(Regions, Jaccard))) + \n",
    "    geom_segment(aes(x = 0, \n",
    "                     xend = Jaccard, \n",
    "                     y = reorder(Regions, Jaccard), \n",
    "                     yend = reorder(Regions, Jaccard)), \n",
    "                     color = 'skyblue') +\n",
    "    geom_point(show.legend = FALSE, color = 'blue') + \n",
    "    labs(title = 'Lolliplot for Jaccard Indices between Regions') +\n",
    "    ylab('Pair of Regions') +\n",
    "    xlab('Jaccard Index')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Thank You! "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}