1.1 How to use this document

This project was produced as part of the University of Pennsylvania’s Master of Urban Spatial Analytics Spring 2018 Practicum (MUSA 801) taught by Ken Steif, Michael Fichman, and Matt Harris. We would like to thank Jonathan Pyle of Philadelphia Legal Assistance for providing feedback and data.

The following document presents an analysis of illegal rentals in Philadelphia and an interactive tool for the Philadelphia Department of Licenses and Inspections (L&I) that will ultimately help proactively prioritize inspections to stymy illegal rentals before they turn into illegal evictions. The tool will provide a risk score for every rental property in the city, indicating each property’s likelihood of being rented illegally. This document includes a detailed overview of the use case, analysis and key findings, as well as the relevant code needed to replicate the analysis.

1.2 Abstract

In Philadelphia, all rental units are required to be certified prior to leasing. Through certification, landlords must obtain a rental license, at a cost of about $55 annually. If a unit is not certified, but is still rented out, it is considered an illegal rental. If a landlord is renting illegally, then any subsequent tenant eviction is also illegal, meaning the tenant has no legal recourse.

With no official record of the rental and the eviction, it is difficult to distill how often landlords are renting out their properties illegally. The only record we have for illegal rentals is the Department of Licenses and Inspections (L&I) Code Violation 9-3902; renting without a license, but this number is far underestimating the problem.

Currently, L&I conducts inspections reactively, according to received 311 complaints. However, L&I only employs 17 inspectors, and there is a backlog of complaints still to be inspected. A statistical model that can predict property-specific illegal rental risk would allow L&I to proactively address illegal rentals before they turn into illegal evictions. This tool predicts the likelihood that a property is being rented to a tenant without a rental license for all non-owner occupied homes in Philadelphia to help L&I prioritize inspections more efficiently.

1.3 Motivation

Between 2005 to 2016, rental prices in Philadelphia went up by 10 percent, while wages only went up by 0.3 percent. As a result, over half of Philadelphia renters are rent-burdened (spending over one-third of income on housing), and nearly one-third of these renters spend over 50 percent of their salary on housing.

Figure 3. The universe of potential rental properties

With our data coming from many sources, the properties in one dataset were not always included in another. As a result, when we joined the different datasets together, these properties, or rows, appear to be missing data. For example, if a property from the OPA dataset does not match with any properties from the violations dataset, based off of its OPA Account Number, we can assume that this property does not have a violation. It is important to note that any assumption will inevitably bias our data. Because across all of our datasets we assumed that missing data meant the absence of a matching property, we bias our data towards underestimating violations, licenses, and tax delinquencies.

Figure 4. Table of properties from each dataset

2.1 Data Stories and Feature Engineering

We sorted our data into four different “stories” that together paint a full picture of illegal rentals. These four stories are used to estimate an illegal rental risk score for each property city wide. This is done to ensure our model is sensitive to a diverse range of variables.

The dependent variable for our analysis - that is the variable we are trying to predict - is illegal rental violations, which we encode in our data as ‘0’ to represent properties that have not received an illegal rental violation and ‘1’ to represent properties that have received an illegal rental violation.

3.1 Where are evictions in Philadelphia, and how has that changed over time?

Understanding where legal evictions are occurring is beneficial for understanding the nature of the rental market in Philadelphia, which ultimately provides insight into illegal rentals.

Between 2016 and 2018, there were about 55,700 evictions in the city. This represents about one fifth of all rental properties in the city. As the below chart shows, this number peaked in 2017.

3.2 Where are illegal rentals in Philadelphia and what types of properties are being illegally rented?

Between 2016-2018, there were around 7,500 citations given out by L&I inspectors to individuals renting their property illegally. It is worth noting that this figure is almost certainly under-reported given the limited inspection resources L&I have at their disposal.

Between 2016-2018, citations were given out for illegally renting a single-family home, illegally renting a multi-family home, and illegally renting a room, with the distribution heavily favoring single family homes.

3.3 What is the relationship between key variables and illegal rentals?

Comparing homes with and without illegal rentals provides a strong basis for understanding the relationship between illegal rentals and explanatory variables. This allows us to see which variables see higher instances of illegal rental violations.

As the below charts show, illegal rentals are more prevalent in properties that have had interior and exterior violations, are tax delinquent, have out of state landlords, and are located on parcels with multiple units. However, they are less prevalent in properties built before 1900.

4.1 Model Building

Illegal rentals are a rare occurrence, comprising less than 0.1 percent of our final dataset. This is in part because illegal rentals are likely undercounted. As such, our models for the entire dataset predicted that almost all properties would not be illegal rentals, which makes sense given the rarity of illegal rentals. To create a more balanced dataset, we randomly sampled our data, so the full dataset had four times as many non-violations than violations.

To build our predictive model we estimated separate models for each of the four stories and then ensembled them into one final predictive model. To determine what variables were most significant, we used stepwise selection. The stepwise model added and deleted variables one at a time so that we could see which variables are significant, but also how the variables correlate with each other. Through this process, we were able to winnow down our over 200 independent variables into 24 variables spread across four stories.

4.2 Validation

The predictive model provides an illegal rental risk score for all rental properties in Philadelphia. This risk score is used to classify properties as ‘worthy of inspection’ or ‘not worthy of inspection’ This ‘worthiness score’ provides a basis for measuring the accuracy and generalizability of our model.

While our continuous relative risk score is useful for L&I, we set a threshold to determine specified outcomes, listed in the table below. These outcomes allow for a more direct comparison of costs between our data driven approach and L&I’s business as usual.

4.3 Cost Benefit Analysis

Ultimately, one of the key decision points in deciding between a data-driven approach and business as usual is the social and economic costs of both approaches. In other words, does our model help L&I use their resources more efficiently and does it reduce the high social cost of an illegal rental violation and an illegal eviction? To determine this we performed a cost benefit analysis, where we calculated a back of the envelope cost for all outcomes and compared this to the current cost of doing business for 1,000 properties.

We estimated the following costs:

• $1,180 for every true positive
• $50 for every true negative
• $50 for every false positive
• An unknown, but high social cost for every family displaced by a false negative

7.1 Data Wrangling, Feature Engineering, and Modeling

First we set up our workspace.

#load libraries

#load libraries
require(ggthemes)
require(IDPmisc)
library(corrplot)
library(stargazer)
library(sjPlot)
library(data.table)
library(ggplot2)
library(ggmap)
library(maps)
library(sp)
library(raster)
library(lattice)
library(latticeExtra)
library(RColorBrewer)
library(rasterVis)
library(dplyr)
library(stringr)
library(tidyverse)
library(tidyr)
library(stats)
library(data.table)
library(grid)
library(gridExtra)
library(rgdal)
library(sf)
library(lattice)
library(caret)
library(tidycensus)
library(FNN)
library(Rcpp)
library(viridisLite)
library(viridis)
library(jsonlite)
library(stringi)
library(downloader)
library(lubridate)
library(spdep)
library(randomForest)
library(kableExtra)


#remove scientific noation
options(scipen=999)

Next we began reading in our datasets, cleaning each dataset, and engineering features.

##### OPA DATA #####

#read data

opa_orig<-fread('opa_properties_public.csv', header=TRUE, sep=",", colClasses='numeric')

#make copy
opafull<-opa_orig

#look for empties and NA's
colSums(opafull=='')
colSums(is.na(opafull))

#first cleaning of unecessary columns
opafull <- subset(opa_orig,select=-c(assessment_date,beginning_point,book_and_page,building_code_description,
                                  category_code_description, cross_reference, geographic_ward,
                                  house_number,market_value_date,objectid, parcel_shape, recording_date,
                                  registry_number, street_designation,street_direction,suffix,year_built_estimate,
                                  date_exterior_condition, mailing_address_1, mailing_address_2,mailing_care_of,
                                  building_code,category_code, owner_2,street_code))

sort(names(opafull))

#distribution of NA values across the datatset
colSums(is.na(opafull))

#distribution of empty values across the dataset
colSums(opafull=='')


##Narrow down opa dataset to possible rental houses

#by zoning code (those that begin with R)
opa1<- opafull %>% filter(str_detect(zoning, '^R'))

#by homestead exempt homes
sum(is.na(opa1$homestead_exemption)) #1270 NA
sum(opa1$homestead_exemption=='') #NA empties
#0 not exempt; 1 is exempt(=homeowners)
table(opa1$homestead_exemption)
#filter by 0, which are the not exempt units; 
#consider NA's which we do not know, but assume is a non-exempt/rental
opa1 <- opa1 %>% filter(opa1$homestead_exemption == 0 | is.na(opa1$homestead_exemption))
table(opa1$homestead_exemption)
sum(is.na(opa1$homestead_exemption))

#OPA dataset cut down from 580k to 510k to 305k


#####################################################################################

#wrangling variables with empties

#basements
table(opa1$basements)

#create variables for basement types
opa1$basement_dummy <- ifelse(opa1$basements=='ZZ',1,0)
opa1$basement_none <- ifelse(opa1$basements=='No_Basement',1,0)
opa1$basement_full <- ifelse(opa1$basements == 'Full' , 1,0)
opa1$basement_partial <- ifelse(opa1$basements == 'Partial', 1,0)
opa1$basement_sizeNotKnown <- ifelse(opa1$basements == 'Size_Unknown', 1,0)

opa1$basement_dummy <- as.factor(opa1$basement_dummy)
opa1$basement_none <- as.factor(opa1$basement_none)
opa1$basement_full <- as.factor(opa1$basement_full)
opa1$basement_partial <- as.factor(opa1$basement_partial)
opa1$basement_sizeNotKnown <- as.factor(opa1$basement_sizeNotKnown)


#central_air
table(opa1$central_air)

#create variables for central air
opa1$central_air_dummy <- ifelse(opa1$central_air_new=='ZZ',1,0)
opa1$central_air_y <- ifelse(opa1$central_air_new==1,1,0)
opa1$central_air_n <- ifelse(opa1$central_air_new==0,1,0)

opa1$central_air_dummy <- as.factor(opa1$central_air_dummy)
opa1$central_air_y <- as.factor(opa1$central_air_y)
opa1$central_air_n <- as.factor(opa1$central_air_n)

#unfinished (u = unfinished/buildings under construction)
table(opa1$unfinished)

#create variables for unfinished buildings
opa1$unfinished_unit<-ifelse(opa1$unfinished=='',0,1)
opa1$unfinished_unit<-as.factor(opa1$unfinished_unit)


#unit (specific unit in condo, building complex, etc, where there are subproperties)
table(opa1$unit)

#create variable for unit type
#1-is a subunit; 0-not a subunit
opa1$unit_new <- ifelse(opa1$unit=='',0,1)
opa1$unit_new <- as.factor(opa1$unit_new)

#utility
table(opa1$utility)

#create variables for utility type
#1 - remodeled/modernized; 0 - no remodeling
opa1$utility_new <- ifelse(opa1$utility_new=='None',0,1)
opa1$utility_new <- as.factor(opa1$utility_new)

#garage type
table(opa1$garage_type)

#create variables for garage type
opa1$garage_type_dummy <- ifelse(opa1$garage_type_new=='ZZ',1,0)
opa1$garage_type_y <- ifelse(opa1$garage_type_new=='Has_Garage',1,0)
opa1$garage_type_n <- ifelse(opa1$garage_type_new=='None',1,0)

opa1$garage_type_dummy <- as.factor(opa1$garage_type_dummy)
opa1$garage_type_y <- as.factor(opa1$garage_type_y)
opa1$garage_type_n <- as.factor(opa1$garage_type_n)


#general construction type
table(opa1$general_construction)

#create variable for construction type
opa1$general_construction_dummy<-ifelse(opa1$general_construction_new=='ZZ',1,0)
opa1$general_construction_brickY<-ifelse(opa1$general_construction_new=='Common_Brick',1,0)
opa1$general_construction_brickN<-ifelse(opa1$general_construction_new=='Not_Common_Brick',1,0)

