Modeling Ames Home Sale Prices: A Multiple Regression & Bootstrap Analysis

1 Introduction

The goal of this project is to model home sale prices in Ames, Iowa using structural, size, and amenity characteristics from the Ames Housing dataset.

We aim to:

• Build an interpretable multiple linear regression model for sale price.
  
• Diagnose violations of linear-model assumptions and apply a Box–Cox transformation when appropriate.
  
• Compare several candidate models (full, Box–Cox, and stepwise AIC) using goodness-of-fit criteria.
  
• Use bootstrap methods to construct confidence intervals for regression coefficients and evaluate the robustness of parametric inference.

Practical questions include:

• Which features (living area, quality, bathrooms, garages, fireplaces, age) have the strongest impact on sale price?
  
• How large are those impacts in dollar terms?
  
• Are bootstrap confidence intervals consistent with the usual parametric intervals?

2 Data Set and Practical Questions

2.1 Data source

We use the make_ames() function from the AmesHousing package. Each observation corresponds to a single residential home sale.

housing <- AmesHousing::make_ames()
dim(housing)

[1] 2930 81

head(housing[, 1:6])

The response variable is

• Sale_Price: sale price of the home (USD).

2.2 Variables used in the analysis

We select predictors related to size, quality, condition, age, and amenities, and construct a few engineered variables.

dat0 <- housing %>%
  dplyr::select(
    Sale_Price,
    Gr_Liv_Area, Total_Bsmt_SF, First_Flr_SF, Second_Flr_SF,
    Overall_Qual, Overall_Cond,
    Year_Built, Year_Remod_Add,
    Full_Bath, Half_Bath, Bsmt_Full_Bath, Bsmt_Half_Bath,
    Garage_Cars, Garage_Area,
    Fireplaces,
    Lot_Area, Mas_Vnr_Area
  ) %>%
  tidyr::drop_na(Sale_Price)

dat <- dat0 %>%
  dplyr::mutate(
    Total_Baths = Full_Bath + 0.5 * Half_Bath +
      Bsmt_Full_Bath + 0.5 * Bsmt_Half_Bath,
    House_Age         = max(Year_Built, na.rm = TRUE) - Year_Built,
    Years_Since_Remod = max(Year_Remod_Add, na.rm = TRUE) - Year_Remod_Add,
    Overall_Qual_S    = as.numeric(Overall_Qual),
    Overall_Cond_S    = as.numeric(Overall_Cond)
  ) %>%
  dplyr::select(
    Sale_Price,
    Gr_Liv_Area, Total_Bsmt_SF, First_Flr_SF, Second_Flr_SF,
    Overall_Qual_S, Overall_Cond_S,
    House_Age, Years_Since_Remod,
    Total_Baths, Garage_Cars, Garage_Area, Fireplaces,
    Lot_Area, Mas_Vnr_Area
  )

dim(dat)

[1] 2930 15

2.2.1 Variable descriptions

• Sale_Price (numeric): sale price of the home (USD).
  
• Gr_Liv_Area (numeric): above-ground living area (sq ft).
  
• Total_Bsmt_SF (numeric): total basement area (sq ft).
  
• First_Flr_SF, Second_Flr_SF (numeric): first- and second-floor areas.
  
• Overall_Qual_S (integer): overall material and finish quality (higher = better).
  
• Overall_Cond_S (integer): overall condition rating.
  
• House_Age (numeric): age of the house at time of sale.
  
• Years_Since_Remod (numeric): years since last remodeling.
  
• Total_Baths (numeric): effective number of bathrooms
  (full = 1, half = 0.5, summed across basement and basement baths).
  
• Garage_Cars (integer): garage car capacity.
  
• Garage_Area (numeric): garage size in sq ft.
  
• Fireplaces (integer): number of fireplaces.
  
• Lot_Area (numeric): lot size in sq ft.
  
• Mas_Vnr_Area (numeric): masonry veneer area (sq ft).

2.2.2 Practical questions

1. Which structural and amenity variables have the strongest impact on sale price?
   
2. Can a linear regression model explain a substantial portion of the variation in sale price while meeting key assumptions?
   
3. Are conclusions based on parametric inference stable under bootstrap resampling?

