{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Calculating Return on Investment for Hot Spots Policing\n", "\n", "By [Andrew Wheeler, PhD](mailto:apwheele@gmail.com) \n", "Website: [andrewpwheeler.com](https://andrewpwheeler.com/)\n", "\n", "So we want to go from the Cohen's D measures in the [Braga hot spot meta-analysis](https://www.tandfonline.com/doi/abs/10.1080/07418825.2012.673632), to a measure more directly related to how much crime is reduced. What I want to do in the end is to recreate the pre and post counts of crime, so I can estimate the direct crime reduction effect of the number of crime reduced.\n", "\n", "First, the Cohen's D measure for these interventions is calculated [as below](https://www.ncbi.nlm.nih.gov/pubmed/7870860):\n", "\n", "$$ D = \\log(OR) \\cdot \\sqrt{ 3 }/\\pi $$\n", "\n", "Where $OR$ is the odds ratio. The odds ratio is calculated from the pre-post counts of crime table. \n", "\n", "| | Pre | Post |\n", "|---------|-------|-------|\n", "| Treated | $t_0$ | $t_1$ |\n", "| Control | $c_0$ | $c_1$ |\n", "\n", "The odds ratio is calculated as:\n", "\n", "$$ OR = \\frac{t_1/t_0}{c_1/c_0} $$\n", "\n", "So we have a few unknowns here to back calculate a return on investment for a particular crime intervention. First, the baseline of crime makes a big difference. Using the typical interpretation of the $OR$ as a percent reduction, say we have a 10% decline (i.e. an $OR = 0.9$). If we start with a baseline of 100 crimes, that is a net reduction of 10 crimes ($100 - OR \\cdot 100 = 10$). But if we start with a baseline of 10 crimes, it is only a reduction of 1 crime ($10 - OR \\cdot 10 = 1$). So those are big differences! The ROI is tied directly to the baseline crime counts. \n", "\n", "Second, this is ignoring the control areas. For simplicity in this analysis, first I am going to assume the control areas do not change, so $c_1 = c_0$. This then reduces the odds ratio formula to simply $OR = t_1/t_0$. Second, I assume that $t_0 = c_0 = c_1$. This is basically saying that the control areas had the same numbers of crime to the treated area at baseline. \n", "\n", "So the formula then to go from the Cohen's D estimates in the Braga meta-analysis to an estimate of $t_1$ and $t_0$ would be:\n", "\n", "$$ \\exp( \\frac{D \\cdot \\pi}{\\sqrt{3}} ) = OR = t_1/t_0 $$ \n", "\n", "So subsquently if I want to calculate $t_1 - t_0 = \\Delta t$, if we fix $t_0$ to some arbitrary value, we then have:\n", "\n", "$$ OR \\cdot t_0 = t_1 $$\n", "\n", "So then:\n", "\n", "$$ \\Delta t = OR \\cdot t_0 - t_0 $$\n", "\n", "What I do below is based on the overall Cohen's D effect size, I plot $\\Delta t$ given a fixed value of $t_0$. Note that the Braga meta-analysis gives the estimates where D is positive (but signifies a reduction in crime). So I actually take the negative of the Cohen's D value. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "meta analysis D of -0.184 to odds ratio of 0.716\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAAF3CAYAAAC123K4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdf3yO9f7A8dfH1BCrUciPQqHG6Z6Z1jpKX6U6HU6dnPwWIiolpV+qo1Q6dVKplIhCJfRDoU6SGMuMmd35UaMvOiiULSYztn2+fxhfp+PH7u2+ruu+3ns/Hw8P22z3/X6dz7fvZ/d93fd1GWstSimllPKfSl4PoJRSSqmy0U1cKaWU8indxJVSSimf0k1cKaWU8indxJVSSimf0k1cKaWU8qnKXg8QqjPPPNM2atTI6zGUUkopV6xcufIXa+1Zx/o3323ijRo1IiMjw+sxwmb16tX84Q9/8HqMsJPYpU3+IbFLYhPI7Ap3kzHmh+P9mz6d7rFdu3Z5PYIjJHZpk39I7JLYBDK73GzSTVwppZTyKd3EPRYIBLwewRESu7TJPyR2SWwCmV1uNukm7rG8vDyvR3CExC5t8g+JXRKbQGaXm026iXts48aNXo/gCIld2uQfErskNoHMLjebdBNXSimlfEo3cY9Jfc+7xC5t8g+JXRKbQGaXm026iXusZs2aXo/gCIld2uQfErskNoHMLjebdBP3WGZmptcjOEJilzb5h8QuiU0gs8vNJt3ElVJKKZ/STdxjsbGxXo/gCIld2uQfErskNoHMLjebjLXWtTsLh8TERCvp3OlKKaXkKCoqorCwkOjo6LDdpjFmpbU28Vj/po/EPZaSkuL1CI6Q2KVN/iGxS2ITyOr68ccf6dChAzfddJNr96mbuMf89kxIaUns0ib/kNglsQnkdM2dO5eLLrqI9PR0LrzwQtfuVzdxjxljvB7BERK7tMk/JHZJbAL/dxUUFDB06FA6depEgwYNWLlyJdddd51r96/HxJVSSqkyWL9+Pd26dWPVqlUMGTKEZ599lipVqoT9fvSYeAQLBoNej+AIiV3a5B8SuyQ2gX+7pk6dSkJCAj/88AOffPIJL7300pEN3M0m3cQ9lpub6/UIjpDYpU3+IbFLYhP4rysvL4/evXvTp08fEhMTCQaD/OUvf/mP73GzSTdxpZRSqhRWrlxJQkIC06ZNY+TIkSxYsIAGDRp4OpNu4h5LSEjwegRHSOzSJv+Q2CWxCfzRVVxczAsvvEBycjIFBQUsWrSIESNGEBUVdczvd7NJN3GP5eTkeD2CIyR2aZN/SOyS2ASR37Vz5046duzIsGHD+POf/0xWVhaXXXbZCX/GzSbdxD22efNmr0dwhMQubfIPiV0SmyCyuxYsWEAgEOCrr77i1Vdf5aOPPirVFcrcbNJNXCmllDrKwYMHefjhh+nQoQOxsbEsX76cO+64IyLf017Z6wEquiZNmng9giMkdmmTf0jsktgEkde1efNmunfvzrJlyxgwYABjxozhtNNOC+k23GzSTdxjNWrU8HoER0js0ib/kNglsQkiq+v999/n1ltvxVrL9OnT6dq1a5lux80mfTrdY3490cHJSOzSJv+Q2CWxCSKja9++fQwcOJAuXbpwwQUXkJWVVeYNHPRkL0oppZQr1qxZQ5s2bZg4cSIPPfQQS5YsoXHjxl6PVWq6iXusVq1aXo/gCIld2uQfErskNoF3XdZaXn/9ddq0aUNOTg7z5s3jH//4B6ecckq5b9vNJr0AiseKi4upVEne71ISu7TJPyR2SWwCb7pyc3MZMGAAH330Eddeey1Tpkyhdu3aYbv9cDfpBVAi2OLFi70ewRESu7TJPyR2SWwC97tSU1MJBALMnj2b0aNH8+mnn4Z1Awd3m3QTV0opJV5RURFPPfUU7dq149RTT2Xp0qUMGzbM989u6FvMPFa5sswlkNilTf4hsUtiE7jTtW3bNnr16sWiRYvo0aMH48aNIyYmxrH7c3Ot9Ji4UkopsebOnUvfvn3Jz8/n1VdfpU+fPhF55rUT0WPiESwzM9PrERwhsUub/ENil8QmcK6roKCAoUOH0qlTJxo2bEhmZiZ9+/Z1ZQN3c61kPj/jI3v27PF6BEdI7NIm/5DYJbEJnOlav3493bp1Y9WqVQwZMoR//vOfREdHh/1+jsfNtdJNXCmllAjWWqZOncrgwYOpUqUKs2fPplOnTl6P5Sg9Ju6xvXv3Ur16da/HCDuJXdrkHxK7JDZB+Lry8vK4/fbbeffdd2nXrh3vvvsu9evXD8OEoQv3Wukx8Qi2Y8cOr0dwhMQubfIPiV0SmyA8XRkZGbRq1Yr33nuPJ554ggULFni2gYO7a6WbuMe2bNni9QiOkNilTf4hsUtiE5Svq7i4mOeff55LL72UAwcOkJKSwt///neioqLCOGHo3FwrPSaulFLKd3bu3EmfPn34/PPP+etf/8rEiROpWbOm12O5Th+Je6xp06Zej+AIiV3a5B8SuyQ2Qdm6vvzySwKBAAsXLuS1117jww8/jKgN3M210k3cY26+7cFNEru0yT8kdklsgtC6Dh48yPDhw7n66quJjY1lxYoV3H777RF38hY310o3cY+tWbPG6xEcIbFLm/xDYpfEJih91+bNm7n88st55plnGDBgABkZGfzhD39weLqycXOt9Ji4UkqpiPb+++9z6623Yq1lxowZdOnSxeuRIoY+EvdYuC+BFykkdmmTf0jsktgEJ+7at28fAwcOpEuXLlx44YVkZWX5YgN3c630ZC8eKywsFHl1Iold2uQfErskNsHxu1avXk3Xrl357rvveOihhxg5ciSnnHKKBxOGLtxrFREnezHGbDbGrDbGZBljMkq+VtMYM98Ys6Hk71i35okUqampXo/gCIld2uQfErskNsF/d1lrGTduHG3atCE3N5cvvviCp59+2jcbOLi7Vm4/nf4/1tr4o36jeAhYYK1tCiwo+VwppVQFlJOTQ+fOnbnjjjto3749wWCQq666yuuxIprXx8SvB6aUfDwFuMHDWTyhbxvxD23yD4ldEpvg/7tSU1OJj49n7ty5jB49mrlz5/r2dQBurpVrx8SNMZuAXMAC4621E4wxv1przzjqe3KttSd8Sl3aMXGllKrIioqKePrpp3n88cdp3Lgx06dPJzHxmId/K6wTHRN381USf7TW/miMqQ3MN8Z8V9ofNMYMBAYC1KtXj0WLFgHQpEkTatSoQTAYBKBWrVq0aNGCxYsXA1C5cmXatm1LZmbmkeu7JiYmsmPHjiPntm3atCnR0dFH3tdXu3ZtmjVrduSYRnR0NMnJyWRkZLB3714AkpKS2Lp1K9u2bQOgefPmREVFsW7dOgDq1q1L48aNSUtLA6Bq1aokJSWRnp5Ofn4+AMnJyWzatImNGzdSrVo14uLiKCoqIjs7G4D69evToEED0tPTAahevTqJiYmkpaVRUFAAQNu2bVm/fj07d+4EoGXLlhQUFLBhwwYAGjZsSJ06dTj8S09MTAwJCQmkpqZSWFgIwOWXX87atWvZtWsXAIFAgLy8PDZu3AhAo0aNqFmz5pGL3MfGxhIIBEhJScFaizGGdu3aEQwGyc3NBSAhIYFVq1Zx+BdECeu0fft29u3bR2Jioq5ThK/T4VnPO+88MeuUk5PDunXrqFatmqh1io2NpV+/fqxatYqrrrqKf/7zn1xwwQVH/n+8H9dp8+bN7Nu3j5YtW4ZtnU7IWuv6H+Bx4D4gGzi75GtnA9kn+9nWrVtbSRYuXOj1CI6Q2KVN/iGxS1rT7Nmzba1atWyVKlXs5MmTbXFxsdcjhU241wrIsMfZE105Jm6MOc0YU+Pwx8DVwBpgNtCn5Nv6AJ+4MY9SSilvFBQUcPfdd/OXv/yFc845h/Hjx9OnT5+IO3WqX7hyTNwY0wSYVfJpZWCatXaUMaYWMBM4B/g3cJO1NudEtyXtmHh+fj5Vq1b1eoywk9ilTf4hsUtCU3Z2Nt26dSMrK4u7776bZ599luLiYt93/V6418rz94lbazdaawMlf1pYa0eVfH2XtfZKa23Tkr9PuIFLtHXrVq9HcITELm3yD4ldfm6y1jJ58mRat27Nli1bmDNnDmPGjCE6OtrXXcfjZpPXbzGr8A6/SEQaiV3a5B8Su/zatGfPHnr16kW/fv1o06YNwWCQjh07Hvl3v3adiJtNuokrpZRyxIoVK0hISGD69Ok8+eSTfPnll9SvX9/rsUTRTdxjzZs393oER0js0ib/kNjlp6bi4mJGjx7NpZdeysGDB1m8eDGPPvooUVFR//W9fuoqLTeb5J1N32eO9X/UEkjs0ib/kNjll6YdO3bQt29fPv/8c2688UYmTpxIbOzxz+Hll65QuNmkj8Q9dvhECdJI7NIm/5DY5Yem+fPnEwgEWLhwIePGjeODDz444QYO/ugKlZtNuokrpZQql4MHDzJ8+HCuueYaatWqxYoVK7jtttv0vd8u0KfTPVa3bl2vR3CExC5t8g+JXZHatGnTJrp37056ejoDBw7kxRdfpFq1aqX++UjtKg83m3QT91jjxo29HsEREru0yT8kdkVi04wZMxg4cCDGGGbOnMlNN90U8m1EYld5udmkT6d77PDFAqSR2KVN/iGxK5Ka9u3bx6233kq3bt1o0aIFWVlZZdrAIbK6wsXNJt3ElVJKldo333xDYmIikyZNYvjw4aSkpNCoUSOvx6qw9Ol0j0k7Z/BhEru0yT8kdnndZK1l3Lhx3HvvvcTGxjJ//nyuvPLKct+u111OcLPJlQughJO0C6AopVSky8nJYcCAAcyaNYs//elPTJ48mdq1a3s9VoXh+QVQ1PGlp6d7PYIjJHZpk39I7PKqacmSJcTHxzN37lyef/555s6dG9YNXNeqfHQT91h+fr7XIzhCYpc2+YfELrebioqKeOKJJ7jiiiuIjo5m6dKl3HvvvVSqFN5tQ9eqfPSYuFJKqf+wdetWevXqRUpKCr169eK1116jRo0aXo+ljkGPiXusoKCA6Ohor8cIO4ld2uQfErvcapo9ezb9+vWjoKCA1157jZtvvtnR+9O1Ojk9Jh7BNm3a5PUIjpDYpU3+IbHL6ab9+/czZMgQrr/+es4991wyMzMd38BB16q8dBP32Pbt270ewRESu7TJPyR2OdmUnZ1NcnIyr7zyCkOHDiUtLY1mzZo5dn9H07UqHz0mrpRSFZS1lsmTJ3PnnXdSrVo15syZQ8eOHb0eS4VAH4l7LC4uzusRHCGxS5v8Q2JXuJv27NlDz549ueWWW0hKSiIYDHqygetalY9u4h4rKiryegRHSOzSJv+Q2BXOphUrVtCqVStmzpzJU089xfz586lXr17Ybj8Uulblo5u4x7Kzs70ewRESu7TJPyR2haOpuLiY0aNHc+mll1JYWEhKSgqPPPIIUVFRYZiwbHStykePiSulVAWwY8cO+vTpw7x58+jcuTNvvPEGsbGxXo+lykk3cY/Vr1/f6xEcIbFLm/xDYld5mubPn0/v3r3ZvXs3r7/++pFrgEcCXavy0afTPdagQQOvR3CExC5t8g+JXWVpOnjwIA899BBXX301tWrVYsWKFQwaNChiNnDQtSov3cQ9JvHk/yCzS5v8Q2JXqE2bNm3isssu49lnn2XgwIGsWLGCli1bOjRd2elalY8+na6UUsLMmDHjyFPmM2fO5KabbvJ6JOUQfSTuserVq3s9giMkdmmTf0jsKk3Tb7/9xoABA+jWrRstWrQgKysr4jfwirpW4aIXQFFKKQG++eYbunbtSnZ2NsOHD+fxxx/nlFNO8XosFQZ6AZQIlpaW5vUIjpDYpU3+IbHreE3WWl599VUuvvhidu/ezfz58xk1apRvNvCKtFZO0GPiHisoKPB6BEdI7NIm/5DYdaymnJwc+vfvz8cff8x1113H5MmTOeusszyYruwqylo5RR+JK6WUDy1ZsoRAIMCnn37KCy+8wJw5c3y3gavy02PiHissLKRyZXlPiEjs0ib/kNh1uKmoqIinnnqKJ554giZNmjB9+nRat27t9XhlJnmtwkWPiUew9evXez2CIyR2aZN/SOxav349W7dupX379jz++OP07NmTzMxMX2/gIHet3KKbuMd27tzp9QiOkNilTf4hsWvWrFkEAgFWrlzJ1KlTmTp1KjVq1PB6rHKTuFZuNsl6DkMppYTZv38/999/P2PHjiUhIYHp06fTtGlTr8dSEUI3cY9F4mkQw0Filzb5h5Su7777jm7duhEMBhk0aBAvvfQS0dHRXo8VVlLW6mhuNukm7jGJb68AmV3a5B9+77LWMnnyZO68806qVavG3LlziY+PF7eBg//X6lj0LWYVyIYNG7wewRESu7TJP/zctXv3bnr27Mktt9xCUlISwWCQP//5z75uOhGJXW426SaulFIRYvny5bRq1YqZM2cyatQo5s+fT7169bweS0Uw3cQ91rBhQ69HcITELm3yD791FRcX89xzz/HHP/6RoqIiFi9ezMMPP0xUVNSR7/FbU2lJ7HKzSY+Je6xOnTpej+AIiV3a5B9+6tqxYwd9+vRh3rx5dO7cmTfeeIPY2Nj/+j4/NYVCYpebTfpI3GOSzj53NIld2uQffun64osvCAQCpKSk8Prrr/P+++8fcwMH/zSFSmKXm026iSullMsOHDjAgw8+yDXXXMOZZ57JihUrGDRoEMYYr0dTPqNPp3ssJibG6xEcIbFLm/wjkrs2btxI9+7dWb58OYMGDeKFF16gWrVqJ/25SG4qD4ldbjbpBVCUUsol06dPZ9CgQVSqVIk33niDv/3tb16PpHwgYi6AYoyJMsasMsbMLfm8sTEm3RizwRgzwxhzqpvzRILU1FSvR3CExC5t8o9I6/rtt9/o378/3bt3p2XLlmRlZYW8gUdaU7hI7HKzye1j4ncD3x71+bPAi9bapkAu0N/leTxXWFjo9QiOkNilTf4RSV3BYJDExETeeustHnnkEVJSUjj33HNDvp1IagoniV1uNrm2iRtjGgB/BiaWfG6A9sAHJd8yBbjBrXmUUspJ1lrGjh1LUlISu3fv5ssvv+Spp54Sd+1s5S3XjokbYz4A/gHUAO4D+gLLrLXnl/x7Q+Bf1toTnjle2jHx4uJiKlWS9yYBiV3a5B9ed+3atYv+/fvzySefcN111zF58mTOOuusct2m101OkdgV7qYTHRN35VdCY0xHYKe1dqUx5orDXz7Gtx7zNwpjzEBgIEC9evVYtGgRAE2aNKFGjRoEg0EAatWqRYsWLVi8eDEAlStXpm3btmRmZrJnzx4AEhMT2bFjB1u2bAGgadOmREdHs2bNGgBq165Ns2bNjhzTiI6OJjk5mYyMDPbu3QtAUlISW7duZdu2bQA0b96cqKgo1q1bB0DdunVp3LgxaWlpAFStWpWkpCTS09PJz88HIDk5mU2bNrFp0yaqVq1KXFwcRUVFZGdnA1C/fn0aNGhAeno6ANWrVycxMZG0tLQjJ9dv27Yt69evP3Lt2pYtW1JQUHDkvL0NGzakTp06R96zGBMTQ0JCAqmpqUee7rn88stZu3Ytu3btAiAQCJCXl8fGjRsBaNSoETVr1iQzMxOA2NjYI+9rtdZijKFdu3YEg0Fyc3MBSEhIYPXq1Rw8eFDMOm3fvp38/Hxat26t6xTh63T4e8455xxP1mnKlCk8+OCD5ObmMmrUKHr27MnatWvLtU45OTl8++23VK1aVdQ6xcXFsXHjRvbv3+/6Ojnx31NOTg6bN28mPz+fFi1ahG2dTsha6/gfDj0C3wpsBrYD+4B3gV+AyiXfkwzMO9lttW7d2kqycOFCr0dwhMQubfIPL7oOHjxoH3vsMVupUiXbtGlTu3LlyrDevq6Vf4S7Cciwx9kTXXkOw1o73FrbwFrbCOgGfGWt7QksBA6/RLMP8Ikb8yilVDht2bKF9u3bM3LkSHr27MnKlStJSEjweixVAXh9IOJB4F5jzPdALWCSx/O4LhAIeD2CIyR2aZN/uNn1ySefEB8fz6pVq5g6dSpTp06lRo0aYb8fXSv/cLPJ9U3cWrvIWtux5OON1tqLrbXnW2tvstbKuzr8SeTl5Xk9giMkdmmTf7jRtX//fu666y5uuOEGGjVqRGZmJr1793bs/nSt/MPNJq8fiVd4h19IIY3ELm3yD6e7vv32W5KSkhg7diz33HMPS5cupWnTpo7ep66Vf7jZFPKr040xpwH7rbVFDsyjlFIRy1rLW2+9xV133UW1atX49NNPue6667weS1VgJ30kboypZIzpYYz51BizE/gO+MkYs9YY85wxxtlfP4Vr1KiR1yM4QmKXNvmHE127d++mR48e9O/fn0suuYRgMOjqBq5r5R9uNpXm6fSFwHnAcKCutbahtbY2cBmwDHjGGNPLwRlFq1mzptcjOEJilzb5R7i7li9fTqtWrXj//fcZNWoUX3zxBfXq1QvrfZyMrpV/uNlUmk38Kmvtk8Bua23x4S9aa3OstR9aazsDMxybULjDJxOQRmKXNvlHuLqKi4v55z//yR//+EeKi4tZvHgxDz/8MFFRUWG5/VDoWvmHm00n3cSttQdLPpz1+38zxlzyu+9RSikRtm/fzrXXXsuDDz7IDTfcQFZWFpdeeqnXYyn1H0pzTLyLMeYZoIYx5kJjzNG/gk5wbrSKITY21usRHCGxS5v8o7xdX3zxBYFAgCVLljB+/HhmzpzJGWecEabpykbXyj/cbDrpBVCMMfWBK4EXgBVAc+BX4CeglrX2EqeHPJq0C6AopSLHgQMHePTRR3nuuedo0aIFM2bMoEWLFl6PpSq4E10ApTRPp2+z1k4FrrfW/sla2wS4ikPnQ18b3lErnpSUFK9HcITELm3yj7J0bdy4kbZt2/Lcc89x2223sWLFiojawHWt/MPNplK/T9xa+7UxJh7oDnQFdgAXAP0dmq1CONkzIX4lsUub/CPUrvfee49BgwYRFRXFBx98QOfOnR2arOx0rfzDzabSHBNvZowZYYz5jkPnNs8BrrDWJpV8rMrBmGNdkdX/JHZpk3+Utuu3337jlltuoUePHvzhD38gKysrIjdw0LXyEzebSnNMvJhDx8L7W2vX/O7fNpY8ve4aPSaulAqHrKwsunXrxvr163nkkUd47LHHqFw55JNYKuW4ch0TBzpz6Drg840xbxtjOhljTgnngBXZ4QvGSyOxS5v840Rd1lpeeeUVkpKSyMvLY8GCBTz55JMRv4FXxLXyKzebSvPCtlnW2q7A+cDnwCBgqzHmLSDG4fnEy83N9XoER0js0ib/OF7Xrl27uOGGGxgyZAgdOnQgKyuL//mf/3F5urKpaGvlZ242lfoqZtba36y175ZcRvRCDp1ydbVjkymlVBilpKQQCAT4/PPPGTNmDHPmzOGss87yeiylyqU0x8SNPck3leZ7wkXaMfE9e/YQEyPvCQ2JXdrkH0d3FRYW8tRTT/Hkk09y3nnnMX36dBISEjyeMHQVYa2kCHdTeY+JLzTG3GWMOed3N3qqMaa9MWYK0Cccg1ZEOTkyX+AvsUub/ONw15YtW2jfvj0jR46kV69erFy50pcbOMhfK0ncbCrNJn4tUAS8Z4z50RizzhizEdjAofeMv2itnezgjKJt3rzZ6xEcIbFLm/xj8+bNfPLJJwQCAVatWsXbb7/NlClTqFGjhtejlZnktZLGzaaTvhzTWrsfeA14reRV6WcC+dbaX50eTimlQrV//35eeuklPv74Y1q3bs17771H06ZNvR5LKUeU+oVtcOhqZdban3QDD58mTVx9m71rJHZpU+T79ttvSUpK4uOPP+bee+9l6dKlYjZwaWt1mMQuN5tC2sRV+Pn56b0TkdilTZHLWsukSZNITEzkp59+YsaMGTz//POceuqpXo8WNlLW6vckdrnZpJu4xySe6ABkdmlTZNq9ezfdu3dnwIABJCcnEwwGqV27ttdjhZ2EtToWiV0RdbIXpZSKVOnp6bRq1YoPPviAp59+mnnz5nH22Wd7PZZSrin1Jm6MuckYU6Pk40eNMR8ZY/z5Xo0IUqtWLa9HcITELm2KHMXFxTz77LO0bduW4uJilixZwvDhw4mKigL823UiEptAZpebTSc92cuRbzTmG2vtRcaYthy6lvho4OGSq5m5RtrJXoqLi6lUSd4TIhK7tCkybN++nZtvvpn58+dz0003MWHCBM4444z/+B4/dp2MxCaQ2RXupvKe7OWwopK//wyMs9Z+Ash51YhHFi9e7PUIjpDYpU3emzdvHoFAgNTUVCZMmMCMGTP+awMH/3WVhsQmkNnlZlMom/g2Y8x4oCvwmTEmOsSfV0qpMjlw4AD3338/1157LbVr1yYjI4Nbb71V5LWolQpFKJtwF2AecE3J+8RrAvc7MlUFEumXPywriV3a5I3//d//pW3btowePZrbb7+d5cuXExcXd8Kf8UNXqCQ2gcwuN5tCOSZugF5AY2vtEyXnUq9rrV3u5IC/J+2YuFLq+KZNm8Ztt91GVFQUkyZN4sYbb/R6JKVcF65j4q8Bl3DofOkAecCr5ZytwsvMzPR6BEdI7NIm9+zdu5d+/frRs2dPLrroIrKyskLawCO1qzwkNoHMLjebQtnEk6y1g4H9ANbaXPSFbeW2Z88er0dwhMQubXJHVlYWiYmJTJkyhb///e8sWrSIc889N6TbiMSu8pLYBDK73GwKZRM/aIyJAiyAMeYsoNiRqZRSFY61lpdffpmkpCTy8vJYsGABTzzxhMhjpkqFSyjHxHty6JXpCcAU4G/Ao9ba950b779JOya+d+9eqlev7vUYYSexS5uc88svv3DLLbcwZ84cOnbsyFtvvcWZZ55Z5tuLlK5wktgEMrvC3RSWY+LW2neBBzh0opefgBvc3sAl2rFjh9cjOEJilzY5Y9GiRcTHxzNv3jzGjBnD7Nmzy7WBQ2R0hZvEJpDZ5WZTqJci/c5a+6q1dqy19lunhqpItmzZ4vUIjpDYpU3hVVhYyIgRI2jfvj3VqlVj2bJl3H333WF577eulX9I7HKzqdQHm4wxicAjwLklP2cAa629yKHZlFJCbdmyhR49epCamkqfPn0YO3asuKdUlXJDKK8YeWNmNTEAACAASURBVJdDJ3dZjb6gLWyaNm3q9QiOkNilTeExa9Ys+vfvz8GDB3nnnXfo2bNn2O9D18o/JHa52RTKJv6ztXa2Y5NUUNHR0V6P4AiJXdpUPvn5+dx333289tprtG7dmunTp3P++ec7cl+6Vv4hscvNplCOiT9mjJlojOlujLnx8B/HJqsg1qxZ4/UIjpDYpU1l9+2335KUlMRrr73GsGHDWLp0qWMbOOha+YnELjebQnkk3g+4ADiF/3863QIfhXsopZQM1lomTZrEkCFDqF69Op999hl/+tOfvB5LKTFC2cQD1to/ODZJBVW7dm2vR3CExC5tCs3u3bsZOHAgM2fO5Morr+Ttt9/m7LPPduz+jqZr5R8Su9xsCuVkL28AL1pr1zk70olJO9lLYWGhyDNSSezSptJbtmwZ3bt3Z8uWLTz11FM88MADVKrk3pWLda38Q2JXuJvCdQGUtkDQGJNtjPnGGLPaGPNNeEasuFJTU70ewRESu7Tp5IqLi3nmmWdo27btkdt/6KGHXN3AD9+vNBKbQGaXm02h/KpwDSXvDXdoFqWUj23fvp3evXvz5Zdf0qVLF8aPH88ZZ5zh9VhKiXbSTdwYk2qtbQus5T838MMbeoxDs1UIEt9eATK7tOn4Pv/8c26++Wb27t3LG2+8Qf/+/cNy5rWy0rXyD4ldbjaV+ph4pJB2TFwpPztw4AAPP/wwzz//PC1btmTGjBnExcV5PZZSopT7mLg5pGE5BqhijFlujAkaY9YaY0aWfL2xMSbdGLPBGDPDGFPhrk8u9RcSiV3a9J++//57/vjHP/L8889zxx13sHz58ojZwHWt/ENil5tNpdrE7aGH6x+X434KgPbW2gAQD1xrjLkEeJZDr3hvCuQC/ctxH760d+9er0dwhMQubfp/7777LgkJCXz//fd89NFHvPrqq1StWjXM05WdrpV/SOxysymUl4wuM8a0Kcud2EMOV51S8scC7YEPSr4+BbihLLevlHLH3r176du3L7169eKiiy4iGAzy17/+1euxlKqwQnmf+DqgGfAD8BshXsXMGBMFrATOB14FngOWWWvPL/n3hsC/rLUtT3Q70o6J5+fnR9QjmHCR2FXRm1atWkW3bt3YsGEDjz76KCNGjIjY9/dW9LXyE4ld4W460THxUP4LLNe5Eq21RUC8MeYMYBZw4bG+7Vg/a4wZCAwEqFevHosWLQKgSZMm1KhRg2AwCECtWrVo0aIFixcvBqBy5cq0bduWzMxM9uzZA0BiYiI7duw4cr3Xpk2bEh0dfeRct7Vr16ZZs2ZH3ucXHR1NcnIyGRkZR54iSUpKYuvWrWzbtg2A5s2bExUVxbp1h86DU7duXRo3bkxaWhoAVatWJSkpifT0dPLz8wFITk5m06ZN/PDDD0RHRxMXF0dRURHZ2dkA1K9fnwYNGpCeng5A9erVSUxMJC0tjYKCAgDatm3L+vXr2blzJwAtW7akoKCADRs2ANCwYUPq1Klz5PhMTEwMCQkJpKamUlhYCMDll1/O2rVr2bVrFwCBQIC8vDw2btwIQKNGjahZsyaZmZkAxMbGEggESElJwVqLMYZ27doRDAbJzc0FICEhge+++459+/aJWaft27dTUFBAq1atKtw6paSk8OGHHzJhwgTOOussXn/99SP/20fiOgGcfvrp1K1bV8w65eTkkJ2dTXR0tJj/ngDi4uLYunXrkTklrNPmzZspKCjgwgsvDNs6nZC19oR/OPTI+Y/H+PplwHkn+/nj3OZjHLqs6S9A5ZKvJQPzTvazrVu3tpIsXLjQ6xEcIbGrIjb9/PPPtlOnThawHTt2tD///LM7g5VTRVwrv5LYFe4mIMMeZ08szTHxMUDeMb6eX/JvJ2WMOavkETjGmKrAVcC3wELgbyXf1gf4pDS3p5Ry3qJFiwgEAsybN4+XXnqJ2bNnc+aZZ3o9llLqKKXZxBtZa//r9KrW2gygUSnv52xgYclpWlcA8621c4EHgXuNMd8DtYBJpbw9MZo3b+71CI6Q2FVRmgoLCxkxYgTt27enevXqLFu2jCFDhnh68pZQVZS1kkBil5tNpTkmXuUE/1aqI/clvwS0OsbXNwIXl+Y2pIqKivJ6BEdI7KoITf/+97/p2bMnqamp9O3bl1deeYXq1at7NF3ZVYS1kkJil5tNpXkkvsIYc+vvv2iM6c+hV5urcjj8ohBpJHZJb5o1axbx8fEEg0Heffdd3nrrLV9u4CB/rSSR2OVmU2keiQ8FZhljevL/m3YicCqgbxBVyufy8/MZNmwY48aNIzExkenTp3Peeed5PZZSqhROuolba3cAlxpj/gc4/B7uT621Xzk6WQVRt25dr0dwhMQuiU179+4lKSmJ1atXc9999zFq1ChOPdX/Zz+WuFYSm0Bml5tNegEUjxUUFIi8io/ELklN1lomTpzI3XffTfXq1Zk6dSrXXnut12OFjaS1OkxiE8jsCndTuS+Aopxz+MQI0kjsktL066+/0rVrVwYOHEhcXBzBYFDUBg5y1upoEptAZpebTZF5zkSllCPS0tLo0aMHW7du5ZlnnqFNmzacffbZXo+llCqjk27ixph7T/Tv1toXwjdOxSPtnMGHSezyc1NxcTHPPvssf//732nYsCFLlizhkksuOXJ6S2n8vFbHI7EJZHa52XTSY+LGmMdKPmwOtAFml3zeCVhsrR3g3Hj/TdoxcaWc9tNPP9G7d28WLFhAly5dGD9+PGeccYbXYymlSqlcx8SttSOttSOBM4EEa+0wa+0woDXQILyjVjxSHwlJ7PJj07/+9S8CgQBLly5l4sSJTJ8+/T82cD82lYbELolNILPLzaZQXth2DnDgqM8PUPrTrqrjOHx1H2kkdvmp6cCBAwwbNozrrruOunXrkpGRQf/+/f/r1Kl+agqFxC6JTSCzy82mUF7Y9jaw3Bgzi0OXDP0rMNWRqZRSZfb999/TrVs3Vq5cyR133MHo0aNFHndUSoX4PnFjTAKHLkEKh46Hr3JkqhOQdkxc4nskQWaXH5reeecdbr/9dk455RQmTZrEX/964pMq+qGpLCR2SWwCmV0R+T5xc+h5uDjgdGvtS8AuY0yFvnhJOGzatMnrERwhsSuSm/bu3Uvfvn3p3bs38fHxZGVlnXQDh8huKg+JXRKbQGaXm02hHBN/DUgGupd8nge8GvaJKpjt27d7PYIjJHZFatOqVato3bo1b7/9NiNGjGDhwoWcc845pfrZSG0qL4ldEptAZpebTaEcE0+y1iYYY1YBWGtzjTH+P8myUj5lreXll1/mgQce4KyzzuKrr76iXbt2Xo+llHJRKJv4QWNMFIde1IYx5iyg2JGpKpC4uDivR3CExK5Iavrll1/o168fc+fOpVOnTrz55puceeaZId9OJDWFk8QuiU0gs8vNplCeTn8ZmAXUNsaMAlKBpx2ZqgIpKiryegRHSOyKlKaFCxcSCAT44osvePnll/nkk0/KtIFD5DSFm8QuiU0gs8vNplJv4tbad4EHgH8APwE3WGvfd2qwiiI7O9vrERwhscvrpsLCQv7+979z5ZVXUqNGDdLT07nrrrv+673fofC6ySkSuyQ2gcwuN5tCugCKtfY74DuHZlFKHccPP/xAz549+frrr+nXrx8vv/wy1atX93ospZTHQnmL2RRjzBlHfR5rjHnTmbEqjvr163s9giMkdnnV9NFHHxEfH88333zDtGnTePPNN8O2gUtcJ5DZJbEJZHa52RTKMfGLrLW/Hv7EWpsLtAr/SBVLgwYyTz8vscvtpvz8fG6//XY6d+7M+eefz6pVq+jevfvJfzAEEtcJZHZJbAKZXW42hbKJVzLGxB7+xBhTE70eeblJPPk/yOxys2nt2rVcfPHFvP7669x///18/fXXnHfeeWG/H4nrBDK7JDaBzC43m0LZhJ8H0owxh1/MdhP66nSlwspayxtvvMHQoUOpUaMGn3/+Oddcc43XYymlIlSpN3Fr7VRjTAbQvuRLN1pr1zkzVsUh9cVJErucbvr1118ZOHAg77//Ph06dGDq1KnUrVvX0fuUuE4gs0tiE8jscrOp1BdAMcaMONbXrbVPhHWik5B2ARSlANLS0ujevTvbtm1j1KhR3HfffVSqFMrRLqWUVGG5AArw21F/ioA/odcTL7e0tDSvR3CExC4nmoqKinj66ae57LLLMMaQmprKAw884NoGLnGdQGaXxCaQ2eVmUyhPpz9/9OfGmNHA7LBPVMEUFBR4PYIjJHaFu+nHH3/k5ptvZsGCBXTt2pXx48dz+umnh/U+TkbiOoHMLolNILPLzabyvLq8GtAkXIMoVZF89tln9OnTh99++42JEydyyy23lOvMa0qpiimUY+KrKbn4CRAFnAU8aa19xaHZjknaMfHCwkIqV5b3Tj2JXeFoOnDgAMOHD+eFF17goosuYvr06Vx44YVhmjB0EtcJZHZJbAKZXeFuCtcx8Y5Ap5I/VwP13N7AJVq/fr3XIzhCYld5mzZs2MCll17KCy+8wODBg0lPT/d0AweZ6wQyuyQ2gcwuN5tOuokbY+41xtwLdD7qT1dgSMnXVTns3LnT6xEcIbGrPE3vvPMOCQkJbNy4kVmzZjF27FiqVKkSxunKRuI6gcwuiU0gs8vNptI83q9R8ndzoA3//2K2TsBiJ4ZSSoq9e/cyePBgpk6dymWXXca7775Lw4YNvR5LKSXESTdxa+1IAGPMF0CCtTav5PPHAb0UaTm1bNnS6xEcIbEr1KbMzEy6devG//7v//LYY4/x6KOPRtyxP4nrBDK7JDaBzC43m0I5Jn4OcOCozw+g7xMvN4lvrwCZXaVtstYyZswYLrnkEvbt28dXX33F448/HnEbOMhcJ5DZJbEJZHa52RTKJv42sNwY87gx5jEgHZjizFgVx4YNG7wewRESu0rT9PPPP9OpUyfuuece/vSnPxEMBmnXrp0L05WNxHUCmV0Sm0Bml5tNoZzsZZQx5l/AZSVf6metXeXMWEr5z8KFC+nZsye7du3ilVdeYfDgwfreb6WUo0J6fs9amwlkOjRLhST1RU4Su47XVFhYyOOPP87TTz9Ns2bN+Oyzz4iPj3d5urKRuE4gs0tiE8jscrOp1E+nm0N6Hb4QijHmHGPMxc6NVjHUqVPH6xEcIbHrWE0//PAD7dq1Y9SoUfTr14+VK1f6ZgMHmesEMrskNoHMLjebQjkm/hqQDHQv+TwPeDXsE1Uwks4+dzSJXb9v+vDDD4mPj2f16tVMmzaNSZMmcdppp3k0XdlIXCeQ2SWxCWR2udkUyiaeZK0dDOwHsNbmAqc6MpVSESw/P5/bbruNv/3tbzRt2pRVq1bRvXv3k/+gUkqFWSib+EFjTBQl5083xpwFFDsyVQUSExPj9QiOkNgVExPD2rVradOmDePHj+eBBx4gNTWV8847z+vRykziOoHMLolNILPLzaZQLoDSk0OnW03g0FvL/gY8aq119YQv0i6AovzBWsuECRMYOnQoMTExTJ06lWuuucbrsZRSFUBYLoBirX0XeAD4B/AjcBP6dHq5paamej2CIyR15ebm0qVLF2677TYuu+wygsGgmA1c0jodTWKXxCaQ2eVmU2kugBJjjBlujBnLobO2vVbyc3OALg7PJ15hYaHXIzhCStfSpUuJj4/n448/ZuDAgXz++efUrVvX67HCRso6/Z7ELolNILPLzabSvE/8bSAXSAMGAPdz6BH49dbaLAdnU8ozRUVFPPvss4wYMYJzzjmH1NRU8vPzqVQplJeRKKWUs056TNwYs9pa+4eSj6OAX4BzDl8IxW3SjokXFxeL3Bj83PXjjz/Su3dvvvrqK7p168brr7/O6aef7uum45HYBDK7JDaBzK5wN5X3mPjBwx9Ya4uATaFu4MaYhsaYhcaYb40xa40xd5d8vaYxZr4xZkPJ37Gh3K4Ea9eu9XoER/i169NPPyUQCLBs2TImTZrEtGnTOP300wH/Np2IxCaQ2SWxCWR2udlUmk08YIzZU/InD7jo8MfGmD2lvJ9CYJi19kLgEmCwMSYOeAhYYK1tCiwo+bxC2bVrl9cjOMJvXQUFBdx777107NiRevXqkZGRwS233PIf5z73W1NpSGwCmV0Sm0Bml5tNpbmeeFR578Ra+xPwU8nHecaYb4H6wPXAFSXfNgVYBDxY3vtTKhQbNmygW7duZGZmcuedd/Lcc89RpUoVr8dSSqmTcv0Cx8aYRkArDl3KtE7JBo+19idjTG235/FaIBDwegRH+KXr7bff5o477uDUU0/l448/5vrrrz/u9/qlKRQSm0Bml8QmkNnlZpOrm7gxpjrwITDUWruntJdpNMYMBAYC1KtXj0WLFgHQpEkTatSoQTAYBKBWrVq0aNGCxYsXA1C5cmXatm1LZmYme/YceuY/MTGRHTt2sGXLFgCaNm1KdHQ0a9asAaB27do0a9bsyPv8oqOjSU5OJiMjg7179wKQlJTE1q1b2bZtGwDNmzcnKiqKdevWAVC3bl0aN25MWloaAFWrViUpKYn09HTy8/MBSE5OZtOmTfz73//m1FNPJS4ujqKiIrKzswGoX78+DRo0ID09HYDq1auTmJhIWlrakQvOt23blvXr17Nz504AWrZsSUFBwZFr2TZs2JA6deocOY9vTEwMCQkJpKamHnkLxOWXX87atWuPPP0TCATIy8tj48aNADRq1IiaNWuSmXno4nWxsbEEAgFSUlKw1mKMoV27dgSDQXJzcwFISEhg06ZNR9YlEtdp3759jB8/ntmzZ3PRRRfxyCOPHHnr2LHWafv27Rw4cID4+Hhdpwj+72n79u0AnHnmmezfv1/MOuXk5LB+/XpOPfVUUesUFxfHjh07jrRIWKfNmzdz4MABLrjggrCt0wlZa135A5wCzAPuPepr2cDZJR+fDWSf7HZat25tJVm4cKHXIzgikrsyMjLs+eefbytVqmQff/xxW1hYWKqfi+SmspLYZK3MLolN1srsCncTkGGPsye68rp+c+gh9yTgW2vtC0f902ygT8nHfYBP3JhHVUzWWl588UWSk5PJz89n4cKFPPbYY0RFlftlH0op5Qm3nk7/I9AbWG2MOXyCmIeBZ4CZxpj+wL85dCrXCqVRo0Zej+CISOv6+eef6du3L5999hnXX389kyZNolatWiHdRqQ1hYPEJpDZJbEJZHa52eTKJm6tTQWOdwD8SjdmiFQ1a9b0egRHRFLXV199Ra9evcjJyWHs2LHccccdlPb1GEeLpKZwkdgEMrskNoHMLjebZJ0mx4cOv3BCmkjoKiws5JFHHuGqq64iJiaG9PR0Bg8eXKYNHCKjKdwkNoHMLolNILPLzSbX32KmlBs2b95Mjx49SEtL45ZbbuHll1/mtNNO83ospZQKK93EPRYbK/NMs152ffDBBwwYMIDi4mLee+89unXrFpbblbhWEptAZpfEJpDZ5WbTSS+AEmmkXQBFhc++ffu45557mDBhAhdffDHvvfceTZo08XospZQql/JeAEU5KCUlxesRHOF215o1a7j44ouZMGECDzzwAKmpqWHfwCWulcQmkNklsQlkdrnZpE+ne8xvz4SUlltd1lomTJjA0KFDiYmJYd68eVx99dWO3Zc0EptAZpfEJpDZ5WaTPhL3WFlfKR3p3OjKzc3lpptu4rbbbuPyyy/nm2++cWwDB5lrJbEJZHZJbAKZXW426TFx5Utff/01PXr04Mcff+Tpp59m2LBhVKqkv5MqpeTRY+IR7PDJ8aVxqquoqIhRo0bRrl07KleuzNdff83999/vygYuca0kNoHMLolNILPLzSY9Ju6xw1fAkcaJrh9//JFevXqxcOFCunfvzrhx4zj99NPDfj/HI3GtJDaBzC6JTSCzy80m3cSVL3z66af07duXffv28eabb9K3b1+Rx9KUUioU+nS6xxISErwewRHh6iooKOCee+6hY8eO1K9fn5UrV9KvXz9PNnCJayWxCWR2SWwCmV1uNukm7rGcnByvR3BEOLrWr19PcnIyY8aM4a677mLZsmVccMEFYZiubCSulcQmkNklsQlkdrnZpJu4xzZv3uz1CI4ob9fUqVNJSEjghx9+4OOPP+bll1+mSpUq4RmujCSulcQmkNklsQlkdrnZpJu4iih5eXn07t2bPn360Lp1a4LBINdff73XYymlVETSTdxjUs/tXZaulStXkpCQwLRp0xg5ciRfffUVDRo0cGC6spG4VhKbQGaXxCaQ2eVmk27iHqtRo4bXIzgilK7i4mJeeOEFkpOT2b9/PwsXLmTEiBFERUU5OGHoJK6VxCaQ2SWxCWR2udmkm7jHJJ7oAErftXPnTjp16sSwYcO47rrryMrK4vLLL3d4urKRuFYSm0Bml8QmkNnlZpNu4sozCxYsIBAIsGDBAsaOHcusWbOoVauW12MppZRv6CbuMamb1om6Dh48yMMPP0yHDh0444wzSE9PZ/DgwRF/8haJayWxCWR2SWwCmV1uNukFUDxWXFws8sIdx+vavHkz3bt3Z9myZfTv35+XXnqJ0047zYMJQydxrSQ2gcwuiU0gsyvcTXoBlAi2ePFir0dwxLG63n//feLj41m3bh3Tp09n4sSJvtnAQeZaSWwCmV0Sm0Bml5tNuokrx+3bt49BgwbRpUsXLrjgAlatWkXXrl29HksppXxPN3GPVa4s8xo0h7vWrFlDmzZtmDBhAg8++CBLlizx7ftCJa6VxCaQ2SWxCWR2udmkx8SVI6y1jB8/nnvuuYfTTz+dt99+mw4dOng9llJK+Y4eE49gmZmZXo8Qdrm5uVx11VXcfvvttGvXjmAwKGIDl7hWEptAZpfEJpDZ5WaTbuIe27Nnj9cjhNXXX39NfHw8KSkpPPfcc3z22WfUqVPH67HCQtpagcwmkNklsQlkdrnZpJu4CouioiKeeuop2rVrR+XKlRk7diz33XefuLeOKKVUJNFj4h7bu3cv1atX93qMctm2bRu9evVi0aJF9OjRg3HjxlGpUiXfd/2ehLX6PYlNILNLYhPI7Ap3kx4Tj2A7duzweoRymTt3LoFAgOXLl/PWW2/xzjvvEBMT4/uuY9Em/5DYJbEJZHa52aSbuMe2bNni9QhlUlBQwNChQ+nUqRMNGjRg5cqV9O3b98ipU/3adSLa5B8SuyQ2gcwuN5t0E1chW79+PcnJybz00ksMGTKEZcuWccEFF3g9llJKVTjy3mXvM02bNvV6hFKz1jJ16lQGDx5MdHQ0n3zyCX/5y1+O+b1+6iotbfIPiV0Sm0Bml5tN+kjcY9HR0V6PUCp5eXn07t2bvn37kpiYSDAYPO4GDv7pCoU2+YfELolNILPLzSbdxD22Zs0ar0c4qYyMDFq1asV7773HyJEjWbBgAQ0aNDjhz/ihK1Ta5B8SuyQ2gcwuN5t0E1fHVVxczPPPP8+ll17KgQMHWLRoESNGjCAqKsrr0ZRSSqHHxD1Xu3Ztr0c4pp07d9KnTx8+//xzbrjhBiZNmkTNmjVL/fOR2lUe2uQfErskNoHMLjeb9GQvHissLIy4q/gsWLCAXr16kZubywsvvMDtt99+5K1jpRWJXeWlTf4hsUtiE8jsCneTnuwlgqWmpno9whEHDx7k4YcfpkOHDsTGxrJ8+XLuuOOOkDdwiKyucNEm/5DYJbEJZHa52STr1x9VZps3b6Z79+4sW7aMAQMGMGbMGE477TSvx1JKKXUCuol7LBLeXvH+++9z6623Yq1l+vTpdO3atdy3GQld4aZN/iGxS2ITyOxys0mPiVdg+/btY+jQobzxxhskJSXx3nvv0bhxY6/HUkopdRQ9Jh7BvPqFZPXq1bRp04aJEyfy0EMPsWTJkrBu4BJ/0dIm/5DYJbEJZHa52aSbuMf27t3r6v1Zaxk3bhwXX3wxOTk5zJs3j3/84x+ccsopYb0ft7vcoE3+IbFLYhPI7HKzSTfxCiQnJ4fOnTtzxx13cMUVVxAMBunQoYPXYymllCojPSbusfz8fKpWrer4/aSmptKjRw9++uknnnnmGe655x4qVXLudzi3utykTf4hsUtiE8jsCneT58fEjTFvGmN2GmPWHPW1msaY+caYDSV/x7oxS6TZunWro7dfVFTEk08+Sbt27Tj11FNZunQpw4YNc3QDB+e7vKBN/iGxS2ITyOxys8mtp9MnA9f+7msPAQustU2BBSWfVzjbtm1z9LavuuoqRowYQbdu3cjMzKRNmzaO3d/v71sabfIPiV0Sm0Bml5tNrmzi1trFQM7vvnw9MKXk4ynADW7MUlHMmTOHQCDA8uXLeeutt3jnnXeIiYnxeiyllFJh5OUL2+pYa38CKPlb3lnwS6F58+Zhvb2CggLuvvtu/vKXv9CwYUMyMzPp27dvmU6dWh7h7ooE2uQfErskNoHMLjebfHHGNmPMQGAgQL169Vi0aBEATZo0oUaNGgSDQQBq1apFixYtWLx4MQCVK1embdu2ZGZmsmfPHgASExPZsWMHW7ZsAaBp06ZER0cfuf5r7dq1adas2ZFz30ZHR5OcnExGRsaRtw0kJSWxdevWI0+ZNG/enKioKNatWwdA3bp1ady4MWlpaQBUrVqVpKQk0tPTyc/PByA5OZlNmzaxdetWsrOziYuLo6ioiOzsbADq169PgwYNSE9PB6B69eokJiaSlpZGQUEBAG3btmX9+vXs3LkTgCpVqjBgwADWrl3LjTfeyJNPPkn9+vWP/O8VExNDQkICqampFBYWAnD55Zezdu1adu3aBUAgECAvL4+NGzcC0KhRI2rWrElmZiYAsbGxBAIBUlJSsNZijKFdu3YEg0Fyc3MBSEhI4JdffjnSImGdtm/fTmFhIVFRUeVep5YtW1JQUMCGDRsAaNiwIXXq1Dny3lJdp/KtExz6/xOAmHXKycnh+++/Jzs7W9Q6xcXFsXv3blHrtHnzZgoLCzl4pMR85QAAEapJREFU8GDY1umErLWu/AEaAWuO+jwbOLvk47OB7NLcTuvWra0kCxcuLPdtFBcX27feesuedtpptlatWnb27NnlH6ycwtEVabTJPyR2SWyyVmZXuJuADHucPdHLp9NnA31KPu4DfOLhLL61Z88eevXqRb9+/UhMTCQYDNKpUyevx1JKKeUCt95i9h6QBjQ3xmw1xvQHngE6GGM2AB1KPq9w6tatW+afzcjIICEhgenTp/PEE0+wYMEC6tevH8bpyq48XZFKm/xDYpfEJpDZ5WaTnuzFYwUFBSFf8aa4uJgXX3yR4cOHU7duXaZNm0bbtm0dmrBsytIV6bTJPyR2SWwCmV3hbvL8ZC/q+A6/CKS0du7cyZ///Gfuu+8+OnbsSFZWVsRt4BB6lx9ok39I7JLYBDK73GzSTdxHvvzySwKBAAsXLuS1117jww8/pGbNml6PpZRSyiO6iXusNOfXPXjwIMOHD+fqq68mNjaWFStWcPvtt7v+3u9QSDsXMmiTn0jsktgEMrvcbNJj4hFu06ZN9OjRg2XLlnHrrbcyZswYqlWr5vVYSimlXKLHxCPY4ZMaHMvMmTOJj49n3bp1zJgxgwkTJvhmAz9Rl19pk39I7JLYBDK73GzSTdxjh89kdLR9+/Zx66230rVrV+Li4sjKyqJLly4eTFd2x+ryO23yD4ldEptAZpebTbqJR5jVq1eTmJjIpEmTGD58OIsXL6Zx48Zej6WUUioC6TFxjx1+P6G1lnHjxnHvvfcSGxvL22+/zVVXXeX1eGWm7/30B4lNILNLYhPI7NL3iVcgmzZtIicnh86dOzN48GDat29PMBj09QYOh7qk0Sb/kNglsQlkdrnZpJu4x7744gvi4+OZO3cuo0ePZu7cudSu7f+rsh6+SpEk2uQfErskNoHMLjebfHEpUomKiooYNWoUI0eOpHHjxixdupTExGM+W6KUUkodk27iHti6dSu9evUiJSWFzp078+abbxITE+P1WGEVFxfn9Qhhp03+IbFLYhPI7HKzSTdxl82ePZt+/fpRUFDA5MmT6dChg7gNHA490yCNNvmHxC6JTSCzy80mPSbukv379zNkyBCuv/56zj33XFauXEmfPn1Yv36916M5Ijs72+sRwk6b/ENil8QmkNnlZpNu4i7Izs4mOTmZV155hbvvvpu0tDSaN2/u9VhKKaV8Tp9Od5C1lilTpnDnnXdSpUoV5syZQ8eOHf/je+rXr+/RdM6S2KVN/iGxS2ITyOxys0kfiTtkz5499OrVi379+tGmTRuCweB/beAADRo08GA650ns0ib/kNglsQlkdrnZpJu4A1asWEGrVq2YPn06Tz75JF9++eVxfzOTePJ/kNmlTf4hsUtiE8js0gug+FRxcTGjR4/m0ksvpbCwkMWLF/Poo48SFRXl9WhKKaUE0mPiYbJjxw769OnDvHnzuPHGG5k4cSKxsbEn/bnq1au7MJ37JHZpk39I7JLYBDK73GzSC6CEwfz58+nduze//vorY8aMYdCgQRhjvB5LKaWUAHoBFIccPHiQhx56iKuvvppatWqxYsUKbrvttpA28LS0NAcn9I7ELm3yD4ldEptAZpebTfp0ehlt2rSJ7t27k56ezsCBA3nxxRepVq1ayLdTUFDgwHTek9ilTf4hsUtiE8jscrNJN/EymDFjBgMHDsQYw8yZM7npppu8HkkppVQFpMfEQ/Dbb79x9913M2nSJJKTk5k2bRqNGjUq120WFhZSubK836UkdmmTf0jsktgEMrvC3aTHxMPgm2++ITExkTfffJPhw4eTkpJS7g0cEHvudIld2uQfErskNoHMLjebdBM/CWstr776KhdffDG//vor8+fP5+mnn+aUU04Jy+3v3LkzLLcTaSR2aZN/SOyS2AQyu9xs0k38BHJycrjxxhu58847ad++PcFgkCuvvNLrsZRSSilAN/HjWrJkCYFAgE8//ZTnn3+euXPnUrt27bDfT8uWLcN+m5FAYpc2+YfELolNILPLzSbdxH+nqKiIkSNHcsUVV1ClShWWLl3KvffeS6VKzvxPJfHtFSCzS5v8Q2KXxCaQ2eVmk27iR9m6dStXXnkljz/+OD169CAzM5PExGO+IDBsNmzY4Ojte0Vilzb5h8QuiU0gs8vNJlmv6y+H2bNn069fPwoKCpgyZQo333yz1yMppZRSJ1ThH4nv37+fIUOGcP3113PuueeSmZnp6gbesGFD1+7LTRK7tMk/JHZJbAKZXW42VehN/LvvvuOSSy7hlVdeYejQoaSlpdGsWTNXZ6hTp46r9+cWiV3a5B8SuyQ2gcwuN5sq9Cb+8MMPs23bNubMmcOLL75IdHS06zNE2hXZwkVilzb5h8QuiU0gs8vNpgp9TPz111/n4MGD1K9f3+tRlFJKqZBV6E3cifd9hyomJsbrERwhsUub/ENil8QmkNnlZpNeAEUppZSKYHoBlAiWmprq9QiOkNilTf4hsUtiE8jscrNJN3GPFRYWej2CIyR2aZN/SOyS2AQyu9xs0k1cKaWU8ik9Ju6x4uJix87L7iWJXdrkHxK7JDaBzK5wN+kx8Qi2du1ar0dwhMQubfIPiV0Sm0Bml5tNuol7bNeuXV6P4AiJXdrkHxK7JDaBzC43m3QTV0oppXxKN3GPBQIBr0dwhMQubfIPiV0Sm0Bml5tNuol7LC8vz+sRHCGxS5v8Q2KXxCaQ2eVmk27iHtu4caPXIzhCYpc2+YfELolNILPLzSbdxJX6v/buP1iqso7j+PsjJqJhJKiZkCCZZJaCP9L80Tg6pYyJmJVlxWSjY9kU49gE0TSO/WGmOdWkWZKjFqVTiVFpQ5KmVqCCIBdBkQEng9DM/IGkCN/+OM/m4bZn2b1775571s9rZuee++zZc77f8+yeZ5/nnD3HzKyiKvc7cUlPA0+UHUc/GgX8s+wgBkA35uWcqqMb8+rGnKA78+rvnPaPiL3qPVG5RrzbSHqw6Ef8VdaNeTmn6ujGvLoxJ+jOvDqZk4fTzczMKsqNuJmZWUW5ES/fj8oOYIB0Y17OqTq6Ma9uzAm6M6+O5eRj4mZmZhXlnriZmVlFuRHvIEljJN0laaWkFZK+lMovkfR3SUvTY3LZsbZC0jpJy1PsD6ayPSX9QdLq9PfNZcfZLEkH5epiqaTnJU2vYj1Jul7SU5J6cmV160aZ70l6XNLDkiaVF3mxgpyukLQqxT1X0ohUPlbS5lydXVte5I0V5FX4npM0M9XVo5I+WE7UjRXkdEsun3WSlqbyStRVg/14OZ+riPCjQw9gX2BSmh4OPAYcDFwCXFx2fG3ktQ4Y1avsW8CMND0DuLzsOPuY2xDgH8D+Vawn4ARgEtCzo7oBJgN3AAKOBhaVHX8LOX0A2DlNX57LaWx+vsH8KMir7nsu7TeWAUOBccAaYEjZOTSTU6/nvw18vUp11WA/Xsrnyj3xDoqIDRGxJE2/AKwE9is3qgEzBbgxTd8InFFiLO04CVgTEZW8wFBE3AP8q1dxUd1MAW6KzEJghKR9OxNp8+rlFBHzI+LV9O9CYHTHA2tTQV0VmQLcHBEvR8Ra4HHgqAELro8a5SRJwEeBn3c0qDY12I+X8rlyI14SSWOBicCiVPSFNNRyfZWGnpMA5ktaLOn8VLZPRGyA7E0P7F1adO05m+13MlWup5qiutkP+Ftuviep5pfMc8l6PjXjJD0k6U+Sji8rqDbUe891Q10dD2yMiNW5skrVVa/9eCmfKzfiJZD0RuBXwPSIeB74ATAeOAzYQDbEVCXHRsQk4FTgQkknlB1Qf5C0C3A68ItUVPV62hHVKavUz1ckzQJeBeakog3A2yJiInAR8DNJe5QVXx8UvecqX1fAx9n+C3Kl6qrOfrxw1jpl/VZXbsQ7TNIbyCp+TkTcChARGyNia0RsA65jEA6LNRIR69Pfp4C5ZPFvrA0Zpb9PlRdhn50KLImIjVD9esopqpsngTG5+UYD6zscW59JmgacBpwT6WBkGm5+Jk0vJjt2/I7yomxNg/dc1etqZ+BM4JZaWZXqqt5+nJI+V27EOygdA/oxsDIirsqV54+PTAV6er92sJK0u6ThtWmyE4x6gHnAtDTbNODX5UTYlu16ClWup16K6mYe8Ol0Nu3RwHO14cHBTtIpwFeA0yPipVz5XpKGpOkDgAOBytz7ssF7bh5wtqShksaR5XV/p+Nrw8nAqoh4slZQlboq2o9T1ueq7DP9Xk8P4DiyYZSHgaXpMRn4CbA8lc8D9i071hZyOoDsLNllwApgViofCSwAVqe/e5Yda4t57QY8A7wpV1a5eiL7ErIB2ELWI/hsUd2QDftdTdYDWg4cUXb8LeT0ONlxx9rn6to074fT+3IZsAT4UNnxt5hX4XsOmJXq6lHg1LLjbzanVH4DcEGveStRVw3246V8rnzFNjMzs4rycLqZmVlFuRE3MzOrKDfiZmZmFeVG3MzMrKLciJuZmVWUG3EzM7OKciNuZmZWUW7EzZogaWu6x/EySUskva+fl/9ibvov/bTMt0i6WdIaSY9Iul1S3ctY9tc624mhjXWMkPT5Fl8zLN1ko3aFsNGSPtZg/l0k3ZMuF2o2aLgRN2vO5og4LCIOBWYClw3UiiKi7S8I6dKQc4G7I2J8RBwMfBXYp/d8knbqj3X2NYZ+MAJoqREnu9PZrRGxNf1/Etl9r+uKiFfIrsJV2NCblcGNuFnr9gCeBZB0W7oF64rabVjT9eR/l3rtPfkenqRPSro/9ep/WOsJ5kl6UdJYSSslXZeWPV/SsBaWcyKwJSKurRVExNKIuDe37GvILm85pjYSkJ5bJWl2in2OpJMl/VnSaklH9UcM6fUXpXX0SJqeW39Pbh0XS7pkB9vjm8D4FMcVjbZ/zjmka1tLOg64CjgrLWNcnfkBbkuvMxs03IibNWdY2sGvAmYD30jl50bE4cARwBcljQROAdZHxKERcQjwewBJ7yTryR0bEYcBW2ncKBwIXB0R7wL+TXZt6WaXcwiwuMGyDwJuioiJEfFEr+feDnwXeA8wAfgE2fWiLybrSbcdg6TDgc8A7wWOBs6TNLFBvFCwPYAZwJo0UvJlCrZ/bt27AAdExDqAiLgPeACYkpaxtmD9PcCRO4jRrKN8fMesOZtTY4WkY4CbJB1C1nBPTfOMIWtolgNXSroc+G2t50k2ZHs48EA20swwGt+idW1ELE3Ti4GxfVxOPU9ExMIG610OIGkFsCAiQtLyfozhOGBuRGxK67kVOJ7sJh9FirZHb0Xbv2YU2ZeAvIPIbiRSuxvfNcArZIcC5gBExFZJr0gaHhEvNJGj2YBzI27Wooj4q6RRwEfIbql4TES8JOluYNeIeCz1NCcDl0maHxGXkt3N6MaImNnkql7OTW8layhpcjkrgLMaPL+pyfVuy/2/jdf2Ge3GoILyV9l+hHDXgrjy22M7DbZ/zeb8ctPoyXMRsSUVnQn8MiJ+I+kWYE7utUOB/xTEbtZxHk43a5GkCcAQskbl2dSATyAbFkbSW4GXIuKnwJW8dsLUArLjrnun+faUtH8fQmhmOX8Ehko6Lxf3kZLe34f1DUQM9wBnSNot9XynAvcCG4G9JY2UNBQ4rYlYXgCG59ZRtP0BiIhngSGSag35OGB9bpbRZLc1hezLQm25I4Gnc429WencEzdrzjBJtaFcAdOAO4ELJD1MNhRbG55+N3CFpG1k91H+HEBEPCLpa8B8STul5y4Eeh+TbqiZ5aTh76nAdyTNIOs9rgOmt5z5AMQQEasl3QDcn2afHREPAUi6FFgErAVWNRHLM+nEux7gDrJ6+b/t38t8siH9O9M6RqXXn0923+vRZPeJznd0TgRu31E8Zp3k+4mb2etOOonuooj4VJ3ndge+T/al477aMfF03H5mRDza0WDNGnBP3MxedyLiIUl3SRqS+6147blNZGfO/086o/02N+A22LgnbmZmVlE+sc3MzKyi3IibmZlVlBtxMzOzinIjbmZmVlFuxM3MzCrKjbiZmVlFuRE3MzOrKDfiZmZmFfVfwcnvXFZu7BcAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import math\n", "\n", "def d_to_or(d):\n", " return math.exp( d*math.pi / math.sqrt(3) )\n", "\n", "#overall meta analysis estimate of D\n", "big_D = -0.184\n", "\n", "#So about a 15% decline\n", "meta_or = d_to_or(big_D)\n", "print( \"meta analysis D of {:1.3f} to odds ratio of {:1.3f}\".format(big_D, meta_or) )\n", "\n", "#plot over a range of 10 to 200\n", "baseline = np.arange(10,201)\n", "reduction = baseline - baseline*meta_or #This is the positive reduction in crimes\n", "\n", "fig, ax = plt.subplots(figsize=(8,6))\n", "ax.plot(baseline,reduction, color='k')\n", "ax.set_xlabel('Baseline Crime Counts ($t_0$)')\n", "ax.set_ylabel('Reduced Crimes ($\\Delta t$)')\n", "ax.grid(True,linestyle='--')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So that was alot of work just to figure out a percent reduction estimate!\n", "\n", "So what I do now is take the specific estimates from Braga's meta-analysis (Table 3), for the overall effects for violent and property crimes. Then I take the cost of law-enforcement estimates from this [Hunt paper](https://www.cambridge.org/core/journals/journal-of-benefit-cost-analysis/article/estimates-of-law-enforcement-costs-by-crime-type-for-benefitcost-analyses/0A1A55F70324FDBAA947FF1F18AA1B74). Those are costs directly related to law enforcement, like how much time it takes to investigate/respond to crime. Finally, I grab the biggest hot spot in my current working paper on [cost of crime hot spots](https://osf.io/preprints/socarxiv/nmq8r/), and get an ROI estimate conditional on the total number of crimes in the hot spot. (In particular it is hot spot 59, the largest in the [center of this map](https://apwheele.github.io/MathPosts/HotSpotMap.html).)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total estimated crime reduction ROI for hot spot 59: $357,289\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
CrimeCostTotalCrimeCrimeReductionROI
0Murder12435300.0000000.000000
1Assault829213035.356277293174.245531
2Robbery22293810.33491223036.518002
3Burglary1185294.0982854856.467525
4Theft102421229.95987530678.912323
5Veh Theft769517.2073285542.435613
\n", "
" ], "text/plain": [ " Crime Cost TotalCrime CrimeReduction ROI\n", "0 Murder 124353 0 0.000000 0.000000\n", "1 Assault 8292 130 35.356277 293174.245531\n", "2 Robbery 2229 38 10.334912 23036.518002\n", "3 Burglary 1185 29 4.098285 4856.467525\n", "4 Theft 1024 212 29.959875 30678.912323\n", "5 Veh Theft 769 51 7.207328 5542.435613" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Now lets go through an example for one of my hot spot areas\n", "#https://apwheele.github.io/MathPosts/HotSpotMap.html\n", "#Choosing the biggest area, ID 59 in the center of map\n", "\n", "#Taken from Braga's Table 3, specific estimates for violent/property\n", "or_v = d_to_or(-0.175)\n", "or_p = d_to_or(-0.084)\n", "\n", "#Hot spot 59 and cost of crime estimates for Texas from Hunt paper\n", "dat = [('v','Murder',124353,0),\n", " ('v','Assault',8292,130),\n", " ('v','Robbery',2229,38),\n", " ('p','Burglary',1185,29),\n", " ('p','Theft',1024,212),\n", " ('p','Veh Theft',769,51)]\n", "\n", "hs_59 = pd.DataFrame(dat, columns=['Type','Crime','Cost','TotalCrime'])\n", "hs_59['OR'] = or_v\n", "hs_59.loc[ hs_59['Type'] == 'p', 'OR'] = or_p\n", "hs_59['CrimeReduction'] = hs_59['TotalCrime'] - hs_59['TotalCrime']*hs_59['OR']\n", "hs_59['ROI'] = hs_59['CrimeReduction']*hs_59['Cost']\n", "\n", "print( \"Total estimated crime reduction ROI for hot spot 59: ${:,.0f}\".format( hs_59['ROI'].sum() ))\n", "\n", "var_order = ['Crime','Cost','TotalCrime','CrimeReduction','ROI']\n", "hs_59[var_order]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some caveats to this of course. One, you can always do a poor job on a hot spots intervention (so the lower bound of ROI is a negative number). Doing planning on the overall effect size though seems reasonable on its face to me. (So your banking on being as effective as the average hot spots experiment.) I wouldn't assume offhand you are going to be more effective than this, at least to justify an intervention from a cost-benefit perspective.\n", "\n", "This also averages together all different crimes, which may not be reasonable given the intervention. For example, one of the most successful hot spots experiments was the [Kansas City gun experiment](https://www.tandfonline.com/doi/abs/10.1080/07418829500096241). An intervention like that may not be reasonable to extrapolate to all crimes. So if the Dallas PD wanted to do something like that, it would probably only make sense to look at the relevant violent crimes, not all of the crimes (so only look at a specific row of the above table).\n", "\n", "These estimates are the 'credits' of conducting a hot spots policing strategy, not the 'debits' of officer time. So if what you want to do costs way more than the ROI in that area (like pay a ton of overtime to officers), it seems hot spots policing is potentially not worth the effort.\n", "\n", "These ROI estimates are also for 1.5 years of data, and the long term effectiveness of hot spots policing has not been established. So that means these estimates are not 'you do this hot spots thing for a short time period crime goes down forever', they are 'if you continue to do this hot spots thing, you should see reduced crimes around this margin'.\n", "\n", "So for this, if you wanted to create a POP officer for just this one hot spot, that ends up being a return of around $150k per year. So that seems to justify that single position. Smaller hot spots will need to be more along the lines of shifting current resources, as oppossed to creating new positions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "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.7.4" } }, "nbformat": 4, "nbformat_minor": 2 }