# If needed, install the regular packages once with:
# install.packages(c("sf", "tidyverse", "tmap", "spdep", "spatialreg",
#                    "rmarkdown", "knitr", "scales"))
# On Windows with R 4.6, install INLA once from the testing binary repository:
# if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
# BiocManager::install(c("graph", "Rgraphviz"), ask = FALSE, update = FALSE)
# install.packages("INLA",
#   repos = c(CRAN = "https://cloud.r-project.org",
#             INLA = "https://inla.r-inla-download.org/R/testing"),
#   dependencies = TRUE,
#   type = "binary")

#load libraries
library(sf)
library(tidyverse)
library(tmap)
library(spdep)
library(spatialreg)
library(INLA)

#tract data
acs_tract_nc <-
st_read("https://drive.google.com/uc?export=download&id=1cnz4xgdDRZvlXzN3IyvCxOETRWfrpw-0")
## Reading layer `acs_tract_data_nc_2020_2024' from data source 
##   `https://drive.google.com/uc?export=download&id=1cnz4xgdDRZvlXzN3IyvCxOETRWfrpw-0' 
##   using driver `GeoJSON'
## Simple feature collection with 2498 features and 47 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32187 ymin: 33.86641 xmax: -75.46062 ymax: 36.58812
## Geodetic CRS:  NAD83

Executive summary

North Carolina’s household income pattern is not spatially random. High-income and low-income census tracts form recognizable geographic clusters, and a conventional regression leaves meaningful spatial structure unexplained. The spatial models therefore give planners a more realistic picture than OLS alone: local socioeconomic conditions matter, but neighboring context and broader regional processes matter too.

Why geography matters

North Carolina has made progress on poverty overall, but that progress is uneven. The North Carolina Department of Commerce reported that the state poverty rate fell from 17.8% in 2013 to 12.8% in 2023, while county rates still ranged from 7.2% in Wake County to 28.6% in Scotland County. This statewide improvement alongside large local gaps shows why averages can hide the places with the greatest need (North Carolina Department of Commerce, 2025).

The longer-run economic context also matters. State labor-market research found that wage inequality intensified during the 1990s and early 2000s and improved little after the Great Recession (North Carolina Department of Commerce, 2022). For regional planners, concentrated low income can signal overlapping needs involving housing affordability, transportation access, workforce development, education, and public services.

This report uses 2020–2024 American Community Survey census-tract estimates. The dependent variable is median household income. The predictors capture poverty, renting, educational attainment, and population density. These are associations, not proof that one factor causes another.

Method-to-data check. These observations are census-tract polygons, so the analysis uses tract adjacency, Moran’s I, LISA, and areal spatial regression. Point-pattern tools such as KDE, average nearest neighbor, and Ripley’s L are not appropriate here because the outcome is a value attached to an area rather than a set of event locations.

Descriptive patterns

The analytic file contains 2,498 census tracts. Tract median household income averages $75,118, with a median of $67,929. The middle half of tracts ranges from $53,430 to $88,167.

Distribution of median household income across North Carolina census tracts
Tracts Mean SD Minimum 25th percentile Median 75th percentile Maximum
2498 $75,118 $32,850 $14,821 $53,430 $67,929 $88,167 $250,001

The map shows that income inequality is both regional and local. Higher-income tracts are especially visible around major metropolitan and suburban areas, while lower-income tracts appear across parts of eastern and southeastern North Carolina and within some cities. Using an unweighted median across tracts within each county, the five highest are Chatham County, Wake County, Union County, Orange County, Camden County; the five lowest are Tyrrell County, Robeson County, Washington County, Bladen County, Lenoir County. These are descriptive tract summaries rather than official county-level median-income estimates.

Global spatial dependence

The null hypothesis for Global Moran’s I is that tract incomes are spatially random under the chosen queen-contiguity structure. A positive, statistically significant value means similar income levels tend to occur near one another.

