Using SD models to understand the differences between population measures at varying levels of geographic disagregation
=======================================================================================================================

In this recipe, we will use data at a national level to infer parameters
for a population aging model. We’ll then try two different ways of using
this trained model to understand variation between the behavior of each
of the states.

About this technique
--------------------

Firstly, we’ll use the parameters fit at the national level to predict
census data at the disaggregated level, and compare these predicted
state-level outputs with the measured values. This gives us a sense for
how different the populations of the states are from what we should
expect given our understanding of the national picture.

Secondly, we’ll fit parameters to a model at each of the state levels
actual measured census data, and then compare the differences in fit
parameters to each other and to the national expectation. This is a
helpful analysis if the parameter itself (and its inter-state variance)
is what we find interesting.

.. code:: ipython3

    %pylab inline
    import pandas as pd
    import pysd
    import scipy.optimize
    import geopandas as gpd


.. parsed-literal::

    Populating the interactive namespace from numpy and matplotlib


##Ingredients

Population data by age cohort
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

We start with data from the decennial census years 2000 and 2010, for
the male population by state and county. We have aggregated the data
into ten-year age cohorts (with the last cohort representing everyone
over 80 years old). The data collection task is described
`here <data/Census/US%20Census%20Data%20Collection.ipynb>`__. In this
analysis we will only use data at the state and national levels.

.. code:: ipython3

    data = pd.read_csv('../../data/Census/Males by decade and county.csv', header=[0,1], index_col=[0,1])
    data.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead tr th {
            text-align: left;
        }
    
        .dataframe thead tr:last-of-type th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr>
          <th></th>
          <th></th>
          <th colspan="9" halign="left">2000</th>
          <th colspan="9" halign="left">2010</th>
        </tr>
        <tr>
          <th></th>
          <th></th>
          <th>dec_1</th>
          <th>dec_2</th>
          <th>dec_3</th>
          <th>dec_4</th>
          <th>dec_5</th>
          <th>dec_6</th>
          <th>dec_7</th>
          <th>dec_8</th>
          <th>dec_9</th>
          <th>dec_1</th>
          <th>dec_2</th>
          <th>dec_3</th>
          <th>dec_4</th>
          <th>dec_5</th>
          <th>dec_6</th>
          <th>dec_7</th>
          <th>dec_8</th>
          <th>dec_9</th>
        </tr>
        <tr>
          <th>state</th>
          <th>county</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th rowspan="5" valign="top">1.0</th>
          <th>1.0</th>
          <td>3375.0</td>
          <td>3630.0</td>
          <td>2461.0</td>
          <td>3407.0</td>
          <td>3283.0</td>
          <td>2319.0</td>
          <td>1637.0</td>
          <td>825.0</td>
          <td>284.0</td>
          <td>3867.0</td>
          <td>4384.0</td>
          <td>3082.0</td>
          <td>3598.0</td>
          <td>4148.0</td>
          <td>3390.0</td>
          <td>2293.0</td>
          <td>1353.0</td>
          <td>454.0</td>
        </tr>
        <tr>
          <th>3.0</th>
          <td>9323.0</td>
          <td>10094.0</td>
          <td>7600.0</td>
          <td>9725.0</td>
          <td>10379.0</td>
          <td>8519.0</td>
          <td>6675.0</td>
          <td>4711.0</td>
          <td>1822.0</td>
          <td>11446.0</td>
          <td>12006.0</td>
          <td>9976.0</td>
          <td>11042.0</td>
          <td>12517.0</td>
          <td>12368.0</td>
          <td>10623.0</td>
          <td>6307.0</td>
          <td>2911.0</td>
        </tr>
        <tr>
          <th>5.0</th>
          <td>2002.0</td>
          <td>2198.0</td>
          <td>2412.0</td>
          <td>2465.0</td>
          <td>2178.0</td>
          <td>1699.0</td>
          <td>1026.0</td>
          <td>689.0</td>
          <td>301.0</td>
          <td>1673.0</td>
          <td>1739.0</td>
          <td>2260.0</td>
          <td>2208.0</td>
          <td>2233.0</td>
          <td>1910.0</td>
          <td>1490.0</td>
          <td>739.0</td>
          <td>324.0</td>
        </tr>
        <tr>
          <th>7.0</th>
          <td>1546.0</td>
          <td>1460.0</td>
          <td>1680.0</td>
          <td>1762.0</td>
          <td>1624.0</td>
          <td>1237.0</td>
          <td>774.0</td>
          <td>475.0</td>
          <td>187.0</td>
          <td>1471.0</td>
          <td>1577.0</td>
          <td>1798.0</td>
          <td>2016.0</td>
          <td>1928.0</td>
          <td>1581.0</td>
          <td>1140.0</td>
          <td>579.0</td>
          <td>211.0</td>
        </tr>
        <tr>
          <th>9.0</th>
          <td>3741.0</td>
          <td>3615.0</td>
          <td>3393.0</td>
          <td>3901.0</td>
          <td>3773.0</td>
          <td>3007.0</td>
          <td>2227.0</td>
          <td>1269.0</td>
          <td>550.0</td>
          <td>3741.0</td>
          <td>4252.0</td>
          <td>3312.0</td>
          <td>3719.0</td>
          <td>4129.0</td>
          <td>3782.0</td>
          <td>3052.0</td>
          <td>1723.0</td>
          <td>652.0</td>
        </tr>
      </tbody>
    </table>
    </div>



