Tutorial: Bivariate Relationships and Spatial Regression — Completed
July 21, 2026
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)Tutorial
This tutorial looks at the relationship between educational attainment and the decrease in workplace mobility during the first months of COVID-19. Correlation, scatterplots, and bivariate maps describe the relationship. OLS and spatial regression models then test the relationship while accounting for spatial dependence.
Read in the data
library(tidyverse)
library(tmap)
library(sf)
library(spatialreg)
library(spdep)
county_covid <- st_read(
"https://drive.google.com/uc?export=download&id=15njmLhQMvg4KNpFthR0XyS1LD3reotus",
quiet = TRUE
) |>
# Convert the original negative value so larger values mean a larger
# decrease in travel to workplaces.
mutate(workplace_change = workplace_change * -1)
tmap_mode("plot")Testing correlation
cor(county_covid$pct_bachelors, county_covid$workplace_change)[1] 0.7981348Q1: What does the correlation tell us about the strength and direction of the relationship?
The correlation is about 0.80, so there is a strong positive relationship. Counties with a higher percentage of residents with a bachelor’s degree generally had a larger decrease in workplace mobility during the early months of COVID-19.
Scatterplot
ggplot(
county_covid,
aes(x = pct_bachelors, y = workplace_change)
) +
geom_point() +
labs(
x = "Population with a bachelor's degree (%)",
y = "Decrease in workplace mobility"
)
The upward pattern matches the correlation result. The relationship is fairly strong, although a few counties sit well above or below the main pattern.
Bivariate map
tm_shape(county_covid) +
tm_polygons(
fill = tm_vars(
c("pct_bachelors", "workplace_change"),
multivariate = TRUE
),
fill.scale = tm_scale_bivariate(
values = "bivario.folk_warmth",
scale1 = tm_scale_intervals(
n = 4,
style = "quantile",
labels = c("Lo", "", "", "Hi")
),
scale2 = tm_scale_intervals(
n = 4,
style = "quantile",
labels = c("Lo", "", "", "Hi")
)
)
)
Q2: What does the bivariate map tell us about where the relationship is strongest?
The strongest high-high pattern appears in major metropolitan counties around the Triangle, Charlotte, the Triad, and Asheville. These counties had both higher bachelor’s-degree attainment and larger decreases in workplace mobility. Several rural eastern and southeastern counties show the opposite low-low pattern. A few counties fall into mixed categories, so the relationship is strong overall but not identical everywhere.
Ordinary least squares regression
ols_model <- lm(
workplace_change ~ pct_bachelors,
data = county_covid
)
summary(ols_model)Call:
lm(formula = workplace_change ~ pct_bachelors, data = county_covid)
Residuals:
Min 1Q Median 3Q Max
-10.231 -2.797 -0.116 2.187 18.263
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.7051 1.1202 15.80 <2e-16
pct_bachelors 0.5452 0.0420 12.98 <2e-16
Residual standard error: 4.254 on 96 degrees of freedom
Multiple R-squared: 0.6370, Adjusted R-squared: 0.6332
F-statistic: 168.5 on 1 and 96 DF, p-value: < 2.2e-16The coefficient is positive and statistically significant. A one-percentage-point increase in bachelor’s-degree attainment is associated with about a 0.55-point larger decrease in workplace mobility. The model explains about 63% of the county-level variation. This is an association and does not show that educational attainment caused mobility to change.
Testing the OLS residuals for spatial autocorrelation
county_covid <- county_covid |>
mutate(olsresid = resid(ols_model))
tm_shape(county_covid) +
tm_polygons(
fill = "olsresid",
fill.scale = tm_scale_intervals(
values = "brewer.bu_pu",
style = "quantile",
n = 5
)
)
nb <- poly2nb(county_covid, queen = TRUE)
nbw <- nb2listw(nb, style = "W", zero.policy = TRUE)
gmoran <- moran.test(
county_covid$olsresid,
nbw,
zero.policy = TRUE
)
gmoranMoran I test under randomisation
Moran I statistic = 0.330622
Expectation = -0.010309
Variance = 0.004134
Standard deviate = 5.3026
p-value = 5.71e-08Q3: Do the results indicate that we should use a spatial model?
Yes. The residual Moran’s I is positive and highly significant, which means the OLS residuals are still spatially clustered. The independence assumption is violated, so a spatial regression model is more appropriate.
Spatial regression
Spatial lag model
fit_lag <- lagsarlm(
workplace_change ~ pct_bachelors,
data = county_covid,
listw = nbw,
zero.policy = TRUE
)
summary(fit_lag)Spatial lag model
Estimate Std. Error z value p-value
(Intercept) 5.65338 2.40201 2.3536 0.01859
pct_bachelors 0.47263 0.03912 12.0814 < 2e-16
Rho: 0.44395, p-value = 7.63e-08
Log likelihood: -268.0403
AIC: 544.08 (OLS AIC: 565.86)
LM test for residual autocorrelation: p-value = 0.35393The lag model still shows a positive association: a one-point increase in bachelor’s-degree attainment is associated with about a 0.47-point increase in workplace change. Rho is positive and significant, meaning a county’s workplace change is related to values in neighboring counties. The AIC improves from 565.86 to 544.08, and the remaining residual spatial autocorrelation is not significant.
Spatial error model
fit_error <- errorsarlm(
workplace_change ~ pct_bachelors,
data = county_covid,
listw = nbw,
zero.policy = TRUE
)
summary(fit_error)Spatial error model
Estimate Std. Error z value p-value
(Intercept) 18.89736 1.29866 14.552 < 2e-16
pct_bachelors 0.50168 0.03924 12.786 < 2e-16
Lambda: 0.59034, p-value = 1.04e-09
Log likelihood: -268.5434
AIC: 545.09 (OLS AIC: 565.86)The error model estimates that a one-point increase in bachelor’s-degree attainment is associated with about a 0.50-point increase in workplace change. Lambda is positive and significant, showing that unmeasured factors affecting mobility are spatially clustered. Its AIC is also much lower than the OLS AIC.
Both spatial models fit much better than OLS and lead to the same basic conclusion about education. The lag model has a slightly smaller AIC, but the difference is only about one point. The error model is a reasonable choice if shared policies, commuting systems, employment patterns, or other unmeasured regional conditions are creating the spatial dependence rather than one county’s mobility directly changing its neighbor’s mobility.
Mini Challenge
acs_tract_nc <- st_read(
"https://drive.google.com/uc?export=download&id=1cnz4xgdDRZvlXzN3IyvCxOETRWfrpw-0",
quiet = TRUE
)I chose percent of the population below the poverty level as the predictor because it has a clear connection to household income while still measuring a different tract characteristic.
Explore the variables
income_map <- tm_shape(acs_tract_nc) +
tm_polygons(
fill = "median_hh_inc",
fill.scale = tm_scale_intervals(style = "quantile", n = 5),
col_alpha = 0
) +
tm_title("Median household income")
poverty_map <- tm_shape(acs_tract_nc) +
tm_polygons(
fill = "pct_under_poverty_level",
fill.scale = tm_scale_intervals(style = "quantile", n = 5),
col_alpha = 0
) +
tm_title("Population below poverty level")
tmap_arrange(income_map, poverty_map, ncol = 2)
The maps show a general inverse pattern. Many tracts with higher poverty have lower median household income, while higher-income areas tend to have lower poverty.
ggplot(
acs_tract_nc,
aes(x = pct_under_poverty_level, y = median_hh_inc)
) +
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", se = FALSE) +
scale_y_continuous(labels = scales::dollar) +
labs(
x = "Population below the poverty level (%)",
y = "Median household income"
)
Fit and interpret the OLS model
mini_ols <- lm(
median_hh_inc ~ pct_under_poverty_level,
data = acs_tract_nc
)
summary(mini_ols)Call:
lm(formula = median_hh_inc ~ pct_under_poverty_level,
data = acs_tract_nc)
Residuals:
Min 1Q Median 3Q Max
-52016 -16166 -4967 9625 149503
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 105088.88 865.66 121.40 <2e-16
pct_under_poverty_level -2125.38 50.06 -42.45 <2e-16
Residual standard error: 25038 on 2496 degrees of freedom
Multiple R-squared: 0.4193, Adjusted R-squared: 0.4191
F-statistic: 1802.2 on 1 and 2496 DF, p-value: < 2.2e-16There is a strong, statistically significant negative relationship. Each one-percentage-point increase in the tract poverty rate is associated with about $2,125 lower median household income. The model predicts an income of about $105,089 when the poverty rate is zero, although that intercept is mainly a reference point. Poverty alone explains about 42% of the variation in tract median household income. This is a bivariate association, not evidence that the poverty rate by itself causes the income difference, and the spatial independence of these tract-level residuals has not been tested here.
Overall takeaway
The county analysis shows a strong relationship between bachelor’s-degree attainment and decreased workplace mobility, but OLS leaves significant spatial structure in its residuals. Both spatial models address that dependence and fit better. In the mini challenge, tract poverty is strongly associated with lower median household income, but the model should still be interpreted as descriptive rather than causal.
install.packages(“rsconnect”) # only needed once
result <- rsconnect::rpubsUpload( title = “Spatial Modeling Tutorial”, contentFile = “C:/Users/sonia/Downloads/vasa_spatial_modeling.html”, originalDoc = NULL )
browseURL(result$continueUrl)