Week 6: Feature Engineering for Linear Regression

Logs, squared terms, categories, interactions, and district fixed effects in R

Published

October 9, 2026

What We Did

Watched the recorded lecture on feature engineering.
Used 2023 ACS five-year data to predict median gross rent for Philadelphia census tracts.
Compared transformations and model specifications using the same ten cross-validation folds.
Reviewed the code and checks needed for Part 5 of the assignment.

Key Concepts

1. Feature engineering and linear models

Feature engineering changes the variables supplied to a model. Logs, squared terms, and interactions let linear regression represent relationships that a straight line in the original variables would miss.

\[ \hat y = \hat\beta_0 + \hat\beta_1 x + \hat\beta_2 x^2 \]

This is still a linear regression: the model is linear in its coefficients, even though its relationship with $x$ is curved.

The lecture evaluates each specification by cross-validation RMSE. A transformation needs to improve predictions on held-out data; a better training fit alone is not enough.

2. Build the predictors and keep the sample consistent

The baseline uses four predictors:

Variable	Meaning	Construction or unit
`income`	Median household income	Dollars
`pct_renter`	Share of households that rent	`100 * hh_renterE / hh_totalE`
`bldg_age`	Age based on median year structure built	`2023 - yr_builtE`
`dist_ch_mi`	Distance from tract geometry to City Hall	Feet divided by 5,280

In the wide ACS output, E identifies the estimate and M identifies the margin of error. The lecture renames rentE to rent and retains rentM as rent_moe.

city_hall <- st_sfc(
  st_point(c(-75.1636, 39.9526)),
  crs = 4326
) %>%
  st_transform(2272)

tracts <- tracts_raw %>%
  mutate(
    rent       = rentE,
    rent_moe   = rentM,
    income     = incomeE,
    bldg_age   = 2023 - yr_builtE,
    pct_renter = 100 * hh_renterE / hh_totalE,
    dist_ch_mi = as.numeric(st_distance(geometry, city_hall)) / 5280
  ) %>%
  select(GEOID, rent, rent_moe, income, bldg_age, pct_renter, dist_ch_mi) %>%
  filter(
    !is.na(rent), !is.na(income),
    !is.na(bldg_age), !is.na(pct_renter)
  )

tracts_raw has already been transformed to EPSG:2272, whose units are feet. The City Hall point starts in EPSG:4326, with coordinates ordered as longitude, latitude, and is transformed to match.

st_distance() here measures the shortest distance between each tract polygon and the point. It does not measure distance from the tract centroid. A tract containing City Hall has distance zero.

The cleaned lecture sample contains 377 tracts. All specifications use these rows. If a later feature introduces missing or invalid values, handle them before comparing models and rebuild the folds if the sample changes.

3. Baseline model and fixed-fold cross-validation

f0 <- rent ~ income + pct_renter + bldg_age + dist_ch_mi
m0 <- lm(f0, data = tracts)

f0 stores the formula. lm() estimates the coefficients, and m0 stores the fitted model. predict(m0, newdata = ...) applies those coefficients to another dataset.

The lecture assigns each row to one of ten folds once:

set.seed(5080)
folds <- sample(rep(1:10, length.out = nrow(tracts)))

cv_rmse <- function(f, data, log_outcome = FALSE) {
  preds <- numeric(nrow(data))

  for (k in 1:10) {
    m <- lm(f, data = data[folds != k, ])
    p <- predict(m, newdata = data[folds == k, ])

    if (log_outcome) p <- exp(p)

    preds[folds == k] <- p
  }

  sqrt(mean((data$rent - preds)^2))
}

folds != k selects the nine training folds.
folds == k selects the held-out fold.
preds[folds == k] <- p puts predictions back into their original row positions.
After ten iterations, every row has one out-of-fold prediction.

\[ RMSE = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_i-\hat y_i)^2} \]

This helper computes one RMSE across all out-of-fold predictions. It does not take an unweighted average of ten separate fold RMSEs. RMSE has the same units as the outcome and gives more weight to large errors than mean absolute error does.

The function uses the external object folds, so every input dataset must preserve the original row count and order. Reordering the data changes which observations belong to each fold.