###A model of an aging population

The model we’ll use to represent the population is a simple aging chain,
with individuals aggregated into stocks by decade, to match the
agregation we used for the above data. Each cohort is promoted with a
timescale of 10 years, and there is some net inmigration, outmigration,
and death subsumed into the ``loss`` flow associated with each cohort.
This loss is controled by some yearly fraction that it will be our task
to understand.

.. image:: ../../../source/models/Aging_Chain/Aging_Chain.png
   :width: 600 px

.. code:: ipython3

    model = pysd.read_vensim('../../models/Aging_Chain/Aging_Chain.mdl')

Our model is initialy parameterized with 10 individuals in each stock,
no births, and a uniform loss rate of 5%. We’ll use data to set the
initial conditions, and infer the loss rates. Estimating births is
difficult, and so for this analysis, we’ll pay attention only to
individuals who have been born before the year 2000.

.. code:: ipython3

    model.run().plot();



.. image:: Exploring_models_across_geographic_scales_files/Exploring_models_across_geographic_scales_8_0.png


Geography Information
~~~~~~~~~~~~~~~~~~~~~

This information comes to us as a shape file ``.shp`` with its
associated ``.dbf`` and ``.shx`` conspirator files. Lets check the
plotting functionality:

.. code:: ipython3

    state_geo = gpd.read_file('../../data/Census/US_State.shp')
    state_geo.set_index('StateFIPSN', inplace=True)
    state_geo.plot();
    state_geo.head(2)




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>OBJECTID</th>
          <th>FIPS</th>
          <th>FIPSNum</th>
          <th>StateFIPS</th>
          <th>StateName</th>
          <th>CensusReg</th>
          <th>CensusDiv</th>
          <th>XCentroid</th>
          <th>YCentroid</th>
          <th>Notes</th>
          <th>geometry</th>
        </tr>
        <tr>
          <th>StateFIPSN</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>2</th>
          <td>1</td>
          <td>02000</td>
          <td>2000</td>
          <td>02</td>
          <td>Alaska</td>
          <td>West</td>
          <td>Pacific</td>
          <td>-1.882092e+06</td>
          <td>2.310348e+06</td>
          <td>None</td>
          <td>MULTIPOLYGON (((-2247528.774 2237995.012, -224...</td>
        </tr>
        <tr>
          <th>53</th>
          <td>9</td>
          <td>53000</td>
          <td>53000</td>
          <td>53</td>
          <td>Washington</td>
          <td>West</td>
          <td>Pacific</td>
          <td>-1.837353e+06</td>
          <td>1.340481e+06</td>
          <td>None</td>
          <td>MULTIPOLYGON (((-2124362.243 1480441.851, -212...</td>
        </tr>
      </tbody>
    </table>
    </div>




.. image:: Exploring_models_across_geographic_scales_files/Exploring_models_across_geographic_scales_10_1.png


Recipe Part A: Predict state-level values from national model fit
-----------------------------------------------------------------

Step 1: Initialize the model using census data
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

We can aggregate the county level data to the national scale by summing
across all geographies. This is relatively straightforward.

.. code:: ipython3

    country = data.sum()
    country




.. parsed-literal::

    2000  dec_1    20332536.0
          dec_2    20909490.0
          dec_3    19485544.0
          dec_4    21638975.0
          dec_5    21016627.0
          dec_6    15115009.0
          dec_7     9536197.0
          dec_8     6946906.0
          dec_9     3060483.0
    2010  dec_1    20703935.0
          dec_2    21878666.0
          dec_3    21645336.0
          dec_4    20033352.0
          dec_5    21597437.0
          dec_6    20451686.0
          dec_7    13926846.0
          dec_8     7424945.0
          dec_9     4083435.0
    dtype: float64



If we run the model using national data from the year 2000 as starting
conditions, we can see how the cohorts develop, given our arbitrary loss
rate values:

.. code:: ipython3

    model.run(return_timestamps=range(2000,2011), 
              initial_condition=(2000, country['2000'])).plot();



.. image:: Exploring_models_across_geographic_scales_files/Exploring_models_across_geographic_scales_14_0.png


Step 2: Fit the national level model to the remaining data
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

We’ve used half of our data (from the year 2000 census) to initialize
our model. Now we’ll use an optimization routine to choose the loss rate
parameters that best predict the census 2010 data. We’ll use the same
basic operations described in the previous recipe: `Fitting with
Optimization <2_1_Fitting_with_Optimization.ipynb>`__.

To make this simple, we’ll write a function that takes a list of
potential model parameters, and returns the model’s prediction in the
year 2010

.. code:: ipython3

    def exec_model(paramlist):
        params = dict(zip(['dec_%i_loss_rate'%i for i in range(1,10)], paramlist)) 
        output = model.run(initial_condition=(2000,country['2000']),
                           params=params, return_timestamps=2010)
        return output

Now we’ll define an error function that calls this executor and
calculates a sum of squared errors value. Remember that because we don’t
have birth information, we’ll only calculate error based upon age
cohorts 2 through 9.

.. code:: ipython3

    def error(paramlist):
        output = exec_model(paramlist)
        errors = output - country['2010']
        #don't tally errors in the first cohort, as we don't have info about births
        return sum(errors.values[0,1:]**2)

Now we can use the minimize function from scipy to find a vector of
parameters which brings the 2010 predictions into alignment with the
data.

.. code:: ipython3

    res = scipy.optimize.minimize(error, x0=[.05]*9,
                                  method='L-BFGS-B')
    country_level_fit_params = dict(zip(['dec_%i_loss_rate'%i for i in range(1,10)], res['x']))
    country_level_fit_params




.. parsed-literal::

    {'dec_1_loss_rate': 0.05,
     'dec_2_loss_rate': 0.05,
     'dec_3_loss_rate': 0.05,
     'dec_4_loss_rate': 0.05,
     'dec_5_loss_rate': 0.05,
     'dec_6_loss_rate': 0.05,
     'dec_7_loss_rate': 0.05,
     'dec_8_loss_rate': 0.05,
     'dec_9_loss_rate': 0.05}



If we run the national model forwards with these parameters, we see
generally good behavior, except for the 0-9yr demographic bracket, from
whom we expect less self-control. (And because we don’t have births
data.)

.. code:: ipython3

    model.run(params=country_level_fit_params,
              return_timestamps=range(2000,2011), 
              initial_condition=(2000, country['2000'])).plot();



.. image:: Exploring_models_across_geographic_scales_files/Exploring_models_across_geographic_scales_22_0.png


###Step 3: Make state-level predictions

If we want to look at the variances between the states and the national
level, we can try making state-level predictions using state-specific
initial conditions, but parameters fit at the national level.

.. code:: ipython3

    states = data.groupby(level=0).sum()
    states.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead tr th {
            text-align: left;
        }
    
        .dataframe thead tr:last-of-type th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr>
          <th></th>
          <th colspan="9" halign="left">2000</th>
          <th colspan="9" halign="left">2010</th>
        </tr>
        <tr>
          <th></th>
          <th>dec_1</th>
          <th>dec_2</th>
          <th>dec_3</th>
          <th>dec_4</th>
          <th>dec_5</th>
          <th>dec_6</th>
          <th>dec_7</th>
          <th>dec_8</th>
          <th>dec_9</th>
          <th>dec_1</th>
          <th>dec_2</th>
          <th>dec_3</th>
          <th>dec_4</th>
          <th>dec_5</th>
          <th>dec_6</th>
          <th>dec_7</th>
          <th>dec_8</th>
          <th>dec_9</th>
        </tr>
        <tr>
          <th>state</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>1.0</th>
          <td>312841.0</td>
          <td>329043.0</td>
          <td>301076.0</td>
          <td>315262.0</td>
          <td>321447.0</td>
          <td>246427.0</td>
          <td>165327.0</td>
          <td>109918.0</td>
          <td>45131.0</td>
          <td>312605.0</td>
          <td>338568.0</td>
          <td>321236.0</td>
          <td>297502.0</td>
          <td>321810.0</td>
          <td>318358.0</td>
          <td>229496.0</td>
          <td>124070.0</td>
          <td>56543.0</td>
        </tr>
        <tr>
          <th>2.0</th>
          <td>50687.0</td>
          <td>53992.0</td>
          <td>42537.0</td>
          <td>51442.0</td>
          <td>56047.0</td>
          <td>35804.0</td>
          <td>14974.0</td>
          <td>7628.0</td>
          <td>2325.0</td>
          <td>53034.0</td>
          <td>52278.0</td>
          <td>58166.0</td>
          <td>47753.0</td>
          <td>51856.0</td>
          <td>54170.0</td>
          <td>29869.0</td>
          <td>10392.0</td>
          <td>4151.0</td>
        </tr>
        <tr>
          <th>4.0</th>
          <td>395110.0</td>
          <td>384672.0</td>
          <td>386486.0</td>
          <td>391330.0</td>
          <td>352471.0</td>
          <td>257798.0</td>
          <td>187193.0</td>
          <td>144837.0</td>
          <td>61148.0</td>
          <td>463808.0</td>
          <td>466275.0</td>
          <td>455170.0</td>
          <td>422447.0</td>
          <td>418398.0</td>
          <td>381076.0</td>
          <td>300553.0</td>
          <td>178849.0</td>
          <td>89247.0</td>
        </tr>
        <tr>
          <th>5.0</th>
          <td>188589.0</td>
          <td>201405.0</td>
          <td>180801.0</td>
          <td>187516.0</td>
          <td>186931.0</td>
          <td>149142.0</td>
          <td>104621.0</td>
          <td>72629.0</td>
          <td>33050.0</td>
          <td>201821.0</td>
          <td>205074.0</td>
          <td>196956.0</td>
          <td>183761.0</td>
          <td>192596.0</td>
          <td>188081.0</td>
          <td>143285.0</td>
          <td>81138.0</td>
          <td>38925.0</td>
        </tr>
        <tr>
          <th>6.0</th>
          <td>2669364.0</td>
          <td>2588761.0</td>
          <td>2556975.0</td>
          <td>2812648.0</td>
          <td>2495158.0</td>
          <td>1692007.0</td>
          <td>1002881.0</td>
          <td>725610.0</td>
          <td>331342.0</td>
          <td>2573619.0</td>
          <td>2780997.0</td>
          <td>2849483.0</td>
          <td>2595717.0</td>
          <td>2655307.0</td>
          <td>2336519.0</td>
          <td>1489395.0</td>
          <td>780576.0</td>
          <td>456217.0</td>
        </tr>
      </tbody>
    </table>
    </div>



We can now generate a prediction by setting the model’s intitial
conditions with state level data, and parameters fit in the national
case. I’ve created a ``model_runner`` helper function to make the code
easier to read, but this could be conducted in a single line if we so
chose.

.. code:: ipython3

    def model_runner(row):
        result = model.run(params=country_level_fit_params, 
                           initial_condition=(2000, row['2000']), 
                           return_timestamps=2010,
                           return_columns=[
                               "dec_1", "dec_2", "dec_3", 
                               "dec_4", "dec_5", "dec_6",
                               "dec_7", "dec_8", "dec_9"
                           ])
        return result.loc[2010]
        
    state_predictions = states.apply(model_runner, axis=1)
    state_predictions.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>dec_1</th>
          <th>dec_2</th>
          <th>dec_3</th>
          <th>dec_4</th>
          <th>dec_5</th>
          <th>dec_6</th>
          <th>dec_7</th>
          <th>dec_8</th>
          <th>dec_9</th>
        </tr>
        <tr>
          <th>state</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>1.0</th>
          <td>68817.239206</td>
          <td>1.425135e+05</td>
          <td>1.752830e+05</td>
          <td>1.856499e+05</td>
          <td>1.905074e+05</td>
          <td>1.766502e+05</td>
          <td>1.430809e+05</td>
          <td>104576.864395</td>
          <td>112780.892489</td>
        </tr>
        <tr>
          <th>2.0</th>
          <td>11149.879343</td>
          <td>2.323983e+04</td>
          <td>2.717858e+04</td>
          <td>2.883606e+04</td>
          <td>3.114106e+04</td>
          <td>2.841718e+04</td>
          <td>2.006445e+04</td>
          <td>11731.576605</td>
          <td>9263.036022</td>
        </tr>
        <tr>
          <th>4.0</th>
          <td>86914.373061</td>
          <td>1.731934e+05</td>
          <td>2.158222e+05</td>
          <td>2.308789e+05</td>
          <td>2.268514e+05</td>
          <td>1.985338e+05</td>
          <td>1.576861e+05</td>
          <td>120469.144700</td>
          <td>140567.692416</td>
        </tr>
        <tr>
          <th>5.0</th>
          <td>41484.889527</td>
          <td>8.658169e+04</td>
          <td>1.061958e+05</td>
          <td>1.115457e+05</td>
          <td>1.128048e+05</td>
          <td>1.048007e+05</td>
          <td>8.660780e+04</td>
          <td>65340.896065</td>
          <td>74727.981971</td>
        </tr>
        <tr>
          <th>6.0</th>
          <td>587193.689179</td>
          <td>1.167877e+06</td>
          <td>1.443928e+06</td>
          <td>1.583683e+06</td>
          <td>1.589019e+06</td>
          <td>1.372822e+06</td>
          <td>1.015240e+06</td>
          <td>699665.480214</td>
          <td>759574.642778</td>
        </tr>
      </tbody>
    </table>
    </div>



Step 4: Compare model predictions with measured data
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Comparing the state level predictions with the actual data, we can see
where our model most under or overpredicts population for each
region/cohort combination.

.. code:: ipython3

    diff = state_predictions-states['2010']
    diff.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>dec_1</th>
          <th>dec_2</th>
          <th>dec_3</th>
          <th>dec_4</th>
          <th>dec_5</th>
          <th>dec_6</th>
          <th>dec_7</th>
          <th>dec_8</th>
          <th>dec_9</th>
        </tr>
        <tr>
          <th>state</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>1.0</th>
          <td>-2.437878e+05</td>
          <td>-1.960545e+05</td>
          <td>-1.459530e+05</td>
          <td>-1.118521e+05</td>
          <td>-1.313026e+05</td>
          <td>-141707.812835</td>
          <td>-86415.080653</td>
          <td>-19493.135605</td>
          <td>56237.892489</td>
        </tr>
        <tr>
          <th>2.0</th>
          <td>-4.188412e+04</td>
          <td>-2.903817e+04</td>
          <td>-3.098742e+04</td>
          <td>-1.891694e+04</td>
          <td>-2.071494e+04</td>
          <td>-25752.819389</td>
          <td>-9804.550877</td>
          <td>1339.576605</td>
          <td>5112.036022</td>
        </tr>
        <tr>
          <th>4.0</th>
          <td>-3.768936e+05</td>
          <td>-2.930816e+05</td>
          <td>-2.393478e+05</td>
          <td>-1.915681e+05</td>
          <td>-1.915466e+05</td>
          <td>-182542.183715</td>
          <td>-142866.899500</td>
          <td>-58379.855300</td>
          <td>51320.692416</td>
        </tr>
        <tr>
          <th>5.0</th>
          <td>-1.603361e+05</td>
          <td>-1.184923e+05</td>
          <td>-9.076017e+04</td>
          <td>-7.221535e+04</td>
          <td>-7.979121e+04</td>
          <td>-83280.336688</td>
          <td>-56677.198680</td>
          <td>-15797.103935</td>
          <td>35802.981971</td>
        </tr>
        <tr>
          <th>6.0</th>
          <td>-1.986425e+06</td>
          <td>-1.613120e+06</td>
          <td>-1.405555e+06</td>
          <td>-1.012034e+06</td>
          <td>-1.066288e+06</td>
          <td>-963697.077944</td>
          <td>-474154.512342</td>
          <td>-80910.519786</td>
          <td>303357.642778</td>
        </tr>
      </tbody>
    </table>
    </div>



This is a little easier to understand if we weight it by the actual
measured population:

.. code:: ipython3

    diff_percent = (state_predictions-states['2010'])/states['2010']
    diff_percent.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>dec_1</th>
          <th>dec_2</th>
          <th>dec_3</th>
          <th>dec_4</th>
          <th>dec_5</th>
          <th>dec_6</th>
          <th>dec_7</th>
          <th>dec_8</th>
          <th>dec_9</th>
        </tr>
        <tr>
          <th>state</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>1.0</th>
          <td>-0.779859</td>
          <td>-0.579070</td>
          <td>-0.454348</td>
          <td>-0.375971</td>
          <td>-0.408013</td>
          <td>-0.445121</td>
          <td>-0.376543</td>
          <td>-0.157114</td>
          <td>0.994604</td>
        </tr>
        <tr>
          <th>2.0</th>
          <td>-0.789760</td>
          <td>-0.555457</td>
          <td>-0.532741</td>
          <td>-0.396141</td>
          <td>-0.399470</td>
          <td>-0.475407</td>
          <td>-0.328252</td>
          <td>0.128905</td>
          <td>1.231519</td>
        </tr>
        <tr>
          <th>4.0</th>
          <td>-0.812607</td>
          <td>-0.628559</td>
          <td>-0.525843</td>
          <td>-0.453473</td>
          <td>-0.457810</td>
          <td>-0.479018</td>
          <td>-0.475347</td>
          <td>-0.326420</td>
          <td>0.575041</td>
        </tr>
        <tr>
          <th>5.0</th>
          <td>-0.794447</td>
          <td>-0.577803</td>
          <td>-0.460814</td>
          <td>-0.392985</td>
          <td>-0.414293</td>
          <td>-0.442790</td>
          <td>-0.395556</td>
          <td>-0.194694</td>
          <td>0.919794</td>
        </tr>
        <tr>
          <th>6.0</th>
          <td>-0.771841</td>
          <td>-0.580051</td>
          <td>-0.493267</td>
          <td>-0.389886</td>
          <td>-0.401568</td>
          <td>-0.412450</td>
          <td>-0.318354</td>
          <td>-0.103655</td>
          <td>0.664942</td>
        </tr>
      </tbody>
    </table>
    </div>



###Step 5: Merge with geo data to plot

I’m using geopandas to manage the shapefiles, and it has its own
plotting functionality. Unfortunately, it is not a particularly well
defined functionality.

.. code:: ipython3

    geo_diff = state_geo.join(diff_percent)
    geo_diff.plot(column='dec_4')
    geo_diff.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>OBJECTID</th>
          <th>FIPS</th>
          <th>FIPSNum</th>
          <th>StateFIPS</th>
          <th>StateName</th>
          <th>CensusReg</th>
          <th>CensusDiv</th>
          <th>XCentroid</th>
          <th>YCentroid</th>
          <th>Notes</th>
          <th>geometry</th>
          <th>dec_1</th>
          <th>dec_2</th>
          <th>dec_3</th>
          <th>dec_4</th>
          <th>dec_5</th>
          <th>dec_6</th>
          <th>dec_7</th>
          <th>dec_8</th>
          <th>dec_9</th>
        </tr>
        <tr>
          <th>StateFIPSN</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>2</th>
          <td>1</td>
          <td>02000</td>
          <td>2000</td>
          <td>02</td>
          <td>Alaska</td>
          <td>West</td>
          <td>Pacific</td>
          <td>-1.882092e+06</td>
          <td>2.310348e+06</td>
          <td>None</td>
          <td>MULTIPOLYGON (((-2247528.774 2237995.012, -224...</td>
          <td>-0.789760</td>
          <td>-0.555457</td>
          <td>-0.532741</td>
          <td>-0.396141</td>
          <td>-0.399470</td>
          <td>-0.475407</td>
          <td>-0.328252</td>
          <td>0.128905</td>
          <td>1.231519</td>
        </tr>
        <tr>
          <th>53</th>
          <td>9</td>
          <td>53000</td>
          <td>53000</td>
          <td>53</td>
          <td>Washington</td>
          <td>West</td>
          <td>Pacific</td>
          <td>-1.837353e+06</td>
          <td>1.340481e+06</td>
          <td>None</td>
          <td>MULTIPOLYGON (((-2124362.243 1480441.851, -212...</td>
          <td>-0.792005</td>
          <td>-0.585096</td>
          <td>-0.510514</td>
          <td>-0.430876</td>
          <td>-0.426070</td>
          <td>-0.459076</td>
          <td>-0.393843</td>
          <td>-0.135940</td>
          <td>0.674966</td>
        </tr>
        <tr>
          <th>23</th>
          <td>10</td>
          <td>23000</td>
          <td>23000</td>
          <td>23</td>
          <td>Maine</td>
          <td>Northeast</td>
          <td>New England</td>
          <td>2.068850e+06</td>
          <td>1.172787e+06</td>
          <td>None</td>
          <td>MULTIPOLYGON (((1951177.135 1127914.539, 19511...</td>
          <td>-0.765509</td>
          <td>-0.556475</td>
          <td>-0.409031</td>
          <td>-0.338041</td>
          <td>-0.430710</td>
          <td>-0.473652</td>
          <td>-0.409398</td>
          <td>-0.164981</td>
          <td>0.726818</td>
        </tr>
        <tr>
          <th>27</th>
          <td>11</td>
          <td>27000</td>
          <td>27000</td>
          <td>27</td>
          <td>Minnesota</td>
          <td>Midwest</td>
          <td>West North Central</td>
          <td>1.310476e+05</td>
          <td>9.821300e+05</td>
          <td>None</td>
          <td>POLYGON ((-91052.168 1282100.079, -30035.391 1...</td>
          <td>-0.786822</td>
          <td>-0.557537</td>
          <td>-0.463915</td>
          <td>-0.369641</td>
          <td>-0.403840</td>
          <td>-0.452106</td>
          <td>-0.335580</td>
          <td>-0.123440</td>
          <td>0.645598</td>
        </tr>
        <tr>
          <th>26</th>
          <td>18</td>
          <td>26000</td>
          <td>26000</td>
          <td>26</td>
          <td>Michigan</td>
          <td>Midwest</td>
          <td>East North Central</td>
          <td>8.425679e+05</td>
          <td>8.094373e+05</td>
          <td>None</td>
          <td>MULTIPOLYGON (((764918.731 786184.882, 772386....</td>
          <td>-0.746814</td>
          <td>-0.547370</td>
          <td>-0.378482</td>
          <td>-0.282202</td>
          <td>-0.364176</td>
          <td>-0.423316</td>
          <td>-0.329958</td>
          <td>-0.074694</td>
          <td>0.747247</td>
        </tr>
      </tbody>
    </table>
    </div>




.. image:: Exploring_models_across_geographic_scales_files/Exploring_models_across_geographic_scales_33_1.png


Recipe Part B: fit state-by-state models
----------------------------------------

Now lets try optimizing the model’s parameters specifically to each
state, and comparing with the national picture.

Step 1: Write the optimization functions to account for the state
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

We’ll start as before with functions that run the model and compute the
error (this time with a parameter for the information about the state in
question) and add a function to optimize and return the best fit
parameters for each state.

.. code:: ipython3

    def exec_model(paramlist, state):
        params = dict(zip(['dec_%i_loss_rate'%i for i in range(1,10)], paramlist)) 
        output = model.run(initial_condition=(2000,state['2000']),
                           params=params, return_timestamps=2010).loc[2010]
        return output
    
    def error(paramlist, state):
        output = exec_model(paramlist, state)
        errors = output - state['2010']
        #don't tally errors in the first cohort, as we don't have info about births
        sse = sum(errors.values[1:]**2)
        return sse

###Step 2: Apply optimization to each state We can wrap the optimizer in
a function that takes census information about each state and returns an
optimized set of parameters for that state. If we apply it to the states
dataframe, we can get out a similar dataframe that includes optimized
parameters.

.. code:: ipython3

    %%capture 
    def optimize_params(row):
        res = scipy.optimize.minimize(lambda x: error(x, row),
                                      x0=[.05]*9,
                                      method='L-BFGS-B');
        return pd.Series(index=['dec_%i_loss_rate'%i for i in range(1,10)], data=res['x'])
        
    state_fit_params = states.apply(optimize_params, axis=1)
    state_fit_params.head()

Step 3: Merge with geographic data
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

As we’re looking at model parameters which themselves are multiplied by
populations to generate actual flows of people, we can look at the
difference between parameters directly without needing to normalize.

.. code:: ipython3

    geo_diff = state_geo.join(state_fit_params)
    geo_diff.plot(column='dec_4_loss_rate')
    geo_diff.head(3)




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>OBJECTID</th>
          <th>FIPS</th>
          <th>FIPSNum</th>
          <th>StateFIPS</th>
          <th>StateName</th>
          <th>CensusReg</th>
          <th>CensusDiv</th>
          <th>XCentroid</th>
          <th>YCentroid</th>
          <th>Notes</th>
          <th>geometry</th>
          <th>dec_1_loss_rate</th>
          <th>dec_2_loss_rate</th>
          <th>dec_3_loss_rate</th>
          <th>dec_4_loss_rate</th>
          <th>dec_5_loss_rate</th>
          <th>dec_6_loss_rate</th>
          <th>dec_7_loss_rate</th>
          <th>dec_8_loss_rate</th>
          <th>dec_9_loss_rate</th>
        </tr>
        <tr>
          <th>StateFIPSN</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>2</th>
          <td>1</td>
          <td>02000</td>
          <td>2000</td>
          <td>02</td>
          <td>Alaska</td>
          <td>West</td>
          <td>Pacific</td>
          <td>-1.882092e+06</td>
          <td>2.310348e+06</td>
          <td>None</td>
          <td>MULTIPOLYGON (((-2247528.774 2237995.012, -224...</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
        </tr>
        <tr>
          <th>53</th>
          <td>9</td>
          <td>53000</td>
          <td>53000</td>
          <td>53</td>
          <td>Washington</td>
          <td>West</td>
          <td>Pacific</td>
          <td>-1.837353e+06</td>
          <td>1.340481e+06</td>
          <td>None</td>
          <td>MULTIPOLYGON (((-2124362.243 1480441.851, -212...</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
        </tr>
        <tr>
          <th>23</th>
          <td>10</td>
          <td>23000</td>
          <td>23000</td>
          <td>23</td>
          <td>Maine</td>
          <td>Northeast</td>
          <td>New England</td>
          <td>2.068850e+06</td>
          <td>1.172787e+06</td>
          <td>None</td>
          <td>MULTIPOLYGON (((1951177.135 1127914.539, 19511...</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
          <td>0.05</td>
        </tr>
      </tbody>
    </table>
    </div>




.. image:: Exploring_models_across_geographic_scales_files/Exploring_models_across_geographic_scales_39_1.png