opa1$general_construction_dummy <- as.factor(opa1$general_construction_dummy)
opa1$general_construction_brickY <- as.factor(opa1$general_construction_brickY)
opa1$general_construction_brickN <- as.factor(opa1$general_construction_brickN)


#house_extension
table(opa1$house_extension)

#create variable for house extesnion
opa1$house_extension_n<-ifelse(opa1$house_extension_new=='No_Extension',1,0)
opa1$house_extension_n <- as.factor(opa1$house_extension_n)


#heater type
table(opa1$type_heater)

#heater Yes/No
opa1$type_heater_yn <- ifelse(opa1$type_heater=='','ZZ', 
                              ifelse(opa1$type_heater=='0','NoHeater','HasHeater'))


#other building
table(opa1$other_building)

#create variable for other building
opa1$other_bldg_yn <- ifelse(opa1$other_building=='Y',1,0)
opa1$other_bldg_yn <- as.factor(opa1$other_bldg_yn)

#wrangling variables with NA's from the OPA dataset

#remove one row with several NA columns throughout - this should remove only one row
opa1 <- opa1[!is.na(opa1$number_of_bedrooms),]

#check
colSums(is.na(opa1))
colSums(opa1=='')


#exempt building; assuming NA is not exempt
table(opa1$exempt_building)

#make variable for exempt building
opa1$exempt_building_dummy <- ifelse(is.na(opa1$exempt_building),1,0)
opa1$exempt_building_y <- ifelse(is.na(opa1$exempt_building),0, ifelse(opa1$exempt_building==0,0,1))

opa1$exempt_building_dummy <- as.factor(opa1$exempt_building_dummy)
opa1$exempt_building_y <- as.factor(opa1$exempt_building_y)


#exempt land; assuming NA is not exempt
table(opa1$exempt_land)

#make variable for exempt land
opa1$exempt_land_dummy <- ifelse(is.na(opa1$exempt_land),1,0)
opa1$exempt_land_y <- ifelse(is.na(opa1$exempt_land),0, ifelse(opa1$exempt_land==0,0,1))

opa1$exempt_land_dummy<-as.factor(opa1$exempt_land_dummy)
opa1$exempt_land_y<-as.factor(opa1$exempt_land_y)


#taxable building; assuming NA is not taxable
table(opa1$taxable_building)

#make variable for taxable building
opa1$taxable_building_dummy <- ifelse(is.na(opa1$taxable_building),1,0)
opa1$taxable_building_y <- ifelse(is.na(opa1$taxable_building),0, ifelse(opa1$taxable_building==0,0,1))

opa1$taxable_building_dummy<-as.factor(opa1$taxable_building_dummy)
opa1$taxable_building_y<-as.factor(opa1$taxable_building_y)


#taxable land; assuming NA is not taxable
table(opa1$taxable_land)

#make variable for taxable land
opa1$taxable_land_dummy <- ifelse(is.na(opa1$taxable_land),1,0)
opa1$taxable_land_y <- ifelse(is.na(opa1$taxable_land),0, ifelse(opa1$taxable_land==0,0,1))

opa1$taxable_land_dummy<-as.factor(opa1$taxable_land_dummy)
opa1$taxable_land_y<-as.factor(opa1$taxable_land_y)

#number of stories
table(opa1$number_stories)

#make variables for number of stories
opa1$story1<-ifelse(opa1$number_stories == "1", 1,0)
opa1$story2<-ifelse(opa1$number_stories == "2", 1,0)
opa1$story3<-ifelse(opa1$number_stories == "3", 1,0)
opa1$story4plus<-ifelse(opa1$number_stories >= "4", 1,0)

opa1$story1<-as.factor(opa1$story1)
opa1$story2<-as.factor(opa1$story2)
opa1$story3<-as.factor(opa1$story3)
opa1$story4plus<-as.factor(opa1$story4plus)

#year built; **remove NA's later (only 3)
class(opa1$year_built)
opa1$year_built_copy<-as.numeric(as.character(opa1$year_built))

#make variables for year built
opa1$built_pre1900<-ifelse(opa1$year_built_copy <= 1900, 1,0) 
opa1$built_btw1900and1950<-ifelse(opa1$year_built_copy > 1900 & opa1$year_built_copy < 1950, 1,0)
opa1$built_post1950<-ifelse(opa1$year_built_copy >= 1950, 1,0)

opa1$built_pre1900<-as.factor(opa1$built_pre1900)
opa1$built_btw1900and1950<-as.factor(opa1$built_btw1900and1950)
opa1$built_post1950<-as.factor(opa1$built_post1950)

#exterior condition; aggregate into avg/not avg;
#assumes NA are not avg (21,814 NA's)
table(opa1$exterior_condition)
sum(is.na(opa1$exterior_condition))

opa1$exterior_condition_new <- opa1$exterior_condition
opa1$exterior_condition_new <- ifelse(is.na(opa1$exterior_condition_new),999,opa1$exterior_condition_new)
opa1$exterior_condition_new <- as.factor(opa1$exterior_condition_new)

#make variables for exterior condition
opa1$extcond_dummy <- ifelse(opa1$exterior_condition_new==999,1,0)
opa1$extcond_avg <- ifelse(opa1$exterior_condition_new==2 | 
                           opa1$exterior_condition_new==3 | 
                           opa1$exterior_condition_new==4, 1,0)

opa1$extcond_dummy <- as.factor(opa1$extcond_dummy)
opa1$extcond_avg <- as.factor(opa1$extcond_avg)


#interior condition
#assumes NA are not avg (22,059 NA's)
table(opa1$interior_condition)
sum(is.na(opa1$interior_condition))

opa1$interior_condition_new <- opa1$interior_condition
opa1$interior_condition_new <- ifelse(is.na(opa1$interior_condition_new),999,opa1$interior_condition_new)
opa1$interior_condition_new <- as.factor(opa1$interior_condition_new)

#make variables for interior condition
opa1$intcond_dummy <- ifelse(opa1$interior_condition_new==999,1,0)
opa1$intcond_avg <- ifelse(opa1$interior_condition_new==2 | 
                           opa1$interior_condition_new==3 | 
                           opa1$exterior_condition_new==4, 1,0)

opa1$intcond_dummy <- as.factor(opa1$intcond_dummy)
opa1$intcond_avg <- as.factor(opa1$intcond_avg)

#garage spaces
class(opa1$garage_spaces)
table(opa1$garage_spaces)

#subdivide garage spaces
opa1$garage_spaces_temp <-opa1$garage_spaces
opa1$garage_spaces_temp <- ifelse(is.na(opa1$garage_spaces_temp),999,
                                  ifelse(opa1$garage_spaces_temp>=4,4,opa1$garage_spaces_temp))

#make variables for garage spaces
opa1$garagespace0<-ifelse(opa1$garage_spaces_temp == 0, 1,0)
opa1$garagespace1<-ifelse(opa1$garage_spaces_temp == 1, 1,0)
opa1$garagespace2<-ifelse(opa1$garage_spaces_temp == 2, 1,0)
opa1$garagespace3<-ifelse(opa1$garage_spaces_temp == 3, 1,0)
opa1$garagespace4plus<-ifelse(opa1$garage_spaces_temp >= 4, 1,0)

opa1$garage_spaces_temp <- as.factor(opa1$garage_spaces_temp)
opa1$garagespace0 <- as.factor(opa1$garagespace0)
opa1$garagespace1 <- as.factor(opa1$garagespace1)
opa1$garagespace2 <- as.factor(opa1$garagespace2)
opa1$garagespace3 <- as.factor(opa1$garagespace3)
opa1$garagespace4plus <- as.factor(opa1$garagespace4plus)

#1-there is garage; 0-no garage, assuming NA is no garage
opa1$garage_yn <- ifelse(is.na(opa1$garage_spaces) | opa1$garage_spaces==0,0,1)
opa1$garage_yn <- as.factor(opa1$garage_yn)

table(opa1$garage_yn)

#Miscellaneous columns

opa1$depth_copy <- as.numeric(opa1$depth)
opa1$depth_dummy <- ifelse(opa1$depth==0,1,0)
opa1$depth_dummy <- as.factor(opa1$depth_dummy)

opa1$frontage_copy <- opa1$frontage
opa1$frontage_dummy <- ifelse(opa1$frontage==0,1,0)
opa1$frontage_dummy <- as.factor(opa1$frontage_dummy)

opa1$total_area_copy<-opa1$total_area
opa1$total_area_dummy<-ifelse(opa1$total_area==0,1,0)
opa1$total_area_dummy <- as.factor(opa1$total_area_dummy)

opa1$total_livable_area_copy <- as.numeric(opa1$total_livable_area)
opa1$total_livable_area_dummy <- ifelse(opa1$total_livable_area==0,1,0)
opa1$total_livable_area_dummy <- as.factor(opa1$total_livable_area_dummy)

opa1$market_value_copy <- opa1$market_value
opa1$market_value_copy <- ifelse(is.na(opa1$market_value_copy),0,opa1$market_value_copy)
opa1$market_value_dummy <- ifelse(opa1$market_value_copy==0,1,0) 
opa1$market_value_dummy <- as.factor(opa1$market_value_dummy)


#MORE OPA FEATURE ENGINEERING

#Create absentee landlord variable (defined as outside of Philadelphia)
sum(is.na(opa1$mailing_city_state))
sum(opa1$mailing_city_state=='') #138666 missing

#0 owner is in city; 1 owner is not in city Philadelphia (absent landlord)
opa1$absenteelandlord<-ifelse(str_detect(opa1$mailing_city_state,"^PHILA"), 0, 
                             ifelse(str_detect(opa1$mailing_city_state,"^PHILD"), 0, 
                                    ifelse(str_detect(opa1$mailing_city_state,"^PHILK"), 0,1)))
opa1$absenteelandlord <- as.factor(opa1$absenteelandlord)
table(opa1$absenteelandlord)


#further break down into people that are not in Pennsylvania or New Jersey
#1-out of state; 0-in state
opa1$outofstateLL<-ifelse(grepl("PA", opa1$mailing_city_state), 0, ifelse(grepl("NJ", opa1$mailing_city_state), 0, 1))
#1-out of philly, 0-in philly
opa1$outofPhillyLL<-ifelse(grepl("PHILA", opa1$mailing_city_state), 0,1)

#make into factors
opa1$outofstateLL<-as.factor(opa1$outofstateLL)
opa1$outofPhillyLL<-as.factor(opa1$outofPhillyLL)


#check final opa dataset
colSums(is.na(opa1))
colSums(opa1=='')



#subset with usefull columns for final data modeling

#by removing edited/feature engineered columns
opa_clean <- subset(opa1,select=-c(basements,central_air,garage_type,unfinished,unit,utility,general_construction,
                                   house_extension,site_type,topography,type_heater,other_building,view_type,
                                   exempt_building,exempt_land,fireplaces,taxable_building,taxable_land,
                                   year_built,depth,frontage,number_of_bedrooms, number_of_rooms,number_stories,
                                   number_of_bathrooms,quality_grade,exterior_condition,interior_condition,
                                   garage_spaces,homestead_exemption,fuel,separate_utilities,sewer,
                                   mailing_street,market_value))

colSums(opa_clean=='')
colSums(is.na(opa_clean))

opa_clean$parcel_number_copy <- as.numeric(opa1$parcel_number)

#copy for knn later
opa_final<-opa_clean


##### VIOLATIONS DATA #####


#read data
violations<-read.csv('li_violations.csv')

#make copy
violations2 <- violations

#check empties and NA
colSums(is.na(violations2))
colSums(violations2=='')

