# combine test & training
min_test_id <- min(test_csv$Id)
test_csv$SalePrice <- NA
all_data <- rbind(train_csv, test_csv)

# mutate variables in train and test data sets
current_year <- as.integer(format(Sys.Date(), "%Y"))

all_data_cleaned <- all_data %>% # change the actual year values to age in unit years
mutate(years_since_built = current_year - YearBuilt, years_since_remodel = current_year - 
    YearRemodAdd, years_since_garage_built = current_year - as.integer(GarageYrBlt)) %>% 
    # Per the data dictionary, the following variables have meaningful NA values
mutate_at(vars(Alley, BsmtQual, BsmtCond, BsmtExposure, BsmtFinType1, BsmtFinType2, 
    FireplaceQu, GarageType, GarageFinish, GarageQual, GarageCond, PoolQC, Fence, 
    MiscFeature), fct_explicit_na, na_level = "None") %>% dplyr::select(-c(YearBuilt, 
    YearRemodAdd, GarageYrBlt))

# Remove 2 outliers due to being partial sales of homes
# https://ww2.amstat.org/publications/jse/v19n3/decock/DataDocumentation.txt
plot(SalePrice ~ GrLivArea, all_data_cleaned)

rows_to_remove <- which(all_data_cleaned$GrLivArea > 4000 & all_data_cleaned$SalePrice < 
    2e+05)
all_data_cleaned <- all_data_cleaned[-rows_to_remove, ]


inspect_metadata <- function(df) {
    ### Takes a data frame & Checks for NAs, incorrect class types, inspects the
    ### unique values & adds summary metrics
    df_len <- nrow(df)
    NA_ct = as.vector(rapply(df, function(x) sum(is.na(x))))
    var_names <- names(df)
    # create dataframe
    df_metadata <- data.frame(position = 1:length(var_names), var_names = var_names, 
        NA_ct = NA_ct, NA_percent = round(NA_ct/df_len, 4), class_type = rapply(df, 
            class), uniq_value_ct = rapply(df, function(x) length(unique(x))), 
        uniq_values = rapply(df, function(x) paste(sort(unique(x)), collapse = ",")))
    metrics <- data.frame(describe(df))
    factor_indices <- which(df_metadata$class_type == "factor")
    metrics[factor_indices, ] <- ""
    metrics <- metrics[, 2:ncol(metrics)]
    df_summary <- cbind(df_metadata, metrics)
    return(df_summary)
}

# View(df_metadata <- inspect_metadata(imputed_Data$ximp))

x_data <- subset(all_data_cleaned, select = -c(Id, SalePrice))

imputed_Data <- missForest(x_data)

# check imputed values
sum(is.na(imputed_Data$ximp))
dim(imputed_Data$ximp)
# check imputation error
imputed_Data$OOBerror

all_data_completed <- cbind(Id = all_data_cleaned$Id, imputed_Data$ximp, SalePrice = all_data_cleaned$SalePrice)

# Split test & training
train_df <- subset(all_data_completed, Id < min_test_id, select = -Id)
test_df <- subset(all_data_completed, Id >= min_test_id, select = -c(Id, SalePrice))

baseline <- lm(SalePrice ~ ., data = train_df)

par(mfrow = c(2, 2))
plot(baseline)

model_summary <- function(model, y_var) {
    ### Summarizes the model's key statistics in one row
    df_summary <- glance(summary(model))
    df_summary$n_vars <- ncol(model$model) - 1
    df_summary$RootMSLE = RMSLE(model$fitted.values, model$model[, y_var])
    df_summary$model_name <- deparse(substitute(model))
    df_summary <- subset(df_summary, select = c(model_name, n_vars, adj.r.squared, 
        statistic, RootMSLE, p.value, sigma, df))
}
mod_sum <- model_summary(baseline, "SalePrice")
kable(all_results <- mod_sum)

train_df_transformed <- baseline$model

PT <- powerTransform(as.formula(baseline$call), family = "bcnPower", data = train_df_transformed)

train_df_transformed$SalePrice <- bcnPower(train_df_transformed$SalePrice, lambda = PT$lambda, 
    gamma = PT$gamma)

baseline_PT <- lm(SalePrice ~ ., data = train_df_transformed)

par(mfrow = c(2, 2))
plot(baseline_PT)

mod_sum <- model_summary(baseline_PT, "SalePrice")

kable(all_results <- rbind(all_results, mod_sum))

# if(!exists('AIC_model'))
AIC_model <- stepAIC(baseline)

par(mfrow = c(2, 2))
plot(AIC_model)

sresid <- studres(AIC_model)
hist(sresid, freq = FALSE, breaks = 30, main = "AIC_model of Studentized Residuals")
xfit <- seq(min(na.omit(sresid)), max(na.omit(sresid)), by = 0.1)
yfit <- dnorm(xfit)
lines(xfit, yfit)

set.seed(40)
durbinWatsonTest(AIC_model)

ncvTest(AIC_model)

mod_sum <- model_summary(AIC_model, "SalePrice")
kable(all_results <- rbind(all_results, mod_sum))

train_df_AIC_PT <- AIC_model$model
PT <- powerTransform(as.formula(AIC_model$call), family = "bcnPower", data = train_df_AIC_PT)
summary(PT)

train_df_AIC_PT$SalePrice <- bcnPower(train_df_AIC_PT$SalePrice, lambda = PT$lambda, 
    gamma = PT$gamma)

AIC_model_PT <- lm(SalePrice ~ ., data = train_df_AIC_PT)
# if(!exists('AIC_model_PT'))
AIC_model_PT <- stepAIC(AIC_model_PT)

par(mfrow = c(2, 2))
plot(AIC_model_PT)

sresid <- studres(AIC_model_PT)
hist(sresid, freq = FALSE, breaks = 20, main = "AIC_model_PT of Studentized Residuals")
xfit <- seq(min(na.omit(sresid)), max(na.omit(sresid)), by = 0.1)
yfit <- dnorm(xfit)
lines(xfit, yfit)

set.seed(40)
durbinWatsonTest(AIC_model_PT)

ncvTest(AIC_model_PT)

mod_sum <- model_summary(AIC_model_PT, "SalePrice")
kable(all_results <- rbind(all_results, mod_sum))

significant_predictors <- intersect(names(AIC_model$model), names(test_df))

AIC_model_fit <- predict(AIC_model, newdata = subset(test_df, select = significant_predictors))

kable(t(round(describe(AIC_model_fit), 2)))

AIC_model_fit_output <- data.frame(Id = test_csv$Id, SalePrice = AIC_model_fit)

write.csv(AIC_model_fit_output, "AIC_model_fit.csv", row.names = FALSE)

y_train <- train_df$SalePrice
x_train_df <- subset(train_df, select = -SalePrice)

(rf_fit <- randomForest(x_train_df, y = y_train, ntree = 500, importance = T))

plot(rf_fit)

rf_model_fit <- predict(rf_fit, test_df)

# Summary of Predicted Output
kable(t(round(describe(rf_model_fit), 2)))

rf_model_fit_output <- data.frame(Id = test_csv$Id, SalePrice = rf_model_fit)

write.csv(rf_model_fit_output, "rf_model_fit.csv", row.names = FALSE)