set.seed(12345)
global_moran <- moran.mc(
  nc$median_hh_inc,
  listw,
  nsim = 999,
  zero.policy = TRUE
)

global_moran_table <- tibble(
  Statistic = "Global Moran's I",
  Estimate = unname(global_moran$statistic),
  `Permutation p-value` = global_moran$p.value,
  Permutations = 999
)

knitr::kable(
  global_moran_table |>
    mutate(
      Estimate = fmt_num(Estimate),
      `Permutation p-value` = ifelse(`Permutation p-value` < .001, "< .001", fmt_num(`Permutation p-value`))
    ),
  align = c("l", "r", "r", "r")
)
Statistic Estimate Permutation p-value Permutations
Global Moran’s I 0.585 0.001 999

Global Moran’s I is 0.585 with a permutation p-value = .001. I reject the null of spatial randomness. Tracts with similar household incomes are located near one another much more often than would be expected by chance. That finding also means a model that assumes every tract is independent is likely incomplete.

Local clusters (LISA)

Global Moran’s I tells us that clustering exists, but not where it occurs. Local Indicators of Spatial Association (LISA) identify four types of local pattern:

  • High–High: a high-income tract surrounded by high-income neighbors.
  • Low–Low: a low-income tract surrounded by low-income neighbors.
  • High–Low: a high-income tract surrounded by lower-income neighbors.
  • Low–High: a low-income tract surrounded by higher-income neighbors.
set.seed(12345)
local_i <- localmoran_perm(
  nc$median_hh_inc,
  listw,
  nsim = 999,
  zero.policy = TRUE,
  alternative = "two.sided",
  iseed = 12345
)

sim_p_name <- if ("Pr(z != E(Ii)) Sim" %in% colnames(local_i)) {
  "Pr(z != E(Ii)) Sim"
} else {
  grep("Pr.*Sim", colnames(local_i), value = TRUE)[1]
}

if (is.na(sim_p_name) || length(sim_p_name) == 0) {
  stop("Could not identify the permutation p-value returned by localmoran_perm().")
}

income_z <- as.numeric(scale(nc$median_hh_inc))
lag_income_z <- lag.listw(listw, income_z, zero.policy = TRUE)
local_p <- local_i[, sim_p_name]

nc$lisa_cluster <- case_when(
  local_p < 0.05 & income_z > 0 & lag_income_z > 0 ~ "High–High",
  local_p < 0.05 & income_z < 0 & lag_income_z < 0 ~ "Low–Low",
  local_p < 0.05 & income_z > 0 & lag_income_z < 0 ~ "High–Low",
  local_p < 0.05 & income_z < 0 & lag_income_z > 0 ~ "Low–High",
  TRUE ~ "Not significant"
)

nc$lisa_cluster <- factor(
  nc$lisa_cluster,
  levels = c("High–High", "Low–Low", "High–Low", "Low–High", "Not significant")
)

lisa_counts <- nc |>
  st_drop_geometry() |>
  count(lisa_cluster, name = "Tracts") |>
  mutate(Percent = 100 * Tracts / sum(Tracts))

Number of tracts in each LISA category
lisa_cluster Tracts Percent
High–High 259 10.4%
Low–Low 306 12.2%
High–Low 12 0.5%
Low–High 17 0.7%
Not significant 1904 76.2%

The LISA map makes the statewide result more specific. High–High clusters are concentrated heavily in and around the state’s largest metropolitan economies, especially Wake and Mecklenburg counties and nearby suburban counties. Low–Low clusters are more geographically widespread, including parts of the northeast and southeast as well as lower-income areas inside some metropolitan counties. The High–Low and Low–High categories are spatial outliers; they may be useful for local follow-up because they differ from their immediate surroundings.

LISA performs a separate test for every tract, so some significant locations can occur by chance. I treat the map as an exploratory planning tool rather than a final designation of neighborhood need.

Ordinary least squares model