3 Exploratory Data Analysis

Min. 1st Qu. Median Mean 3rd Qu. Max. 12789 129500 160000 180796 213500 755000

The median home sells for around the mid $160,000s, with substantial variation across the market.

We examine relationships among sale price and several key predictors using a manageable subsample.

set.seed(321)
small <- dat %>% dplyr::sample_n(min(1000, nrow(dat)))

GGally::ggpairs(
  dplyr::select(small,
                Sale_Price, Gr_Liv_Area, Total_Bsmt_SF,
                Overall_Qual_S, Total_Baths, Garage_Cars)
)

Interpretation.
Sale price is strongly and positively associated with living area, basement size, overall quality, and total bathrooms. The relationships appear roughly linear, with increasing spread for larger, more expensive homes, suggesting mild heteroscedasticity.

4 Modeling Strategy

We start with the following full multiple linear regression model:

full_formula <- Sale_Price ~ Gr_Liv_Area + Total_Bsmt_SF +
  First_Flr_SF + Second_Flr_SF +
  Overall_Qual_S + Overall_Cond_S +
  House_Age + Years_Since_Remod +
  Total_Baths + Garage_Cars + Garage_Area +
  Fireplaces + Lot_Area + Mas_Vnr_Area

Plan:

1. Fit the full model and check multicollinearity and residual diagnostics.
   
2. Fit a Box–Cox transformed model to improve normality and variance stability.
   
3. Use stepwise AIC to obtain a parsimonious model.
   
4. Compare models with \(R^2\), adjusted \(R^2\), Mallows’ \(C_p\), and PRESS.
   
5. Select a final model for interpretation and bootstrap inference.

4.1 Full multiple regression model

m_full <- lm(full_formula, data = dat)
summary(m_full)

Call: lm(formula = full_formula, data = dat)

Residuals: Min 1Q Median 3Q Max -518019 -17759 -2451 14208 296657

Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.770e+04 6.299e+03 -12.335 < 2e-16 Gr_Liv_Area 2.189e+01 1.403e+01 1.560 0.1189
Total_Bsmt_SF 1.949e+01 2.639e+00 7.385 1.97e-13 First_Flr_SF 2.543e+01 1.440e+01 1.767 0.0774 .
Second_Flr_SF 1.443e+01 1.420e+01 1.016 0.3097
Overall_Qual_S 1.817e+04 7.712e+02 23.559 < 2e-16 Overall_Cond_S 4.595e+03 6.830e+02 6.728 2.06e-11 House_Age -2.637e+02 3.736e+01 -7.059 2.08e-12 Years_Since_Remod -2.433e+02 4.472e+01 -5.442 5.72e-08 Total_Baths 7.264e+03 1.168e+03 6.220 5.67e-10 Garage_Cars 2.841e+03 1.982e+03 1.434 0.1518
Garage_Area 3.202e+01 6.810e+00 4.702 2.70e-06 Fireplaces 8.071e+03 1.174e+03 6.874 7.63e-12 Lot_Area 5.961e-01 8.903e-02 6.696 2.56e-11 Mas_Vnr_Area 3.817e+01 4.181e+00 9.128 < 2e-16 *** — Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘.’ 0.1 ’ ’ 1

Residual standard error: 34630 on 2915 degrees of freedom Multiple R-squared: 0.813, Adjusted R-squared: 0.8121 F-statistic: 905.3 on 14 and 2915 DF, p-value: < 2.2e-16

The full model explains over 80% of the variation in sale price, with many predictors (especially quality, baths, and fireplaces) highly significant.

4.1.1 Multicollinearity

vif_vals <- car::vif(m_full)
sort(round(vif_vals, 2), decreasing = TRUE)

  Gr_Liv_Area     Second_Flr_SF      First_Flr_SF       Garage_Cars 
       122.87             90.44             77.75              5.56 
  Garage_Area     Total_Bsmt_SF         House_Age    Overall_Qual_S 
         5.25              3.31              3.12              2.89 
  Total_Baths Years_Since_Remod    Overall_Cond_S        Fireplaces 
         2.17              2.13              1.41              1.41 
 Mas_Vnr_Area          Lot_Area 
         1.36              1.20 

