Simulated_assemblages.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparing real and simulated assemblages\n",
    "## Assessing typicality and deviation from expectations empirically\n",
    "\n",
    "\n",
    "With domestic assemblages, it is unclear in many cases whether the handful of sherds and other objects are meaningfully deposited rather than reflecting depositional and post-depositional processes. One way of considering this problem is to construct alternative assemblages via bootstrapping and other simulation techniques that let us glimpse larger universes of archaeological sites and compare them to our own.\n",
    "\n",
    "This code attempts to do so by taking the empirical frequency and density of artifacts separated by categories across multiple assemblages and generating statistically similar assemblages in order to estimate what \"typical\" assemblages with these properties would look like. These generated assmeblages can be contrasted with individual real assemblages or the set of all observed assemblages to quantify how close the observed archaeological record is to a random categorical distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A quick note about this code: this is a Jupyter Notebook. You may be running this in Binder or be new to Jupyter; if so, just click \"Run\" on the toolbar above to run each cell in turn, and then you can see the output step-by-step. Otherwise, go to Cell > Run All, which will run the entire the entire notebook.\n",
    "\n",
    "This code was written and runs in Python 3, and hasn't been tested in Python 2. It might work! I just haven't checked it yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the packages we need\n",
    "import numpy as np\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt\n",
    "import random as rd\n",
    "import pandas as pd\n",
    "#import the actual code, stored in a separate module.\n",
    "from simulate_assemblages import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Although not in the source code, let's visualize the results using some nice plots\n",
    "## this requires sklearn for the pca; you're welcome to not use this \n",
    "## and just switch the PCA over to something else\n",
    "import sklearn.preprocessing\n",
    "import sklearn.decomposition\n",
    "from sklearn.pipeline import Pipeline\n",
    "\n",
    "def plot_results(result):\n",
    "    \"\"\"Take the result object from either function and return a nice plot of the distribution.\"\"\"\n",
    "    #Create two plots: one a distance histogram, one a PCA scatter plot\n",
    "    fig, (ax1, ax2) = plt.subplots(1,2, figsize = (10,5))\n",
    "    fig.suptitle('Simulated assemblages (N = {0})'.format(result['N']))\n",
    "    \n",
    "    # First, the histogram of distances\n",
    "    ax1.hist(result['distance_simulation'], bins = 20)\n",
    "    ax1.plot([result['distance_observed']]*2, [0,result['N']/5], 'r-')\n",
    "    mean_dist = result['distance_simulation'].mean()\n",
    "    ax1.plot([mean_dist]*2, [0,result['N']/5], 'k-')\n",
    "    ax1.set_title('Distance from center (p = {0})'.format(result['p-value']))\n",
    "    ax1.set_xlabel('Distance from center of distribution')\n",
    "    \n",
    "    #PCA plot: create a 2-d plot of the (large) dimensional distribution array\n",
    "    collector_flat = result['raw_simulation'].reshape(result['N'], -1)\n",
    "    data_flat = result['raw_observed'].reshape(1, -1)\n",
    "    scaler = sklearn.preprocessing.StandardScaler().fit(collector_flat)\n",
    "    pca = sklearn.decomposition.PCA(n_components = 2)\n",
    "    pca.fit(scaler.transform(collector_flat))\n",
    "    workflow = Pipeline(steps = [('scale',scaler),('pca',pca)])\n",
    "    collector_transformed = workflow.transform(collector_flat)\n",
    "    observed_transformed = workflow.transform(data_flat)\n",
    "    mean_transformed = workflow.transform(collector_flat.mean(axis = 0).reshape(1,-1))\n",
    "    ax2.plot(collector_transformed[:,0], collector_transformed[:,1],'b.', alpha = .15)\n",
    "    ax2.plot(observed_transformed[:,0], observed_transformed[:,1],'r.', markersize = 10)\n",
    "    ax2.plot(mean_transformed[:,0], mean_transformed[:,1],'.', color = 'gray', markersize = 10)\n",
    "    ax2.set_xlabel('Normalized distance in projection')\n",
    "    ax2.set_title('2-d projection of {0}-d distribution ({1:.0f}% of variance)'.format(collector_flat.shape[1],round(100*sum(workflow['pca'].explained_variance_ratio_),0)))    #pca.explained_variance_ratio_ #note: this should just be used for visualization, NOT for classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example usage\n",
    "\n",
    "Our example data record the distribution of cooking pot fragments from Terrace House C at Monastiraki Katalimata, East Crete, Greece based on three main parts of the house: the East living quarters (\"East House\"), the West living quarters (\"West House\"), and the open-air area to the east (\"East Area\"). (Nowicki 2008)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = {'dish' : [1,0,2],\n",
    "     'tray' : [1,0,0],\n",
    "     'pot' : [3,3,1]}\n",
    "df = pd.DataFrame(data = dt, index = ['East House','West House','East Area'])\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Note that these are example data, and in practice you would likely want to separate vessels represented by single sherds vs. mendable or largely restorable vessels. The main takeaway is that this method works on any discrete count data.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = simulate_assemblages_collection(df, 10000, 2, 'calculated')\n",
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1000 simulated sites, consisting of three assemblages with three different types of vessels were generated using the empirical frequency and density of objects. Distances from the calculated center of the distribution were recorded as a measure of difference, and can be compared with the difference between the simulated assemblages and the observed assemblages to assess how close to an expected set of assemblages the distribution is. \n",
    "\n",
    "The output object, here stored as `r`, includes all generated assemblages, the observed assemblages, calculated distances for both the generated and observed assemblages, the p-value for measuring whether our observed data can be classified as unusually far from the center, and the parameters of the model. We can visualize these using the `plot_results` function defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_results(r)\n",
    "print('{0:2.1f}% of simulated sites were further from the center than the observed assemblages.'.format(100*r['p-value']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Overall, Monastiraki Katalimata is very close to typical, or what we would expect given the density and frequency of objects; it is typically within the closest 20-25% of simulated sites to the mean frequencies or expected distribution. This means that we shouldn't be \"surprised\" at the distribution of objects across the different assemblages: the frequency and density of artifacts explains the observed distribution pretty well. Whether the frequency of objects is itself meaningful would require comparison with another site.\n",
    "\n",
    "We can also check a single assemblage against the properties of all assemblages under consideration. This method addresses whether a single assemblage within a larger site may be less typical or may require an explanation greater than random chance. Let's look at the East Area, an outdoor cooking area of Terrace House C."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = simulate_assemblage_single(df, df.index.get_loc('East Area'), N = 10000, center = 'observed')\n",
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Despite the overall house having a close-to-expected distribution of different cooking vessels based on the density and frequency of different types, the East Area itself is somewhat less typical - only ~20% of all simulated assemblages are further from the expected distribution. While this may not be \"significant\" in a hypothesis testing framework, it is notable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_results(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The power of such an analysis is the ability to run it across scales and datasets, checking how typicality exists across and within sites."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, these same scores can be calculated for all assemblages in a site, enabling more \"typical\" or \"expected\" assemblages to be separated from \"atypical\" or \"surprising\" assemblages based on the relative variation between the frequency and density of different categories of objects between assemblages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create a series of the p-values in order\n",
    "pvs = []\n",
    "for i in df.index:\n",
    "    r = simulate_assemblage_single(df, df.index.get_loc(i), N = 10000, center = 'observed')\n",
    "    pvs.append(r)\n",
    "p = pd.Series([r['p-value'] for r in pvs], name = 'p-value', index = df.index)\n",
    "# concatenate the series and sort by p-value so the lowest (and therefore most exceptional) are at the top\n",
    "df_p = pd.concat((df,p), axis = 1).sort_values('p-value')\n",
    "df_p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These types of comparisons can easily be applied at scale to entire sites or regions, providing an intuitive, empirical tool for comparing expected and observed distributions of objects across assemblages."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bibliography \n",
    "\n",
    "Nowicki, Krzysztof. 2008. Monastiraki Katalimata: Excavation of a Cretan Refuge Site, 1993-2000. Prehistory Monographs 24. Philadelphia, PA: INSTAP Academic Press.\n"
   ]
  },
  {
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}