model_log <- tibble(
  model = "m0",
  added = "baseline",
  cv_rmse = cv_rmse(f0, tracts)
)

model_log is a record of the models tried. The word log in this name has nothing to do with a logarithmic transformation.

4. Log transformations and coefficient interpretation

A log compresses a long right tail. It also makes equal proportional changes in a predictor have equal distances on the transformed scale. In R, log() uses the natural logarithm, and exp() reverses it.

f_logx <- rent ~ log(income) + pct_renter + bldg_age + dist_ch_mi
m_logx <- lm(f_logx, data = tracts)

For this model, the coefficient on log(income) is about 502.18. Holding the other predictors constant, a 1% increase in income corresponds to about $5.02 higher predicted rent.

Model	Formula	Interpretation of $\beta$, holding other variables constant
Level-level	`y ~ x`	+1 unit of X corresponds to +$\beta$ units of Y
Level-log	`y ~ log(x)`	+1% in X corresponds to approximately +$\beta/100$ units of Y
Log-level	`log(y) ~ x`	+1 unit of X corresponds to approximately $100\beta$% change in Y
Log-log	`log(y) ~ log(x)`	+1% in X corresponds to approximately $\beta$% change in Y

For level-log models, a proportional increase $r$ in X gives the exact predicted change:

\[ \Delta\hat y = \beta\ln(1+r) \]

b <- coef(m_logx)["log(income)"]
b / 100         # About $5.02 for a 1% income increase
b * log(1.10)   # About $47.86 for a 10% income increase

For log-level models, the exact percentage change in the exponentiated prediction for a one-unit increase in X is:

\[ 100\left(e^\beta-1\right)\% \]

The log-log income coefficient in the lecture is about 0.363. A 1% increase in income corresponds to approximately 0.363% higher predicted rent. This coefficient is an elasticity.

Zeros and `log1p()`

log(0)     # -Inf
log1p(0)   # 0
log1p(x)   # log(1 + x)

Counts such as crimes or 311 calls can include meaningful zeros. log1p() preserves those observations when transforming nonnegative counts. Its interpretation concerns x + 1, so the usual percentage interpretation for log(x) is less accurate near zero. Neither function makes arbitrary negative values suitable for logging.

Logging Y changes the prediction scale

f_logy <- log(rent) ~ log(income) + pct_renter + bldg_age + dist_ch_mi
cv_rmse(f_logy, tracts, log_outcome = TRUE)

Predictions from this model are on the log-rent scale. Apply exp() to each prediction before comparing it with observed rent in dollars. A log-scale RMSE of 0.188 cannot be compared with a dollar RMSE of 265, and exponentiating the final RMSE does not fix this.

Direct exponentiation does not generally recover the conditional mean of rent; it can underpredict that mean. The lecture uses the simple back-transformation for model comparisons and notes that bias corrections exist.

A skewed predictor does not violate a requirement that X be normally distributed. Linear regression has no such requirement. Use plots and CV to decide whether a log is useful.

5. Squared terms and changing slopes

f_sq <- rent ~ log(income) + pct_renter +
  bldg_age + I(bldg_age^2) + dist_ch_mi

m_sq <- lm(f_sq, data = tracts)

Keep both the original variable and its square. I() tells R to perform arithmetic inside the formula; ^ otherwise has a special formula meaning.

For the age component:

\[ \hat y = \cdots + \beta_1 age + \beta_2 age^2 \]

The slope depends on age:

\[ \frac{d\hat y}{d age} = \beta_1 + 2\beta_2 age \]

The lecture estimates $\beta_1\approx-28.004$ and $\beta_2\approx0.1823$. The positive squared coefficient gives an upward-opening curve. The linear coefficient alone does not mean that rent falls by $28 for every additional year of age throughout the sample.

The turning point is:

\[ age^* = -\frac{\beta_1}{2\beta_2} \]

b1 <- coef(m_sq)["bldg_age"]
b2 <- coef(m_sq)["I(bldg_age^2)"]

-b1 / (2 * b2)         # About 76.8 years
range(tracts$bldg_age) # 5 to 85 years

