<jupyter_start><jupyter_text>Spatial AnalysisVisualize statistical findings.  Setup<jupyter_code>%load_ext autoreload
%autoreload 2<jupyter_output><empty_output><jupyter_text>Imports<jupyter_code># standard 
import sys
import warnings

# data science
import pandas as pd

# geospatial
import geopandas as gpd

# visualization
import seaborn as sns

# module
sys.path.append("../src/")
import viz, stats
warnings.filterwarnings('ignore')<jupyter_output><empty_output><jupyter_text>Data Read<jupyter_code>df = gpd.read_file("../data/processed.csv",
                   GEOM_POSSIBLE_NAMES="geometry", KEEP_GEOM_COLUMNS="NO")
df.head(1)<jupyter_output><empty_output><jupyter_text>Spatial Autocorrelation AnalysisSpatial autocorrelation refers to the combination of two types of similarity: 1. spatial similarity  2. attribute similarityThese two similarities are combined to get an idea of how much a target feature changes over in a region with respect to itself (ie autocorrelation). The fundamental objective is to visualize and statistically quantify a relationship between the target feature and geography.<jupyter_code>target_feature = "median_hh_income" # using this feature for this example
df[target_feature] = df[target_feature].astype(float)<jupyter_output><empty_output><jupyter_text>Spatial Weights/Similarity First I will calculate spatial similarity by using the queen contiguity weight in particular, which checks if an adjacency relationship exists between two polygons (if a pair shares an edge or vertex).<jupyter_code>spatial_weights = stats.get_spatial_weights(df)<jupyter_output><empty_output><jupyter_text>Target Feature Similarity (Spatial Lag) Now that we have measured similarity with respect to the spatial region, we are ready to measure it with respect to the target feature. This is calculated (for each district $i$) as:$y_{\text{lag }i} = \sum_j w_{i, j} y_j$where $w_{i, j}$ is a matrix of the spatial weights (row-standardized if continuous) that we just calculated above.Mathematically, this is the weighted sum of neighboring-district target feature values.<jupyter_code>y = df[target_feature]
df["y_lag"] = stats.get_lagged_y(spatial_weights, y)<jupyter_output><empty_output><jupyter_text>Visualizing Baseline/Lagged Target Feature Values<jupyter_code>viz.lagged_feature(df, target_feature, df["y_lag"], "Median Household Income")<jupyter_output><empty_output><jupyter_text>We should expect to see a change in the range (due to calculations when lagging the raw values) and a shift in the colors of some districts. It is expected that any high values (red in color) will 'spread' to the neighboring districts (think of it as ironing them out).    Statistically Comparing the Target with the Spatially Lagged Target Above is a nice visual, but we still need to compare the two versions of the target (raw and spatially lagged) statistically to draw more reliable conclusions. To recap, we now have calculated the raw target values versus the weighted sum of the target with respect to each area's neighborhood (the lagged target values). To see if the geospatial aspect of this data affects the target variable distribution, we statistically compare these two versions of the target feature.  A Note on Population Size When the target variable is a rate based on populations of different sizes, the Moran's I values need to be adjusted (during both global and local calculations). This is done in my code.   Global Autocorrelation To start, we compare the target feature over the scope of the entire map, with the ambition of discovering any geographical clusters that exist.In this case, we have null hypothesis: - $H_0$: the data is completely spatially random (CSR)/no spatial autocorrelation. I will use Moran's I test to measure the correlation between a feature (raw target or lagged target) and the geographical space.<jupyter_code>morans_i_test, _ = stats.moran_test(y, spatial_weights)
round(morans_i_test.I, 3) # the test statistic<jupyter_output><empty_output><jupyter_text>This statistic $i \ \epsilon \ [-1, 1]$, where -1 is perfect dispersion (think of a chess board) and perfect separation (one half white, the other black, for example). A value of 0 is random.  SimulatingWe will take a computational approach (I found multiple sources describing that, due to the variability in geographical structures, analytical approaches to spatial autocorrelation are very unreliable). Here we assume the target feature to be distributed CSR. We then calculate a Moran's I statistic, and build up a distribution of statistics by repeating this process. Then we see where our actual calculated statistic falls relative to this distribution to calculated a significance value respective to the CSR hypothesis.<jupyter_code>viz.moran_test_simulation(morans_i_test, target_feature_name="Median Household Income")
# The simulated p value (recalling the distribution above is simulated). 
morans_i_test.p_sim<jupyter_output><empty_output><jupyter_text>If we reject the null hypothesis/with a significant p-value, we can conclude the target feature distribution is not spatially random, and that there is a signficant spatial association relative to Dodoma region districts.Recall this result applies to the map as a whole, as it is the global spatial autocorrelation.   Local AutocorrelationWhereas global autocorrelation can be interpreted as finding overall clustering of the target feature distribution relative to a region, the local autocorrelation looks to identify specific local clusters themselves. It is very similar to the above, only we perform a similar procedure over every individual geometric polygon (district) and get a significance measure for each. We can further put significant findings in context by turning to a local analysis where the attention shifts to detection of spatially correlated districts and spatial outliers.Mathematically, the local Moran's I is simply the product of the target feature value and its spatial lag at a particular location (district). As a result, the global Moran's I we calculated above is the average rate of change of the local statistics. As before, here I take the computational approach (as the analytical is very unreliable for a small sample size like we have), which looks to use conditional permutations for each location. The conditional part here means I hold out the target feature observation at a particular location $i$ (to avoid additional spatial dependency), and then randomly shuffle the remaining values before computing the local Moran value. This builds up a simulated reference distribution, which I do for each individual location.From here, we get determine if each district is locally significant in terms of the target feature distribution, relative to the other districts. To make this result more interpretable, I group findings according to the quadrant of the Moran's I scatterplot they come from; this is explained in much more detail below.To recap and put this all in context, the main objective here is to identify specific areas in the region that are driving overall patterns observed in the target feature.  Moran's i Scatterplot<jupyter_code>viz.morans_i_scatter(y, df["y_lag"], target_feature_name="Median Household Income")<jupyter_output><empty_output><jupyter_text>The interpretation of the above is as follows:We have essentially taken an array of local Moran's i statistics instead of a single global I like we did earlier (which would be the slope of the red regression line above).By using the mean dividing lines of each axis (the raw feature value and spatially lagged feature values), we divide the two dimensional space into four quadrants. Each quadrant has a characteristic we may be interested in characterizing a district by:- 1 (top right):    - there is a spatial association here (+, +)    - points represent districts with high average target feature value surrounded by districts that also have high average    - for simplicity, these will be referenced as 'diamond diamond'- 2 (top left):     - there is no spatial association here    - points represent districts with low average target feature value, but surrounded by districts with high average    - for simplicity, these will be referenced as 'rough in the diamonds'- 3 (bottom left):    - there is a spatial association here (-, -)    - points represent districts with low average target feature value surrounded by districts that also have a low average    - for simplicity, these will be referenced as 'rough rough'- 4 (bottom right):    - there is no spatial association here    - points represent districts with high average target feature value, but are surrounded by districts with low average    - for simplicity, these will be referenced as 'diamonds in the rough'    In technical terms, we can group these four quadrants into two supgroups:1. spatial clusters: quadrants one and three, where there is a spatial correlation (positive or negative)2. spatial outliers: quadrants two and four  Doing this for all Target Features Can calculate items over all features using the following code block:<jupyter_code># calculating moran's features for each feature of interest
df["population"] = df["population"].astype(float)
for feature in ["hs_diploma_rate", "vacancy_rate", "median_hh_income"]:
    df[feature] = df[feature].astype(float)
    df = stats.get_moran_features_for_col(df, feature, .05)<jupyter_output><empty_output><jupyter_text>Comparing Results with Custom DCI Metric<jupyter_code>quad_col = target_feature+"_moran_quadrant_descriptive"
sig_df = df[df[quad_col]!="No Significance"][["dci", quad_col]]
viz.dci_bar(sig_df, target_feature, "Median Household Income")<jupyter_output><empty_output><jupyter_text>Comparing Community Distress Metric with Moran's I Quadrant for Household Income<jupyter_code>viz.compare_dci_quadrants(df, target_feature, "Median Household Income")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.054201088051581706, 'fraction_numerical': 0.0018134193028410237, 'mean_word_length': 5.675857431069267, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '259061', 'license': 'MIT', 'n_tokens_mistral': 2588, 'n_tokens_neox': 2420, 'n_words': 1256, 'path': 'notebooks/spatial_analysis.ipynb', 'repo_name': 'sawezo/moran_dc', 'size': '9926'}