Big_Data_Analytics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "rn-_KJ3eqxVy"
   },
   "source": [
    "# The Woof Factor in Zürich\n",
    "---\n",
    "\n",
    "Through use of data freely available from *Stadt Zürich Open Data*, we analyse dogs registered in Zurich and combine this data with other information available about the districts (kreis) that make up Zurich. We perform exploratory analysis using a number of visualisations and a model is developed in which we predict dog breed given related data.\n",
    "\n",
    "*Project Group 23:* Corey Bothwell, Nicoletta Farabullini, Jacob Gelling, Andris Prokofjevs, Qasim Warraich <br>\n",
    "Link to Repo: https://github.com/nfarabullini/Big-Data-Analytics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Configuration and Dependencies\n",
    "---\n",
    "## R Kernel\n",
    "This Jupyter notebook requires use of the R Kernel, which can be installed with a very simple 2 step installation process.\n",
    "\n",
    "The detailed instructions can be found here: https://irkernel.github.io/installation/\n",
    "\n",
    "## Browser\n",
    "We reccomend a **chromium based browser** to view this notebook. We have discovered an **issue** rendering one of the visualisations **in firefox**.\n",
    "\n",
    "## Libraries\n",
    "The required packages must first be installed. On Linux, this requires the system to have Curl (for communicating with the Yandex Translate API) and GDAL (for generating the map vizualisations) installed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HYonn-cUpjzv"
   },
   "source": [
    "# R Packages\n",
    "---\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 625
    },
    "colab_type": "code",
    "id": "8uvzQkZ1pf5R",
    "outputId": "0102d1c9-c4a0-4188-cef8-abae10a9b2d0",
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "install.packages(\"readr\")\n",
    "install.packages(\"data.table\")\n",
    "install.packages(\"RYandexTranslate\")\n",
    "install.packages(\"plyr\")\n",
    "install.packages(\"dplyr\")\n",
    "install.packages(\"stringr\")\n",
    "install.packages(\"class\")\n",
    "install.packages(\"e1071\")\n",
    "install.packages(\"ggplot2\")\n",
    "install.packages(\"reshape2\")\n",
    "install.packages(\"leaflet\")\n",
    "install.packages(\"rgdal\")\n",
    "install.packages(\"stringi\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Libraries used throughout the analysis are imported here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(readr)\n",
    "library(data.table)\n",
    "library(RYandexTranslate)\n",
    "library(plyr)\n",
    "library(dplyr)\n",
    "library(stringr)\n",
    "library(class)\n",
    "library(e1071)\n",
    "library(ggplot2)\n",
    "library(reshape2)\n",
    "library(leaflet)\n",
    "library(rgdal)\n",
    "library(stringi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Aquisition\n",
    "---\n",
    "\n",
    "All data was aquired from Stadt Zürich Open Data. The below table details data sources used and their respective URL. This data is processed, cleaned and merged in the next section.\n",
    "\n",
    "| Content           | URL                                                                                      |\n",
    "|-------------------|------------------------------------------------------------------------------------------|\n",
    "| Dog Registrations | https://data.stadt-zuerich.ch/dataset/sid_stapo_hundebestand                             |\n",
    "| Wealth Data       | https://data.stadt-zuerich.ch/dataset/fd_median_einkommen_quartier_od1003                |\n",
    "| Income Data       | https://data.stadt-zuerich.ch/dataset/fd_median_vermoegen_quartier_od1004                |\n",
    "| Education Data    | https://data.stadt-zuerich.ch/dataset/bfs_bev_bildungsstand_statquartier_seit1970_od1012 |\n",
    "| Home Type Data    | https://data.stadt-zuerich.ch/dataset/bau_best_geb_whg_bev_gebaeudeart_quartier_seit2008 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "fxWlXtTDj5EE"
   },
   "source": [
    "# Data Processing and Cleaning\n",
    "---\n",
    "\n",
    "This section is used for importing the raw data from CSV files into R data tables. Appropriate data merging and cleaning is performed.\n",
    "\n",
    "Also performed is automatic translations from German to English for some of the column content that would otherwise be too laborious to do by hand."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "NDSSGT4R4aca"
   },
   "source": [
    "Rather than running data cleaning each time we run the notebook during development, which can be slow due to the Yandex Translate API, an R data image of the alreadly cleaned data provided with the Notebook file can be loaded instead by running the following block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "DYmD8dDM4QoG"
   },
   "outputs": [],
   "source": [
    "load(\"dogs.Rdata\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5BqJj3pYm1jr"
   },
   "source": [
    "## Import Dog Data\n",
    "\n",
    "The dog registration data for 2020, for which the source is available at [Stadt Zürich Open Data](https://data.stadt-zuerich.ch/dataset/sid_stapo_hundebestand/resource/a05e2101-7997-4bb5-bed8-c5a61cfffdcf), is imported and converted to a R data table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 510
    },
    "colab_type": "code",
    "id": "ynxU6g3GnNLJ",
    "outputId": "85b79030-11ba-41ef-83f2-2b51ad2d62db"
   },
   "outputs": [],
   "source": [
    "dogs2020 <- data.table(read_csv(\"data_sources/20200306_hundehalter.csv\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data cleaning is now performed. Due to what is likely a typographical error, a dog with a district value of 8 is removed (non-existant stadtquartier).\n",
    "\n",
    "Data we are uninterested in and rows that are incomplete are also removed, however for the 2020 dog dataset we find that 0 rows are ommited."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dog w. district \"8\", typographical error? We remove it.\n",
    "dogs2020 <- dogs2020[(dogs2020$STADTQUARTIER!=8), ]\n",
    "\n",
    "# Removing unnececesary columns\n",
    "dogs2020[, c(\"RASSE1_MISCHLING\", \"RASSE2\", \"RASSE2_MISCHLING\"):=NULL]\n",
    "\n",
    "# If a row has a NA entry in one of the cells, remove the entire row\n",
    "dogs2020 <- na.omit(dogs2020)\n",
    "# 0 rows ommited, still leave the code in place in case we change data basis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to the international makeup of our team members, we translate the column names from German to English."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rename columns\n",
    "setnames(dogs2020,\n",
    "         old = c(\"HALTER_ID\", \"ALTER\", \"GESCHLECHT\", \"STADTQUARTIER\", \"STADTKREIS\",   \"RASSE1\", \"RASSENTYP\",  \"GEBURTSJAHR_HUND\", \"GESCHLECHT_HUND\", \"HUNDEFARBE\"),\n",
    "         new = c(\"OWNER_ID\",  \"AGE\",   \"SEX\",        \"DISTRICT\",      \"DISTRICT_BIG\", \"BREED\",  \"BREED_TYPE\", \"YOB_DOG\",          \"SEX_DOG\",         \"COLOR_DOG\")\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because stadtquartier is a more granular district segmentation than stadtkreis, and not all datasets can be merged by stadtkreis, we used it instead. We found that it is hard to view on a map when the dataset does not contain normal district names, so we extract the district names from the wealth dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import district names\n",
    "district_names <- data.table(read_csv(\"data_sources/wir100od1004.csv\"))\n",
    "district_names <- unique(district_names[, QuarSort, QuarLang])\n",
    "\n",
    "# Rename columns\n",
    "setnames(district_names,\n",
    "         old = c(\"QuarSort\", \"QuarLang\"),\n",
    "         new = c(\"DISTRICT\", \"DISTRICT_NAME\")\n",
    "        )\n",
    "\n",
    "dogs2020 <- merge(dogs2020, district_names , by = \"DISTRICT\", all.x = T)\n",
    "\n",
    "# Remove unused variable\n",
    "rm(district_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "qMHV7jrUm-E9"
   },
   "source": [
    "## Import Wealth Data\n",
    "\n",
    "Next the wealth data, for which the source is available at [Stadt Zürich Open Data](https://data.stadt-zuerich.ch/dataset/fd_median_einkommen_quartier_od1003), is imported and converted to a R data table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 221
    },
    "colab_type": "code",
    "id": "mx1asHDGnDuh",
    "outputId": "ba02bb01-58c5-4962-fa24-240d9a2304f1"
   },
   "outputs": [],
   "source": [
    "wealth <- data.table(read_csv(\"data_sources/wir100od1004.csv\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The latest data from this source is for 2017, so we select only this wealth data, store the average and aggregate the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We don't have data for 2020, the freshest data is for 2017\n",
    "wealth <- wealth[SteuerJahr == 2017,]\n",
    "\n",
    "# Create a new column to store the average\n",
    "wealth[, wealth50 := SteuerVermoegen_p50]\n",
    "\n",
    "# Replace the family values in this column with same divided by 2 for normalization\n",
    "wealth$wealth50[wealth$SteuerTarifSort == 1] <- wealth$wealth50[wealth$SteuerTarifSort == 1]/2\n",
    "\n",
    "# Aggregate data for family status\n",
    "# new table = old table, select mean of wealth50 (ignore NA), aggregate it by quartal\n",
    "wealth_merge <- wealth[,mean(wealth50, na.rm = T), by=QuarSort]\n",
    "\n",
    "# Leaving only quartals that have dogs in them\n",
    "wealth_merge <- wealth_merge[wealth_merge$QuarSort %in% dogs2020$DISTRICT]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We again translate the column names from German to English and then merge the data with the dogs data table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rename columns\n",
    "setnames(wealth_merge,\n",
    "         old = c(\"V1\",           \"QuarSort\"),\n",
    "         new = c(\"WEALTH_T_CHF\", \"DISTRICT\")\n",
    "        )\n",
    "\n",
    "# Merge wealth into dogs dataset\n",
    "dogs2020 <- merge(dogs2020, wealth_merge, by = \"DISTRICT\", all.x = T)\n",
    "# ATTENTION! FYI: not all districts have wealth values! \n",
    "\n",
    "# Remove unused variables\n",
    "rm(wealth, wealth_merge)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "e9zuqYq6nW-J"
   },
   "source": [
    "## Import Income Data\n",
    "\n",
    "Next the income data, for which the source is available at [Stadt Zürich Open Data](https://data.stadt-zuerich.ch/dataset/fd_median_vermoegen_quartier_od1004), is imported and converted to a R data table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 221
    },
    "colab_type": "code",
    "id": "KGQ9qrWsnZoK",
    "outputId": "d36d841c-4300-4d53-b62c-d56ffa5acd7a"
   },
   "outputs": [],
   "source": [
    "income <- data.table(read_csv(\"data_sources/wir100od1003.csv\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, the latest data from this source is for 2017, so we select only this income data, store the average and aggregate the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We don't have data for 2020, the freshest data is for 2017\n",
    "income <- income[SteuerJahr == 2017,]\n",
    "\n",
    "# Create a new column to store the average\n",
    "income[, incomep50 := SteuerEInkommen_p50]\n",
    "\n",
    "# Replace the family values in this column with same divided by 2 for normalization\n",
    "income$incomep50[income$SteuerTarifSort == 1] <- income$incomep50[income$SteuerTarifSort == 1]/2\n",
    "\n",
    "# Aggregate data for family status\n",
    "# new table = old table, select mean of incomep50 (ignore NA), aggregate it by quartal\n",
    "income_merge <- income[,mean(incomep50, na.rm = T), by=QuarSort]\n",
    "\n",
    "# Leaving only quartals that have dogs in them\n",
    "income_merge <- income_merge[income_merge$QuarSort %in% dogs2020$DISTRICT]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We again translate the column names from German to English and then merge the data with the dogs data table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rename columns\n",
    "setnames(income_merge,\n",
    "         old = c(\"V1\",           \"QuarSort\"),\n",
    "         new = c(\"INCOME_T_CHF\", \"DISTRICT\")\n",
    "        )\n",
    "\n",
    "# Merge income into dogs dataset\n",
    "dogs2020 <- merge(dogs2020, income_merge, by = \"DISTRICT\", all.x = T)\n",
    "# ATTENTION! FYI: not all districts have income values! \n",
    "\n",
    "# Remove unused variables\n",
    "rm(income, income_merge)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "I24YVJ_VncAY"
   },
   "source": [
    "## Import Education Data\n",
    "\n",
    "The import of education data is straightforward, requiring only a long to wide reshape of the data. The data source is available at [Stadt Zürich Open Data](https://data.stadt-zuerich.ch/dataset/bfs_bev_bildungsstand_statquartier_seit1970_od1012)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 187
    },
    "colab_type": "code",
    "id": "4NssSsYpneo0",
    "outputId": "105009b9-96f9-4f5a-b71c-18c1cd01ad35"
   },
   "outputs": [],
   "source": [
    "education <- data.table(read_csv(\"data_sources/bil101od1012 (2).csv\"))\n",
    "\n",
    "# Long to wide education reshape\n",
    "education <- dcast(education, RaumSort ~ Bildungsstand, value.var = \"AntBev\")\n",
    "\n",
    "# Rename columns\n",
    "setnames(education,\n",
    "         old = c(\"RaumSort\", \"Obligatorische Schule\",   \"Sekundarstufe II\",     \"Tertiärstufe\"),\n",
    "         new = c(\"DISTRICT\", \"BASIC_SCHOOL_PERCENTAGE\", \"GYMNASIUM_PERCENTAGE\", \"UNIVERSITY_PERCENTAGE\")\n",
    "        )\n",
    "\n",
    "# Merge income into dogs dataset\n",
    "dogs2020 <- merge(dogs2020, education, by = \"DISTRICT\", all.x = T)\n",
    "\n",
    "# Remove unused variable\n",
    "rm(education)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "aHUmAB4pnjK5"
   },
   "source": [
    "## Import Home Type Data\n",
    "\n",
    "Home type data, such as number of resident buildings in a district, is now imported. The data source is available at [Stadt Zürich Open Data](https://data.stadt-zuerich.ch/dataset/bau_best_geb_whg_bev_gebaeudeart_quartier_seit2008/resource/3850add1-264c-4993-98cd-d8a9ba87ee25). 2019 is the latest data available. Each building type is summed and a long to wide reshape of the data is performed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 374
    },
    "colab_type": "code",
    "id": "C80_N4ZgnmDm",
    "outputId": "45d43f21-1696-4bca-f310-6e27c77ed11b"
   },
   "outputs": [],
   "source": [
    "home_type <- data.table(read_csv(\"data_sources/bau_best_geb_whg_bev_gebaeudeart_quartier_seit2008.csv\")) \n",
    "\n",
    "# We don't have data for 2020, the freshest data is for 2019\n",
    "home_type <- home_type[Jahr == 2019,]\n",
    "\n",
    "# Count the sum of every type of building\n",
    "home_type <- home_type[, sum(AnzGeb), by = list(QuarSort,GbdArtPubName)]\n",
    "\n",
    "# Rename columns\n",
    "setnames(home_type,\n",
    "         old = c(\"QuarSort\", \"GbdArtPubName\", \"V1\"), \n",
    "         new = c(\"DISTRICT\", \"Hometype\",      \"Number_homes\")\n",
    "        )\n",
    "\n",
    "# Long to wide hometype reshape\n",
    "home_type <- data.table::dcast(home_type, DISTRICT ~ Hometype, value.var = \"Number_homes\")\n",
    "\n",
    "# Translate hometype columns\n",
    "setnames(home_type,\n",
    "         old = c(\"Produktions- und Lagergebäude\", \"Mehrfamilienhäuser\", \"Einfamilienhäuser\",   \"Infrastrukturgebäude\",     \"Kleingebäude\",    \"Kommerzielle Gebäude\", \"Spezielle Wohngebäude\"), \n",
    "         new = c(\"FACTORIES_AND_WAREHOUSES\",      \"APARTMENTS\",         \"SINGLE_FAMILY_HOMES\", \"INFRASTRUCTURE_BUILDINGS\", \"SMALL_BUILDINGS\", \"COMMERCIAL_BUILDINGS\", \"SPECIAL_ACCOMODATION\")\n",
    "        )\n",
    "\n",
    "# Remove unnececesary column\n",
    "home_type[, Unbekannt:=NULL]\n",
    "\n",
    "# Merge home types into dogs dataset\n",
    "dogs2020 <- merge(dogs2020, home_type, by = \"DISTRICT\", all.x = T)\n",
    "\n",
    "# Remove unused variable\n",
    "rm(home_type)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "55CKOQRFnrCr"
   },
   "source": [
    "## Import Population Data\n",
    "\n",
    "Next population data is imported. The data source is available at [Stadt Zürich Open Data](https://www.stadt-zuerich.ch/prd/de/index/statistik/themen/bevoelkerung/bevoelkerungsentwicklung/kreise-und-quartiere.html#daten). The source CSV file is a little messy, so some additional data cleaning is required. When importing, the first 8 lines are skipped for the parser to reach the CSV header. As one of the columns are unlabelled, a warning is produced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 204
    },
    "colab_type": "code",
    "id": "RLQ9mG4vntd5",
    "outputId": "30696c37-6fa7-47ea-d850-ef27c42ab3b4"
   },
   "outputs": [],
   "source": [
    "# Warning shown due to formatting of souce CSV header, can be ignored\n",
    "pop_per_district <- data.table(read_csv(\"data_sources/2019-Table_1.csv\", col_names = TRUE, skip = 8))\n",
    "\n",
    "# Set column names\n",
    "setnames(pop_per_district, old = \"X1\", new = \"DISTRICT_NAME\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data also contains hidden characters acting as whitespace which the standard regex matching pattern ```\\\\s``` does not always catch. For this reason we remove all non-numeric characters with a regex function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This Reg-Ex Matching Function Removes all Whitespaces and not only conventional ones (bad data design)\n",
    "pop_per_district[,2:5] <- data.table(apply(pop_per_district[,2:5], 2, function(x) gsub('[^0-9.]', '', x)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As normal, we continue by merging the population data into the dogs data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Join by district name (perfect match with no NAs)\n",
    "dogs2020 <- merge(dogs2020, pop_per_district, by = \"DISTRICT_NAME\", all.x = T)\n",
    "\n",
    "# Remove unused variable\n",
    "rm(pop_per_district)\n",
    "\n",
    "# Rename columns\n",
    "setnames(dogs2020,\n",
    "         old = c(\"Total\",            \"Schweizer/-innen\", \"Ausländer/-innen\",   \"Anteil ausländische\\nBevölkerung (%)\"),\n",
    "         new = c(\"TOTAL_POPULATION\", \"SWISS_POPULATION\", \"FOREIGN_POPULATION\", \"FOREIGN_POPULATION_PERCENTAGE\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The columns are then cast as numeric data types for use in mathematical operations later in the analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cast columns as numeric\n",
    "dogs2020$TOTAL_POPULATION <- as.numeric(dogs2020$TOTAL_POPULATION)\n",
    "dogs2020$SWISS_POPULATION   <- as.numeric(dogs2020$SWISS_POPULATION)\n",
    "dogs2020$FOREIGN_POPULATION   <- as.numeric(dogs2020$FOREIGN_POPULATION  )\n",
    "dogs2020$FOREIGN_POPULATION_PERCENTAGE   <- as.numeric(dogs2020$FOREIGN_POPULATION_PERCENTAGE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nvMuam6BnzHX"
   },
   "source": [
    "## Translate Dataframes\n",
    "Using the Yandex Translate API, values in columns are automatically translated from German to English."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3UqXvple8DvM"
   },
   "source": [
    "Firstly the Yandex translate package is fixed, as per the [solution found on GitHub](https://github.com/mukul13/RYandexTranslate/issues/2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "7Wt8fMK68L8c"
   },
   "outputs": [],
   "source": [
    "translate = function (api_key, text = \"\", lang = \"\") \n",
    "{\n",
    "  url = \"https://translate.yandex.net/api/v1.5/tr.json/translate?\"\n",
    "  url = paste(url, \"key=\", api_key, sep = \"\")\n",
    "  if (text != \"\") {\n",
    "    url = paste(url, \"&text=\", text, sep = \"\")\n",
    "  }\n",
    "  if (lang != \"\") {\n",
    "    url = paste(url, \"&lang=\", lang, sep = \"\")\n",
    "  }\n",
    "  url = gsub(pattern = \" \", replacement = \"%20\", x = url)\n",
    "  d = RCurl::getURL(url, ssl.verifyhost = 0L, ssl.verifypeer = 0L, .encoding = \"UTF-8\")\n",
    "  d = jsonlite::fromJSON(d)\n",
    "  d$code = NULL # Remove unused variables\n",
    "  d\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "rsjgeXUK8Yhq"
   },
   "source": [
    "As the Yandex Translate API requires an API key, we define this here. It is connected to Andris' personal account, so please do not disclose it elsewhere."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Andris' personal Yandex Translate API key\n",
    "# PLEASE DO NOT DISCLOSE ELSEWHERE\n",
    "api_key <- \"trnsl.1.1.20200515T134653Z.f9fb709ac3e94036.783aefa609692b463a79b5827d5c0e7f2d037a8c\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define the list of columns for which we want to translate their contents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "column_list <- c(\"BREED\", \"COLOR_DOG\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then loop through columns that we wish to translate and pass these to the Yandex Translate API. As translation can be slow, we pass only unique values from the columns, preventing a translation call for every row, even if the value was translated before. The output is saved to the dogs data table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "id": "NQfpAlpIn1xI",
    "outputId": "1eff4904-573a-4253-ce95-6ed13b2af689"
   },
   "outputs": [],
   "source": [
    "# Loop through columns in column_list\n",
    "for (column_names in column_list) {\n",
    "  # Get all unique values in the coresponding column\n",
    "  unique_values <- unique(dogs2020[, get(column_names)])\n",
    "    \n",
    "  # replacing special character (apostrophy) in one of the breed names, to avoid encoding error in jupyter\n",
    "  unique_values <- str_replace(unique_values, \"\\032\" , \"\")\n",
    "    \n",
    "  # Loop through unique values to be translated\n",
    "  for (unique_num in 1:length(unique_values)) {\n",
    "    \n",
    "    # Print progress for debugging\n",
    "    print(column_names)\n",
    "    print(paste(unique_num, \" out of \", length(unique_values)))\n",
    "      \n",
    "    # data.table synthaxis code\n",
    "    dogs2020[\n",
    "      # Select all rows where value is equal to current unique value\n",
    "      # left part \"get\" is needed in order to use dynamic column_names\n",
    "      dogs2020[,get(column_names)] == unique_values[unique_num],\n",
    "      # Replace old value with translated one\n",
    "      # Left part is \"eval\" is needed in order to use dynamic column_names\n",
    "      # Right part is a Yandex Translate function, providing api_key and text to be translated\n",
    "      eval(column_names) := translate(api_key, text=stringi::stri_unescape_unicode(unique_values[unique_num]),\n",
    "      # Specify language pair and translation direction\n",
    "      lang=\"de-en\"\n",
    "      # Extract only output text from the return given by Yandex\n",
    "      )$text]\n",
    "  \n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We manually replace sex and sex_dog variable to match the English language conventions for sex abbreviations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dogs2020[SEX == \"w\", SEX := \"f\"]\n",
    "dogs2020[SEX_DOG == \"w\", SEX_DOG := \"f\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HfomTPUf4s2X"
   },
   "source": [
    "## Save RData Image\n",
    "An R data image is saved here so that instead of running the entire data cleaning code block again in future runs, the image can be loaded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "l0-a80wa43XP"
   },
   "outputs": [],
   "source": [
    "save.image(\"dogs.Rdata\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Column Overview\n",
    "\n",
    "**DISTRICT_NAME** = String representation of the district's name <br>\n",
    "**DISTRICT** = Numeric representation of a sub district <br>\n",
    "**OWNER_ID** = Unique owner ID <br>\n",
    "**AGE** = Age group of owner (Grouped by decade)  <br>\n",
    "**SEX** = Gender of owner  <br>\n",
    "**DISTRICT_BIG** = Parent of subdistrict  <br>\n",
    "**BREED** = String representation of breed of dog  <br>\n",
    "**BREED_TYPE** = 3 possible values (K,I,II) K = Small Breeds, I = Unrestricted Breeds, II = Restricted Breeds.<br>\n",
    "**YOB_DOG** = Year of dog's birth <br>\n",
    "**SEX_DOG** = Gender of Dog (M or F) <br>\n",
    "**COLOR_DOG** = String representation of dog colour <br>\n",
    "**WEALTH_T_CHF** = Avg district resident wealth in thousands of CHF. <br>\n",
    "**INCOME_T_CHF** = Avg district resident income in thousands of CHF. <br>\n",
    "**BASIC_SCHOOL_PERCENTAGE** = Avg district resident's who have <= Basic school education. <br>\n",
    "**GYMNASIUM_PERCENTAGE** = Avg district resident's who have <= Gymnasium level education. <br>\n",
    "**UNIVERSITY_PERCENTAGE** = Avg district resident's who have <= University level education. <br>\n",
    "**SINGLE_FAMILY_HOMES** = Number of single family homes in a given sub-district.  <br>\n",
    "**INFRASTRUCTURE_BUILDINGS** = Number of infrastructure buildings in a given sub-district. <br>\n",
    "**SMALL_BUILDINGS** = Number of small buildings in a given sub-district.  <br>\n",
    "**COMMERCIAL_BUILDINGS** = Number of commercial buildings in a given sub-district.<br>\n",
    "**APARTMENTS** = Number of apartment buildings in a given sub-district.<br>\n",
    "**FACTORIES_AND_WAREHOUSES** = Number of factories and warehouses in a given sub-district.<br>\n",
    "**SPECIAL_ACCOMODATION** = Number of special accomodations (eg: dorms, carehomes) in a given sub-district.<br>\n",
    "**TOTAL_POPULATION** = Total population of a given sub-district.<br>\n",
    "**SWISS_POPULATION** = Total Swiss population of a given sub-district. <br>\n",
    "**FOREIGN_POPULATION** = Total foreign population of a given sub-district.<br>\n",
    "**FOREIGN_POPULATION_PERCENTAGE** = Total foreign population percentage of a given sub-district.<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Limitations\n",
    "\n",
    "It is important to discuss some limitations inherent in the data. To wit, the import data often cannot be matched up exactly year to year. Concretely, sources in our data imports maintain years between 2017 to 2020. This could skew the results in an inaccurate manner if large variations in the measured data (for instance, wealth) were present in the unmeasured years. A critical assumption in our analysis is therefore that these variables do not vary drastically from year to year, and that an appropriate 'general' model can be inferred with years that do not match one-to-one.\n",
    "\n",
    "Futhermore, within our income & wealth datasets, some districts do not maintain a value across one or more of the three categorical dimensions (family type) of the appropriate wealth or income class. To account for this, we have computed an average across all available categorical dimensions. If these missing values, which function as omitted variables, were to significantly alter the underlying distribution, our results may not hold in practice.\n",
    "\n",
    "Finally, it is important to note that the data is effectively split into two sets of granularities. The first are those characteristics of individual dogs and owners themselves, (Sex of Dog, Age of Owner, etc.) while the second is are those characteristics computed across districts in aggregate. These include features such as Average Income, Average Wealth, etc. and those features can be utilised as proxies on individual owners in those districts. While less severe to underlying stastistical distributions than the other limitations mentioned in this paragraph, it behooves practicioners to keep this distinction in granularity in mind when conducting analysis. We account for this important distinction by effectively splitting our analysis into two classes, that which applies to owners & dogs individually, and that which applies to districts in aggregate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "XlVYqq67tNF1"
   },
   "source": [
    "## Backup Output\n",
    "Save a backup of the dogs data set in another variable. This is useful as the data is modified during the investigation. **This backup is not optional, please execute this block below.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "KdoFnYpZtVpH"
   },
   "outputs": [],
   "source": [
    "dogs2020_backup = copy(dogs2020)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary Statistics\n",
    "\n",
    "Finally, before we continue, let us compute some basic summary statistics from the dataset that will give us an overview of the distribution of dogs and their owners in Zuerich and that will help us later on.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_gender <- unique(dogs2020$SEX)\n",
    "total_gender <- lapply(unique_gender, function(x){length(which(dogs2020$SEX == x))})\n",
    "percent_gender <- round(unlist(total_gender)/sum(unlist(total_gender))*100)\n",
    "\n",
    "unique_age <- unique(dogs2020$AGE)\n",
    "total_age <- lapply(unique_age, function(x){length(which(dogs2020$AGE == x))})\n",
    "percent_age <- round(unlist(total_age)/sum(unlist(total_age))*100)\n",
    "\n",
    "# Table illustrates the percentage of males and females owning a dog\n",
    "table_gender <- data.table(\"Owner Gender\" = unique_gender, \"Percentage %\" = percent_gender)\n",
    "\n",
    "# Table illustrates the percentage of dog ownership for each age group\n",
    "table_age <- data.table(\"Owner Age Group\" = unique_age, \"Percentage %\" = percent_age)\n",
    "setorder(table_age, \"Owner Age Group\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at those tables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table_gender\n",
    "table_age"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wow! So Females own dogs in a much higher amount than males for our sample. (population of Zuerich) This will play an important role in our analysis later on. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "N--2CUsM0CDE"
   },
   "source": [
    "\n",
    "# Visualisation\n",
    "\n",
    "### Maps and Charts\n",
    "In order to effectively perform some exploratory data analysis, we decide to build a visualisation over the districts of Zürich. We are interested in providing a mental map for how wealth is distributed across the city. A lot of the analysis presented in this notebook hinges on economic analysis so we thought this would be a fitting visualisation to accompany the rest of the material present in the notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Here we group the subdistricts of Zurich together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "TFBq2zS1Ewp8"
   },
   "outputs": [],
   "source": [
    "zh_rg <- readOGR(\"./data_sources/stzh.adm_stadtkreise_v.json\")\n",
    "\n",
    "# group sub-districts together\n",
    "list_districts <- list(\n",
    "  District_1 <- c(\n",
    "    \"Rathaus\",\n",
    "    \"Hochschulen\",\n",
    "    \"Lindenhof\",\n",
    "    \"City\"\n",
    "  ),\n",
    "  District_2\t<- c(\n",
    "    \"Wollishofen\",\n",
    "    \"Leimbach\",\n",
    "    \"Enge\"\n",
    "  ),\n",
    "  District_3 <- c(\n",
    "    \"Alt-Wiedikon\",\n",
    "    \"Friesenberg\",\n",
    "    \"Sihlfeld\"\n",
    "  ),\n",
    "  District_4 <- c(\n",
    "    \"Werd\",\n",
    "    \"Langstrasse\",\n",
    "    \"Hard\"\n",
    "  ),\n",
    "  District_5 <- c(\n",
    "    \"Gewerbeschule\",\n",
    "    \"Escher Wyss\"\n",
    "  ),\n",
    "  District_6\t<- c(\n",
    "    \"Unterstrass\",\n",
    "    \"Oberstrass\",\n",
    "    \"Unterstrass\"\n",
    "  ),\n",
    "  District_7 <- c(\t\t\n",
    "    \"Fluntern\",\n",
    "    \"Hottingen\",\n",
    "    \"Hirslanden\",\n",
    "    \"Witikon\"\n",
    "  ),\n",
    "  District_8 <- c(\n",
    "    \"Seefeld\",\n",
    "    \"Mühlebach\",\n",
    "    \"Weinegg\"\n",
    "  ),\n",
    "  District_9 <- c(\t\n",
    "    \"Albisrieden\",\n",
    "    \"Altstetten\"\n",
    "  ),\n",
    "  District_10 <- c(\t\n",
    "    \"Höngg\",\n",
    "    \"Wipkingen\"\n",
    "  ),\n",
    "  District_11 <- c(\t\n",
    "    \"Affoltern\",\n",
    "    \"Oerlikon\",\n",
    "    \"Seebach\"\n",
    "  ),\n",
    "  District_12 <- c(\n",
    "    \"Saatlen\",\n",
    "    \"Schwamendingen-Mitte\",\n",
    "    \"Hirzenbach\"\n",
    "  )\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### We then calculate the average wealth per city districts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_wealth <- unlist(lapply(seq_len(length(list_districts)), function (z) {\n",
    "  mean(unlist(lapply(list_districts[[z]], function (y) {\n",
    "    mean(unlist(lapply(y, function (x) {\n",
    "      dogs2020[which(x == dogs2020$DISTRICT_NAME),]$WEALTH_T_CHF\n",
    "    })))\n",
    "  })))\n",
    "}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### This block handles the colouring and generation of a leaflet map that will show a district wise wealth map of Zürich."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "qM9dsO-aE8e0"
   },
   "outputs": [],
   "source": [
    "bins <- c(0, 25, 30, 50, 60, 75, 100, 120, 150, 175)\n",
    "pal <- colorBin(\"Greens\", domain = avg_wealth, bins = bins)\n",
    "\n",
    "map <- leaflet(zh_rg) %>%\n",
    "  addPolygons(fillColor = ~pal(unlist(avg_wealth)), weight = 2, fillOpacity = 0.9, \n",
    "              opacity = 1) %>%\n",
    "  addTiles() %>% \n",
    "  addLegend(colors = pal(unlist(avg_wealth)), labels = zh_rg$kname, title = \"Zurich Districts\", opacity = 1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "## Leaflet map\n",
    "This map shows a district wise breakdown of wealth in the City of Zürich. The darker colours represent lower wealth where as the lighter shades represent wealthier districts. In other words, Kreis 12 is the poorest district, whereas Districts 1 and 2 are the wealthiest. <br>\n",
    "**Note: You may need to execute the block below twice to actually get the map to appear.** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "d1-vswSaoGvZ"
   },
   "source": [
    "# Modelling and Analysis\n",
    "---\n",
    "This section of code explores predicting dog breed, starting from a sophisticated KNN Algorithm (both a naive approach and moving onwards) and then moves on to a more traditional Naive Bayes approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## K-Nearest Neighbours Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our problem is a classification problem because we have labled data and our goal is to predict a class of new data based on known attributes. There are multiple appropriate classification algorithms for such a task, we arbritraily chose KNN to serve as benchmark for the dataset to see what weaknesses we could identify in this model. This was in an attempt to narrow down on a more appropriate model for our data analysis. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5qnDRGFdpwjR"
   },
   "source": [
    "### Model Preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dogs2020$COLOR_DOG[1:20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we recode color to make bigger chunks. **Note: Clearing variable i from the enviroment can produce a warning on first-run of the code, however is necessary on subsequent runs and the error can be safely ignored.** As you can see some dogs are multicoloured so too many distinct colour groups that do not represent the differences between the dogs are formed. To combat this we reduce the colours to only the first listed or primary colour of the dog. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Double-colored dogs will be recoded as single-colored with first color taken as main\n",
    "a <- str_split(dogs2020$COLOR_DOG, \"/\")\n",
    "new_colors <- a[[1]][1]\n",
    "rm(i)\n",
    "for (i in 2:length(a)) {\n",
    "  new_colors <- append(new_colors, a[[i]][1])\n",
    "}\n",
    "\n",
    "# Put new colors into the dogs data table\n",
    "dogs2020$COLOR_DOG <- new_colors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we are going to be putting the data into a knn model, we must first refactor the data into numeric values. This is done for all values we are interested in. The associated labels are saved for breed and breed type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 51
    },
    "colab_type": "code",
    "id": "VET_xDP0p0JH",
    "outputId": "10dce23d-7930-47b9-afe8-6677b2181bcc"
   },
   "outputs": [],
   "source": [
    "# Refactoring text attributes to numeric, so that knn can work with it\n",
    "dogs2020$DISTRICT_NAME <- as.numeric(as.factor(dogs2020$DISTRICT_NAME))\n",
    "dogs2020$AGE <- as.numeric(as.factor(dogs2020$AGE))\n",
    "dogs2020$SEX <- as.numeric(as.factor(dogs2020$SEX))\n",
    "dogs2020$SEX_DOG <- as.numeric(as.factor(dogs2020$SEX_DOG))\n",
    "dogs2020$COLOR_DOG <- as.numeric(as.factor(dogs2020$COLOR_DOG))\n",
    "dogs2020$YOB_DOG <- as.numeric(as.factor(dogs2020$YOB_DOG))\n",
    "\n",
    "# Saving labels in order to see later what numbers represent what breedtypes\n",
    "breedtype_labels <- as.factor(dogs2020$BREED_TYPE)\n",
    "dogs2020$BREED_TYPE <- as.numeric(as.factor(dogs2020$BREED_TYPE))\n",
    "\n",
    "# Saving labels in order to see later what numbers represent what breeds\n",
    "breed_labels <- as.factor(dogs2020$BREED)\n",
    "dogs2020$BREED <- as.numeric(as.factor(dogs2020$BREED))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to train and test the data, we create training (70%) and testing (30%) datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data separation, so that we have testing and training group\n",
    "set.seed(1) # We set this to ensure reproducibility and result consistencey.\n",
    "\n",
    "\n",
    "dat.d <- sample(1:nrow(dogs2020),size=nrow(dogs2020)*0.7,replace = FALSE) # Randomly select 70% of the data.\n",
    "\n",
    "# This is the only setting where we are not \"cheating\" by using dog associated attributes\n",
    "train.dogs <- dogs2020[dat.d, c(\"DISTRICT_NAME\", \"AGE\", \"SEX\")] # 70% training data\n",
    "test.dogs <- dogs2020[-dat.d, c(\"DISTRICT_NAME\", \"AGE\", \"SEX\")] # 30% testing data\n",
    "\n",
    "# Provide labels for training and testing datasets\n",
    "train.dogs_labels <- dogs2020[dat.d,BREED]\n",
    "test.dogs_labels <-dogs2020[-dat.d,BREED]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JlCqPfaRp_Xa"
   },
   "source": [
    "## \"Naive\" KNN\n",
    "Here we perform a simple k-nearest neighbors classification. We try to directly predict breeds by providing all available parameters to the model (excluding population and economic data because they are district level variables and not granular enough to be used on a per breed basis). The accuracy of this direct prediction is roughly 9.22%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the table below you can see the output of the KNN model prediction for breed vs the real breed values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "id": "ZQqHvoD5qCKc",
    "outputId": "0c234121-6626-4423-af11-999789995ed8"
   },
   "outputs": [],
   "source": [
    "# Calculate an approx. k value\n",
    "k_value <-  sqrt(nrow(dogs2020))\n",
    "\n",
    "# Perform knn classification using given training data and testing data\n",
    "knn.test <- knn(train=train.dogs, test=test.dogs, cl=train.dogs_labels, k=k_value)\n",
    "\n",
    "# Calculate accuracy\n",
    "ACC.test <- 100 * sum(test.dogs_labels == knn.test) / NROW(test.dogs_labels)\n",
    "\n",
    "# Print results\n",
    "table(knn.test ,test.dogs_labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Accuracy test output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ACC.test "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ZPZgLDmuqL94"
   },
   "source": [
    "## KNN with Breed Type Grouping (Increasing Model Accuracy)\n",
    "\n",
    "Based on the relatively low accuracy of ~9% we next decided to try and further separate the dataset into smaller chunks with the idea that it would decrease chunk level variety and potentially increase accuracy. As a chunk seperator we used breed types because it's a natural seperator that existed in our dataset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is the point where we exclude all population and income level data due to the granularity issues discussed \n",
    "# above.\n",
    "\n",
    "train.dogs <- dogs2020[dat.d, c(\"DISTRICT_NAME\", \"AGE\", \"SEX\")] # 70% training data\n",
    "test.dogs <- dogs2020[-dat.d,c(\"DISTRICT_NAME\", \"AGE\", \"SEX\")] # Remaining 30% test data\n",
    "\n",
    "train.dogs_labels <- dogs2020[dat.d,BREED_TYPE]\n",
    "test.dogs_labels <-dogs2020[-dat.d,BREED_TYPE]\n",
    "\n",
    "#Approximate k value\n",
    "k_value <-  sqrt(nrow(dogs2020))\n",
    "\n",
    "knn.test <- knn(train=train.dogs, test=test.dogs, cl=train.dogs_labels, k=k_value)\n",
    "\n",
    "ACC.test <- 100 * sum(test.dogs_labels == knn.test)/NROW(test.dogs_labels)\n",
    "\n",
    "table(knn.test ,test.dogs_labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Output of accuracy test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ACC.test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### This apporach yields ~62% accuracy in breed type prediction now we want to predict exact breed in order to do so we need to save breed type prediction. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we save predicted breedtype, by applying algorithm to whole dataset. (this can open potential for overfitting)\n",
    "train.dogs$PREDICTED_BREED_TYPE <- knn(train=train.dogs, test=train.dogs, cl=train.dogs_labels, k=k_value)\n",
    "test.dogs$PREDICTED_BREED_TYPE <- knn(train=train.dogs[,c(\"DISTRICT_NAME\", \"AGE\", \"SEX\")], test=test.dogs, cl=train.dogs_labels, k=k_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Perform KNN on breed directly\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm(i) # Safety measures\n",
    "results <- 0 # Safety measures\n",
    "for (i in 1:length(unique(test.dogs_labels))) # Loops through all unique dog breed types.\n",
    "{ \n",
    "  if (any(train.dogs$PREDICTED_BREED_TYPE == i) == F) {  \n",
    "    next()\n",
    "  }\n",
    "  \n",
    "  train.dogs_labels <- dogs2020[dat.d,BREED] # Label assignment for testing and training datasets,\n",
    "  test.dogs_labels <- dogs2020[-dat.d,BREED] # because it was previously used for breed types not breeds.\n",
    "  \n",
    "  # Calculating the model per breed type.   \n",
    "  knn.test <- knn(train=train.dogs[train.dogs$PREDICTED_BREED_TYPE == i], \n",
    "                  test=test.dogs[test.dogs$PREDICTED_BREED_TYPE == i], \n",
    "                  cl=train.dogs_labels[train.dogs$PREDICTED_BREED_TYPE == i], k=k_value)\n",
    "  \n",
    "  ACC.test <- 100 *\n",
    "    sum(test.dogs_labels[test.dogs$PREDICTED_BREED_TYPE == i] == knn.test)/\n",
    "    NROW(test.dogs_labels[test.dogs$PREDICTED_BREED_TYPE == i])\n",
    "  \n",
    "  table(knn.test ,test.dogs_labels[test.dogs$PREDICTED_BREED_TYPE == i])\n",
    "  \n",
    "  results[i] <- ACC.test\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculating weighted accuracy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(length(test.dogs_labels[test.dogs$PREDICTED_BREED_TYPE == 1]) * results[1]/100\n",
    "  + length(test.dogs_labels[test.dogs$PREDICTED_BREED_TYPE == 3]) * results[3]/100 ) / \n",
    "  (length(test.dogs_labels[test.dogs$PREDICTED_BREED_TYPE == 1]) + length(test.dogs_labels[test.dogs$PREDICTED_BREED_TYPE == 3]))*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The result of ~9%, is not substantially better than approach without breed type grouping."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## \"Naive\" KNN with Breed Frequency Filtering.\n",
    "As a second attempt to improve accuracy we remove certain breeds that are not highly represented in the data in order to reduce the possible number of outcomes. This approach does not compromise the strength of the prediction as certain breeds appear an insufficient (< 20) amount of times in the dataset.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### In order to test it, we have to calculate how many people own each breed:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "breed_filter <- data.table(aggregate(OWNER_ID ~ BREED, data = dogs2020, FUN = function(x){NROW(x)}))\n",
    "# Renaming\n",
    "setnames(breed_filter, old = c(\"OWNER_ID\")\n",
    "         , new = c(\"NUMBER_OWNERS\"))\n",
    "\n",
    "# Sorting by number of owners\n",
    "breed_filter <- breed_filter[order(-rank(NUMBER_OWNERS),)]\n",
    "breed_filter <- breed_filter[breed_filter[, NUMBER_OWNERS > 20]]\n",
    "\n",
    "# Leaving in dataset only desired breeds\n",
    "dogs2020 <- dogs2020[BREED %in% breed_filter[, BREED]]\n",
    "\n",
    "# 1342 out of 7839 rows are out. (6497 left)\n",
    "# 64 breeds out of 313 breeds left"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Recalculation of the naive approach but with the filtering of breeds described above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Now we rerun an exact copy of code in the first try\n",
    "\n",
    "dat.d <- sample(1:nrow(dogs2020),size=nrow(dogs2020)*0.7,replace = FALSE) #random selection of 70% data.\n",
    "\n",
    "# This is the only setting where we are not \"cheating\" by using dog associated attributes\n",
    "train.dogs <- dogs2020[dat.d, c(\"DISTRICT_NAME\", \"AGE\", \"SEX\")] # 70% training data\n",
    "test.dogs <- dogs2020[-dat.d,c(\"DISTRICT_NAME\", \"AGE\", \"SEX\")] # Remaining 30% test data\n",
    "\n",
    "train.dogs_labels <- dogs2020[dat.d,BREED]\n",
    "test.dogs_labels <-dogs2020[-dat.d,BREED]\n",
    "\n",
    "#Aprox k value\n",
    "k_value <-  sqrt(nrow(dogs2020))\n",
    "\n",
    "knn.test <- knn(train=train.dogs, test=test.dogs, cl=train.dogs_labels, k=k_value)\n",
    "\n",
    "ACC.test <- 100 * sum(test.dogs_labels == knn.test)/NROW(test.dogs_labels)\n",
    "\n",
    "table(knn.test ,test.dogs_labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Result of accuracy test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ACC.test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " #### Accuracy of direct prediction is ~11%. Result is better, but still not significantly improven."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## KNN with  Frequency Based Grouping\n",
    "\n",
    "Based on the previous KNN data we go on to gorup breeds with similar frequencies into groups inorder to produce more chunks with breeds that are more comparable in terms of frequency to avoid potential overpreferencing with KNN. It is a similar idea to the above Breed Type based grouping but with a different form. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculating values for grouping\n",
    "owners_sum <- sum(breed_filter$NUMBER_OWNERS)\n",
    "number_groups <- 6\n",
    "groupsize <- round(owners_sum / number_groups)\n",
    "\n",
    "# Adding grouping variable\n",
    "breed_filter[,group:=0]\n",
    "\n",
    "# In any situation first value is group 1\n",
    "breed_filter$group[1] <- 1\n",
    "\n",
    "# Assigning groups\n",
    "for (i in 2:nrow(breed_filter)) {\n",
    "  if (sum(breed_filter[group == breed_filter$group[i-1], NUMBER_OWNERS]) > groupsize){\n",
    "    breed_filter$group[i] <- breed_filter$group[i-1]+1\n",
    "  } else {\n",
    "    breed_filter$group[i] <- breed_filter$group[i-1]\n",
    "  }\n",
    "}\n",
    "breed_filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Preperation for KNN and Execution of KNN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transfer grouping to original dataset\n",
    "dogs2020 <- merge(dogs2020, breed_filter[,c(\"BREED\", \"group\")], by = \"BREED\", all.x = T)\n",
    "\n",
    "# Rewrite train and test datsets.\n",
    "# This is the only setting where we are not \"cheating\" by using dog associated attributes\n",
    "train.dogs <- dogs2020[dat.d, c(\"DISTRICT_NAME\", \"AGE\", \"SEX\", \"group\")] # 70% training data\n",
    "test.dogs <- dogs2020[-dat.d,c(\"DISTRICT_NAME\", \"AGE\", \"SEX\",\"group\")] # remaining 30% test data\n",
    "\n",
    "train.dogs_labels <- dogs2020[dat.d,BREED]\n",
    "test.dogs_labels <-dogs2020[-dat.d,BREED]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ####  Here we once more execute code very similar to the above part with breed_type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm(i)\n",
    "results <- 0\n",
    "for (i in 1:length(unique(breed_filter$group)))\n",
    "{\n",
    "  if (any(train.dogs$group == i) == F) {\n",
    "    next()\n",
    "  }\n",
    "  \n",
    "  train.dogs_labels <- dogs2020[dat.d,BREED]\n",
    "  test.dogs_labels <- dogs2020[-dat.d,BREED]\n",
    "  \n",
    "  knn.test <- knn(train=train.dogs[train.dogs$group == i], \n",
    "                  test=test.dogs[test.dogs$group == i], \n",
    "                  cl=train.dogs_labels[train.dogs$group == i], k=k_value)\n",
    "  \n",
    "  ACC.test <- 100 *\n",
    "    sum(test.dogs_labels[test.dogs$group == i] == knn.test)/\n",
    "    NROW(test.dogs_labels[test.dogs$group == i])\n",
    "  \n",
    "  table(knn.test ,test.dogs_labels[test.dogs$group == i])\n",
    "  \n",
    "  results[i] <- ACC.test\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculating weighted accuracy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(length(test.dogs_labels[test.dogs$group == 1]) * results[1]/100 \n",
    "  + length(test.dogs_labels[test.dogs$group == 2]) * results[2]/100\n",
    "  + length(test.dogs_labels[test.dogs$group == 3]) * results[3]/100\n",
    "  + length(test.dogs_labels[test.dogs$group == 4]) * results[4]/100\n",
    "  + length(test.dogs_labels[test.dogs$group == 5]) * results[5]/100\n",
    "  + length(test.dogs_labels[test.dogs$group == 6]) * results[6]/100) / \n",
    "    (length(test.dogs_labels[test.dogs$group == 1]) \n",
    "  + length(test.dogs_labels[test.dogs$group == 2])\n",
    "  + length(test.dogs_labels[test.dogs$group == 3])\n",
    "  + length(test.dogs_labels[test.dogs$group == 4])\n",
    "  + length(test.dogs_labels[test.dogs$group == 5])\n",
    "  + length(test.dogs_labels[test.dogs$group == 6])) * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The accuracy is ~29%, which is substantially better than 9% we got in the beginning. This is the best result we could get. However, the knn basically still follows the distribution of data. Because of that, one might consider invoking Occam's Razor, that is to say, perhaps a simpler more robust approach might lead to a comparable result with less complexity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "29F9YbZ-umjf"
   },
   "source": [
    "## Restore Backup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "eCjIyjrSumyX"
   },
   "outputs": [],
   "source": [
    "dogs2020 = copy(dogs2020_backup)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "E3zqXtrFv-f6"
   },
   "source": [
    "---\n",
    "# A Second Approach - Naive Bayes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "eubb1WX9xCjp"
   },
   "source": [
    "##  Who's the master anyway?\n",
    "Faced with some of the considerations of the KNN model, let's revisit our problem with a different approach.\n",
    "\n",
    "In this section we start again from square one with our basic question: Do any characteristics of Dog Owners influence their choice of Dog Breed, and if so, which ones?\n",
    "\n",
    "Using a Naive Bayes model we calculate frequency values for unique dog breeds against certain charateristics of the owners, specifically, age and gender. We also plot these frequency calculations on a per breed basis to better visualise the output of the model.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "VoJo7RGHxM2u"
   },
   "source": [
    "### Subset Data\n",
    "First we subset the main dataset into those specific characteristics of individual dogs and owners.\n",
    "We then remove less common breeds in order to present a tractable result.\n",
    "\n",
    "It is very important to note that this is the same numerical approach utilised while performing the 3rd KNN model above (namely Naive KNN with Breed Freq. Filtering). In doing so we illuminate some of the advantages and tradeoffs inherent in model selection, especially the necessary consideration that must be given between more sophisticated with inherent fragility and more traditional models.\n",
    "\n",
    "Finally we declare a helper variable of these common breeds to be used in the analysis.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "QKiNNNVRxPcY"
   },
   "outputs": [],
   "source": [
    "# Subset our Data\n",
    "dog_owner_chars <- subset(dogs2020, select=c(\"BREED\", \"BREED_TYPE\", \"YOB_DOG\", \"SEX_DOG\", \"COLOR_DOG\", \n",
    "                                             \"OWNER_ID\", \"AGE\", \"SEX\", \"WEALTH_T_CHF\", \"INCOME_T_CHF\", \n",
    "                                             \"DISTRICT_NAME\"))\n",
    "\n",
    "# Remove Outlier Breeds (<20 Entries)\n",
    "dog_owner_chars <- ddply(dog_owner_chars, \"BREED\", function(d) {if(nrow(d)>19) d else NULL})\n",
    "\n",
    "\n",
    "\n",
    "# A helper data frame with Unique Breeds\n",
    "unique_common_breeds <- unique(subset(dog_owner_chars, select=c(\"BREED\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Nskij5cBxQwd"
   },
   "source": [
    "### Naive Bayes Model\n",
    "Using the naiveBayes() function from the e1071 library we create a dataframe with relevant characteristics, where our independent variable is the breed, with our characteristics of owners acting as dependent variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "setDT(dog_owner_chars)\n",
    "\n",
    "dog_owner_chars <- subset(dog_owner_chars, select=c(\"BREED\", \"AGE\", \"SEX\"))\n",
    "\n",
    "# Naive Bayes Implementation\n",
    "nb <- naiveBayes(BREED ~ ., data=dog_owner_chars, laplace = 0, na.action = na.pass)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Printing the model to the output allows us to view a-priori probabilities as well as conditional probabilities of each breed on Age & Gender respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A quick glance at the data indicates that of these common breeds, the incidence of a male owner predicts a certain breed with a higher probability in only three cases! (English Bulldog, Rottweiler, Siberian Husky)\n",
    "\n",
    "It is important to remember our underlying distribution when considering this result. The preponderence of female dog owners in our underlying distribution means that it is very likely that the majority of breeds will be owned by a female owner, and indeed this is what we see. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Frequency vs Age\n",
    "We now use a simple for loop to construct histograms of the conditional probabilities of breed type on age.\n",
    "While most follow a \"normal\" distribution (placed in quotes as the distribution of population is ceratainly not completely normal, and there are other aspects of age that determine willingness to own a dog overall, for instance many small children are not able to care for a dog themselves, and similarly for the very old, a dog may be 'too much work'), some do not.\n",
    "\n",
    "A quick look at the data reveals a couple of highlights, to include that Border Collies seem to skew towards younger owners, while Collies tend towards older owners, while Rottweilers are almost dominated by owners in the 21-30 range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert nb into a data frame\n",
    "nb_df_age <- as.data.frame(nb$tables$AGE)\n",
    "\n",
    "for (i in 1:length(unique_common_breeds$BREED)) {\n",
    "  breed_name <- unique_common_breeds$BREED[i]\n",
    "  breed <- which(nb_df_age$Y == breed_name)\n",
    "  # create data frame for breed values\n",
    "  d_age <- data.frame(x = nb_df_age[breed,]$AGE, y = nb_df_age[breed,]$Freq)\n",
    "  # create plot\n",
    "  plot <- ggplot(d_age, aes(x = x, y = y, group = 1)) + geom_point() + geom_line() + \n",
    "    labs(x = \"Age group\", y = \"Frequency\", title = paste(breed_name, \"-- Frequency Across Age Groups\"))\n",
    "  print(plot)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Frequency vs Gender\n",
    "We now repeat the process for the Sex (or Gender) variable, this time utilising a bar graph.\n",
    "Astute readers will of course notice the predominance of female owners. This is likely due to the overall distribution of dog owners in general. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Similar approach for sex\n",
    "nb_df_sex <- as.data.frame(nb$tables$SEX)\n",
    "\n",
    "for (i in 1:length(unique_common_breeds$BREED)) {\n",
    "  breed_name <- unique_common_breeds$BREED[i]\n",
    "  breed <- which(nb_df_sex$Y == breed_name)\n",
    "  # Create data frame for sex values\n",
    "  # geom_bar() in ggplot2 takes all of the values that one wants to plot in the bar plot and automatically calculates teh frquency. \n",
    "  # Here the frequency is already given. \n",
    "  # A workaround is to create a data frame that replicates the frequency value with respect the number given for frequency, such that the right amount of values is output\n",
    "  rep_freq_f <- cbind(\"x\" = rep(nb_df_sex[breed,]$Freq[1], nb_df_sex[breed,]$Freq[1]*100),\n",
    "                      \"y\" = rep(\"f\", nb_df_sex[breed,]$Freq[1]*100))\n",
    "  rep_freq_m <- cbind(\"x\" = rep(nb_df_sex[breed,]$Freq[2], nb_df_sex[breed,]$Freq[2]*100),\n",
    "                      \"y\" = rep(\"m\", nb_df_sex[breed,]$Freq[2]*100))\n",
    "  freq <- as.data.frame(rbind(rep_freq_f, rep_freq_m))\n",
    "  # Create plot\n",
    "  plot <- ggplot(freq, aes(x = y)) + geom_bar() + \n",
    "    labs(x = \"Sex\", y = \"Frequency\", title = paste(breed_name, \"-- Frequency Across Genders\")) \n",
    "  print(plot)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "6ryLgVKdwWwm"
   },
   "source": [
    "---\n",
    "# Investigation Part 2 : Location, location, location?\n",
    "\n",
    "## An Analysis in Aggregate\n",
    "\n",
    "Let's now switch gears. We have performed some categorical analysis on the dataset to predict dog breed selection conditional on characteristics of owners. \n",
    "\n",
    "Now we turn our investigation to more general anaylsis of districts. For instance, we may ask ourselves: Does wealth or education influence preponderance of dog ownership? Does the proportion of single family homes vs. apartments play a predictive role? What characteristics of neighbourhoods predict dog ownership, if any? And what sort of policy levers may be used by policy leaders to promote or discourage dog ownership? (If they so choose of course!)\n",
    "\n",
    "To such matters we now turn our attention. We will begin by investigating the data in aggregate across district. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "NUK0HexcweRb"
   },
   "source": [
    "### Subset and Aggregate Data\n",
    "First, we subset the data into those characteristics that are relevant for district level analysis. We then take an aggregate of unique owner ids to produce a figure that represents the total number of dog owners in Zürich. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First we will need to Subset the Data and Take aggregates to get District Level Statistics\n",
    "# Relevant Columns w. Duplicates Removed\n",
    "district_dog <- unique(subset(dogs2020, select=c(\"OWNER_ID\", \"DISTRICT_NAME\", \"WEALTH_T_CHF\", \"INCOME_T_CHF\",\n",
    "                                                 \"BASIC_SCHOOL_PERCENTAGE\", \"GYMNASIUM_PERCENTAGE\", \"UNIVERSITY_PERCENTAGE\",\n",
    "                                                 \"SINGLE_FAMILY_HOMES\", \"INFRASTRUCTURE_BUILDINGS\", \"SMALL_BUILDINGS\", \"COMMERCIAL_BUILDINGS\",\n",
    "                                                 \"APARTMENTS\", \"FACTORIES_AND_WAREHOUSES\", \"SPECIAL_ACCOMODATION\", \"TOTAL_POPULATION\",\n",
    "                                                 \"FOREIGN_POPULATION_PERCENTAGE\")))\n",
    "\n",
    "# Aggregation - Count total unique owners and preserve select district data\n",
    "ddc <- aggregate(OWNER_ID ~ DISTRICT_NAME + WEALTH_T_CHF + INCOME_T_CHF + BASIC_SCHOOL_PERCENTAGE\n",
    "                 + GYMNASIUM_PERCENTAGE + UNIVERSITY_PERCENTAGE + SINGLE_FAMILY_HOMES\n",
    "                 + INFRASTRUCTURE_BUILDINGS + SMALL_BUILDINGS + COMMERCIAL_BUILDINGS\n",
    "                 + APARTMENTS + FACTORIES_AND_WAREHOUSES + SPECIAL_ACCOMODATION + TOTAL_POPULATION + FOREIGN_POPULATION_PERCENTAGE, \n",
    "                 data = district_dog, FUN = function(x){NROW(x)})\n",
    "\n",
    "# Rename Column(s) and Remove Redundant Data\n",
    "colnames(ddc)[16] <- \"TOTAL_UNIQUE_OWNERS\"\n",
    "rm(district_dog)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we compute some additional figures for the analysis that may be interesting. These additional attributes per district are described via the following: \n",
    "\n",
    "**NON_DOG_OWNERS** = By removing the amount of dog owners from the total population, we can compute a figure for non-dog owners <br>\n",
    "**PERCENT_DOG_OWNERS** = A simple ratio of dog owners to total population <br>\n",
    "**TOTAL_BUILDINGS** = We can combine all types of buildings to produce an aggregate <br>\n",
    "**TOTAL_RESIDENCES** = Similar to TOTAL_BUILDINGS, but for residential type buildings only. This is important to try and distinguish and isolate the effect of a type of building from the overall distribution of total buildings, on dog ownership. <br>\n",
    "\n",
    "It's important to notice that we have taken a percentage of each home type both over total buildings, as well as total residences. Again, this is to try and isolate the effect of a home type within residences from its frequency in the underlying distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute Difference\n",
    "ddc$NON_DOG_OWNERS <- ddc$TOTAL_POPULATION - ddc$TOTAL_UNIQUE_OWNERS\n",
    "ddc$PERCENT_DOG_OWNERS <- ddc$TOTAL_UNIQUE_OWNERS / ddc$TOTAL_POPULATION\n",
    "\n",
    "# Compute Total Building Count for Percentages % of Total Buildings and Residential Buildings\n",
    "ddc$TOTAL_BUILDINGS <- ddc$SINGLE_FAMILY_HOMES + ddc$INFRASTRUCTURE_BUILDINGS + ddc$SMALL_BUILDINGS\n",
    "                       + ddc$COMMERCIAL_BUILDINGS + ddc$APARTMENTS + ddc$FACTORIES_AND_WAREHOUSES + ddc$SPECIAL_ACCOMODATION\n",
    "ddc$TOTAL_RESIDENCES <- ddc$SINGLE_FAMILY_HOMES + ddc$APARTMENTS + ddc$SPECIAL_ACCOMODATION\n",
    "\n",
    "ddc$SINGLE_FAMILY_HOMES_PERCENTAGE <- ddc$SINGLE_FAMILY_HOMES / ddc$TOTAL_BUILDINGS\n",
    "ddc$SINGLE_FAMILY_HOMES_RESIDENCE_PERCENTAGE <- ddc$SINGLE_FAMILY_HOMES / ddc$TOTAL_RESIDENCES\n",
    "\n",
    "ddc$APARTMENTS_PERCENTAGE <- ddc$APARTMENTS / ddc$TOTAL_BUILDINGS\n",
    "ddc$APARTMENTS_RESIDENCE_PERCENTAGE <- ddc$APARTMENTS / ddc$TOTAL_RESIDENCES\n",
    "\n",
    "ddc$SPECIAL_ACCOMODATION_PERCENTAGE <- ddc$SPECIAL_ACCOMODATION / ddc$TOTAL_BUILDINGS\n",
    "ddc$SPECIAL_ACCOMODATION_RESIDENCE_PERCENTAGE <- ddc$SPECIAL_ACCOMODATION / ddc$TOTAL_RESIDENCES"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we display a summary of our table to begin our investigation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Summary Statistics\n",
    "summary(ddc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Model\n",
    "Let's continue to investigate. First, let's compute scatter plots of dog ownership percentage per district on individual district characteristics, including population, foreign population percentage, wealth & income, and various measures of building composition. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute Scatter Plots for Preliminary Investigation of Independent Variables\n",
    "scatter.smooth(x=ddc$TOTAL_POPULATION, y=ddc$PERCENT_DOG_OWNERS, main=\"Dog Ownership ~ Population\")\n",
    "scatter.smooth(x=ddc$FOREIGN_POPULATION_PERCENTAGE, y=ddc$PERCENT_DOG_OWNERS, main=\"Dog Ownership ~ Foreign Population %\")\n",
    "\n",
    "scatter.smooth(x=ddc$WEALTH_T_CHF, y=ddc$PERCENT_DOG_OWNERS, main=\"Dog Ownership ~ Wealth\")\n",
    "scatter.smooth(x=ddc$INCOME_T_CHF, y=ddc$PERCENT_DOG_OWNERS, main=\"Dog Ownership ~ Income\")\n",
    "\n",
    "scatter.smooth(x=ddc$SINGLE_FAMILY_HOMES_PERCENTAGE, y=ddc$PERCENT_DOG_OWNERS, main=\"Dog Ownership ~ Single Fam. Home %\")\n",
    "scatter.smooth(x=ddc$SINGLE_FAMILY_HOMES_RESIDENCE_PERCENTAGE, y=ddc$PERCENT_DOG_OWNERS, main=\"Dog Ownership ~ Single Fam. Home Res. %\")\n",
    "\n",
    "scatter.smooth(x=ddc$APARTMENTS_PERCENTAGE, y=ddc$PERCENT_DOG_OWNERS, main=\"Dog Ownership ~ Apartments %\")\n",
    "scatter.smooth(x=ddc$APARTMENTS_RESIDENCE_PERCENTAGE, y=ddc$PERCENT_DOG_OWNERS, main=\"Dog Ownership ~ Apartments Res. %\")\n",
    "\n",
    "scatter.smooth(x=ddc$SPECIAL_ACCOMODATION_PERCENTAGE, y=ddc$PERCENT_DOG_OWNERS, main=\"Dog Ownership ~ Spec. Acc. %\")\n",
    "scatter.smooth(x=ddc$SPECIAL_ACCOMODATION_RESIDENCE_PERCENTAGE, y=ddc$PERCENT_DOG_OWNERS, main=\"Dog Ownership ~ Spec. Acc. Res. %\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hmm. That's not very satisfying in the majority of the cases. We are seeing some pretty scattered distributions, with questionable correlations for the population and building composition characteristics. We should be very wary of trying to draw any conclusions from those features. \n",
    "\n",
    "Fortunately, our wealth and income graphs have produced plots with trends that could be viewed definitively as 'Significant'. How fortunate! Let's continue to investigate this path of our analysis.\n",
    "\n",
    "First let's compute the correlations for each of these two features with dog ownership percentage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Based upon the plots we decide to investigate the quality of wealth and income on Dog Ownership %\n",
    "cor_wealth <- cor(ddc$PERCENT_DOG_OWNERS, ddc$WEALTH_T_CHF)\n",
    "print('Wealth Corr: ')\n",
    "print(cor_wealth)\n",
    "cor_income <- cor(ddc$PERCENT_DOG_OWNERS, ddc$INCOME_T_CHF)\n",
    "print('Income Corr: ')\n",
    "print(cor_income)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Those are interesting correlations, suggesting that an increase in Wealth or Income is highly correlated with the amount of dog ownership for districts in aggregate. This is a very good foundation for a linear regression. Let's proceed. \n",
    "\n",
    "Let's keep in mind some of the fundamental assumptions that allow an ordinary least square (OLS) regression to function as the Best Linear Unbiased Estimator (BLUE), namely that the errors of each y term encompassed in the linear model are uncorrelated with one another, with mean 0 and homoscedastic finite variances. Additionally, there should be strong theoretical/foundational underpinning for a linear relationship.\n",
    "\n",
    "Let's proceed. \n",
    "\n",
    "We compute three regressions, two with a single variable, and one with multivariable input features of both wealth and income."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Two Independent Regressions\n",
    "linearModWealth <- lm(PERCENT_DOG_OWNERS ~ WEALTH_T_CHF, data=ddc)\n",
    "linearModIncome <- lm(PERCENT_DOG_OWNERS ~ INCOME_T_CHF, data=ddc)\n",
    "linearModCombi <- lm(PERCENT_DOG_OWNERS ~ WEALTH_T_CHF + INCOME_T_CHF, data=ddc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute summaries of each for an overview of our investigation. The computer performs all necessary caculations for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary(linearModWealth)\n",
    "summary(linearModIncome)\n",
    "summary(linearModCombi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Investigating the data, we are confronted with some facts.\n",
    "\n",
    "First let us speak to statistical significance. In both of the single-variable models, the predictors are statistically significant, with t-values small enough to trigger the categorical significance code '0', as shown by the outputs. \n",
    "\n",
    "With the 'combi' model, the t values explode into insignificant figures. How can this be? We will revisit this shortly. \n",
    "\n",
    "Now, in terms of significance of magnitude, it appears both Income and Wealth do indeed have an effect on the percentage of dog ownership in a district. This is great news! However, in terms of magnitude, only Income seems to have an effect that causes us to 'Raise our eyebrows'. \n",
    "\n",
    "In plain English: An additional 1000 CHF in average wealth for a district only predicts a 0.0055 % increase (please note this is percent) in dog ownership. It is stastistically significant, but not significant in magnitude. \n",
    "\n",
    "However and additional 1000 CHF in average income for a district predicts a 0.0269 % increase in dog ownership. Given the low magnitude of our numbers overall (percentage of dog ownership overall is in the low single digits), we can conclude that this is indeed a reliable result.\n",
    "\n",
    "So income of a district does indeed predict a higher rate of dog ownership, as stated: each 1000 CHF of additional average income predicts a 0.0269% increase in dog ownership.\n",
    "\n",
    "To revisit our final question, why then does the combined model break down? If both wealth and income effect dog ownership, shouldn't they both predict dog ownership when combined?\n",
    "\n",
    "The answer lies in a concept called multi-collinearity. In effect, it is almost certain that income and wealth are correlate <em>with each other</em> as well. When this happens, the variance in dog ownership cannot be reliably discerned from the variance in income and wealth independently, as they are so highly correlated with each other! This can be verified by the F-stastistic for the combined model. What this F statistic says (with a p value of 0.003158, so very significant statistically) is \"Does <em>some</em> variable in our idenpendent variables have an effect on our independent variable?\" A glance at the value (0.0003158) reveals that indeed, something is affecting the independent variable, dog ownership percentage.\n",
    "\n",
    "To think analogiously, we know that wealth or income predicts dog ownership, but we don't know which one (or both!), as they almost track each other!\n",
    "\n",
    "\n",
    "In our case we have made the sensible choice of income to represent both of these variables, which serves as an instrument for the overall depedence of dog ownership on financial condition. The income variable effectively 'encapsulates' the wealth variable in this case. \n",
    "\n",
    "Thus we conclude this project and our overall investigation with our most reliable and significant result: Dog Ownership is Conditional on Financial Condition, with Ownership Rates rising with increases in Income and Wealth. \n",
    "\n",
    "Thanks for a wonderful course and all the best.<br>\n",
    "Signed, \n",
    "\n",
    "Corey Bothwell<br/>\n",
    "Nicolleta Farabullini<br/>\n",
    "Jacob Gelling<br/>\n",
    "Andris Prokofjevs<br/>\n",
    "Qasim Warraich<br/>\n",
    "\n",
    "***"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Big Data Analytics.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}