IS605 Final

# find percent missing
pMiss <- function(x){sum(is.na(x))/length(x)*100}

missing_train <- data.frame(cbind(apply(train,2,pMiss)))
missing_train$Feature <- rownames(missing_train)
colnames(missing_train)[1] <- 'Value'

missing_test <- data.frame(cbind(apply(test,2,pMiss)))
missing_test$Feature <- rownames(missing_test)
colnames(missing_test)[1] <- 'Value'

missing <- cbind(missing_test[,1], missing_train[-81,1])
rownames(missing) <- rownames(missing_test)
colnames(missing) <- c('Test', 'Train')
missing <- data.frame(missing)

pander(missing[with(missing, order(-Train)), ][1:5,], caption='Most missing data (Pct)')

drops <- names(train) %in% c('PoolQC', 'MiscFeature', 'Alley', 'Fence')
train <- train[!drops]
drops <- names(test) %in% c('PoolQC', 'MiscFeature', 'Alley', 'Fence')
test <- test[!drops]

# get names of categorical, numeric variables
numerics <- colnames(train[, sapply(train, is.numeric)])
categoricals <- colnames(train[, sapply(train, is.factor)])
categoricals <- categoricals[-1]

corPlot(train[,numerics[1:17]], main='Correlation Plot I')

corPlot(train[,numerics[18:34]], main='Correlation Plot II')

# training variables to drop 
drops <- names(train) %in% c('LotFrontage', 'MasVrnArea', 'GarageYrBlt', 'BedroomAbvGr', 'KitchenAbvGr', 'EnclosedPorch', '3SsnPorch', 'ScreenPorch', 'PoolArea MiscVal', 'YrSold')

train <- train[!drops]

# test variables to drop
drops <- names(test) %in% c('LotFrontage', 'MasVrnArea', 'GarageYrBlt', 'BedroomAbvGr', 'KitchenAbvGr', 'EnclosedPorch', '3SsnPorch', 'ScreenPorch', 'PoolArea MiscVal', 'YrSold')

test <- test[!drops]

# plots for categorical variables

par(mfrow=c(6,2))
for(i in 31:42) {
  plot(train[, categoricals[i]], main=categoricals[i], col='lightblue')
}

# training variables to drop 
drops <- names(train) %in% c("Street", "LotShape","LandContour", "Utilities", "LotConfig", "LandSlope", "Condition1", "Condition2", "OverallCond", "RoofMatl", "Exterior1st", "Exterior2nd", "MasVnrType", "ExterQual", "ExterCond", "Foundation", "BsmtQual", "BsmtCond", "BsmtExposure", "BsmtFinType1", "BsmtFinType2", "Heating", "HeatingQC", "Electrical", "Functional", "FireplaceQu", "GarageFinish",  "GarageQual", "GarageCond",  "PavedDrive", "SaleType", "SaleCondition", "KitchenQual")

train <- train[!drops]

# test variables to drop
drops <- names(test) %in% c("Street", "LotShape","LandContour", "Utilities", "LotConfig", "LandSlope", "Condition1", "Condition2", "OverallCond", "RoofMatl", "Exterior1st", "Exterior2nd", "MasVnrType", "ExterQual", "ExterCond", "Foundation", "BsmtQual", "BsmtCond", "BsmtExposure", "BsmtFinType1", "BsmtFinType2", "Heating", "HeatingQC", "Electrical", "Functional", "FireplaceQu", "GarageFinish",  "GarageQual", "GarageCond",  "PavedDrive", "SaleType", "SaleCondition", "KitchenQual")

test <- test[!drops]

# add variable for total square feet of living space
train$TotSqFt <- train$BsmtFinSF1 + train$BsmtFinSF2 + train$`1stFlrSF` + train$`2ndFlrSF`
test$TotSqFt <- test$BsmtFinSF1 + test$BsmtFinSF2 + test$`1stFlrSF` + test$`2ndFlrSF`

# drop the components and others  
drops <- names(train) %in% c("1stFlrSF", "2ndFlrSF","MSSubClass", "RoofStyle","MasVnrArea","BsmtFinSF1", "BsmtFinSF2", "BsmtUnfSF", "TotalBsmtSF","LowQualFinSF", "GrLivArea", "BsmtFullBath", "BsmtHalfBath", "GarageCars", "WoodDeckSF", "MiscVal", "PoolArea")
                             
train <- train[!drops]