#first cleaning of unecessary columns
violations2<-subset(violations,select=-c(geocode_x,geocode_y,objectid,casenumber,apfailkey,
                                         the_geom_webmercator))

unlist(lapply(violations2,class))

#limit violations to 2016/17/18
violations2$year<-substring(violations2$violationdate,1,4)
violations_clean<-subset(violations2, violations2$year=="2016" | violations2$year=='2017' | violations2$year=='2018')

#cuts down violations 1,319,904 to 340,684

#create binary for illegal rental -
#### 9-3902(1) = 1 and 2 family homes, 9-3902(2) = Multi-family homes, 9-3902(3) = rooming unit or house, 9-3902(4) = hotel
violations_clean$illegalRental <- ifelse(violations_clean$violationtype == '9-3902 (1)' | 
                                           violations_clean$violationtype == '9-3902 (2)' |
                                           violations_clean$violationtype == '9-3902 (3)',1,0)

#create binary for single family illegal rentals
violations_clean$SFIllegalRent<-ifelse(violations_clean$violationtype == '9-3902 (1)', 1,0)
violations_clean$SFIllegalRent<-as.factor(violations_clean$SFIllegalRent)

#create binary for multi-family illegal rentals
violations_clean$MFIllegalRent<-ifelse(violations_clean$violationtype == '9-3902 (2)', 1,0)
violations_clean$MFIllegalRent<-as.factor(violations_clean$MFIllegalRent)

#create binary for rooming unit or house illegal rental
violations_clean$RoomIllegalRent<-ifelse(violations_clean$violationtype == '9-3902 (3)', 1,0)
violations_clean$RoomIllegalRent<-as.factor(violations_clean$RoomIllegalRent)

#violation status
#table(violations_clean$status)
#table(violations_clean$casestatus)
violations_clean$issue_fixed<-ifelse(violations_clean$status=="COMPLIED",1, ifelse(violations_clean$status=="RESOLVE", 1,0))
violations_clean$active_case<-ifelse(violations_clean$casestatus=="CLOSED", 0, 1)

violations_clean$issue_fixed<-as.factor(violations_clean$issue_fixed)
violations_clean$active_case<-as.factor(violations_clean$active_case)


#Number of Previous Violations

violation_count<-violations_clean%>% 
  group_by(opa_account_num)%>%
  summarize(viol_count=n())

violation_count_other<-violations_clean %>%
  filter(violationtype != '9-3902 (1)')%>%
  filter(violationtype != '9-3902 (2)')%>%
  filter(violationtype != '9-3902 (3)')%>%
  group_by(opa_account_num) %>%
  summarize(other_violations=n())

#violation_count_other<-violation_count_other[!is.na(violation_count_other$opa_account_num),]

violations_clean$opa_account_num_copy<-as.numeric(as.character(violations_clean$opa_account_num))
violation_count$opa_account_num_copy<-as.numeric(as.character(violation_count$opa_account_num))
violation_count_other$opa_account_num_copy<-as.numeric(as.character(violation_count_other$opa_account_num))

violations_clean <- left_join(violations_clean,violation_count, by='opa_account_num_copy')
violations_clean <- left_join(violations_clean,violation_count_other, by='opa_account_num_copy')

violations_clean$other_violations<-as.numeric(violations_clean$other_violations)
violations_clean$other_violations<-ifelse(!is.na(violations_clean$other_violations),
                                          violations_clean$other_violations, 0)


#Landlord with previous violations
landlord_count<-violations_clean%>% 
  group_by(ownername)%>%
  summarize(owner_count=n())
landlord_count <- landlord_count[-c(1), ]

violations_clean <- left_join(violations_clean,landlord_count, by='ownername')

#assumes NA (no owner name on record) is 1 owner only
violations_clean$owner_count<-ifelse(is.na(violations_clean$owner_count),1,violations_clean$owner_count)
violations_clean$owner_count_dummy<-ifelse(is.na(violations_clean$owner_count),1,0)
violations_clean$owner_count_dummy<-as.factor(violations_clean$owner_count_dummy)

#Type of violation
#table(violations_clean$violationdescription)
violations_clean$exterior<-ifelse(grepl("EXT", violations_clean$violationdescription), 1,0)
violations_clean$interior<-ifelse(grepl("INT", violations_clean$violationdescription), 1,0)


#Length of time since last violation (based on Jan 2019 date)

#make copy
violations_clean2<-violations_clean
#calculate time since last violation
violations_clean2$now<-as.Date('2019-01-31')
violations_clean2$violationsdate2<-as.Date(substring(violations_clean2$violationdate,1,10))

violations_clean3 <-
  violations_clean2 %>%
  mutate(now2 = ymd(now),
         then2 = ymd(violationsdate2))

#time length in days
violations_clean3$lentim_sincelast_viol <- violations_clean3$now2-violations_clean3$then2
class(violations_clean3$lentim_sincelast_viol)
violations_clean3$lentim_sincelast_viol <- as.numeric(as.character(violations_clean3$lentim_sincelast_vio))

#sum exterior, interior and illegalRental by opaaccount
violations_summed <-
  violations_clean3 %>%
  group_by(opa_account_num_copy) %>%
  summarise(had_extvio = sum(exterior),
            had_intvio = sum(interior),
            had_illegalRent = sum(illegalRental))
violations_summed <- violations_summed[!is.na(violations_summed$opa_account_num_copy), ]


violations_clean4<-left_join(violations_clean3,violations_summed, by='opa_account_num_copy')

#history of illegal rental and interior and exterior violations
#could be biased if a property was never cited for illegal rental

#hist_ columns as binary
#1-had prior history of illegal rental; 0-no history of illegal rental
violations_clean4$hist_illegalRent<-ifelse(violations_clean4$had_illegalRent==0,0,1)
#1-had prior history of external violation; 0-no history of external violation
violations_clean4$hist_extviol <- ifelse(violations_clean4$had_extvio==0,0,1)
#1-had prior history of internal violation; 0-no history of internal violation
violations_clean4$hist_intviol <- ifelse(violations_clean4$had_intvio==0,0,1)


#make numeric columns as factors since binary
violations_clean4$illegalRental<-as.factor(violations_clean4$illegalRental)
violations_clean4$exterior<-as.factor(violations_clean4$exterior)
violations_clean4$interior<-as.factor(violations_clean4$interior)
violations_clean4$hist_illegalRent<-as.factor(violations_clean4$hist_illegalRent)
violations_clean4$hist_extviol<-as.factor(violations_clean4$hist_extviol)
violations_clean4$hist_intviol<-as.factor(violations_clean4$hist_intviol)

#make binary for multiple units
#update status
violations_clean4$status2 <- sub("^$", "OPEN", violations_clean4$status)

violations_clean4$multipleunits<-ifelse(violations_clean4$unit=="",0,1)

##make binary for has had violation (0 if only 1 violation; +1 if has had previous violations)
violations_clean4$had_anyviolation<-ifelse(violations_clean4$other_violations >0, 1,0)
table(violations_clean4$had_anyviolation)


#remove unecessary columns
sort(names(violations_clean4))
colSums(is.na(violations_clean4))


vio_clean_full<-subset(violations_clean4, select=-c(the_geom, addresskey, zip, censustract, 
                                                    aptype, casegroup, casepriority, now, now2, 
                                                    then2, violationsdate2,prioritydesc,
                                                    caseresolutioncode,casestatus, status,unit))

#make column class accordingly
unlist(lapply(vio_clean_full,class))

vio_clean_full$year<-as.factor(vio_clean_full$year)
vio_clean_full$status2<-as.factor(vio_clean_full$status2)
vio_clean_full$had_anyviolation<-as.factor(vio_clean_full$had_anyviolation) #not really needed
vio_clean_full$multipleunits<-as.factor(vio_clean_full$multipleunits)

vio_clean_full$mostrecentinsp2<-as.Date(vio_clean_full$mostrecentinsp)

#keep only one record of unique opa;
violations_clean_unique <-
  vio_clean_full %>%
  group_by(opa_account_num.x) %>%
  slice(which.max(mostrecentinsp2))

names(violations_clean_unique)
colSums(is.na(violations_clean_unique)) #1 NA - can remove later
#violations_clean_unique<- violations_clean_unique[!is.na(violations_clean_unique$opa_account_num),]

#copy for knn later
violations_final <- violations_clean_unique


##############################################################################################

##### BUSINESS LICENSES #####

#read data
business_orig<- st_read('li_business_licenses.shp')

#make copy
business<-business_orig
sort(names(business))


#create ACTIVE / INACTIVE column
business$ActiveInactive <- ifelse(business$licensesta=='Active','Active','Inactive,Expired,Closed,Revoked')

#subset RENTALS / VACANT RESIDENTIALS / HIGH RISE
business2 <- subset(business, business$licensetyp == 'Rental' | business$licensetyp == 'High Rise')

#make LAT and LNG coordinates from geometry
business2$lng <- st_coordinates(business2)[,1]
business2$lat <- st_coordinates(business2)[,2]

#check type
class(business2) #sf;data.frame

#check for NA's
sum(is.na(business2$opaaccount)) #4237 NA's
#cut out NA's
business2<-business2[!is.na(business2$opaaccount),]

class(business2$opaaccount) #factor
business2$opaaccount_Copy <- as.numeric(as.character(business2$opaaccount))
class(business2$opaaccount_Copy) #numeric

class(business2$fulladdres) #factor
business2$fulladdres_Copy <- as.character(business2$fulladdres)
class(business2$fulladdres_Copy) #character

class(business2) #sf; data.frame
business2_df <- as.data.frame(business2)
class(business2_df) #data.frame

table(droplevels(business2_df$licensetyp)) # 207596 rentals; 613 high rise
table(business2_df$ActiveInactive)         # 883639 active; 124570 inactive (~60% inactive)
table(business2_df$legalentit)             # 66433 company; 141577 individual

#feature engineering
#Number of previous licenses
business_temp <- 
  business2_df %>%
  group_by(fulladdres) %>%
  summarise(num_prev_lic=n())

business2_df<- left_join(business2_df,business_temp, by='fulladdres')

#Active or inactive based on most recent license date
business3 <- business2_df
business3$mostrecent2<- as.Date(as.character(business3$mostrecent))
business3$status_mostrecent <- ifelse (business3$mostrecent2 < as.Date('2018-01-31'), 'inactive','active')

#Length of time since last license if most recent is inactive
business3$license_year <- format(as.Date(business3$mostrecent, format="%d/%m/%Y"),"%Y")

class(business3$license_year) #character
business3$license_year <- as.integer(business3$license_year)
class(business3$license_year) #integer

business3$time_since_last <- ifelse(business3$num_prev_lic == 1 & business3$status_mostrecent == 'active', "only one active",
                                    ifelse (business3$num_prev_lic == 1 & business3$status_mostrecent == 'inactive', 2018-business3$license_year, 
                                            ifelse(business3$num_prev_lic > 1 & business3$status_mostrecent == 'active', 'newlicense',
                                                   ifelse(business3$num_prev_lic > 1 & business3$status_mostrecent == 'inactive', 2018-business3$license_year, 00000))))

business3_noDup_mostrec <-
  business3 %>%
  group_by(opaaccount_Copy) %>%
  slice(which.max(as.Date(mostrecent2)))

#check for consistency
#temp<-subset(business3_noDup_mostrec,select=c('ActiveInactive','status_mostrecent'))
table(business3_noDup_mostrec$ActiveInactive)
table(business3_noDup_mostrec$status_mostrecent)

names(business3_noDup_mostrec)

business3_noDup_mostrec$activitystat<-ifelse(business3_noDup_mostrec$ActiveInactive=="Active",1,0)
business3_noDup_mostrec$activitystat<-as.factor(business3_noDup_mostrec$activitystat)
class(business3_noDup_mostrec$activitystat)