Several size-related variables (Gr_Liv_Area, First_Flr_SF, Second_Flr_SF) have large VIFs, indicating notable multicollinearity.

4.1.2 Residual diagnostics

op <- par(mfrow = c(2, 2))
plot(m_full)

par(op)

Residual plots show mild skewness and hints of non-constant variance; Q–Q plots indicate slightly heavy upper tails.

5 Box–Cox Transformation

y <- dat$Sale_Price
eps <- ifelse(any(y <= 0), 1, 0)

bc <- MASS::boxcox(m_full, lambda = seq(-1, 1, by = 0.1), plotit = TRUE)

(lambda_hat <- bc$x[which.max(bc$y)])

[1] 0.1515152

The estimated lambda is around 0.15, suggesting a mild power transformation.

# Recompute Box–Cox inside this chunk so lambda_hat always exists
y <- dat$Sale_Price
eps <- ifelse(any(y <= 0), 1, 0)

# Box–Cox profile based on full model
bc <- MASS::boxcox(m_full, lambda = seq(-1, 1, by = 0.1))

lambda_hat <- bc$x[which.max(bc$y)]
lambda_hat

[1] 0.1515152

# Transform the response using lambda_hat
y_bc <- if (abs(lambda_hat) < 1e-8) {
  log(y + eps)
} else {
  ((y + eps)^lambda_hat - 1) / lambda_hat
}

# Fit Box–Cox transformed model
m_bc <- lm(
  y_bc ~ Gr_Liv_Area + Total_Bsmt_SF + First_Flr_SF + Second_Flr_SF +
    Overall_Qual_S + Overall_Cond_S + House_Age + Years_Since_Remod +
    Total_Baths + Garage_Cars + Garage_Area + Fireplaces +
    Lot_Area + Mas_Vnr_Area,
  data = dat
)

summary(m_bc)

Call: lm(formula = y_bc ~ Gr_Liv_Area + Total_Bsmt_SF + First_Flr_SF + Second_Flr_SF + Overall_Qual_S + Overall_Cond_S + House_Age + Years_Since_Remod + Total_Baths + Garage_Cars + Garage_Area + Fireplaces + Lot_Area + Mas_Vnr_Area, data = dat)

Residuals: Min 1Q Median 3Q Max -15.1262 -0.4469 0.0052 0.4908 3.2168

Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.549e+01 1.681e-01 151.637 < 2e-16 Gr_Liv_Area 8.797e-04 3.744e-04 2.349 0.018874
Total_Bsmt_SF 6.307e-04 7.044e-05 8.954 < 2e-16 First_Flr_SF 5.040e-04 3.842e-04 1.312 0.189683
Second_Flr_SF 2.632e-04 3.791e-04 0.694 0.487479
Overall_Qual_S 5.684e-01 2.058e-02 27.614 < 2e-16 Overall_Cond_S 3.081e-01 1.823e-02 16.904 < 2e-16 House_Age -1.622e-02 9.971e-04 -16.267 < 2e-16 Years_Since_Remod -6.867e-03 1.193e-03 -5.754 9.62e-09 Total_Baths 2.759e-01 3.117e-02 8.854 < 2e-16 Garage_Cars 2.375e-01 5.289e-02 4.490 7.40e-06 Garage_Area 6.377e-04 1.817e-04 3.509 0.000457 Fireplaces 3.425e-01 3.134e-02 10.929 < 2e-16 Lot_Area 1.835e-05 2.376e-06 7.722 1.56e-14 ** Mas_Vnr_Area 1.986e-04 1.116e-04 1.779 0.075309 .
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘.’ 0.1 ’ ’ 1

Residual standard error: 0.9242 on 2915 degrees of freedom Multiple R-squared: 0.8671, Adjusted R-squared: 0.8664 F-statistic: 1358 on 14 and 2915 DF, p-value: < 2.2e-16

op <- par(mfrow = c(2, 2))
plot(m_bc)

par(op)

The Box–Cox model improves residual normality and slightly increases adjusted \(R^2\), although interpretation becomes less direct due to the transformed response.

6 Model Selection and Goodness-of-Fit