The OLS model relates tract median household income to four tract characteristics: poverty rate, renter share, bachelor’s-degree attainment, and the log of population density. The density measure is logged because the original distribution is strongly right-skewed. Each percentage variable is measured in percentage points.

ols_formula <- median_hh_inc ~
  pct_under_poverty_level +
  pct_renter +
  pct_bachelors_or_more +
  log_pop_density

ols_model <- lm(ols_formula, data = st_drop_geometry(nc))

ols_table <- as.data.frame(coef(summary(ols_model))) |>
  rownames_to_column("Term") |>
  as_tibble() |>
  rename(
    Estimate = Estimate,
    `Standard error` = `Std. Error`,
    `t value` = `t value`,
    `p value` = `Pr(>|t|)`
  ) |>
  mutate(
    Term = recode(
      Term,
      `(Intercept)` = "Intercept",
      pct_under_poverty_level = "Poverty rate (percentage points)",
      pct_renter = "Renter share (percentage points)",
      pct_bachelors_or_more = "Bachelor's degree or higher (percentage points)",
      log_pop_density = "Log population density"
    )
  )

knitr::kable(
  ols_table |>
    mutate(
      across(c(Estimate, `Standard error`), ~ fmt_dollar(.x)),
      `t value` = fmt_num(`t value`, 2),
      `p value` = ifelse(`p value` < .001, "< .001", fmt_num(`p value`))
    ),
  align = c("l", "r", "r", "r", "r"),
  caption = paste0("OLS estimates (adjusted R-squared = ", fmt_num(summary(ols_model)$adj.r.squared), ")")
)
OLS estimates (adjusted R-squared = 0.731)
Term Estimate Standard error t value p value
Intercept $50,897 $1,671 30.47 < .001
Poverty rate (percentage points) -$621 $47 -13.36 < .001
Renter share (percentage points) -$611 $24 -25.40 < .001
Bachelor’s degree or higher (percentage points) $945 $23 40.38 < .001
Log population density $3,365 $321 10.49 < .001

Holding the other included variables constant, a one-percentage-point increase in the poverty rate is associated with -$621 in median household income. A one-point increase in renter share is associated with -$611, while a one-point increase in bachelor’s-degree attainment is associated with $945. These estimates describe conditional relationships within this dataset; they should not be interpreted as causal effects.

Do the OLS residuals remain spatial?

set.seed(12345)
ols_resid_moran <- moran.mc(
  residuals(ols_model),
  listw,
  nsim = 999,
  zero.policy = TRUE
)

knitr::kable(
  tibble(
    Statistic = "Moran's I for OLS residuals",
    Estimate = fmt_num(unname(ols_resid_moran$statistic)),
    `Permutation p-value` = ifelse(
      ols_resid_moran$p.value < .001,
      "< .001",
      fmt_num(ols_resid_moran$p.value)
    )
  ),
  align = c("l", "r", "r")
)
Statistic Estimate Permutation p-value
Moran’s I for OLS residuals 0.271 0.001

The OLS residual Moran’s I is 0.271 with p = .001. The residuals are still positively clustered, so the predictors explain a large share of income differences but do not satisfy the OLS independence assumption. Standard OLS inference may therefore be too confident, and the model is missing spatial structure.

Choosing a traditional spatial model

A spatial lag model (SLM) is designed for dependence in the outcome itself: a tract’s income is related to neighboring tract incomes. A spatial error model (SEM) is designed for spatial dependence in omitted or unmeasured processes. Robust Rao’s score diagnostics and model fit are used together to choose between them.

# Current spdep versions use lm.RStests(); the fallback supports older versions.
if ("lm.RStests" %in% getNamespaceExports("spdep")) {
  rs_tests <- lm.RStests(
    ols_model,
    listw,
    test = c("RSerr", "RSlag", "adjRSerr", "adjRSlag"),
    zero.policy = TRUE
  )
  robust_error_name <- "adjRSerr"
  robust_lag_name <- "adjRSlag"
} else {
  rs_tests <- lm.LMtests(
    ols_model,
    listw,
    test = c("LMerr", "LMlag", "RLMerr", "RLMlag"),
    zero.policy = TRUE
  )
  robust_error_name <- "RLMerr"
  robust_lag_name <- "RLMlag"
}

