# Required once if any package is missing:
# install.packages(c("sf", "tidyverse", "spdep", "spatialreg",
#                    "gridExtra", "knitr", "rmarkdown"))
# Optional, for map scale bars and north arrows:
# install.packages("ggspatial")

library(sf)
library(tidyverse)
library(spdep)
library(spatialreg)
library(gridExtra)

# The source file is cached locally on first knit. Reading the Google Drive
# URL on every knit is fragile: Drive intermittently returns an HTML
# virus-scan interstitial instead of the GeoJSON, which breaks the build for
# anyone re-running this document.
data_url  <- "https://drive.google.com/uc?export=download&id=1hxz7coJASXKz-yGC1Z7yVQPsqf8mIBzE"
data_path <- file.path("data", "segregation_tracts.geojson")

if (!file.exists(data_path)) {
  dir.create("data", showWarnings = FALSE)
  download.file(data_url, data_path, mode = "wb", quiet = TRUE)
}

data <- st_read(data_path, quiet = TRUE)

# METRO_NAME is truncated at 30 characters in the source file (for example
# "Minneapolis-St. Paul-Bloomingt"). The truncation is in the supplied data,
# not a typo, so the filters below match the truncated strings exactly.
metro_san_diego   <- "San Diego-Carlsbad, CA"
metro_minneapolis <- "Minneapolis-St. Paul-Bloomingt"
# ---------------------------------------------------------------------
# 1. Connected-component detection (iterative; avoids the recursion depth
#    problem in older versions of spdep::n.comp.nb).
# ---------------------------------------------------------------------
nb_components <- function(nb) {
  n <- length(nb)
  membership <- integer(n)
  component_number <- 0L

  for (start in seq_len(n)) {
    if (membership[start] != 0L) next

    component_number <- component_number + 1L
    membership[start] <- component_number

    queue <- integer(n)
    queue[1] <- start
    first <- 1L
    last  <- 1L

    while (first <= last) {
      current <- queue[first]
      first <- first + 1L

      neighbors <- as.integer(nb[[current]])
      neighbors <- unique(
        neighbors[!is.na(neighbors) & neighbors >= 1L & neighbors <= n]
      )

      unseen <- neighbors[membership[neighbors] == 0L]

      if (length(unseen) > 0L) {
        membership[unseen] <- component_number
        queue[(last + 1L):(last + length(unseen))] <- unseen
        last <- last + length(unseen)
      }
    }
  }

  list(nc = component_number, comp.id = membership)
}

# ---------------------------------------------------------------------
# 2. Queen-contiguity weights, with disconnected fragments joined to the
#    LARGEST component (not an arbitrary one) by shortest centroid
#    distance. Every added link is logged so the report can disclose how
#    many artificial neighbour relationships the weights matrix contains
#    and how far they span.
# ---------------------------------------------------------------------
build_weights <- function(x) {
  nb <- poly2nb(
    x,
    queen = TRUE,
    row.names = as.character(x$GEOID20),
    snap = 50
  )

  # Some spdep versions encode an isolate as 0L rather than integer(0).
  n <- length(nb)
  for (i in seq_len(n)) {
    neighbors <- as.integer(nb[[i]])
    neighbors <- neighbors[!is.na(neighbors) & neighbors >= 1L & neighbors <= n]
    nb[[i]] <- sort(unique(neighbors))
  }
  attr(nb, "sym") <- TRUE

  centers <- suppressWarnings(st_centroid(st_geometry(x)))
  added <- list()

  repeat {
    components <- nb_components(nb)
    if (components$nc <= 1L) break

    sizes  <- tabulate(components$comp.id, nbins = components$nc)
    anchor <- which.max(sizes)

    inside  <- which(components$comp.id != anchor)
    outside <- which(components$comp.id == anchor)

    distance_matrix <- st_distance(centers[inside], centers[outside])
    closest <- arrayInd(which.min(distance_matrix), dim(distance_matrix))
    i <- inside[closest[1]]
    j <- outside[closest[2]]

    added[[length(added) + 1L]] <- tibble(
      from_geoid  = as.character(x$GEOID20[i]),
      to_geoid    = as.character(x$GEOID20[j]),
      distance_km = as.numeric(min(distance_matrix)) / 1000
    )

    nb[[i]] <- sort(unique(as.integer(c(nb[[i]], j))))
    nb[[j]] <- sort(unique(as.integer(c(nb[[j]], i))))
    attr(nb, "sym") <- TRUE
  }

  stopifnot(nb_components(nb)$nc == 1L, is.symmetric.nb(nb))

  listw <- nb2listw(nb, style = "W", zero.policy = TRUE)

  if (any(card(nb) < 1L)) stop("A zero-neighbour tract remained after linking.")
  weight_values <- unlist(listw$weights, use.names = FALSE)
  if (length(weight_values) == 0L || any(!is.finite(weight_values))) {
    stop("Non-finite spatial weights were created.")
  }

  list(
    nb    = nb,
    listw = listw,
    added = if (length(added) == 0L) tibble(
      from_geoid = character(), to_geoid = character(), distance_km = numeric()
    ) else bind_rows(added)
  )
}

# ---------------------------------------------------------------------
# 3. LISA classification. Returns a factor with the five standard
#    categories, using whichever vector of p-values is supplied (raw or
#    FDR-adjusted).
# ---------------------------------------------------------------------
lisa_classify <- function(values, listw, p_values, alpha = .05) {
  z      <- as.numeric(scale(values))
  lag_z  <- lag.listw(listw, z, zero.policy = TRUE)

  out <- case_when(
    p_values < alpha & z > 0 & lag_z > 0 ~ "High–High",
    p_values < alpha & z < 0 & lag_z < 0 ~ "Low–Low",
    p_values < alpha & z > 0 & lag_z < 0 ~ "High–Low",
    p_values < alpha & z < 0 & lag_z > 0 ~ "Low–High",
    TRUE ~ "Not significant"
  )

  factor(
    out,
    levels = c("High–High", "Low–Low", "High–Low", "Low–High", "Not significant")
  )
}

# ---------------------------------------------------------------------
# 4. Coefficient extraction that works for both lm and sarlm objects, so
#    OLS and SEM results can be compared with one code path.
# ---------------------------------------------------------------------
coef_table <- function(model) {
  m <- if (inherits(model, "Sarlm") || inherits(model, "sarlm")) {
    summary(model)$Coef
  } else {
    coef(summary(model))
  }
  tibble(
    Term     = rownames(m),
    Estimate = as.numeric(m[, 1]),
    SE       = as.numeric(m[, 2]),
    p        = as.numeric(m[, 4])
  )
}

