Predicting Housing Prices

#Present a correlation matrix

# Numeric variable selecting for correlation matrix 

numericVars <- 
  select_if(st_drop_geometry(Mklg.sf), is.numeric) %>% na.omit()

# Correlation matrix

testRes = cor.mtest(numericVars, conf.level = 0.95)

c2 <- numericVars %>% dplyr::select(Price = price, HouseArea = shape_Area, HouseIncome  = MedHHInc, PctWhite = pctWhite, Beadroom = bedrooms, HouseAge = Age,
                                                HeatedArea  = heatedarea, FullBaths= fullbaths, Neighborhood = neighborhood, treeCanopy, FirePlace =  numfirepla,
                                                  testElementary, raceAf, housingSize, foodService,
                                                  violent, ageDeath, Grocery,
                                                  Parks, college, wifi, daycare, Firestaion, charger,
                                                  buildingAge, internet, schoolAttendance,
                                                  pedSignal,busStop, lots, transit,shopping, goldLine,library,
                                                  row, art, storm, historic2,
                                                  bicycle, drivingCommuters, impSurface, recreation,
                                                  middleProficiency, voters, Landuse , Incidents, request311, pctWorkableAge) %>%
  cor() %>%
  corrplot(sub = "Figure 3-3. Correlation Matrix of Variables", method = 'circle', p.mat = testRes$p, insig='blank', col=COL1('Blues', 10), order = 'FPC', type = 'lower', tl.cex = 0.5, tl.col = greyPalette5[4],  addgrid.col = 'transparent', diag = TRUE, tl.srt = 45, cl.cex = 0.6, cl.offset = 0.2, bg = "transparent", sig.level = 0.1)

#Present 4 home price correlation scatter plots that you think are of interest. I’m going to look for interesting open data that you’ve integrated with the home sale observations.

#log transform

Mklg.sf <- 
  Mklg.sf %>%
  mutate(lg_shape_Area = log(shape_Area+1), lg_Age = log(Age+1), lg_MedHHInc = log(MedHHInc+1))


css1 <-
  ggplot(Mklg.sf, aes(x = price, y = lg_shape_Area))+
  geom_point(colour = greyPalette5[5], size = 0.8) +
  geom_smooth(method = "lm", se = FALSE, colour = bluePalette5[5], size = 1)+
  stat_poly_eq(label.y = 0.95, label.x = 0.05, size = 3, color = bluePalette5[5])+
  labs(subtitle = "Plot1 : logged House Area",
       x = "Price",
       y = "logged House Area") +
  corTheme() 

css2 <-
  ggplot(Mklg.sf, aes(x = price, y = lg_Age))+
  geom_point(colour = greyPalette5[5], size = 0.8) +
  geom_smooth(method = "lm", se = FALSE, colour = bluePalette5[5], size = 1)+
  stat_poly_eq(label.y = 0.95,label.x = 0.05, size = 3, color = bluePalette5[5])+
  labs(subtitle = "Plot2 : logged House Age",
       x = "Price",
       y = "logged House Age") +
  corTheme()

css3 <-
  ggplot(Mklg.sf, aes(x = price, y = bedrooms))+
  geom_point(colour = greyPalette5[5], size = 0.8) +
  geom_smooth(method = "lm", se = FALSE, colour = bluePalette5[5], size = 1) +
  stat_poly_eq(label.y = 0.95, label.x = 0.05, size = 3, color = bluePalette5[5])+
  labs(subtitle = "Plot3 : Number of Bedrooms",
       x = "Price",
       y = "Bedrooms") +
  corTheme()

css4 <-
  ggplot(Mklg.sf, aes(x = price, y = MedHHInc))+
  geom_point(colour = greyPalette5[5], size = 0.8) +
  geom_smooth(method = "lm", se = FALSE, colour = bluePalette5[5], size = 1)+
  stat_poly_eq(label.y = 0.95, label.x = 0.05, size = 3 ,color = bluePalette5[5]) +
  labs(subtitle = "Plot4 : Median Household Income",
       x = "Price",
       y = "Household Income") +
  corTheme()

# Calculate tract's average value of dependent variable