# test variables to drop
drops <- names(test) %in% c("1stFlrSF", "2ndFlrSF","MSSubClass", "RoofStyle","MasVnrArea","BsmtFinSF1", "BsmtFinSF2", "BsmtUnfSF", "TotalBsmtSF","LowQualFinSF", "GrLivArea", "BsmtFullBath", "BsmtHalfBath", "GarageCars", "WoodDeckSF", "MiscVal", "PoolArea")

test <- test[!drops]

# find percent missing
pMiss <- function(x){sum(is.na(x))/length(x)*100}

missing_train <- data.frame(cbind(apply(train,2,pMiss)))
missing_train$Feature <- rownames(missing_train)
colnames(missing_train)[1] <- 'Value'

missing_test <- data.frame(cbind(apply(test,2,pMiss)))
missing_test$Feature <- rownames(missing_test)
colnames(missing_test)[1] <- 'Value'

missing <- cbind(missing_test[,1], missing_train[-20,1])
rownames(missing) <- rownames(missing_test)
colnames(missing) <- c('Test', 'Train')
missing <- data.frame(missing)

pander(missing[with(missing, order(-Test, -Train)), ][1:5,], caption='Percent missing')

tempvars <- c('MSZoning', 'GarageType', 'GarageArea', 'TotSqFt')
summary(test[, tempvars])

summary(train$GarageType)

# impute Garage area
# mean(test$GarageArea, na.rm=T) # 473
# which(is.na(test$GarageArea)) # get NA values
# test$GarageArea[1117]  # check
test$GarageArea[1117] <- 473 # update
# test$GarageArea[1117]  # check

# impute Garage area
# mean(test$TotSqFt, na.rm=T) # 1974.728
# which(is.na(test$TotSqFt)) # get NA values
# test$TotSqFt[661]  # check
test$TotSqFt[661] <- 1975 # update
# test$TotSqFt[661]  # check

# impute Zoning
# mode is 'RL'
# which(is.na(test$MSZoning)) # get NA values
# test$MSZoning[c(456, 757, 791, 1445)]
test$MSZoning[c(456, 757, 791, 1445)] <- 'RL'
# test$MSZoning[c(456, 757, 791, 1445)]

# impute Garage type by randomly assigning proportional values
train_rows_to_impute <- rownames(train[which(is.na(train$GarageType)), ])
test_rows_to_impute <- rownames(test[which(is.na(test$GarageType)), ])

# training data first
# get proportions
# summary(train$GarageType)
# trgarage <- c(6, 870, 19, 88, 9, 387)
# round(prop.table(trgarage)*81,1)
# sum(0, 51, 1, 5, 1, 23)

# randomly select fills
Attchd <- sample(train_rows_to_impute, 51)
temp <- train_rows_to_impute[!train_rows_to_impute %in% Attchd] 
Detchd <- sample(temp, 23)
temp <- temp[!temp %in% Detchd] 
BuiltIn <- sample(temp, 5)
temp <- temp[!temp %in% BuiltIn] 
Basment <- sample(temp, 1)
CarPort <- temp[!temp %in% Basment] 

# make changes
train[Attchd, "GarageType"] <- 'Attchd'
train[Detchd, "GarageType"] <- 'Detchd'
train[BuiltIn, "GarageType"] <- 'BuiltIn'
train[Basment, "GarageType"] <- 'Basment'
train[CarPort, "GarageType"] <- 'CarPort'

# check for NAs
# length(which(is.na(train$GarageType)==TRUE))

# test data, same technique
# tsgarage<- c(17, 853, 17, 98, 6, 392)
# round(prop.table(tsgarage)*76,1)
# sum(1, 47, 1, 5, 0, 22)

# randomly select fills
Attchd <- sample(test_rows_to_impute, 47)
temp <- test_rows_to_impute[!test_rows_to_impute %in% Attchd] 
Detchd <- sample(temp, 22)
temp <- temp[!temp %in% Detchd] 
BuiltIn <- sample(temp, 5)
temp <- temp[!temp %in% BuiltIn] 
Basment <- sample(temp, 1)
twoTypes <- temp[!temp %in% Basment] 

# make changes
test[Attchd, "GarageType"] <- '2Types'
test[Detchd, "GarageType"] <- 'Detchd'
test[BuiltIn, "GarageType"] <- 'BuiltIn'
test[Basment, "GarageType"] <- 'Basment'
test[twoTypes, "GarageType"] <- '2Types'