# The strongest process is the largest ABSOLUTE coefficient among terms
# that are statistically distinguishable from zero. A large but
# non-significant estimate should never be reported as "most influential".
strongest_term <- function(model,
                           terms = c("redlining_z", "home_value_z", "environment_z"),
                           alpha = .05) {
  tab <- coef_table(model) |> filter(Term %in% terms, p < alpha)
  if (nrow(tab) == 0) return(NA_character_)
  tab$Term[which.max(abs(tab$Estimate))]
}

# Plain-language description of a single coefficient that respects BOTH
# sign and significance.
describe_coef <- function(model, term, alpha = .05) {
  tab <- coef_table(model)
  row <- tab[tab$Term == term, ]
  if (nrow(row) == 0) return("not estimated")
  if (row$p >= alpha) {
    return(paste0(
      "statistically indistinguishable from zero (b = ", fmt_num(row$Estimate),
      ", p ", fmt_p(row$p), ")"
    ))
  }
  paste0(
    "significantly ", if (row$Estimate > 0) "positive" else "negative",
    " (b = ", fmt_num(row$Estimate), ", p ", fmt_p(row$p), ")"
  )
}

# ---------------------------------------------------------------------
# 5. Variance inflation factors, computed directly so the report does not
#    depend on the car package.
# ---------------------------------------------------------------------
compute_vif <- function(model) {
  preds <- attr(terms(model), "term.labels")
  mf <- model.frame(model)
  vals <- vapply(preds, function(p) {
    others <- setdiff(preds, p)
    r2 <- summary(lm(reformulate(others, response = p), data = mf))$r.squared
    1 / (1 - r2)
  }, numeric(1))
  tibble(Term = preds, VIF = unname(vals))
}

# ---------------------------------------------------------------------
# 6. Lagrange multiplier / Rao score tests for choosing between a spatial
#    error and a spatial lag specification. spdep renamed these functions,
#    so both names are attempted.
# ---------------------------------------------------------------------
run_spatial_tests <- function(ols, listw) {
  res <- tryCatch(
    spdep::lm.RStests(
      ols, listw,
      test = c("RSerr", "RSlag", "adjRSerr", "adjRSlag"),
      zero.policy = TRUE
    ),
    error = function(e) {
      spdep::lm.LMtests(
        ols, listw,
        test = c("LMerr", "LMlag", "RLMerr", "RLMlag"),
        zero.policy = TRUE
      )
    }
  )

  tibble(
    Test        = names(res),
    Statistic   = vapply(res, function(z) unname(z$statistic), numeric(1)),
    `p value`   = vapply(res, function(z) unname(z$p.value), numeric(1))
  )
}

# ---------------------------------------------------------------------
# 7. Main pipeline. Two analytic samples are built on purpose:
#      * ice_data   - every tract with a valid ICE value. Used for the
#                     descriptive table, the ICE map, Moran's I and LISA,
#                     none of which require the predictors.
#      * model_data - complete cases on all four variables. Used only for
#                     the regressions.
#    Filtering everything down to the regression sample would needlessly
#    discard tracts from the descriptive and autocorrelation analysis and
#    would put artificial holes in the contiguity graph.
# ---------------------------------------------------------------------
prepare_city <- function(source, metro_value, short_name) {
  model_vars <- c(
    "index_concentration_extremes", "redlining",
    "median_home_value", "env_deprivation_index"
  )

  metro <- source |>
    filter(METRO_NAME == metro_value) |>
    st_make_valid() |>
    st_transform(5070)

  if (nrow(metro) == 0) stop(paste("No tracts matched", metro_value))

  missing_summary <- tibble(
    City     = short_name,
    Variable = model_vars,
    Missing  = vapply(
      model_vars,
      function(v) sum(!is.finite(metro[[v]])),
      integer(1)
    )
  ) |>
    mutate(`Percent missing` = 100 * Missing / nrow(metro))

  metro <- metro |>
    mutate(
      in_ice_sample   = is.finite(index_concentration_extremes),
      in_model_sample = if_all(all_of(model_vars), is.finite)
    )

  ice_data <- metro |> filter(in_ice_sample)

  model_data <- metro |>
    filter(in_model_sample) |>
    mutate(
      city          = short_name,
      redlining_z   = as.numeric(scale(redlining)),
      home_value_z  = as.numeric(scale(median_home_value)),
      environment_z = as.numeric(scale(env_deprivation_index))
    )

  if (nrow(model_data) < 30) {
    stop(paste("Too few complete tracts to model", short_name))
  }

  ice_w   <- build_weights(ice_data)
  model_w <- build_weights(model_data)

  # --- Global and local autocorrelation in ICE -------------------------
  set.seed(12345)
  moran_ice <- moran.mc(
    ice_data$index_concentration_extremes,
    ice_w$listw, nsim = 999, zero.policy = TRUE
  )

  set.seed(12345)
  local_i <- localmoran_perm(
    ice_data$index_concentration_extremes,
    ice_w$listw, nsim = 999, zero.policy = TRUE,
    alternative = "two.sided", iseed = 12345
  )

  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 (length(p_name) == 0L || is.na(p_name)) {
    stop("The LISA permutation p-value column was not found.")
  }

  local_p     <- local_i[, p_name]
  local_p_fdr <- p.adjust(local_p, method = "fdr")

  ice_data$lisa_cluster <- lisa_classify(
    ice_data$index_concentration_extremes, ice_w$listw, local_p
  )
  ice_data$lisa_cluster_fdr <- lisa_classify(
    ice_data$index_concentration_extremes, ice_w$listw, local_p_fdr
  )

  # --- Regression ------------------------------------------------------
  model_formula <- index_concentration_extremes ~
    redlining_z + home_value_z + environment_z

  ols <- lm(model_formula, data = st_drop_geometry(model_data))

  set.seed(12345)
  moran_ols <- moran.mc(
    residuals(ols), model_w$listw, nsim = 999, zero.policy = TRUE
  )

  spatial_tests <- run_spatial_tests(ols, model_w$listw)

  sem <- errorsarlm(
    model_formula,
    data = st_drop_geometry(model_data),
    listw = model_w$listw,
    zero.policy = TRUE,
    method = "Matrix",
    interval = c(-0.99, 0.99),
    tol.solve = 1e-12
  )

  set.seed(12345)
  moran_sem <- moran.mc(
    residuals(sem), model_w$listw, nsim = 999, zero.policy = TRUE
  )

  list(
    name             = short_name,
    metro            = metro,
    ice_data         = ice_data,
    model_data       = model_data,
    missing_summary  = missing_summary,
    ice_w            = ice_w,
    model_w          = model_w,
    moran_ice        = moran_ice,
    local_p          = local_p,
    local_p_fdr      = local_p_fdr,
    ols              = ols,
    vif              = compute_vif(ols),
    spatial_tests    = spatial_tests,
    moran_ols        = moran_ols,
    sem              = sem,
    moran_sem        = moran_sem
  )
}
san_diego <- prepare_city(data, metro_san_diego, "San Diego–Carlsbad")
minneapolis <- prepare_city(data, metro_minneapolis, "Minneapolis–St. Paul")

