linear models.R

set.seed(2626)

# Create toy dataset

x1 <- rnorm(200, mean = 10, sd =2)
x2 <- rnorm(200, mean = 50, sd =7)
y <- 10 + 4*x1 + 0.5*x2 - 0.1*x1*x2 + rnorm(200, mean = 0, sd = 5)
data <- data.frame(x1, x2, y)

# Spliting the data

train <- sample(1:nrow(data), 100)
training_set <- data[train,]
testing_set <- data[-train,]

# Fitting linear models

lm.1 <- lm(y ~ x1 + x2, data = training_set)
summary(lm.1)

lm.2 <- lm(y ~ x1 + x2 + x1:x2, data = training_set)
summary(lm.2)

lm.3 <- lm(y ~ x1*x2 + I(x1^2) + I(x2^2), data = training_set)
summary(lm.3)

loess.1 <- loess(y ~ x1*x2, span = 0.5, degree = 2, data = training_set)
summary(loess.1)

# Retrieving in-sample RMSE of each model

RMSE <- function(residuals) sqrt(mean(residuals^2))

rmse.lm.1 <- RMSE(lm.1$residuals)
rmse.lm.2 <- RMSE(lm.2$residuals)
rmse.lm.3 <- RMSE(lm.3$residuals)
rmse.loess.1 <- RMSE(loess.1$residuals)

# Calculating the out-of-sample RMSE of each model

test_RMSE <- function(model) RMSE(testing_set$y - predict(model, data = testing_set))
test_rmse.lm.1 <- test_RMSE(lm.1)
test_rmse.lm.2 <- test_RMSE(lm.2)
test_rmse.lm.3 <- test_RMSE(lm.3)
test_rmse.loess.1 <- test_RMSE(loess.1)

training_rmse <- c(rmse.lm.1, rmse.lm.2, rmse.lm.3, rmse.loess.1)
test_rmse <- c(test_rmse.lm.1, test_rmse.lm.2, test_rmse.lm.3, test_rmse.loess.1)
data.frame(training_rmse, test_rmse)