#clean columns again
license_final<-subset(business3_noDup_mostrec, select=-c(objectid,licensenum,revenuecod,
                                                         censustrac,legalfirst,legallastn,
                                                         business_n,eclipse_ad,licensesta,
                                                         fulladdres,zip,ActiveInactive))


#variable classes to change
class(license_final$num_prev_lic)
class(license_final$mostrecent2)
class(license_final$status_mostrecent)
class(license_final$license_year)
class(license_final$time_since_last)
class(license_final$numberofun)

license_final$status_mostrecent<-as.factor(license_final$status_mostrecent)
license_final$time_since_last<-as.factor(license_final$time_since_last)
license_final$numberofun<-as.factor(license_final$numberofun)

#fixed effect for unkown inactive date
license_final$ZZinactivedate<-ifelse(is.na(license_final$inactiveda), 1,0)
license_final$ZZinactivedate<-as.factor(license_final$ZZinactivedate)

license_final #cleaned business license data

#####   TAX DELINQUENCY   #####
tax_delinquency <- read.csv('tax_delinquency_clean.csv')
names(tax_delinquency)

#clean tax delinquency, remove columns that we don't need
tax_delinquency_clean<-subset(tax_delinquency,select=-c(index,objectid,zip_code,zip_4,co_owner,mailing_address,
                                                        mailing_city,mailing_state,mailing_zip))

names(tax_delinquency_clean)

#check class of columns
class(tax_delinquency_clean$owes)
class(tax_delinquency_clean$num_years_owed)
class(tax_delinquency_clean$total_due)
tax_delinquency_clean$owes<-as.factor(tax_delinquency_clean$owes)

#####   CENSUS DATA   #####

#read in census data
demos<-read.csv('censuscombined.csv')
names(demos)
unlist(lapply(demos,class))

#remove tracts without population
demos <- subset(demos, demos$pct_white!= 0 & demos$pct_black!=0)
#remove tracts where total population is 0
demos <- subset(demos, tot_population > 0)


#recalculate pct of each race accordingly
demos$pct_white <- demos$num_white/demos$tot_population
demos$pct_black <- demos$num_black/demos$tot_population
demos$pct_AmIndian <- demos$num_AmIndian/demos$tot_population
demos$pct_asian <- demos$num_asian/demos$tot_population
demos$pct_nativeHawaiian <- demos$num_nativeHawaiian/demos$tot_population
demos$pct_other <- demos$num_other/demos$tot_population
demos$pct_mixedrace <- demos$num_mixedrace/demos$tot_population
demos$pct_notwhte <- demos$num_notwhite/demos$tot_population
demos$pct_foreignborn <- demos$num_foreignborn/demos$tot_population
demos$pct_naturalizedcitizen <- demos$num_naturalizedcitizen/demos$tot_population
demos$pct_noncitizen <- demos$num_noncitizen/demos$tot_population
demos$pct_female_hh <- demos$tot_female_hh/demos$tot_population
demos$per_elevated <- demos$num_elevated/demos$tot_population
#numeric columns
demos$Pop_pov_Est <- as.numeric(demos$Pop_pov_Est) #integer to numeric
demos$Pop_below_pov <- as.numeric(demos$Pop_below_pov) #integer to numeric
demos$pct_pop_below <- as.numeric(demos$pct_pop_below) #factor to numeric
demos$pct_rent_Inc <- as.numeric(demos$pct_rent_of_Inc) #null to numeric
demos$rent_burden <- as.numeric(demos$rent_burden) #integer to numeric
demos$num_screened <- as.numeric(demos$num_screened) #integer to numeric
demos$num_elevated <- as.numeric(demos$num_elevated) #integer to numeric

unlist(lapply(demos,class))

#remove unencessary columns
demos <- subset(demos, select=-c(num_white,num_black,num_AmIndian,num_asian,
                                 num_nativeHawaiian,num_other,num_mixedrace,
                                 num_notwhite,num_foreignborn,num_naturalizedcitizen,
                                 num_noncitizen,tot_female_hh))

summary(demos)
unlist(lapply(demos,class))


#check where NAs were introduced by coercion
colSums(is.na(demos))

#fixed effect for lead NAs
demos$ZZleadelev<-ifelse(is.na(demos$num_elevated), 1,0)
demos$ZZleadelev<-as.factor(demos$ZZleadelev)

#rename census demos
#names(demos)[1]<-"FIPS"
demos$FIPS<-demos$tract #numeric


##################DOCKET############################

#join in docket data
phillyshape<-st_read('phillyproj2010.shp')
phillyshape = phillyshape%>%st_set_crs(4326)
st_transform(phillyshape, 4326)

evictionsshape<-st_read('evictions_clip.shp')
evictionsshape = evictionsshape%>%st_set_crs(4326)
st_transform(evictionsshape, 4326)
st_intersection(phillyshape, evictionsshape)

evictionsdf2<-read.csv('docket_clean.csv')
names(evictionsshape)
names(evictionsdf2)

#join evictions shape and docket by id

class(evictionsshape$id) #factor
evictionsshape$id_copy<-as.character(evictionsshape$id)
class(evictionsshape$id_copy) #character

class(evictionsdf2$id)   #factor
evictionsdf2$id_copy<-as.character(evictionsdf2$id)
class(evictionsdf2$id_copy) #character

evictions_join<-left_join(evictionsshape,evictionsdf2, by = "id_copy")
names(evictions_join)


#CLEAN DOCKET DATA
evictions_clean<-subset(evictions_join, select=c(ARC_Addres, year.x,month.x,plaintiff.y,X,Y, 
                                                 eviction_count.y,case_count))
evictions_clean<-as.data.frame((evictions_clean))

#rename counts
names(evictions_clean)[2]<-"evict_year"
names(evictions_clean)[3]<-"evict_month"
names(evictions_clean)[7]<-"home_evict_count"
names(evictions_clean)[8]<-"LL_eviction_count"

sum(is.na(evictions_clean$LL_eviction_count))
sum(is.na(evictions_clean$home_evict_count))

#filter to just evict_year
evictions_clean2<-(subset (evictions_clean, evictions_clean$evict_year== 2016 | 
                             evictions_clean$evict_year== 2017 | 
                             evictions_clean$evict_year== 2018))

#remove complete duplicates
evictions_clean3 <- evictions_clean2[!duplicated(evictions_clean2), ]

#group by address and find eviction sums for each address
docket_count<-evictions_clean3%>% 
  group_by(ARC_Addres)%>%
  summarize(sum_home_evict = sum(home_evict_count),
            sum_LL_evict = sum(LL_eviction_count))

#join back to evictions_clean3
evictions_clean4 <- left_join(evictions_clean3,docket_count, by='ARC_Addres')

class(evictions_clean4$evict_year) #factor
evictions_clean4$evict_year_copy <- as.numeric(as.character(evictions_clean4$evict_year))

#keep unique address copies
evictions_unique <- 
  evictions_clean4 %>%
  group_by(ARC_Addres) %>%
  slice(which.max(evict_year_copy))


####KNN violations and licenses###
#KNN FUNCTION
nn_function <- function(measureFrom,measureTo,k) {
  
  nn <-   
    get.knnx(measureTo, measureFrom, k)$nn.dist
  
  output <-
    as.data.frame(nn) %>%
    rownames_to_column(var = "thisPoint") %>%
    gather(points, point_distance, V1:ncol(.)) %>%
    arrange(as.numeric(thisPoint)) %>%
    group_by(thisPoint) %>%
    summarize(pointDistance = mean(point_distance)) %>%
    arrange(as.numeric(thisPoint)) %>% 
    dplyr::select(-thisPoint)
  
  return(output)  
}

#KNN
#create clean sf without nas in geometry
license_final2<- license_final[!is.na(license_final$lat), ]
license.sf<-st_as_sf(license_final2, coords=c('lng','lat'),crs=4326)

rentals_geo<-opa_final[!is.na(opa_final$lat), ]
names(rentals_geo)
sum(is.na(rentals_geo$lat))
rentals.sf<-st_as_sf(rentals_geo, coords=c('lng','lat'), crs=4326)

violations_geo<-violations_final[!is.na(violations_final$lat), ]
violations.sf<-st_as_sf(violations_geo, coords=c('lng','lat'), crs=4326)

#pull the coordinates
rentals.xy <- 
  rentals.sf %>%
  cbind(.,st_coordinates(rentals.sf))  %>%
  st_set_geometry(NULL) %>%
  dplyr::select(X,Y) %>%
  as.matrix()

license.xy <- 
  license.sf %>%
  cbind(.,st_coordinates(license.sf))  %>%
  st_set_geometry(NULL) %>%
  dplyr::select(X,Y) %>%
  as.matrix()

violations.xy <- 
  violations.sf %>%
  cbind(.,st_coordinates(violations.sf))  %>%
  st_set_geometry(NULL) %>%
  dplyr::select(X,Y) %>%
  as.matrix()

#add KNN as columns
knn_licenses_5<-nn_function(rentals.xy,license.xy,5)
rentals_geo <- cbind(rentals_geo, knn_licenses_5) 
names(rentals_geo)
colnames(rentals_geo)[colnames(rentals_geo)=="pointDistance"] <- "knn_lic_5"

knn_violations_5<-nn_function(rentals.xy,violations.xy,5)
rentals_geo<-cbind(rentals_geo, knn_violations_5)
colnames(rentals_geo)[colnames(rentals_geo)=="pointDistance"] <- "knn_viol_5"

names(rentals_geo)
opa_final2<-as.data.frame(rentals_geo)



opa_final2        #opa - housing characteristics
license_final     #business licenses

Following the process of data cleaning, we joined the datasets together and did some additional feature engineering.

### JOIN #1: OPA + BUSINESS ###

#check common column
names(opa_final2) #parcel_number_copy
names(license_final2) #opaaccount_copy

#check correct class before join
class(opa_final2$parcel_number_copy) #numeric
class(license_final2$opaaccount_Copy) #numeric

#make join
join1 <- left_join(opa_final2,license_final2,by=c('parcel_number_copy'='opaaccount_Copy'))
names(join1)

#see how many didn't join
sum(is.na(join1$num_prev_lic)) #188,834 didn't join

#check correct join
join1 %>%
  dplyr::select(location,fulladdres_Copy) %>%
  head()


######### Additional FE column ##########

#1-rental but no license; 0- rental with license
#rental but no license <- anything in opa that didnt join with business
join1$rental_butNolic <- ifelse(is.na(join1$licensetyp), 1,0)
join1$rental_butNolic <- as.factor(join1$rental_butNolic)


##############################################################################################

### JOIN #2: OPA/BUSINESS + VIOLATIONS ###

#check common column
sort(names(join1)) #parcel_number_copy
names(violations_final) #opa_account_num_copy

#check correct class before join
class(join1$parcel_number_copy) #numeric
class(violations_final$opa_account_num_copy) #numeric

#make join
join2 <- left_join(join1,violations_final, by=c('parcel_number_copy' = 'opa_account_num_copy'))
names(join2)

#see how many didn't join
sum(is.na(join2$viol_count)) #247,043 didn't join



##############################################################################################

### JOIN #3: OPA/BUSINESS/VIOLATIONS + TAX DELINQUENCY ###

#check common column
sort(names(join2)) #parcel_number_copy
names(tax_delinquency_clean) #opa_number

#check correct class before join
class(join2$parcel_number_copy) #numeric
class(tax_delinquency_clean$opa_number) #integer
tax_delinquency_clean$opa_number_copy <- as.numeric(tax_delinquency_clean$opa_number)
class(tax_delinquency_clean$opa_number_copy) #numeric

#join by parcel_number_copy and opa_number_copy

#make join
join3 <- left_join(join2, tax_delinquency_clean, by=c('parcel_number_copy'= 'opa_number_copy'))
names(join3)

#see how many didn't join
sum(is.na(join3$owes)) #280,005 didn't join