cities <- list(san_diego, minneapolis)
city_a <- san_diego$name
city_b <- minneapolis$name

# --- Coverage / missingness -------------------------------------------
coverage_table <- map_dfr(cities, function(r) {
  tibble(
    City                  = r$name,
    `Tracts in supplied file` = nrow(r$metro),
    `Valid ICE`           = nrow(r$ice_data),
    `Complete cases`      = nrow(r$model_data),
    `Percent complete`    = 100 * nrow(r$model_data) / nrow(r$metro)
  )
})

missing_table <- map_dfr(cities, ~ .x$missing_summary)

link_table <- map_dfr(cities, function(r) {
  tibble(
    City                   = r$name,
    `Sample`               = c("ICE sample", "Regression sample"),
    `Artificial links`     = c(nrow(r$ice_w$added), nrow(r$model_w$added)),
    `Longest link (km)`    = c(
      if (nrow(r$ice_w$added)) max(r$ice_w$added$distance_km) else NA_real_,
      if (nrow(r$model_w$added)) max(r$model_w$added$distance_km) else NA_real_
    )
  )
})

# --- Descriptives ------------------------------------------------------
desc_table <- map_dfr(cities, function(r) {
  y <- r$ice_data$index_concentration_extremes
  tibble(
    City                = r$name,
    Tracts              = length(y),
    Mean                = mean(y),
    SD                  = sd(y),
    Minimum             = min(y),
    `25th percentile`   = unname(quantile(y, .25)),
    Median              = median(y),
    `75th percentile`   = unname(quantile(y, .75)),
    Maximum             = max(y)
  )
})

# --- LISA counts (raw and FDR-adjusted) --------------------------------
lisa_table <- map_dfr(cities, function(r) {
  r$ice_data |>
    st_drop_geometry() |>
    count(lisa_cluster, .drop = FALSE, name = "Tracts") |>
    mutate(City = r$name, Percent = 100 * Tracts / sum(Tracts)) |>
    select(City, `LISA cluster` = lisa_cluster, Tracts, Percent)
})

lisa_robustness <- map_dfr(cities, function(r) {
  d <- st_drop_geometry(r$ice_data)
  tibble(
    City                 = r$name,
    `Significant (raw)`  = sum(d$lisa_cluster != "Not significant"),
    `Significant (FDR)`  = sum(d$lisa_cluster_fdr != "Not significant"),
    `Retained`           = 100 * sum(d$lisa_cluster_fdr != "Not significant") /
      max(1, sum(d$lisa_cluster != "Not significant"))
  )
})

# --- Model results in one tidy frame -----------------------------------
tidy_model <- function(r, which_model) {
  m <- if (which_model == "OLS") r$ols else r$sem
  coef_table(m) |>
    mutate(City = r$name, Model = which_model, .before = 1)
}

all_model_results <- bind_rows(
  map_dfr(cities, tidy_model, which_model = "OLS"),
  map_dfr(cities, tidy_model, which_model = "Spatial error")
)

ols_fit <- map_dfr(cities, function(r) {
  tibble(
    City                 = r$name,
    Tracts               = nrow(r$model_data),
    `Adjusted R-squared` = summary(r$ols)$adj.r.squared,
    AIC                  = AIC(r$ols)
  )
})

sem_fit <- map_dfr(cities, function(r) {
  tibble(
    City      = r$name,
    Lambda    = unname(r$sem$lambda),
    `SEM AIC` = AIC(r$sem),
    `OLS AIC` = AIC(r$ols)
  )
})

vif_table <- map_dfr(cities, function(r) mutate(r$vif, City = r$name, .before = 1))

correlation_table <- map_dfr(cities, function(r) {
  d <- r$model_data |>
    st_drop_geometry() |>
    select(
      ICE = index_concentration_extremes,
      redlining_z, home_value_z, environment_z
    )
  as.data.frame(cor(d)) |>
    rownames_to_column("Variable") |>
    mutate(City = r$name, .before = 1)
})

spatial_test_table <- map_dfr(cities, function(r) {
  mutate(r$spatial_tests, City = r$name, .before = 1)
})

residual_table <- map_dfr(cities, function(r) {
  tibble(
    City                  = r$name,
    Model                 = c("OLS", "Spatial error"),
    `Residual Moran's I`  = c(
      unname(r$moran_ols$statistic), unname(r$moran_sem$statistic)
    ),
    `Permutation p-value` = c(r$moran_ols$p.value, r$moran_sem$p.value)
  )
})

# --- Formal cross-city coefficient comparison --------------------------
# The two cities are separate, non-overlapping samples, so their estimates
# are independent and the standard error of the difference is simply the
# root of the summed variances. This is the test the nonstationarity claim
# actually requires; comparing point estimates by eye is not sufficient.
coefficient_difference <- function(model_name) {
  a <- all_model_results |>
    filter(Model == model_name, City == city_a, Term != "(Intercept)") |>
    arrange(Term)
  b <- all_model_results |>
    filter(Model == model_name, City == city_b, Term != "(Intercept)") |>
    arrange(Term)

  difference <- a$Estimate - b$Estimate
  se_diff    <- sqrt(a$SE^2 + b$SE^2)
  z          <- difference / se_diff

  out <- tibble(
    Model        = model_name,
    Term         = a$Term,
    A            = a$Estimate,
    B            = b$Estimate,
    Difference   = difference,
    `SE of difference` = se_diff,
    z            = z,
    `p value`    = 2 * pnorm(-abs(z))
  )
  names(out)[3:4] <- c(city_a, city_b)
  out
}

difference_table <- bind_rows(
  coefficient_difference("OLS"),
  coefficient_difference("Spatial error")
)

n_sig_differences <- sum(
  difference_table$`p value`[difference_table$Model == "Spatial error"] < .05
)

sd_strongest   <- strongest_term(san_diego$ols)
mpls_strongest <- strongest_term(minneapolis$ols)

name_or_none <- function(term) {
  if (is.na(term)) "no predictor reached significance" else unname(pretty_term[term])
}

Executive summary