Mklg.sf.DV <- Mklg.sf %>%
  group_by(GEOID) %>%
  summarise(Avg.HP = mean(price),
            AVg.AG = mean(Age),
            Avg.AR = mean(shape_Area)) %>%
  st_drop_geometry() %>%
  mutate_all(~replace(., is.na(.), 0))

Mklg.Nhoods <- Nhoods %>%
  left_join(Mklg.sf.DV, by = "GEOID")

#Develop 1 map of your dependent variable (sale price)

M35 <- ggplot() +
  geom_sf(data = Mklg.Nhoods, aes(fill = q5(Avg.HP)), color = greyPalette5[1], size = 0.2) +
  scale_fill_manual(values = purplePalette5,
                    labels = round(as.numeric(qBr(Mklg.Nhoods, "Avg.HP"))),
                    name = "House Price($)\n(Quintile)",
                    na.translate=FALSE) +
  labs(title = "Map of Mean House Price", 
       subtitle = "Mecklenberg County",
       caption = 'Figure 3-5') +
  mapThememin()

M35

#Develop 3 maps of 3 of your most interesting independent variables.


M352 <- ggplot() +
  geom_sf(data = Mklg.Nhoods, aes(fill = q5(Avg.HP)), color = greyPalette5[1], size = 0.2) +
  scale_fill_manual(values = purplePalette5,
                    labels = round(as.numeric(qBr(Mklg.Nhoods, "Avg.HP"))/1000),
                    name = "House Price\n(1,000$)",
                    na.translate=FALSE) +
  labs(title = "Mean House Price",
       caption = 'Figure 3-6-1')  +
  mapThememin()

#1IDV : House Age
M361 <- ggplot() +
  geom_sf(data = Mklg.Nhoods, aes(fill = rev(q5(AVg.AG))), color = greyPalette5[1], size = 0.2) +
  scale_fill_manual(values = bluePalette5,
                    labels = rev(qBr(Mklg.Nhoods, "AVg.AG")),
                    name = "Age(years)\n(Quintile)",
                    na.translate=FALSE) +
  labs(title = "Mean House Age", 
       caption = 'Figure 3-6-1') +
  mapThememin()


#2IDV : House Area
M362 <- ggplot()+
  geom_sf(data = Mklg.Nhoods, aes(fill = q5(Avg.AR)), color = greyPalette5[1], size = 0.2) +
  scale_fill_manual(values = bluePalette5,
                    labels = round(as.numeric(qBr(Mklg.Nhoods, "Avg.AR"))/100),
                    name = "Area\n(100f²)",
                    na.translate=FALSE) +
  labs(title = "Mean House Shape Area",
       caption = 'Figure 3-6-2') +
  mapThememin()


#3IDV : Median Household Income
M363 <- ggplot() +
  geom_sf(data = Mklg.Nhoods, aes(fill = q5(MedHHInc)), color = greyPalette5[1], size = 0.2) +
  scale_fill_manual(values = bluePalette5,
                    labels = round(as.numeric(qBr(Mklg.Nhoods, "MedHHInc"))/1000),
                    name = "Income\n(1,000$)",
                    na.translate=FALSE) +
  labs(title = "Median Household Income",
       caption = 'Figure 3-6-3') +
  mapThememin()

#4IDV : percent of White
M364 <- ggplot() +
  geom_sf(data = Mklg.Nhoods, aes(fill = q5(pctWhite)), color = greyPalette5[1], size = 0.2) +
  scale_fill_manual(values = bluePalette5,
                    labels = qBr(Mklg.Nhoods, "pctWhite"),
                    name = "White(%)\n(Quintile)",
                    na.translate=FALSE) +
  labs(title = "Percent of White",
       caption = 'Figure 3-6-4') +
  mapThememin()

#5IDV : percent of Vacant Housing
M365 <- ggplot() +
  geom_sf(data = Mklg.Nhoods, aes(fill = q5(pctVacnant)), color = greyPalette5[1], size = 0.2) +
  scale_fill_manual(values = bluePalette5,
                    labels = qBr(Mklg.Nhoods, "pctVacnant"),
                    name = "Vacant(%)\n(Quintile)",
                    na.translate=FALSE) +
  labs(title = "Percent of Vacant Housing",
       caption = 'Figure 3-6-5') +
  mapThememin()