##############################################################################################

### JOIN #4: OPA/BUSINESS/VIOLATIONS/TAX DELINQUENCY + CENSUS ###

# SPATIAL JOIN: join3 + census tracts #

#open census tracts shapefile
census_tracts<-st_read('Census_Tracts_2010.shp') 

#test
ggplot() + geom_sf(data=census_tracts)

#make join3 data frame in sf object
join3_geo<-st_as_sf(join3, coords=c('lng.x','lat.x'), crs=4326)
join3_geo

#check
head(join3$lat.x)
head(join3$lng.x)
head(join3_geo$geometry)

#spatial join census_tracts to join3_geo
join4<-st_join(join3_geo, census_tracts, join=st_within)


##############################################################################################

### JOIN #5: OPA/BUSINESS/VIOLATIONS/TAX DELINQUENCY/CENSUS_TRACTS + DEMOGRAPHIC INFO ###

#check common column
join4$FIPS_geoidCopy <- as.numeric(as.character(join4$GEOID10))
sort(names(join4)) #FIPS_geoidCopy
names(demos) #FIPS

#check correct class before join
class(join4$FIPS_geoidCopy) #numeric
class(demos$FIPS) #numeric

#make join
join5<-left_join(join4, demos, by=c("FIPS_geoidCopy" = "FIPS"))
names(join5)

#see how many didn't join
sum(is.na(join5$pct_white)) # 2,347 didn't join



##############################################################################################

### JOIN #6: OPA/BUSINESS/VIOLATIONS/TAX DELINQUENCY/CENSUS_TRACTS/DEMOGRAPHIC INFO + DOCKET ###

#check common column
sort(names(join5)) #location
names(evictions_unique) #ARC_Addres

#check correct class before join
class(join5$location) #factor
join5$location_copy <- as.character(join5$location)
class(join5$location_copy) #character

class(evictions_unique$ARC_Addres) #factor
evictions_unique$ARC_Addres_copy <- as.character(evictions_unique$ARC_Addres)
class(evictions_unique$ARC_Addres_copy) #character

#make join
join6 <-left_join(join5, evictions_unique, by=c("location_copy"="ARC_Addres_copy"))
names(join6)

#see how many didn't join
sum(is.na(join6$evict_year)) # 294,060  didn't join


##############################################################################################

### JOIN #7: OPA/BUSINESS/VIOLATIONS/TAX DELINQUENCY/CENSUS_TRACTS/DEMOGRAPHIC INFO/DOCKET ###
###          + PHILLY NEIGHBORHOODS

#Add in Philly neighborhoods shapefile
phillyneighbs<-st_read('Neighbs_4326.shp')

#clean
phillyneighbs<-subset(phillyneighbs,select=c(MAPNAME, geometry))

#test
ggplot() + geom_sf(data=phillyneighbs)

#remove unused geometry columns
join6_copy<-subset(join6,select=-c(geometry.y,geometry))
#reassign active geometry column
st_geometry(join6_copy)<-'geometry.x'
#remove unused geometry columns
join6_copy<-subset(join6_copy,select=-c(geometry))

#check
head(join6_copy$geometry.x)

#spatial join census_tracts to join3_geo
join7<-st_join(join6_copy, phillyneighbs, join=st_within)

Next we cleaned the joined dataset to prepare it for the exploratory analysis and model building.

##### CLEANING #####
colSums(is.na(join7))

#removing unecessary rows first
data <- subset(join7, select=-c(geocode_y,geocode_x,street_add,opaaccount,business_m,lng.y,lat.y,fulladdres_Copy,
                                lng,opa_account_num.x,address,ownername,lat.x.x,opa_account_num.y,opa_number,
                                lat.y.y,lon,OBJECTID,STATEFP10,COUNTYFP10,TRACTCE10,NAME10,NAMELSAD10,MTFCC10,
                                FUNCSTAT10,ALAND10,AWATER10,LOGRECNO,FIPS_geoidCopy,tract,location_copy,ARC_Addres,
                                plaintiff.y,X,Y,street_address))


colSums(is.na(data))

#assume NA for illegalRental is 0; 
data$illegalRental_dv <- data$illegalRental
data$illegalRental_dv <- ifelse(is.na(data$illegalRental) | data$illegalRental=='0',0,1)

table(data$illegalRental_dv)
#0-302586;1-2225

colSums(is.na(data))
#make copy
data2<-data

#checked that all na's in this are 0 illegalRental;
data2<-data2[!is.na(data2$off_street_open),]
data2<-data2[!is.na(data2$GEOID10),]

#second filtering of columns
data2<-subset(data2,select=-c(illegalRental,SFIllegalRent,MFIllegalRent,RoomIllegalRent,
                              mostrecent,initialiss,numberofun,legalname,evict_year,licensetyp))


colSums(is.na(data2))

data2$tot_pop2<-data2$tot_population
data2$tot_pop2<-ifelse(is.na(data2$tot_pop2),0,data2$tot_pop2)

data2$pct_white2<-data2$pct_white
data2$pct_white2<-ifelse(is.na(data2$pct_white2),0,data2$pct_white2)

data2$pct_black2<-data2$pct_black
data2$pct_black2<-ifelse(is.na(data2$pct_black2),0,data2$pct_black2)

data2$pct_AmIndian2<-data2$pct_AmIndian
data2$pct_AmIndian2<-ifelse(is.na(data2$pct_AmIndian2),0,data2$pct_AmIndian2)

data2$pct_asian2<-data2$pct_asian
data2$pct_asian2<-ifelse(is.na(data2$pct_asian2),0,data2$pct_asian2)

data2$pct_nativeHawaiian2<-data2$pct_nativeHawaiian
data2$pct_nativeHawaiian2<-ifelse(is.na(data2$pct_nativeHawaiian2),0,data2$pct_nativeHawaiian2)

data2$pct_other2<-data2$pct_other
data2$pct_other2<-ifelse(is.na(data2$pct_other2),0,data2$pct_other2)

data2$pct_mixedrace2<-data2$pct_mixedrace
data2$pct_mixedrace2<-ifelse(is.na(data2$pct_mixedrace2),0,data2$pct_mixedrace2)

data2$pct_notwhte2<-data2$pct_notwhte
data2$pct_notwhte2<-ifelse(is.na(data2$pct_notwhte2),0,data2$pct_notwhte2)

data2$pct_foreignborn2<-data2$pct_foreignborn
data2$pct_foreignborn2<-ifelse(is.na(data2$pct_foreignborn2),0,data2$pct_foreignborn2)

data2$pct_naturalizedcitizen2<-data2$pct_naturalizedcitizen
data2$pct_naturalizedcitizen2<-ifelse(is.na(data2$pct_naturalizedcitizen2),0,data2$pct_naturalizedcitizen2)

data2$pct_noncitizen2<-data2$pct_noncitizen
data2$pct_noncitizen2<-ifelse(is.na(data2$pct_noncitizen2),0,data2$pct_noncitizen2)

data2$Pop_pov_Est2<-data2$Pop_pov_Est
data2$Pop_pov_Est2<-ifelse(is.na(data2$Pop_pov_Est2),0,data2$Pop_pov_Est2)

data2$Pop_below_pov2<-data2$Pop_below_pov
data2$Pop_below_pov2<-ifelse(is.na(data2$Pop_below_pov2),0,data2$Pop_below_pov2)

data2$pct_pop_below2<-data2$pct_pop_below
data2$pct_pop_below2<-ifelse(is.na(data2$pct_pop_below2),0,data2$pct_pop_below2)

data2$pct_rent_Inc2<-data2$pct_rent_Inc
data2$pct_rent_Inc2<-ifelse(is.na(data2$pct_rent_Inc2),0,data2$pct_rent_Inc2)

data2$rent_burden2<-data2$rent_burden
data2$rent_burden2<-ifelse(is.na(data2$rent_burden2),0,data2$rent_burden2)

data2$pct_female_hh2<-data2$pct_female_hh
data2$pct_female_hh2<-ifelse(is.na(data2$pct_female_hh2),0,data2$pct_female_hh2)

data2$per_elevated2<-data2$per_elevated
data2$per_elevated2<-ifelse(is.na(data2$per_elevated2),0,data2$per_elevated2)

data2$pct_rent_Inc2<-data2$pct_rent_Inc
data2$pct_rent_Inc2<-ifelse(is.na(data2$pct_rent_Inc2),0,data2$pct_rent_Inc2)

data2$ZZleadelev2<-data2$ZZleadelev
data2$ZZleadelev2<-ifelse(is.na(data2$ZZleadelev2) | data2$ZZleadelev2=='0','0','1')
data2$ZZleadelev2<-as.factor(data2$ZZleadelev2)

data2$home_evict_count2<-data2$home_evict_count
data2$home_evict_count2<-ifelse(is.na(data2$home_evict_count2),0,data2$home_evict_count2)

data2$LL_eviction_count2<-data2$LL_eviction_count
data2$LL_eviction_count2<-ifelse(is.na(data2$LL_eviction_count2),0,data2$LL_eviction_count2)
data2$LL_eviction_count2<-as.numeric(data2$LL_eviction_count2)

data2$sum_home_evict2_dummy<-ifelse(is.na(data2$sum_home_evict),1,0)
data2$sum_home_evict2_dummy<-as.factor(data2$sum_home_evict2_dummy)
data2$sum_home_evict2<-data2$sum_home_evict
data2$sum_home_evict2<-ifelse(is.na(data2$sum_home_evict2),0,data2$sum_home_evict2)

data2$sum_LL_evict2_dummy<-ifelse(is.na(data2$sum_LL_evict),1,0)
data2$sum_LL_evict2_dummy<-as.factor(data2$sum_LL_evict2_dummy)
data2$sum_LL_evict2<-data2$sum_LL_evict
data2$sum_LL_evict2<-ifelse(is.na(data2$sum_LL_evict2),0,data2$sum_LL_evict2)


#remove columns
data3<- subset(data2, select=-c(tot_population,pct_white,pct_black,pct_AmIndian,pct_asian,
                                pct_nativeHawaiian,pct_other,pct_other,pct_mixedrace, pct_notwhte,
                                pct_foreignborn,pct_naturalizedcitizen,pct_noncitizen, 
                                Pop_pov_Est,Pop_below_pov,pct_pop_below,pct_rent_Inc, 
                                rent_burden,pct_female_hh,num_elevated,num_screened,
                                per_elevated,pct_rent_Inc,pct_rent_of_Inc,home_evict_count,LL_eviction_count,
                                sum_home_evict,sum_LL_evict,ZZleadelev))


colSums(is.na(data3))

data3$owes2<-data3$owes
data3$owes2<-ifelse(is.na(data3$owes2),0,data3$owes2)
data3$owes2<-as.factor(data3$owes2)

data3$hist_extviol2<-data3$hist_extviol
data3$hist_extviol2<-ifelse(is.na(data3$hist_extviol2) | data3$hist_extviol2=='0','0','1')
data3$hist_extviol2<-as.factor(data3$hist_extviol2)

data3$hist_intviol2<-data3$hist_intviol
data3$hist_intviol2<-ifelse(is.na(data3$hist_intviol2) | data3$hist_intviol2=='0','0','1')
data3$hist_intviol2<-as.factor(data3$hist_intviol2)

data3$hist_illegalRent2<-data3$hist_illegalRent
data3$hist_illegalRent2<-ifelse(is.na(data3$hist_illegalRent2) | data3$hist_illegalRent2=='0','0','1')
data3$hist_illegalRent2<-as.factor(data3$hist_illegalRent2)

data3$had_anyviolation2<-data3$had_anyviolation
data3$had_anyviolation_dummy<-ifelse(is.na(data3$had_anyviolation2),1,0)
data3$had_anyviolation_dummy<-as.factor(data3$had_anyviolation_dummy)
data3$had_anyviolation2<-ifelse(is.na(data3$had_anyviolation2),'0',
                                ifelse(data3$had_anyviolation2=='1','1','0'))