rs_table <- bind_rows(lapply(names(rs_tests), function(nm) {
  test <- rs_tests[[nm]]
  tibble(
    Test = nm,
    Statistic = unname(test$statistic),
    `p value` = test$p.value
  )
}))

slm_model <- lagsarlm(
  ols_formula,
  data = st_drop_geometry(nc),
  listw = listw,
  method = "Matrix",
  zero.policy = TRUE
)

sem_model <- errorsarlm(
  ols_formula,
  data = st_drop_geometry(nc),
  listw = listw,
  method = "Matrix",
  zero.policy = TRUE
)

p_lag <- rs_tests[[robust_lag_name]]$p.value
p_error <- rs_tests[[robust_error_name]]$p.value

# If only one robust diagnostic is significant, follow that diagnostic.
# If both (or neither) are significant, use the lower AIC as the tie-breaker.
if (p_lag < 0.05 && p_error >= 0.05) {
  selected_model_name <- "Spatial Lag Model (SLM)"
  spatial_model <- slm_model
} else if (p_error < 0.05 && p_lag >= 0.05) {
  selected_model_name <- "Spatial Error Model (SEM)"
  spatial_model <- sem_model
} else if (AIC(slm_model) <= AIC(sem_model)) {
  selected_model_name <- "Spatial Lag Model (SLM)"
  spatial_model <- slm_model
} else {
  selected_model_name <- "Spatial Error Model (SEM)"
  spatial_model <- sem_model
}

knitr::kable(
  rs_table |>
    mutate(
      Statistic = fmt_num(Statistic, 2),
      `p value` = ifelse(`p value` < .001, "< .001", fmt_num(`p value`))
    ),
  align = c("l", "r", "r"),
  caption = "Rao's score diagnostics for spatial dependence"
)
Rao’s score diagnostics for spatial dependence
Test Statistic p value
RSerr 499.36 < .001
RSlag 474.04 < .001
adjRSerr 122.56 < .001
adjRSlag 97.24 < .001 
selection_table <- tibble(
  Candidate = c("Spatial Lag Model (SLM)", "Spatial Error Model (SEM)"),
  AIC = c(AIC(slm_model), AIC(sem_model)),
  `Spatial parameter` = c(slm_model$rho, sem_model$lambda)
)

knitr::kable(
  selection_table |>
    mutate(
      AIC = fmt_num(AIC, 1),
      `Spatial parameter` = fmt_num(`Spatial parameter`)
    ),
  align = c("l", "r", "r"),
  caption = "Traditional spatial-model comparison"
)
Traditional spatial-model comparison
Candidate AIC Spatial parameter
Spatial Lag Model (SLM) 55,335.8 0.380
Spatial Error Model (SEM) 55,367.4 0.535

The robust lag-test p-value is < .001, and the robust error-test p-value is < .001. When both tests are significant, neither diagnostic alone provides a clean choice, so I use AIC as a tie-breaker. The selected model is the Spatial Lag Model (SLM) because it has the stronger combined diagnostic and fit evidence.

spatial_summary <- summary(spatial_model)
spatial_coef <- as.data.frame(spatial_summary$Coef) |>
  rownames_to_column("Term") |>
  as_tibble()

colnames(spatial_coef)[2:5] <- c("Estimate", "Standard error", "z value", "p value")

spatial_coef <- spatial_coef |>
  mutate(
    Term = recode(
      Term,
      `(Intercept)` = "Intercept",
      pct_under_poverty_level = "Poverty rate (percentage points)",
      pct_renter = "Renter share (percentage points)",
      pct_bachelors_or_more = "Bachelor's degree or higher (percentage points)",
      log_pop_density = "Log population density"
    )
  )

