{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"library(tidyverse)\nlibrary(rpart)\nlibrary(randomForest)\nlibrary(caret)\nlibrary(rpart)\nlibrary(rpart.plot)\nlibrary(knitr)\nlibrary(ggplot2)\nlibrary(plyr)\nlibrary(dplyr)\nlibrary(corrplot)\nlibrary(gridExtra)\nlibrary(scales)\nlibrary(Rmisc)\nlibrary(ggrepel)\nlibrary(psych)\nlibrary(xgboost)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:55:25.138858Z","iopub.execute_input":"2022-07-25T00:55:25.191235Z","iopub.status.idle":"2022-07-25T00:55:25.279865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train <- read.csv('../input/house-prices-advanced-regression-techniques/train.csv')","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:55:26.650019Z","iopub.execute_input":"2022-07-25T00:55:26.655688Z","iopub.status.idle":"2022-07-25T00:55:26.750226Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"head(train)\ncolnames(train)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:55:27.906741Z","iopub.execute_input":"2022-07-25T00:55:27.909985Z","iopub.status.idle":"2022-07-25T00:55:27.986286Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"DATA CLEANING","metadata":{}},{"cell_type":"code","source":"#make paved street a binary variable\ntrain$paved <- ifelse(train$Street == 'Pave', 1, 0)\ntable(train$Street)\ntable(train$paved)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:55:30.188634Z","iopub.execute_input":"2022-07-25T00:55:30.190192Z","iopub.status.idle":"2022-07-25T00:55:30.218531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#make lot size a binary variable of regular (1) or irregular (0)\ntrain$LotReg <- ifelse(train$LotShape == 'Reg', 1, 0)\ntable(train$LotShape)\ntable(train$LotReg)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:55:32.516145Z","iopub.execute_input":"2022-07-25T00:55:32.517645Z","iopub.status.idle":"2022-07-25T00:55:32.544677Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"CORRELATION TABLES\n\nTo be added for further analysis/data exploration","metadata":{}},{"cell_type":"markdown","source":"DATA CLEANUP AND PREPARATION","metadata":{}},{"cell_type":"code","source":"which(colSums(is.na(train))> 0)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:55:35.339418Z","iopub.execute_input":"2022-07-25T00:55:35.341157Z","iopub.status.idle":"2022-07-25T00:55:35.361136Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"unique(train$Fence)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:55:37.045204Z","iopub.execute_input":"2022-07-25T00:55:37.047003Z","iopub.status.idle":"2022-07-25T00:55:37.063509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Lot frontage with NA values likely mean that the frontage is actually 0\ntrain$LotFrontage[is.na(train$LotFrontage)] <- 0\n\n#likewise MAsontry Veneer area being NA would be when there is no Masonry Veneer\ntrain$MasVnrArea[is.na(train$MasVnrArea)] <- 0\ntrain$MasVnrType[is.na(train$MasVnrType)] <- 'other'\n\n#the last data where NA values are an issue is Year Garage built. We can either set the NA values to 0, or we could simply remove the column.\n#For now, we are going to remove it. Further analysis will be necessary here\n\ndrop <- c('GarageTrBlt')\ntrain <- train[!(names(train) %in% drop)]\n\n\n#for the rest of the NA values, by the description, they are categorical variables where NA is equivalent to None\n\nwhich(colSums(is.na(train))> 0)\ntrain[is.na(train)] <- 'NA'","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:55:38.091993Z","iopub.execute_input":"2022-07-25T00:55:38.093648Z","iopub.status.idle":"2022-07-25T00:55:38.128328Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#split the data into test and train values\nResponse <- train$SalePrice\n\npartition <- createDataPartition(y = Response, p = 0.8, list = F)\n\ntraining <- train[partition,]\ntesting <- train[-partition,]\n\ndim(training)\ndim(testing)\n","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:55:39.402870Z","iopub.execute_input":"2022-07-25T00:55:39.404402Z","iopub.status.idle":"2022-07-25T00:55:39.448832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Basic Linear Regression","metadata":{}},{"cell_type":"code","source":"model.basicLM <- lm(formula = SalePrice ~ . , data = training)\nsummary(model.basicLM)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:55:41.003037Z","iopub.execute_input":"2022-07-25T00:55:41.004575Z","iopub.status.idle":"2022-07-25T00:55:41.171741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"RANDOM FOREST","metadata":{}},{"cell_type":"code","source":"model.RF <- randomForest(SalePrice ~. , data = training)\nsummary(model.RF)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:55:42.687969Z","iopub.execute_input":"2022-07-25T00:55:42.689517Z","iopub.status.idle":"2022-07-25T00:55:50.930705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"REGRESSION TREE","metadata":{"execution":{"iopub.status.busy":"2022-07-05T00:41:23.876894Z","iopub.execute_input":"2022-07-05T00:41:23.878629Z","iopub.status.idle":"2022-07-05T00:41:23.913738Z"}}},{"cell_type":"code","source":"model.tree1 <- rpart(SalePrice ~., data = training, method = 'anova')\nsummary(model.tree1)\nrpart.plot(model.tree1)\nplotcp(model.tree1)\n","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:55:58.635411Z","iopub.execute_input":"2022-07-25T00:55:58.636795Z","iopub.status.idle":"2022-07-25T00:55:59.707081Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cp <- 0.06193911\nmodel.tree1.prune <- prune(model.tree1, cp = 0.03375881)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:01.681010Z","iopub.execute_input":"2022-07-25T00:56:01.682744Z","iopub.status.idle":"2022-07-25T00:56:01.699769Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rpart.plot(model.tree1.prune)\nplotcp(model.tree1.prune)\nsummary(model.tree1.prune)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:03.790012Z","iopub.execute_input":"2022-07-25T00:56:03.791577Z","iopub.status.idle":"2022-07-25T00:56:04.147833Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"put training and testing together to get to know the average house price by data","metadata":{}},{"cell_type":"code","source":"together <- rbind(training, testing)\ndim(together)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:06.424850Z","iopub.execute_input":"2022-07-25T00:56:06.426394Z","iopub.status.idle":"2022-07-25T00:56:06.448873Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Using the geom_histogram in ggplot to draw the visual report of house price","metadata":{}},{"cell_type":"code","source":"library(scales)\nggplot(data=together[!is.na(together$SalePrice),], aes(x=SalePrice)) +\n                                                     geom_histogram(fill=\"yellow\", binwidth = 15000) +\n                                                     scale_x_continuous(breaks= seq(0, 800000, by=100000), labels = comma)\n\nsummary (together$SalePrice)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:08.310073Z","iopub.execute_input":"2022-07-25T00:56:08.311521Z","iopub.status.idle":"2022-07-25T00:56:08.648308Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The histogram about the house price are right skewed. That means most people usually spent their money to buy cheap house and just few of them can pay for expensive fancy houses.\nThe summary values tell us about the minimum, median, maximum values of house price and the median price is $163000.","metadata":{}},{"cell_type":"code","source":"num_variables <- which(sapply(together, is.numeric))\nnum_variables_names <- names(num_variables)\nlength(num_variables)\n\nnum_var <- together[, num_variables]\ncorr_variables <- cor(num_var, use=\"pairwise.complete.obs\")\n\nsorted_corr <- as.matrix(sort(corr_variables[,'SalePrice'], decreasing = TRUE))\nhigh_rate <- names(which(apply(sorted_corr, 1, function(x) abs(x) > 0.5)))\ncorr_variables <- corr_variables[high_rate, high_rate]\n                               \ncorrplot.mixed(corr_variables, tl.col=\"black\", tl.pos = \"lt\")\n\n","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:10.755217Z","iopub.execute_input":"2022-07-25T00:56:10.756777Z","iopub.status.idle":"2022-07-25T00:56:10.940628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Seleced columns which contains numeric values. There are 38 columns in this data. \nAnd, selected the columns which is more than 0.5 correlation ratio. \n\n\nWhich house features/criteria have the largest magnitude in price change? \n\nThe correlation between GarageCars and GarageArea is very high (0.89), and both have similar (high) correlations with SalePrice. \nThe other 6 six variables with a correlation higher than 0.5 with \nSalePrice are: -TotalBsmtSF: Total square feet of basement area \n-1stFlrSF: First Floor square feet \n-FullBath: Full bathrooms above grade \n-TotRmsAbvGrd: Total rooms above grade (does not include bathrooms) \n-YearBuilt: Original construction date \n-YearRemodAdd: Remodel date (same as construction date if no remodeling or additions)","metadata":{}},{"cell_type":"code","source":"ggplot(data=together[!is.na(together$SalePrice),], aes(x=factor(OverallQual), y=SalePrice))+\n        geom_boxplot(col='red') + labs(x='Overall Quality') +\n        scale_y_continuous(breaks= seq(0, 1000000, by=150000), labels = comma)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:13.415000Z","iopub.execute_input":"2022-07-25T00:56:13.416673Z","iopub.status.idle":"2022-07-25T00:56:13.803758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The positive correlation is certainly there indeed, and seems to be a slightly upward curve. \nRegarding outliers, I do not see any extreme values. \nIf there is a candidate to take out as an outlier later on, it seems to be the expensive house with grade 4.","metadata":{}},{"cell_type":"code","source":"ms1 <- ggplot(together[!is.na(together$SalePrice),], aes(x=YrSold, y=SalePrice)) +\n        geom_bar(stat='summary', fun.y = \"median\", fill='blue') +\n        theme(axis.text.x = element_text(angle = 45, hjust = 1)) +\n        scale_y_continuous(breaks= seq(0, 800000, by=50000), labels = comma) +\n        geom_label(stat = \"count\", aes(label = ..count.., y = ..count..), size=3) +\n        geom_hline(yintercept=163000, linetype=\"dashed\", color = \"red\")\n\n\ngrid.arrange(ms1)\n\neight <- filter(together, YrSold==2008)\nnrow(together)\nnrow(eight)\nq1 <- ggplot(data=eight, aes(x=as.factor(OverallQual))) +\n        geom_histogram(stat='count')\n\nlayout <- matrix(c(1,2,8,3,4,8,5,6,7),3,3,byrow=TRUE)\nmultiplot(q1, layout=layout)\n","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:16.523062Z","iopub.execute_input":"2022-07-25T00:56:16.524643Z","iopub.status.idle":"2022-07-25T00:56:17.136321Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Using the GrLivArea to find out the outlier as 524,1229 and remove it our final model","metadata":{}},{"cell_type":"code","source":"ggplot(data = together[!is.na(together$SalePrice),], aes(x=GrLivArea, y=SalePrice))+\n        geom_point(col='blue') + geom_smooth(method = \"lm\", se=FALSE, color=\"black\", aes(group=1)) +\n        scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma) +\n        geom_text_repel(aes(label = ifelse(together$GrLivArea[!is.na(together$SalePrice)]>4500, rownames(together), '')))","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:19.842416Z","iopub.execute_input":"2022-07-25T00:56:19.843971Z","iopub.status.idle":"2022-07-25T00:56:20.292345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"remove highest correlation between two columns and outliers. \nand using the qqplot to how it looks like it. ","metadata":{}},{"cell_type":"code","source":"drop <- c('YearRemodAdd', 'GarageYrBlt', 'GarageArea', 'GarageCond', 'TotalBsmtSF', 'TotalRmsAbvGrd', 'BsmtFinSF1')\n\ntogether <- together[,!(names(together) %in% drop)]\n\ntogether <- together[-c(524, 1299),]\n\nhead(together)\ncolnames(together)\nskew(together$SalePrice)\n\nqqnorm(together$SalePrice)\nqqline(together$SalePrice)\n\n\n","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:22.314422Z","iopub.execute_input":"2022-07-25T00:56:22.315979Z","iopub.status.idle":"2022-07-25T00:56:22.518719Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"1.87 value is too high, it is not normally distributed. using the log of sales price to fix it.\nI change to 1.12 and qq plot look so better than first one. ","metadata":{}},{"cell_type":"code","source":"together$SalePrice <- log(together$SalePrice)\nskew(together$SalePrice)\n\nqqnorm(together$SalePrice)\nqqline(together$SalePrice)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:25.014067Z","iopub.execute_input":"2022-07-25T00:56:25.015666Z","iopub.status.idle":"2022-07-25T00:56:25.162030Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"length(together)\nnum_variables <- which(sapply(together, is.numeric))\nnum_variables_names <- names(num_variables)\nlength(num_variables)\n\nmodel <- together[, num_variables]\nhead(model)\n\nidx <- sample(1:nrow(model), size = nrow(model)*0.8, replace=FALSE)\ntraining <- model[idx,]\ntesting <- model[-idx,]\n\nhead(testing)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:27.478952Z","iopub.execute_input":"2022-07-25T00:56:27.480510Z","iopub.status.idle":"2022-07-25T00:56:27.551252Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nmy_control <-trainControl(method=\"cv\", number=5)\nlassoGrid <- expand.grid(alpha = 1, lambda = seq(0.001,0.1,by = 0.0005))\nlasso_mod <- train(x=training, y=training$SalePrice, method='glmnet', trControl= my_control, tuneGrid=lassoGrid) \nlasso_mod$bestTune\nmin(lasso_mod$results$RMSE)\n\nlassoVarImp <- varImp(lasso_mod,scale=F)\nlassoImportance <- lassoVarImp$importance\n\nLassoPred <- predict(lasso_mod, training)\npredictions_lasso <- exp(LassoPred) #need to reverse the log to the real values\nhead(predictions_lasso)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:29.059221Z","iopub.execute_input":"2022-07-25T00:56:29.060770Z","iopub.status.idle":"2022-07-25T00:56:30.342839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xgb_grid = expand.grid(\nnrounds = 1000,\neta = c(0.1, 0.05, 0.01),\nmax_depth = c(2, 3, 4, 5, 6),\ngamma = 0,\ncolsample_bytree=1,\nmin_child_weight=c(1, 2, 3, 4 ,5),\nsubsample=1\n)\n\n#xgb.caret <- train(x=training, y=training$SalePrice, method = 'xgbTree', trControl = my_control, tuneGrid = xgb_grid)\n#xgb.caret$bestTune","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:31.967921Z","iopub.execute_input":"2022-07-25T00:56:31.969569Z","iopub.status.idle":"2022-07-25T00:56:31.985080Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"the best parameters are nrounds = 1000, max_depth = 3, eta = 0.01, gamma = 0, colsample_bytree = 1, min_child_weight = 1, subsample = 1","metadata":{}},{"cell_type":"code","source":"training_xgb <- xgb.DMatrix(data = as.matrix(training), label = training$SalePrice)\nset.seed(123)\nxgb.model <- xgb.cv(data = training_xgb, \n                    params = list(booster = \"gbtree\", eta = 0.01, gamma = 0, max_depth = 3,\n                                  min_child_weight = 1, subsample = 1, colsample_bytree = 1,\n                                  objective = \"reg:squarederror\", eval_metric = \"rmse\"),\n                    nrounds = 500, early_stopping_rounds = 10, nfold = 5,\n                    showsd = T, stratified = T, print_every_n = 40, maximize = F)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:34.852053Z","iopub.execute_input":"2022-07-25T00:56:34.853796Z","iopub.status.idle":"2022-07-25T00:56:38.570566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"best nrow is 500","metadata":{}},{"cell_type":"code","source":"xgb.model <- xgb.train(data = training_xgb,\n                       params = list(booster = \"gbtree\", eta = 0.01, gamma = 0, max_depth = 3,\n                                     min_child_weight = 1, subsample = 1, colsample_bytree = 1,\n                                     objective = \"reg:squarederror\", eval_metric = \"rmse\"),\n                       nrounds = 500)\n\n\ntesting_xgb <- xgb.DMatrix(data = as.matrix(testing), label = testing$SalePrice)\n\ntest.predict.xgboost.log <- predict(xgb.model, testing_xgb)\ntest.predict.xgboost <- exp(test.predict.xgboost.log)   \nhead(test.predict.xgboost)\n","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:40.873841Z","iopub.execute_input":"2022-07-25T00:56:40.875352Z","iopub.status.idle":"2022-07-25T00:56:41.452551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub.avg <- data.frame(SalePrice = predictions_lasso)\nsummary(sub.avg)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:56:43.708993Z","iopub.execute_input":"2022-07-25T00:56:43.710756Z","iopub.status.idle":"2022-07-25T00:56:43.732076Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"write.csv(sub.avg, file = 'saleprice.csv', row.names = F)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T00:58:12.402971Z","iopub.execute_input":"2022-07-25T00:58:12.408274Z","iopub.status.idle":"2022-07-25T00:58:12.430451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"boost.avg <- data.frame(SalePrice = test.predict.xgboost)\nsummary(boost.avg)\n\nwrite.csv(boost.avg, file = 'saleprice_boost.csv', row.names = F)","metadata":{"execution":{"iopub.status.busy":"2022-07-25T01:04:03.225196Z","iopub.execute_input":"2022-07-25T01:04:03.226729Z","iopub.status.idle":"2022-07-25T01:04:03.261132Z"},"trusted":true},"execution_count":null,"outputs":[]}]}