6.1 Stepwise AIC model

m_step <- MASS::stepAIC(m_full, trace = FALSE)
summary(m_step)

Call: lm(formula = Sale_Price ~ Gr_Liv_Area + Total_Bsmt_SF + First_Flr_SF + Overall_Qual_S + Overall_Cond_S + House_Age + Years_Since_Remod + Total_Baths + Garage_Cars + Garage_Area + Fireplaces + Lot_Area + Mas_Vnr_Area, data = dat)

Residuals: Min 1Q Median 3Q Max -517086 -17863 -2425 14259 297387

Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.772e+04 6.299e+03 -12.338 < 2e-16 Gr_Liv_Area 3.599e+01 2.063e+00 17.443 < 2e-16 Total_Bsmt_SF 1.938e+01 2.637e+00 7.350 2.56e-13 First_Flr_SF 1.114e+01 3.063e+00 3.639 0.000279 Overall_Qual_S 1.818e+04 7.711e+02 23.576 < 2e-16 Overall_Cond_S 4.626e+03 6.823e+02 6.781 1.44e-11 House_Age -2.670e+02 3.723e+01 -7.172 9.34e-13 Years_Since_Remod -2.423e+02 4.471e+01 -5.420 6.46e-08 Total_Baths 7.324e+03 1.166e+03 6.280 3.89e-10 Garage_Cars 2.887e+03 1.981e+03 1.457 0.145096
Garage_Area 3.199e+01 6.810e+00 4.698 2.75e-06 Fireplaces 8.092e+03 1.174e+03 6.893 6.67e-12 Lot_Area 5.982e-01 8.900e-02 6.721 2.16e-11 Mas_Vnr_Area 3.841e+01 4.175e+00 9.200 < 2e-16 *** — Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘.’ 0.1 ’ ’ 1

Residual standard error: 34630 on 2916 degrees of freedom Multiple R-squared: 0.8129, Adjusted R-squared: 0.8121 F-statistic: 974.8 on 13 and 2916 DF, p-value: < 2.2e-16

round(car::vif(m_step), 2)

  Gr_Liv_Area     Total_Bsmt_SF      First_Flr_SF    Overall_Qual_S 
         2.66              3.30              3.52              2.89 

Overall_Cond_S House_Age Years_Since_Remod Total_Baths 1.40 3.10 2.12 2.17 Garage_Cars Garage_Area Fireplaces Lot_Area 5.55 5.25 1.41 1.20 Mas_Vnr_Area 1.36

The stepwise procedure removes some collinear or weak predictors while maintaining strong fit and reducing VIF values.

6.2 Goodness-of-fit comparison

Goodness-of-fit summary for candidate models.

Name	SSE	R2	R2_Adj	Cp	PRESS	p
Box–Cox	2490	0.867	0.866	-2900	2589	15
Stepwise (orig)	3.497e+12	0.813	0.812	29.04	3.626e+12	14
Full	3.495e+12	0.813	0.812	30	3.629e+12	15

Interpretation.

• The Box–Cox model has the best raw \(R^2\), but uses a transformed response.
  
• The stepwise model achieves similar adjusted \(R^2\) and PRESS with fewer predictors and lower multicollinearity.

We select the stepwise original model as the final model due to its balance of fit and interpretability.

7 Final Model and Interpretation

m_final <- m_step
summary(m_final)

Call: lm(formula = Sale_Price ~ Gr_Liv_Area + Total_Bsmt_SF + First_Flr_SF + Overall_Qual_S + Overall_Cond_S + House_Age + Years_Since_Remod + Total_Baths + Garage_Cars + Garage_Area + Fireplaces + Lot_Area + Mas_Vnr_Area, data = dat)