knitr::kable(
  spatial_coef |>
    mutate(
      across(c(Estimate, `Standard error`), ~ fmt_dollar(.x)),
      `z value` = fmt_num(`z value`, 2),
      `p value` = ifelse(`p value` < .001, "< .001", fmt_num(`p value`))
    ),
  align = c("l", "r", "r", "r", "r"),
  caption = paste("Coefficient estimates from the selected", selected_model_name)
)
Coefficient estimates from the selected Spatial Lag Model (SLM)
Term Estimate Standard error z value p value
Intercept $32,804 $1,711 19.18 < .001
Poverty rate (percentage points) -$445 $43 -10.41 < .001
Renter share (percentage points) -$547 $22 -24.99 < .001
Bachelor’s degree or higher (percentage points) $680 $24 28.06 < .001
Log population density $2,408 $293 8.23 < .001

The estimated spatial-lag parameter (rho) is 0.380. Its positive value means that, even after adjusting for the included tract characteristics, median household income remains connected to income in neighboring tracts. The regression coefficients are conditional/direct associations; feedback through neighboring tracts means the total regional association can be larger than a coefficient read in isolation.

Bayesian BYM2 model

The Bayesian BYM2 model separates residual variation into two pieces: a spatially structured component shared by neighboring tracts and an unstructured tract-specific component. This avoids forcing all unexplained variation into only one type of error and produces stabilized estimates for each tract. Income is modeled in $10,000 units for numerical stability; all reported fixed and spatial effects are converted back to dollars.

graph_file <- tempfile(fileext = ".adj")
nb2INLA(graph_file, nb)
nc_graph <- inla.read.graph(graph_file)

bym2_formula <- income_10k ~
  pct_under_poverty_level +
  pct_renter +
  pct_bachelors_or_more +
  log_pop_density +
  f(
    tract_id,
    model = "bym2",
    graph = nc_graph,
    scale.model = TRUE,
    constr = TRUE,
    adjust.for.con.comp = TRUE
  )

bym2_model <- inla(
  bym2_formula,
  data = st_drop_geometry(nc),
  family = "gaussian",
  control.predictor = list(compute = TRUE),
  control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE),
  verbose = FALSE
)

bym2_fixed <- bym2_model$summary.fixed |>
  as.data.frame() |>
  rownames_to_column("Term") |>
  as_tibble() |>
  select(Term, mean, sd, `0.025quant`, `0.975quant`) |>
  mutate(
    across(c(mean, sd, `0.025quant`, `0.975quant`), ~ .x * 10000),
    Term = recode(
      Term,
      `(Intercept)` = "Intercept",
      pct_under_poverty_level = "Poverty rate (percentage points)",
      pct_renter = "Renter share (percentage points)",
      pct_bachelors_or_more = "Bachelor's degree or higher (percentage points)",
      log_pop_density = "Log population density"
    )
  )

knitr::kable(
  bym2_fixed |>
    mutate(across(-Term, ~ fmt_dollar(.x))) |>
    rename(
      `Posterior mean` = mean,
      `Posterior SD` = sd,
      `2.5%` = `0.025quant`,
      `97.5%` = `0.975quant`
    ),
  align = c("l", "r", "r", "r", "r"),
  caption = "BYM2 fixed effects and 95% credible intervals"
)
BYM2 fixed effects and 95% credible intervals
Term Posterior mean Posterior SD 2.5% 97.5%
Intercept $64,226 $2,319 $59,691 $68,786
Poverty rate (percentage points) -$445 $44 -$531 -$359
Renter share (percentage points) -$604 $23 -$649 -$559
Bachelor’s degree or higher (percentage points) $873 $29 $816 $929
Log population density $1,229 $398 $448 $2,009 
hyper <- bym2_model$summary.hyperpar |>
  as.data.frame() |>
  rownames_to_column("Parameter") |>
  as_tibble()

