Predicting County-Level Small Business Formation

This RMarkdown is the main analysis for the County-Level Small Business Formation project. It reads the modeling table produced by 01_data_pull.Rmd (saved at data/panel_clean.csv), performs exploratory analysis, fits three regression models (LASSO, Random Forest, kNN), compares them on a held-out test set, and produces the figures and tables that feed the written report.

1. Reading the modeling data

The data acquisition step (downloading from BDS and ACS, joining on FIPS, deriving the target and the six predictors) lives in the companion file 01_data_pull.Rmd. That workbook writes its output to data/panel_clean.csv, which we read in here.

panel_clean <- read_csv("data/panel_clean.csv", show_col_types = FALSE)

cat("Rows:", nrow(panel_clean), "  Columns:", ncol(panel_clean), "\n")

## Rows: 3069   Columns: 14

glimpse(panel_clean)

## Rows: 3,069
## Columns: 14
## $ fips                    <chr> "01001", "01003", "01005", "01007", "01009", "…
## $ name                    <chr> "Autauga County, Alabama", "Baldwin County, Al…
## $ startup_rate            <dbl> 1.1301341, 2.6089312, 1.0502080, 1.1737089, 1.…
## $ estabs_entry            <dbl> 67, 626, 26, 26, 84, 10, 33, 167, 34, 36, 61, …
## $ estabs                  <dbl> 783, 5827, 388, 297, 724, 105, 419, 2113, 554,…
## $ firms                   <dbl> 727, 5000, 351, 269, 660, 102, 383, 1747, 517,…
## $ emp                     <dbl> 9025, 71881, 6304, 4131, 7450, 1641, 6047, 356…
## $ total_pop               <dbl> 59285, 239945, 24757, 22152, 59292, 10157, 188…
## $ median_income           <dbl> 69841, 75019, 44290, 51215, 61096, 36723, 4488…
## $ pct_bachelors_or_higher <dbl> 28.282680, 32.797637, 11.464715, 11.468207, 15…
## $ pct_foreign_born        <dbl> 2.5790672, 3.8075392, 3.1344670, 1.2639942, 4.…
## $ lfp_rate                <dbl> 58.97954, 58.33333, 44.85755, 51.57901, 57.374…
## $ mean_commute            <dbl> 27.13771, 26.28998, 25.84178, 30.69759, 35.257…
## $ pct_25_44               <dbl> 26.15670, 23.20448, 26.29155, 27.85302, 23.849…

1.1 Locking the predictor set

Every model in this workbook uses the same six predictors. Storing them in a single named vector means no model can accidentally see a different feature set, and we can reference one canonical name throughout.

predictor_vars <- c(
  "median_income",
  "pct_bachelors_or_higher",
  "pct_foreign_born",
  "lfp_rate",
  "mean_commute",
  "pct_25_44"
)

1.2 Sanity check on missingness

The acquisition step already filtered out counties with any missing modeling value. We re-confirm here that the round trip through the CSV preserved that.

panel_clean |>
  select(startup_rate, all_of(predictor_vars)) |>
  summarize(across(everything(), ~ sum(is.na(.)))) |>
  kable(caption = "NA counts in modeling columns (should all be zero)")

NA counts in modeling columns (should all be zero)
startup_rate	median_income	pct_bachelors_or_higher	pct_foreign_born	lfp_rate	mean_commute	pct_25_44
0	0	0	0	0	0	0

2. Exploratory Data Analysis

2.1 Numerical summary

A skim() of the target plus the six predictors gives us the essential distributional facts (mean, sd, percentiles, histogram sparkline) in one frame.

panel_clean |>
  select(startup_rate, all_of(predictor_vars)) |>
  skim()