Social segregation is strongly clustered in both metropolitan areas. In each metro, historical redlining is negatively associated with the Index of Concentration at the Extremes and median home value is positively associated with it, and in both cases substantial spatial structure remains in the OLS residuals, which a spatial error model absorbs. The most striking result is not how differently the two metros behave but how similarly: after standardising predictors within each city, the coefficients are close enough that 0 of three cross-city differences in the spatial error models are statistically distinguishable from zero. Spatial nonstationarity is not detectable in the three measured coefficients; the clearer cross-city difference is in the amount of unexplained spatial structure.

1. Research context and hypotheses

Why these cities?

San Diego–Carlsbad is in the assignment’s large coastal/high-density group, while Minneapolis–St. Paul is in the mid-sized/smaller metro group. The pairing is useful because both metros were shaped by discriminatory housing policy, but their housing markets, growth histories, climates, and regional forms differ substantially. San Diego is the primary city.

The outcome is the Index of Concentration at the Extremes (ICE), supplied in the course dataset at the census-tract level. It ranges from −1 to 1: values near 1 indicate concentration of high-income, non-Hispanic White households and values near −1 indicate concentration of low-income households of colour. A value near zero indicates a tract that is mixed on these dimensions, which is not the same thing as integration.

The standard race–income ICE construction contrasts non-Hispanic White households with income of at least $100,000 against households of colour with income of no more than $25,000, divided by all households in the tract (IPUMS CDOH, 2024). The supplied GeoJSON does not identify the ACS vintage used to create its precomputed ICE field, so the analysis does not assign an unsupported year to the measure.

San Diego: high housing costs layered over older divisions

San Diego’s own fair-housing assessment describes continuing regional and local patterns of segregation by race, ethnicity, disability, and income (City of San Diego, 2021). The historical geography is unusually visible. A 1936 HOLC map rated much of southeastern San Diego poorly while rating places such as La Jolla and Coronado more favourably; local housing experts argue that those lines still resemble present socioeconomic divisions, even though some neighbourhoods have changed through redevelopment and gentrification (KPBS, 2018).

Contemporary housing pressure reinforces this pattern. SANDAG reports substantial differences in overcrowding and housing conditions across racial and tenure groups, with renters more likely than owners to live in overcrowded housing (SANDAG, 2023). Home value is therefore more than a market indicator: it also captures who can gain access to high-opportunity neighbourhoods in a very expensive region.

Environmental inequality provides the third process. The City of San Diego’s Environmental Justice Element recognises that some communities carry disproportionate health risks from development patterns, pollution exposure, transportation infrastructure, and limited access to healthy public space (City of San Diego, 2024). Statewide research using CalEnviroScreen likewise found that communities of colour face greater cumulative environmental health burdens (Cushing et al., 2015).

Minneapolis–St. Paul: covenants, redlining, and uneven environmental privilege

Minneapolis–St. Paul’s history cannot be explained by HOLC redlining alone. The University of Minnesota’s Mapping Prejudice project documented racial covenants that barred non-White residents from owning or occupying property, and its maps connect those covenants to present neighbourhoods of colour and later freeway construction (Mapping Prejudice, n.d.). The Metropolitan Council similarly identifies redlining as a barrier to homeownership and intergenerational wealth whose economic, environmental, and health effects remain visible across the region (Metropolitan Council, n.d.).

Walker, Keeler, and Derickson found that areas protected by racial covenants have more tree cover and lower temperatures today, linking present environmental privilege to historical housing discrimination (Walker et al., 2024). Minnesota climate officials also report a nearly 11 °F surface-temperature gap between formerly redlined and non-redlined parts of Minneapolis in one multi-city study (Minnesota State Climatology Office, n.d.).

Hypotheses

The assignment asks for hypotheses about relative strength across the two cities, not only about direction. The table below states each prediction before the models are estimated, and Section 5 returns to it.

Pre-registered expectations, stated before estimation
Term Expected direction Expected strength (San Diego) Expected cross-city comparison
Historical redlining (z) Negative Moderate to strong Weaker in San Diego than in Minneapolis–St. Paul: postwar growth, redevelopment and gentrification have partly overwritten the small 1930s HOLC footprint, whereas the Twin Cities’ graded core remains closer to its historical form.
Median home value (z) Positive Strongest of the three Stronger in San Diego than in Minneapolis–St. Paul: extreme and sustained price pressure sorts households by income and wealth more sharply than in a lower-cost market.
Environmental deprivation (z) Negative Moderate Weaker in San Diego than in Minneapolis–St. Paul: the Twin Cities literature ties environmental amenity directly to covenant geography, whereas San Diego’s burdens are concentrated in a few corridors.

Interpretation note. Redlining values run from about 1 (mostly A-graded) to 4 (mostly D-graded), so a negative coefficient means that historically downgraded tracts have lower ICE today. A higher environmental-deprivation score means worse conditions, so that coefficient is also expected to be negative. The expected home-value coefficient is positive.

2. Data coverage and analytic samples

The supplied dataset contains the HOLC-linked tract footprint rather than a complete enumeration of either metropolitan area. Every spatial statistic below therefore describes the tracts supplied for that footprint, not the full modern MSA.

Analytic sample sizes within the tract footprint supplied for each metro
City Tracts in supplied file Valid ICE Complete cases Percent complete
San Diego–Carlsbad 131 131 131 100.0%
Minneapolis–St. Paul 188 188 188 100.0%
Missingness by variable within each metro
City Variable Missing Percent missing
San Diego–Carlsbad index_concentration_extremes 0 0.0%
San Diego–Carlsbad redlining 0 0.0%
San Diego–Carlsbad median_home_value 0 0.0%
San Diego–Carlsbad env_deprivation_index 0 0.0%
Minneapolis–St. Paul index_concentration_extremes 0 0.0%
Minneapolis–St. Paul redlining 0 0.0%
Minneapolis–St. Paul median_home_value 0 0.0%
Minneapolis–St. Paul env_deprivation_index 0 0.0%
Figure 1. Tracts in the supplied HOLC-linked footprint. Dark tracts are complete regression cases; any light tracts are excluded because at least one analysis field is missing.

Figure 1. Tracts in the supplied HOLC-linked footprint. Dark tracts are complete regression cases; any light tracts are excluded because at least one analysis field is missing.

Of the 131 San Diego tracts supplied, 131 have a valid ICE value and 131 are complete on all four analysis variables. The corresponding Minneapolis–St. Paul figures are 188, 188, and 188. These counts show whether any observations are lost within the supplied file; they should not be mistaken for the total number of census tracts in either MSA. The source footprint is limited to tracts linked to the areas HOLC mapped in the 1930s, which largely represent the pre-war urbanised core.