phi_row <- grep("^Phi for", hyper$Parameter)
phi_mean <- if (length(phi_row) > 0) hyper$mean[phi_row[1]] else NA_real_

knitr::kable(
  hyper |>
    select(Parameter, mean, sd, `0.025quant`, `0.975quant`) |>
    mutate(across(-Parameter, ~ fmt_num(.x))) |>
    rename(
      `Posterior mean` = mean,
      `Posterior SD` = sd,
      `2.5%` = `0.025quant`,
      `97.5%` = `0.975quant`
    ),
  align = c("l", "r", "r", "r", "r"),
  caption = "BYM2 hyperparameters"
)
Interpretable BYM2 spatial-structure result
Measure Posterior mean 95% credible interval
Share of residual tract variation that is spatially structured (phi) 53.1% 43.6% to 62.2% 
# INLA stores the n combined BYM2 tract effects first, followed by n internal
# component values. Retain the combined effect aligned to the observed tracts.
nc$bym2_spatial_effect <-
  bym2_model$summary.random$tract_id$mean[seq_len(nrow(nc))] * 10000

The posterior mean of phi is 53.1% (95% credible interval: 43.6% to 62.2%). This means that about half of the remaining tract-level variation is geographically structured and about half is tract-specific. The map displays their combined effect.

Comparing the three approaches

predictor_names <- c(
  "pct_under_poverty_level",
  "pct_renter",
  "pct_bachelors_or_more",
  "log_pop_density"
)

term_labels <- c(
  pct_under_poverty_level = "Poverty rate",
  pct_renter = "Renter share",
  pct_bachelors_or_more = "Bachelor's degree or higher",
  log_pop_density = "Log population density"
)

comparison <- tibble(
  Term = unname(term_labels[predictor_names]),
  OLS = unname(coef(ols_model)[predictor_names]),
  Traditional_spatial = unname(coef(spatial_model)[predictor_names]),
  BYM2 = bym2_model$summary.fixed[predictor_names, "mean"] * 10000
) |>
  mutate(across(-Term, ~ fmt_dollar(.x))) |>
  rename(`Selected traditional spatial model` = Traditional_spatial)

knitr::kable(
  comparison,
  align = c("l", "r", "r", "r"),
  caption = "Coefficient comparison across OLS, the selected traditional spatial model, and BYM2"
)
Coefficient comparison across OLS, the selected traditional spatial model, and BYM2
Term OLS Selected traditional spatial model BYM2
Poverty rate -$621 -$445 -$445
Renter share -$611 -$547 -$604
Bachelor’s degree or higher $945 $680 $873
Log population density $3,365 $2,408 $1,229 
rmse <- function(actual, predicted) sqrt(mean((actual - predicted)^2, na.rm = TRUE))

fit_comparison <- tibble(
  Model = c("OLS", selected_model_name, "Bayesian BYM2"),
  `In-sample RMSE` = c(
    rmse(nc$median_hh_inc, fitted(ols_model)),
    rmse(nc$median_hh_inc, fitted(spatial_model)),
    rmse(nc$median_hh_inc, bym2_model$summary.fitted.values$mean * 10000)
  ),
  `Model-specific fit statistic` = c(
    paste0("AIC ", fmt_num(AIC(ols_model), 1)),
    paste0("AIC ", fmt_num(AIC(spatial_model), 1)),
    paste0("WAIC ", fmt_num(bym2_model$waic$waic, 1))
  )
)

knitr::kable(
  fit_comparison |>
    mutate(`In-sample RMSE` = fmt_dollar(`In-sample RMSE`)),
  align = c("l", "r", "l"),
  caption = "Overall model comparison"
)
Model comparison using interpretable, framework-appropriate evidence
Model Key evidence What it shows
OLS Adjusted R-squared = 0.731; residual Moran’s I = 0.271 Strong covariate fit, but significant spatial structure remains.
Spatial Lag Model (SLM) rho = 0.380; AIC = 55,335.8 (433.5 lower than OLS) Better accounts for dependence between neighboring tract incomes.
Bayesian BYM2 Phi = 53.1%; 95% credible interval 43.6% to 62.2% Separates remaining variation into spatially structured and tract-specific portions.

