{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Global Spatial Autocorrelation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import mapclassify as mc\n",
    "import libpysal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A well-used functionality in PySAL is the use of PySAL to conduct exploratory spatial data analysis. This notebook will provide an overview of ways to conduct exploratory spatial analysis in Python. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import geopandas as gpd\n",
    "shp_link = \"data/texas.shp\"\n",
    "tx = gpd.read_file(shp_link)\n",
    "hr10 = mc.Quantiles(tx.HR90, k=10)\n",
    "f, ax = plt.subplots(1, figsize=(9, 9))\n",
    "tx.assign(cl=hr10.yb).plot(column='cl', categorical=True, \\\n",
    "        k=10, cmap='OrRd', linewidth=0.1, ax=ax, \\\n",
    "        edgecolor='white', legend=True)\n",
    "ax.set_axis_off()\n",
    "plt.title(\"HR90 Deciles\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spatial Autocorrelation\n",
    "\n",
    "Visual inspection of the map pattern for HR90 deciles allows us to search for spatial structure. If the spatial distribution of the rates was random, then we should not see any clustering of similar values on the map. However, our visual system is drawn to the darker clusters in the south west as well as the east, and a concentration of the lighter hues (lower homicide rates) moving north to the pan handle.\n",
    "\n",
    "Our brains are very powerful pattern recognition machines. However, sometimes they can be too powerful and lead us to detect false positives, or patterns where there are no statistical patterns. This is a particular concern when dealing with visualization of irregular polygons of differning sizes and shapes.\n",
    "\n",
    "The concept of *spatial autocorrelation* relates to the combination of two types of similarity: spatial similarity and attribute similarity. Although there are many different measures of spatial autocorrelation, they all combine these two types of simmilarity into a summary measure.\n",
    "\n",
    "Let's use PySAL to generate these two types of similarity measures.\n",
    "\n",
    "### Spatial Similarity\n",
    "\n",
    "We have already encountered spatial weights in a previous notebook. In spatial autocorrelation analysis, the spatial weights are used to formalize the notion of spatial similarity. As we have seen there are many ways to define spatial weights, here we will use queen contiguity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = libpysal.weights.Queen.from_shapefile(\"data/texas.shp\")\n",
    "W.transform = 'r'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Attribute Similarity\n",
    "\n",
    "So the spatial weight between counties $i$ and $j$ indicates if the two counties are neighbors (i.e., geographically similar). What we also need is a measure of attribute similarity to pair up with this concept of spatial similarity.\n",
    "The **spatial lag** is a derived variable that accomplishes this for us. For county $i$ the spatial lag is defined as:\n",
    "$$HR90Lag_i = \\sum_j w_{i,j} HR90_j$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HR90Lag = libpysal.weights.lag_spatial(W, data.HR90)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HR90LagQ10 = mc.Quantiles(HR90Lag, k=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(1, figsize=(9, 9))\n",
    "tx.assign(cl=HR90LagQ10.yb).plot(column='cl', categorical=True, \\\n",
    "        k=10, cmap='OrRd', linewidth=0.1, ax=ax, \\\n",
    "        edgecolor='white', legend=True)\n",
    "ax.set_axis_off()\n",
    "plt.title(\"HR90 Spatial Lag Deciles\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The decile map for the spatial lag tends to enhance the impression of value similarity in space. However, we still have the challenge of visually associating the value of the homicide rate in a county with the value of the spatial lag of rates for the county. The latter is a weighted average of homicide rates in the focal county's neighborhood.\n",
    "\n",
    "To complement the geovisualization of these associations we can turn to formal statistical measures of spatial autocorrelation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HR90 = data.HR90\n",
    "b,a = np.polyfit(HR90, HR90Lag, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(1, figsize=(9, 9))\n",
    "\n",
    "plt.plot(HR90, HR90Lag, '.', color='firebrick')\n",
    "\n",
    " # dashed vert at mean of the last year's PCI\n",
    "plt.vlines(HR90.mean(), HR90Lag.min(), HR90Lag.max(), linestyle='--')\n",
    " # dashed horizontal at mean of lagged PCI\n",
    "plt.hlines(HR90Lag.mean(), HR90.min(), HR90.max(), linestyle='--')\n",
    "\n",
    "# red line of best fit using global I as slope\n",
    "plt.plot(HR90, a + b*HR90, 'r')\n",
    "plt.title('Moran Scatterplot')\n",
    "plt.ylabel('Spatial Lag of HR90')\n",
    "plt.xlabel('HR90')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Global Spatial Autocorrelation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In PySAL, commonly-used analysis methods are very easy to access. For example, if we were interested in examining the spatial dependence in `HR90` we could quickly compute a Moran's $I$ statistic:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import esda\n",
    "I_HR90 = esda.Moran(data.HR90.values, W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I_HR90.I, I_HR90.p_sim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, the $I$ statistic is $0.859$ for this data, and has a very small $p$ value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b # note I is same as the slope of the line in the scatterplot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualize the distribution of simulated $I$ statistics using the stored collection of simulated statistics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I_HR90.sim[0:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simple way to visualize this distribution is to make a KDEplot (like we've done before), and add a rug showing all of the simulated points, and a vertical line denoting the observed value of the statistic:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.kdeplot(I_HR90.sim, shade=True)\n",
    "plt.vlines(I_HR90.sim, 0, 0.5)\n",
    "plt.vlines(I_HR90.I, 0, 10, 'r')\n",
    "plt.xlim([-0.15, 0.15])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead, if our $I$ statistic were close to our expected value, `I_HR90.EI`, our plot might look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.kdeplot(I_HR90.sim, shade=True)\n",
    "plt.vlines(I_HR90.sim, 0, 1)\n",
    "plt.vlines(I_HR90.EI+.01, 0, 10, 'r')\n",
    "plt.xlim([-0.15, 0.15])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result of applying Moran's I is that we conclude the map pattern is not spatially random, but instead there is a signficant spatial association in homicide rates in Texas counties in 1990.\n",
    "\n",
    "This result applies to the map as a whole, and is sometimes referred to as \"global spatial autocorrelation\". Later we turn to a local analysis where the attention shifts to detection of hot spots, cold spots and spatial outliers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "1. Repeat the global analysis for the years 1960, 70, 80 and compare the results to what we found in 1990.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(12345)\n",
    "I_pv = [esda.Moran(data[var].values, W).p_sim for var in ['HR60', 'HR70', 'HR80', 'HR90']]\n",
    "I_pv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}