This has three consequences that are carried through the rest of the report:

  1. The supplied footprint represents the historically graded core, not the metro. Statements about suburban rings or outlying growth areas cannot be supported by these models, and are avoided below.
  2. The contiguity graph may contain disconnected fragments. Ungraded areas outside the supplied footprint do not enter the graph, so the analysis does not represent adjacency across the complete contemporary metropolitan area.
  3. Disconnected fragments must be joined artificially. The table below records how many such links were added and how far the longest one spans.
Artificial neighbour links required to produce a single connected weights graph
City Sample Artificial links Longest link (km)
San Diego–Carlsbad ICE sample 5 2.9
San Diego–Carlsbad Regression sample 5 2.9
Minneapolis–St. Paul ICE sample 2 4.3
Minneapolis–St. Paul Regression sample 2 4.3

Because the predictors are not needed for the descriptive, Moran’s I, or LISA analyses, those use every tract with valid ICE; only the regressions require complete cases on all four fields. If the missingness table reports zero missing values, the two samples are identical.

3. Describing social segregation

Descriptive statistics

Distribution of the Index of Concentration at the Extremes (ICE sample)
City Tracts Mean SD Minimum 25th percentile Median 75th percentile Maximum
San Diego–Carlsbad 131 0.193 0.252 -0.297 -0.009 0.193 0.384 0.699
Minneapolis–St. Paul 188 0.199 0.241 -0.425 0.025 0.183 0.385 0.721

The two distributions are remarkably alike. The means differ by 0.006 and the standard deviations by 0.011 — differences far smaller than the within-city spread. Both distributions sit well above zero (means of 0.193 and 0.199), which reflects the fact that ICE is measured over tracts rather than over people and that both metros contain many predominantly White, moderate-to-high-income tracts.

Two features deserve emphasis. First, neither metro contains a tract anywhere near the theoretical minimum of −1: the lowest observed values are -0.297 and -0.425, with the lower of the two in Minneapolis–St. Paul. “Concentrated disadvantage” in this analysis is therefore relative rather than absolute — no tract in either sample is close to being entirely low-income and non-White on this measure, while the maxima (0.699 and 0.721) approach the concentrated-advantage end far more closely. The index is asymmetric in practice even though it is symmetric by construction. Second, the interquartile ranges are 0.393 and 0.359, with the wider spread in San Diego–Carlsbad.

The substantive implication is the motivation for everything that follows: if the two metros have nearly identical distributions of ICE, then whatever distinguishes them must be spatial — how those values are arranged, not which values occur. A non-spatial summary cannot detect that difference, which is exactly what Moran’s I and LISA are for.

ICE maps

Figure 2. ICE by census tract, on a shared colour scale. Any light grey tracts lack valid ICE. The maps show only the tract footprints supplied for the two metros.

Figure 2. ICE by census tract, on a shared colour scale. Any light grey tracts lack valid ICE. The maps show only the tract footprints supplied for the two metros.

Within the mapped footprint, San Diego shows a clear divide between areas of concentrated advantage and concentrated disadvantage rather than a random mix of tract values. Higher-ICE tracts are concentrated farther west, closer to the coast, while many of the lowest-ICE tracts form a compact group farther inland and to the southeast. This description is limited to the shaded HOLC-linked footprint and does not describe North County as a whole. The broad orientation is consistent with the 1936 HOLC geography and the city’s current fair-housing findings (City of San Diego, 2021), though it is not a reproduction of the 1930s map: redevelopment and gentrification have changed some neighbourhoods.

Minneapolis–St. Paul also shows clear spatial sorting within its mapped core: similar ICE values form adjoining groups instead of alternating randomly from tract to tract. The shaded footprint does not cover the full suburban ring, so the map supports a comparison of clusters within the historically mapped core, not a claim about the entire modern metropolitan area. Mapping Prejudice helps explain why the observed pattern cannot be attributed to HOLC grades alone: racial covenants created a wider geography of exclusion, reserving land, housing wealth, and environmental amenities for White households.

Early comparison. The distributions are nearly identical (Section 3), so the interesting difference is arrangement. San Diego’s pattern is stretched along a coastal axis, while Minneapolis–St. Paul’s clustering is contained within a different, twin-city historical footprint. The same ICE value therefore sits inside a different urban history in each metro.

4. Global and local spatial structure

Global Moran’s I

Queen contiguity defines tracts as neighbours when they share an edge or a corner, and the weights are row-standardised. Moran’s I is evaluated against 999 random permutations, so the smallest attainable p-value is 0.001; a reported value of 0.001 means “at the permutation floor”, not an exact probability.

Global spatial autocorrelation in ICE
City Tracts Moran’s I Permutation p-value
San Diego–Carlsbad 131 0.765 0.001 (floor)
Minneapolis–St. Paul 188 0.585 0.001 (floor)

San Diego’s Moran’s I is 0.765 (p = .001) and Minneapolis–St. Paul’s is 0.585 (p = .001). Both reject the null hypothesis of spatial randomness decisively. The larger statistic belongs to San Diego–Carlsbad, indicating stronger overall clustering under the same neighbour definition.

That difference should be read with the coverage results in mind. Both graphs contain gaps created by missing tracts, and gaps attenuate measured autocorrelation, so these values are conservative estimates of the true clustering in each metro rather than exact quantities.

Local Indicators of Spatial Association (LISA)

Figure 3. Local ICE clusters from 999 permutations at p < .05, unadjusted for multiple testing. Light grey tracts are outside the ICE sample.

Figure 3. Local ICE clusters from 999 permutations at p < .05, unadjusted for multiple testing. Light grey tracts are outside the ICE sample.

LISA classification by metropolitan area (unadjusted p < .05)
City LISA cluster Tracts Percent
San Diego–Carlsbad High–High 22 16.8%
San Diego–Carlsbad Low–Low 29 22.1%
San Diego–Carlsbad High–Low 0 0.0%
San Diego–Carlsbad Low–High 1 0.8%
San Diego–Carlsbad Not significant 79 60.3%
Minneapolis–St. Paul High–High 28 14.9%
Minneapolis–St. Paul Low–Low 22 11.7%
Minneapolis–St. Paul High–Low 0 0.0%
Minneapolis–St. Paul Low–High 1 0.5%
Minneapolis–St. Paul Not significant 137 72.9%
Robustness of LISA results to false-discovery-rate correction
City Significant (raw) Significant (FDR) Retained
San Diego–Carlsbad 52 33 63.5%
Minneapolis–St. Paul 51 25 49.0%

San Diego has 22 High–High tracts and 29 Low–Low tracts; Minneapolis–St. Paul has 28 and 22. In both metros the clusters are large, contiguous blocks rather than scattered individual tracts, which is the local expression of the global statistics above.