Data summary
Name	select(panel_clean, start…
Number of rows	3069
Number of columns	7
_______________________
Column type frequency:
numeric	7
________________________
Group variables	None

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
startup_rate	1	2.03	0.95	0.00	1.45	1.85	2.38	15.01	▇▁▁▁▁
median_income	1	66194.40	17444.07	25425.00	55082.00	63786.00	73855.00	178707.00	▃▇▁▁▁
pct_bachelors_or_higher	1	24.12	10.20	4.40	17.01	21.45	28.75	79.73	▆▇▂▁▁
pct_foreign_born	1	4.89	5.77	0.00	1.42	2.86	6.05	54.33	▇▁▁▁▁
lfp_rate	1	58.31	7.58	18.58	53.64	58.93	63.71	84.31	▁▁▆▇▁
mean_commute	1	24.44	5.70	5.02	20.61	24.19	27.95	49.41	▁▆▇▂▁
pct_25_44	1	23.86	3.27	6.93	22.00	23.62	25.58	47.81	▁▇▇▁▁

2.2 Distribution of the target

ggplot(panel_clean, aes(x = startup_rate)) +
  geom_histogram(bins = 60, fill = "steelblue", color = "white") +
  labs(title = "Distribution of county startup rate, 2023",
       subtitle = paste0("n = ", nrow(panel_clean), " counties"),
       x = "Startups per 1,000 residents",
       y = "Count of counties") +
  theme_minimal()

2.3 Predictor distributions

panel_clean |>
  select(all_of(predictor_vars)) |>
  pivot_longer(everything(), names_to = "var", values_to = "value") |>
  ggplot(aes(value)) +
  geom_histogram(bins = 40, fill = "steelblue", color = "white") +
  facet_wrap(~ var, scales = "free", ncol = 3) +
  labs(title = "Distribution of each predictor",
       x = NULL, y = "Count") +
  theme_minimal()

2.4 Predictor versus target

For each predictor, we plot it against the startup rate and overlay a LOESS smoother to surface any nonlinear shape. Each panel uses its own x-axis scale so the predictor’s natural range is preserved.

panel_clean |>
  select(startup_rate, all_of(predictor_vars)) |>
  pivot_longer(-startup_rate, names_to = "var", values_to = "value") |>
  ggplot(aes(value, startup_rate)) +
  geom_point(alpha = 0.2, size = 0.5, color = "gray40") +
  geom_smooth(method = "loess", color = "firebrick", se = FALSE, linewidth = 0.8) +
  facet_wrap(~ var, scales = "free_x", ncol = 3) +
  labs(title = "Each predictor vs startup rate",
       x = NULL,
       y = "Startups per 1,000 residents") +
  theme_minimal()

2.5 Geographic pattern

A county-level choropleth gives us the spatial story: where in the United States is business formation most concentrated?

plot_usmap(data = panel_clean,
           values = "startup_rate",
           regions = "counties",
           color = NA) +
  scale_fill_continuous(low = "white", high = "firebrick",
                        name = "Startups per 1,000",
                        na.value = "grey90") +
  labs(title = "County-level startup rate, 2023") +
  theme(legend.position = "right")

3. Train / Test Split

We use a single 80/20 random split for the entire workbook so every model is evaluated on the same held-out test set. The split is unstratified (a simple random sample of row indices); cross-validation for any model that needs hyperparameter tuning is configured locally inside each model section, not globally here.

set.seed(123)

n <- nrow(panel_clean)
test_index <- sample(seq_len(n), size = round(0.20 * n))
test_set   <- panel_clean[test_index, ]
train_set  <- panel_clean[-test_index, ]

cat("Train rows:", nrow(train_set),
    " Test rows:",  nrow(test_set), "\n")

## Train rows: 2455  Test rows: 614

# Confirm the random split did not produce a noticeably skewed test set
bind_rows(
  tibble(set = "train", startup_rate = train_set$startup_rate),
  tibble(set = "test",  startup_rate = test_set$startup_rate)
) |>
  group_by(set) |>
  summarize(n = n(),
            mean   = mean(startup_rate),
            median = median(startup_rate),
            sd     = sd(startup_rate)) |>
  kable(digits = 3, caption = "Target distribution: train vs test")

Target distribution: train vs test
set	n	mean	median	sd
test	614	2.005	1.808	1.037
train	2455	2.039	1.861	0.930

4. LASSO Regression

LASSO is the linear, shrinkage-based baseline. It penalizes the sum of absolute coefficients, which both regularizes against overfitting and performs automatic variable selection (some coefficients shrink exactly to zero). Tuning is done via 10-fold cross-validation on the training set; predictions on the held-out test set are then scored with postResample. Marginal-effect plots are produced for the top predictors by absolute standardized coefficient.

4.1 Fit

# Build model matrices (exclude intercept column)
x_train <- model.matrix(~ . - 1,
                        data = train_set |> select(all_of(predictor_vars)))
y_train <- train_set$startup_rate

x_test  <- model.matrix(~ . - 1,
                        data = test_set |> select(all_of(predictor_vars)))
y_test  <- test_set$startup_rate

set.seed(123)
cv_lasso <- cv.glmnet(x_train, y_train, alpha = 1)

plot(cv_lasso, main = "LASSO: CV MSE vs lambda")

cat("lambda.min =", cv_lasso$lambda.min, "\n")

## lambda.min = 0.01263224

cat("lambda.1se =", cv_lasso$lambda.1se, "\n")

## lambda.1se = 0.141901

4.2 Coefficients at lambda.min

coef(cv_lasso, s = "lambda.min")

## 7 x 1 sparse Matrix of class "dgCMatrix"
##                            lambda.min
## (Intercept)              2.187951e+00
## median_income            6.200927e-06
## pct_bachelors_or_higher  3.118455e-02
## pct_foreign_born         8.111984e-03
## lfp_rate                 .           
## mean_commute            -3.267528e-02
## pct_25_44               -2.316523e-02

4.3 Test-set performance

yhat_lasso <- as.numeric(predict(cv_lasso, newx = x_test, s = "lambda.min"))
lasso_perf <- postResample(yhat_lasso, y_test)
print(lasso_perf)

##      RMSE  Rsquared       MAE 
## 0.8963876 0.2584143 0.5321753

4.4 Marginal-effect plots

Hold every other predictor at its training-set mean and sweep one predictor across its 5th to 95th percentile, plotting the predicted startup rate. Showing the top 4 predictors by IQR-weighted impact (|coefficient| × IQR(predictor)), since raw coefficients are not comparable across predictors with very different scales (median_income is in tens of thousands while the percentages range 0 to 100).

# Identify top 4 predictors by IQR-weighted impact at lambda.min
lasso_coefs <- as.numeric(coef(cv_lasso, s = "lambda.min"))[-1]
names(lasso_coefs) <- predictor_vars

iqr_impact <- sapply(predictor_vars, function(v) {
  abs(lasso_coefs[[v]]) * IQR(train_set[[v]])
})

top4_lasso <- names(sort(iqr_impact, decreasing = TRUE))[1:4]
cat("Top 4 LASSO predictors by IQR-weighted impact:\n")

## Top 4 LASSO predictors by IQR-weighted impact:

print(top4_lasso)

## [1] "pct_bachelors_or_higher" "mean_commute"           
## [3] "median_income"           "pct_25_44"

# Predictor means for the marginal-effect grid
mean_vals <- train_set |> summarize(across(all_of(predictor_vars), mean))

for (var in top4_lasso) {
  var_range <- seq(quantile(train_set[[var]], 0.05),
                   quantile(train_set[[var]], 0.95),
                   length.out = 60)
  grid <- mean_vals |> slice(rep(1, length(var_range)))
  grid[[var]] <- var_range
  grid$yhat <- as.numeric(predict(cv_lasso,
                                  newx = as.matrix(grid[, predictor_vars]),
                                  s = "lambda.min"))

  p <- ggplot(grid, aes(.data[[var]], yhat)) +
    geom_line(color = "steelblue", linewidth = 1.2) +
    labs(title = paste("LASSO: marginal effect of", var),
         y = "Predicted startup rate (per 1,000)",
         x = var) +
    theme_minimal()
  print(p)
}

5. Random Forest

Random Forest is the tree-based ensemble. It grows many decision trees on bootstrap samples of the training data, randomly subsets the predictor pool at each split, and averages predictions across the forest. This decorrelates the trees and reduces variance. The main tuning parameter is mtry (the number of predictors considered at each split), tuned here with 10-fold cross-validation over four candidate values.

5.1 Fit and tune

set.seed(123)
rf_grid <- expand.grid(mtry = c(2, 3, 4, 5))

rf_caret <- train(
  reformulate(predictor_vars, "startup_rate"),
  data      = train_set,
  method    = "rf",
  trControl = trainControl(method = "cv", number = 10),
  tuneGrid  = rf_grid,
  ntree     = 500
)

ggplot(rf_caret, highlight = TRUE) +
  labs(title = "Random Forest: 10-fold CV RMSE vs mtry") +
  theme_minimal()

cat("Best mtry:\n")

## Best mtry:

print(rf_caret$bestTune)

##   mtry
## 1    2

5.2 Test-set performance

yhat_rf <- predict(rf_caret, newdata = test_set)
rf_perf <- postResample(yhat_rf, test_set$startup_rate)
print(rf_perf)

##      RMSE  Rsquared       MAE 
## 0.8791867 0.2826909 0.5083562

5.3 Variable importance

varImp(rf_caret) |>
  plot(main = "Random Forest: variable importance")

5.4 Partial-dependence plots

For each of the top 4 predictors by importance, sweep that predictor across its observed range while integrating out the others. Partial dependence is the model’s view of the marginal effect after accounting for predictor interactions, the natural counterpart to the LASSO marginal-effect plots in Section 4.

top4_rf <- varImp(rf_caret)$importance |>
  tibble::rownames_to_column("var") |>
  arrange(desc(Overall)) |>
  slice(1:4) |>
  pull(var)

cat("Top 4 RF predictors by importance:\n")

## Top 4 RF predictors by importance:

print(top4_rf)

## [1] "pct_bachelors_or_higher" "mean_commute"           
## [3] "pct_25_44"               "median_income"

for (v in top4_rf) {
  p <- partial(rf_caret, pred.var = v, train = train_set,
               plot = TRUE, plot.engine = "ggplot2",
               rug = TRUE) +
    labs(title = paste("RF partial dependence:", v),
         y = "Partial dependence (predicted startup rate)") +
    theme_minimal()
  print(p)
}

6. kNN Regression

kNN is the nonparametric, distance-based comparison. Each prediction is the average of the \(k\) nearest training counties in (centered and scaled) predictor space. Because kNN is sensitive to the scale of predictors (median income would otherwise dominate any distance calculation), we center and scale every predictor before fitting. The single tuning parameter is \(k\), the neighborhood size; we tune it via 10-fold cross-validation across a wide grid.

6.1 Fit and tune k

set.seed(123)
knn_caret <- train(
  reformulate(predictor_vars, "startup_rate"),
  data       = train_set,
  method     = "knn",
  trControl  = trainControl(method = "cv", number = 10),
  preProcess = c("center", "scale"),
  tuneGrid   = data.frame(k = seq(5, 80, by = 5))
)

ggplot(knn_caret, highlight = TRUE) +
  labs(title = "kNN: 10-fold CV RMSE vs k") +
  theme_minimal()

cat("Best k:\n")

## Best k:

print(knn_caret$bestTune)

##    k
## 6 30

6.2 Test-set performance

yhat_knn <- predict(knn_caret, newdata = test_set)
knn_perf <- postResample(yhat_knn, test_set$startup_rate)
print(knn_perf)

##      RMSE  Rsquared       MAE 
## 0.8779608 0.2979948 0.4922012

6.3 Marginal-effect plots

kNN does not provide coefficients or built-in variable-importance scores in the same form as LASSO and RF. We construct marginal-effect plots the same way as the LASSO section (sweep one predictor across its 5th to 95th percentile while holding the others at their training-set mean), and pick the top 4 predictors by their estimated 25th-to-75th-percentile prediction impact.

mean_vals <- train_set |> summarize(across(all_of(predictor_vars), mean))

# IQR-based impact for kNN: predict at 25th vs 75th percentile, hold others at mean
knn_impact <- sapply(predictor_vars, function(v) {
  q25 <- mean_vals; q25[[v]] <- quantile(train_set[[v]], 0.25)
  q75 <- mean_vals; q75[[v]] <- quantile(train_set[[v]], 0.75)
  abs(predict(knn_caret, newdata = q75) - predict(knn_caret, newdata = q25))
})
names(knn_impact) <- predictor_vars

top4_knn <- names(sort(knn_impact, decreasing = TRUE))[1:4]
cat("Top 4 kNN predictors by 25th-to-75th-percentile impact:\n")

## Top 4 kNN predictors by 25th-to-75th-percentile impact:

print(top4_knn)

## [1] "pct_bachelors_or_higher" "pct_25_44"              
## [3] "pct_foreign_born"        "lfp_rate"

for (var in top4_knn) {
  var_range <- seq(quantile(train_set[[var]], 0.05),
                   quantile(train_set[[var]], 0.95),
                   length.out = 60)
  grid <- mean_vals |> slice(rep(1, length(var_range)))
  grid[[var]] <- var_range
  grid$yhat <- predict(knn_caret, newdata = grid)

  p <- ggplot(grid, aes(.data[[var]], yhat)) +
    geom_line(color = "coral", linewidth = 1.2) +
    labs(title = paste("kNN: marginal effect of", var),
         y = "Predicted startup rate (per 1,000)",
         x = var) +
    theme_minimal()
  print(p)
}

7. Model Comparison

All three models were fit on the same 80/20 train-test split and tuned with 10-fold cross-validation on the training set only. Test-set performance is therefore directly comparable. Below we put RMSE, R-squared, and MAE next to each other, then plot predicted-vs-actual and residuals for each model.

7.1 Performance table

comparison <- tibble(
  Model    = c("LASSO", "Random Forest", "kNN"),
  RMSE     = c(lasso_perf["RMSE"],     rf_perf["RMSE"],     knn_perf["RMSE"]),
  Rsquared = c(lasso_perf["Rsquared"], rf_perf["Rsquared"], knn_perf["Rsquared"]),
  MAE      = c(lasso_perf["MAE"],      rf_perf["MAE"],      knn_perf["MAE"])
)
kable(comparison, digits = 4,
      caption = "Test-set performance: 80/20 hold-out, 10-fold CV for tuning")

Test-set performance: 80/20 hold-out, 10-fold CV for tuning
Model	RMSE	Rsquared	MAE
LASSO	0.8964	0.2584	0.5322
Random Forest	0.8792	0.2827	0.5084
kNN	0.8780	0.2980	0.4922

7.2 Predicted vs actual

Each panel shows test-set predicted startup rate against actual startup rate; the dashed red line is y = x (perfect prediction). Points hugging the diagonal indicate accurate predictions; points scattered far from it are model errors.

preds_long <- tibble(
  actual = rep(test_set$startup_rate, 3),
  pred   = c(yhat_lasso, yhat_rf, yhat_knn),
  model  = factor(rep(c("LASSO", "Random Forest", "kNN"),
                      each = nrow(test_set)),
                  levels = c("LASSO", "Random Forest", "kNN"))
)

ggplot(preds_long, aes(actual, pred)) +
  geom_point(alpha = 0.3, size = 0.7) +
  geom_abline(slope = 1, intercept = 0,
              color = "firebrick", linetype = "dashed", linewidth = 0.8) +
  facet_wrap(~ model) +
  labs(title = "Predicted vs actual startup rate (test set)",
       x = "Actual (per 1,000)",
       y = "Predicted (per 1,000)") +
  theme_minimal()

7.3 Residuals vs predicted

Residuals (actual minus predicted) plotted against predicted values surface systematic over- or under-prediction; for an unbiased model the points should scatter symmetrically around zero across the full range of predictions.

preds_long |>
  mutate(residual = actual - pred) |>
  ggplot(aes(pred, residual)) +
  geom_point(alpha = 0.3, size = 0.7) +
  geom_hline(yintercept = 0,
             color = "firebrick", linetype = "dashed", linewidth = 0.8) +
  facet_wrap(~ model) +
  labs(title = "Residuals vs predicted (test set)",
       x = "Predicted startup rate (per 1,000)",
       y = "Residual (actual - predicted)") +
  theme_minimal()

8. Worked-example predictions for specific counties

Picking a handful of named counties to walk through what each model predicts and how that compares to the actual 2023 startup rate. Counties are chosen to span the range: a known urban-tech metro (Travis County, TX, home to Austin), a large industrial metro (Cook County, IL, home to Chicago), a county that sits very near the demographic medians of the dataset, and a rural midwestern county at the low end of the startup-rate distribution.

# Hand-picked counties spanning a range of profiles
named_picks <- c("48453",   # Travis County, TX  (Austin, tech-heavy metro)
                 "17031",   # Cook County, IL    (Chicago, large industrial metro)
                 "17175")   # Stark County, IL   (rural midwestern)

# Plus a programmatically-chosen "median demographics" county from the test set:
# pick the test-set county whose vector of standardized predictors is closest to the origin
predictor_means <- sapply(predictor_vars, function(v) mean(train_set[[v]]))
predictor_sds   <- sapply(predictor_vars, function(v) sd(train_set[[v]]))

z_scores <- as.matrix(test_set[, predictor_vars])
z_scores <- sweep(z_scores, 2, predictor_means, "-")
z_scores <- sweep(z_scores, 2, predictor_sds,   "/")
distances <- sqrt(rowSums(z_scores^2))
median_county_fips <- test_set$fips[which.min(distances)]
cat("Programmatically-chosen median-demographic county FIPS:",
    median_county_fips, "\n")

## Programmatically-chosen median-demographic county FIPS: 01103

example_fips <- c(named_picks, median_county_fips)

examples <- panel_clean |>
  filter(fips %in% example_fips)

# Display predictor values
examples |>
  select(fips, name, startup_rate, all_of(predictor_vars)) |>
  kable(digits = 2,
        caption = "Selected counties: actual startup rate and predictor values")

Selected counties: actual startup rate and predictor values
fips	name	startup_rate	median_income	pct_bachelors_or_higher	pct_foreign_born	lfp_rate	mean_commute	pct_25_44
01103	Morgan County, Alabama	1.62	64858	23.24	5.50	58.32	23.97	24.57
17031	Cook County, Illinois	2.31	81797	41.95	21.40	66.20	31.85	29.76
17175	Stark County, Illinois	2.82	62284	19.27	2.08	58.12	26.66	21.77
48453	Travis County, Texas	3.99	97169	55.51	17.26	73.19	25.40	36.32

# Build prediction inputs in the same shape each model expects
ex_x_lasso <- model.matrix(~ . - 1,
                           data = examples |> select(all_of(predictor_vars)))

examples_pred <- examples |>
  transmute(
    fips, name,
    actual = round(startup_rate, 2),
    LASSO  = round(as.numeric(predict(cv_lasso, newx = ex_x_lasso,
                                      s = "lambda.min")), 2),
    RF     = round(predict(rf_caret, newdata = examples), 2),
    kNN    = round(predict(knn_caret, newdata = examples), 2)
  )

kable(examples_pred,
      caption = "Predicted vs actual startup rate (per 1,000) for selected counties")

Predicted vs actual startup rate (per 1,000) for selected counties
fips	name	actual	LASSO	RF	kNN
01103	Morgan County, Alabama	1.62	2.01	1.76	1.85
17031	Cook County, Illinois	2.31	2.45	2.57	2.66
17175	Stark County, Illinois	2.82	1.82	1.58	1.62
48453	Travis County, Texas	3.99	2.99	3.84	3.41

9. Limitations and possible enhancements

Several limitations of this analysis are worth naming explicitly before the report draws conclusions from it.

Spatial autocorrelation. Counties next to each other are not statistically independent, but every model in this workbook treats each county as an independent observation. A spatial random-effects model (e.g., a geographically weighted regression or an ICAR random effect by state) would explicitly account for the fact that, say, a Sun Belt county is likely correlated with its neighbors.

Single-year cross-section. The analysis uses 2023 BDS only. Counties that are rapidly gaining or losing population, or are riding a multi-year boom or bust trajectory, look identical to a snapshot model. A panel approach pooling 2018–2023 would give us within-county variation and meaningfully more predictive power.

Small-population county noise. As the worked example for Stark County, IL shows, a small-population county where one or two new establishments register can post a high per-capita startup rate that none of our models will predict from demographics alone. We saw this also in the choropleth as bright red dots in sparsely populated rural areas. A multi-year aggregation, or modeling raw counts instead of rates with a population offset, would dampen this noise.

Heavy-tail compression. All three models compress predictions toward the mean. None of them predict above ~5 startups per 1,000 even when actual values reach 15. The Travis County example shows the closest the models came to a high-end prediction (RF at 3.84 vs actual 3.99), but the dataset’s rare extreme counties remain inherently hard to learn from a six-predictor model.

Predictor scope. Six predictors is conservative. Plausible additions, all of which are publicly available and could plug into the same modeling pipeline, include: SBA 7(a) loan volume by county (access to capital), FCC broadband speed (not just adoption), state-level tax climate, and a categorical region or state effect.

Modeling alternatives. Boosted trees (gbm or xgboost) would likely beat all three models on raw RMSE at the cost of more tuning. A generalized additive model (GAM) with smoothers on each predictor would produce smoother, more interpretable nonlinearities than RF while still capturing the J-shapes that LASSO smoothed away.

Analysis: County-Level Small Business Formation

Lennox Magak

2026-04-20