#6IDV : Gas Consumption
M366 <- ggplot() +
  geom_sf(data = Mklg.Nhoods, aes(fill = q5(gasConsumption)), color = greyPalette5[1], size = 0.2) +
  scale_fill_manual(values = bluePalette5,
                    labels = qBr(Mklg.Nhoods, "gasConsumption"),
                    name = "Vacant(%)\n(Quintile)",
                    na.translate=FALSE) +
  labs(title = "Gas Consumption",
       caption = 'Figure 3-6-6') +
  mapThememin()

#7IDV : Crime Violent
M367 <- ggplot() +
  geom_sf(data = Mklg.Nhoods, aes(fill = q5(violent)), color = greyPalette5[1], size = 0.2) +
  scale_fill_manual(values = bluePalette5,
                    labels = qBr(Mklg.Nhoods, "violent"),
                    name = "Vacant(%)\n(Quintile)",
                    na.translate=FALSE) +
  labs(title = "Crime Violent",
       caption = 'Figure 3-6-7') +
  mapThememin()

#8IDV : Elementary Proficiency
M368 <- ggplot() +
  geom_sf(data = Mklg.Nhoods, aes(fill = q5(testElementary)), color = greyPalette5[1], size = 0.2) +
  scale_fill_manual(values = bluePalette5,
                    labels = qBr(Mklg.Nhoods, "testElementary"),
                    name = "Proficiency(%)\n(Quintile)",
                    na.translate=FALSE) +
  labs(title = "Elementary Proficiency",
       caption = 'Figure 3-6-8') +
  mapThememin()

#9IDV : Percent of Poverty
M369 <- ggplot()+
  geom_sf(data = Mklg.Nhoods, aes(fill = q5(pctPoverty)), color = greyPalette5[1], size = 0.2) +
  scale_fill_manual(values = bluePalette5,
                    labels = qBr(Mklg.Nhoods, "pctPoverty"),
                    name = "Poverty(%)\n(Quintile)",
                    na.translate=FALSE) +
  labs(title = "Percent of Poverty",
       caption = 'Figure 3-6-9') +
  mapThememin()

#Split the ‘toPredict’ == “MODELLING” into a training and test set. Provide a polished table of your training set lm summary results (coefficients, R2 etc)

Mklg.Modelling <- Mklg.sf %>% filter(toPredict == "MODELLING")
Mklg.Challenge <- Mklg.sf %>% filter(toPredict == "CHALLENGE")

set.seed(127)

inTrain <- createDataPartition(y = paste(Mklg.Modelling$lot_num), p = 0.8, list = FALSE)

Mklg.training <- Mklg.Modelling[inTrain,]
Mklg.test <- Mklg.Modelling[-inTrain,]


reg.training <- lm(price ~ ., data = st_drop_geometry(Mklg.training) %>% 
                                    dplyr::select(price, shape_Area, MedHHInc, pctWhite, bedrooms, Age,
                                                  heatedarea, fullbaths, neighborhood, treeCanopy, numfirepla,
                                                  testElementary, raceAf, housingSize, foodService,
                                                  violent, ageDeath, Grocery,
                                                  Parks, college, wifi, daycare, Firestaion, charger,
                                                  buildingAge, internet, schoolAttendance,
                                                  pedSignal,busStop, lots, transit,shopping, goldLine,library,
                                                  row, art, storm, historic2,
                                                  bicycle, drivingCommuters, impSurface, recreation,
                                                  middleProficiency, voters, Landuse , Incidents, request311, pctWorkableAge)) # Incidents, request311,

stargazer(reg.training, title = "Table 5-1. Summary Statistics of Training Set", style = "apsr", type = "html", single.row = TRUE,  notes.align = "l", notes.label = "")

#Provide a polished table of mean absolute error and MAPE for a single test set. Check out the “kable” function for markdown to create nice tables

Mklg.test <-
  Mklg.test %>%
  mutate(price.Predict = predict(reg.training, Mklg.test),
         price.Error = price.Predict - price,
         price.AbsError = abs(price.Predict - price),
         price.APE = (abs(price.Predict - price)) / price.Predict)