# check for NAs
length(which(is.na(train[,])==TRUE)); length(which(is.na(test[,])==TRUE))

# drop ID
test <- test[, -1]
train <- train[,-1]

pander(head(train,3))

pander(head(test,3))

numerics <- colnames(train[, sapply(train, is.numeric)])
describe(train[,numerics])[11]

# make copies for the modeling; make a small train and test set based on earlier testing
df_train <- train
df_test <- test
df_small <- df_train[, -c(5,9,14,17)]
df_test_small <- df_test[, -c(5,9,14,17)]

# make sure that OverallQual is a properly ordered factor
df_small$OverallQual <- factor(df_small$OverallQual, ordered=TRUE, levels=c(1:10))
df_test_small$OverallQual <- factor(df_test_small$OverallQual, ordered=TRUE, levels=c(1:10))

# using caret package, set up cross validation and our comparison metric
trainControl <- trainControl(method="repeatedcv", number=10, repeats=5)
metric <- "RMSE"

# build the four test models
# 
# LM
set.seed(100)
fit2.lm <- train(SalePrice~., data=df_small, method="lm", metric=metric, preProc=c("center", "scale", "BoxCox"), trControl=trainControl)

# Random Forest
set.seed(100)
fit2.rf <- train(SalePrice~., data=df_small, method="rf", metric=metric,
                preProc=c("BoxCox"),
                trControl=trainControl)

# Stochastic Gradient Boosting
set.seed(100)
fit2.gbm <- train(SalePrice~., data=df_small, method="gbm", metric=metric,
                 preProc=c("BoxCox"),
                 trControl=trainControl, verbose=FALSE)

# Cubist
set.seed(100)
fit2.cubist <- train(SalePrice~., data=df_small, method="cubist", metric=metric, preProc=c("BoxCox"), trControl=trainControl)

# Compare algorithms
ensembleResults <- resamples(list(LM=fit2.lm, GBM=fit2.gbm, RF=fit2.rf, CUBIST=fit2.cubist))

summary(ensembleResults)

# Tune the Cubist algorithm
trainControl <- trainControl(method="repeatedcv", number=10, repeats=3)
metric <- "RMSE"

set.seed(100)
grid <- expand.grid(.committees=seq(5, 15, by=1), .neighbors=c(3, 5, 7))

tune.cubist <- train(SalePrice~., data=df_small, method="cubist", metric=metric,
                  preProc=c("BoxCox"), tuneGrid=grid, trControl=trainControl)

plot(tune.cubist)

print(tune.cubist)

cubist.predictions <- predict(tune.cubist)

plot.new()
par(mfrow=c(1,2))
hist(cubist.predictions, freq=F,
      main='Cubist Sales Price',
      col='lightblue', breaks=30,
      xlim=c(50000, 600000))
lines(density(cubist.predictions), col='red', lwd=3)

hist(df_small$SalePrice, freq=F,
      main='Training Sale Price',
      col='lightblue', breaks=30,
      xlim=c(50000, 600000))
lines(density(df_small$SalePrice), col='red', lwd=3)

cubist_test.predictions <- predict(tune.cubist, df_test_small)

plot.new()
par(mfrow=c(1,2))
hist(cubist_test.predictions, freq=F,
      main='Cubist Test Predictions',
      col='lightblue', breaks=30,
      xlim=c(50000, 600000))
lines(density(cubist_test.predictions), col='red', lwd=3)

hist(df_small$SalePrice, freq=F,
      main='Training Sales Price',
      col='lightblue', breaks=30,
      xlim=c(50000, 600000))
lines(density(df_small$SalePrice), col='red', lwd=3)

# root mean squared logarithmic error
# msle(df_small$SalePrice, cubist.predictions)^0.5

# variable importance
varImp(tune.cubist)

# prepare the data transforms
set.seed(100)
x <- df_small[,c(1:13,15)]
y <- df_small[,14]
preprocessParams <- preProcess(x, method=c("BoxCox"))
transX <- predict(preprocessParams, x)
# train the final model
finalModel <- cubist(x=transX, y=y, committees=13)

print(finalModel)

# transform the test data and predict
set.seed(100)
valX <- df_test_small[, 1:14]
trans_valX <- predict(preprocessParams, valX)

# use final model to make predictions on the validation dataset
cubist_final.predictions <- predict(finalModel, newdata=trans_valX, neighbors=3)

# save the predictions for submission
write.csv(cubist_final.predictions, "kaggle_cubist.predictions.csv")