---
title: "Lecture 2: Introduction to Prediction - Code Examples"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Code Examples from Lecture 2

## Setup

```{r load-packages}
library(mosaicData)
library(tidyverse)
```

```{r data-prep}
# Read in Saratoga County houses data and select six columns with renaming
sh = SaratogaHouses %>%
  select(price,
         livingArea,
         age,
         bedrooms,
         bathrooms,
         heating,
         new = newConstruction)
```

## Idea 1: The Sample Mean

```{r sample-mean}
# Compute the sample mean
c = mean(sh$price)
c

# Compute MSE (Mean Squared Error)
mse = mean((sh$price - c)^2)
mse

# Compute RMSE (Root Mean Squared Error)
rmse = sqrt(mse)
rmse
```

## Idea 2: Group Means

### Group by Heating Type

```{r group-by-heating}
# Compute sample means for all three heating types
c_hotair = mean(sh$price[sh$heating == "hot air"])
c_hotwater = mean(sh$price[sh$heating == "hot water/steam"])
c_electric = mean(sh$price[sh$heating == "electric"])

# Display group means
cat("c_hotair (mean for hot air heating):", c_hotair, "\n")
cat("c_hotwater (mean for hot water/steam heating):", c_hotwater, "\n")
cat("c_electric (mean for electric heating):", c_electric, "\n")

# Compute MSE and RMSE for the three heating type means
# Predicted values based on heating type
predicted_values = case_when(
  sh$heating == "hot air" ~ c_hotair,
  sh$heating == "hot water/steam" ~ c_hotwater,
  sh$heating == "electric" ~ c_electric
)
mse_grouped = mean((sh$price - predicted_values)^2)
rmse_grouped = sqrt(mse_grouped)

cat("MSE for three heating type means:", mse_grouped, "\n")
cat("RMSE for three heating type means:", rmse_grouped, "\n")
```

### Binning Living Area

```{r bin-living-area}
# Bin livingArea by rounding to nearest multiple of 200
sh$livingArea_bin = round(sh$livingArea / 200) * 200

# Compute sample mean by bin
bin_means = sh %>%
  group_by(livingArea_bin) %>%
  summarize(bin_mean = mean(price), .groups = "drop")

# Merge bin means back to original data
sh = sh %>%
  left_join(bin_means, by = "livingArea_bin")

# Compute MSE and RMSE using bin-based predictions
mse_binned = mean((sh$price - sh$bin_mean)^2)
rmse_binned = sqrt(mse_binned)

cat("MSE for bin means:", mse_binned, "\n")
cat("RMSE for bin means:", rmse_binned, "\n")
```

## Idea 3: Ordinary Least Squares (OLS) Linear Regression

### Example 1: Single Variable (Living Area)

```{r ols-single-variable}
# Simple linear regression against 'livingArea' variable
lm_area = lm(price ~ livingArea, data = sh)
summary(lm_area)

# Compute predictions and error metrics
sh$lm_area_predictions = predict(lm_area, sh)
mse_lm_area = mean((sh$price - sh$lm_area_predictions)^2)
rmse_lm_area = sqrt(mse_lm_area)

cat("MSE for regression on 'livingArea':", mse_lm_area, "\n")
cat("RMSE for regression on 'livingArea':", rmse_lm_area, "\n")

# Plot the model
ggplot(data = sh, aes(x = livingArea, y = price)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)
```

### Example 2: Categorical Variable (Heating)

```{r ols-categorical}
# Simple linear regression against 'heating' variable with intercept
lm_heating = lm(price ~ 1 + heating, data = sh)
summary(lm_heating)

# Compute predictions and error metrics
sh$lm_heating_predictions = predict(lm_heating, sh)
mse_lm_heating = mean((sh$price - sh$lm_heating_predictions)^2)
rmse_lm_heating = sqrt(mse_lm_heating)

cat("MSE for regression on 'heating':", mse_lm_heating, "\n")
cat("RMSE for regression on 'heating':", rmse_lm_heating, "\n")
```

