pacman::p_load(httr,jsonlite,sf, spdep, progress, tmap,
tidyverse,knitr, spdep,GWmodel, rsample,
SpatialML,Metrics,kableExtra, randomForestExplainer,
performance,lubridate,colorspace,ggdist, ggridges, ggthemes,
plotly,treemap,ggrepel,scales)Take-home Exercise 3 - Modelling Spatial Heterogeneity in HDB Resale Prices via GWR
1 Overview
Housing represents one of the most significant forms of wealth in Singapore, and the resale market for Housing and Development Board (HDB) flats is a critical barometer of housing affordability and accessibility. The resale price of an HDB flat is influenced by a complex interplay of structural characteristics (e.g., floor area, remaining lease, floor level) and locational amenities (e.g., proximity to public transport, schools, and parks).
Traditional global regression models, such as Ordinary Least Squares (OLS), assume that the relationships between these predictors and price are constant across space. However, housing markets are inherently spatial, and this assumption of spatial stationarity is often violated. To address this, our study employs Geographically Weighted Regression (GWR), a local statistical technique that allows us to model how these relationships vary geographically across Singapore, providing a more nuanced and spatially explicit explanation of HDB resale prices.
2 Objective
This study aims to calibrate a spatially explicit explanatory model to identify and evaluate the factors influencing 4-room HDB resale flat prices in Singapore for the period January to September 2025. Moving beyond global averages, the analysis employs Geographically Weighted Regression (GWR) to investigate the spatial heterogeneity of these factors. Specifically, this research seeks to address three core questions:
Which structural attributes and locational amenities exert a statistically significant influence on 4-room HDB resale prices?
How do the magnitudes, directions, and significances of these relationships exhibit spatial non-stationarity across Singapore’s planning regions?
To what extent does the geographically weighted modelling framework improve explanatory power and model fit over a traditional, aspatial Ordinary Least Squares (OLS) benchmark?
3 Dataset
The study integrates both aspatial and geospatial datasets to model HDB resale prices in Singapore. The variables were selected based on hedonic price theory and prior literature on Singapore’s housing market, categorized into structural attributes and locational amenities that theoretically influence property values.
| Dataset | Description | Format | Source |
|---|---|---|---|
| HDB Resale Prices | Resale flat prices based on registration date from Jan 2017 onwards | CSV | data.gov.sg |
| Singapore Master Plan 2019 | URA 2019 Master Plan Planning Subzone boundary data | KML | data.gov.sg |
| Bus Stop Locations | Land Transport Authority bus stop locations across Singapore | SHP | LTA DataMall |
| MRT/LRT Station Exits | Locations of all MRT and LRT station exits in Singapore | KML | data.gov.sg |
| Hawker Centres | Directory of hawker centres managed by National Environment Agency | KML | data.gov.sg |
| Supermarkets | Locations of supermarkets regulated by Singapore Food Agency | KML | data.gov.sg |
| Park Facilities | Representative locations of facilities in parks under National Parks Board | GEOJSON | data.gov.sg |
| Eldercare Services | Directory of eldercare services provided by Ministry of Health | KML | data.gov.sg |
| Childcare Centres | Locations of licensed childcare centers in Singapore | KML | data.gov.sg |
| Kindergartens | Locations of registered kindergartens in Singapore | KML | data.gov.sg |
| Primary Schools | List of all primary schools in Singapore from Ministry of Education | CSV | MOE |
| Shopping Malls | Comprehensive list of shopping malls in Singapore with coordinates | CSV | Wikipedia |
| CHAS Clinics | Locations of Community Health Assist Scheme clinics | KML | data.gov.sg |
| Community Clubs | Locations of People’s Association community clubs | KML | data.gov.sg |
| Sports Facilities | ActiveSG sports centers, gyms, swimming pools, and stadiums | CSV | ActiveSG |
| Hospitals | List of general hospitals in Singapore | CSV | Wikipedia |
| HDB HIP/MUP Status | Manual compilation of Home Improvement Programme and Main Upgrading Programme status for HDB blocks | CSV | HDB Portal |
Proximity measures will be calculated as Euclidean distances from each HDB block to the nearest facility, while accessibility will be measured as count of facilities within specified buffer zones (e.g., 350m for kindergartens, 1km for primary schools) based on common walking distances.
4 Installing and Loading the R packages
The packages used for this analysis are outlined below. This study primarily relies on the sf ecosystem for geospatial data handling and GWmodel for geographically weighted regression analysis.
| Package | Description |
|---|---|
| sf | Provides functions to manage, process, and manipulate Simple Features, a formal geospatial data standard for storing and accessing spatial geometries like points, lines, and polygons. |
| tidyverse | A collection of packages for data science tasks such as importing, tidying, wrangling, and visualizing data. |
| tmap | Enables cartographic-quality static maps and interactive maps using the leaflet API. |
| GWmodel | Provides tools for calibrating geographically weighted models, including Geographically Weighted Regression (GWR). |
| olsrr | Facilitates building OLS regression models and performing diagnostic tests. |
| ggplot2 | A powerful graphics system for creating elegant plots based on the Grammar of Graphics. |
| spdep | Offers functions to create spatial weights matrices from polygon contiguities and test for spatial autocorrelation. |
| gtsummary | Creates publication-ready summary tables for regression models with elegant and flexible formatting. |
| sfdep | Builds on spdep with a tidyverse-friendly interface, list-columns, and enhanced spatial dependence analysis tools. |
| httr | Provides tools for working with HTTP requests, useful for accessing web APIs like OneMap SG for geocoding. |
| jsonlite | A fast and flexible JSON parser and generator optimized for statistical and web data. |
| corrplot | Enables visual exploration of correlation matrices using various plotting methods. |
| performance | Offers comprehensive tools for assessing regression model performance, including R², RMSE, and diagnostics. |
5 Data Import and Preperation
Before importing the data, we must establish our spatial framework and core dataset. We begin by loading the Master Plan 2019 Subzone (MPSZ) boundaries, which provide the essential geographical context for Singapore’s planning regions. Simultaneously, we import the raw HDB resale transaction data. The initial step involves filtering this dataset to focus exclusively on 4-room flats and the relevant time period (Jan-Sep 2025) as defined by our objective. This filtering is critical to ensure our model is calibrated on a homogeneous and temporally consistent property type.
The code reads the MPSZ polygon data and checks its Coordinate Reference System (CRS) to ensure it uses the local projected system, SVY21 (EPSG:3414), for accurate spatial measurements. It then imports the CSV file, filters the rows based on flat_type and month, and provides a summary of the data structure to confirm the successful load and filtering.
5.1 MPSZ19 Geospatial Data
mpsz19 <- read_rds("data/mpsz.rds")
st_crs(mpsz19)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]5.2 HDB Resale Data
As we are building a predictive model for 4‑Room HDB resale prices, the code below filters for flats of type 4‑Room and resale records between July 2024 and September 2025
HDBresale <- read_csv("data/rawdata/ResaleflatpricesbasedonregistrationdatefromJan2017onwards.csv")%>%
filter(flat_type=="4 ROOM")%>%
filter(month >="2025-01" & month < "2025-10")Rows: 218966 Columns: 11
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (8): month, town, flat_type, block, street_name, storey_range, flat_mode...
dbl (3): floor_area_sqm, lease_commence_date, resale_price
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Let’s check the data structure:
kable(str(HDBresale))spc_tbl_ [8,679 × 11] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
$ month : chr [1:8679] "2025-01" "2025-01" "2025-01" "2025-02" ...
$ town : chr [1:8679] "ANG MO KIO" "ANG MO KIO" "ANG MO KIO" "ANG MO KIO" ...
$ flat_type : chr [1:8679] "4 ROOM" "4 ROOM" "4 ROOM" "4 ROOM" ...
$ block : chr [1:8679] "337" "310C" "337" "226" ...
$ street_name : chr [1:8679] "ANG MO KIO AVE 1" "ANG MO KIO AVE 1" "ANG MO KIO AVE 1" "ANG MO KIO AVE 1" ...
$ storey_range : chr [1:8679] "13 TO 15" "10 TO 12" "10 TO 12" "07 TO 09" ...
$ floor_area_sqm : num [1:8679] 91 94 91 91 96 92 98 97 94 87 ...
$ flat_model : chr [1:8679] "New Generation" "Model A" "New Generation" "New Generation" ...
$ lease_commence_date: num [1:8679] 1982 2012 1982 1978 2012 ...
$ remaining_lease : chr [1:8679] "56 years 03 months" "86 years 09 months" "56 years 02 months" "52 years" ...
$ resale_price : num [1:8679] 582000 865000 530000 560000 1038000 ...
- attr(*, "spec")=
.. cols(
.. month = col_character(),
.. town = col_character(),
.. flat_type = col_character(),
.. block = col_character(),
.. street_name = col_character(),
.. storey_range = col_character(),
.. floor_area_sqm = col_double(),
.. flat_model = col_character(),
.. lease_commence_date = col_double(),
.. remaining_lease = col_character(),
.. resale_price = col_double()
.. )
- attr(*, "problems")=<externalptr> After filtering, there are total 8,679 observations for the study.
The code below check wether data has missing values:
any(is.na(HDBresale))[1] FALSEThe output, specifically the st_crs(mpsz19) result, must be verified to confirm it is indeed EPSG:3414. The str() output of the HDB data should be checked to ensure the filtering worked correctly and that there are no obvious data import errors, such as incorrect column types or a unexpected number of observations.
5.2.1 Geocode HDB resale data
The raw HDB resale data contains valuable information like block numbers and street names, but it lacks the precise geographic coordinates (latitude and longitude) required for any spatial analysis. To transform this from a simple table into a mappable, spatially-aware dataset for Geographically Weighted Regression, we must convert these text addresses into precise locations on the map. This process, called geocoding, is the foundational step that enables all subsequent distance calculations and spatial modeling.
This code performs two main operations. First, it cleans the street names by replacing “St.” with “SAINT” using gsub() to prevent geocoding errors where the abbreviation might be misinterpreted. Second, it defines a custom geocode() function that queries Singapore’s OneMap API for each “Block + Street Name” combination. The function sends an HTTP request, parses the JSON response, and safely extracts the coordinates. A loop then iterates through all 8,679 HDB transactions, calling this function for each address and storing the results in new LATITUDE and LONGITUDE columns. Finally, the completed dataset is saved to disk and reloaded to ensure data persistence.
HDBresale$street_name <- gsub("ST\\.","SAINT",HDBresale$street_name) Next, since we need the exact longitude and latitude for further geospatial modeling, we use OneMap’s API to retrieve these coordinates for each observation:
geocode <- function(block, streetname) {
base_url <- "https://www.onemap.gov.sg/api/common/elastic/search"
address <- paste(block, streetname)
query <- list(
searchVal = address,
returnGeom = "Y",
getAddrDetails = "N",
pageNum = "1"
)
res <- GET(base_url, query = query)
restext <- content(res, as = "text", encoding = "UTF-8")
json <- fromJSON(restext, simplifyVector = FALSE)
# Handle missing results safely
if (is.null(json$results) || length(json$results) == 0) {
return(data.frame(LATITUDE = NA, LONGITUDE = NA))
}
result <- json$results[[1]]
# extract LATITUDE and LONGITUDE w/o warning messages
lat <- suppressWarnings(as.numeric(result$LATITUDE))
lon <- suppressWarnings(as.numeric(result$LONGITUDE))
return(data.frame(LATITUDE = lat, LONGITUDE = lon))
}HDBresale$LATITUDE <- 0
HDBresale$LONGTITUDE <- 0
for (i in seq_len(nrow(HDBresale))) {
if (HDBresale$LATITUDE[i] == 0 ||
HDBresale$LONGTITUDE[i] == 0) {
temp_output <- geocode(HDBresale[i, 4], HDBresale[i, 5])
HDBresale$LATITUDE[i] <- temp_output$LATITUDE
HDBresale$LONGTITUDE[i] <- temp_output$LONGITUDE
}
}write_rds(HDBresale,"data/HDBresale.rds")HDBresale <- read_rds("data/HDBresale.rds") HDBresale <- HDBresale%>%
filter(month >="2025-01" & month < "2025-10")The successful execution of this code creates a significantly enhanced dataset. The new HDBresale dataframe now contains numeric LATITUDE and LONGITUDE columns for every transaction record. This transformation from textual addresses to precise coordinates means we now have a georeferenced dataset ready for conversion to a spatial format. The absence of error messages indicates that the OneMap API handled most queries successfully, though some records might contain NA values if addresses couldn’t be matched. The final filtered dataset containing only 2025 transactions with coordinates forms the essential foundation for all subsequent spatial operations, distance calculations, and ultimately, the Geographically Weighted Regression analysis.
5.2.2 Prepare Data
This section focuses on the crucial data preparation phase for the structural attributes that fundamentally influence HDB resale prices. The raw dataset contains variables in a format unsuitable for direct statistical modeling, requiring transformation into meaningful numerical and categorical predictors. Specifically, the textual storey_range (e.g., “10 TO 12”) must be converted into an ordinal variable to capture the non-linear price premium associated with higher floors. Similarly, the lease_commence_date needs to be used to calculate the unit_age, and the text-based remaining_lease must be parsed into a numerical value representing years. By engineering these variables—floor_level, unit_age, and a numerical remaining_lease—we create a robust analytical foundation that allows the subsequent Geographically Weighted Regression model to accurately quantify how these intrinsic structural characteristics, from floor height to lease decay, contribute to a property’s value across different locations in Singapore.
Storey Range
Before getting started, let’s visualize the storey range first:
ggplot(HDBresale, aes(x = storey_range)) +
geom_bar(fill = "#93A4CC", color = "grey30") +
geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5, size = 3, color = "grey30") +
labs(title = "Distribution of HDB Storey Ranges",
x = "Storey Range",
y = "Count") +
theme_minimal() +
theme(
plot.title = element_text(size = 14, hjust = 0.5,
margin = ggplot2::margin(b = 15)),
panel.background = element_rect(fill = "#f6f6f6", color = NA),
plot.background = element_rect(fill = "#f6f6f6", color = NA),
panel.grid = element_blank(),
axis.ticks = element_line(color = "grey30"),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.x = element_line(color = "grey30"),
axis.text.x = element_text(angle = 90, hjust = 1, size = 8,
margin = ggplot2::margin(t = 10)),
axis.title.x = element_text(size = 10,
margin = ggplot2::margin(t = 10))
)Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(count)` instead.
The distribution of storey ranges for 4-room HDB units exhibits a pronounced right skew, revealing a clear market concentration in lower to mid-level floors. Specifically, over three-quarters (77.5%) of all transactions occur in units situated between the 1st and 12th floors. To transform this textual data into a robust analytical variable for our model, we adopt a five-tier classification system (Very Low to High), grounded in established industry standards from PropertyGuru and the Singapore Public Housing Wiki. This categorization not only mitigates potential overfitting by reducing granularity but also captures the expected non-linear price premiums associated with height, where the value difference between floors 2 and 3 may be negligible compared to that between floors 21 and 22. The final classification—Very Low (1-3), Low (4-6), Mid (7-9), Mid-High (10-21), and High (22+)—strategically aligns common HDB conventions with the specific distribution observed in our dataset, ensuring our floor_level variable is both interpretable and statistically sound for spatial analysis.
The code below extracts the maximum floor level for each storey range and label it based on classification above:
HDBresale$storey_max <- sapply(strsplit(HDBresale$storey_range, " TO "), function(x) max(as.numeric(x)))
HDBresale$floor_level <- cut(
HDBresale$storey_max,
breaks = c(0, 3, 6, 9, 21, Inf),
labels = c(1, 2, 3, 4, 5),
include.lowest = TRUE
)
table(HDBresale$floor_level)
1 2 3 4 5
1429 1832 1788 3213 417 Next, we plot the graph again:
ggplot(HDBresale, aes(x = floor_level)) +
geom_bar(fill = "#93A4CC", color = "grey30") +
geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5, size = 3, color = "grey30") +
labs(title = "Distribution of HDB Floor Level",
x = "Floor Level",
y = "Count") +
theme_minimal() +
theme(
plot.title = element_text(size = 14, hjust = 0.5,
margin = ggplot2::margin(b = 15)),
panel.background = element_rect(fill = "#f6f6f6", color = NA),
plot.background = element_rect(fill = "#f6f6f6", color = NA),
panel.grid = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.ticks = element_line(color = "grey30"),
axis.line.x = element_line(color = "grey30"),
axis.text.x = element_text(hjust = 0.5, size = 8,
margin = ggplot2::margin(t = 10)),
axis.title.x = element_text(size = 10,
margin = ggplot2::margin(t = 10))
)
After combining the storey ranges, we can see that mid-high floors (10–21) have the highest count, with 5,406 units
Unit Age
HDB’s Minimum Occupation Period (MOP) policy regulates how long a HDB buyer must live in the flat before allowed to sell or rent it out entirely. Here’s the summary:
| Mode of Purchase | MOP Duration | Notes |
|---|---|---|
| BTO / SBF (Standard flats) | 5 years | Typical case |
| Plus / Prime flats | 10 years | For highly desirable locations |
| DBSS flats | 5 years | Built by private developers |
| SERS (before 7 Apr 2022) | 7 years from the date of selection 5 years from the date of collecting the keys |
Transitional rule |
| SERS (after 7 Apr 2022) | 5 years | Same as standard flats |
| Resale flat (Standard/Unclassified) | 5 years | |
| Resale flat (Plus / Prime) | 10 years |
Next, we calculate the age of each unit by subtracting the lease_commence_date from the year of the resale date and visualize the distribution:
HDBresale <- HDBresale%>%
mutate(resale_date = as.Date(paste0(month, "-01")),
unit_age = year(resale_date)-lease_commence_date)
summary(HDBresale$unit_age) Min. 1st Qu. Median Mean 3rd Qu. Max.
4.00 9.00 23.00 22.85 37.00 58.00 This code performs two critical data transformations using dplyr::mutate(). First, it converts the month character variable (e.g., “2025-01”) into a proper Date object (resale_date) by appending “-01” to represent the first day of each transaction month. Second, it calculates unit_age by subtracting the lease_commence_date from the year component of the resale_date, giving us the exact age of each flat in years at the time of its resale transaction. The summary() function then provides a statistical overview of this newly created variable.
The summary indicates that the minimum unit age is 4 years. To verify whether any flats were sold before completing the 5-year Minimum Occupation Period (MOP), the code checks for cases with lease commencement dates from 2022 onward:
kable(HDBresale %>% filter(unit_age<5 & lease_commence_date>=2022))| month | town | flat_type | block | street_name | storey_range | floor_area_sqm | flat_model | lease_commence_date | remaining_lease | resale_price | LATITUDE | LONGTITUDE | storey_max | floor_level | resale_date | unit_age |
|---|
The validation check confirms full compliance with HDB’s Minimum Occupation Period (MOP) policy, as the query returned no results (NA). This indicates that all transactions from 2022 onward involve flats that are at least 5 years old, validating the dataset’s integrity for our analysis period.
To systematically examine the age distribution of resale units, the visualization below employs a histogram with consistent parameters: the x-axis is bounded between 0 and 60 years to capture the full lease cycle, while a 5-year bin width provides optimal granularity for identifying cohort patterns in the housing stock.
ggplot(HDBresale, aes(x = unit_age)) +
geom_histogram(boundary = 0, binwidth = 5,fill = "#93A4CC", color = "grey30") +
geom_text(stat = 'bin', binwidth = 5, boundary = 0,
aes(label = ..count..), vjust = -0.7, size = 3, color = "grey30") +
labs(title = "Distribution of HDB Unit Age",
x = "Unit Age (years)",
y = "Count") +
scale_x_continuous(
limits = c(0, 60),
breaks = seq(0, 60, by = 5)
) +
theme_minimal() +
theme(
plot.title = element_text(size = 14, hjust = 0.5,
margin = ggplot2::margin(b = 15)),
panel.background = element_rect(fill = "#f6f6f6", color = NA),
plot.background = element_rect(fill = "#f6f6f6", color = NA),
panel.grid = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.x = element_line(color = "grey30"),
axis.ticks = element_line(color = "grey30"),
axis.text.x = element_text(hjust = 0.5, size = 8,
margin = ggplot2::margin(t = 10)),
axis.title.x = element_text(size = 10,
margin = ggplot2::margin(t = 10))
)
The distribution of unit age for resale HDBs is multimodal, with peaks at 5–10 years (1,907 units), 25–30 years (1,089 units), and 35–40 years (1,397 units).
Remaining Lease
The current remaining_lease attribute is stored as character text, which prevents quantitative analysis of how lease decay influences resale prices. To enable statistical modeling, we transform this textual data into a continuous numerical variable representing the remaining lease duration in decimal years. The custom function below systematically parses the character strings, extracting both years and months components before converting them into a standardized numeric format suitable for regression analysis.
parse_remaining_lease <- function(x) {
years <- str_extract(x, "(\\d+)\\s*years?") %>% str_extract("\\d+") %>% as.numeric()
months <- str_extract(x, "(\\d+)\\s*months?") %>% str_extract("\\d+") %>% as.numeric()
if (is.na(years)) years <- 0
if (is.na(months)) months <- 0
years + months/12
}
HDBresale <- HDBresale %>%
mutate(remaining_lease = round(sapply(remaining_lease, parse_remaining_lease), 1))Price Per sqm
HDBresale <- HDBresale %>%
mutate(price_per_sqm = round(resale_price / floor_area_sqm, 1))5.2.3 Convert Data into sf Format
Following the geocoding process, we transform the tabular HDB data into a formal spatial dataset using the st_as_sf() function, which creates a Simple Features (sf) dataframe. This conversion explicitly incorporates the latitude and longitude coordinates as geometric point data, establishing the fundamental geographic structure required for all subsequent spatial analysis and geographically weighted modeling.
HDBresale_sf <- st_as_sf(HDBresale,
coords = c("LONGTITUDE", "LATITUDE"),
crs = 4326)
HDBresale_sfSimple feature collection with 8679 features and 16 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.6859 ymin: 1.272099 xmax: 103.9649 ymax: 1.457071
Geodetic CRS: WGS 84
# A tibble: 8,679 × 17
month town flat_type block street_name storey_range floor_area_sqm
* <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
1 2025-01 ANG MO KIO 4 ROOM 337 ANG MO KIO AV… 13 TO 15 91
2 2025-01 ANG MO KIO 4 ROOM 310C ANG MO KIO AV… 10 TO 12 94
3 2025-01 ANG MO KIO 4 ROOM 337 ANG MO KIO AV… 10 TO 12 91
4 2025-02 ANG MO KIO 4 ROOM 226 ANG MO KIO AV… 07 TO 09 91
5 2025-02 ANG MO KIO 4 ROOM 310C ANG MO KIO AV… 28 TO 30 96
6 2025-02 ANG MO KIO 4 ROOM 226 ANG MO KIO AV… 04 TO 06 92
7 2025-03 ANG MO KIO 4 ROOM 304 ANG MO KIO AV… 10 TO 12 98
8 2025-03 ANG MO KIO 4 ROOM 303 ANG MO KIO AV… 04 TO 06 97
9 2025-03 ANG MO KIO 4 ROOM 310B ANG MO KIO AV… 16 TO 18 94
10 2025-03 ANG MO KIO 4 ROOM 310B ANG MO KIO AV… 16 TO 18 87
# ℹ 8,669 more rows
# ℹ 10 more variables: flat_model <chr>, lease_commence_date <dbl>,
# remaining_lease <dbl>, resale_price <dbl>, storey_max <dbl>,
# floor_level <fct>, resale_date <date>, unit_age <dbl>, price_per_sqm <dbl>,
# geometry <POINT [°]>Next, we transform the sf data from WGS84 coordinate system to SVY21 for modelling purpose:
HDBresale_sf <- HDBresale_sf %>%
st_transform(3414)After transformation, let’s quickly plot the resale HDB distribution using tmap:
tm_shape(mpsz19)+
tm_polygons(fill="REGION_N",fill_alpha=0.5, lwd =0.2)+
tm_shape(HDBresale_sf)+
tm_dots(size = 0.05)+
tm_layout(bg.color = "#f6f6f6",outer.bg.color = "#f6f6f6",frame = FALSE)+
tm_legend(frame = FALSE,bg.color = "#f6f6f6")+
tm_title("HDB Resale Unit Locations",
position = tm_pos_out("center", "top", "center"))── tmap v3 code detected ───────────────────────────────────────────────────────[v3->v4] `tm_legend()`: use 'tm_legend()' inside a layer function, e.g.
'tm_polygons(..., fill.legend = tm_legend())'
This message is displayed once every 8 hours.
5.3 Housing Price Factor Data
To operationalize the locational determinants of HDB resale values, we now enrich our spatial dataset with proximity and accessibility metrics to key urban amenities. This involves calculating both Euclidean distances to major facilities and density counts within pedestrian catchment areas, capturing fundamental aspects of location utility. Specifically, we compute: straight-line distances to the CBD, MRT stations, hawker centres, shopping malls, supermarkets, parks, and primary schools; alongside counts of kindergartens, childcare centres, and bus stops within a 350-meter walking radius, and primary schools within a 1-kilometer catchment area.
5.3.1 Proximity to CBD
Centrality remains a cornerstone of urban economic theory, with proximity to the Central Business District (CBD) historically commanding significant price premiums. Following URA’s planning framework and established geographical definitions, we delineate the CBD boundary by amalgamating eight contiguous subzones within the Downtown Core planning area—Anson, Cecil, Clifford Pier, Maxwell, Phillip, Raffles Place, Tanjong Pagar, and Marina Centre. Using spatial operations, we merge these subzones into a single polygon representing the CBD and subsequently calculate the Euclidean distance from each HDB transaction point to this unified boundary, creating our primary measure of metropolitan accessibility.
CBD <- mpsz19 %>%
filter(SUBZONE_N %in% c("ANSON", "CECIL", "CLIFFORD PIER", "MAXWELL",
"PHILLIP", "RAFFLES PLACE", "TANJONG PAGAR", "MARINA CENTRE"))%>%
st_union()HDBresale_sf <- HDBresale_sf %>%
mutate(
dist_to_CBD = st_distance(., CBD) %>%
as.numeric() %>%
round(1)
)The CBD area is shown in green on the map:
tm_shape(mpsz19) +
tm_polygons(col = "lightgrey", alpha = 0.2, border.col = "grey50") +
tm_shape(CBD) +
tm_polygons(col = "#8EB859", alpha = 0.6, border.col = "#8EB859") +
tm_shape(HDBresale_sf) +
tm_dots(size = 0.05) +
tm_layout(bg.color = "#f6f6f6",
outer.bg.color = "#f6f6f6",
frame = FALSE) +
tm_legend(frame = FALSE, bg.color = "#f6f6f6") +
tm_title("HDB Resale Unit Locations x CBD",
position = tm_pos_out("center", "top", "center"))── tmap v3 code detected ───────────────────────────────────────────────────────[v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
(instead of 'col'), and 'col' for the outlines (instead of 'border.col').
[v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.
[v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.
This message is displayed once every 8 hours.
5.3.2 Proximity to MRT / Count of MRT within 600m
Sys.setenv(OGR_GEOMETRY_ACCEPT_UNCLOSED_RING = "NO")
MRT <- st_read(dsn = "data/rawdata/TrainStation_Aug2025",
layer = "RapidTransitSystemStation") %>%
st_buffer(0) %>%
st_transform(3414)%>%
filter(!st_is_empty(geometry)) %>%
st_centroid() Reading layer `RapidTransitSystemStation' from data source
`D:\ssinha8752\ISSS608-VAA\Take-home_Ex\Take-home_Ex03\data\rawdata\TrainStation_Aug2025'
using driver `ESRI Shapefile'Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
Error 1: Non closed ring detected.Simple feature collection with 231 features and 7 fields (with 1 geometry empty)
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 6068.209 ymin: 27478.44 xmax: 45377.5 ymax: 47913.58
Projected CRS: SVY21Warning: st_centroid assumes attributes are constant over geometriesdist_matrix <- st_distance(HDBresale_sf, MRT)
within_600m <- st_is_within_distance(HDBresale_sf, MRT, dist=600)
HDBresale_sf <- HDBresale_sf %>%
mutate(
dist_to_MRT = apply(dist_matrix, 1, min) %>%
as.numeric() %>%
round(1),
MRT_600m = map_int(within_600m, length)
)tm_shape(mpsz19) +
tm_polygons(col = "lightgrey", fill_alpha = 0.2, border.col = "grey50") +
tm_shape(HDBresale_sf) +
tm_dots(size = 0.1,alpha=0.4) +
tm_shape(MRT) +
tm_dots(size = 0.5, fill="orange",shape =18) +
tm_layout(bg.color = "#f6f6f6",
outer.bg.color = "#f6f6f6",
frame = FALSE) +
tm_title("HDB Resale Unit Locations x MRT Stations",
position = tm_pos_out("center", "top", "center"))── tmap v3 code detected ───────────────────────────────────────────────────────[v3->v4] `tm_dots()`: use `fill_alpha` instead of `alpha`.
5.3.3 Top 50 Primary School
primary_school_rank <- read_csv("data/rawdata/PrimarySchool_Rank.csv")Rows: 91 Columns: 10
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (6): School, GEP, SAP, Gender, Area, Admission
dbl (4): Rank (2024), Vacancy, Applicant, Popularity
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.primary_school <- read_csv("data/rawdata/Generalinformationofschools.csv") %>%
left_join(primary_school_rank, by = c("school_name" ="School"))%>%
select(c("school_name","postal_code","address",`Rank (2024)`))%>%
filter(`Rank (2024)`<= 50)%>%
arrange(`Rank (2024)`)Rows: 337 Columns: 31
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (30): school_name, url_address, address, telephone_no, telephone_no_2, f...
dbl (1): postal_code
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.geocode <- function(address) {
base_url <- "https://www.onemap.gov.sg/api/common/elastic/search"
query <- list(
searchVal = address,
returnGeom = "Y",
getAddrDetails = "N",
pageNum = "1"
)
res <- GET(base_url, query = query)
restext <- content(res, as = "text", encoding = "UTF-8")
json <- fromJSON(restext, simplifyVector = FALSE)
# Handle missing results safely
if (is.null(json$results) || length(json$results) == 0) {
return(data.frame(LATITUDE = NA, LONGITUDE = NA))
}
result <- json$results[[1]]
# extract LATITUDE and LONGITUDE w/o warning messages
lat <- suppressWarnings(as.numeric(result$LATITUDE))
lon <- suppressWarnings(as.numeric(result$LONGITUDE))
return(data.frame(LATITUDE = lat, LONGITUDE = lon))
}primary_school$LATITUDE <- 0
primary_school$LONGITUDE <- 0
for (i in seq_len(nrow(primary_school))) {
if (primary_school$LATITUDE[i] == 0 ||
primary_school$LONGTITUDE[i] == 0) {
temp_output <- geocode(primary_school[i, 3, drop = TRUE])
primary_school$LATITUDE[i] <- temp_output$LATITUDE
primary_school$LONGITUDE[i] <- temp_output$LONGITUDE
}
}primary_school <- primary_school%>%
st_as_sf(coords = c("LONGITUDE", "LATITUDE"),crs = 4326)%>%
st_transform(3414) Warning in min(cc[[1]], na.rm = TRUE): no non-missing arguments to min;
returning InfWarning in min(cc[[2]], na.rm = TRUE): no non-missing arguments to min;
returning InfWarning in max(cc[[1]], na.rm = TRUE): no non-missing arguments to max;
returning -InfWarning in max(cc[[2]], na.rm = TRUE): no non-missing arguments to max;
returning -Infwrite_rds(primary_school,"data/primary_school.rds")primary_school <- read_rds("data/primary_school.rds")5.3.4 Proximity to Other facilities
Beyond major amenities, proximity to daily necessities and community facilities—including parks, supermarkets, hawker centres, kindergartens, childcare centres, and clinics—significantly influences housing preferences and resale values through enhanced liveability and convenience. To systematically quantify these locational advantages while ensuring computational efficiency and methodological consistency, we implement a standardized spatial processing pipeline. This approach calculates the nearest facility distance for each amenity type through a reproducible four-step workflow:
Data Standardization: Import each facility dataset and transform its coordinate reference system to SVY21 (EPSG:3414) to ensure geometric alignment with our HDB point data.
Distance Computation: Calculate the complete distance matrix between all HDB units and facility locations using st_distance(), which computes Euclidean distances in meters.
Proximity Extraction: For each HDB unit, identify the closest facility by extracting the minimum value from its corresponding row in the distance matrix.
Variable Creation: Round the resulting distances to one decimal place and append them as new predictor variables to the main analytical dataset, creating features like dist_park and dist_supermarket.sf`.
add_nearest_distance <- function(hdb_sf, facility_sf, new_col_name) {
facility_sf <- st_transform(facility_sf, 3414)
dist_matrix <- st_distance(hdb_sf, facility_sf)
nearest_dist <- apply(dist_matrix, 1, min) %>%
as.numeric() %>%
round(1)
hdb_sf <- hdb_sf %>%
mutate(!!new_col_name := nearest_dist) #tidy evaluation
return(hdb_sf)
}
add_nearby_count <- function(hdb_sf, facility_sf,new_col_name,radius = 350) {
facility_sf <- st_transform(facility_sf, 3414)
nearby_list <- st_is_within_distance(hdb_sf, facility_sf, dist = radius)
nearby_count <- lengths(nearby_list)
hdb_sf <- hdb_sf %>%
mutate(!!new_col_name := nearby_count)
return(hdb_sf)
}# BusStop
BusStop <- st_read(dsn = "data/rawdata/BusStop_Aug2025",
layer = "BusStop") %>%st_transform(3414)Reading layer `BusStop' from data source
`D:\ssinha8752\ISSS608-VAA\Take-home_Ex\Take-home_Ex03\data\rawdata\BusStop_Aug2025'
using driver `ESRI Shapefile'
Simple feature collection with 5172 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 3970.122 ymin: 26482.1 xmax: 48285.52 ymax: 52983.82
Projected CRS: SVY21HDBresale_sf <- add_nearest_distance(HDBresale_sf, BusStop, "dist_busStop")
# Mall
mall <- read_csv("data/rawdata/singapore_malls_pois.csv")%>%
st_as_sf(coords = c("lon", "lat"),crs = 4326)%>%
st_transform(3414)Rows: 249 Columns: 8
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (5): name, category, address, website, phone
dbl (2): lat, lon
lgl (1): brand
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.HDBresale_sf <- add_nearest_distance(HDBresale_sf, mall, "dist_mall")
# Hawker
hawker <- st_read("data/rawdata/HawkerCentresGEOJSON.geojson") %>% st_transform(3414)Reading layer `SSOT_HAWKERCENTRES' from data source
`D:\ssinha8752\ISSS608-VAA\Take-home_Ex\Take-home_Ex03\data\rawdata\HawkerCentresGEOJSON.geojson'
using driver `GeoJSON'
Simple feature collection with 129 features and 20 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.6972 ymin: 1.272761 xmax: 103.9882 ymax: 1.448263
Geodetic CRS: WGS 84HDBresale_sf <- add_nearest_distance(HDBresale_sf, hawker, "dist_hawker")
# Supermarket
supermarket <- st_read("data/rawdata/SupermarketsGEOJSON.geojson") %>% st_transform(3414)Reading layer `SupermarketsGEOJSON' from data source
`D:\ssinha8752\ISSS608-VAA\Take-home_Ex\Take-home_Ex03\data\rawdata\SupermarketsGEOJSON.geojson'
using driver `GeoJSON'
Simple feature collection with 526 features and 2 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 103.6258 ymin: 1.24715 xmax: 104.0036 ymax: 1.461526
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84HDBresale_sf <- add_nearest_distance(HDBresale_sf, supermarket, "dist_supermarket")
# Park
park <- st_read("data/rawdata/Parks.geojson") %>% st_transform(3414)Reading layer `NATIONALPARKS' from data source
`D:\ssinha8752\ISSS608-VAA\Take-home_Ex\Take-home_Ex03\data\rawdata\Parks.geojson'
using driver `GeoJSON'
Simple feature collection with 450 features and 6 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.6929 ymin: 1.263958 xmax: 104.0538 ymax: 1.461944
Geodetic CRS: WGS 84HDBresale_sf <- add_nearest_distance(HDBresale_sf, park, "dist_park")
# Kindergartens
kindergartens <- st_read("data/rawdata/Kindergartens.geojson") %>% st_transform(3414)Reading layer `KINDERGARTENS' from data source
`D:\ssinha8752\ISSS608-VAA\Take-home_Ex\Take-home_Ex03\data\rawdata\Kindergartens.geojson'
using driver `GeoJSON'
Simple feature collection with 448 features and 16 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.6887 ymin: 1.247759 xmax: 103.9717 ymax: 1.455452
Geodetic CRS: WGS 84HDBresale_sf <- add_nearest_distance(HDBresale_sf, kindergartens, "dist_kindergarten")
# ChildCareServices
childcare <- st_read("data/rawdata/ChildCareServices.geojson") %>% st_transform(3414)Reading layer `CHILDCARE' from data source
`D:\ssinha8752\ISSS608-VAA\Take-home_Ex\Take-home_Ex03\data\rawdata\ChildCareServices.geojson'
using driver `GeoJSON'
Simple feature collection with 1925 features and 16 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
Geodetic CRS: WGS 84HDBresale_sf <- add_nearest_distance(HDBresale_sf, childcare, "dist_childcare")
# Preschool
preschool <- st_read("data/rawdata/PreSchoolsLocation.kml") %>% st_transform(3414)Reading layer `PRESCHOOLS_LOCATION' from data source
`D:\ssinha8752\ISSS608-VAA\Take-home_Ex\Take-home_Ex03\data\rawdata\PreSchoolsLocation.kml'
using driver `KML'
Simple feature collection with 2290 features and 2 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84HDBresale_sf <- add_nearest_distance(HDBresale_sf, preschool, "dist_preschool")
# Clinics
clinics <- st_read("data/rawdata/CHASClinics.geojson") %>% st_transform(3414)Reading layer `CHASClinics' from data source
`D:\ssinha8752\ISSS608-VAA\Take-home_Ex\Take-home_Ex03\data\rawdata\CHASClinics.geojson'
using driver `GeoJSON'
Simple feature collection with 1193 features and 2 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 103.5818 ymin: 1.016264 xmax: 103.9903 ymax: 1.456037
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84HDBresale_sf <- add_nearest_distance(HDBresale_sf, clinics, "dist_clinics")
# Eldercare
eldercare <- st_read(dsn = "data/rawdata/EldercareServicesSHP",
layer = "ELDERCARE") %>% st_transform(3414)Reading layer `ELDERCARE' from data source
`D:\ssinha8752\ISSS608-VAA\Take-home_Ex\Take-home_Ex03\data\rawdata\EldercareServicesSHP'
using driver `ESRI Shapefile'
Simple feature collection with 133 features and 18 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 14481.92 ymin: 28218.43 xmax: 41665.14 ymax: 46804.9
Projected CRS: SVY21HDBresale_sf <- add_nearest_distance(HDBresale_sf, eldercare, "dist_elderlycare")HDBresale_sf <- add_nearby_count(HDBresale_sf, park,"park_600m",600)
HDBresale_sf <- add_nearby_count(HDBresale_sf, preschool, "preschool_350m",350)
HDBresale_sf <- add_nearby_count(HDBresale_sf, BusStop, "busStop_350m",350)
HDBresale_sf <- add_nearby_count(HDBresale_sf, childcare, "childcare_350m",350)
HDBresale_sf <- add_nearby_count(HDBresale_sf, kindergartens, "kindergartens_350m",350)
HDBresale_sf <- add_nearby_count(HDBresale_sf, primary_school, "top_primary_sch_1km",1000)
HDBresale_sf <- add_nearby_count(HDBresale_sf, supermarket, "supermarket_500m",500)write_rds(HDBresale_sf,"data/HDBresale_sf.rds") HDBresale_sf <- read_rds("data/HDBresale_sf.rds")6 Exploratory Data Analysis
6.1 Visualize the distribution for derived variables
df <- HDBresale_sf[, c(7, 10, 11, 13, 15,16,18:20)] %>% st_drop_geometry()
df$floor_level <- as.numeric(df$floor_level)
HDBresale_sf_num <- df %>%
pivot_longer(cols = everything(),
names_to = "variable",
values_to = "value")
ggplot(HDBresale_sf_num, aes(x = value)) +
geom_histogram(bins = 30, fill = "#93A4CC", color = "grey30") +
facet_wrap(~ variable, scales = "free") +
labs(title = "Histograms of All Numeric Variables",
x = "", y = "Frequency") +
scale_x_continuous(labels = comma) +
theme_minimal()+
theme(plot.title = element_text(size = 14, hjust = 0.5),
panel.background = element_rect(fill = "#f6f6f6"),
plot.background = element_rect(fill = "#f6f6f6",color = NA),
panel.grid = element_blank()
)
Observation
This observed right-skewness carries important implications for our modeling approach. First, it suggests that logarithmic transformations may be warranted for certain variables to normalize their distributions and improve model stability. Second, the skewness reflects the inherent urban structure of Singapore, where most residential developments are strategically planned with reasonable access to basic amenities, creating a natural clustering effect. Third, from a hedonic pricing perspective, this pattern implies diminishing marginal utility for proximity—the price premium for reducing distance from 2km to 1km is likely much larger than from 5km to 4km. Understanding these distributional characteristics is crucial for both interpreting our GWR results and for considering potential variable transformations that could enhance model performance and the validity of our spatial inferences.
6.2 Average Resale Price per sqm vs Total Resale Units by Subzone
To analyze spatial patterns in the HDB resale market, we first aggregate transaction data at the subzone level, calculating both transaction volume and average resale price for each planning area. This spatial summary only includes subzones with active resale activity during our study period, ensuring our analysis reflects actual market behavior rather than imputed values. Building on this foundation, we apply the Jenks natural breaks classification method to identify significant price thresholds within the market, creating five distinct price tiers that maximize homogeneity within groups and heterogeneity between them. The highest price tier, defined by a threshold above SGD 1,028,435, captures the most premium subzones where market values substantially exceed the city-wide average, highlighting areas of exceptional housing demand and locational advantage.
HDBresale_sf <- HDBresale_sf %>%
mutate(price_per_sqm = resale_price / floor_area_sqm)
price_VS_units_sold <- st_join(mpsz19, HDBresale_sf) %>%
group_by(PLN_AREA_N,SUBZONE_N) %>%
summarise(
count = sum(!is.na(resale_price)),
avg_price = mean(resale_price, na.rm = TRUE),
avg_price_sqm = mean(price_per_sqm, na.rm = TRUE),
avg_remaining_lease = mean(remaining_lease, na.rm = TRUE),
.groups = "drop"
)tmap_mode("view")ℹ tmap modes "plot" - "view"
ℹ toggle with `tmap::ttm()`tm_shape(price_VS_units_sold) +
tm_polygons(alpha=0.5,lwd=0.2) +
tm_bubbles(
size = "count",
col = "avg_price_sqm",
col.legend = tm_legend("Avg Price"),
palette = "-hcl.heat",
alpha = 1,
scale = 1.5,
border.col = "#f6f6f6",
title.size = tm_legend(title = "Resale Units")) +
tm_layout(bg.color = "#f6f6f6",
outer.bg.color = "#f6f6f6",
frame = FALSE) +
tm_title("Average HDB Resale Price vs Resale Units by Subzone",
position = tm_pos_out("center", "top", "center"))+
tm_view(set_zoom_limits = c(11,14))
── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.[v3->v4] `tm_bubbles()`: use `fill_alpha` instead of `alpha`.[v3->v4] `tm_tm_bubbles()`: migrate the argument(s) related to the scale of the
visual variable `size` namely 'scale' (rename to 'values.scale') to size.scale
= tm_scale_continuous(<HERE>).
ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
a list: 'size.scale = list(<scale1>, <scale2>, ...)'tmap_mode("plot")ℹ tmap modes "plot" - "view"
Observation
The visualization reveals a distinct spatial gradient in housing values, characterized by a clear distance-decay pattern radiating outward from the Central Business District (CBD). Properties in the central core command significant price premiums per square meter, with values systematically decreasing as distance from the urban center increases, reflecting classic urban economic theories of accessibility and land value. However, this price distribution exhibits an inverse relationship with transaction volume—while the central areas demonstrate the highest price points, they concurrently show lower transaction frequencies compared to the suburban regions. This pattern suggests a constrained supply and potentially different market dynamics in the city center, whereas the suburban areas, despite their lower per-square-meter prices, demonstrate substantially higher market activity and liquidity.
6.3 Top 10 Subzones by Average Resale Price & Factors
To investigate the complex interplay between locational context and housing values, we employ comparative distribution analysis across Singapore’s highest-priced subzones. By systematically mapping different housing attributes—including flat model types, floor levels, and regional classifications—against the fill aesthetic in our visualizations, we can visually decode how various factors contribute to price stratification within premium markets. The ggdist package’s stat_halfeye() plot provides an ideal visualization tool for this purpose, as it simultaneously displays the density distribution, quartile ranges, and median values within a single compact geometry. This multi-layered approach allows us to identify not just central tendencies but also the variability and shape of price distributions for each factor, revealing whether certain attributes create consistent price premiums or exhibit context-dependent effects across different geographical contexts.
Explanation of Enha
top10_subzones <- price_VS_units_sold %>%
arrange(desc(avg_price_sqm)) %>%
slice_head(n = 10)
HDBresale_sf <- st_join(HDBresale_sf, mpsz19[, c("SUBZONE_N", "REGION_N")], left = TRUE)
top10_resale <- HDBresale_sf %>%
filter(SUBZONE_N %in% top10_subzones$SUBZONE_N)ggplot(top10_resale, aes(
x = price_per_sqm,
y = reorder(SUBZONE_N, price_per_sqm, FUN = median)
)) +
stat_halfeye(adjust = 1,
fill="#93A4CC",
justification = -0.2,
.width = 0,
point_colour = NA)+
geom_boxplot(width = 0.2, fill = "white", outlier.size = 1) +
labs(
title = "Top 10 Subzones by Average Resale Price Per sqm",
x = "Resale Price Per sqm (SGD)",
y = "Subzone"
) +
scale_x_continuous(labels = comma) +
theme_ridges() +
theme(
plot.title = element_text(size = 12, hjust = 0.5,
margin = ggplot2::margin(b = 15)),
panel.background = element_rect(fill = "#f6f6f6", color = NA),
plot.background = element_rect(fill = "#f6f6f6", color = NA),
axis.line.x = element_line(color = "grey30"),
axis.ticks = element_line(color = "grey30"),
axis.title.y = element_text(size = 12),
axis.text.y = element_text(hjust = 0.5, size = 8),
axis.text.x = element_text(hjust = 0.5, size = 8,
margin = ggplot2::margin(t = 10)),
axis.title.x = element_text(size = 10,
margin = ggplot2::margin(t = 10)),
legend.position = "right",
legend.title = element_text(size = 8),
legend.key.size = unit(0.3, "cm"),
legend.text = element_text(size = 6),
legend.background = element_rect(fill = "#f6f6f6")
)
ggplot(top10_resale, aes(
x = price_per_sqm,
y = reorder(SUBZONE_N, price_per_sqm, FUN = median),
fill = flat_model,
)) +
stat_halfeye(adjust = 2,
alpha = 0.8,
justification = -0.2,
.width = 0,
point_colour = NA)+
geom_boxplot(width = 0.1, fill = "white", outlier.size = 1) +
labs(
title = "Top 10 Subzones by Average Resale Price Per sqm & Flat Model",
x = "Resale Price Per sqm (SGD)",
y = "Subzone"
) +
scale_x_continuous(labels = comma) +
theme_ridges() +
theme(
plot.title = element_text(size = 14, hjust = 0.5,
margin = ggplot2::margin(b = 15)),
panel.background = element_rect(fill = "#f6f6f6", color = NA),
plot.background = element_rect(fill = "#f6f6f6", color = NA),
axis.line.x = element_line(color = "grey30"),
axis.ticks = element_line(color = "grey30"),
axis.title.y = element_text(size = 12),
axis.text.y = element_text(hjust = 0.5, size = 8),
axis.text.x = element_text(hjust = 0.5, size = 8,
margin = ggplot2::margin(t = 10)),
axis.title.x = element_text(size = 10,
margin = ggplot2::margin(t = 10)),
legend.position = "right",
legend.title = element_text(size = 8),
legend.key.size = unit(0.3, "cm"),
legend.text = element_text(size = 8),
legend.background = element_rect(fill = "#f6f6f6")
)
ggplot(top10_resale, aes(
x = price_per_sqm,
y = reorder(SUBZONE_N, price_per_sqm, FUN = median),
fill = floor_level,
)) +
stat_halfeye(adjust = 2,
alpha = 0.8,
justification = -0.2,
.width = 0,
point_colour = NA)+
geom_boxplot(width = 0.1, fill = "white", outlier.size = 1) +
scale_x_continuous(labels = comma) +
labs(
title = "Top 10 Subzones by Average Resale Price Per sqm & Floor Level",
x = "Resale Price Per sqm (SGD)",
y = "Subzone"
) +
theme_ridges() +
theme(
plot.title = element_text(size = 14, hjust = 0.5,
margin = ggplot2::margin(b = 15)),
panel.background = element_rect(fill = "#f6f6f6", color = NA),
plot.background = element_rect(fill = "#f6f6f6", color = NA),
axis.line.x = element_line(color = "grey30"),
axis.ticks = element_line(color = "grey30"),
axis.title.y = element_text(size = 12),
axis.text.y = element_text(hjust = 0.5, size = 8),
axis.text.x = element_text(hjust = 0.5, size = 8,
margin = ggplot2::margin(t = 10)),
axis.title.x = element_text(size = 10,
margin = ggplot2::margin(t = 10)),
legend.position = "right",
legend.title = element_text(size = 8),
legend.key.size = unit(0.3, "cm"),
legend.text = element_text(size = 8),
legend.background = element_rect(fill = "#f6f6f6")
)
7 Number of Top 50 Primary School within 1 KM
ggplot(top10_resale, aes(
x = price_per_sqm,
y = reorder(SUBZONE_N, price_per_sqm, FUN = median),
fill = preschool_350m,
)) +
stat_halfeye(adjust = 2,
alpha = 0.8,
justification = -0.2,
.width = 0,
point_colour = NA)+
geom_boxplot(width = 0.1, fill = "white", outlier.size = 1) +
scale_x_continuous(labels = comma) +
labs(
title = "Top 10 Subzones by Average Resale Price Per sqm \n & Number of Top Primary School within 1 KM ",
x = "Resale Price Per sqm (SGD)",
y = "Subzone"
) +
theme_ridges() +
theme(
plot.title = element_text(size = 14, hjust = 0.5,
margin = ggplot2::margin(b = 15)),
panel.background = element_rect(fill = "#f6f6f6", color = NA),
plot.background = element_rect(fill = "#f6f6f6", color = NA),
axis.line.x = element_line(color = "grey30"),
axis.ticks = element_line(color = "grey30"),
axis.title.y = element_text(size = 12),
axis.text.y = element_text(hjust = 0.5, size = 8),
axis.text.x = element_text(hjust = 0.5, size = 8,
margin = ggplot2::margin(t = 10)),
axis.title.x = element_text(size = 10,
margin = ggplot2::margin(t = 10)),
legend.position = "right",
legend.title = element_text(size = 8),
legend.key.size = unit(0.3, "cm"),
legend.text = element_text(size = 8),
legend.background = element_rect(fill = "#f6f6f6")
)
7.1 Scatter Plot between Resale Prices and Independent Variables
df <- df %>%
mutate(price_per_sqm = resale_price / floor_area_sqm)
vars_to_compare <- setdiff(names(df), c("resale_price","price_per_sqm"))
df_plot <- df %>%
mutate(across(all_of(vars_to_compare))) %>%
pivot_longer(
cols = all_of(vars_to_compare),
names_to = "variable",
values_to = "value"
)
ggplot(df_plot, aes(x = value, y = price_per_sqm)) +
geom_point(alpha = 0.3, color="#93A4CC") +
geom_smooth(method = "lm", color = "black") +
facet_wrap(~ variable, scales = "free_x", ncol = 4) +
labs(
title = "Scatter Plots of Derived Variables vs Resale Price",
x = "Value",
y = "Resale Price (SGD)"
) +
theme_minimal(base_size = 10) +
theme(
plot.title = element_text(size = 14, hjust = 0.5),
plot.background = element_rect(fill = "#f6f6f6", color = NA),
panel.background = element_rect(fill = "#f6f6f6", color = NA),
strip.background = element_rect(fill = "#e1e1e1"),
axis.text.x = element_text(angle = 0, hjust = 1, size = 8),
panel.grid = element_blank()
)`geom_smooth()` using formula = 'y ~ x'
7.1.1 Remaining Lease
ggplot(HDBresale_sf %>% st_drop_geometry(),
aes(x = remaining_lease,
y = resale_price,
color = REGION_N)) +
geom_point(alpha = 0.6, size = 2) +
geom_text_repel(
aes(label = SUBZONE_N),
size = 3,
box.padding = 0.3,
point.padding = 0.5,
max.overlaps = 8
) +
scale_color_brewer(palette = "Set2") +
labs(
title = "Remaining Lease vs Resale Price",
x = "Remaining Lease (years)",
y = "HDB Resale Price (SGD)",
color = "Planning Region"
) +
theme_minimal()+
theme(plot.title = element_text(size = 14, hjust = 0.5),
panel.background = element_rect(fill = "#f6f6f6"),
plot.background = element_rect(fill = "#f6f6f6",color = NA),
panel.grid = element_blank()
)Warning: ggrepel: 8657 unlabeled data points (too many overlaps). Consider
increasing max.overlaps
7.2 Treemap by Planning Region, Subzone and Flatmodel
HDBresale_sf <- HDBresale_sf %>%
mutate(model_group = case_when(
flat_model %in% c("Standard", "Simplified", "Improved", "Model A2", "New Generation") ~ 1,
flat_model %in% c("Model A", "Premium Apartment", "Adjoined flat") ~ 2,
flat_model %in% c("DBSS", "Terrace") ~ 3,
flat_model %in% c("Premium Apartment Loft", "Type S1") ~ 4,
TRUE ~ NA
))flat_summary <- HDBresale_sf %>%
st_drop_geometry() %>%
group_by(flat_model) %>%
summarise(
median_price = median(resale_price, na.rm = TRUE),
units_sold = n(),
.groups = "drop"
) %>%
arrange(median_price) %>%
mutate(flat_model = factor(flat_model, levels = flat_model)) # order factor
ggplot(HDBresale_sf %>% st_drop_geometry(),
aes(x = factor(flat_model, levels = flat_summary$flat_model),
y = resale_price)) +
geom_boxplot(alpha = 0.7, outlier.size = 1, fill="#93A4CC") +
# units sold
geom_point(data = flat_summary,
aes(x = flat_model, y = median_price, size = units_sold),
color = "black", alpha = 0.6) +
scale_size_continuous(name = "Units Sold") +
scale_y_continuous(labels = comma) +
labs(
title = "HDB Resale Price by Flat Model (Ordered by Median Price)",
x = "Flat Model",
y = "Resale Price (SGD)",
fill = "Flat Model"
) +
theme_classic()+
theme(
plot.title = element_text(size = 12, hjust = 0.5),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.title = element_text(size = 8),
legend.key.size = unit(0.3, "cm"),
legend.text = element_text(size = 8),
panel.background = element_rect(fill = "#f6f6f6", color = NA),
plot.background = element_rect(fill = "#f6f6f6", color = NA),
legend.background = element_rect(fill = "#f6f6f6")
)Ignoring unknown labels:
• fill : "Flat Model"
treemap_data <- HDBresale_sf %>%
st_drop_geometry() %>% # remove geometry for speed
group_by(REGION_N,SUBZONE_N,flat_model) %>%
summarise(
count = n(),
med_resale_price= median(resale_price, na.rm = TRUE),
.groups = "drop"
)par(bg = "#f6f6f6")
treemap(treemap_data,
index=c("flat_model","REGION_N","SUBZONE_N"),
vSize="count",
vColor="med_resale_price",
sortID = "REGION_N",
type="value",
title="4-Room HDB Resale Transactions by Planning Region, Flat Model & Subzone: Jul 2024 – Jun 2025",
title.legend = "Median Resale Price"
)Warning in treemap(treemap_data, index = c("flat_model", "REGION_N",
"SUBZONE_N"), : sortID must be a numeric variable
Observation
From both box-plot and tree map, We can see Model A is the most common model across all the region. Based on the treemap, we can observe that although overall resale price is mainly affected by location, different model like Type S1 vs standard also has huge difference, hence, the model type should also considered as variable for the predictive model.
The model will be split into 4 groups based on their box plot distribution:
- Low (1): Standard, Simplified, Improved, Model A2, New Genegration
- Mid (2): Model A, Premium Apartment, Adjoined Flat
- Mid-High (3): DBSS, Terrace
- High (4): Premium Apartment Loft, Type S1
HDBresale_sf <- HDBresale_sf %>%
mutate(model_group = case_when(
flat_model %in% c("Standard", "Simplified", "Improved", "Model A2", "New Generation") ~ 1,
flat_model %in% c("Model A", "Premium Apartment", "Adjoined flat") ~ 2,
flat_model %in% c("DBSS", "Terrace") ~ 3,
flat_model %in% c("Premium Apartment Loft", "Type S1") ~ 4,
TRUE ~ NA
))8 Calibrate Explanatory Model
8.1 Computing Correlation Matrix
mdata <- HDBresale_sf %>%
arrange(resale_date)%>%
select(price_per_sqm, everything(),
-c(month,town,flat_type,block,street_name,storey_max,storey_range,flat_model,
lease_commence_date,SUBZONE_N,REGION_N,resale_price,resale_date))par(bg="#f6f6f6")
mdata_nogeo <- mdata %>%
st_drop_geometry() %>%
mutate(across(where(is.factor), ~ as.numeric(.)))
corrplot::corrplot(cor(mdata_nogeo[, 2:23]),
diag = FALSE,
tl.pos = "td",
tl.cex = 1,
method = "number",
type = "upper")
From the correlation matrix above, we observe that three pairs of independent variables are highly correlated:
Remaining lease vs. unit_age (-1): These variables capture the same concept and are interchangeable. Therefore,
unit_agewill be dropped.dist_childcare vs. dist_preschool (0.99): Since preschool has more data points,
dist_childcarewill be removed.preschool_350m vs. childcare_350m (0.97): For the same reason as above,
childcare_350mwill be dropped.mdata <- mdata %>% select(-c(unit_age,dist_childcare,childcare_350m))
par(bg="#f6f6f6")
mdata_nogeo <- mdata %>%
st_drop_geometry() %>%
mutate(across(where(is.factor), ~ as.numeric(.)))
corrplot::corrplot(cor(mdata_nogeo[, 2:23]),
diag = FALSE,
tl.pos = "td",
tl.cex = 1,
method = "number",
type = "upper")Warning in cor(mdata_nogeo[, 2:23]): the standard deviation is zero
mdata <- mdata %>%
mutate(
model_group = as.factor(model_group),
floor_level = as.factor(floor_level) # your 5-level bin
)8.2 Building Hedonic Pricing Models using GWmodel
8.2.1 Compute fixed bandwith
bw_fixed_CV <- bw.gwr(
formula = price_per_sqm ~ floor_area_sqm + remaining_lease + floor_level +
dist_to_CBD + dist_to_MRT + dist_busStop +
dist_mall + dist_hawker + dist_supermarket + dist_park +
dist_kindergarten + dist_preschool + dist_clinics + dist_elderlycare +
MRT_600m + park_600m + preschool_350m + busStop_350m + kindergartens_350m +
supermarket_500m + model_group,
data = mdata,
approach = "CV",
kernel = "gaussian",
adaptive = FALSE,
longlat = FALSE
)Take a cup of tea and have a break, it will take a few minutes.
-----A kind suggestion from GWmodel development group
Fixed bandwidth: 22980.96 CV score: 4655095701
Fixed bandwidth: 14205.85 CV score: 4497948661
Fixed bandwidth: 8782.54 CV score: 4105167018
Fixed bandwidth: 5430.749 CV score: 3306972173
Fixed bandwidth: 3359.228 CV score: 2311279853
Fixed bandwidth: 2078.957 CV score: 1577651920
Fixed bandwidth: 1287.706 CV score: 1246997073
Fixed bandwidth: 798.6866 CV score: 1058482760
Fixed bandwidth: 496.4557 CV score: 1168585756
Fixed bandwidth: 985.4755 CV score: 1128770780
Fixed bandwidth: 683.2446 CV score: 1021057895
Fixed bandwidth: 611.8976 CV score: 1063255449
Fixed bandwidth: 727.3395 CV score: 1033030202
Fixed bandwidth: 655.9925 CV score: 1022702738
Fixed bandwidth: 700.0874 CV score: 1024604293
Fixed bandwidth: 672.8352 CV score: 1020223906
Fixed bandwidth: 666.4019 CV score: 1020504712
Fixed bandwidth: 676.8113 CV score: 1020377948
Fixed bandwidth: 670.3779 CV score: 1020247931
Fixed bandwidth: 674.354 CV score: 1020255944
Fixed bandwidth: 671.8966 CV score: 1020221690
Fixed bandwidth: 671.3165 CV score: 1020227284
Fixed bandwidth: 672.2551 CV score: 1020220911
Fixed bandwidth: 672.4767 CV score: 1020221438
Fixed bandwidth: 672.1182 CV score: 1020220969
Fixed bandwidth: 672.3398 CV score: 1020221022
Fixed bandwidth: 672.2028 CV score: 1020220898
Fixed bandwidth: 672.1705 CV score: 1020220912
Fixed bandwidth: 672.2228 CV score: 1020220898
Fixed bandwidth: 672.2352 CV score: 1020220901
Fixed bandwidth: 672.2152 CV score: 1020220897
Fixed bandwidth: 672.2105 CV score: 1020220897
Fixed bandwidth: 672.2181 CV score: 1020220897
Fixed bandwidth: 672.2134 CV score: 1020220897
Fixed bandwidth: 672.2163 CV score: 1020220897
Fixed bandwidth: 672.2145 CV score: 1020220897
Fixed bandwidth: 672.2141 CV score: 1020220897
Fixed bandwidth: 672.2138 CV score: 1020220897
Fixed bandwidth: 672.2142 CV score: 1020220897
Fixed bandwidth: 672.2143 CV score: 1020220897
Fixed bandwidth: 672.2142 CV score: 1020220897
Fixed bandwidth: 672.2143 CV score: 1020220897
Fixed bandwidth: 672.2142 CV score: 1020220897 The fixed bandwidth is 672.2142 using CV method
write_rds(bw_fixed_CV,"data/bw_fixed_CV.rds")bw_fixed_CV<- read_rds("data/bw_fixed_CV.rds")
gwr_fixed_CV <- gwr.basic(
formula = price_per_sqm ~ floor_area_sqm + remaining_lease + floor_level +
dist_to_CBD + dist_to_MRT + dist_busStop +
dist_mall + dist_hawker + dist_supermarket + dist_park +
dist_kindergarten + dist_preschool + dist_clinics + dist_elderlycare +
park_600m + preschool_350m + busStop_350m + kindergartens_350m +
supermarket_500m + model_group,
data=mdata,
bw=bw_fixed_CV,
kernel = 'gaussian',
longlat = FALSE)write_rds(gwr_fixed_CV,"data/gwr_fixed_CV.rds")gwr_fixed_CV <- read_rds("data/gwr_fixed_CV.rds")
gwr_fixed_CV ***********************************************************************
* Package GWmodel *
***********************************************************************
Program starts at: 2025-11-17 23:37:44.19304
Call:
gwr.basic(formula = price_per_sqm ~ floor_area_sqm + remaining_lease +
floor_level + dist_to_CBD + dist_to_MRT + dist_busStop +
dist_mall + dist_hawker + dist_supermarket + dist_park +
dist_kindergarten + dist_preschool + dist_clinics + dist_elderlycare +
park_600m + preschool_350m + busStop_350m + kindergartens_350m +
supermarket_500m + model_group, data = mdata, bw = bw_fixed_CV,
kernel = "gaussian", longlat = FALSE)
Dependent (y) variable: price_per_sqm
Independent variables: floor_area_sqm remaining_lease floor_level dist_to_CBD dist_to_MRT dist_busStop dist_mall dist_hawker dist_supermarket dist_park dist_kindergarten dist_preschool dist_clinics dist_elderlycare park_600m preschool_350m busStop_350m kindergartens_350m supermarket_500m model_group
Number of data points: 8679
***********************************************************************
* Results of Global Regression *
***********************************************************************
Call:
lm(formula = formula, data = data)
Residuals:
Min 1Q Median 3Q Max
-3364.2 -490.2 15.1 445.1 4297.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.013e+03 2.278e+02 26.399 < 2e-16 ***
floor_area_sqm -1.605e+01 1.899e+00 -8.449 < 2e-16 ***
remaining_lease 6.499e+01 1.081e+00 60.115 < 2e-16 ***
floor_level2 2.419e+02 2.658e+01 9.102 < 2e-16 ***
floor_level3 4.140e+02 2.674e+01 15.479 < 2e-16 ***
floor_level4 5.455e+02 2.444e+01 22.320 < 2e-16 ***
floor_level5 1.968e+03 4.525e+01 43.491 < 2e-16 ***
dist_to_CBD -2.044e-01 2.416e-03 -84.589 < 2e-16 ***
dist_to_MRT -4.364e-01 2.360e-02 -18.491 < 2e-16 ***
dist_busStop 2.871e-01 1.646e-01 1.744 0.0811 .
dist_mall -2.803e-01 2.672e-02 -10.490 < 2e-16 ***
dist_hawker -2.690e-01 1.918e-02 -14.024 < 2e-16 ***
dist_supermarket 3.360e-01 5.580e-02 6.022 1.79e-09 ***
dist_park 6.534e-02 2.569e-02 2.543 0.0110 *
dist_kindergarten -2.632e-01 5.616e-02 -4.686 2.82e-06 ***
dist_preschool -6.024e-02 1.240e-01 -0.486 0.6270
dist_clinics 3.787e-02 8.198e-02 0.462 0.6441
dist_elderlycare -1.036e-01 1.447e-02 -7.159 8.82e-13 ***
park_600m 1.134e+02 8.987e+00 12.621 < 2e-16 ***
preschool_350m -4.684e+01 4.272e+00 -10.964 < 2e-16 ***
busStop_350m 3.376e+00 3.255e+00 1.037 0.2997
kindergartens_350m 1.296e+02 1.160e+01 11.180 < 2e-16 ***
supermarket_500m 1.008e+02 7.648e+00 13.181 < 2e-16 ***
model_group2 -1.124e+01 3.938e+01 -0.285 0.7753
model_group3 1.494e+03 7.686e+01 19.437 < 2e-16 ***
model_group4 2.679e+03 1.321e+02 20.270 < 2e-16 ***
---Significance stars
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 751.4 on 8653 degrees of freedom
Multiple R-squared: 0.8281
Adjusted R-squared: 0.8276
F-statistic: 1668 on 25 and 8653 DF, p-value: < 2.2e-16
***Extra Diagnostic information
Residual sum of squares: 4884965275
Sigma(hat): 750.319
AIC: 139600.5
AICc: 139600.7
BIC: 131357.3
***********************************************************************
* Results of Geographically Weighted Regression *
***********************************************************************
*********************Model calibration information*********************
Kernel function: gaussian
Fixed bandwidth: 672.2142
Regression points: the same locations as observations are used.
Distance metric: Euclidean distance metric is used.
****************Summary of GWR coefficient estimates:******************
Min. 1st Qu. Median 3rd Qu.
Intercept -3.5843e+03 3.8903e+03 6.4167e+03 8.6772e+03
floor_area_sqm -4.2606e+02 -3.8063e+01 -2.9129e+01 -2.0980e+01
remaining_lease -3.7643e+02 4.6611e+01 5.7164e+01 7.3462e+01
floor_level2 -4.0920e+02 1.5423e+02 2.2683e+02 3.0761e+02
floor_level3 -1.7619e+01 3.2404e+02 4.2516e+02 5.3522e+02
floor_level4 1.9430e+02 4.6377e+02 6.2161e+02 7.3640e+02
floor_level5 -9.3204e+37 4.4364e+02 9.9035e+02 1.4994e+03
dist_to_CBD -1.7805e+00 -2.9932e-01 -1.0267e-01 6.1503e-02
dist_to_MRT -3.4208e+00 -8.7592e-01 -4.9720e-01 -1.6733e-01
dist_busStop -4.2855e+00 -5.9561e-01 -1.2874e-01 6.9800e-01
dist_mall -3.3404e+00 -7.0561e-01 -2.9149e-01 7.7960e-02
dist_hawker -5.7742e+00 -3.1459e-01 6.4994e-03 3.1578e-01
dist_supermarket -7.2036e+00 -5.1646e-01 -1.4130e-01 2.6524e-01
dist_park -7.1020e+00 -2.0538e-01 -1.7320e-04 2.3653e-01
dist_kindergarten -2.8872e+00 -1.8935e-01 1.0024e-01 5.7194e-01
dist_preschool -1.8827e+01 -4.4825e-01 -1.1806e-01 2.4189e-01
dist_clinics -5.9629e+00 -3.1940e-01 1.9226e-02 3.3525e-01
dist_elderlycare -2.8077e+00 -2.2201e-01 -3.3851e-02 1.6265e-01
park_600m -7.3192e+02 -4.9963e+01 2.9536e+00 6.4593e+01
preschool_350m -5.9658e+02 -2.4785e+01 -2.2362e+00 2.8998e+01
busStop_350m -3.9683e+02 -1.3531e+01 4.2461e+00 1.5123e+01
kindergartens_350m -1.2491e+03 -7.3308e+01 -1.3861e+01 5.3213e+01
supermarket_500m -7.6839e+02 -6.5539e+00 2.4385e+01 8.6818e+01
model_group2 -3.0426e+03 -1.3500e+02 7.7710e+00 1.4505e+02
model_group3 -4.4349e+03 9.9500e+02 1.5079e+03 1.8523e+03
model_group4 -7.4898e+35 1.7541e+03 4.1780e+03 6.4699e+03
Max.
Intercept 7.1238e+04
floor_area_sqm 4.0660e+01
remaining_lease 1.6295e+02
floor_level2 8.1150e+02
floor_level3 1.6558e+03
floor_level4 1.7468e+03
floor_level5 9.0775e+36
dist_to_CBD 1.0791e+00
dist_to_MRT 5.8927e+00
dist_busStop 9.3336e+00
dist_mall 1.8475e+00
dist_hawker 4.0307e+00
dist_supermarket 3.8384e+00
dist_park 2.5008e+00
dist_kindergarten 6.6986e+00
dist_preschool 4.5883e+00
dist_clinics 4.3097e+00
dist_elderlycare 2.0773e+00
park_600m 7.7157e+02
preschool_350m 7.6399e+02
busStop_350m 1.6250e+02
kindergartens_350m 8.7900e+02
supermarket_500m 3.2524e+02
model_group2 1.8389e+04
model_group3 2.0389e+04
model_group4 1.7473e+37
************************Diagnostic information*************************
Number of data points: 8679
Effective number of parameters (2trace(S) - trace(S'S)): 1410.186
Effective degrees of freedom (n-2trace(S) + trace(S'S)): 7268.814
AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 125340.7
AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 123872.2
BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 124297.3
Residual sum of squares: 704789520
R-square value: 0.9752025
Adjusted R-square value: 0.9703909
***********************************************************************
Program stops at: 2025-11-17 23:38:58.482437 As we removed dist_busStop, dist_park, dist_preschool, dist_clinics, top_primary_sch_1km and then run again
bw_fixed_CV <- bw.gwr(
formula = price_per_sqm ~ floor_area_sqm + remaining_lease + floor_level +
dist_to_CBD + dist_to_MRT +
dist_mall + dist_hawker + dist_supermarket +
dist_kindergarten + dist_elderlycare +
park_600m + preschool_350m + busStop_350m + kindergartens_350m +
supermarket_500m + model_group,
data = mdata,
approach = "CV",
kernel = "gaussian",
adaptive = FALSE,
longlat = FALSE
)Take a cup of tea and have a break, it will take a few minutes.
-----A kind suggestion from GWmodel development group
Fixed bandwidth: 22980.96 CV score: 4810016253
Fixed bandwidth: 14205.85 CV score: 4639942412
Fixed bandwidth: 8782.54 CV score: 4220229509
Fixed bandwidth: 5430.749 CV score: 3379881814
Fixed bandwidth: 3359.228 CV score: 2359720175
Fixed bandwidth: 2078.957 CV score: 1635393930
Fixed bandwidth: 1287.706 CV score: 1294539880
Fixed bandwidth: 798.6866 CV score: 1095119153
Fixed bandwidth: 496.4557 CV score: 1722347028
Fixed bandwidth: 985.4755 CV score: 1171279710
Fixed bandwidth: 683.2446 CV score: 1054155227
Fixed bandwidth: 611.8976 CV score: 1142472270
Fixed bandwidth: 727.3395 CV score: 1066292401
Fixed bandwidth: 655.9925 CV score: 1060552390
Fixed bandwidth: 700.0874 CV score: 1057247174
Fixed bandwidth: 672.8352 CV score: 1054299229
Fixed bandwidth: 689.678 CV score: 1054962966
Fixed bandwidth: 679.2686 CV score: 1053966817
Fixed bandwidth: 676.8113 CV score: 1053992475
Fixed bandwidth: 680.7873 CV score: 1054007079
Fixed bandwidth: 678.33 CV score: 1053962916
Fixed bandwidth: 677.7499 CV score: 1053968864
Fixed bandwidth: 678.6885 CV score: 1053962459
Fixed bandwidth: 678.9101 CV score: 1053963385
Fixed bandwidth: 678.5516 CV score: 1053962347
Fixed bandwidth: 678.4669 CV score: 1053962454
Fixed bandwidth: 678.6039 CV score: 1053962348
Fixed bandwidth: 678.5192 CV score: 1053962372
Fixed bandwidth: 678.5715 CV score: 1053962341
Fixed bandwidth: 678.5839 CV score: 1053962341
Fixed bandwidth: 678.5639 CV score: 1053962342
Fixed bandwidth: 678.5763 CV score: 1053962341
Fixed bandwidth: 678.5792 CV score: 1053962341
Fixed bandwidth: 678.5745 CV score: 1053962341
Fixed bandwidth: 678.5774 CV score: 1053962341
Fixed bandwidth: 678.5756 CV score: 1053962341
Fixed bandwidth: 678.5767 CV score: 1053962341
Fixed bandwidth: 678.5769 CV score: 1053962341
Fixed bandwidth: 678.5765 CV score: 1053962341
Fixed bandwidth: 678.5768 CV score: 1053962341
Fixed bandwidth: 678.5766 CV score: 1053962341
Fixed bandwidth: 678.5767 CV score: 1053962341
Fixed bandwidth: 678.5767 CV score: 1053962341 write_rds(bw_fixed_CV,"data_bw_fixed_CV2.rds")gwr_fixed_CV <- gwr.basic(
formula = price_per_sqm ~ floor_area_sqm + remaining_lease + floor_level +
dist_to_CBD + dist_to_MRT +
dist_mall + dist_hawker + dist_supermarket +
dist_kindergarten + dist_elderlycare +
park_600m + preschool_350m + busStop_350m + kindergartens_350m +
supermarket_500m + model_group,
data=mdata,
bw=bw_fixed_CV,
kernel = 'gaussian',
longlat = FALSE)
write_rds(bw_fixed_CV,"bw_fixed_CV.rds")gwr_fixed_CV ***********************************************************************
* Package GWmodel *
***********************************************************************
Program starts at: 2025-11-18 00:01:35.411795
Call:
gwr.basic(formula = price_per_sqm ~ floor_area_sqm + remaining_lease +
floor_level + dist_to_CBD + dist_to_MRT + dist_mall + dist_hawker +
dist_supermarket + dist_kindergarten + dist_elderlycare +
park_600m + preschool_350m + busStop_350m + kindergartens_350m +
supermarket_500m + model_group, data = mdata, bw = bw_fixed_CV,
kernel = "gaussian", longlat = FALSE)
Dependent (y) variable: price_per_sqm
Independent variables: floor_area_sqm remaining_lease floor_level dist_to_CBD dist_to_MRT dist_mall dist_hawker dist_supermarket dist_kindergarten dist_elderlycare park_600m preschool_350m busStop_350m kindergartens_350m supermarket_500m model_group
Number of data points: 8679
***********************************************************************
* Results of Global Regression *
***********************************************************************
Call:
lm(formula = formula, data = data)
Residuals:
Min 1Q Median 3Q Max
-3355.9 -486.3 16.1 442.4 4338.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.066e+03 2.239e+02 27.096 < 2e-16 ***
floor_area_sqm -1.579e+01 1.897e+00 -8.325 < 2e-16 ***
remaining_lease 6.509e+01 1.062e+00 61.285 < 2e-16 ***
floor_level2 2.424e+02 2.658e+01 9.120 < 2e-16 ***
floor_level3 4.130e+02 2.674e+01 15.444 < 2e-16 ***
floor_level4 5.453e+02 2.444e+01 22.308 < 2e-16 ***
floor_level5 1.968e+03 4.524e+01 43.492 < 2e-16 ***
dist_to_CBD -2.031e-01 2.356e-03 -86.211 < 2e-16 ***
dist_to_MRT -4.359e-01 2.359e-02 -18.476 < 2e-16 ***
dist_mall -2.775e-01 2.670e-02 -10.396 < 2e-16 ***
dist_hawker -2.689e-01 1.912e-02 -14.059 < 2e-16 ***
dist_supermarket 3.368e-01 5.371e-02 6.271 3.75e-10 ***
dist_kindergarten -2.584e-01 5.561e-02 -4.647 3.41e-06 ***
dist_elderlycare -1.086e-01 1.432e-02 -7.581 3.77e-14 ***
park_600m 1.027e+02 7.852e+00 13.074 < 2e-16 ***
preschool_350m -4.622e+01 4.215e+00 -10.966 < 2e-16 ***
busStop_350m 1.579e+00 3.094e+00 0.510 0.610
kindergartens_350m 1.292e+02 1.156e+01 11.182 < 2e-16 ***
supermarket_500m 1.018e+02 7.631e+00 13.340 < 2e-16 ***
model_group2 -1.132e+01 3.926e+01 -0.288 0.773
model_group3 1.479e+03 7.664e+01 19.293 < 2e-16 ***
model_group4 2.673e+03 1.321e+02 20.238 < 2e-16 ***
---Significance stars
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 751.6 on 8657 degrees of freedom
Multiple R-squared: 0.8279
Adjusted R-squared: 0.8275
F-statistic: 1984 on 21 and 8657 DF, p-value: < 2.2e-16
***Extra Diagnostic information
Residual sum of squares: 4890447019
Sigma(hat): 750.7399
AIC: 139602.3
AICc: 139602.4
BIC: 131294.4
***********************************************************************
* Results of Geographically Weighted Regression *
***********************************************************************
*********************Model calibration information*********************
Kernel function: gaussian
Fixed bandwidth: 678.5767
Regression points: the same locations as observations are used.
Distance metric: Euclidean distance metric is used.
****************Summary of GWR coefficient estimates:******************
Min. 1st Qu. Median 3rd Qu.
Intercept -5.3878e+03 3.9657e+03 6.5644e+03 9.0640e+03
floor_area_sqm -3.4518e+02 -3.8538e+01 -2.9359e+01 -2.0393e+01
remaining_lease -2.2872e+02 4.6336e+01 5.7063e+01 7.1713e+01
floor_level2 -4.8812e+02 1.5566e+02 2.2730e+02 3.0003e+02
floor_level3 -9.8179e+01 3.2202e+02 4.2280e+02 5.3271e+02
floor_level4 2.2354e+02 4.6416e+02 6.0743e+02 7.2737e+02
floor_level5 -6.1019e+36 5.3711e+02 1.0207e+03 1.4587e+03
dist_to_CBD -9.5249e-01 -3.2285e-01 -1.1960e-01 7.4781e-02
dist_to_MRT -3.4236e+00 -8.6092e-01 -5.2542e-01 -2.2625e-01
dist_mall -4.6694e+00 -6.2933e-01 -2.3513e-01 7.0526e-02
dist_hawker -3.7817e+00 -2.8361e-01 1.0859e-03 2.5780e-01
dist_supermarket -3.1855e+00 -5.0994e-01 -1.8297e-01 2.1538e-01
dist_kindergarten -4.3882e+00 -1.9435e-01 3.4679e-02 4.2540e-01
dist_elderlycare -3.7126e+00 -2.7857e-01 -5.5558e-02 1.4540e-01
park_600m -1.2001e+03 -6.6256e+01 5.1819e+00 6.3417e+01
preschool_350m -2.3301e+02 -2.2718e+01 1.0884e+00 2.9436e+01
busStop_350m -2.0409e+02 -1.5170e+01 3.3228e+00 1.3517e+01
kindergartens_350m -7.0783e+02 -7.5257e+01 -1.3503e+01 5.2554e+01
supermarket_500m -5.4379e+02 -1.3116e+01 2.6995e+01 9.0547e+01
model_group2 -2.8500e+03 -1.3684e+02 2.3279e+01 1.6494e+02
model_group3 -2.5272e+03 9.8036e+02 1.5172e+03 1.8389e+03
model_group4 -3.2303e+36 1.5959e+03 3.9684e+03 6.8291e+03
Max.
Intercept 4.3673e+04
floor_area_sqm 4.2616e+01
remaining_lease 1.6069e+02
floor_level2 9.6840e+02
floor_level3 1.6497e+03
floor_level4 1.4380e+03
floor_level5 2.7933e+36
dist_to_CBD 1.2553e+00
dist_to_MRT 1.7633e+00
dist_mall 2.2999e+00
dist_hawker 2.7559e+00
dist_supermarket 4.7624e+00
dist_kindergarten 3.3142e+00
dist_elderlycare 2.2852e+00
park_600m 5.1669e+02
preschool_350m 4.4785e+02
busStop_350m 2.2887e+02
kindergartens_350m 3.3286e+02
supermarket_500m 3.8812e+02
model_group2 1.3471e+04
model_group3 1.5263e+04
model_group4 5.7557e+35
************************Diagnostic information*************************
Number of data points: 8679
Effective number of parameters (2trace(S) - trace(S'S)): 1191.027
Effective degrees of freedom (n-2trace(S) + trace(S'S)): 7487.973
AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 125622.3
AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 124441.8
BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 123397.1
Residual sum of squares: 768562842
R-square value: 0.9729586
Adjusted R-square value: 0.9686569
***********************************************************************
Program stops at: 2025-11-18 00:02:26.644254 bw_fixed_AIC <- bw.gwr(
formula = price_per_sqm ~ floor_area_sqm + remaining_lease + floor_level +
dist_to_CBD + dist_to_MRT + dist_busStop +
dist_mall + dist_hawker + dist_supermarket + dist_kindergarten + dist_elderlycare +
park_600m + preschool_350m + kindergartens_350m + supermarket_500m + model_group,
data = mdata,
approach = "AIC",
kernel = "gaussian",
adaptive = FALSE,
longlat = FALSE
)Take a cup of tea and have a break, it will take a few minutes.
-----A kind suggestion from GWmodel development group
Fixed bandwidth: 22980.96 AICc value: 139414.2
Fixed bandwidth: 14205.85 AICc value: 139104.1
Fixed bandwidth: 8782.54 AICc value: 138279.1
Fixed bandwidth: 5430.749 AICc value: 136342.9
Fixed bandwidth: 3359.228 AICc value: 133192.2
Fixed bandwidth: 2078.957 AICc value: 129860.8
Fixed bandwidth: 1287.706 AICc value: 127682.2
Fixed bandwidth: 798.6866 AICc value: 126050.8
Fixed bandwidth: 496.4557 AICc value: 124967.5
Error in eval(expr, envir) : inv(): matrix is singular
Fixed bandwidth: 309.6667 AICc value: Inf
Fixed bandwidth: 611.8976 AICc value: 125342.4
Error in eval(expr, envir) : inv(): matrix is singular
Fixed bandwidth: 425.1087 AICc value: Inf
Fixed bandwidth: 540.5506 AICc value: 125092.5
Fixed bandwidth: 469.2035 AICc value: Inf
Fixed bandwidth: 513.2984 AICc value: 125011.7
Fixed bandwidth: 486.0463 AICc value: 124942.7
Fixed bandwidth: 479.6129 AICc value: 124928.3
Fixed bandwidth: 475.6369 AICc value: 124919.8
Fixed bandwidth: 473.1796 AICc value: Inf
Fixed bandwidth: 477.1556 AICc value: 124923
Fixed bandwidth: 474.6983 AICc value: 124917.8
Fixed bandwidth: 474.1182 AICc value: 124916.6
Fixed bandwidth: 473.7597 AICc value: Inf
Fixed bandwidth: 474.3398 AICc value: 124917.1
Fixed bandwidth: 473.9813 AICc value: Inf
Fixed bandwidth: 474.2028 AICc value: 124916.8
Fixed bandwidth: 474.0659 AICc value: Inf
Fixed bandwidth: 474.1505 AICc value: 124916.7
Fixed bandwidth: 474.0982 AICc value: 124916.6
Fixed bandwidth: 474.0859 AICc value: Inf
Fixed bandwidth: 474.1058 AICc value: 124916.6
Fixed bandwidth: 474.0935 AICc value: Inf
Fixed bandwidth: 474.1011 AICc value: 124916.6
Fixed bandwidth: 474.0964 AICc value: 124916.6
Fixed bandwidth: 474.0953 AICc value: 124916.6
Fixed bandwidth: 474.0946 AICc value: 124916.6
Fixed bandwidth: 474.0942 AICc value: 124916.6
Fixed bandwidth: 474.0939 AICc value: 124916.6
Fixed bandwidth: 474.0938 AICc value: 124916.6
Fixed bandwidth: 474.0937 AICc value: 124916.6
Fixed bandwidth: 474.0936 AICc value: 124916.6
Fixed bandwidth: 474.0936 AICc value: 124916.6 write_rds(bw_fixed_AIC,"bw_fixed_AIC2.rds")write_rds(bw_fixed_AIC,"data/bw_fixed_AIC.rds")The fixed bandwidth is 474 using CV method
gwr_fixed_AIC <- gwr.basic(
formula = price_per_sqm ~ floor_area_sqm + remaining_lease + floor_level +
dist_to_CBD + dist_to_MRT + dist_busStop +
dist_mall + dist_hawker + dist_supermarket + dist_kindergarten + dist_elderlycare +
park_600m + preschool_350m + kindergartens_350m + supermarket_500m + model_group,
data=mdata,
bw=bw_fixed_AIC,
kernel = 'gaussian',
longlat = FALSE)write_rds(gwr_fixed_AIC,"data/gwr_fixed_AIC.rds")gwr_fixed_AIC <- read_rds("data/gwr_fixed_AIC.rds")
gwr_fixed_AIC ***********************************************************************
* Package GWmodel *
***********************************************************************
Program starts at: 2025-11-18 00:34:42.774258
Call:
gwr.basic(formula = price_per_sqm ~ floor_area_sqm + remaining_lease +
floor_level + dist_to_CBD + dist_to_MRT + dist_busStop +
dist_mall + dist_hawker + dist_supermarket + dist_kindergarten +
dist_elderlycare + park_600m + preschool_350m + kindergartens_350m +
supermarket_500m + model_group, data = mdata, bw = bw_fixed_AIC,
kernel = "gaussian", longlat = FALSE)
Dependent (y) variable: price_per_sqm
Independent variables: floor_area_sqm remaining_lease floor_level dist_to_CBD dist_to_MRT dist_busStop dist_mall dist_hawker dist_supermarket dist_kindergarten dist_elderlycare park_600m preschool_350m kindergartens_350m supermarket_500m model_group
Number of data points: 8679
***********************************************************************
* Results of Global Regression *
***********************************************************************
Call:
lm(formula = formula, data = data)
Residuals:
Min 1Q Median 3Q Max
-3336.7 -487.6 16.5 442.1 4321.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.057e+03 2.226e+02 27.212 < 2e-16 ***
floor_area_sqm -1.582e+01 1.896e+00 -8.345 < 2e-16 ***
remaining_lease 6.505e+01 1.060e+00 61.362 < 2e-16 ***
floor_level2 2.417e+02 2.658e+01 9.092 < 2e-16 ***
floor_level3 4.130e+02 2.673e+01 15.449 < 2e-16 ***
floor_level4 5.445e+02 2.444e+01 22.283 < 2e-16 ***
floor_level5 1.968e+03 4.524e+01 43.510 < 2e-16 ***
dist_to_CBD -2.027e-01 2.339e-03 -86.682 < 2e-16 ***
dist_to_MRT -4.338e-01 2.357e-02 -18.403 < 2e-16 ***
dist_busStop 2.304e-01 1.562e-01 1.475 0.140
dist_mall -2.845e-01 2.615e-02 -10.877 < 2e-16 ***
dist_hawker -2.708e-01 1.904e-02 -14.221 < 2e-16 ***
dist_supermarket 3.309e-01 5.375e-02 6.156 7.78e-10 ***
dist_kindergarten -2.557e-01 5.540e-02 -4.616 3.97e-06 ***
dist_elderlycare -1.077e-01 1.433e-02 -7.510 6.49e-14 ***
park_600m 1.020e+02 7.847e+00 12.996 < 2e-16 ***
preschool_350m -4.535e+01 4.159e+00 -10.906 < 2e-16 ***
kindergartens_350m 1.285e+02 1.155e+01 11.122 < 2e-16 ***
supermarket_500m 1.014e+02 7.633e+00 13.281 < 2e-16 ***
model_group2 -1.192e+01 3.925e+01 -0.304 0.761
model_group3 1.485e+03 7.662e+01 19.374 < 2e-16 ***
model_group4 2.686e+03 1.320e+02 20.345 < 2e-16 ***
---Significance stars
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 751.5 on 8657 degrees of freedom
Multiple R-squared: 0.828
Adjusted R-squared: 0.8276
F-statistic: 1984 on 21 and 8657 DF, p-value: < 2.2e-16
***Extra Diagnostic information
Residual sum of squares: 4889364830
Sigma(hat): 750.6568
AIC: 139600.4
AICc: 139600.5
BIC: 131292.5
***********************************************************************
* Results of Geographically Weighted Regression *
***********************************************************************
*********************Model calibration information*********************
Kernel function: gaussian
Fixed bandwidth: 474.0936
Regression points: the same locations as observations are used.
Distance metric: Euclidean distance metric is used.
****************Summary of GWR coefficient estimates:******************
Min. 1st Qu. Median 3rd Qu.
Intercept -4.2499e+04 3.0270e+03 6.0830e+03 9.4979e+03
floor_area_sqm -2.1775e+03 -4.0017e+01 -3.0245e+01 -1.9190e+01
remaining_lease -1.7573e+03 4.5047e+01 5.6610e+01 7.2157e+01
floor_level2 -5.3191e+02 1.4550e+02 2.3095e+02 3.2052e+02
floor_level3 -2.8664e+02 3.0832e+02 4.0921e+02 5.5928e+02
floor_level4 9.4862e+01 4.5141e+02 6.1325e+02 7.5180e+02
floor_level5 -1.0084e+87 2.1375e+02 1.1017e+03 1.6346e+03
dist_to_CBD -3.1434e+00 -3.6464e-01 -8.1324e-02 1.6157e-01
dist_to_MRT -2.6991e+01 -9.2770e-01 -5.1137e-01 -1.6757e-01
dist_busStop -3.4142e+01 -6.0290e-01 -2.0562e-01 6.5297e-01
dist_mall -2.6852e+01 -7.7525e-01 -2.4515e-01 1.0852e-01
dist_hawker -1.6034e+01 -4.0433e-01 -3.9552e-02 3.8581e-01
dist_supermarket -3.9498e+00 -5.7176e-01 -1.6231e-01 3.0413e-01
dist_kindergarten -2.6437e+01 -2.7762e-01 -9.2100e-03 4.2589e-01
dist_elderlycare -9.4409e+00 -3.1364e-01 -3.5153e-02 2.4378e-01
park_600m -3.0788e+03 -6.2231e+01 1.1352e+01 1.0375e+02
preschool_350m -4.7861e+02 -2.3821e+01 2.8942e+00 2.5947e+01
kindergartens_350m -3.9948e+03 -8.3659e+01 -1.6577e+01 5.3759e+01
supermarket_500m -4.1333e+02 -1.6446e+01 2.0093e+01 6.1090e+01
model_group2 -2.4200e+03 -1.8903e+02 -2.4600e+01 2.0402e+02
model_group3 -9.1300e+03 8.9739e+02 1.5741e+03 2.1804e+03
model_group4 -2.6227e+86 7.4872e+02 4.6558e+03 8.8517e+03
Max.
Intercept 2.7228e+05
floor_area_sqm 7.9484e+01
remaining_lease 1.6597e+02
floor_level2 1.4454e+03
floor_level3 2.0969e+03
floor_level4 2.1917e+03
floor_level5 1.1719e+88
dist_to_CBD 1.2536e+01
dist_to_MRT 1.2168e+01
dist_busStop 1.3305e+01
dist_mall 4.0919e+00
dist_hawker 9.7659e+00
dist_supermarket 3.0568e+01
dist_kindergarten 5.4412e+00
dist_elderlycare 7.1871e+00
park_600m 6.0257e+02
preschool_350m 2.6442e+03
kindergartens_350m 2.5239e+03
supermarket_500m 4.6737e+02
model_group2 9.0806e+04
model_group3 8.9977e+04
model_group4 1.3702e+88
************************Diagnostic information*************************
Number of data points: 8679
Effective number of parameters (2trace(S) - trace(S'S)): 1801.348
Effective degrees of freedom (n-2trace(S) + trace(S'S)): 6877.652
AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 124916.6
AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 122876.4
BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 125918.7
Residual sum of squares: 605344500
R-square value: 0.9787014
Adjusted R-square value: 0.9731222
***********************************************************************
Program stops at: 2025-11-18 00:35:33.23478 bw_adaptive_AIC <- bw.gwr(
formula = price_per_sqm ~ floor_area_sqm + remaining_lease + floor_level +
dist_to_CBD + dist_to_MRT + dist_busStop +
dist_mall + dist_hawker + dist_supermarket + dist_kindergarten + dist_elderlycare +
park_600m + preschool_350m + kindergartens_350m + supermarket_500m + model_group,
data=mdata,
approach = "AIC",
kernel = "gaussian",
adaptive = TRUE,
longlat = FALSE
)Take a cup of tea and have a break, it will take a few minutes.
-----A kind suggestion from GWmodel development group
Adaptive bandwidth (number of nearest neighbours): 5371 AICc value: 138898.6
Adaptive bandwidth (number of nearest neighbours): 3327 AICc value: 138314
Adaptive bandwidth (number of nearest neighbours): 2063 AICc value: 137235
Adaptive bandwidth (number of nearest neighbours): 1282 AICc value: 135390.3
Adaptive bandwidth (number of nearest neighbours): 799 AICc value: 133376.4
Adaptive bandwidth (number of nearest neighbours): 501 AICc value: 131469.3
Adaptive bandwidth (number of nearest neighbours): 316 AICc value: 129829.7
Error in eval(expr, envir) : inv(): matrix is singular
Adaptive bandwidth (number of nearest neighbours): 202 AICc value: Inf
Adaptive bandwidth (number of nearest neighbours): 386 AICc value: 130537.5
Adaptive bandwidth (number of nearest neighbours): 272 AICc value: 129334.5
Adaptive bandwidth (number of nearest neighbours): 245 AICc value: 128992.3
Error in eval(expr, envir) : inv(): matrix is singular
Adaptive bandwidth (number of nearest neighbours): 228 AICc value: Inf
Adaptive bandwidth (number of nearest neighbours): 255 AICc value: 129106.9
Adaptive bandwidth (number of nearest neighbours): 238 AICc value: 128910.7
Adaptive bandwidth (number of nearest neighbours): 234 AICc value: 128860.4
Error in eval(expr, envir) : inv(): matrix is singular
Adaptive bandwidth (number of nearest neighbours): 231 AICc value: Inf
Adaptive bandwidth (number of nearest neighbours): 235 AICc value: 128871
Error in eval(expr, envir) : inv(): matrix is singular
Adaptive bandwidth (number of nearest neighbours): 232 AICc value: Inf
Adaptive bandwidth (number of nearest neighbours): 233 AICc value: 128844.5
Adaptive bandwidth (number of nearest neighbours): 235 AICc value: 128871
Adaptive bandwidth (number of nearest neighbours): 234 AICc value: 128860.4
Adaptive bandwidth (number of nearest neighbours): 234 AICc value: 128860.4
Adaptive bandwidth (number of nearest neighbours): 233 AICc value: 128844.5 gwr_adp_AIC <- gwr.basic(
formula = price_per_sqm ~ floor_area_sqm + remaining_lease + floor_level +
dist_to_CBD + dist_to_MRT + dist_busStop +
dist_mall + dist_hawker + dist_supermarket + dist_kindergarten + dist_elderlycare +
park_600m + preschool_350m + kindergartens_350m + supermarket_500m + model_group,
data=mdata,
bw=bw_adaptive_AIC,
adaptive=TRUE,
kernel = 'gaussian',
longlat = FALSE)
write_rds(gwr_adp_AIC,"data/gwr_adp_AIC.rds")bw_adaptive_CV <- bw.gwr(
formula = price_per_sqm ~ floor_area_sqm + remaining_lease + floor_level +
dist_to_CBD + dist_to_MRT + dist_busStop +
dist_mall + dist_hawker + dist_supermarket + dist_kindergarten + dist_elderlycare +
park_600m + preschool_350m + kindergartens_350m + supermarket_500m + model_group,
data=mdata,
approach = "CV",
kernel = "gaussian",
adaptive = TRUE,
longlat = FALSE
)Take a cup of tea and have a break, it will take a few minutes.
-----A kind suggestion from GWmodel development group
Adaptive bandwidth: 5371 CV score: 4533592483
Adaptive bandwidth: 3327 CV score: 4238072332
Adaptive bandwidth: 2063 CV score: 3743398092
Adaptive bandwidth: 1282 CV score: 3022416694
Adaptive bandwidth: 799 CV score: 2394258363
Adaptive bandwidth: 501 CV score: 1921989904
Adaptive bandwidth: 316 CV score: 1600392492
Adaptive bandwidth: 202 CV score: Inf
Adaptive bandwidth: 386 CV score: 1731368590
Adaptive bandwidth: 272 CV score: 1513529634
Adaptive bandwidth: 245 CV score: 1456872851
Adaptive bandwidth: 228 CV score: Inf
Adaptive bandwidth: 255 CV score: 1475755679
Adaptive bandwidth: 238 CV score: 1443619486
Adaptive bandwidth: 234 CV score: 1435471300
Adaptive bandwidth: 231 CV score: Inf
Adaptive bandwidth: 235 CV score: 1437160574
Adaptive bandwidth: 232 CV score: Inf
Adaptive bandwidth: 233 CV score: 1432912308
Adaptive bandwidth: 235 CV score: 1437160574
Adaptive bandwidth: 234 CV score: 1435471300
Adaptive bandwidth: 234 CV score: 1435471300
Adaptive bandwidth: 233 CV score: 1432912308 gwr_adp_CV <- gwr.basic(
formula = price_per_sqm ~ floor_area_sqm + remaining_lease + floor_level +
dist_to_CBD + dist_to_MRT + dist_busStop +
dist_mall + dist_hawker + dist_supermarket + dist_kindergarten + dist_elderlycare +
park_600m + preschool_350m + kindergartens_350m + supermarket_500m + model_group,
data=mdata,
adaptive=TRUE,
bw=bw_adaptive_CV,
kernel = 'gaussian',
longlat = FALSE)
write_rds(gwr_adp_CV,"data/gwr_adp_CV.rds")gwr_adp_CV ***********************************************************************
* Package GWmodel *
***********************************************************************
Program starts at: 2025-11-18 01:04:08.633028
Call:
gwr.basic(formula = price_per_sqm ~ floor_area_sqm + remaining_lease +
floor_level + dist_to_CBD + dist_to_MRT + dist_busStop +
dist_mall + dist_hawker + dist_supermarket + dist_kindergarten +
dist_elderlycare + park_600m + preschool_350m + kindergartens_350m +
supermarket_500m + model_group, data = mdata, bw = bw_adaptive_CV,
kernel = "gaussian", adaptive = TRUE, longlat = FALSE)
Dependent (y) variable: price_per_sqm
Independent variables: floor_area_sqm remaining_lease floor_level dist_to_CBD dist_to_MRT dist_busStop dist_mall dist_hawker dist_supermarket dist_kindergarten dist_elderlycare park_600m preschool_350m kindergartens_350m supermarket_500m model_group
Number of data points: 8679
***********************************************************************
* Results of Global Regression *
***********************************************************************
Call:
lm(formula = formula, data = data)
Residuals:
Min 1Q Median 3Q Max
-3336.7 -487.6 16.5 442.1 4321.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.057e+03 2.226e+02 27.212 < 2e-16 ***
floor_area_sqm -1.582e+01 1.896e+00 -8.345 < 2e-16 ***
remaining_lease 6.505e+01 1.060e+00 61.362 < 2e-16 ***
floor_level2 2.417e+02 2.658e+01 9.092 < 2e-16 ***
floor_level3 4.130e+02 2.673e+01 15.449 < 2e-16 ***
floor_level4 5.445e+02 2.444e+01 22.283 < 2e-16 ***
floor_level5 1.968e+03 4.524e+01 43.510 < 2e-16 ***
dist_to_CBD -2.027e-01 2.339e-03 -86.682 < 2e-16 ***
dist_to_MRT -4.338e-01 2.357e-02 -18.403 < 2e-16 ***
dist_busStop 2.304e-01 1.562e-01 1.475 0.140
dist_mall -2.845e-01 2.615e-02 -10.877 < 2e-16 ***
dist_hawker -2.708e-01 1.904e-02 -14.221 < 2e-16 ***
dist_supermarket 3.309e-01 5.375e-02 6.156 7.78e-10 ***
dist_kindergarten -2.557e-01 5.540e-02 -4.616 3.97e-06 ***
dist_elderlycare -1.077e-01 1.433e-02 -7.510 6.49e-14 ***
park_600m 1.020e+02 7.847e+00 12.996 < 2e-16 ***
preschool_350m -4.535e+01 4.159e+00 -10.906 < 2e-16 ***
kindergartens_350m 1.285e+02 1.155e+01 11.122 < 2e-16 ***
supermarket_500m 1.014e+02 7.633e+00 13.281 < 2e-16 ***
model_group2 -1.192e+01 3.925e+01 -0.304 0.761
model_group3 1.485e+03 7.662e+01 19.374 < 2e-16 ***
model_group4 2.686e+03 1.320e+02 20.345 < 2e-16 ***
---Significance stars
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 751.5 on 8657 degrees of freedom
Multiple R-squared: 0.828
Adjusted R-squared: 0.8276
F-statistic: 1984 on 21 and 8657 DF, p-value: < 2.2e-16
***Extra Diagnostic information
Residual sum of squares: 4889364830
Sigma(hat): 750.6568
AIC: 139600.4
AICc: 139600.5
BIC: 131292.5
***********************************************************************
* Results of Geographically Weighted Regression *
***********************************************************************
*********************Model calibration information*********************
Kernel function: gaussian
Adaptive bandwidth: 233 (number of nearest neighbours)
Regression points: the same locations as observations are used.
Distance metric: Euclidean distance metric is used.
****************Summary of GWR coefficient estimates:******************
Min. 1st Qu. Median 3rd Qu.
Intercept -1.3085e+02 4.8560e+03 6.7204e+03 8.0603e+03
floor_area_sqm -6.1338e+01 -3.6759e+01 -2.6839e+01 -2.0794e+01
remaining_lease 3.0338e+01 4.7350e+01 5.6395e+01 7.6595e+01
floor_level2 1.0040e+02 1.8831e+02 2.3785e+02 3.1622e+02
floor_level3 2.0715e+02 3.3940e+02 4.4139e+02 5.2263e+02
floor_level4 3.2709e+02 4.8655e+02 6.3054e+02 7.2906e+02
floor_level5 -1.6063e+32 5.7569e+02 1.0076e+03 1.5006e+03
dist_to_CBD -6.1596e-01 -2.4767e-01 -1.1512e-01 -9.3120e-03
dist_to_MRT -2.8322e+00 -1.1376e+00 -5.4537e-01 -2.4390e-01
dist_busStop -1.4849e+00 -5.1144e-01 3.5358e-02 8.2475e-01
dist_mall -1.6161e+00 -4.2553e-01 -1.9301e-01 -2.7310e-02
dist_hawker -6.5848e-01 -3.0984e-01 -1.0289e-01 1.4240e-01
dist_supermarket -1.2167e+00 -4.4819e-01 -1.2668e-01 2.3656e-01
dist_kindergarten -1.3721e+00 -1.3381e-01 1.2038e-01 3.6664e-01
dist_elderlycare -1.1500e+00 -2.4834e-01 -7.1580e-02 2.9221e-02
park_600m -2.0665e+02 -3.7811e+01 6.7861e+00 4.8876e+01
preschool_350m -6.3518e+01 -1.8432e+01 -6.2757e-01 2.2349e+01
kindergartens_350m -2.6952e+02 -4.8270e+01 1.0605e+00 4.5349e+01
supermarket_500m -1.2502e+02 8.3366e+00 4.0508e+01 1.1433e+02
model_group2 -7.0829e+02 -1.0376e+02 4.0179e+01 1.3359e+02
model_group3 3.0897e+02 1.1358e+03 1.4752e+03 1.7917e+03
model_group4 -9.6081e+32 2.3366e+03 3.7418e+03 5.6745e+03
Max.
Intercept 1.4001e+04
floor_area_sqm 2.3735e+01
remaining_lease 1.1118e+02
floor_level2 5.2489e+02
floor_level3 8.4398e+02
floor_level4 9.8722e+02
floor_level5 3.1131e+27
dist_to_CBD 2.5500e-01
dist_to_MRT 6.2420e-01
dist_busStop 3.8852e+00
dist_mall 8.5520e-01
dist_hawker 1.3993e+00
dist_supermarket 1.7628e+00
dist_kindergarten 1.8805e+00
dist_elderlycare 4.9720e-01
park_600m 3.1384e+02
preschool_350m 8.5157e+01
kindergartens_350m 2.8352e+02
supermarket_500m 3.0201e+02
model_group2 5.2912e+02
model_group3 5.5998e+03
model_group4 1.0333e+27
************************Diagnostic information*************************
Number of data points: 8679
Effective number of parameters (2trace(S) - trace(S'S)): 492.9679
Effective degrees of freedom (n-2trace(S) + trace(S'S)): 8186.032
AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 128844.5
AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 128427.1
BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 122816.8
Residual sum of squares: 1298417379
R-square value: 0.9543161
Adjusted R-square value: 0.9515646
***********************************************************************
Program stops at: 2025-11-18 01:05:05.658468 gwr_adp_AIC ***********************************************************************
* Package GWmodel *
***********************************************************************
Program starts at: 2025-11-18 00:52:37.916063
Call:
gwr.basic(formula = price_per_sqm ~ floor_area_sqm + remaining_lease +
floor_level + dist_to_CBD + dist_to_MRT + dist_busStop +
dist_mall + dist_hawker + dist_supermarket + dist_kindergarten +
dist_elderlycare + park_600m + preschool_350m + kindergartens_350m +
supermarket_500m + model_group, data = mdata, bw = bw_adaptive_AIC,
kernel = "gaussian", adaptive = TRUE, longlat = FALSE)
Dependent (y) variable: price_per_sqm
Independent variables: floor_area_sqm remaining_lease floor_level dist_to_CBD dist_to_MRT dist_busStop dist_mall dist_hawker dist_supermarket dist_kindergarten dist_elderlycare park_600m preschool_350m kindergartens_350m supermarket_500m model_group
Number of data points: 8679
***********************************************************************
* Results of Global Regression *
***********************************************************************
Call:
lm(formula = formula, data = data)
Residuals:
Min 1Q Median 3Q Max
-3336.7 -487.6 16.5 442.1 4321.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.057e+03 2.226e+02 27.212 < 2e-16 ***
floor_area_sqm -1.582e+01 1.896e+00 -8.345 < 2e-16 ***
remaining_lease 6.505e+01 1.060e+00 61.362 < 2e-16 ***
floor_level2 2.417e+02 2.658e+01 9.092 < 2e-16 ***
floor_level3 4.130e+02 2.673e+01 15.449 < 2e-16 ***
floor_level4 5.445e+02 2.444e+01 22.283 < 2e-16 ***
floor_level5 1.968e+03 4.524e+01 43.510 < 2e-16 ***
dist_to_CBD -2.027e-01 2.339e-03 -86.682 < 2e-16 ***
dist_to_MRT -4.338e-01 2.357e-02 -18.403 < 2e-16 ***
dist_busStop 2.304e-01 1.562e-01 1.475 0.140
dist_mall -2.845e-01 2.615e-02 -10.877 < 2e-16 ***
dist_hawker -2.708e-01 1.904e-02 -14.221 < 2e-16 ***
dist_supermarket 3.309e-01 5.375e-02 6.156 7.78e-10 ***
dist_kindergarten -2.557e-01 5.540e-02 -4.616 3.97e-06 ***
dist_elderlycare -1.077e-01 1.433e-02 -7.510 6.49e-14 ***
park_600m 1.020e+02 7.847e+00 12.996 < 2e-16 ***
preschool_350m -4.535e+01 4.159e+00 -10.906 < 2e-16 ***
kindergartens_350m 1.285e+02 1.155e+01 11.122 < 2e-16 ***
supermarket_500m 1.014e+02 7.633e+00 13.281 < 2e-16 ***
model_group2 -1.192e+01 3.925e+01 -0.304 0.761
model_group3 1.485e+03 7.662e+01 19.374 < 2e-16 ***
model_group4 2.686e+03 1.320e+02 20.345 < 2e-16 ***
---Significance stars
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 751.5 on 8657 degrees of freedom
Multiple R-squared: 0.828
Adjusted R-squared: 0.8276
F-statistic: 1984 on 21 and 8657 DF, p-value: < 2.2e-16
***Extra Diagnostic information
Residual sum of squares: 4889364830
Sigma(hat): 750.6568
AIC: 139600.4
AICc: 139600.5
BIC: 131292.5
***********************************************************************
* Results of Geographically Weighted Regression *
***********************************************************************
*********************Model calibration information*********************
Kernel function: gaussian
Adaptive bandwidth: 233 (number of nearest neighbours)
Regression points: the same locations as observations are used.
Distance metric: Euclidean distance metric is used.
****************Summary of GWR coefficient estimates:******************
Min. 1st Qu. Median 3rd Qu.
Intercept -1.3085e+02 4.8560e+03 6.7204e+03 8.0603e+03
floor_area_sqm -6.1338e+01 -3.6759e+01 -2.6839e+01 -2.0794e+01
remaining_lease 3.0338e+01 4.7350e+01 5.6395e+01 7.6595e+01
floor_level2 1.0040e+02 1.8831e+02 2.3785e+02 3.1622e+02
floor_level3 2.0715e+02 3.3940e+02 4.4139e+02 5.2263e+02
floor_level4 3.2709e+02 4.8655e+02 6.3054e+02 7.2906e+02
floor_level5 -1.6063e+32 5.7569e+02 1.0076e+03 1.5006e+03
dist_to_CBD -6.1596e-01 -2.4767e-01 -1.1512e-01 -9.3120e-03
dist_to_MRT -2.8322e+00 -1.1376e+00 -5.4537e-01 -2.4390e-01
dist_busStop -1.4849e+00 -5.1144e-01 3.5358e-02 8.2475e-01
dist_mall -1.6161e+00 -4.2553e-01 -1.9301e-01 -2.7310e-02
dist_hawker -6.5848e-01 -3.0984e-01 -1.0289e-01 1.4240e-01
dist_supermarket -1.2167e+00 -4.4819e-01 -1.2668e-01 2.3656e-01
dist_kindergarten -1.3721e+00 -1.3381e-01 1.2038e-01 3.6664e-01
dist_elderlycare -1.1500e+00 -2.4834e-01 -7.1580e-02 2.9221e-02
park_600m -2.0665e+02 -3.7811e+01 6.7861e+00 4.8876e+01
preschool_350m -6.3518e+01 -1.8432e+01 -6.2757e-01 2.2349e+01
kindergartens_350m -2.6952e+02 -4.8270e+01 1.0605e+00 4.5349e+01
supermarket_500m -1.2502e+02 8.3366e+00 4.0508e+01 1.1433e+02
model_group2 -7.0829e+02 -1.0376e+02 4.0179e+01 1.3359e+02
model_group3 3.0897e+02 1.1358e+03 1.4752e+03 1.7917e+03
model_group4 -9.6081e+32 2.3366e+03 3.7418e+03 5.6745e+03
Max.
Intercept 1.4001e+04
floor_area_sqm 2.3735e+01
remaining_lease 1.1118e+02
floor_level2 5.2489e+02
floor_level3 8.4398e+02
floor_level4 9.8722e+02
floor_level5 3.1131e+27
dist_to_CBD 2.5500e-01
dist_to_MRT 6.2420e-01
dist_busStop 3.8852e+00
dist_mall 8.5520e-01
dist_hawker 1.3993e+00
dist_supermarket 1.7628e+00
dist_kindergarten 1.8805e+00
dist_elderlycare 4.9720e-01
park_600m 3.1384e+02
preschool_350m 8.5157e+01
kindergartens_350m 2.8352e+02
supermarket_500m 3.0201e+02
model_group2 5.2912e+02
model_group3 5.5998e+03
model_group4 1.0333e+27
************************Diagnostic information*************************
Number of data points: 8679
Effective number of parameters (2trace(S) - trace(S'S)): 492.9679
Effective degrees of freedom (n-2trace(S) + trace(S'S)): 8186.032
AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 128844.5
AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 128427.1
BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 122816.8
Residual sum of squares: 1298417379
R-square value: 0.9543161
Adjusted R-square value: 0.9515646
***********************************************************************
Program stops at: 2025-11-18 00:53:34.897985 # Assuming your GWR model is called 'gwr_model'
# Extract the SpatialPointsDataFrame with local coefficients
gwr_spdf <- gwr_fixed_CV$SDF
# Convert to sf
gwr_sf <- st_as_sf(gwr_spdf)
# Plot coefficient map for dist_park
ggplot(gwr_sf) +
geom_sf(aes(color = busStop_350m), size = 3) +
scale_color_viridis_c() +
theme_minimal() +
labs(title = "GWR Coefficients: Distance to Park")
The diagnostic metrics reveal a substantial improvement in model performance when incorporating spatial context. The global OLS model explains approximately 73% of the variance in HDB resale prices (R² = 0.73), whereas the geographically weighted approach achieves a remarkable 97% explanatory power (R² = 0.97). This dramatic enhancement demonstrates that accounting for spatial heterogeneity is not merely beneficial but essential for accurately modeling housing prices in Singapore’s complex urban landscape.
While predictive accuracy is important, our primary focus as an explanatory model requires deeper statistical evaluation beyond R-squared alone. For robust model comparison between global and local specifications, we must examine multiple diagnostic criteria: statistical significance (p-values) to identify reliable predictors, information-theoretic measures (Akaike Information Criterion) to assess model parsimony and fit quality, and goodness-of-fit indicators (R-squared) to quantify variance explanation. This comprehensive diagnostic approach allows us to determine whether the increased complexity of the GWR framework yields meaningful explanatory advantages over the conventional global model.
bw_simple <- bw.gwr(
formula = price_per_sqm ~ floor_area_sqm + dist_to_CBD + dist_to_MRT,
data = mdata,
approach = "CV",
kernel = "gaussian",
adaptive = FALSE,
longlat = FALSE
)Take a cup of tea and have a break, it will take a few minutes.
-----A kind suggestion from GWmodel development group
Fixed bandwidth: 22980.96 CV score: 15713868486
Fixed bandwidth: 14205.85 CV score: 15492812959
Fixed bandwidth: 8782.54 CV score: 14940873025
Fixed bandwidth: 5430.749 CV score: 13908054524
Fixed bandwidth: 3359.228 CV score: 12709502259
Fixed bandwidth: 2078.957 CV score: 11348181278
Fixed bandwidth: 1287.706 CV score: 9877848053
Fixed bandwidth: 798.6866 CV score: 8150874748
Fixed bandwidth: 496.4557 CV score: 6304944185
Fixed bandwidth: 309.6667 CV score: 4362703888
Fixed bandwidth: 194.2248 CV score: 4807284631
Fixed bandwidth: 381.0138 CV score: 5180864868
Fixed bandwidth: 265.5719 CV score: 5446367198
Fixed bandwidth: 336.9189 CV score: 4680418149
Fixed bandwidth: 292.824 CV score: 4167776456
Fixed bandwidth: 282.4146 CV score: 4617912492
Fixed bandwidth: 299.2573 CV score: 4241884607
Fixed bandwidth: 288.848 CV score: 4144515394
Fixed bandwidth: 286.3906 CV score: 4130732476
Fixed bandwidth: 284.8719 CV score: 4455664026
Fixed bandwidth: 287.3292 CV score: 4455624529
Fixed bandwidth: 285.8105 CV score: 4804523753
Fixed bandwidth: 286.7492 CV score: 4145642675
Fixed bandwidth: 286.1691 CV score: 150076430449
Fixed bandwidth: 286.5276 CV score: 4107187709
Fixed bandwidth: 286.6122 CV score: 4273587001
Fixed bandwidth: 286.4753 CV score: 4145749169
Fixed bandwidth: 286.5599 CV score: 5749300726
Fixed bandwidth: 286.5076 CV score: 7747469567
Fixed bandwidth: 286.5399 CV score: 12672433275
Fixed bandwidth: 286.5199 CV score: 5824962005
Fixed bandwidth: 286.5323 CV score: 4246346534
Fixed bandwidth: 286.5247 CV score: 6265764393
Fixed bandwidth: 286.5294 CV score: 4166272162
Fixed bandwidth: 286.5265 CV score: 4535058577
Fixed bandwidth: 286.5283 CV score: 4753056701
Fixed bandwidth: 286.5272 CV score: 4702701575
Fixed bandwidth: 286.5278 CV score: 4177748689
Fixed bandwidth: 286.5274 CV score: 9741788250
Fixed bandwidth: 286.5277 CV score: 4292132827
Fixed bandwidth: 286.5275 CV score: 4455701010
Fixed bandwidth: 286.5276 CV score: 4988757336
Fixed bandwidth: 286.5276 CV score: 7324350889 gwr_simple <- gwr.basic(
formula = price_per_sqm ~ floor_area_sqm + dist_to_CBD + dist_to_MRT,
data = mdata,
bw = bw_simple,
kernel = "gaussian",
adaptive = FALSE,
longlat = FALSE
)
gwr_spdf <- gwr_simple$SDF
gwr_sf <- st_as_sf(gwr_spdf)
ggplot(gwr_sf) +
geom_sf(aes(color = dist_to_CBD), size = 0.5) +
scale_color_viridis_c() +
theme_minimal() +
labs(title = "GWR Coefficients: Distance to CBD")
9 Conclusion
The analysis successfully demonstrates that HDB resale prices are driven by a complex interplay of structural and locational factors whose influences vary significantly across Singapore. The Geographically Weighted Regression (GWR) framework substantially outperforms the traditional global OLS model, providing compelling evidence that spatial non-stationarity is a critical feature of Singapore’s housing market.
10 Key Findings
- Significant Factors Influencing HDB Resale Prices The analysis identified multiple statistically significant determinants:
Structural Factors: floor_area_sqm, remaining_lease, floor_level, and model_group all showed strong significant relationships with price, with higher floors and premium flat models commanding substantial premiums.
Locational Factors: dist_to_CBD emerged as the most consistently significant variable, followed by dist_to_MRT, dist_mall, and dist_hawker. Accessibility measures like MRT_600m and busStop_350m also demonstrated significant impacts.
- Evidence of Spatial Non-Stationarity The GWR outputs reveal pronounced spatial variation in coefficient estimates:
Varying Magnitudes: The local coefficients for key variables like dist_to_CBD ranged from -0.616 to 0.255, indicating that the distance penalty varies dramatically across locations.
Direction Changes: Some areas showed positive coefficients for variables that were negative globally, suggesting contextual reversals in relationship directions.
Regional Patterns: The influence of amenities like parks, schools, and public transport showed clear spatial patterning across planning regions.
- Superior Performance of GWR over OLS The model comparison provides strong justification for using geographically weighted methods:
Explanatory Power: GWR achieved an R² of 0.97 compared to OLS’s 0.73, representing a 33% improvement in variance explanation.
Model Fit: The GWR model with adaptive bandwidth (233 neighbors) achieved the lowest AICc (128,844.5) compared to the global model (139,600.5), indicating better fit with appropriate complexity.
Diagnostic Superiority: All GWR variants consistently outperformed the global OLS across multiple metrics including residual sum of squares and adjusted R².
11 Reference
Here are 4 key references formatted as individual points:
Prof TS Kam Geospatial Course - https://isss626-ay2025-26aug.netlify.app/
Gollini, I., Lu, B., Charlton, M., Brunsdon, C., & Harris, P. (2015). GWmodel: An R Package for Exploring Spatial Heterogeneity using Geographically Weighted Models. Journal of Statistical Software, 63(17), 1–50.
Pebesma, E. (2018). Simple Features for R: Standardized Support for Spatial Vector Data. The R Journal, 10(1), 439–446.
Rosen, S. (1974). Hedonic Prices and Implicit Markets: Product Differentiation in Pure Competition. Journal of Political Economy, 82(1), 34–55.
Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002). Geographically weighted regression: the analysis of spatially varying relationships. John Wiley & Sons.