Check whether the turning point lies inside the observed range. A turning point outside that range does not establish an observed reversal. A quadratic also keeps curving outside the training range, so extrapolated predictions can become unreasonable.

The old-building pile-up in the lecture plot partly comes from Census coding of early construction years. Check coding limits before interpreting a concentration of observations as a housing-market pattern.

6. Categories, dummy variables, and reference groups

Use a categorical predictor for existing groups, such as zoning classes, or when a continuous relationship is better represented by distinct groups.

tracts <- tracts %>%
  mutate(
    ring = case_when(
      dist_ch_mi < 2 ~ "Core",
      dist_ch_mi < 5 ~ "Inner",
      TRUE ~ "Outer"
    ),
    ring = factor(ring, levels = c("Core", "Inner", "Outer"))
  )

case_when() uses the first matching condition. The Inner group therefore contains distances from 2 up to, but not including, 5 miles. Core is the reference because it is the first specified factor level.

For numeric breakpoints, use cut():

tracts <- tracts %>%
  mutate(renter_band = cut(
    pct_renter,
    breaks = c(0, 33, 66, 100),
    labels = c("Mostly owners", "Mixed", "Mostly renters"),
    include.lowest = TRUE
  ))

With an intercept and the usual treatment coding, R represents $k$ categories using $k-1$ dummy variables:

Ring	`ringInner`	`ringOuter`
Core	0	0
Inner	1	0
Outer	0	1

model.matrix(~ ring, data = st_drop_geometry(tracts))

In the lecture model, ringInner is about -122 and ringOuter is about -164. Holding the other predictors constant, predicted rent in Inner tracts is about $122 lower than in Core tracts, and predicted rent in Outer tracts is about $164 lower than in Core tracts. These coefficients are differences from Core, not rent levels.

tracts_r <- tracts %>%
  mutate(ring = relevel(ring, ref = "Outer"))

Changing the reference changes the coefficient comparisons while leaving predictions unchanged. Numeric IDs also need factor() when their numbers are labels; treating district codes 1 through 18 as a continuous predictor would impose a numerical slope across arbitrary codes.

Binning removes information. Distances of 2.1 and 4.9 miles both become Inner. In this example, replacing continuous distance with rings increased CV RMSE.

7. Interactions: different slopes by group

An additive category term shifts the intercept. An interaction allows the relationship with another predictor to vary by group.

f_int <- rent ~ log(income) * ring + pct_renter +
  bldg_age + I(bldg_age^2)

In an R formula:

x * group
# Expands to: x + group + x:group

Writing $x=\ln(income)$, the relevant terms are:

\[ \hat y = \beta_0 + \beta_1x + \gamma_I D_I + \gamma_O D_O + \delta_I xD_I + \delta_O xD_O + \cdots \]

Ring	Intercept, with other predictors set to zero	Slope on log income
Core	$\beta_0$	$\beta_1$
Inner	$\beta_0+\gamma_I$	$\beta_1+\delta_I$
Outer	$\beta_0+\gamma_O$	$\beta_1+\delta_O$

The lecture coefficients give these log-income slopes:

Core: 653.70.
Inner: $653.70-219.10=434.60$.
Outer: $653.70-244.42=409.28$.

The main log-income coefficient is the Core slope. Each interaction coefficient is a difference from that slope. With an interaction, the Outer-Core rent difference also depends on income: $\gamma_O+\delta_Ox$.

Use interactions when the relationship or a plot suggests different slopes. Each additional group slope requires another coefficient estimated from the same sample.

8. District fixed effects

The lecture uses a district factor to allow separate intercepts for Philadelphia’s planning districts:

\[ y_{ig} = \alpha_g + \beta x_{ig} + \epsilon_{ig} \]

$\alpha_g$ is the intercept for district $g$. It captures district-level differences left unexplained by the measured predictors, without identifying their individual causes.

tracts <- tracts %>%
  st_join(districts %>% select(district = dist_name), largest = TRUE) %>%
  mutate(district = factor(district))

f_fe <- rent ~ log(income) + pct_renter +
  bldg_age + I(bldg_age^2) + district