Residuals: Min 1Q Median 3Q Max -517086 -17863 -2425 14259 297387

Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.772e+04 6.299e+03 -12.338 < 2e-16 Gr_Liv_Area 3.599e+01 2.063e+00 17.443 < 2e-16 Total_Bsmt_SF 1.938e+01 2.637e+00 7.350 2.56e-13 First_Flr_SF 1.114e+01 3.063e+00 3.639 0.000279 Overall_Qual_S 1.818e+04 7.711e+02 23.576 < 2e-16 Overall_Cond_S 4.626e+03 6.823e+02 6.781 1.44e-11 House_Age -2.670e+02 3.723e+01 -7.172 9.34e-13 Years_Since_Remod -2.423e+02 4.471e+01 -5.420 6.46e-08 Total_Baths 7.324e+03 1.166e+03 6.280 3.89e-10 Garage_Cars 2.887e+03 1.981e+03 1.457 0.145096
Garage_Area 3.199e+01 6.810e+00 4.698 2.75e-06 Fireplaces 8.092e+03 1.174e+03 6.893 6.67e-12 Lot_Area 5.982e-01 8.900e-02 6.721 2.16e-11 Mas_Vnr_Area 3.841e+01 4.175e+00 9.200 < 2e-16 *** — Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘.’ 0.1 ’ ’ 1

Residual standard error: 34630 on 2916 degrees of freedom Multiple R-squared: 0.8129, Adjusted R-squared: 0.8121 F-statistic: 974.8 on 13 and 2916 DF, p-value: < 2.2e-16

confint(m_final)

                      2.5 %        97.5 %

(Intercept) -9.007652e+04 -6.537308e+04 Gr_Liv_Area 3.194066e+01 4.003102e+01 Total_Bsmt_SF 1.421075e+01 2.455175e+01 First_Flr_SF 5.138295e+00 1.714841e+01 Overall_Qual_S 1.666798e+04 1.969192e+04 Overall_Cond_S 3.288468e+03 5.964038e+03 House_Age -3.399672e+02 -1.939860e+02 Years_Since_Remod -3.299438e+02 -1.546299e+02 Total_Baths 5.037435e+03 9.610963e+03 Garage_Cars -9.971247e+02 6.771933e+03 Garage_Area 1.863786e+01 4.534281e+01 Fireplaces 5.790373e+03 1.039430e+04 Lot_Area 4.236867e-01 7.727211e-01 Mas_Vnr_Area 3.022111e+01 4.659273e+01

Parametric estimates and confidence intervals for the final model. (continued below)

term	estimate	std.error	statistic	p.value	conf.low
(Intercept)	-77725	6299	-12.34	0	-90077
Gr_Liv_Area	35.99	2.063	17.44	0	31.94
Total_Bsmt_SF	19.38	2.637	7.35	0	14.21
Overall_Qual_S	18180	771.1	23.58	0	16668
Overall_Cond_S	4626	682.3	6.781	0	3288
House_Age	-267	37.23	-7.172	0	-340
Years_Since_Remod	-242.3	44.71	-5.42	0	-329.9
Total_Baths	7324	1166	6.28	0	5037
Garage_Area	31.99	6.81	4.698	0	18.64
Fireplaces	8092	1174	6.893	0	5790
Lot_Area	0.5982	0.089	6.721	0	0.4237
Mas_Vnr_Area	38.41	4.175	9.2	0	30.22
First_Flr_SF	11.14	3.063	3.639	3e-04	5.138
Garage_Cars	2887	1981	1.458	0.1451	-997.1

conf.high
-65373
40.03
24.55
19692
5964
-194
-154.6
9611
45.34
10394
0.7727
46.59
17.15
6772

Key findings:

• Overall quality has one of the largest positive effects; each one-unit increase in quality rating is associated with a sizable increase in sale price, holding other variables fixed.
  
• Total bathrooms and garage capacity/area have strong positive coefficients, indicating that additional bathrooms and parking capacity yield meaningful price premiums.
  
• Fireplaces add value, though on a smaller scale than quality and size.
  
• House age and years since remodel generally have negative coefficients, meaning older and less recently updated homes sell for less.

8 8 Bootstrap Regression

We use bootstrap methods to assess the robustness of parametric inference for the final model.

8.1 Bootstrapping cases

set.seed(321)

X <- model.matrix(m_final)
y_final <- model.response(model.frame(m_final))
p <- length(coef(m_final))
coef_names <- names(coef(m_final))
B <- 1000

boot_beta_case <- matrix(NA_real_, nrow = B, ncol = p)
colnames(boot_beta_case) <- coef_names

