5SSG2060_Practical_9_Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5SSG2060 Practical Week 10 :  Multivariate Regression & Spatial Regression\n",
    "<a href=\"#This Week's Overview\">This Week's Overview</a>\n",
    " \n",
    " <a href=\"#Learn Outcomes\">Learn Outcomes</a> \n",
    " \n",
    " <a href='#Get organized'>Get organized</a>\n",
    "- <a href='#Import libraries'>Import libraries</a>\n",
    "- <a href='#Prepare data'>Prepare data</a>\n",
    "\n",
    "<a href='#Aspatial Multivariate Regression'>Aspatial Multivariate Regression</a>\n",
    "\n",
    "\n",
    "<a href='#Spatial Regression'>Spatial Regression</a>\n",
    "  - <a href='#Spatial Lag model'>Spatial Lag model<a/>\n",
    "  - <a href='#Spatial Error model'>Spatial Error model<a/>\n",
    "  - <a href='#Prediction performance of spatial models'>Prediction performance of spatial models</a>\n",
    "  - <a href='#GWR Prediction'>GWR Prediction<a/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id=\"This Week's Overview\">This Week's Overview</a>\n",
    "This week, we will mainly use `PySAL` and `Numpy` to realize both nonspatial and spatial regression purposes. The spatial regression models will be introduced with both Spatial Lag model and Spatial Error model for your comparison, and GWR (geographical weighted regression) prediction will also follow as the end of this practical. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id=\"Learn Outcomes\">Learn Outcomes</a>\n",
    "By the end of this practical, you would be able to conduct and interpret\n",
    "- Aspatial Multivariate Regression Model\n",
    "- Spatial Lag Model\n",
    "- Spatial Error Model\n",
    "\n",
    "And be able to do corresponding prediction using geographical weighted regression model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id='Get organized'>Get organized</a>\n",
    "## <a id='Import libraries'>Import libraries</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import urllib\n",
    "import zipfile\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import pysal as ps\n",
    "from pysal.spreg import ols\n",
    "from pysal.spreg import ml_error\n",
    "from pysal.spreg import ml_lag\n",
    "\n",
    "import seaborn as sns\n",
    "sns.set(style=\"whitegrid\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id='Prepare data'>Prepare data</a>\n",
    "Please copy the **shapefile data \"borough_airbnb_housing\"** we produced last week into your \"data\" folder for this week."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in your data and get the headers presented\n",
    "gdf=gpd.read_file('data/borough_airbnb_housing.shp')\n",
    "gdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the labels changed as expected (10 letters limit, remember?), we need to rename them back into some human language again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#change column names here. \n",
    "cnames = {'Borough' : 'borough', \n",
    "          'Sector' : 'sector', \n",
    "          'Population': 'population', \n",
    "          'Households' : 'households', \n",
    "          'Dwellings' : 'dwellings', \n",
    "          'person_dwe' : 'per_dwell', \n",
    "          'house_stoc' : 'house_growth', \n",
    "          'house_st_1' : 'house_stock',\n",
    "          'net_homes' : 'net_home', \n",
    "          'private_re': 'private_rent_price',\n",
    "          'house_pric': 'house_price',\n",
    "          'ratio_pric': 'ratio_price_earn',\n",
    "          'price': 'listing_price',\n",
    "          'number_of_': 'reviews',\n",
    "          'Number of': 'airbnb_amount',\n",
    "          'geometry': 'geometry' \n",
    "         }\n",
    "gdf.rename(columns=cnames, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id='Aspatial Multivariate Regression'>Aspatial Multivariate Regression<a/>\n",
    "We will use `PySAL` to analyse the relations between `Airbnb listing price` and other housing and earning variables by borough in London. Before introducing explicitly spatial methods next week, we will run a simple ordinary least squares (OLS) linear regression model to analyze the relationship between these variables. Hence providing us the capability to interpret regression results on which the spatial models will build on; and enabling us to evaluate the meaningfulness of spatial regression models upon comparions. \n",
    "\n",
    "Normally, an OLS linear regression is to explain a given variable $y$, as a linear function of a set of other variables $X$:\n",
    "\n",
    "$y_i$=α+β$X_i$+$ϵ_i$\n",
    "\n",
    "where $i$ is the boroughs in our example, and $\\epsilon_i$ is the error term for our probabilistic model. \n",
    "\n",
    "Now let's start with read in the dependent variable $y$ and independent variables $X$:\n",
    "\n",
    "We will start directly with more than 3 independant variables for today, but if in your own project, the X is a 2 dimensional array, which means only 2 independant variables, then the ordinary least squares (OLS) regression is performed with [pysal.spreg](https://pysal.readthedocs.io/en/v1.11.0/library/spreg/index.html), with outputs callable through its `summary` feature and I leave more to explore by yourself. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read the .dbf file from your shapefile data\n",
    "f = ps.open('data/borough_airbnb_housing.dbf','r')\n",
    "# Read in the listing_price (dependent variable) into an array y\n",
    "y = np.array(f.by_col['price'])\n",
    "y.shape = (len(y),1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# value for independent variables into a one dimmensional array X. \n",
    "# You can feel free to change the independant variables\n",
    "X= []\n",
    "X.append(f.by_col['number_of_']) # number_of_reviews\n",
    "X.append(f.by_col['private_re']) # private_rent_price\n",
    "X.append(f.by_col['ratio_pric']) # ratio_price_earning\n",
    "X.append(f.by_col['person_dwe']) # person_dwelling\n",
    "X = np.array(X).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m1 = ols.OLS(y, X, name_y = 'price', name_x = ['number_of_reviews', 'private_rent','ratio_price_earning', 'person_dwelling'], name_ds = 'borough_airbnb_housing')\n",
    "print(m1.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 1\">Task 1<a/>\n",
    "Interpret the regression summary, and get your answers to:\n",
    "- Is current OLS model well enough? What is the adjusted R-squared value and how to interpret it?\n",
    "- Which variable(s) have significant influence on airbnb listing price, and how they are going to affect the pricing?\n",
    "- Which are the parameters telling us the model fit?\n",
    "- Is multicollinearity existing? If so, do you have better solutions? If you are going to replace some variables for improvement, which variables will be considered?\n",
    "\n",
    "However, you may already found this method is not good enough! We call the array by using pysal wrapper $f$, so we can only define **X** with the original cutted labels for variables, which is not easy to interpret. How if we want to keep using the geodataframe to do the regression? We can use `statsmodels` (you've seen it last week in optional section) to realize the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "x = ['airbnb_amount', 'private_rent_price', 'ratio_price_earn', 'per_dwell'] # define input/independant variables\n",
    "X = gdf.loc[:, x].dropna()\n",
    "y = np.array(gdf['listing_price']) # define output/dependant variable\n",
    "# Note the difference in argument order\n",
    "model = sm.OLS(y, X).fit()\n",
    "predictions = model.predict(X) # make the predictions by the model\n",
    "\n",
    "# Print out the statistics\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What have you got this time? The output summary is very different from the previous one.\n",
    "\n",
    "\"Ordinary Least Squares (OLS)\" method means that we’re trying to fit a regression line that would minimize the square of distance from the regression line. Date and Time are the precise moment you got this summary, Number of observations is the count of rows, and in our dataset it equals to the number of boroughs. $Df$ of residuals and models relates to the degrees of freedom (\"number of values in the final calculation of a statistic that are free to vary\").\n",
    "\n",
    "Take the private_rent_price variable as an example, the coefficient of 0.1103 means that as the private rent price increases by 1, the predicted value of airbnb listing price increases by 11%. Similar explanations could be made on other coefficients with a signifant $p$ value.\n",
    "\n",
    "Other important measures include: \n",
    "\n",
    "- R-squared: how much percentage of variance the model can explain; \n",
    "Std err (standard error): standard deviation of the sampling distribution of a statistic, e.g. mean; \n",
    "\n",
    "- t scores and p-values: hypothesis test. If private_rent_price variable has statistically significant p-value; there is a 95% confidence intervals for it.\n",
    "\n",
    "If we want to add a constant intercept to the model, we can call the following:\n",
    "### <a id=\"Task 2\">Task 2<a/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = sm.add_constant(X) ## let's add an intercept (beta_0) to our model\n",
    "\n",
    "# Note the difference in argument order\n",
    "model_c = ??? ## sm.OLS(output, input)\n",
    "predictions = model_c.predict(X)\n",
    "\n",
    "# Print out the statistics\n",
    "model_c.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the constant term the coefficients are different now. The previous model without constant term, we forced our model to go through the origin data, but now we have a y-intercept at 47.1647, the slope of the private rent price variable also decreased from 0.11 to 0.09.\n",
    "\n",
    "The model explanation capability decreased from 98.1% of the variance in dependent variable, airbnb listing price, to 84.8%, reading from the R-squared values of summaries. Think of another scenario, how if we add other variables to this regression model, will the R² get higher? Discuss with your neighbours if you want to share your opinion.\n",
    "\n",
    "We can see both \"private rent price\" and \"average persons living in each dwelling\" are statistically significant in predicting (or estimating) the airbnb listing price. As the private rent price increases by 1, the airbnb listing price will increase by 0.09, which is almost 10 times relationship; when number of people living in each dwelling increases by 1 (people), the airbnb listing price will drop by 26.45 (sterlings). However, with an added constant, the number of airbnbs in borough begun to be insignificant in predicting the price.\n",
    "\n",
    "In all, `Statsmodels` is functionally sufficient to conduct multivariate regression, and you can add up to 13 variables for your research！In addition, it also provide follow-on methods to diagnose your regression result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Regression diagnostics'>Regression diagnostics<a/>\n",
    "\n",
    "The rationale is to check the normality of the residuals from our model. The meanings of statistics could be googled or discuss with your neighbors in class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.compat import lzip\n",
    "import statsmodels.stats.api as sms\n",
    "# Omni test\n",
    "name = ['Chi^2', 'Two-tail probability']\n",
    "test = sms.omni_normtest(model.resid)\n",
    "lzip(name, test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 3\">Task 3<a/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Jarque-Bera test\n",
    "name = ['Jarque-Bera', 'Chi^2 two-tail prob.', 'Skew', 'Kurtosis']\n",
    "test = sms.jarque_bera(?.resid)\n",
    "lzip(name, test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multicollinearity on condition number\n",
    "np.linalg.cond(model.model.exog)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Heteroskedasticity tests\n",
    "# Breush-Pagan test\n",
    "name = ['Lagrange multiplier statistic', 'p-value',\n",
    "        'f-value', 'f p-value']\n",
    "test = sms.het_breuschpagan(model.resid, model.model.exog)\n",
    "lzip(name, ?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Goldfeld-Quandt test\n",
    "name = ['F statistic', 'p-value']\n",
    "test = sms.het_goldfeldquandt(model.resid, ?)\n",
    "lzip(?, test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Partial Regression Plots'>Partial Regression Plots<a/>\n",
    "Partial regression plots show us the relationships between each independent variable and dependent varaiable individually, with a linear fit line added. However, you may find partical regression take natural log transformed data as default variables, the axis is with $ e(X) $ specification, and the scales for axis had changed even with negative values. But for our analysis, we didn't do data transformation beforehand for consistence, hence the axis have changed into Xth power of e. It doesn't affect the result for this practical, we only take it as an additional supportive part, helping you better interpret your regression summary intuitively. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12,8))\n",
    "fig = sm.graphics.plot_partregress_grid(model, fig=fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "fig = sm.graphics.plot_partregress('listing_price', 'private_rent_price', ['airbnb_amount', 'ratio_price_earn', 'per_dwell'],  ax=ax, data=gdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Influence Plots'>Influence Plots<a/>\n",
    "Influence plots could be used to signify the influences on regression result from each row (borough for our data, and 33 in total)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualize the influence from each borough\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "fig = sm.graphics.influence_plot(model, ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Influence tests\n",
    "from statsmodels.stats.outliers_influence import OLSInfluence\n",
    "test_class = OLSInfluence(model)\n",
    "test_class.dfbetas[:5,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It will return a ndarray ($m*n$), where $m$ is the number you defined (here we give 5, but you may change it by yourself), $n$ is the number of independent variables (4 for our model), the value is studentized residuals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Leverage-Resid2 Plot'>Leverage-Resid2 Plot<a/>\n",
    "Closely related to the influence_plot is the leverage-resid2 plot. Let's have a look of the leverage plot first. In order to interpret the diagnostic plot better, you are encourage to explore the [leverage](http://sphweb.bumc.bu.edu/otlt/MPH-Modules/BS/R/R5_Correlation-Regression/R5_Correlation-Regression7.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "fig = sm.graphics.plot_leverage_resid2(model, ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot_regress_exog function is a convenience function that gives a 2x2 plot containing the dependent variable and fitted values with confidence intervals vs. the independent variable chosen, the residuals of the model vs. the chosen independent variable, a partial regression plot, and a CCPR plot. This function can be used for quickly checking modeling assumptions with respect to a single regressor. Let's take the most significant variable \"private_rent_price\" as an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12,8))\n",
    "fig = sm.graphics.plot_regress_exog(model, 'private_rent_price', fig=fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot_fit function plots the fitted values versus a chosen independent variable. It includes prediction confidence intervals and optionally plots the true dependent variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "fig = sm.graphics.plot_fit(model, \"private_rent_price\", ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More functions to explore from `statsmodel` at http://www.statsmodels.org/stable/index.html.\n",
    "## <a id='Spatial Regression'>Spatial Regression</a>\n",
    "We need to use Moran's I tool to check for spatial autocorrelation firstly before conducting spatial regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_queen = ps.weights.Queen.from_dataframe(gdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mi = ps.Moran(gdf['listing_price'], w_queen, two_tailed=False)\n",
    "print(\"The Statistic Moran's I is: \"+str(\"%.4f\"%mi.I),\n",
    "      \"\\nThe Expected Value for Statistic I is: \"+str(\"%.4f\"%mi.EI),\n",
    "      \"\\nThe Significance Test Value is: \"+str(\"%.4f\"%mi.p_norm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Spatial Lag model'>Spatial Lag model</a>\n",
    "In a similar way to how we have included the spatial lag, one could think the airbnb listing prices surrounding a given property also enter its own price function. Recall your memory on Spatial Lag model from lecture, and use `ML_Lag` class in `PySAL` to estimate this model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sl = ml_lag.ML_Lag(y, X, w_queen, name_y='price', name_x=['number_of_reviews', 'private_rent','ratio_price_earning', 'person_dwelling'], name_w='w_queen', name_ds='borough_airbnb_housing')\n",
    "print(sl.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, results are again very similar in all the other variable. It is also very clear that the estimate of the spatial lag of price is statistically significant. This points to evidence that there are processes of spatial interaction between property owners when they set their price."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Spatial Error model'>Spatial Error model</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "se = ml_error.ML_Error(y, X, w_queen, name_y='price', name_x=['number_of_reviews', 'private_rent','ratio_price_earning', 'person_dwelling'], name_w='w_queen', name_ds='borough_airbnb_housing')\n",
    "print(se.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Task 4'>Task 4</a>\n",
    "Discuss you interpretation on the results from Spatial Lag Model and Spatial Error Model with your neighbors, comparing them with the OLS regression result last week, and summarize your conclusions below:\n",
    "\n",
    "**----------------------------------------------------------------**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Prediction performance of spatial models'>Prediction performance of spatial models</a>\n",
    "We can use the mean squared error (MSE), a standard metric of accuracy in the machine learning literature, to evaluate whether explicitly spatial models are better than traditional, non-spatial ones:\n",
    "### <a id='Task 5'>Task 5</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# repeatly to build up your aspatial ordinary least square regression model\n",
    "m1 = ???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error as mse\n",
    "\n",
    "mses = pd.Series({'OLS': mse(y, m1.predy.flatten()), \n",
    "                  'SL': mse(y, sl.predy.flatten()), \n",
    "                  'SE': mse(y, se.predy.flatten())\n",
    "                    })\n",
    "mses.sort_values()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inclusion of the spatial lag of price marginally reduces the MSE, however, does a better job at improving the accuracy of the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='GWR Prediction'>GWR Prediction</a>\n",
    "Geographically weighted regression (**GWR**) can fit Gaussian, Poisson, and logistic models to estimate a 'GWR Results' object. Now let's compare the maps $before$ and $after$ applying the GWR prediction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pysal.contrib.gwr.sel_bw import Sel_BW\n",
    "from pysal.contrib.gwr.gwr import GWR\n",
    "from pysal.contrib.glm.family import Gaussian\n",
    "from scipy.stats import pearsonr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vmin, vmax = np.min(gdf['listing_price']), np.max(gdf['listing_price']) \n",
    "ax = gdf.plot('listing_price', vmin=vmin, vmax=vmax, figsize=(8,8), cmap='Reds')\n",
    "ax.set_title('listing_price'+' t-vals')\n",
    "fig = ax.get_figure()\n",
    "cax = fig.add_axes([1.0, 0.3, 0.02, 0.4]) # the position and size of colormap legend bar\n",
    "sm_price = plt.cm.ScalarMappable(norm=plt.Normalize(vmin=vmin, vmax=vmax), cmap='Reds')\n",
    "sm_price._A = []\n",
    "fig.colorbar(sm_price, cax=cax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prep data into design matrix and coordinates\n",
    "\n",
    "# Dependent variable\n",
    "y = gdf.listing_price\n",
    "y = np.array(y).reshape(-1,1) # make array change the list format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Design matrix - covariates - intercept added automatically\n",
    "amou = np.array(gdf.airbnb_amount).reshape(-1,1)\n",
    "rent = np.array(gdf.private_rent_price).reshape(-1,1)\n",
    "ratio = np.array(gdf.ratio_price_earn).reshape(-1,1)\n",
    "pdwell = np.array(gdf.per_dwell).reshape(-1,1)\n",
    "\n",
    "X = np.hstack([amou, rent, ratio, pdwell])\n",
    "labels = ['Intercept', 'amount', 'rent', 'ratio', 'dwellresid']\n",
    "\n",
    "# standardization\n",
    "X_s = (X - X.mean(axis=0)) / X.std(axis=0)\n",
    "y_s = (y - y.mean(axis=0)) / y.std(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Coordinates for calibration points\n",
    "\n",
    "def getXY(pt):\n",
    "    return (pt.x, pt.y)\n",
    "centroidseries = gdf['geometry'].centroid\n",
    "u,v = [list(t) for t in zip(*map(getXY, centroidseries))]\n",
    "\n",
    "coords = list(zip(u,v))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Find optimal bandwidth using golden section search to minimize AICc\n",
    "\n",
    "#Instantiate bandwidth selection class - bisquare NN (adaptive)\n",
    "selector = Sel_BW(coords, y_s, X_s, kernel='bisquare', fixed=False)\n",
    "\n",
    "#Find optimal bandwidth by minimizing AICc using golden section search algorithm\n",
    "bw = selector.search(search='golden_section', criterion='AICc')\n",
    "\n",
    "print (bw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Instantiate GWR model and then estimate parameters and diagnostics using fit method\n",
    "model = GWR(coords, y, X, bw, family=Gaussian(), fixed=False, kernel='bisquare')\n",
    "results = model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Results in a set of mappable results \n",
    "results.params.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Map Parameter estimates and T-vals for each covariate\n",
    "for param in range(results.params.shape[1]):\n",
    "    gdf[str(param)] = results.params[:,param]\n",
    "    vmin, vmax = np.min(gdf[str(param)]), np.max(gdf[str(param)]) \n",
    "    ax = gdf.plot(str(param), vmin=vmin, vmax=vmax, figsize=(8,8), cmap='YlOrRd')\n",
    "    ax.set_title(labels[param] + ' Estimates')\n",
    "    fig = ax.get_figure()\n",
    "    cax = fig.add_axes([1.0, 0.3, 0.03, 0.4])\n",
    "    sm = plt.cm.ScalarMappable(norm=plt.Normalize(vmin=vmin, vmax=vmax), cmap='YlOrRd')\n",
    "    sm._A = []\n",
    "    fig.colorbar(sm, cax=cax)\n",
    "    \n",
    "    gdf[str(param)] = results.tvalues[:,param]\n",
    "    vmin, vmax = np.min(gdf[str(param)]), np.max(gdf[str(param)]) \n",
    "    ax = gdf.plot(str(param), vmin=vmin, vmax=vmax, figsize=(8,8), cmap='Greys')\n",
    "    ax.set_title(labels[param] + ' t-vals')\n",
    "    fig = ax.get_figure()\n",
    "    cax = fig.add_axes([1.0, 0.3, 0.02, 0.4])\n",
    "    sm = plt.cm.ScalarMappable(norm=plt.Normalize(vmin=vmin, vmax=vmax), cmap='Greys')\n",
    "    sm._A = []\n",
    "    fig.colorbar(sm, cax=cax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Map local R-square values which is a weighted R-square at each observation location\n",
    "\n",
    "gdf['localR2'] = results.localR2\n",
    "vmin, vmax = np.min(gdf['localR2']), np.max(gdf['localR2']) \n",
    "ax = gdf.plot('localR2', vmin=vmin, vmax=vmax, figsize=(8,8), cmap='PuBuGn')\n",
    "ax.set_title('Local R-Squared')\n",
    "fig = ax.get_figure()\n",
    "cax = fig.add_axes([1.0, 0.3, 0.02, 0.4])\n",
    "sm = plt.cm.ScalarMappable(norm=plt.Normalize(vmin=vmin, vmax=vmax), cmap='PuBuGn')\n",
    "sm._A = []\n",
    "fig.colorbar(sm, cax=cax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Credits！"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Geocomp",
   "language": "python",
   "name": "gsa2018"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}