data3$had_anyviolation2<-as.factor(data3$had_anyviolation2)

data3$other_violations2<-data3$other_violations
data3$other_violations2<-ifelse(is.na(data3$other_violations2),0,data3$other_violations2+1)

data3$issue_fixed_dummy<-ifelse(is.na(data3$issue_fixed),1,0)
data3$issue_fixed_dummy<-as.factor(data3$issue_fixed_dummy)

data3$issue_fixed_y<-ifelse(is.na(data3$issue_fixed),'0',
                            ifelse(data3$issue_fixed=='1','1','0'))
data3$issue_fixed_y<-as.factor(data3$issue_fixed_y)

data3$issue_fixed_n<-ifelse(is.na(data3$issue_fixed),'0',
                            ifelse(data3$issue_fixed=='0','1','0'))
data3$issue_fixed_n<-as.factor(data3$issue_fixed_n)


data3$active_case_dummy<-ifelse(is.na(data3$active_case),1,0)
data3$active_case_dummy<-as.factor(data3$active_case_dummy)

data3$active_case_y<-ifelse(is.na(data3$active_case),'0',
                            ifelse(data3$active_case=='1','1','0'))
data3$active_case_y<-as.factor(data3$active_case_y)

data3$active_case_n<-ifelse(is.na(data3$active_case),'0',
                            ifelse(data3$active_case=='0','1','0'))
data3$active_case_n<-as.factor(data3$active_case_n)


data3$total_due_dummy<-ifelse(is.na(data3$total_due),1,0)
data3$total_due_dummy<-as.factor(data3$total_due_dummy)
data3$total_due2<-data3$total_due 
data3$total_due2<-ifelse(is.na(data3$total_due2),0,data3$total_due2)

data3$num_years_owed_dummy<-ifelse(is.na(data3$num_years_owed),1,0)
data3$num_years_owed_dummy<-as.factor(data3$num_years_owed_dummy)

#0 because there is nothing due
data3$num_years_owed2<-data3$num_years_owed 
data3$num_years_owed2<-ifelse(is.na(data3$num_years_owed2),0,data3$num_years_owed2)

data3$lentim_sincelast_viol_dummy<-ifelse(is.na(data3$lentim_sincelast_viol),1,0)
data3$lentim_sincelast_viol_dummy<-as.factor(data3$lentim_sincelast_viol_dummy)

#0 because never had violation
data3$lentim_sincelast_viol2<-data3$lentim_sincelast_viol 
data3$lentim_sincelast_viol2<-ifelse(is.na(data3$lentim_sincelast_viol2),0,data3$lentim_sincelast_viol2)

data3$viol_count_dummy<-ifelse(is.na(data3$viol_count),1,0)
data3$viol_count_dummy<-as.factor(data3$viol_count_dummy)

#0 because never had violation
data3$viol_count2<-data3$viol_count
data3$viol_count2<-ifelse(is.na(data3$viol_count2),0,data3$viol_count2)


data4 <- subset(data3, select=-c(owes,had_extvio,had_intvio,had_anyviolation,had_illegalRent,
                                 hist_illegalRent,hist_extviol,hist_intviol,
                                 exterior,interior,issue_fixed,active_case,owner,
                                 total_due,num_years_owed,lentim_sincelast_viol,
                                 viol_count,violationtype,other_violations))

colSums(is.na(data4))


#Length of time since last license if most recent is inactive
#NA's we assume the license is active
data4$time_since_last_dummy<-ifelse(is.na(data4$time_since_last),1,0)
data4$time_since_last_dummy<-as.factor(data4$time_since_last_dummy)
data4$time_since_last2<-data4$time_since_last
data4$time_since_last2<-ifelse(is.na(data4$time_since_last2),'only one active',data4$time_since_last2)
data4$time_since_last2<-as.factor(data4$time_since_last2)

data4$owner_count_dummy<-ifelse(is.na(data4$owner_count_dummy),1,0)
data4$owner_count_dummy<-as.factor(data4$owner_count_dummy)

#assumes NA is at least 1 owner
data4$owner_count<-ifelse(is.na(data4$owner_count),1,data4$owner_count)

#assumes NA is not a multiple unit
data4$multipleunits_dummy<-ifelse(is.na(data4$multipleunits),1,0)
data4$multipleunits_dummy<-as.factor(data4$multipleunits_dummy)

data4$multipleunits_y<-ifelse(is.na(data4$multipleunits),'0',
                              ifelse(data4$multipleunits=='1','1','0'))
data4$multipleunits_y<-as.factor(data4$multipleunits_y)


data4$ZZinactivedate<-ifelse(is.na(data4$ZZinactivedate),1,0)
data4$ZZinactivedate<-as.factor(data4$ZZinactivedate)

#number of previous business license
#dummy for unknown number of previous licenses
data4$num_prev_lic_dummy<-ifelse(is.na(data4$num_prev_lic),1,0)
data4$num_prev_lic_dummy<-as.factor(data4$num_prev_lic_dummy)
#assumes at least 1 if NA
data4$num_prev_lic<-ifelse(is.na(data4$num_prev_lic),1,data4$num_prev_lic)


data4$status_mostrecent_dummy<-ifelse(is.na(data4$status_mostrecent),1,0)
data4$status_mostrecent_dummy<-as.factor(data4$status_mostrecent_dummy)

data4$status_mostrecent_active<-ifelse(is.na(data4$status_mostrecent),'0',
                                       ifelse(data4$status_mostrecent=='active','1','0'))
data4$status_mostrecent_active<-as.factor(data4$status_mostrecent_active)
data4$status_mostrecent_inactive<-ifelse(is.na(data4$status_mostrecent),'0',
                                         ifelse(data4$status_mostrecent=='inactive','1','0'))
data4$status_mostrecent_inactive<-as.factor(data4$status_mostrecent_inactive)

#individual or company
data4$legalentit_dummy<-ifelse(is.na(data4$legalentit),1,0)
data4$legalentit_dummy<-as.factor(data4$legalentit_dummy)

data4$legalentit_indiv<-ifelse(is.na(data4$legalentit),'0',
                               ifelse(data4$legalentit=='Individual','1','0'))
data4$legalentit_indiv<-as.factor(data4$legalentit_indiv)

data4$legalentit_compny<-ifelse(is.na(data4$legalentit),'0',
                                ifelse(data4$legalentit=='Company','1','0'))
data4$legalentit_compny<-as.factor(data4$legalentit_compny)


data5<-subset(data4,select=-c(violationdescription,caseresolutiondate,violationdate, 
                              caseaddeddate,organization,time_since_last,mostrecent2,
                              inactiveda,license_year,councildis,mostrecentinsp,
                              owner_count_dummy,year,mostrecentinsp2,status2,legalentit,
                              status_mostrecent,multipleunits))

colSums(is.na(data5))

data6<-subset(data5,select=-c(evict_month,evict_year_copy,activitystat))

colSums(is.na(data6))

Next we begain building the models using a stepwise regression that helped with model selection.

#the below shows the model building process for the violation story only. The process is the same for all 4 stories. 

###stepwise regression
library(MASS)
library(leaps)

violationsMod<-glm(illegalRental_dv~.,data=trainingViolations)


violationsStepMod<-stepAIC(violationsMod, direction = 'both', trace=FALSE)
summary(violationsStepMod)


##tuning the stepwise
violationModels <- regsubsets(illegalRental_dv~., data = trainingViolations, nvmax = 5, method = "seqrep") 
summary(violationModels)


# Set seed for reproducibility 
set.seed(123) 

# Set up repeated k-fold cross-validation 
train.control <- trainControl(method = "cv", number = 10) 

# Train the model 
Violation.step.model <- train(illegalRental_dv ~., data = trainingViolations, method = "leapSeq", tuneGrid = data.frame(nvmax = 1:5), trControl = train.control ) 
Violation.step.model$results


#Find best nvmax
Violation.step.model$bestTune

#Summary of the best model
summary(Violations.step.model$finalModel)

#find regression coefficients, meaning run a model with only the recommended variables
finalViolationModel<-glm(illegalRental_dv~issue_fixed_y+issue_fixed_dummy+viol_count2+hist_intviol2+hist_extviol2+owner_count+lentim_sincelast_viol2+active_case_y, data=trainingViolations, family=binomial)


summary(finalViolationModel)

####train and validate stories individually

trainingViolations_sample <-
  mutate(trainingViolations_sample, illegalRental_dv2 = ifelse(illegalRental_dv == 1, "illegal", "not_illegal"))


##caret train

##sensitivity = true positive
##specificity = true negative

ctrl <- trainControl(method = "cv", number = 50, classProbs=TRUE, summaryFunction=twoClassSummary) 

cvFitViolations<-train(illegalRental_dv2 ~ ., data = trainingViolations_sample %>%
                         dplyr::select(illegalRental_dv2,issue_fixed_y,issue_fixed_dummy,viol_count2,hist_intviol2,hist_extviol2,owner_count,
                                       lentim_sincelast_viol2,active_case_y), 
                       method='glm', family='binomial', metric='ROC', na.action=na.omit, trControl=ctrl)

cvFitViolations


#cross validation outputs
cvFitViolations$resample %>%
  dplyr::select(-Resample) %>%
  gather() %>%
  ggplot(., aes(value)) + 
  geom_histogram(bins=35) +
  facet_wrap(~key) +
  scale_x_continuous(limits = c(0, 1)) +
  labs(x="Goodness of Fit",
       y="Count",
       title="CV Goodness of Fit Metrics - Violations Story")

Next, because the models were unable to predict for illegal rentals because they are a rare event, we took a random sample of the data and used that smaller dataset to test the model.

####make subsets for each of the storie

####create a smaller sample where non-illegal rentals are 4X illegal rentals
##create a subset of the data 
data7_IR<-subset(data7, data7$illegalRental_dv == 1)
table(data7_IR$illegalRental_dv)

data7_noIR<-subset(data7, data7$illegalRental_dv == 0)
table(data7_noIR$illegalRental_dv)

#select 8900 random rows from temp 0 (2225*4=8900)
random_0 <- data7_noIR[sample(nrow(data7_noIR), 8900), ]

#rbind 8900 rows of 0 to the 2225 rows of 1
newsample <- rbind(data7_IR,random_0)

#shuffle rows - generate a random ordering
set.seed(1) ## make reproducible here, but not if generating many random samples
rand <- sample(nrow(newsample))
newsample<-newsample[rand,]

houseCharacterStory_sample<-subset(newsample, select=c(illegalRental_dv,quality_grade_copy,central_air_dummy,central_air_y,central_air_n,
                                            unfinished_unit,general_construction_dummy,general_construction_brickY,
                                            general_construction_brickN,house_extension_n,type_heater_yn,type_heater_dummy,
                                            other_bldg_yn,exempt_building_dummy,exempt_building_y,exempt_land_dummy,
                                            exempt_land_y,taxable_building_dummy,taxable_building_y,taxable_land_dummy,
                                            taxable_land_y,bedrooms0,bedrooms1,bedrooms2,bedrooms3,bedrooms4,bedrooms5plus,
                                            rooms0to2,rooms3to5,rooms6to8,rooms9to11,bathrooms1,bathrooms2to3,bathrooms4plus,story1,
                                            story2,story3,story4plus,year_built_copy,extcond_dummy,extcond_avg,intcond_dummy,
                                            intcond_avg,depth_copy,depth_dummy,frontage_copy,frontage_dummy,total_area_copy,
                                            total_area_dummy,market_value_copy,market_value_dummy,multipleunits_dummy,     multipleunits_n, multipleunits_y))