Spatial outliers are almost absent: 1 in San Diego and 1 in Minneapolis–St. Paul. That absence is substantively informative. A High–Low tract would be an enclave of concentrated advantage embedded in disadvantage, and a Low–High tract the reverse. Their near-total absence means segregation in both metros operates at a scale larger than the individual tract: advantage and disadvantage come in districts, not islands. Policy that treats a single tract as the unit of intervention is working below the scale at which the pattern exists.

Because LISA runs one hypothesis test per tract, roughly five percent of tracts would be flagged by chance alone at α = .05. The robustness table above re-runs the classification using false-discovery-rate-adjusted p-values. The share of significant tracts that survive adjustment is 63.5% in San Diego and 49.0% in Minneapolis–St. Paul, so the cluster cores reported above are not artefacts of multiple testing.

5. Modelling the processes

Specification and diagnostics

For each city the outcome is ICE and the predictors are historical redlining, median home value, and environmental deprivation, each standardised within that city. A coefficient is therefore the expected change in ICE associated with a one-standard-deviation increase in that predictor within that metro, holding the other two constant.

# This is the model estimated inside prepare_city() above; it is shown here
# only for reference.
index_concentration_extremes ~ redlining_z + home_value_z + environment_z

Before interpreting coefficients, the predictors need to be checked for collinearity: if they overlap heavily, the partial coefficients will be unstable and “which process matters most” becomes unanswerable.

Correlations among ICE and the three standardised predictors, by city
City Variable ICE redlining_z home_value_z environment_z
San Diego–Carlsbad ICE 1.00 -0.72 0.78 -0.18
San Diego–Carlsbad Redlining (z) -0.72 1.00 -0.72 0.33
San Diego–Carlsbad Home value (z) 0.78 -0.72 1.00 -0.12
San Diego–Carlsbad Environmental deprivation (z) -0.18 0.33 -0.12 1.00
Minneapolis–St. Paul ICE 1.00 -0.63 0.66 -0.15
Minneapolis–St. Paul Redlining (z) -0.63 1.00 -0.52 0.33
Minneapolis–St. Paul Home value (z) 0.66 -0.52 1.00 -0.25
Minneapolis–St. Paul Environmental deprivation (z) -0.15 0.33 -0.25 1.00
Variance inflation factors
City Term VIF
San Diego–Carlsbad Historical redlining (z) 2.39
San Diego–Carlsbad Median home value (z) 2.17
San Diego–Carlsbad Environmental deprivation (z) 1.16
Minneapolis–St. Paul Historical redlining (z) 1.45
Minneapolis–St. Paul Median home value (z) 1.39
Minneapolis–St. Paul Environmental deprivation (z) 1.13

The largest variance inflation factor is 2.39, so the coefficients below are not destabilised by collinearity. The correlation table matters for a second reason: it reports each predictor’s bivariate relationship with ICE alongside the multivariate coefficient. Environmental deprivation correlates with ICE at -0.18 in San Diego and -0.15 in Minneapolis–St. Paul on its own. Comparing those figures with the adjusted coefficients shows how much of the environmental association is shared with redlining and home value rather than independent of them.

Choosing between a spatial error and a spatial lag specification

Lagrange multiplier / Rao score tests on the OLS residuals
City Test Statistic p value
San Diego–Carlsbad RSerr 63.28 < .001
San Diego–Carlsbad RSlag 60.10 < .001
San Diego–Carlsbad adjRSerr 9.76 0.002
San Diego–Carlsbad adjRSlag 6.58 0.010
Minneapolis–St. Paul RSerr 36.59 < .001
Minneapolis–St. Paul RSlag 44.52 < .001
Minneapolis–St. Paul adjRSerr 1.44 0.231
Minneapolis–St. Paul adjRSlag 9.36 0.002

These tests distinguish an error process (spatially clustered omitted variables) from a lag process (ICE in one tract directly influencing ICE in its neighbours). The robust versions are the ones to read when both simple tests are significant, since each simple test has power against the other alternative. The substantive argument also favours the error specification here: zoning, school-district boundaries, transit investment, racial covenants, and neighbourhood reputation are all spatially clustered influences that the three predictors do not measure, and there is no strong theoretical reason to think a tract’s ICE mechanically causes its neighbour’s.

San Diego–Carlsbad: Both robust tests are significant, so the diagnostics do not uniquely distinguish a lag from an error process. The SEM is used because the omitted influences identified above are spatially clustered and there is no strong mechanism by which one tract’s ICE directly causes its neighbour’s ICE.

Minneapolis–St. Paul: The robust lag test is significant while the robust error test is not. The diagnostics favour a lag model; the SEM below is therefore interpreted as an assignment-directed sensitivity model rather than a uniquely preferred specification.

OLS results

OLS models of ICE using city-standardised predictors
City Model Term Estimate Standard error p value
San Diego–Carlsbad OLS Intercept 0.193 0.013 < .001
San Diego–Carlsbad OLS Historical redlining (z) -0.083 0.020 < .001
San Diego–Carlsbad OLS Median home value (z) 0.136 0.019 < .001
San Diego–Carlsbad OLS Environmental deprivation (z) -0.002 0.014 0.881
Minneapolis–St. Paul OLS Intercept 0.199 0.012 < .001
Minneapolis–St. Paul OLS Historical redlining (z) -0.102 0.014 < .001
Minneapolis–St. Paul OLS Median home value (z) 0.114 0.014 < .001
Minneapolis–St. Paul OLS Environmental deprivation (z) 0.025 0.012 0.043
OLS model fit
City Tracts Adjusted R-squared AIC
San Diego–Carlsbad 131.000 0.651 -121.613
Minneapolis–St. Paul 188.000 0.557 -148.558

San Diego

Holding the other predictors constant, redlining is significantly negative (b = -0.083, p < .001), median home value is significantly positive (b = 0.136, p < .001), and environmental deprivation is statistically indistinguishable from zero (b = -0.002, p = .881). Among the predictors that reach significance, the largest absolute coefficient belongs to Median home value (z).

This supports the primary San Diego hypothesis: the contemporary housing market is the strongest of the three measured processes, consistent with a region where price pressure sorts households by income and wealth. The environmental hypothesis is not supported. The environmental-deprivation coefficient is -0.002 with a standard error of 0.014 and p = .881 — statistically indistinguishable from zero, not merely small. A moderate negative relationship was predicted and was not found. Two readings are consistent with the literature cited earlier. The composite index averages ozone, particulate matter, diesel, air-toxics risk, and proximity to rail and airports across a whole tract, which can wash out the corridor-scale burdens that San Diego’s Environmental Justice Element identifies. Alternatively, the environmental burden may be almost entirely shared with redlining and home value, in which case it has little left to explain once they are in the model — the bivariate correlation reported above distinguishes these cases.