for (b in 1:B) {
  idx <- sample.int(nrow(X), replace = TRUE)
  fit_b <- lm(y_final[idx] ~ X[idx, -1])
  boot_beta_case[b, ] <- coef(fit_b)
}

ci_case_mat <- apply(boot_beta_case, 2,
                     quantile, probs = c(0.025, 0.975), na.rm = TRUE)

ci_case_tbl <- tibble::tibble(
  term  = colnames(ci_case_mat),
  `2.5%` = as.numeric(ci_case_mat[1, ]),
  `97.5%`= as.numeric(ci_case_mat[2, ])
)

pander::pander(ci_case_tbl,
               caption = "95% bootstrap confidence intervals (case resampling).")

95% bootstrap confidence intervals (case resampling).

term	2.5%	97.5%
(Intercept)	-93979	-62811
Gr_Liv_Area	25.08	46.42
Total_Bsmt_SF	4.742	31.93
First_Flr_SF	2.983	19.83
Overall_Qual_S	15904	20600
Overall_Cond_S	3282	6082
House_Age	-340.4	-186.2
Years_Since_Remod	-309	-176.6
Total_Baths	4413	10224
Garage_Cars	-4657	11975
Garage_Area	5.562	54.25
Fireplaces	5292	10937
Lot_Area	0.363	1.063
Mas_Vnr_Area	23.79	53.63

8.2 Bootstrapping residuals

fit_hat   <- fitted(m_final)
resid_hat <- resid(m_final)

boot_beta_res <- matrix(NA_real_, nrow = B, ncol = p)
colnames(boot_beta_res) <- coef_names

for (b in 1:B) {
  y_star <- fit_hat + sample(resid_hat, replace = TRUE)
  fit_b  <- lm(y_star ~ X[, -1])
  boot_beta_res[b, ] <- coef(fit_b)
}

ci_resid_mat <- apply(boot_beta_res, 2,
                      quantile, probs = c(0.025, 0.975), na.rm = TRUE)

ci_resid_tbl <- tibble::tibble(
  term  = colnames(ci_resid_mat),
  `2.5%` = as.numeric(ci_resid_mat[1, ]),
  `97.5%`= as.numeric(ci_resid_mat[2, ])
)

pander::pander(ci_resid_tbl,
               caption = "95% bootstrap confidence intervals (residual resampling).")

95% bootstrap confidence intervals (residual resampling).

term	2.5%	97.5%
(Intercept)	-90205	-65312
Gr_Liv_Area	32.02	40
Total_Bsmt_SF	14.1	24.98
First_Flr_SF	4.273	17.05
Overall_Qual_S	16719	19599
Overall_Cond_S	3271	5969
House_Age	-336	-199.7
Years_Since_Remod	-330.9	-149.2
Total_Baths	5015	9557
Garage_Cars	-1096	6664
Garage_Area	18.6	45.32
Fireplaces	5949	10589
Lot_Area	0.4408	0.7659
Mas_Vnr_Area	29.71	46.29

8.3 Comparison with parametric intervals

Across major predictors, the bootstrap intervals from both methods are very similar to the parametric confidence intervals, and all three procedures agree on which coefficients are significantly different from zero. This indicates that the sample size is large enough for classical linear-model inference to be reliable, even in the presence of mild deviations from ideal assumptions.

9 Conclusions

The final stepwise regression model explains a large portion of the variation in Ames home sale prices and highlights that buyers pay the highest premiums for:

• Higher overall quality and condition
  
• Greater living area and bathroom counts
  
• Larger garages and additional amenities such as fireplaces

Older homes and those with long intervals since the last remodel tend to sell for lower prices, holding other factors constant.

Bootstrap inference corroborates the parametric results, suggesting that the key conclusions about drivers of sale price are robust.

10 Limitations and Future Work

• Important factors such as neighborhood school quality, proximity to workplaces, and broader economic trends are not included.
  
• Nonlinearities and interaction effects were not fully explored.
  
• The analysis is specific to Ames and may not generalize to other housing markets.

Future work could incorporate additional location-based variables, explore flexible nonlinear models (splines, tree-based methods), and validate model performance using out-of-sample predictions or time-based cross-validation.