violationsStory_sample<-subset(newsample, select=c(illegalRental_dv, issue_fixed_y, issue_fixed_dummy, viol_count2, viol_count_dummy,
                                        hist_intviol2, hist_extviol2, owner_count, lentim_sincelast_viol2, active_case_y))


legalStory_sample<-subset(newsample, select=c(illegalRental_dv,absenteelandlord, outofPhillyLL, outofstateLL, owes2, status_mostrecent_inactive, status_mostrecent_active,  
                                   num_prev_lic,num_prev_lic_dummy, time_since_last2, time_since_last_dummy,rental_butNolic,
                                   ZZinactivedate))



demosGeoStory_sample<-subset(newsample, select=c(illegalRental_dv,pct_white2,pct_black2,pct_AmIndian2,pct_asian2,pct_nativeHawaiian2,
                                      pct_other2,pct_mixedrace2,pct_notwhte2, pct_foreignborn2,
                                      pct_naturalizedcitizen2,pct_noncitizen2,
                                      pct_pop_below2,rent_burden2,pct_female_hh2,per_elevated2,
                                      knn_lic_5, knn_viol_5, MAPNAME))

##the below shows the prediction steps for the violation story only, however the process is the same for all four stories. 

###create training and test set for new sample
inTrain <- createDataPartition(
  y = violationsStory_sample$illegalRental_dv, 
  p = .75, list = FALSE)
trainingViolations_sample <- violationsStory_sample[ inTrain,] #the new training set
testViolations_sample <- violationsStory_sample[-inTrain,]  #the new test set

##final modes with sample data
finalViolationModel<-glm(illegalRental_dv~issue_fixed_y+issue_fixed_dummy+viol_count2+hist_intviol2+hist_extviol2+owner_count+lentim_sincelast_viol2+active_case_y, data=trainingViolations_sample, family=binomial)
finalHouseCharacterModel<-glm(illegalRental_dv~extcond_avg+year_built_copy+story3+multipleunits_dummy+multipleunits_n, data=trainingHouseCharacter_sample, family=binomial)
finalLegalModel<-glm(illegalRental_dv~outofstateLL+owes2+num_prev_lic_dummy+num_prev_lic+time_since_last2+status_mostrecent_active,data=trainingLegal_sample, family=binomial)
finalDemosGeo<-glm(illegalRental_dv~knn_viol_5+knn_lic_5+per_elevated2+pct_black2, data=trainingDemosGeo_sample, family=binomial)

#### predict on sample data
#sample data
sample_test_violation <- predict(finalViolationModel, testViolations_sample, type='response')
testViolations_sample$classProbs <- sample_test_violation

sampletestProbs <- data.frame(Class=testViolations_sample$illegalRental_dv,
                              Probs=sample_test_violation)
head(sampletestProbs)

#distribution plot 
ggplot(sampletestProbs, aes(x = Probs, fill=as.factor(Class))) + 
  geom_density() +
  facet_grid(Class ~ .) +
  scale_fill_manual(values = c("dodgerblue4", "darkgreen"),
                    name = "",
                    labels = c("Not Illegal Rental","Illegal Rental")) +
  labs(x="Illegal Rental", y="Density of probabilities") +
  theme(strip.text.x = element_text(size = 18))


# distribution of the prediction score grouped by known outcome
ggplot( sampletestProbs, aes( Probs, color = as.factor(Class) ) ) + 
  geom_density( size = 1 ) +
  ggtitle( "Test Set's Predicted Score (Sample Data)" ) + 
  scale_fill_manual(values = c('dodgerblue4', 'darkgreen'),
                    name = "data", 
                    labels = c( "Not Illegal", "Illegal"))

The below code runs the ensemble model using a Random Forest regression.

#####ENSEMBLE 1 - Random Forest MODEL######

#violationsStory
#legalStory
#demosGeoStory
#houseCharacterStory

sum(is.na(data7))
data7$time_since_last2<-as.numeric(data7$time_since_last2)
newsample$year_built_copy<-as.factor(newsample$year_built_copy)

story.Violations <- newsample %>% dplyr::select(illegalRental_dv,issue_fixed_y,issue_fixed_dummy,viol_count2,hist_intviol2,hist_extviol2,owner_count,lentim_sincelast_viol2,active_case_y)
story.Legal <- newsample %>% dplyr::select(illegalRental_dv,outofstateLL,owes2,num_prev_lic_dummy,num_prev_lic,time_since_last2,status_mostrecent_active)
story.housingCharacter <- newsample %>% dplyr::select(illegalRental_dv,extcond_avg,year_built_copy,story3,multipleunits_dummy,multipleunits_n)
story.DemosGeo <-newsample %>% dplyr::select(illegalRental_dv,knn_viol_5,knn_lic_5,per_elevated2,pct_black2)

ensembleFunction <- function(inputStory, iterations, sampleSize) {
  #create an empty dataframe
  endList <- data.frame(matrix(NA, nrow = nrow(inputStory), ncol = 0))
  #build n models with a m% test set
  for (i in 1: iterations ) {
    sample <- sample.int(n = nrow(inputStory), size = floor( sampleSize *nrow(inputStory)), replace = F)
    train <- inputStory[sample, ]
    #train the model
    thisTrain <- randomForest(train$illegalRental_dv ~ ., data = train)
    #create a vector of predictions for this model.
    thisPrediction = exp(predict(thisTrain, newsample, type='response'))
    #append the predictions to the data frame         
    endList <- cbind(endList, thisPrediction)
  }
  #label each prediction column 1 through n
  colnames(endList) <- paste("prediction", 1:ncol(endList), sep = "")
  #create a data frame of average, median and sd predictions. I have turned off the latter two for now
  return(data.frame(mean.Prediction = rowMeans(endList), na.rm=T))
}

allStories.predictions_RF <- 
  data.frame(
    cbind(
      story1.mean.predictions = ensembleFunction(story.Violations,100,.4)[,1],
      story2.mean.predictions = ensembleFunction(story.Legal,100,.4)[,1],
      story3.mean.predictions = ensembleFunction(story.housingCharacter,100,.4)[,1],
      story4.mean.predictions = ensembleFunction(story.DemosGeo,100,.4)[,1]
    ))

head(allStories.predictions_RF)

#cbind all story predictions back to newsample dataset
names(newsample)
IR<-subset(newsample, select=c(illegalRental_dv))
allStories.predictions2_rf<-cbind(allStories.predictions_RF, IR)
class(allStories.predictions_RF)

##run random forest regression on predictions
finalReg_rf <- randomForest(newsample$illegalRental_dv ~ .,data=allStories.predictions2_rf )
summary(finalReg_rf)


#final predictions
finalRegTest_rf <- predict(finalReg_rf, allStories.predictions2_rf, type='response')
allStories.predictions2_rf$classProbs <- finalRegTest_rf

finalRegProbs_rf <- data.frame(Class=allStories.predictions2_rf$illegalRental_dv,
                            Probs=finalRegTest_rf)
head(finalRegProbs_rf)

#confusion matrix with selected threshold
finalRegProbs_rf$predClass  = ifelse(finalRegProbs_rf$Probs > .65 ,1,0)
head(finalRegProbs_rf)
caret::confusionMatrix(factor(finalRegProbs_rf$Class), 
                       factor(finalRegProbs_rf$predClass), positive='1')


####ROC/AUC curves for sample RF ensemble
install.packages('plotROC')
library(plotROC)
names(finalRegProbs_rf2)

palette<-c("#c51b8a")
ggplot(finalRegProbs_rf2, aes(d = Class, m = Probs, colour = "#c51b8a")) + 
  geom_roc(n.cuts = 50, labels = FALSE) + 
  style_roc(theme = theme_grey) +
  geom_abline(slope = 1, intercept = 0, size = 1.5, color = '#feebe2') +
  labs(title = 'Final Model ROC', 
       subtitle = 'AUC = 0.97', 
       xlab = "False positive fraction", 
       ylab = "True positive fraction") +
  scale_colour_manual(values=palette)+
  theme_minimal()+
  theme(plot.title=element_text(hjust=0.5), legend.position='none')


#save AUC
sample_auc <- auc(finalRegProbs_rf2$Class, finalRegProbs_rf2$Probs) ##AUC =0.97


###How well does our model predict across neighborhoods

neighbs<-subset(newsample, select = c(MAPNAME))
head(neighbs)
finalRegProbs_rf_neighbs<-cbind(finalRegProbs_rf, neighbs)
head(finalRegProbs_rf_neighbs)

finalRegProbs_rf_neighbs2<-
  data.frame(class = finalRegProbs_rf_neighbs$Class, 
             nhood = finalRegProbs_rf_neighbs$MAPNAME, 
             probs = finalRegProbs_rf_neighbs$Probs, 
             predClass = finalRegProbs_rf_neighbs$predClass)


head(finalRegProbs_rf_neighbs2)  



options(yardstick.event_first = FALSE)


finalRegProbs_rf_neighbs2 %>%
  dplyr::select(-probs) %>%
  gather(Variable, predClass, -nhood, -class) #%>%
  group_by(nhood) %>%
  summarize(Sensitivity = round(yardstick::sens_vec(factor(class),factor(predClass)),2),
            Specificity = round(yardstick::spec_vec(factor(class),factor(predClass)),2),
            Accuracy = round(yardstick::accuracy_vec(factor(class),factor(predClass)),2),
            AUC = round(Metrics::auc(factor(class),factor(predClass)),2))
  
trueFalsePos_withNhood<-
  finalRegProbs_rf.withNhood %>%
    dplyr::select(-Probs) %>%
    mutate(truePositive = ifelse(predClass == 1 & Class ==1, 1,0),
           trueNegative = ifelse(predClass == 0 & Class ==0, 1,0)) %>%
    group_by(nhood) %>%
    summarize(sumTruePositive = sum(truePositive),
              sumTrueNegative = sum(trueNegative),
              totalObservations = n()) %>%
    mutate(truePositiveRate = sumTruePositive / totalObservations,
           trueNegativeRate = sumTrueNegative / totalObservations) %>%
    left_join(phillyneighbs, by=c("nhood"="MAPNAME")) %>%
    st_sf() 

head(trueFalsePos_withNhood)
trueFalsePos_withNhood<-subset(trueFalsePos_withNhood, trueFalsePos_withNhood$sumTruePositive>=1)
trueFalsePos_withNhood<-subset(trueFalsePos_withNhood, trueFalsePos_withNhood$totalObservations>=6)
head(trueFalsePos_withNhood)

7.2 Visuals

The visuals presented here largely, but not exclusively, feed into the exploratory analysis. To begin, we set up our visuals workspace.

#maptheme
mapTheme <- function(base_size = 12) {
  theme(
    text = element_text( color = "black"),
    plot.title = element_text(size = 14,colour = "black"),
    plot.subtitle=element_text(face="italic"),
    plot.caption=element_text(hjust=0),
    axis.ticks = element_blank(),
    panel.background = element_blank(),axis.title = element_blank(),
    axis.text = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=2)
  )
}

#Here is a function that will create quintile breaks. Note that the inclusion of `round` might not be advantageous if you are dealing with say, percentages.
quintileBreaks <- function(df,variable) {
  as.character(round(quantile(df[[variable]],
                              c(.01,.2,.4,.6,.8),na.rm=T)))
}
#Next, we can create a custom palette
palette5 <- c("#ffffcc","#a1dab4","#41b6c4","#2c7fb8","#253494")
palette5.2 <- c("#feebe2","#fbb4b9","#f768a1","#c51b8a","#7a0177")


#read in shapes
phillyshape2<-st_read('C:/Users/Fay/Desktop/GIS/Philly Shape Files/phillyproj2010.shp')%>%st_transform(crs=4326)

Next we make maps and graphics.

