Welcome to the Lincoln Institute of Land Policy GWR Analysis Tool. We hope you find this R script useful.

Equity is a critical component of a property tax system. An equitable assessment system requires that all properties are assessed accurately and consistently, regardless of their price, location, size, style and any other attribute of a property. Since location is a key driver of real estate value, understanding how well assessments capture locational influences is critical to determining whether assessments are equitable.

This is a spatial analysis tool for assessors to identify locations with potential assessment inequity. Such areas may be flagged for additional quality control review, including data correction, updating, or even valuation adjustments.

Please note: the analysis provided by this script does not prove or confirm the existence of assessment inequity. Additionally, the algorithm is susceptible to extreme values (outliers) that may distort results. Assessment inequity should always be tested using ratio studies consistent with the standards provided by the International Association of Assessing Officers.

# Load necessary libraries
library(leaflet)
library(sp)
library(GWmodel)
## Loading required package: robustbase
## Loading required package: Rcpp
## Welcome to GWmodel version 2.4-2.
library(htmlwidgets)
library(viridis)
## Warning: package 'viridis' was built under R version 4.4.3
## Loading required package: viridisLite
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# Set working directory and load data
setwd("C:/Users/pbida/OneDrive/Desktop/R")
data <- read.csv("ward 20 with x-y.csv")

# Create SpatialPointsDataFrame from coordinates
coords <- cbind(data$Longitude, data$Latitude)
spData <- SpatialPointsDataFrame(coords, data = data)

# Define palette for A.S.Ratio and build map for Ratio
pal <- colorNumeric(palette = "RdYlGn", domain = data$A.S.Ratio)
map_ratio <- leaflet(spData) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addCircleMarkers(
    lng = ~coords[,1],
    lat = ~coords[,2],
    color = ~pal(A.S.Ratio),
    popup = ~paste(
      "<b>Parcel ID:</b>", data$Parcel.Id, "<br>",
      "<b>Assessed Value:</b> $", formatC(data$Current.Assessed.Value, big.mark = ",", format = "f", digits = 0), "<br>",
      "<b>Sale Price:</b> $", formatC(data$Sale.Price, big.mark = ",", format = "f", digits = 0), "<br>",
      "<b>Ratio (Assessed/Sale):</b>", round(A.S.Ratio, 2)
    ),
    layerId = ~data$Parcel.Id,
    fillOpacity = 0.8,
    stroke = FALSE
  ) %>%
  addLegend(pal = pal, values = ~A.S.Ratio, title = "Ratio (Assessed/Sale)", opacity = 1) %>%
  onRender("
    function(el, x) {
      var map = this;
      map.eachLayer(function(layer) {
        if (layer instanceof L.CircleMarker) {
          layer.on('click', function(e) {
            map.setView(e.latlng, 20);
          });
        }
      });
    }
  ")

# Print map for Ratio
map_ratio
# Run GWR model
DM <- gw.dist(dp.locat = coordinates(spData))
bw1 <- bw.gwr(Sale.Price ~ 0 + Current.Assessed.Value,
              data = spData,
              adaptive = TRUE, 
              kernel = "exponential", 
              dMat = DM)
## Adaptive bandwidth: 338 CV score: 1.776027e+12 
## Adaptive bandwidth: 217 CV score: 1.771337e+12 
## Adaptive bandwidth: 141 CV score: 1.767606e+12 
## Adaptive bandwidth: 95 CV score: 1.762855e+12 
## Adaptive bandwidth: 65 CV score: 1.757158e+12 
## Adaptive bandwidth: 48 CV score: 1.751124e+12 
## Adaptive bandwidth: 36 CV score: 1.737922e+12 
## Adaptive bandwidth: 30 CV score: 1.737462e+12 
## Adaptive bandwidth: 25 CV score: 1.7351e+12 
## Adaptive bandwidth: 23 CV score: 1.734449e+12 
## Adaptive bandwidth: 20 CV score: 1.735896e+12 
## Adaptive bandwidth: 23 CV score: 1.734449e+12
gwr.res1 <- gwr.basic(Sale.Price ~ 0 + Current.Assessed.Value, 
                      data = spData, 
                      adaptive = TRUE, 
                      bw = bw1, 
                      kernel = "exponential", 
                      dMat = DM)

# Convert GWR results to a data frame and filter insignificant observations 
gwr.df <- as.data.frame(gwr.res1$SDF)
gwr.df_filt <- gwr.df %>% filter(abs(Current.Assessed.Value_TV) > 1.99)

# Build map for local R² using filtered GWR results
pal_r2 <- colorNumeric(palette = viridis(100), domain = gwr.df_filt$Local_R2)
map_local_r2 <- leaflet(gwr.df_filt) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addCircleMarkers(
    lng = ~coords.x1, 
    lat = ~coords.x2,
    color = ~pal_r2(Local_R2),
    popup = ~paste(
      "<b>Parcel ID:</b>", formatC(data$Parcel.Id, format = "f", digits = 0), "<br>",
      "<b>Local R²:</b>", formatC(Local_R2, format = "f", digits = 2), "<br>",
      "<b>Assessment Coefficient:</b>", round(Current.Assessed.Value, 2)
    ),
    fillOpacity = 0.8,
    stroke = FALSE,
    layerId = ~data$Parcel.Id
  ) %>%
  addLegend(pal = pal_r2, values = ~Local_R2, title = "Local R²", opacity = 1) %>%
  onRender("
    function(el, x) {
      var map = this;
      map.eachLayer(function(layer) {
        if (layer instanceof L.CircleMarker) {
          layer.on('click', function(e) {
            map.setView(e.latlng, 20);
          });
        }
      });
    }
  ")

# Print map for Local R²
map_local_r2
# Build map for Current.Assessed.Value using RdYlGn palette
pal_assessed <- colorNumeric(palette = "RdYlGn", domain = gwr.df_filt$Current.Assessed.Value)
map_assessed <- leaflet(gwr.df_filt) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addCircleMarkers(
    lng = ~coords.x1, 
    lat = ~coords.x2,
    color = ~pal_assessed(Current.Assessed.Value),
    popup = ~paste(
      "<b>Parcel ID:</b>", formatC(data$Parcel.Id, format = "f", digits = 0), "<br>",
      "<b>Assessment Coefficient:</b>", round(Current.Assessed.Value, 2), "<br>",
      "<b>Local R²:</b>", formatC(Local_R2, format = "f", digits = 2)
    ),
    fillOpacity = 0.8,
    stroke = FALSE,
    layerId = ~data$Parcel.Id
  ) %>%
  addLegend(pal = pal_assessed, values = ~Current.Assessed.Value, title = "Current Assessed Value", opacity = 1) %>%
  onRender("
    function(el, x) {
      var map = this;
      map.eachLayer(function(layer) {
        if (layer instanceof L.CircleMarker) {
          layer.on('click', function(e) {
            map.setView(e.latlng, 20);
          });
        }
      });
    }
  ")

# Print map for Current.Assessed.Value
map_assessed