### Example 3: All Features

```{r ols-all-features}
# Linear regression with all features and intercept
lm_model = lm(price ~ 1 + livingArea + age + bedrooms +
                bathrooms + heating + new, data = sh)
summary(lm_model)

# Get predictions from the linear model
sh$lm_predictions = predict(lm_model, sh)

# Compute MSE and RMSE for linear regression
mse_lm = mean((sh$price - sh$lm_predictions)^2)
rmse_lm = sqrt(mse_lm)

cat("MSE for linear regression:", mse_lm, "\n")
cat("RMSE for linear regression:", rmse_lm, "\n")
```

## Idea 4: Feature Engineering

### Higher Order Terms

```{r higher-order-terms}
# Linear regression with livingArea and its second power
lm_area2 = lm(price ~ 1 + livingArea + age + bedrooms +
                bathrooms + heating + new + I(livingArea^2), data = sh)

# Compute predictions and error metrics
sh$lm_area2_predictions = predict(lm_area2, sh)
mse_lm_area2 = mean((sh$price - sh$lm_area2_predictions)^2)
rmse_lm_area2 = sqrt(mse_lm_area2)

cat("MSE for regression with livingArea^2:", mse_lm_area2, "\n")
cat("RMSE for regression with livingArea^2:", rmse_lm_area2, "\n")

# Display R-squared
summary(lm_area2)$r.squared
```

### Even Higher Order Terms

```{r higher-order-234}
# Linear regression with livingArea and powers 2 through 4
lm_area234 = lm(price ~ 1 + livingArea + age + bedrooms +
                  bathrooms + heating + new + I(livingArea^2) +
                  I(livingArea^3) + I(livingArea^4), data = sh)

# Compute predictions and error metrics
sh$lm_area234_predictions = predict(lm_area234, sh)
mse_lm_area234 = mean((sh$price - sh$lm_area234_predictions)^2)
rmse_lm_area234 = sqrt(mse_lm_area234)

cat("MSE for regression with livingArea^2-4:", mse_lm_area234, "\n")
cat("RMSE for regression with livingArea^2-4:", rmse_lm_area234, "\n")

# Display R-squared
summary(lm_area234)$r.squared
```

## Additional Examples from Lecture

### Simple Regression with 'new' Variable

```{r ols-new-variable}
# Simple linear regression against 'new' variable
lm_new = lm(price ~ new, data = sh)
summary(lm_new)

sh$lm_new_predictions = predict(lm_new, sh)
mse_lm_new = mean((sh$price - sh$lm_new_predictions)^2)
rmse_lm_new = sqrt(mse_lm_new)

cat("MSE for regression on 'new':", mse_lm_new, "\n")
cat("RMSE for regression on 'new':", rmse_lm_new, "\n")
```

### Collinearity Example

```{r collinearity-example}
# Create a copy of livingArea to demonstrate collinearity
sh$livingArea_copy = sh$livingArea

# Try to fit a model with perfectly collinear variables
fm_collinear = lm(data = sh, price ~ 1 + livingArea + livingArea_copy)
coef(fm_collinear)
```

## Comparison of Model Performance

```{r model-comparison}
# Create a summary table of all models
model_summary = data.frame(
  Model = c("Sample Mean", "Heating Groups", "Living Area Bins",
            "Linear: Living Area", "Linear: Heating", "Linear: All Features",
            "Linear: + Area^2", "Linear: + Area^2-4"),
  RMSE = c(rmse, rmse_grouped, rmse_binned, rmse_lm_area,
           rmse_lm_heating, rmse_lm, rmse_lm_area2, rmse_lm_area234),
  R_squared = c(0,
                1 - mse_grouped/mse,
                1 - mse_binned/mse,
                summary(lm_area)$r.squared,
                summary(lm_heating)$r.squared,
                summary(lm_model)$r.squared,
                summary(lm_area2)$r.squared,
                summary(lm_area234)$r.squared)
)

print(model_summary)
```