Mklg.training <-
  Mklg.training %>%
  mutate(price.Predict = predict(reg.training, Mklg.training),
         price.Error = price.Predict - price,
         price.AbsError = abs(price.Predict - price),
         price.APE = (abs(price.Predict - price)) / price.Predict)

Mklg.Challenge <-
  Mklg.Challenge %>%
  mutate(price = predict(reg.training, Mklg.Challenge),
         price.Predict = predict(reg.training, Mklg.Challenge),
         price.Error = price.Predict - price,
         price.AbsError = abs(price.Predict - price),
         price.APE = (abs(price.Predict - price)) / price.Predict)

Mklg.Challenge.csv <- Mklg.Challenge %>% dplyr::select(musaID, price.Predict) %>% st_drop_geometry()

write_csv(Mklg.Challenge.csv, file =  "Thwitches.csv")

R = cor(Mklg.test$price, Mklg.test$price.Predict, method = "pearson", use = "complete.obs")
R2 = R*R
MAE = mean(Mklg.test$price.AbsError, na.rm = T)
MAPE = mean(Mklg.test$price.APE, na.rm = T)

Mklg.test.summary <- as.data.frame(cbind(R2, MAE, MAPE))

colnames(Mklg.test.summary)<-c("R Squre", "MAE", "MAPE")

kable(Mklg.test.summary, caption = "<span style='font-size:14px'><b>Table 5-2. Summary table of Test set results</b></span>", align = "c") %>% 
  kable_minimal(full_width = T, html_font = "Arial", font_size = 14) %>%
  footnote(general = "", fixed_small_size = T, general_title = "")

#Provide the results of your cross-validation tests. This includes mean and standard deviation MAE. Do 100 folds and plot your cross-validation MAE as a histogram. Is your model generalizable to new data?


fitControl <- trainControl(method = "cv", number = 500)

set.seed(123)

reg.cv <- 
  train(price ~ ., data = st_drop_geometry(Mklg.sf) %>% 
                                dplyr::select(price, shape_Area, MedHHInc, pctWhite, bedrooms, Age,
                                                  heatedarea, fullbaths, neighborhood, treeCanopy, numfirepla,
                                                  testElementary, raceAf, housingSize, foodService,
                                                  violent, ageDeath, Grocery,
                                                  Parks, college, wifi, daycare, Firestaion, charger,
                                                  buildingAge, internet, schoolAttendance,
                                                  pedSignal,busStop, lots, transit,shopping, goldLine,library,
                                                  row, art, storm, historic2,
                                                  bicycle, drivingCommuters, impSurface, recreation,
                                                  middleProficiency, voters, Landuse , Incidents, request311, pctWorkableAge), 
     method = "lm", trControl = fitControl, na.action = na.pass)

#reg.cv
#reg.cv$resample[1:5,]
#reg.cv$resample[,3]

MAE.cv = data.table(MAE.cv= reg.cv$resample[,3])

#MAE.cv

MAE.sd <- as.character(round(sd(reg.cv$resample[,3])))
MAE.mean <- as.character(round(mean(reg.cv$resample[,3])))

caption <- paste("Sd =",MAE.sd, "  Mean =",MAE.mean, " K = 500")

ggplot(MAE.cv, aes(x = MAE.cv))+
  geom_histogram(fill=bluePalette5[3], color = greyPalette5[1], size = 2, xlim=c(20,50))+
  labs(title = "Cross-Validation MAE histogram",
       subtitle = caption,
       x="MAE",
       y="Count",
       caption = "Figure 5-3") +
  corTheme2()  

#Plot predicted prices as a function of observed prices

Mklg.training$set <- "Training"
Mklg.test$set <- "Test"
Mklg.Challenge$set <- "challenge"

Mklg.All <-
  rbind(Mklg.training, Mklg.test) 

Mklg.All.Challenge <-
  rbind(Mklg.All ,Mklg.Challenge) 

#head(Mklg.All)