These are conditional associations in cross-sectional data, not causal effects.

Minneapolis–St. Paul

Redlining is significantly negative (b = -0.102, p < .001), median home value is significantly positive (b = 0.114, p < .001), and environmental deprivation is significantly positive (b = 0.025, p = .043). The largest significant absolute coefficient belongs to Median home value (z).

The redlining coefficient is -0.102 (p < .001). The HOLC measure therefore retains a clear independent association with ICE, even though Mapping Prejudice shows that covenants operated across a wider area than the graded neighbourhoods. The measured effect should be read as a lower bound on the influence of discriminatory housing policy in this metro, since the covenant geography is not in the model at all. The environmental result is the most surprising finding in the report. The coefficient is 0.025 with p = .043 — statistically significant and positive, the opposite of the predicted sign. Taken at face value it says that, among tracts alike in redlining grade and home value, greater environmental deprivation is associated with higher ICE. The most likely explanation is compositional rather than substantive. The composite index is weighted toward proximity to transport infrastructure and air-toxics exposure, both of which are elevated near the region’s freeway and rail corridors and near the airport — locations that also include established, higher-income, predominantly White areas. Once redlining and home value absorb the disadvantage gradient, what remains of the environmental index may be tracking infrastructure adjacency rather than deprivation. The Walker et al. finding concerns tree canopy and surface temperature specifically, which this composite does not include, so the result does not contradict that literature so much as fail to measure the same thing.

Residual spatial autocorrelation

Spatial autocorrelation in residuals before and after the spatial error model
City Model Residual Moran’s I Permutation p-value
San Diego–Carlsbad OLS 0.464 0.001 (floor)
San Diego–Carlsbad Spatial error -0.025 0.599
Minneapolis–St. Paul OLS 0.282 0.001 (floor)
Minneapolis–St. Paul Spatial error -0.012 0.544

San Diego’s OLS residual Moran’s I is 0.464 (p = .001) and Minneapolis–St. Paul’s is 0.282 (p = .001). Both models leave substantial clustered structure unexplained, which violates the OLS independence assumption and means the OLS standard errors are too small. A spatial error model is required, not optional.

The two values are not equal, and the gap is itself a finding. San Diego retains the larger residual statistic (0.464 against 0.282), meaning that the same three predictors leave more geographically organised variation unaccounted for in San Diego than in Minneapolis–St. Paul.

Spatial error models

Spatial error models of ICE
City Model Term Estimate Standard error p value
San Diego–Carlsbad Spatial error Intercept 0.199 0.030 < .001
San Diego–Carlsbad Spatial error Historical redlining (z) -0.072 0.020 < .001
San Diego–Carlsbad Spatial error Median home value (z) 0.099 0.021 < .001
San Diego–Carlsbad Spatial error Environmental deprivation (z) -0.009 0.020 0.664
Minneapolis–St. Paul Spatial error Intercept 0.200 0.023 < .001
Minneapolis–St. Paul Spatial error Historical redlining (z) -0.087 0.015 < .001
Minneapolis–St. Paul Spatial error Median home value (z) 0.095 0.016 < .001
Minneapolis–St. Paul Spatial error Environmental deprivation (z) 0.003 0.016 0.879
Spatial error dependence and model fit
City Lambda SEM AIC OLS AIC
San Diego–Carlsbad 0.669 -171.300 -121.613
Minneapolis–St. Paul 0.558 -180.375 -148.558

The spatial-error parameter lambda is 0.669 in San Diego and 0.558 in Minneapolis–St. Paul, and AIC falls in both cities relative to OLS. After the SEM, residual Moran’s I drops to -0.025 and -0.012 respectively.

Neither remaining residual pattern is statistically significant, so the SEM adequately absorbs the detected spatial structure in both cities.

Figure 4. Standardised coefficients with 95% intervals, OLS versus spatial error, by city. Overlapping intervals across the two panels indicate that a coefficient does not differ detectably between metros.

Figure 4. Standardised coefficients with 95% intervals, OLS versus spatial error, by city. Overlapping intervals across the two panels indicate that a coefficient does not differ detectably between metros.

Comparing the two model families within each city: coefficients that shrink after spatial adjustment were partly absorbing broader geographic context under OLS, while coefficients that hold steady are robust to the modelled dependence. In both cities the OLS standard errors are the smaller ones, which is the expected consequence of ignoring positive residual autocorrelation — a further reason to report the SEM estimates as the headline results.

6. Cross-city comparison and spatial nonstationarity

Spatial nonstationarity means that a process does not hold the same relationship with the outcome everywhere. Establishing it requires showing that coefficients differ by more than sampling error, not simply that the point estimates are unequal. Because the two cities are separate, non-overlapping samples, their estimates are independent and the difference can be tested directly.

Formal test of cross-city differences in standardised coefficients
Model Term San Diego–Carlsbad Minneapolis–St. Paul Difference SE of difference z p value
OLS Environmental deprivation (z) -0.002 0.025 -0.028 0.019 -1.469 0.142
OLS Median home value (z) 0.136 0.114 0.022 0.024 0.944 0.345
OLS Historical redlining (z) -0.083 -0.102 0.019 0.025 0.768 0.442
Spatial error Environmental deprivation (z) -0.009 0.003 -0.011 0.026 -0.433 0.665
Spatial error Median home value (z) 0.099 0.095 0.005 0.026 0.177 0.860
Spatial error Historical redlining (z) -0.072 -0.087 0.014 0.025 0.568 0.570

In the spatial error models, 0 of 3 cross-city coefficient differences are statistically distinguishable from zero at p < .05. This is the central result of the comparison, and it is not the one the hypotheses predicted. Despite different climates, densities, growth histories, and housing markets, none of the three standardised relationships with ICE differs detectably between these samples. This is a failure to find coefficient nonstationarity, not proof that the processes are identical.

That finding is still substantively meaningful. Within the uncertainty of this two-city analysis, the associations between historical redlining and ICE and between housing values and ICE are similar across two different urban contexts. A one-standard-deviation worse redlining grade is associated with about the same reduction in ICE in both samples, although the design has limited power to rule out smaller differences.

Where the two metros do differ is in what the model fails to explain:

