---
title: "In-Class Exercise 08"
date: "6 March 2023"
date-modified: "`r Sys.Date()`"
author: "Xinyee How"
format: html
editor: visual
execute: 
  echo: true
  eval: true
  warning: false
---

## Installing packages

```{r}
pacman::p_load(olsrr, corrplot, ggpubr, sf, spdep, GWmodel, tmap, tidyverse, gtsummary)
```

### Importing geospatial data

```{r}
mpsz = st_read(dsn = "data/geospatial", layer = "MP14_SUBZONE_WEB_PL")
```

Setting CRS

```{r}
mpsz_svy21 <- st_transform(mpsz, 3414)
```

```{r}
st_crs(mpsz_svy21)
```

```{r}
st_bbox(mpsz_svy21) #view extent
```

### Importing aspatial data

```{r}
condo_resale = read_csv("data/aspatial/Condo_resale_2015.csv")
```

Inspecting aspatial data

```{r}
glimpse(condo_resale)
```

Inspecting coord column

```{r}
head(condo_resale$LONGITUDE) #see the data in XCOORD column
```

```{r}
head(condo_resale$LATITUDE) #see the data in YCOORD column
```

Display summary statistics

```{r}
summary(condo_resale)
```

## Converting aspatial data frame into a sf object

```{r}
condo_resale.sf <- st_as_sf(condo_resale,
                            coords = c("LONGITUDE", "LATITUDE"),
                            crs=4326) %>%
  st_transform(crs=3414)
```

**Inspecting new sf layer**

```{r}
head(condo_resale.sf)
```

## Exploratory Data Analysis (EDA)

**Plotting the distribution of *SELLING_PRICE***

```{r}
ggplot(data=condo_resale.sf, aes(x=`SELLING_PRICE`)) +
  geom_histogram(bins=20, color="black", fill="light blue")
```

Graph looks right-skewed. Next, we can do a log transformation to reduce the skewness

```{r}
condo_resale.sf <- condo_resale.sf %>%
  mutate(`LOG_SELLING_PRICE` = log(SELLING_PRICE))
```

**Plotting chart**

```{r}
ggplot(data=condo_resale.sf, aes(x=`LOG_SELLING_PRICE`)) +
  geom_histogram(bins=20, color="black", fill="light blue")
```

Looks better.

**Plotting multiple histograms by variables**

```{r}
AREA_SQM <- ggplot(data=condo_resale.sf, aes(x= `AREA_SQM`)) + 
  geom_histogram(bins=20, color="black", fill="light blue")

AGE <- ggplot(data=condo_resale.sf, aes(x= `AGE`)) +
  geom_histogram(bins=20, color="black", fill="light blue")

PROX_CBD <- ggplot(data=condo_resale.sf, aes(x= `PROX_CBD`)) +
  geom_histogram(bins=20, color="black", fill="light blue")

PROX_CHILDCARE <- ggplot(data=condo_resale.sf, aes(x= `PROX_CHILDCARE`)) + 
  geom_histogram(bins=20, color="black", fill="light blue")

PROX_ELDERLYCARE <- ggplot(data=condo_resale.sf, aes(x= `PROX_ELDERLYCARE`)) +
  geom_histogram(bins=20, color="black", fill="light blue")

PROX_URA_GROWTH_AREA <- ggplot(data=condo_resale.sf, 
                               aes(x= `PROX_URA_GROWTH_AREA`)) +
  geom_histogram(bins=20, color="black", fill="light blue")

PROX_HAWKER_MARKET <- ggplot(data=condo_resale.sf, aes(x= `PROX_HAWKER_MARKET`)) +
  geom_histogram(bins=20, color="black", fill="light blue")

PROX_KINDERGARTEN <- ggplot(data=condo_resale.sf, aes(x= `PROX_KINDERGARTEN`)) +
  geom_histogram(bins=20, color="black", fill="light blue")

PROX_MRT <- ggplot(data=condo_resale.sf, aes(x= `PROX_MRT`)) +
  geom_histogram(bins=20, color="black", fill="light blue")

PROX_PARK <- ggplot(data=condo_resale.sf, aes(x= `PROX_PARK`)) +
  geom_histogram(bins=20, color="black", fill="light blue")

PROX_PRIMARY_SCH <- ggplot(data=condo_resale.sf, aes(x= `PROX_PRIMARY_SCH`)) +
  geom_histogram(bins=20, color="black", fill="light blue")

PROX_TOP_PRIMARY_SCH <- ggplot(data=condo_resale.sf, 
                               aes(x= `PROX_TOP_PRIMARY_SCH`)) +
  geom_histogram(bins=20, color="black", fill="light blue")

ggarrange(AREA_SQM, AGE, PROX_CBD, PROX_CHILDCARE, PROX_ELDERLYCARE, 
          PROX_URA_GROWTH_AREA, PROX_HAWKER_MARKET, PROX_KINDERGARTEN, PROX_MRT,
          PROX_PARK, PROX_PRIMARY_SCH, PROX_TOP_PRIMARY_SCH,  
          ncol = 3, nrow = 4)
```