ggplot(Mklg.All, aes(x = price.Predict, y = price))+
  geom_point(size =0.8, alpha = 0.8, color = bluePalette5[3])+
  geom_abline(color = greyPalette5[4])+
  geom_smooth(method = "lm", se = FALSE, colour = bluePalette5[4], size = 1)+
  labs(title = "Predicted prices as a function of Observed prices",
       subtitle = "Grey line : Prefect prediction     Blue line : Average Prediction",
       x="Predicted price",
       y="Observed price",
       caption = "Figure 5-4") +
  corTheme2()  

coords <- st_coordinates(Mklg.sf) 
neighborList <- knn2nb(knearneigh(coords, 5))
spatialWeights <- nb2listw(neighborList, style="W")
Mklg.sf$lagPrice <- lag.listw(spatialWeights, Mklg.sf$price)

coords.test <-  st_coordinates(Mklg.test) 
neighborList.test <- knn2nb(knearneigh(coords.test, 5))
spatialWeights.test <- nb2listw(neighborList.test, style="W")
Mklg.test$lagPrice <- lag.listw(spatialWeights.test, Mklg.test$price)
Mklg.test$lagPriceError <- lag.listw(spatialWeights.test, Mklg.test$price.Error) 


#Spatial lag : Price as a function of the spatial lag of price

SL1 <- ggplot(Mklg.test, aes(x = lagPrice, y = price))+
  geom_point(size =0.8, alpha = 0.8, color = bluePalette5[3])+
  geom_abline(color = greyPalette5[4])+
  geom_smooth(method = "lm", se = FALSE, colour = bluePalette5[4], size = 1)+
  labs(title = "Prices as a function of Spatial Lag Price",
       subtitle = "Grey line : Perfect     Blue line : Average",
       x="Spatial Lag Price",
       y="Observed price",
       caption = "Figure 5-6-1") +
  corTheme3()  

#Spatial lag : Error as a function of the spatial lag of price

SL2 <- ggplot(Mklg.test, aes(x = lagPriceError, y = price.Error))+
  geom_point(size =0.8, alpha = 0.8, color = bluePalette5[3])+
  geom_abline(color = greyPalette5[4])+
  geom_smooth(method = "lm", se = FALSE, colour = bluePalette5[4], size = 1)+
  scale_y_continuous(position = "right")+
  labs(title = "Errors as a function of Spatial Lag Errors",
       subtitle = "Grey line : Perfect     Blue line : Average",
       x="Spatial Lag Error",
       y="Error",
       caption = "Figure 5-6-2") +
  corTheme4()  

moranTest <- moran.mc(Mklg.test$price.Error, 
                      spatialWeights.test, nsim = 999)


MT <- ggplot(as.data.frame(moranTest$res[c(1:999)]), aes(moranTest$res[c(1:999)])) +
  geom_histogram(binwidth = 0.01) +
  geom_vline(aes(xintercept = moranTest$statistic), colour = bluePalette5[5],size=1) +
  scale_x_continuous(limits = c(-1, 1))+
  labs(title = "Observed and permuted Moran's I",
       subtitle = "Observed Moran's I = 0.1958 : Blueline",
       x="Moran's I",
       y="Count",
       caption = "Figure 5-6-3")  +
  corTheme2()  

SL1 + SL2

#Provide a map of your predicted values for where ‘toPredict’ is both “MODELLING” and “CHALLENGE”

PP <- ggplot() +
  geom_sf(data = Nhoods, fill = greyPalette5[2], size = 0.3, color = greyPalette5[1]) +
  geom_sf(data = Mklg.All.Challenge, aes(colour = q5(price.Predict)), alpha = 0.4, size = 0.5, show.legend = "point")+
  labs(title = "Map of Predicted House Price", 
       subtitle = "Entire data set",
       caption = 'Figure 5-7-1')+
  scale_colour_manual(values = bluePalette5,
                      labels=round(as.numeric(qBr(Mklg.All,"price.Predict"))/1000),
                      name="Predicted\nPrice\n(1,000$)",
                      na.translate=FALSE)+
  mapThememin()

