This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.
SPATIAL FOR WHEN I MERGE THINGS
library(sf)## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUElibrary(spdep)## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(spatialreg)## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
##
## Attaching package: 'spatialreg'
##
## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOptiontract_sf <- st_read("C:/Users/evett/Downloads/school_work/IPRO/Chi_nbhd.shp")## Reading layer `Chi_nbhd' from data source
## `C:\Users\evett\Downloads\school_work\IPRO\Chi_nbhd.shp' using driver `ESRI Shapefile'
## Simple feature collection with 98 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -87.94011 ymin: 41.64454 xmax: -87.52414 ymax: 42.02304
## Geodetic CRS: WGS 84chicago_sf <- st_make_valid(tract_sf)
# Create neighborhood adjacency
nb <- poly2nb(chicago_sf, queen = TRUE)
weights <- nb2listw(nb, style = "W", zero.policy = TRUE)library(readxl)
finalidk = spatialpred=read.csv("C:/Users/evett/Downloads/school_work/IPRO/assessor_rm_sli_merged.csv")
finalidk$Private.Service.Line.Material = as.character(finalidk$Private.Service.Line.Material)
finalidk$Public.Service.Line.Material = as.character(finalidk$Public.Service.Line.Material)
finalidk$ward = as.character(finalidk$ward)
finalidk$property_zip <- sub("-.*", "", finalidk$property_zip)
finalidk$property_zip[finalidk$property_zip == 0] <- NA
finalidk$property_zip = as.character(finalidk$property_zip)
finalidk$Age=as.numeric(finalidk$Age)
finalidk$tract_pop=as.numeric(finalidk$tract_pop)
finalidk$tract_white_perc=as.numeric(finalidk$tract_white_perc)
finalidk$tract_black_perc=as.numeric(finalidk$tract_black_perc)
finalidk$tract_asian_perc=as.numeric(finalidk$tract_asian_perc)
finalidk$tract_his_perc=as.numeric(finalidk$tract_his_perc)
finalidk$tract_other_perc=as.numeric(finalidk$tract_other_perc)
finalidk$Rooms=as.numeric(finalidk$Rooms)
finalidk$Bedrooms=as.numeric(finalidk$Bedrooms)
finalidk$ohare_noise = as.character(finalidk$ohare_noise)
finalidk$Wall.Material = as.character(finalidk$Wall.Material)
finalidk$Roof.Material = as.character(finalidk$Roof.Material)
finalidk$Garage.1.Attached = ifelse(finalidk$Garage.1.Attached %in% c("false", "true"), 1, 0)
finalidk$Porch = as.character(finalidk$Porch)
finalidk$Use = as.character(finalidk$Use)
finalidk$Land.Square.Feet=as.numeric(finalidk$Land.Square.Feet)colnames(finalidk)## [1] "Property.Address" "Assessment.Triad"
## [3] "Property.Class" "Modeling.Group"
## [5] "Proportion.Assessed" "Large.Home.Indicator"
## [7] "Wall.Material" "Roof.Material"
## [9] "Basement" "Basement.Finish"
## [11] "Central.Heating" "Central.Air.Conditioning"
## [13] "Fireplaces" "Attic.Type"
## [15] "Design.Plan" "Cathedral.Ceiling"
## [17] "Garage.1.Size" "Garage.1.Material"
## [19] "Garage.1.Attached" "Repair.Condition"
## [21] "Use" "Type.of.Residence"
## [23] "Attic.Finish" "Renovation"
## [25] "Porch" "FEMA.Floodplain"
## [27] "Flood.Risk.Factor" "geoposition"
## [29] "Sale.Price" "Land.Square.Feet"
## [31] "Age" "Apartments"
## [33] "Rooms" "Bedrooms"
## [35] "Full.Baths" "Half.Baths"
## [37] "Garage.1.Area" "Building.Square.Feet"
## [39] "Tract.Median.Income" "Tax.Rate"
## [41] "property_address" "property_city"
## [43] "property_zip" "township"
## [45] "township_name" "nbhd"
## [47] "tract_geoid" "municipality"
## [49] "reps_dist" "senate_dist"
## [51] "ward" "ssa_name"
## [53] "ssa_no" "ohare_noise"
## [55] "floodplain" "fs_flood_risk_direction"
## [57] "withinmr100" "withinmr101300"
## [59] "school_elem_district" "school_hs_district"
## [61] "indicator_has_latlon" "longitude"
## [63] "latitude" "tract_pop"
## [65] "tract_white_perc" "tract_black_perc"
## [67] "tract_asian_perc" "tract_his_perc"
## [69] "tract_other_perc" "tract_midincome"
## [71] "Address" "Private.Service.Line.Material"
## [73] "Public.Service.Line.Material" "Cleaned.Address"
## [75] "sli_count"spatialpred= spatialpred %>% drop_na(longitude, latitude)
predictors_sf <- st_as_sf(spatialpred, coords = c("longitude", "latitude"), crs = 4326) # Ensure correct CRScolnames(chicago_sf)## [1] "pri_neigh" "sec_neigh" "shape_area" "shape_len" "geometry"merged_data <- st_join(predictors_sf, chicago_sf, join = st_within)merged_data <- merged_data %>%
mutate(
Public.Service.Line.Material = case_when(
Public.Service.Line.Material == "LEAD" ~ 1,
Public.Service.Line.Material == "SUSPECTED LEAD" ~ NA,
Public.Service.Line.Material == "NOT LEAD" ~ 0,
Public.Service.Line.Material == "GALVANIZED REQUIRING REPLACEMENT" ~ 1,
TRUE ~ NA_real_
),
Private.Service.Line.Material = case_when(
Private.Service.Line.Material == "LEAD" ~ 1,
Private.Service.Line.Material == "SUSPECTED LEAD" ~ NA,
Private.Service.Line.Material == "NOT LEAD" ~ 0,
Private.Service.Line.Material == "GALVANIZED REQUIRING REPLACEMENT" ~ 1,
TRUE ~ NA_real_
)
)%>%
drop_na(Public.Service.Line.Material, Private.Service.Line.Material)colnames(merged_data)## [1] "Property.Address" "Assessment.Triad"
## [3] "Property.Class" "Modeling.Group"
## [5] "Proportion.Assessed" "Large.Home.Indicator"
## [7] "Wall.Material" "Roof.Material"
## [9] "Basement" "Basement.Finish"
## [11] "Central.Heating" "Central.Air.Conditioning"
## [13] "Fireplaces" "Attic.Type"
## [15] "Design.Plan" "Cathedral.Ceiling"
## [17] "Garage.1.Size" "Garage.1.Material"
## [19] "Garage.1.Attached" "Repair.Condition"
## [21] "Use" "Type.of.Residence"
## [23] "Attic.Finish" "Renovation"
## [25] "Porch" "FEMA.Floodplain"
## [27] "Flood.Risk.Factor" "geoposition"
## [29] "Sale.Price" "Land.Square.Feet"
## [31] "Age" "Apartments"
## [33] "Rooms" "Bedrooms"
## [35] "Full.Baths" "Half.Baths"
## [37] "Garage.1.Area" "Building.Square.Feet"
## [39] "Tract.Median.Income" "Tax.Rate"
## [41] "property_address" "property_city"
## [43] "property_zip" "township"
## [45] "township_name" "nbhd"
## [47] "tract_geoid" "municipality"
## [49] "reps_dist" "senate_dist"
## [51] "ward" "ssa_name"
## [53] "ssa_no" "ohare_noise"
## [55] "floodplain" "fs_flood_risk_direction"
## [57] "withinmr100" "withinmr101300"
## [59] "school_elem_district" "school_hs_district"
## [61] "indicator_has_latlon" "tract_pop"
## [63] "tract_white_perc" "tract_black_perc"
## [65] "tract_asian_perc" "tract_his_perc"
## [67] "tract_other_perc" "tract_midincome"
## [69] "Address" "Private.Service.Line.Material"
## [71] "Public.Service.Line.Material" "Cleaned.Address"
## [73] "sli_count" "pri_neigh"
## [75] "sec_neigh" "shape_area"
## [77] "shape_len" "geometry"convert to logit
# Compute smoothed probabilities (rolling mean)
merged_data$y <- cumsum(merged_data$Public.Service.Line.Material + 0.5) / (seq_along(merged_data$Public.Service.Line.Material) + 1)
merged_data$y1 <- cumsum(merged_data$Private.Service.Line.Material + 0.5) / (seq_along(merged_data$Private.Service.Line.Material) + 1)
# Logit transformation
merged_data$logitpub <- log(merged_data$y / (1 - merged_data$y))## Warning in log(merged_data$y/(1 - merged_data$y)): NaNs produced# Logit transformation
merged_data$logitpriv <- log(merged_data$y1 / (1 - merged_data$y1))## Warning in log(merged_data$y1/(1 - merged_data$y1)): NaNs producedmerged_data <- merged_data %>%
filter(
!is.na(logitpriv), # Remove NA values
!is.na(logitpub),
!is.nan(logitpriv), # Remove NaN values
!is.nan(logitpub),
is.finite(logitpriv), # Remove Inf and -Inf values
is.finite(logitpub)
)colnames(merged_data)## [1] "Property.Address" "Assessment.Triad"
## [3] "Property.Class" "Modeling.Group"
## [5] "Proportion.Assessed" "Large.Home.Indicator"
## [7] "Wall.Material" "Roof.Material"
## [9] "Basement" "Basement.Finish"
## [11] "Central.Heating" "Central.Air.Conditioning"
## [13] "Fireplaces" "Attic.Type"
## [15] "Design.Plan" "Cathedral.Ceiling"
## [17] "Garage.1.Size" "Garage.1.Material"
## [19] "Garage.1.Attached" "Repair.Condition"
## [21] "Use" "Type.of.Residence"
## [23] "Attic.Finish" "Renovation"
## [25] "Porch" "FEMA.Floodplain"
## [27] "Flood.Risk.Factor" "geoposition"
## [29] "Sale.Price" "Land.Square.Feet"
## [31] "Age" "Apartments"
## [33] "Rooms" "Bedrooms"
## [35] "Full.Baths" "Half.Baths"
## [37] "Garage.1.Area" "Building.Square.Feet"
## [39] "Tract.Median.Income" "Tax.Rate"
## [41] "property_address" "property_city"
## [43] "property_zip" "township"
## [45] "township_name" "nbhd"
## [47] "tract_geoid" "municipality"
## [49] "reps_dist" "senate_dist"
## [51] "ward" "ssa_name"
## [53] "ssa_no" "ohare_noise"
## [55] "floodplain" "fs_flood_risk_direction"
## [57] "withinmr100" "withinmr101300"
## [59] "school_elem_district" "school_hs_district"
## [61] "indicator_has_latlon" "tract_pop"
## [63] "tract_white_perc" "tract_black_perc"
## [65] "tract_asian_perc" "tract_his_perc"
## [67] "tract_other_perc" "tract_midincome"
## [69] "Address" "Private.Service.Line.Material"
## [71] "Public.Service.Line.Material" "Cleaned.Address"
## [73] "sli_count" "pri_neigh"
## [75] "sec_neigh" "shape_area"
## [77] "shape_len" "geometry"
## [79] "y" "y1"
## [81] "logitpub" "logitpriv"merged_data <- merged_data %>%
st_drop_geometry() %>% # Drop the POINT geometries
group_by(pri_neigh) %>%
summarise(
across(where(is.numeric), ~ mean(.x, na.rm = TRUE))
) %>%
left_join(select(chicago_sf, pri_neigh, geometry), by = "pri_neigh") %>%
# Reattach polygon geometries
st_as_sf() # Convert back to sf
# Ensure geometries are MULTIPOLYGON
merged_data <- st_cast(merged_data, "MULTIPOLYGON")colnames(merged_data)## [1] "pri_neigh" "Property.Class"
## [3] "Proportion.Assessed" "Fireplaces"
## [5] "FEMA.Floodplain" "Flood.Risk.Factor"
## [7] "Sale.Price" "Land.Square.Feet"
## [9] "Age" "Rooms"
## [11] "Bedrooms" "Full.Baths"
## [13] "Half.Baths" "Garage.1.Area"
## [15] "Building.Square.Feet" "Tract.Median.Income"
## [17] "Tax.Rate" "township"
## [19] "nbhd" "tract_geoid"
## [21] "reps_dist" "senate_dist"
## [23] "ward" "ohare_noise"
## [25] "floodplain" "fs_flood_risk_direction"
## [27] "withinmr100" "withinmr101300"
## [29] "indicator_has_latlon" "tract_pop"
## [31] "tract_white_perc" "tract_black_perc"
## [33] "tract_asian_perc" "tract_his_perc"
## [35] "tract_other_perc" "tract_midincome"
## [37] "Private.Service.Line.Material" "Public.Service.Line.Material"
## [39] "sli_count" "shape_area"
## [41] "shape_len" "y"
## [43] "y1" "logitpub"
## [45] "logitpriv" "geometry"run models
Public_ols_model <- lm(logitpub ~ Age + ward + tract_white_perc, data = merged_data) # Replace with actual variables
summary(Public_ols_model)##
## Call:
## lm(formula = logitpub ~ Age + ward + tract_white_perc, data = merged_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.41831 -0.21368 -0.10364 0.09284 1.66822
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.6939497 0.1282007 13.213 <2e-16 ***
## Age 0.0009647 0.0015630 0.617 0.5386
## ward -0.0072835 0.0031204 -2.334 0.0218 *
## tract_white_perc 0.1220599 0.1465838 0.833 0.4072
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3475 on 90 degrees of freedom
## Multiple R-squared: 0.06173, Adjusted R-squared: 0.03045
## F-statistic: 1.974 on 3 and 90 DF, p-value: 0.1236Private_ols_model <- lm(logitpriv ~ Age + senate_dist + ward, data = merged_data)
# Replace with actual variables
summary(Private_ols_model)##
## Call:
## lm(formula = logitpriv ~ Age + senate_dist + ward, data = merged_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.49034 -0.15810 -0.02499 0.08558 1.02396
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.324075 0.091567 14.460 < 2e-16 ***
## Age 0.000336 0.001062 0.316 0.75237
## senate_dist 0.013697 0.004262 3.214 0.00182 **
## ward -0.003942 0.001785 -2.209 0.02972 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2365 on 90 degrees of freedom
## Multiple R-squared: 0.1518, Adjusted R-squared: 0.1235
## F-statistic: 5.368 on 3 and 90 DF, p-value: 0.001908merged_data <- merged_data %>%
filter(!st_is_empty(geometry)) %>% # Remove rows where geometry is missing
st_cast("MULTIPOLYGON") # Ensure correct geometry type
nb <- poly2nb(merged_data, queen = TRUE)
weights <- nb2listw(nb, style = "W", zero.policy = TRUE)
Public_ols_model <- lm(logitpub ~ tract_white_perc + ward + senate_dist + reps_dist + tract_black_perc, data = merged_data) # Replace with actual variables
Private_ols_model <- lm(logitpriv ~ tract_white_perc + ward + senate_dist + tract_black_perc + reps_dist, data = merged_data)
# Replace with actual variables
#moran.test(residuals(Public_ols_model), weights)
AIC(Public_ols_model)## [1] 56.0155AIC(Private_ols_model)## [1] -7.092734set.seed(123)
# Set of allowed predictors
predictor_vars <- setdiff(names(merged_data), c("logitpub", "logitpriv", "geometry", "pri_neigh", "y", "y1", "Private.Service.Line.Material", "Public.Service.Line.Material", "indicator_has_latlon", "sli_count", "shape_area", "shape_len"))
model_sig_results <- list()
model_count <- 0
attempts <- 0
while (model_count < 100 && attempts < 500) {
attempts <- attempts + 1
num_vars <- sample(2:7, 1)
selected_vars <- sample(predictor_vars, num_vars)
formula_str <- paste("logitpub ~", paste(selected_vars, collapse = " + "))
suppressWarnings(
tryCatch({
model <- lagsarlm(as.formula(formula_str), data = merged_data, listw = weights, Durbin = TRUE)
model_sum <- summary(model)
coef_table <- model_sum$Coef
row_names <- rownames(coef_table)
# Remove intercept & lag terms
is_valid <- !(grepl("Intercept", row_names) | grepl("^lag\\.", row_names))
valid_pvals <- coef_table[is_valid, "Pr(>|z|)"]
sig_count <- sum(valid_pvals < 0.01)
if (sig_count >= 2) {
model_sig_results[[formula_str]] <- sig_count
model_count <- model_count + 1
}
}, error = function(e) {})
)
}
# Sort by number of significant predictors (descending)
sorted_models <- sort(unlist(model_sig_results), decreasing = TRUE)
# Print the formula and number of sig predictors
#print(noquote(formatC(names(sorted_models), width = 100, flag = "-")))
print(sorted_models)## logitpub ~ Land.Square.Feet + Flood.Risk.Factor + Bedrooms + FEMA.Floodplain + Age + township + tract_other_perc
## 3
## logitpub ~ tract_white_perc + reps_dist + Building.Square.Feet + tract_pop + Tax.Rate + Land.Square.Feet
## 2
## logitpub ~ Full.Baths + FEMA.Floodplain + Land.Square.Feet
## 2
## logitpub ~ Building.Square.Feet + floodplain + tract_black_perc + Age
## 2
## logitpub ~ floodplain + Tax.Rate + Full.Baths + ohare_noise + tract_black_perc + withinmr101300 + Half.Baths
## 2
## logitpub ~ Tract.Median.Income + Garage.1.Area + Building.Square.Feet + tract_black_perc + ohare_noise + nbhd
## 2
## logitpub ~ Garage.1.Area + Building.Square.Feet + FEMA.Floodplain + tract_asian_perc
## 2
## logitpub ~ floodplain + Flood.Risk.Factor + Building.Square.Feet + Property.Class
## 2
## logitpub ~ Building.Square.Feet + Land.Square.Feet + tract_asian_perc
## 2
## logitpub ~ Fireplaces + nbhd + Rooms + Property.Class + Land.Square.Feet + Tract.Median.Income + township
## 2
## logitpub ~ tract_black_perc + FEMA.Floodplain + tract_other_perc + Full.Baths
## 2
## logitpub ~ tract_pop + tract_asian_perc + Half.Baths + Full.Baths + floodplain + tract_other_perc + tract_black_perc
## 2
## logitpub ~ Full.Baths + Half.Baths + Flood.Risk.Factor + floodplain + tract_black_perc + FEMA.Floodplain
## 2
## logitpub ~ tract_black_perc + Full.Baths + tract_pop
## 2
## logitpub ~ FEMA.Floodplain + nbhd + Full.Baths + reps_dist + tract_white_perc
## 2
## logitpub ~ tract_white_perc + Bedrooms + reps_dist + Land.Square.Feet + tract_asian_perc
## 2set.seed(123)
# Set of allowed predictors
predictor_vars <- setdiff(names(merged_data), c("logitpub", "logitpriv", "geometry", "pri_neigh", "y", "y1", "Private.Service.Line.Material", "Public.Service.Line.Material", "indicator_has_latlon", "sli_count", "shape_area", "shape_len"))
model_sig_results <- list()
model_count <- 0
attempts <- 0
while (model_count < 100 && attempts < 500) {
attempts <- attempts + 1
num_vars <- sample(2:7, 1)
selected_vars <- sample(predictor_vars, num_vars)
formula_str <- paste("logitpriv ~", paste(selected_vars, collapse = " + "))
suppressWarnings(
tryCatch({
model <- lagsarlm(as.formula(formula_str), data = merged_data, listw = weights, Durbin = TRUE)
model_sum <- summary(model)
coef_table <- model_sum$Coef
row_names <- rownames(coef_table)
# Remove intercept & lag terms
is_valid <- !(grepl("Intercept", row_names) | grepl("^lag\\.", row_names))
valid_pvals <- coef_table[is_valid, "Pr(>|z|)"]
sig_count <- sum(valid_pvals < 0.01)
if (sig_count >= 2) {
model_sig_results[[formula_str]] <- sig_count
model_count <- model_count + 1
}
}, error = function(e) {})
)
}
# Sort by number of significant predictors (descending)
sorted_models <- sort(unlist(model_sig_results), decreasing = TRUE)
# Print the formula and number of sig predictors
#print(noquote(formatC(names(sorted_models), width = 100, flag = "-")))
print(sorted_models)## logitpriv ~ Land.Square.Feet + Flood.Risk.Factor + Bedrooms + FEMA.Floodplain + Age + township + tract_other_perc
## 3
## logitpriv ~ Fireplaces + nbhd + Rooms + Property.Class + Land.Square.Feet + Tract.Median.Income + township
## 3
## logitpriv ~ tract_white_perc + Bedrooms + reps_dist + Land.Square.Feet + tract_asian_perc
## 3
## logitpriv ~ tract_asian_perc + Sale.Price + Land.Square.Feet
## 2
## logitpriv ~ floodplain + tract_pop + Fireplaces + Building.Square.Feet + Age
## 2
## logitpriv ~ tract_white_perc + reps_dist + Building.Square.Feet + tract_pop + Tax.Rate + Land.Square.Feet
## 2
## logitpriv ~ Full.Baths + FEMA.Floodplain + Land.Square.Feet
## 2
## logitpriv ~ Building.Square.Feet + floodplain + tract_black_perc + Age
## 2
## logitpriv ~ Tax.Rate + withinmr100 + tract_midincome + Proportion.Assessed + tract_asian_perc + Building.Square.Feet
## 2
## logitpriv ~ Garage.1.Area + Building.Square.Feet + FEMA.Floodplain + tract_asian_perc
## 2
## logitpriv ~ township + Land.Square.Feet + Garage.1.Area + Tract.Median.Income + Half.Baths + Building.Square.Feet
## 2
## logitpriv ~ floodplain + Flood.Risk.Factor + Building.Square.Feet + Property.Class
## 2
## logitpriv ~ Half.Baths + Bedrooms + tract_black_perc + tract_white_perc + tract_his_perc + tract_pop + Property.Class
## 2
## logitpriv ~ Building.Square.Feet + Land.Square.Feet + tract_asian_perc
## 2
## logitpriv ~ township + Proportion.Assessed + withinmr101300 + tract_asian_perc + tract_black_perc + Full.Baths
## 2
## logitpriv ~ tract_pop + tract_asian_perc + Half.Baths + Full.Baths + floodplain + tract_other_perc + tract_black_perc
## 2
## logitpriv ~ Land.Square.Feet + Garage.1.Area + reps_dist + Rooms + senate_dist + Half.Baths
## 2
## logitpriv ~ Proportion.Assessed + nbhd + Full.Baths + tract_midincome + tract_pop
## 2
## logitpriv ~ Rooms + Flood.Risk.Factor + Age + tract_pop + tract_black_perc + nbhd
## 2
## logitpriv ~ tract_asian_perc + Age + Flood.Risk.Factor + Half.Baths + reps_dist + tract_black_perc + tract_his_perc
## 2
## logitpriv ~ tract_black_perc + Full.Baths + tract_pop
## 2
## logitpriv ~ tract_asian_perc + withinmr100 + Building.Square.Feet + tract_his_perc + reps_dist + Land.Square.Feet + Rooms
## 2
## logitpriv ~ Land.Square.Feet + tract_asian_perc + Half.Baths + FEMA.Floodplain + tract_his_perc + township + Bedrooms
## 2
## logitpriv ~ FEMA.Floodplain + nbhd + Full.Baths + reps_dist + tract_white_perc
## 2Public
set.seed(123)
# Exclude response and unwanted variables
predictor_vars <- setdiff(names(merged_data), c("logitpub", "logitpriv", "geometry", "pri_neigh", "y", "y1"))
# Storage
all_models_lm <- list()
aic_values <- numeric()
failed_models <- list()
i <- 0
model_attempts <- 0
while (i < 100 && model_attempts < 500) {
model_attempts <- model_attempts + 1
# Randomly select 1–5 predictors
num_vars <- sample(2:7, 1)
random_vars <- sample(predictor_vars, num_vars, replace = FALSE)
formula_str <- paste("logitpub ~", paste(random_vars, collapse = " + "))
suppressWarnings(
tryCatch({
#Change the type of models
model <- lagsarlm(as.formula(formula_str), data = merged_data, listw = weights, Durbin = TRUE)
# Skip numerical Hessian models
if (!is.null(model$ase) && model$ase == "Numerical Hessian used") {
failed_models[[formula_str]] <- "Numerical Hessian used"
} else {
all_models_lm[[formula_str]] <- summary(model)
aic_values[formula_str] <- AIC(model)
i <- i + 1
}
}, error = function(e) {
failed_models[[formula_str]] <- "Error"
})
)
}
# Sort and extract top 5
sorted_aic_values <- sort(aic_values)
top_5_models <- all_models_lm[names(sorted_aic_values[1:5])]
top_5_aics <- sorted_aic_values[1:5]
top_5_aics## logitpub ~ senate_dist + Rooms + tract_other_perc + Full.Baths + FEMA.Floodplain + Land.Square.Feet
## -14.9187164
## logitpub ~ Bedrooms + Land.Square.Feet + sli_count + Private.Service.Line.Material + ohare_noise + tract_other_perc + floodplain
## -12.2031273
## logitpub ~ fs_flood_risk_direction + Private.Service.Line.Material + Fireplaces + Land.Square.Feet + FEMA.Floodplain
## -2.3058166
## logitpub ~ Building.Square.Feet + Public.Service.Line.Material + floodplain + tract_white_perc
## 0.4674571
## logitpub ~ fs_flood_risk_direction + Tax.Rate + floodplain + Full.Baths + reps_dist + Fireplaces + tract_pop
## 7.1577523Private
set.seed(123)
# Exclude response and unwanted variables
predictor_vars <- setdiff(names(merged_data), c("logitpub", "logitpriv", "geometry", "pri_neigh", "y", "y1"))
# Storage
all_models_lm <- list()
aic_values <- numeric()
failed_models <- list()
i <- 0
model_attempts <- 0
while (i < 100 && model_attempts < 500) {
model_attempts <- model_attempts + 1
# Randomly select 1–5 predictors
num_vars <- sample(2:7, 1)
random_vars <- sample(predictor_vars, num_vars, replace = FALSE)
formula_str <- paste("logitpriv ~", paste(random_vars, collapse = " + "))
suppressWarnings(
tryCatch({
model <- lagsarlm(as.formula(formula_str), data = merged_data, listw = weights, Durbin = TRUE)
# Skip numerical Hessian models
if (!is.null(model$ase) && model$ase == "Numerical Hessian used") {
failed_models[[formula_str]] <- "Numerical Hessian used"
} else {
all_models_lm[[formula_str]] <- summary(model)
aic_values[formula_str] <- AIC(model)
i <- i + 1
}
}, error = function(e) {
failed_models[[formula_str]] <- "Error"
})
)
}
# Sort and extract top 5
sorted_aic_values <- sort(aic_values)
top_5_models <- all_models_lm[names(sorted_aic_values[1:5])]
top_5_aics <- sorted_aic_values[1:5]
top_5_aics## logitpriv ~ senate_dist + Rooms + tract_other_perc + Full.Baths + FEMA.Floodplain + Land.Square.Feet
## -74.02306
## logitpriv ~ Bedrooms + Land.Square.Feet + sli_count + Private.Service.Line.Material + ohare_noise + tract_other_perc + floodplain
## -68.09521
## logitpriv ~ fs_flood_risk_direction + Private.Service.Line.Material + Fireplaces + Land.Square.Feet + FEMA.Floodplain
## -61.29176
## logitpriv ~ Building.Square.Feet + Public.Service.Line.Material + floodplain + tract_white_perc
## -60.71233
## logitpriv ~ fs_flood_risk_direction + Tax.Rate + floodplain + Full.Baths + reps_dist + Fireplaces + tract_pop
## -57.33318PUB - SDM
Pub_SAR_model <- lagsarlm(logitpub ~ Land.Square.Feet + Flood.Risk.Factor + Bedrooms + FEMA.Floodplain + Age + township + tract_white_perc, data = merged_data, listw = weights)
summary(Pub_SAR_model)##
## Call:lagsarlm(formula = logitpub ~ Land.Square.Feet + Flood.Risk.Factor +
## Bedrooms + FEMA.Floodplain + Age + township + tract_white_perc,
## data = merged_data, listw = weights)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.509815 -0.147700 -0.021794 0.073496 1.353493
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.5245e+00 1.0789e+00 1.4130 0.1577
## Land.Square.Feet -4.2977e-08 3.0258e-05 -0.0014 0.9989
## Flood.Risk.Factor -6.6034e-03 2.1827e-02 -0.3025 0.7622
## Bedrooms 2.8343e-02 3.2568e-02 0.8703 0.3842
## FEMA.Floodplain -8.0694e+00 1.1225e+01 -0.7189 0.4722
## Age -5.7513e-04 1.0939e-03 -0.5258 0.5991
## township -1.7097e-02 1.3995e-02 -1.2217 0.2218
## tract_white_perc 1.0498e-01 9.2339e-02 1.1369 0.2556
##
## Rho: 0.77904, LR test value: 46.61, p-value: 8.6634e-12
## Asymptotic standard error: 0.067534
## z-value: 11.535, p-value: < 2.22e-16
## Wald statistic: 133.07, p-value: < 2.22e-16
##
## Log likelihood: -6.526481 for lag model
## ML residual variance (sigma squared): 0.056652, (sigma: 0.23802)
## Number of observations: 93
## Number of parameters estimated: 10
## AIC: 33.053, (AIC for lm: 77.663)
## LM test for residual autocorrelation
## test value: 0.1316, p-value: 0.71678Pub_SEM_model <- errorsarlm(logitpub ~ Land.Square.Feet + Flood.Risk.Factor + Bedrooms + FEMA.Floodplain + Age + township + tract_white_perc, data = merged_data, listw = weights)
summary(Pub_SEM_model)##
## Call:errorsarlm(formula = logitpub ~ Land.Square.Feet + Flood.Risk.Factor +
## Bedrooms + FEMA.Floodplain + Age + township + tract_white_perc,
## data = merged_data, listw = weights)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.4281100 -0.1220387 -0.0047681 0.0689265 1.1193451
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.8692e+00 1.7785e+00 1.0510 0.293253
## Land.Square.Feet -8.8382e-05 3.7917e-05 -2.3310 0.019756
## Flood.Risk.Factor -7.7267e-03 2.0410e-02 -0.3786 0.705001
## Bedrooms 7.0496e-02 3.1123e-02 2.2651 0.023508
## FEMA.Floodplain -3.0080e+01 9.2919e+00 -3.2372 0.001207
## Age -1.0729e-03 9.6609e-04 -1.1106 0.266749
## township -2.8057e-03 2.3972e-02 -0.1170 0.906830
## tract_white_perc 5.4720e-01 1.5780e-01 3.4676 0.000525
##
## Lambda: 0.8957, LR test value: 63.293, p-value: 1.7764e-15
## Asymptotic standard error: 0.041021
## z-value: 21.835, p-value: < 2.22e-16
## Wald statistic: 476.78, p-value: < 2.22e-16
##
## Log likelihood: 1.815261 for error model
## ML residual variance (sigma squared): 0.04309, (sigma: 0.20758)
## Number of observations: 93
## Number of parameters estimated: 10
## AIC: 16.369, (AIC for lm: 77.663)Pub_SDM_model <- lagsarlm(logitpub ~ Land.Square.Feet + Flood.Risk.Factor + Bedrooms + FEMA.Floodplain + Age + township + tract_white_perc, data = merged_data, listw = weights, Durbin = TRUE)
summary(Pub_SDM_model)##
## Call:lagsarlm(formula = logitpub ~ Land.Square.Feet + Flood.Risk.Factor +
## Bedrooms + FEMA.Floodplain + Age + township + tract_white_perc,
## data = merged_data, listw = weights, Durbin = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.553142 -0.083872 -0.008002 0.079308 0.577838
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.8235e+00 1.2881e+00 -1.4157 0.156867
## Land.Square.Feet -1.4532e-04 3.5436e-05 -4.1009 4.116e-05
## Flood.Risk.Factor -1.1192e-02 1.8902e-02 -0.5921 0.553781
## Bedrooms 8.8497e-02 2.9188e-02 3.0320 0.002429
## FEMA.Floodplain -8.7491e+00 9.5740e+00 -0.9138 0.360801
## Age 4.2347e-04 9.4397e-04 0.4486 0.653716
## township 1.0792e-02 2.3880e-02 0.4519 0.651312
## tract_white_perc 6.3368e-01 1.5065e-01 4.2064 2.595e-05
## lag.Land.Square.Feet 2.2096e-04 4.9731e-05 4.4431 8.868e-06
## lag.Flood.Risk.Factor -7.6335e-02 4.3937e-02 -1.7374 0.082319
## lag.Bedrooms -8.2404e-02 5.1418e-02 -1.6026 0.109014
## lag.FEMA.Floodplain 1.4397e+02 2.3680e+01 6.0799 1.203e-09
## lag.Age 1.6097e-03 1.9505e-03 0.8253 0.409223
## lag.township 1.9127e-02 3.1726e-02 0.6029 0.546600
## lag.tract_white_perc -7.3911e-01 1.8499e-01 -3.9954 6.460e-05
##
## Rho: 0.62943, LR test value: 36.121, p-value: 1.854e-09
## Asymptotic standard error: 0.084491
## z-value: 7.4497, p-value: 9.3481e-14
## Wald statistic: 55.498, p-value: 9.3592e-14
##
## Log likelihood: 21.20113 for mixed model
## ML residual variance (sigma squared): 0.033591, (sigma: 0.18328)
## Number of observations: 93
## Number of parameters estimated: 17
## AIC: -8.4023, (AIC for lm: 25.719)
## LM test for residual autocorrelation
## test value: 3.0131, p-value: 0.082593PRI - SDM
Priv_SAR_model <- lagsarlm(logitpriv ~ Land.Square.Feet + Flood.Risk.Factor + Bedrooms + FEMA.Floodplain + Age + township + tract_white_perc, data = merged_data, listw = weights)
summary(Priv_SAR_model)##
## Call:lagsarlm(formula = logitpriv ~ Land.Square.Feet + Flood.Risk.Factor +
## Bedrooms + FEMA.Floodplain + Age + township + tract_white_perc,
## data = merged_data, listw = weights)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.346571 -0.114883 -0.011503 0.068667 0.848440
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0354e+00 7.5998e-01 1.3624 0.1731
## Land.Square.Feet 1.0654e-06 2.1233e-05 0.0502 0.9600
## Flood.Risk.Factor -2.7035e-03 1.5312e-02 -0.1766 0.8599
## Bedrooms 2.6138e-02 2.2861e-02 1.1433 0.2529
## FEMA.Floodplain 1.4878e+00 7.8738e+00 0.1890 0.8501
## Age -3.7078e-04 7.6767e-04 -0.4830 0.6291
## township -1.1391e-02 9.8187e-03 -1.1601 0.2460
## tract_white_perc 6.8784e-02 6.4709e-02 1.0630 0.2878
##
## Rho: 0.78411, LR test value: 49.477, p-value: 2.0075e-12
## Asymptotic standard error: 0.066459
## z-value: 11.799, p-value: < 2.22e-16
## Wald statistic: 139.21, p-value: < 2.22e-16
##
## Log likelihood: 26.26215 for lag model
## ML residual variance (sigma squared): 0.027899, (sigma: 0.16703)
## Number of observations: 93
## Number of parameters estimated: 10
## AIC: -32.524, (AIC for lm: 14.952)
## LM test for residual autocorrelation
## test value: 4.9339e-05, p-value: 0.9944Priv_SEM_model <- errorsarlm(logitpriv ~ Land.Square.Feet + Flood.Risk.Factor + Bedrooms + FEMA.Floodplain + Age + township + tract_white_perc, data = merged_data, listw = weights)
summary(Priv_SEM_model)##
## Call:errorsarlm(formula = logitpriv ~ Land.Square.Feet + Flood.Risk.Factor +
## Bedrooms + FEMA.Floodplain + Age + township + tract_white_perc,
## data = merged_data, listw = weights)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3210619 -0.0889912 0.0013805 0.0652986 0.7352221
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4046e+00 1.2666e+00 1.1089 0.2674568
## Land.Square.Feet -5.7810e-05 2.7123e-05 -2.1314 0.0330558
## Flood.Risk.Factor -2.6802e-03 1.4612e-02 -0.1834 0.8544680
## Bedrooms 5.3297e-02 2.2284e-02 2.3917 0.0167697
## FEMA.Floodplain -1.2730e+01 6.6567e+00 -1.9123 0.0558390
## Age -7.7675e-04 6.9195e-04 -1.1226 0.2616283
## township -6.2593e-04 1.7083e-02 -0.0366 0.9707714
## tract_white_perc 3.8084e-01 1.1268e-01 3.3798 0.0007255
##
## Lambda: 0.89015, LR test value: 63.24, p-value: 1.7764e-15
## Asymptotic standard error: 0.042558
## z-value: 20.916, p-value: < 2.22e-16
## Wald statistic: 437.48, p-value: < 2.22e-16
##
## Log likelihood: 33.14406 for error model
## ML residual variance (sigma squared): 0.022096, (sigma: 0.14865)
## Number of observations: 93
## Number of parameters estimated: 10
## AIC: -46.288, (AIC for lm: 14.952)Priv_SDM_model <- lagsarlm(logitpriv ~ Land.Square.Feet + Flood.Risk.Factor + Bedrooms + FEMA.Floodplain + Age + township + tract_white_perc, data = merged_data, listw = weights, Durbin = TRUE)
summary(Priv_SDM_model)##
## Call:lagsarlm(formula = logitpriv ~ Land.Square.Feet + Flood.Risk.Factor +
## Bedrooms + FEMA.Floodplain + Age + township + tract_white_perc,
## data = merged_data, listw = weights, Durbin = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3385658 -0.0669832 0.0006872 0.0588156 0.4373693
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1566e+00 9.4386e-01 -1.2254 0.2204201
## Land.Square.Feet -9.5817e-05 2.5942e-05 -3.6935 0.0002212
## Flood.Risk.Factor -5.5927e-03 1.3838e-02 -0.4042 0.6861007
## Bedrooms 6.8485e-02 2.1369e-02 3.2049 0.0013510
## FEMA.Floodplain 1.4735e-01 6.9319e+00 0.0213 0.9830407
## Age 2.9132e-04 6.9123e-04 0.4215 0.6734223
## township 1.1364e-02 1.7476e-02 0.6503 0.5155210
## tract_white_perc 4.5527e-01 1.1028e-01 4.1283 3.654e-05
## lag.Land.Square.Feet 1.4541e-04 3.6376e-05 3.9974 6.405e-05
## lag.Flood.Risk.Factor -5.6204e-02 3.2161e-02 -1.7476 0.0805356
## lag.Bedrooms -5.7359e-02 3.7598e-02 -1.5256 0.1271140
## lag.FEMA.Floodplain 8.6678e+01 1.7344e+01 4.9975 5.809e-07
## lag.Age 1.5626e-03 1.4285e-03 1.0939 0.2740118
## lag.township 8.4048e-03 2.3222e-02 0.3619 0.7174024
## lag.tract_white_perc -5.3677e-01 1.3549e-01 -3.9617 7.443e-05
##
## Rho: 0.64001, LR test value: 37.505, p-value: 9.1171e-10
## Asymptotic standard error: 0.084159
## z-value: 7.6048, p-value: 2.8644e-14
## Wald statistic: 57.833, p-value: 2.8533e-14
##
## Log likelihood: 50.01726 for mixed model
## ML residual variance (sigma squared): 0.018001, (sigma: 0.13417)
## Number of observations: 93
## Number of parameters estimated: 17
## AIC: -66.035, (AIC for lm: -30.529)
## LM test for residual autocorrelation
## test value: 3.5862, p-value: 0.058263AIC(Public_ols_model)## [1] 56.0155for public
# Load the libraries
library(spdep)
library(sf)
library(tmap)
library(ggplot2)
library(RColorBrewer)
chicago_sf <- st_read("C:/Users/evett/Downloads/school_work/IPRO/Chi_nbhd.shp")## Reading layer `Chi_nbhd' from data source
## `C:\Users\evett\Downloads\school_work\IPRO\Chi_nbhd.shp' using driver `ESRI Shapefile'
## Simple feature collection with 98 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -87.94011 ymin: 41.64454 xmax: -87.52414 ymax: 42.02304
## Geodetic CRS: WGS 84# Ensure valid geometries
chicago_sf <- st_make_valid(chicago_sf)
merged_sf <- st_as_sf(merged_data, coords = c("longitude", "latitude"), crs = 4326)
merged_sf <- st_transform(merged_sf, crs = st_crs(chicago_sf))
nb_merged <- poly2nb(merged_sf, queen = TRUE)
weights_merged <- nb2listw(nb_merged, style = "W", zero.policy = TRUE)
# Add predictions to the spatial data
#Public
merged_data$PUB_SDM_pred <- as.numeric(predict(Pub_SDM_model, newdata = merged_sf, listw = weights_merged))## Warning in predict.Sarlm(Pub_SDM_model, newdata = merged_sf, listw =
## weights_merged): only legacy.mixed=TRUE is supported for pred.type='TS' and
## mixed models. legacy.mixed is forcedmerged_data$PUB_SAR_pred <- as.numeric(predict(Pub_SAR_model, newdata = merged_data, listw = weights))
merged_data$PUB_SEM_pred <- as.numeric(predict(Pub_SEM_model, newdata = merged_data, listw = weights))
#Private
merged_data$PRIV_SDM_pred <- as.numeric(predict(Priv_SDM_model, newdata = merged_sf, listw = weights_merged))## Warning in predict.Sarlm(Priv_SDM_model, newdata = merged_sf, listw =
## weights_merged): only legacy.mixed=TRUE is supported for pred.type='TS' and
## mixed models. legacy.mixed is forcedmerged_data$PRIV_SAR_pred <- as.numeric(predict(Priv_SAR_model, newdata = merged_data, listw = weights))
merged_data$PRIV_SEM_pred <- as.numeric(predict(Priv_SEM_model, newdata = merged_data, listw = weights))
# Define inverse logit function
inv_logit <- function(x) {
exp(x) / (1 + exp(x))
}
# Predict logit values and transform back to probabilities
merged_data$Pred_Probpubsar <- inv_logit(merged_data$PUB_SAR_pred)
merged_data$Pred_Probpubsem <- inv_logit(merged_data$PUB_SEM_pred)
merged_data$Pred_Probpubsdm <- inv_logit(merged_data$PUB_SDM_pred)
merged_data$Pred_Probprivsar <- inv_logit(merged_data$PRIV_SAR_pred)
merged_data$Pred_Probprivsem <- inv_logit(merged_data$PRIV_SEM_pred)
merged_data$Pred_Probprivsdm <- inv_logit(merged_data$PRIV_SDM_pred)
# Add residuals
#merged_data$SAR_resid <- as.numeric(residuals(sar_model))
#merged_data$SEM_resid <- as.numeric(residuals(sem_model))
#merged_data$SDM_resid <- as.numeric(residuals(sdm_model))library(sf)
merged_sf <- st_as_sf(merged_data, coords = c("longitude", "latitude"), crs = 4326)public
library(tmap)
# Map the predictions
tmap_mode("view")## ℹ tmap mode set to "view".tm_shape(merged_sf) +
tm_polygons(
"Pred_Probpubsar",
fill.legend = tm_legend(title = "SAR Predictions"),
palette = "brewer.PuBuGn"
) +
tm_layout(
title = "Public SAR Model Predictions",
legend.outside = TRUE
) +
tm_borders()##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "brewer.PuBuGn" is
## named "pu_bu_gn" (in long format "brewer.pu_bu_gn")# SEM Predictions Map
tmap_mode("view")## ℹ tmap mode set to "view".tm_shape(merged_sf) +
tm_polygons(
fill = "Pred_Probpubsem",
fill.legend = tm_legend(title = "SEM Predictions"),
palette = "brewer.YlOrRd"
) +
tm_layout(
title = "Public SEM Model Predictions",
legend.outside = TRUE
) +
tm_borders()##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "brewer.YlOrRd" is
## named "yl_or_rd" (in long format "brewer.yl_or_rd")# SDM Predictions Map
library(spData)
library(tmap)
tmap_mode("view")## ℹ tmap mode set to "view".tm = tm_shape(merged_sf) +
tm_polygons("Pred_Probpubsdm",
fill.scale = tm_scale_intervals(
values = "brewer.purples"),
fill.legend = tm_legend(title = "Predicted Probability")) +
tm_borders() +
tm_title("Public SDM Model Predictions")
tm#tmap_save(tm, filename = "Public_SR.html", selfcontained = TRUE)Private
# SAR Prediction Map
tmap_mode("view")## ℹ tmap mode set to "view".tm_shape(merged_sf) +
tm_polygons(
fill = "Pred_Probprivsar",
fill.legend = tm_legend(title = "SAR Predictions"),
palette = "brewer.PuBuGn"
) +
tm_layout(
title = "Private SAR Model Predictions",
legend.outside = TRUE
) +
tm_borders()##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "brewer.PuBuGn" is
## named "pu_bu_gn" (in long format "brewer.pu_bu_gn")