largest = TRUE assigns a tract crossing district boundaries to the district with the largest overlapping area. It assigns attributes to the tract rather than splitting it into pieces.

After adding district fixed effects, the log-income coefficient changes from about 519.74 to 427.39. The income relationship now uses variation within districts after controlling for their intercepts. Some of the earlier association reflected differences between districts.

Fixed effects do not account for every unmeasured difference within a district, and they do not turn an association into a causal effect. Residual maps can show whether geographic clusters remain after adding them.

Small groups and unseen factor levels

tracts %>%
  st_drop_geometry() %>%
  count(district) %>%
  arrange(n)

One lecture district has only one tract. When that tract is held out, the training set contains no observation from its district, so prediction encounters a factor level the model has never learned.

The lecture handles groups with fewer than ten observations by assigning their tracts to the district of the nearest tract in a sufficiently large group. It keeps the original district label and creates district_cv for the CV specification. This changes group membership; it does not dissolve or combine tract geometries.

Other options discussed include a coarser geography or pooling small groups into an Other category. A geographically sensible merge gives the resulting category a clearer meaning. About ten observations per group is a classroom rule of thumb for ten-fold CV, not a guarantee that every training fold will contain every category.

9. R syntax to remember

Code	Purpose or check
`log(x)`	Natural log; requires positive values for finite results
`log1p(x)`	Compute `log(1 + x)`; accommodates zero counts
`exp(pred)`	Back-transform a log prediction
`x + I(x^2)`	Include a linear term and its square
`factor(x)`	Treat labels as categories
`relevel(x, ref = "...")`	Choose a reference category
`case_when()` / `cut()`	Construct categories from conditions or breakpoints
`x * group`	Include both main effects and their interaction
`model.matrix()`	Inspect the encoded predictors
`count(group)`	Check group sizes before CV
`st_join(..., largest = TRUE)`	Assign the spatial match with the largest overlap
`add_row()`	Record another model and its CV RMSE