OP <- ggplot() +
  geom_sf(data = Nhoods, fill = greyPalette5[2], size = 0.3, color = greyPalette5[1]) +
  geom_sf(data = Mklg.All.Challenge, aes(colour = q5(price)), alpha = 0.4, size = 0.5, show.legend = "point")+
  labs(title = "Map of Observed House Price", 
       subtitle = "Entire data set",
       caption = 'Figure 5-7-2')+
  scale_colour_manual(values = purplePalette5,
                      labels=round(as.numeric(qBr(Mklg.All,"price"))/1000),
                      name="Observed\nPrice\n(1,000$)",
                      na.translate=FALSE)+
  mapThememin()

NT <- ggplot() +
  geom_sf(data = Nhoods, fill = greyPalette5[2], size = 0.3, color = greyPalette5[1]) +
  geom_sf(data = Mklg.All.Challenge %>% filter(price.Predict < 0) , aes(colour = q5(price)), alpha = 0.4, size = 0.5, show.legend = "point")+
  labs(title = "Map of Observed House Price", 
       subtitle = "Entire data set",
       caption = 'Figure 5-7-2')+
  scale_colour_manual(values = purplePalette5,
                      labels=round(as.numeric(qBr(Mklg.All,"price"))/1000),
                      name="Observed\nPrice\n(1,000$)",
                      na.translate=FALSE)+
  mapThememin()

#NT

OP + PP + plot_layout(nrow = 1)

#Using the test set predictions, provide a map of mean absolute percentage error (MAPE) by neighborhood

Mklg.test.MAPE <- Mklg.test %>%
  group_by(GEOID) %>%
  summarise(meanMAPE = mean(price.APE), meanPrice = mean(price)) %>%
  st_drop_geometry()

Mklg.test.MAPE <- Nhoods %>% 
  left_join(Mklg.test.MAPE, by = "GEOID")

ggplot() +
  geom_sf(data = Mklg.test.MAPE, aes(fill = q5(meanMAPE)), color = greyPalette5[1], size = 0.2) +
  scale_fill_manual(values = bluePalette5 ,
                    labels = round(as.numeric(qBr2(Mklg.test.MAPE, "meanMAPE"))*100)/100,
                    name = "MAPE\n(Quintile)",
                    na.translate=FALSE) +
  labs(title = "Map of MAPE by neighborhood", 
       subtitle = "Test data set",
       caption = 'Figure 5-8') +
  mapThememin()

#Provide a scatterplot plot of MAPE by neighborhood as a function of mean price by neighborhood

ggplot(Mklg.test.MAPE %>% filter(meanMAPE < 2000) %>% filter(meanMAPE > -80), aes(x = meanMAPE, y = meanPrice))+
  geom_point(size=1)+
  geom_smooth(method = "lm", se = FALSE, colour = bluePalette5[4], size = 1)+
  labs(title = "Scatterplot of MAPE as a function of mean price",
       subtitle = "by neighborhood",
       x="MAPE",
       y="price",
       caption = "Figure 5-9") +
  corTheme2()  

#Using tidycensus, split your study area into two groups (perhaps by race or income) and test your model’s generalizability. Is your model generalizable

Mklg.Nhoods.split <- Mklg.Nhoods %>%
  mutate(raceContext = ifelse(pctWhite > 50, "White", "Non-White"),
         incomeContext = ifelse(MedHHInc > mean(MedHHInc,na.rm = T), "High", "Low"))

CT1 <- ggplot() +
  geom_sf(data = Mklg.Nhoods.split, aes(fill = raceContext), color = greyPalette5[2], size = 0.2) +
  scale_fill_manual(values = c(bluePalette5[1],bluePalette5[4]),
                    name = "RaceContext",
                    na.translate=FALSE) +
  labs(title = "Race Context", 
       subtitle = "Majority White = Percent of White > 50%",
       caption = 'Figure 5-10-2') +
  mapThememin()

CT2 <- ggplot() +
  geom_sf(data = Mklg.Nhoods.split, aes(fill = incomeContext), color = greyPalette5[2], size = 0.2) +
  scale_fill_manual(values = c(bluePalette5[4],bluePalette5[1]),
                    name = "IncomeContext",
                    na.translate=FALSE) +
  labs(title = "Income Context", 
       subtitle = "High Income = Income > Median Income",
       caption = 'Figure 5-10-3') +
  mapThememin()

PP + labs(caption = 'Figure 5-10-1')