Where the two metros genuinely differ: spatial structure, not coefficients
City Moran’s I (ICE) Residual Moran’s I (OLS) Lambda (SEM) Adjusted R-squared (OLS)
San Diego–Carlsbad 0.765 0.464 0.669 0.651
Minneapolis–St. Paul 0.585 0.282 0.558 0.557
  1. Unexplained spatial structure is the real difference. The same three predictors leave residual Moran’s I of 0.464 in San Diego against 0.282 in Minneapolis–St. Paul, with lambda of 0.669 against 0.558. Whatever organises segregation beyond redlining, home values, and environmental burden is more strongly clustered in San Diego. Physical geography is a plausible candidate: the coast, canyon topography, military installations, and the international border impose large-scale constraints on the San Diego housing market that have no Twin Cities equivalent.

  2. Historical policy is expressed through different institutions even when the coefficient is similar. San Diego’s HOLC geography maps onto a coastal axis; in Minneapolis–St. Paul the HOLC variable measures only part of the historical process, because racial covenants created exclusion across a wider area than the graded neighbourhoods. Similar coefficients on a variable that means somewhat different things in each place is itself worth stating carefully.

  3. The environmental process behaves differently, but the difference is not significant. The environmental coefficients differ in sign across the two metros in the OLS models, which is the sort of result that invites a nonstationarity story — but the formal test above shows the difference is within sampling error, and both estimates are close to zero. The honest conclusion is that the composite index does not capture the environmental dimension the literature identifies in either city.

Hypotheses against results

Spatial error estimates against the pre-stated expectations (* = p < .05)
Term San Diego (SEM) Minneapolis (SEM) Difference detectable?
Historical redlining (z) -0.072* -0.087* No
Median home value (z) 0.099* 0.095* No
Environmental deprivation (z) -0.009 (n.s.) 0.003 (n.s.) No

Read against the hypothesis table in Section 1: the directional predictions for redlining and home value hold in both metros, while the environmental prediction fails in both. Every prediction about relative strength across the two cities is rejected, because none of the cross-city differences can be distinguished from sampling error.

Bottom line. Both metros contain geographically concentrated advantage and disadvantage, and neither can be explained as the aggregate of household preference. Historical housing discrimination shaped where investment and homeownership were possible; contemporary housing markets reproduce those advantages. What this comparison adds is a caution about how quickly nonstationarity is invoked: unequal point estimates are not evidence of different processes, and here the measured coefficients do not differ detectably across the two regions. The clearest geographic difference lies in the unexplained residual structure, which is what should be investigated next.

7. Limitations

The most important limitation is coverage. The supplied source footprint is restricted to tracts linked to a HOLC redlining score, which means the historically graded urban core rather than the metropolitan area. Conclusions do not generalise to suburban or exurban tracts, and omitted areas plus the artificial connecting links may change the autocorrelation statistics by an unknown amount. Results also depend on the queen-contiguity definition and on the artificial links used to connect disconnected fragments, both documented in Section 2.

This is a cross-sectional, observational analysis, so coefficients are not causal effects. ICE compresses race, ethnicity, and income into one index and cannot describe every dimension of segregation. The HOLC variable omits racial covenants and exclusionary zoning, which matters especially in Minneapolis–St. Paul. Median home value is plausibly an outcome of segregation as well as a process reinforcing it, so its coefficient should not be read as one-directional. The environmental composite does not include tree canopy or surface temperature, the dimensions most directly implicated by the Twin Cities literature. LISA involves many local tests and is exploratory even after FDR adjustment. Finally, tracts vary in population and area, and their boundaries need not correspond to residents’ own sense of neighbourhood.

A two-city comparison also has limited power to detect nonstationarity. Failing to reject equality is not proof of equality, and the null differences reported in Section 6 are consistent with modest true differences that this design cannot resolve.

8. Conclusion

San Diego and Minneapolis–St. Paul both show strong geographic sorting on ICE. Their ICE distributions are nearly identical, so the interesting variation is spatial rather than distributional; Moran’s I and LISA confirm that advantage and disadvantage occur in large contiguous districts rather than isolated tracts in either metro. The regression models show that historical redlining is negatively associated with ICE and median home value positively associated with it in both cities, that substantial spatial structure remains in the OLS residuals, and that a spatial error model absorbs it.

The cross-city comparison produced the opposite of the expected result. Tested formally, the standardised coefficients are not distinguishable between the two metros. The analysis therefore finds no detectable coefficient nonstationarity and locates the clearest geographic difference in the unexplained residual structure, which is larger in San Diego. The methodological lesson is that nonstationarity is a testable claim rather than an interpretive default, and that comparing coefficients without comparing their uncertainty will manufacture differences that the data do not contain.

References

City of San Diego. (2021). Assessment of fair housing. https://www.sandiego.gov/sites/default/files/he_appa_assessmentfairhousing_final.pdf

City of San Diego. (2024). Equity Forward: Environmental Justice Element. https://www.sandiego.gov/planning/environmental-justice-element

Cushing, L., Faust, J., August, L. M., Cendak, R., Wieland, W., & Alexeeff, G. (2015). Racial/ethnic disparities in cumulative environmental health impacts in California: Evidence from a statewide environmental justice screening tool. American Journal of Public Health, 105(11), 2341–2348. https://doi.org/10.2105/AJPH.2015.302643

IPUMS Contextual Determinants of Health. (2024). Income inequity measure documentation. University of Minnesota. https://cdoh.ipums.org/download

KPBS. (2018, April 5). Redlining’s mark on San Diego persists 50 years after housing protections. https://www.kpbs.org/news/midday-edition/2018/04/05/redlinings-mark-on-san-diego-persists

Mapping Prejudice. (n.d.). Maps and data: Racial covenants. University of Minnesota Libraries. https://mappingprejudice.umn.edu/racial-covenants/maps-data

Metropolitan Council. (n.d.). Repairing historic and ongoing harm. Imagine 2050. https://imagine2050.metrocouncil.org/reference-materials/housing/repairing-historic-and-ongoing-harm/

Minnesota State Climatology Office. (n.d.). Disproportionate heat risks for communities of color. https://climate.state.mn.us/disproportionate-heat-risks

San Diego Association of Governments. (2023). Social equity baseline report. https://www.sandag.org/-/media/SANDAG/Documents/PDF/regional-plan/social-equity-in-planning/social-equity-baseline-report-2023-04-01.pdf

Walker, R. H., Keeler, B. L., & Derickson, K. D. (2024). The impacts of racially discriminatory housing policies on the distribution of intra-urban heat and tree canopy: A comparison of racial covenants and redlining in Minneapolis, MN. Landscape and Urban Planning, 246, 105019. https://doi.org/10.1016/j.landurbplan.2024.105019

## Report generated with R version 4.6.1 (2026-06-24 ucrt).
## Key packages: sf 1.1.1, spdep 1.4.2, spatialreg 1.4.3.