#Read in point data
#licenses
rental_license_sf<- st_read('li_business_licenses.shp')%>%st_transform(crs=4326)
rental_license_sf<-rental_license_sf%>%
  filter(rental_license_sf$licensesta=='Active','Active','Inactive,Expired,Closed,Revoked')%>%
  filter(rental_license_sf, business$licensetyp == 'Rental' | rental_license_sf$licensetyp == 'High Rise')
rental_license_sf$lng <- st_coordinates(rental_license_sf)[,1]
rental_license_sf$lat <- st_coordinates(rental_license_sf)[,2]

#EVICTIONS
evictionsshape<-st_read('evictions_4326.shp')%>%st_transform(crs=4326)%>%st_intersect(phillyshape)
evictionsshape$lng <- st_coordinates(evictionsshape)[,1]
evictionsshape$lat <- st_coordinates(evictionsshape)[,2]

#Illegal Rentals
violations2<-violations

#create sf
illegalRentalsf<-st_as_sf(violations_clean, coords=c('lng','lat'),crs=4326)
illegalRentalsf$lng <- st_coordinates(illegalRentalsf)[,1]
illegalRentalsf$lat <- st_coordinates(illegalRentalsf)[,2]

#IMPORT TAX DELINQUENCY
tax_delinquency_4326<-tax_delinquency_4326%>%
  st_transform(crs=4326)%>%
  st_intersection(phillyshape2)
names(tax_delinquency_4326)
tax_delinquency_4326$lng <- st_coordinates(tax_delinquency_4326)[,1]
tax_delinquency_4326$lat <- st_coordinates(tax_delinquency_4326)[,2]
class(tax_delinquency_4326$lat)
tax_delinquency_4326 <- tax_delinquency_4326[!is.na(tax_delinquency_4326$lat),]

#bar plot of kind of IR
names(illegal_rentals)
palette3<-c("#f768a1","#c51b8a","#7a0177")
ggplot(violations_clean, aes(x=factor(IRtype)))+
  geom_bar(stat="count", width=0.9, fill=palette3)+
  theme_minimal()+
  labs(x="Kind of Illegal Rental",
       y="Count",
       title="Illegal Rentals by Type")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(hjust = 0.5))

#*******************************************************************************************************************

#Density maps for illegal rentals and evictions

city <- st_read("http://data.phl.opendata.arcgis.com/datasets/405ec3da942d4e20869d4e1449a2be48_0.geojson")

#creating a 'windown' aka 'owin' in spatstat  

#create a bounding box for city
cityBox <- st_convex_hull(city)

#Convert `city` into a window we need for spatstat. Note it's class.
city.window <-
  as.owin(c(
    #xrange
    min(st_coordinates(cityBox)[,1]),
    max(st_coordinates(cityBox)[,1]),
    #yrange
    min(st_coordinates(cityBox)[,2]),
    max(st_coordinates(cityBox)[,2]))) 


#read in data
evictionpts <- st_read("evictions_4326.shp")

#note the latter fits in the former
plot(city.window)
plot(city,add=T)

#take a sample of your points, remove NAs, us st_difference to get rid of duplicate geometries
evictionpts2 <- evictionpts [1:1000,] %>% 
  mutate(x = st_coordinates(.)[,1],
         y = st_coordinates(.)[,2]) %>%
  filter(!is.na(x)) %>%
  st_difference()

#create the ppp and note the window
evictionpts2.ppp <- as.ppp(cbind(evictionpts2$x,evictionpts2$y), city.window)

#create the density 
evictionDensity <- raster(density(evictionpts2.ppp,.005, eps = .001))
plot(evictionDensity)

#to map - 
#Start by clipping the raster. The mask has to be a 'SpatialPolygons' 

evictionDensity.clip <- mask(evictionDensity, as(as(city,'Spatial'),'SpatialPolygons'))

#plot
palette5 <- c("#ffffcc","#a1dab4","#41b6c4","#2c7fb8","#253494")
palette5.2 <- c("#feebe2","#fbb4b9","#f768a1","#c51b8a","#7a0177")


gplot(evictionDensity.clip) +
  geom_raster(aes(fill = value))+
  scale_fill_distiller(palette="RdPu")+
  labs(title="Density of Evictions")+
  mapTheme()+
  theme(legend.position = "none") 

###repeat the process for illegal rentals

#**************************************************************************************************************#
#CONFUSION MATRICES
#read in data
confusionmatrix<-read.csv('rf0.65.csv')
confusionmatrix<-read.csv('rf0.65_2.csv')

confusionmatrix$Rate<-as.numeric(confusionmatrix$Rate)
names(confusionmatrix)[1]<-"Threshold"

##plot the confusion matrix for 0.4

confusionmatrix$variable <- factor(confusionmatrix$variable, levels = confusionmatrix$variable[order(-confusionmatrix$Rate)])

ggplot(confusionmatrix, aes(x=variable, y=Rate))+
  geom_bar(stat="identity", width=0.9, fill=palette4)+
  theme_minimal()+
  labs(x="",
       y="Rate (%)",
       title="Threshold 0.65",
       subtitle="Sensitivity: 0.8, Specificity: 0.93, Accuracy: 0.91")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(hjust = 0.5))


#**************************************************************************************************************#
#BAR CHART FOR INPUT ROWS

#read in csv
rowcount<-read.csv('properties_bar_chart_dat2.csv')

palette4<-c("#fbb4b9","#f768a1","#c51b8a","#7a0177")

ggplot(tips2, aes(x = reorder(day, -perc), y = perc)) + geom_bar(stat = "identity")
#create barplot
#bar plot of kind of IR
ggplot(rowcount, aes(x=reorder(Source,-Count), y=Count))+
  geom_bar(stat="identity", width=0.9, fill=palette4)+
  theme_minimal()+
  labs(x="",
       y="Count",
       title="Rows by Property Type")+
  theme(plot.title = element_text(hjust = 0.5))


#*************************************************************************************************************

#annual rent
annual_rent<-read.csv('annual_rent.csv')
names(annual_rent)[1]<-"Year"

#PLOT
ggplot(annual_rent, aes(x=Year, y=Median.Rent, group=1)) +
  geom_point(color="#c51b8a")+geom_line(color="#c51b8a")+
  labs(title="Median Rent in Philadelphia",
       x="Year",
       y="Rent (USD)")+
  theme(plot.title = element_text(hjust = 0.5),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA))+
  theme_minimal()


#PLOT ACCURACY BY NEIGHBORHOOD
#cbind sensitivity and specificty to phillyneighbs3

#sensitivity
ggplot() +
  geom_sf(data=phillyneighbs3, aes(fill='gray', colour='gray'))+
  geom_sf(data=phillyneighbs3, aes(fill=factor(ntile(sensitivity,5))))+
  scale_fill_manual(values=palette5.2,
                    labels = quintileBreaks(phillyneighbs3,"sensitivity")) +
  labs(title="Sensitivity by Neighborhood",
       fill="Sensitivity \nby Quantile")+
  mapTheme()

#Specificity
ggplot() +
  geom_sf(data=phillyneighbs3, aes(fill=factor(ntile(specificity,5))))+
  scale_fill_manual(values=palette5.2,
                    labels = quintileBreaks(phillyneighbs3,"specificity")) +
  labs(title="Specificity by Neighborhood",
       fill="Specificity by Quantile")+
  mapTheme()

A final step was developing correlation plots for key variables.

##correlation plots for exploratory analysis - for each chart, show % illegal rental for 0 and 1 for each variable (chart should have variable on x axis and %IR on the y) - how are illegal rentals distributed across variables?

IR<-subset(data7, data7$illegalRental_dv == '1')
NIR<-subset(data7, data7$illegalRental_dv == '0')

#history of interior violations
table(IR$hist_intviol2) #208 - 9%
table(NIR$hist_intviol2) #17618 - 6%
208/2225*100
17618/302539*100

#history of exterior violations
table(IR$hist_extviol2) #494 - 22%
table(NIR$hist_extviol2) #35835 - 12%
494/2225*100
35835/302539*100

#tax delinquency
table(IR$owes2) #383 - 17%
table(NIR$owes2) #24423 - 8%
383/2225*100
24423/302539*100

#out of state LL
table(IR$outofstateLL) #689 - 31%
table(NIR$outofstateLL) #152278 - 50%
689/2225*100
152278/302539*100

table(IR$outofPhillyLL) #1244 - 36%
table(NIR$outofPhillyLL) #198511 - 66%
1244/2225*100
198511/302539*100

#multiple units
table(IR$multipleunits_y) #33 - 1%
table(NIR$multipleunits_y) #271 - .09%
33/2225*100
271/302539*100

#built pre 1900
table(IR$built_pre1900) #138 - 6%
table(NIR$built_pre1900) #50126 - 16%
138/2225*100
50126/302539*100

correlations<-read.csv('correlations.csv')
correlations$PctIR<-as.numeric(correlations$PctIR)
correlations$PctIR.1<-as.numeric(correlations$PctIR.1)
correlations$PctIR.2<-as.numeric(correlations$PctIR.2)
correlations$PctIR.3<-as.numeric(correlations$PctIR.3)
correlations$PctIR.4<-as.numeric(correlations$PctIR.4)
correlations$PctIR.5<-as.numeric(correlations$PctIR.5)


###plots
palette5.2 <- c("#feebe2","#fbb4b9","#f768a1","#c51b8a","#7a0177")

histIntViol<-ggplot(correlations, aes(x =IRStat , y = PCTIntViol)) + 
  geom_bar(stat='identity', wideth =0.9, fill='#fbb4b9') +
  theme_minimal() +
  ylim(0,60) +
  labs(y = "% interior violations", 
       title = 'Properties with interior violations') +
  theme(axis.title.x=element_blank()) +
  theme(plot.title = element_text(hjust = 0.5))
histIntViol

histExtViol<-ggplot(correlations, aes(x =IRStat , y = PCTExtViol)) + 
  geom_bar(stat='identity', wideth =0.9, fill='#fbb4b9') +
  theme_minimal() +
  ylim(0,60) +
  labs(y = "% exterior violations", 
       title = 'Properties with exterior violations') +
  theme(axis.title.x=element_blank()) +
  theme(plot.title = element_text(hjust = 0.5))
histExtViol

taxdel<-ggplot(correlations, aes(x =IRStat , y = PCTTD)) + 
  geom_bar(stat='identity', wideth =0.9, fill='#f768a1') +
  theme_minimal() +
  ylim(0,60) +
  labs(y = "% tax delinquent", 
       title = 'Properties with Tax Delinquency') +
  theme(axis.title.x=element_blank()) +
  theme(plot.title = element_text(hjust = 0.5))
taxdel


LL<-ggplot(correlations, aes(x =IRStat , y = PCTLL)) + 
  geom_bar(stat='identity', wideth =0.9, fill='#f768a1') +
  theme_minimal() +
  ylim(0,60) +
  labs(y = "% out of state landlords", 
       title = 'Properties with out of state landlords') +
  theme(axis.title.x=element_blank()) +
  theme(plot.title = element_text(hjust = 0.5))
LL

pre19<-ggplot(correlations, aes(x =IRStat , y = PCTPre1900)) + 
  geom_bar(stat='identity', wideth =0.9, fill='#c51b8a') +
  theme_minimal() +
  ylim(0,60) +
  labs(y = "% built before 1900", 
       title = 'Properties built before 1900')+
  theme(axis.title.x=element_blank()) +
  theme(plot.title = element_text(hjust = 0.5))
pre19

MU<-ggplot(correlations, aes(x =IRStat , y = PCTMU)) + 
  geom_bar(stat='identity', wideth =0.9, fill='#c51b8a') +
  theme_minimal() +
  ylim(0,60) +
  labs(y = "% multiple units", 
       title = 'Properties with multiple units')+
  theme(axis.title.x=element_blank()) +
  theme(plot.title = element_text(hjust = 0.5))
MU

grid.arrange(histIntViol, histExtViol, taxdel, LL, pre19, MU, nrow = 3)