**Drawing Statistical Point Map**

```{r}
tmap_mode("view")
```

```{r}
tmap_options(check.and.fix = TRUE)
tm_shape(mpsz_svy21)+
  tm_polygons() +
tm_shape(condo_resale.sf) +  
  tm_dots(col = "SELLING_PRICE",
          alpha = 0.6,
          style="quantile") +
  tm_view(set.zoom.limits = c(11,14))
```

```{r}
tmap_mode("plot")
```

### Hedonic Pricing Modelling in R

Simple Linear Regression Method

```{r}
condo.slr <- lm(formula=SELLING_PRICE ~ AREA_SQM, data = condo_resale.sf)
```

Examining layer

```{r}
summary(condo.slr)
```

**Visualising and plotting best fit line**

```{r}
ggplot(data=condo_resale.sf,  
       aes(x=`AREA_SQM`, y=`SELLING_PRICE`)) +
  geom_point() +
  geom_smooth(method = lm)
```

### Multiple Linear Regression Method

**Finding correlationships between variables**

we will be using residual for calculation

```{r}
corrplot(cor(condo_resale[, 5:23]), diag = FALSE, order = "AOE",
         tl.pos = "td", tl.cex = 0.5, method = "number", type = "upper")
```

### Building a hedonic pricing model using multiple linear regression method

```{r}
condo.mlr <- lm(formula = SELLING_PRICE ~ AREA_SQM + AGE    + 
                  PROX_CBD + PROX_CHILDCARE + PROX_ELDERLYCARE +
                  PROX_URA_GROWTH_AREA + PROX_HAWKER_MARKET + PROX_KINDERGARTEN + 
                  PROX_MRT  + PROX_PARK + PROX_PRIMARY_SCH + 
                  PROX_TOP_PRIMARY_SCH + PROX_SHOPPING_MALL + PROX_SUPERMARKET + 
                  PROX_BUS_STOP + NO_Of_UNITS + FAMILY_FRIENDLY + FREEHOLD, 
                data=condo_resale.sf)
summary(condo.mlr)
```

### Preparing Publication Quality Table: olsrr method

Tidy way to generate report

```{r}
condo.mlr1 <- lm(formula = SELLING_PRICE ~ AREA_SQM + AGE + 
                   PROX_CBD + PROX_CHILDCARE + PROX_ELDERLYCARE +
                   PROX_URA_GROWTH_AREA + PROX_MRT  + PROX_PARK + 
                   PROX_PRIMARY_SCH + PROX_SHOPPING_MALL    + PROX_BUS_STOP + 
                   NO_Of_UNITS + FAMILY_FRIENDLY + FREEHOLD,
                 data=condo_resale.sf)
ols_regress(condo.mlr1)
```

### Preparing Publication Quality Table: gtsummary method

```{r}
tbl_regression(condo.mlr1, intercept = TRUE)
```

**Preparing report**

```{r}
tbl_regression(condo.mlr1, 
               intercept = TRUE) %>% 
  add_glance_source_note(
    label = list(sigma ~ "\U03C3"),
    include = c(r.squared, adj.r.squared, 
                AIC, statistic,
                p.value, sigma))
```

**Checking for multicolinearity**

```{r}
ols_vif_tol(condo.mlr1)
```

Test for non-linearity

```{r}
ols_plot_resid_fit(condo.mlr1)
```

Test for Normality Assumption

```{r}
ols_plot_resid_hist(condo.mlr1)
```