The three models answer related but slightly different questions:

  • OLS gives the simplest associations, but its residuals are spatially clustered. It treats neighboring tracts as if they were independent, which the diagnostic test rejects.
  • The selected traditional spatial model explicitly represents one form of spatial dependence. Its coefficients should be read after accounting for the neighboring outcome pattern (SLM) or spatially correlated omitted influences (SEM).
  • BYM2 allows both structured and unstructured residual variation. Its fixed effects can be compared with OLS for direction and approximate magnitude, while its tract effects show where unmeasured geographic context remains.

The coefficient table should not be used to declare one estimate “correct” solely because it is smaller or larger. The main substantive question is whether the direction and practical meaning of the associations remain consistent after spatial structure is acknowledged. AIC is used only to compare OLS with the traditional spatial models, while phi and its credible interval summarize the BYM2 result. A direct predictive comparison would require holding out data or using cross-validation.

Policy implications

Bottom line for planners: Income inequality in North Carolina is spatially concentrated. Statewide policy matters, but uniform statewide delivery is unlikely to be enough. Regional and neighborhood targeting should complement universal programs.

  1. Target clusters, not isolated tracts alone. Low–Low clusters cross tract boundaries. Housing, transportation, workforce, education, and health investments should be coordinated across neighboring jurisdictions rather than assigned one tract at a time.

  2. Protect access in rapidly changing high-income clusters. High–High areas can bring jobs and tax capacity, but they may also create housing-cost pressure. Planners should connect economic growth with affordable housing, transit access, and anti-displacement strategies.

  3. Investigate spatial outliers locally. Low–High tracts may represent pockets of disadvantage within prosperous regions, while High–Low tracts may be locally advantaged but disconnected from surrounding communities. Both patterns can be missed in county averages.

  4. Use the BYM2 residual map to guide follow-up. Large remaining spatial effects suggest that the available ACS variables do not tell the whole story. Those locations are priorities for adding information on employment centers, housing costs, transportation, school access, rurality, and historical investment.

  5. Do not interpret the coefficients causally. The results identify patterns useful for planning and hypothesis generation. They do not show that changing one tract characteristic by itself would automatically produce the estimated income change.

Limitations

ACS tract estimates contain sampling error, and the supplied dataset does not include margins of error. Census tracts also vary in population and land area. Results depend partly on the queen-contiguity definition of a neighbor, and the minimum number of nearest-neighbor links was added to connect otherwise disconnected coastal components so every tract could remain in the analysis. LISA results involve many local tests and should be treated as exploratory. Finally, the analysis is cross-sectional and observational, so it cannot establish causation or show how income patterns changed over time.

Conclusion

Median household income is strongly spatially patterned across North Carolina. Both the global and local statistics show concentrated advantage and disadvantage, and the OLS residual test confirms that observable tract characteristics do not remove this geographic dependence. The traditional spatial model improves on OLS by representing neighboring structure directly, while the BYM2 model shows that residual variation contains both geographic and tract-specific components. For planners, the practical lesson is that income inequality is not only about individual tracts—it is also a regional system that calls for coordinated, place-aware policy.

References

North Carolina Department of Commerce. (2022, September 13). Wage growth and wage inequality in North Carolina. https://www.commerce.nc.gov/blog/2022/09/13/wage-growth-and-wage-inequality-north-carolina

North Carolina Department of Commerce. (2025, February 20). North Carolina’s decade of progress in poverty reduction. https://www.commerce.nc.gov/news/the-lead-feed/nc-poverty-reduction

U.S. Census Bureau. (n.d.). American Community Survey (ACS). https://www.census.gov/programs-surveys/acs.html

## Report generated with R 4.5.3.