--- title: "Week 6: Feature Engineering for Linear Regression" description: "Logs, squared terms, categories, interactions, and district fixed effects in R" date: 2026-10-09 --- ## What We Did - Watched the recorded lecture on feature engineering. - Used 2023 ACS five-year data to predict median gross rent for Philadelphia census tracts. - Compared transformations and model specifications using the same ten cross-validation folds. - Reviewed the code and checks needed for Part 5 of the assignment. ## Key Concepts ### 1. Feature engineering and linear models Feature engineering changes the variables supplied to a model. Logs, squared terms, and interactions let linear regression represent relationships that a straight line in the original variables would miss. $$ \hat y = \hat\beta_0 + \hat\beta_1 x + \hat\beta_2 x^2 $$ This is still a linear regression: the model is linear in its coefficients, even though its relationship with $x$ is curved. The lecture evaluates each specification by cross-validation RMSE. A transformation needs to improve predictions on held-out data; a better training fit alone is not enough. ### 2. Build the predictors and keep the sample consistent The baseline uses four predictors: | Variable | Meaning | Construction or unit | |---|---|---| | `income` | Median household income | Dollars | | `pct_renter` | Share of households that rent | `100 * hh_renterE / hh_totalE` | | `bldg_age` | Age based on median year structure built | `2023 - yr_builtE` | | `dist_ch_mi` | Distance from tract geometry to City Hall | Feet divided by 5,280 | In the wide ACS output, `E` identifies the estimate and `M` identifies the margin of error. The lecture renames `rentE` to `rent` and retains `rentM` as `rent_moe`. ```r city_hall <- st_sfc( st_point(c(-75.1636, 39.9526)), crs = 4326 ) %>% st_transform(2272) tracts <- tracts_raw %>% mutate( rent = rentE, rent_moe = rentM, income = incomeE, bldg_age = 2023 - yr_builtE, pct_renter = 100 * hh_renterE / hh_totalE, dist_ch_mi = as.numeric(st_distance(geometry, city_hall)) / 5280 ) %>% select(GEOID, rent, rent_moe, income, bldg_age, pct_renter, dist_ch_mi) %>% filter( !is.na(rent), !is.na(income), !is.na(bldg_age), !is.na(pct_renter) ) ``` `tracts_raw` has already been transformed to EPSG:2272, whose units are feet. The City Hall point starts in EPSG:4326, with coordinates ordered as longitude, latitude, and is transformed to match. `st_distance()` here measures the shortest distance between each tract polygon and the point. It does not measure distance from the tract centroid. A tract containing City Hall has distance zero. The cleaned lecture sample contains 377 tracts. All specifications use these rows. If a later feature introduces missing or invalid values, handle them before comparing models and rebuild the folds if the sample changes. ### 3. Baseline model and fixed-fold cross-validation ```r f0 <- rent ~ income + pct_renter + bldg_age + dist_ch_mi m0 <- lm(f0, data = tracts) ``` `f0` stores the formula. `lm()` estimates the coefficients, and `m0` stores the fitted model. `predict(m0, newdata = ...)` applies those coefficients to another dataset. The lecture assigns each row to one of ten folds once: ```r set.seed(5080) folds <- sample(rep(1:10, length.out = nrow(tracts))) cv_rmse <- function(f, data, log_outcome = FALSE) { preds <- numeric(nrow(data)) for (k in 1:10) { m <- lm(f, data = data[folds != k, ]) p <- predict(m, newdata = data[folds == k, ]) if (log_outcome) p <- exp(p) preds[folds == k] <- p } sqrt(mean((data$rent - preds)^2)) } ``` - `folds != k` selects the nine training folds. - `folds == k` selects the held-out fold. - `preds[folds == k] <- p` puts predictions back into their original row positions. - After ten iterations, every row has one out-of-fold prediction. $$ RMSE = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_i-\hat y_i)^2} $$ This helper computes one RMSE across all out-of-fold predictions. It does not take an unweighted average of ten separate fold RMSEs. RMSE has the same units as the outcome and gives more weight to large errors than mean absolute error does. The function uses the external object `folds`, so every input dataset must preserve the original row count and order. Reordering the data changes which observations belong to each fold. ```r model_log <- tibble( model = "m0", added = "baseline", cv_rmse = cv_rmse(f0, tracts) ) ``` `model_log` is a record of the models tried. The word `log` in this name has nothing to do with a logarithmic transformation. ### 4. Log transformations and coefficient interpretation A log compresses a long right tail. It also makes equal proportional changes in a predictor have equal distances on the transformed scale. In R, `log()` uses the natural logarithm, and `exp()` reverses it. ```r f_logx <- rent ~ log(income) + pct_renter + bldg_age + dist_ch_mi m_logx <- lm(f_logx, data = tracts) ``` For this model, the coefficient on `log(income)` is about 502.18. Holding the other predictors constant, a 1% increase in income corresponds to about $5.02 higher predicted rent. | Model | Formula | Interpretation of $\beta$, holding other variables constant | |---|---|---| | Level-level | `y ~ x` | +1 unit of X corresponds to +$\beta$ units of Y | | Level-log | `y ~ log(x)` | +1% in X corresponds to approximately +$\beta/100$ units of Y | | Log-level | `log(y) ~ x` | +1 unit of X corresponds to approximately $100\beta$% change in Y | | Log-log | `log(y) ~ log(x)` | +1% in X corresponds to approximately $\beta$% change in Y | For level-log models, a proportional increase $r$ in X gives the exact predicted change: $$ \Delta\hat y = \beta\ln(1+r) $$ ```r b <- coef(m_logx)["log(income)"] b / 100 # About $5.02 for a 1% income increase b * log(1.10) # About $47.86 for a 10% income increase ``` For log-level models, the exact percentage change in the exponentiated prediction for a one-unit increase in X is: $$ 100\left(e^\beta-1\right)\% $$ The log-log income coefficient in the lecture is about 0.363. A 1% increase in income corresponds to approximately 0.363% higher predicted rent. This coefficient is an elasticity. #### Zeros and `log1p()` ```r log(0) # -Inf log1p(0) # 0 log1p(x) # log(1 + x) ``` Counts such as crimes or 311 calls can include meaningful zeros. `log1p()` preserves those observations when transforming nonnegative counts. Its interpretation concerns `x + 1`, so the usual percentage interpretation for `log(x)` is less accurate near zero. Neither function makes arbitrary negative values suitable for logging. #### Logging Y changes the prediction scale ```r f_logy <- log(rent) ~ log(income) + pct_renter + bldg_age + dist_ch_mi cv_rmse(f_logy, tracts, log_outcome = TRUE) ``` Predictions from this model are on the log-rent scale. Apply `exp()` to each prediction before comparing it with observed rent in dollars. A log-scale RMSE of 0.188 cannot be compared with a dollar RMSE of 265, and exponentiating the final RMSE does not fix this. Direct exponentiation does not generally recover the conditional mean of rent; it can underpredict that mean. The lecture uses the simple back-transformation for model comparisons and notes that bias corrections exist. A skewed predictor does not violate a requirement that X be normally distributed. Linear regression has no such requirement. Use plots and CV to decide whether a log is useful. ### 5. Squared terms and changing slopes ```r f_sq <- rent ~ log(income) + pct_renter + bldg_age + I(bldg_age^2) + dist_ch_mi m_sq <- lm(f_sq, data = tracts) ``` Keep both the original variable and its square. `I()` tells R to perform arithmetic inside the formula; `^` otherwise has a special formula meaning. For the age component: $$ \hat y = \cdots + \beta_1 age + \beta_2 age^2 $$ The slope depends on age: $$ \frac{d\hat y}{d age} = \beta_1 + 2\beta_2 age $$ The lecture estimates $\beta_1\approx-28.004$ and $\beta_2\approx0.1823$. The positive squared coefficient gives an upward-opening curve. The linear coefficient alone does not mean that rent falls by $28 for every additional year of age throughout the sample. The turning point is: $$ age^* = -\frac{\beta_1}{2\beta_2} $$ ```r b1 <- coef(m_sq)["bldg_age"] b2 <- coef(m_sq)["I(bldg_age^2)"] -b1 / (2 * b2) # About 76.8 years range(tracts$bldg_age) # 5 to 85 years ``` Check whether the turning point lies inside the observed range. A turning point outside that range does not establish an observed reversal. A quadratic also keeps curving outside the training range, so extrapolated predictions can become unreasonable. The old-building pile-up in the lecture plot partly comes from Census coding of early construction years. Check coding limits before interpreting a concentration of observations as a housing-market pattern. ### 6. Categories, dummy variables, and reference groups Use a categorical predictor for existing groups, such as zoning classes, or when a continuous relationship is better represented by distinct groups. ```r tracts <- tracts %>% mutate( ring = case_when( dist_ch_mi < 2 ~ "Core", dist_ch_mi < 5 ~ "Inner", TRUE ~ "Outer" ), ring = factor(ring, levels = c("Core", "Inner", "Outer")) ) ``` `case_when()` uses the first matching condition. The Inner group therefore contains distances from 2 up to, but not including, 5 miles. Core is the reference because it is the first specified factor level. For numeric breakpoints, use `cut()`: ```r tracts <- tracts %>% mutate(renter_band = cut( pct_renter, breaks = c(0, 33, 66, 100), labels = c("Mostly owners", "Mixed", "Mostly renters"), include.lowest = TRUE )) ``` With an intercept and the usual treatment coding, R represents $k$ categories using $k-1$ dummy variables: | Ring | `ringInner` | `ringOuter` | |---|---|---| | Core | 0 | 0 | | Inner | 1 | 0 | | Outer | 0 | 1 | ```r model.matrix(~ ring, data = st_drop_geometry(tracts)) ``` In the lecture model, `ringInner` is about -122 and `ringOuter` is about -164. Holding the other predictors constant, predicted rent in Inner tracts is about $122 lower than in Core tracts, and predicted rent in Outer tracts is about $164 lower than in Core tracts. These coefficients are differences from Core, not rent levels. ```r tracts_r <- tracts %>% mutate(ring = relevel(ring, ref = "Outer")) ``` Changing the reference changes the coefficient comparisons while leaving predictions unchanged. Numeric IDs also need `factor()` when their numbers are labels; treating district codes 1 through 18 as a continuous predictor would impose a numerical slope across arbitrary codes. Binning removes information. Distances of 2.1 and 4.9 miles both become Inner. In this example, replacing continuous distance with rings increased CV RMSE. ### 7. Interactions: different slopes by group An additive category term shifts the intercept. An interaction allows the relationship with another predictor to vary by group. ```r f_int <- rent ~ log(income) * ring + pct_renter + bldg_age + I(bldg_age^2) ``` In an R formula: ```r x * group # Expands to: x + group + x:group ``` Writing $x=\ln(income)$, the relevant terms are: $$ \hat y = \beta_0 + \beta_1x + \gamma_I D_I + \gamma_O D_O + \delta_I xD_I + \delta_O xD_O + \cdots $$ | Ring | Intercept, with other predictors set to zero | Slope on log income | |---|---|---| | Core | $\beta_0$ | $\beta_1$ | | Inner | $\beta_0+\gamma_I$ | $\beta_1+\delta_I$ | | Outer | $\beta_0+\gamma_O$ | $\beta_1+\delta_O$ | The lecture coefficients give these log-income slopes: - Core: 653.70. - Inner: $653.70-219.10=434.60$. - Outer: $653.70-244.42=409.28$. The main log-income coefficient is the Core slope. Each interaction coefficient is a difference from that slope. With an interaction, the Outer-Core rent difference also depends on income: $\gamma_O+\delta_Ox$. Use interactions when the relationship or a plot suggests different slopes. Each additional group slope requires another coefficient estimated from the same sample. ### 8. District fixed effects The lecture uses a district factor to allow separate intercepts for Philadelphia's planning districts: $$ y_{ig} = \alpha_g + \beta x_{ig} + \epsilon_{ig} $$ $\alpha_g$ is the intercept for district $g$. It captures district-level differences left unexplained by the measured predictors, without identifying their individual causes. ```r tracts <- tracts %>% st_join(districts %>% select(district = dist_name), largest = TRUE) %>% mutate(district = factor(district)) f_fe <- rent ~ log(income) + pct_renter + bldg_age + I(bldg_age^2) + district ``` `largest = TRUE` assigns a tract crossing district boundaries to the district with the largest overlapping area. It assigns attributes to the tract rather than splitting it into pieces. After adding district fixed effects, the log-income coefficient changes from about 519.74 to 427.39. The income relationship now uses variation within districts after controlling for their intercepts. Some of the earlier association reflected differences between districts. Fixed effects do not account for every unmeasured difference within a district, and they do not turn an association into a causal effect. Residual maps can show whether geographic clusters remain after adding them. #### Small groups and unseen factor levels ```r tracts %>% st_drop_geometry() %>% count(district) %>% arrange(n) ``` One lecture district has only one tract. When that tract is held out, the training set contains no observation from its district, so prediction encounters a factor level the model has never learned. The lecture handles groups with fewer than ten observations by assigning their tracts to the district of the nearest tract in a sufficiently large group. It keeps the original district label and creates `district_cv` for the CV specification. This changes group membership; it does not dissolve or combine tract geometries. Other options discussed include a coarser geography or pooling small groups into an Other category. A geographically sensible merge gives the resulting category a clearer meaning. About ten observations per group is a classroom rule of thumb for ten-fold CV, not a guarantee that every training fold will contain every category. ### 9. R syntax to remember | Code | Purpose or check | |---|---| | `log(x)` | Natural log; requires positive values for finite results | | `log1p(x)` | Compute `log(1 + x)`; accommodates zero counts | | `exp(pred)` | Back-transform a log prediction | | `x + I(x^2)` | Include a linear term and its square | | `factor(x)` | Treat labels as categories | | `relevel(x, ref = "...")` | Choose a reference category | | `case_when()` / `cut()` | Construct categories from conditions or breakpoints | | `x * group` | Include both main effects and their interaction | | `model.matrix()` | Inspect the encoded predictors | | `count(group)` | Check group sizes before CV | | `st_join(..., largest = TRUE)` | Assign the spatial match with the largest overlap | | `add_row()` | Record another model and its CV RMSE |