# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ## Measuring Neighbourhood Change with OLS

# Loading in the preliminaries

import numpy as np
import pandas as pd
import geopandas as gpd
import seaborn as sns
import shapely
import pysal
import matplotlib.pyplot as plt
import sklearn
# %matplotlib inline

# +
df = pd.read_csv("data/all_data_2016_CTs_V3.csv", dtype={'ctuid': str})
df["w"] = np.log(0.1 * df["h_n_hhld"] + 1) # weight field

df = df[(df["cma"] == 535)]

dfses = pd.read_csv("data/PCA_SES_Toronto.csv")

df = pd.merge(df,dfses,on = "ct_year", how = "left")

# +
sns.distplot(df["SESi"]);

print(df["SESi"].max())
# -

# This code computes rates of change of a variable in a neighbourhood (census tract) over time. e.g. average rate of change in population in a neighbourhood from 1991 to 2016

# +
# output setup
dfo = pd.DataFrame(df.ctuid.unique())
dfo.columns = ['ctuid']
dfo

# accessibility and SES models
variables = [
    'Ai_jobs_transit',
    'Ai_activities_transit',
    'Ai_jobs_auto',
    'Ai_activities_auto',
    'Ai_jobs_ratio',
    'Ai_activities_ratio',
    'Ai_combine',
    'SESi',
    'inc_avg_inc_16',
    'pop',
    "dw",
    
    "pop_elderly_prop",
    "pop_child_prop",
    "pop_immig",
    "pop_nodegree",
    "pop_unemployed",
    "inc_lico_p",
    "inc_total30p",
    "dw_majorrep",
    "dw_lone_prop",
    "dw_no_car",
]
for variable in variables:
    output = []
    out_row = {"ctuid": "x", variable + "_beta": 0}
    output.append(out_row)
    for ctuid in (df.ctuid.unique()):
        try:
            temp = df[[ 'ctuid','year',variable]]
            temp = temp.loc[(temp.ctuid == ctuid)]
            temp = temp.dropna()
            
            x = temp["year"].values.reshape(len(temp), 1)
            y = temp[variable].values.reshape(len(temp), 1)
            reg = sklearn.linear_model.LinearRegression()
            reg.fit(x,y)
            beta = reg.coef_[0][0]

            # output data
            out_row = {
                "ctuid": ctuid,
                variable + "_beta": beta
            }
            output.append(out_row)
        except:
            out_row = {
                "ctuid": ctuid,
                variable + "_beta": 0
            }
            output.append(out_row)
            
    output = pd.DataFrame(output)
    dfo = dfo.merge(output, on='ctuid', how = 'outer')

    
# travel behaviour
variables = [
    ["mean_nonworktrip_time_adult","w"],
    ["mean_commute_time_adult","w"],
    ["R_activities_per_day","w"],
]
for variable in variables:
    output = []
    out_row = {"ctuid": "x",variable[0] + "_beta": 0}
    output.append(out_row)
    for ctuid in (df.ctuid.unique()):
        try:
            temp = df[[ 'ctuid','year',variable[0],variable[1] ]]
            temp = temp.loc[(temp.ctuid == ctuid)]
            temp = temp.dropna()
            temp["X"] = temp[variable[0]]

            x = temp["year"].values.reshape(len(temp), 1)
            y = temp["X"].values.reshape(len(temp), 1)
            w = temp["w"].values.reshape(len(temp), 1)
            reg = sklearn.linear_model.LinearRegression()
            reg.fit(x,y, sample_weight = temp["w"])
            beta = reg.coef_[0][0]

            # output data
            out_row = {
                "ctuid": ctuid,
                variable[0] + "_beta": beta
            }
            output.append(out_row)
        except:
            out_row = {
                "ctuid": ctuid,
                variable[0] + "_beta": 0
            }
            output.append(out_row)
    output = pd.DataFrame(output)
    dfo = dfo.merge(output, on='ctuid', how = 'outer')
# -

dfo

# +
# write to csv (this is then used to link to boundaries in QGIS to make some maps!)
# dfo.to_csv("data/var_betas.csv")
# -

# Simple example of the final result

plt.scatter(x, y,  color='black')
plt.plot(x, reg.predict(x), color='blue', linewidth=3)
plt.xticks(())
plt.yticks(())
plt.show()

# Exmining scatter plots for relationship between change in one variable compared to change in another, at this neighbourhood level (each dot is a neighbourhood, values represent slope over time)

sns.scatterplot(x="Ai_combine_beta", y="mean_commute_time_adult_beta", data=dfo)

# Correlatiaon matrix of a few key variables

# +
import matplotlib.pyplot as plt

dfs = dfo[["Ai_combine_beta","dw_no_car_beta","mean_commute_time_adult_beta","R_activities_per_day_beta","SESi_beta"]]

plt.matshow(dfs.corr())
plt.show()
print(dfs.corr())


StarcoderdataJupyterScriptsDedupFiltered