{"cells":[{"metadata":{},"cell_type":"markdown","source":"# 1 Introduction","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## 1.1 Business Goal\nAssessing the value of energy efficiency improvements can be challenging as there's no way to truly know how much energy a building would have used without the improvements. The best we can do is to build counterfactual models. Once a building is overhauled the new (lower) energy consumption is compared against modeled values for the original building to calculate the savings from the retrofit. More accurate models could support better market incentives and enable lower cost financing.\n\n## 1.2 Technical Goal\nIn the ASHRAE - GREAT ENERGY PREDICTION III dataset we are required to use the 4 types of energy cost records (meter readings) of 1448 buildings across 16 sites in three years, with their building information and environment parameters, to build a regression model to predict the energy cost (meter reading) of a certain building with specific time and environment conditions. This is an important reference to check if the retrofits of a building can really lower the energy consumption.\n\n## 1.3 Datasets and Variables\nDatasets contain 5 files, 3 types:\n* **train/test.csv**: meter reading records of the buildings of each meter type.\n* **building_meta.csv**: information of the buildings.\n* **weather_[train/test].csv**: the weather parameters of each sites where the buildings are located.\n\nAs the test dataset does not contain the true meter reading records (which is calculated by kaggle after submission), so in this mission we only use the data from the training files. We will split the data into training, validation and testing datasets randomly.\n\nA general overview of the variables each file contains:\n\n![](https://drive.google.com/uc?id=1o9z4uT0YhSyfsd2OaPgN08On-dtAEhbB)","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# 2 Preparation","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Libraries that are needed:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"library(lubridate)\nlibrary(tidyverse)\nlibrary(scales)\nlibrary(ggridges)\nlibrary(corrgram)\nlibrary(data.table)\nlibrary(caret)\nlibrary(xgboost)\nset.seed(15)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Common theme setting for plotting:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# set plotting theme baseline\ntheme_set(theme_minimal() +\n            theme(axis.title.y=element_text(size=15), \n                  axis.title.x=element_text(size=15),\n                  plot.title=element_text(hjust=0.5, size=18, color = \"grey25\", face = \"bold\"), \n                  plot.subtitle=element_text(hjust=0.5, size=16, color=\"grey44\"), \n                  axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 12),\n                  legend.title=element_text(hjust=0.9, size=12), legend.text = element_text(size=12),\n                  panel.grid.major = element_line(linetype = 2),\n                  panel.grid.minor = element_line(linetype = 2)))\n\ncol_pal <- c(\"#5EB296\", \"#4E9EBA\", \"#F29239\", \"#C2CE46\", \"#FF7A7F\", \"#4D4D4D\")\n# mint green; bean blue; orange; wasabi yellow; water red; dark gray","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 3 Data Analysis, Visualization and Preprocessing","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## 3.1 Load and A Glimpse of the Data\nFirst let's load and have a glimpse of the data.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# read in data\ntrain <- read.csv(\"../input/ashrae-energy-prediction/train.csv\")\nbuilding <- read.csv(\"../input/ashrae-energy-prediction/building_metadata.csv\")\nweather_train <- read.csv(\"../input/ashrae-energy-prediction/weather_train.csv\")\n\ndf <- train %>%\n  left_join(building, by=\"building_id\") %>%\n  left_join(weather_train, by=c(\"site_id\", \"timestamp\"))\n\nrm(train, building, weather_train)\nhead(df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.2 Datetime Format\nTimestamp is formatted and we extract date, month, day_number, weekday, hour, is_holiday variables from the timestamp.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"holidays <- as_datetime(\n  c(\"2016-01-01\", \"2016-01-18\", \"2016-02-15\", \n    \"2016-05-30\", \"2016-07-04\", \"2016-09-05\", \n    \"2016-10-10\", \"2016-11-11\", \"2016-11-24\",\n    \"2016-12-26\", \"2017-01-02\", \"2017-01-16\",\n    \"2017-02-20\", \"2017-05-29\", \"2017-07-04\", \n    \"2017-09-04\", \"2017-10-09\", \"2017-11-10\", \n    \"2017-11-23\", \"2017-12-25\"))\n\ndf <- mutate(df, \n             date = ymd(gsub( \" .*$\", \"\", timestamp)),\n             hour = hour(timestamp),\n             month = month(date),\n             weekday = wday(date),\n             day_number = day(date),\n             is_holiday = (as_datetime(timestamp) %in% holidays) + 0\n             )\n\nhead(df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.3 NA Values\nLet's take a look at the NA values of each variable and see if we can deal with specific variables.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"colSums(is.na(df))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.feat <- df %>% \n  select(building_id, meter, timestamp, meter_reading, site_id, primary_use, square_feet, year_built, floor_count, air_temperature, \n         cloud_coverage, dew_temperature, precip_depth_1_hr, sea_level_pressure, wind_direction, wind_speed) %>% \n  as.data.frame()\n\ndf.nonna <- colSums(!is.na(df.feat)) %>% \n  as.data.frame()\nnames(df.nonna)[1] <- \"nums\"\n\ndf.na <- colSums(is.na(df.feat)) %>%\n  as.data.frame()\nnames(df.na)[1] <- \"nums\"\n\ndf.nonna <- mutate(df.nonna, var=row.names(df.nonna), status=\"non_na\")\ndf.na <- mutate(df.na, var=row.names(df.na), status=\"na\")\ntotal.num <- rbind(df.nonna, df.na)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=8, repr.plot.height=8)\n\nggplot(total.num, aes(fill=status, y=nums, x=var)) +\n  geom_bar(position=\"fill\", stat=\"identity\") +\n  scale_fill_manual(values = col_pal) +\n  labs(y=\"percentage\") +\n  ggtitle(\"Percentage of NA Values and Non-NA Values in Train\") +\n  coord_flip() +\n  theme(axis.title.x=element_blank())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"rm(df.feat, df.na, df.nonna, total.num)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.4 Target Variable Analysis\nThe target variable is meter_reading, so we want to investigate and visualize the relationship of the target variable with other variables.\n### 3.4.1 meter_reading by date\nWe can see from the annual trend of meter reading that the energy usage is significantly larger overall in summer than in other seasons, which may indicate temperature has a significant impact on energy consumption","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=8, repr.plot.height=6)\ndf %>% \n  group_by(date) %>% \n  summarise(median_reading=median(meter_reading, na.rm=T)) %>% \n  ggplot(aes(x=date, y=median_reading)) +\n  geom_line(color=col_pal[1]) +\n  geom_smooth(se=F, color=col_pal[3]) +\n  ggtitle(\"Readings Really Peaked from May to October\", subtitle=\"Readings are the median reading per each reading\") +\n  labs(x=\"Reading Date\", y=\"Median Reading\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For each building type, most of them are with the overall meter_reading trend, except for manufacturing. And utility and warehouse have flatter curves and no big differences across seasons.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=16, repr.plot.height=12)\ndf %>% \n  group_by(date, primary_use) %>% \n  summarise(median_reading=median(meter_reading, na.rm=T)) %>% \n  ggplot(aes(x=date, y=median_reading)) +\n  geom_line(color=col_pal[1]) +\n  geom_smooth(se=F, color=col_pal[3]) +\n  ggtitle(\"Manufacturing Distinguishes from the General Trend\", subtitle=\"While some other business types remained fairly stable over the period\") +\n  labs(x=\"Reading Date\", y=\"Median Reading\") +\n  facet_wrap(~ primary_use, scales=\"free\") +\n  theme(axis.text.x=element_text(size=6, angle = 45),\n        axis.text.y=element_text(size=6))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3.4.2 Meter_reading by hours of a day\nThe overall trend of meter reading over the day shows the energy usage rises in working hours and decreases at night.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=8, repr.plot.height=6)\ndf %>% \n  group_by(hour) %>% \n  summarise(median_reading=median(meter_reading, na.rm=T)) %>% \n  ggplot(aes(x=hour, y=median_reading)) +\n  geom_line(color=col_pal[1], size=1) +\n  ggtitle(\"Readings wrt Hours\", subtitle=\"The highest point is in the middle of the day,\\nthe lowest point is around 4am\") +\n  labs(x=\"Reading Hour\", y=\"Median Reading\") +\n  theme(panel.grid.major.y=element_blank())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Although different building types present different curves, most of them are with the overall trend, except for Parking and Utility. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=16, repr.plot.height=12)\ndf %>% \n  group_by(hour, primary_use) %>% \n  summarise(median_reading=median(meter_reading, na.rm=T)) %>% \n  ggplot(aes(x=hour, y=median_reading)) +\n  geom_line(color=col_pal[1], size=1) +\n  ggtitle(\"Time Readings Are the Highest Changes Based on \\nBuilding Type\", subtitle = \"Residential, and Warehouse peak later in day,\\nUtility earlier\") +\n  labs(x=\"Reading Hour\", y=\"Median Reading\") +\n  facet_wrap(~ primary_use, scales=\"free\") +\n  theme(panel.grid.major.y=element_blank())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3.4.3 Meter_reading by primary_use\nThe daily meter reading values of utility, healthcare and education are large. While religious is generally small, it is reasonable to infer the lower frequency of religious activity than others. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=8)\ndf %>% \n  ggplot(aes(x=as.character(primary_use), y=log(meter_reading + 1))) +\n  geom_boxplot(alpha=0.5, fill=col_pal[1], color=col_pal[6]) +\n  ggtitle(\"Utilities and Healthcare Have the Highest Readings\", subtitle = \"Religious worship places the lowest\") +\n  labs(y= \"log(Meter Reading)\") +\n  coord_flip() +\n  theme(axis.title.y = element_blank())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3.4.4 Meter_reading of different type of buildings by weekdays\nNext is the weekly trend of the meter reading of each building type.  Judging from Q1 in the box chart, there was not much difference, with only religious which the reading on the weekend is slightly higher than the weekday. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=8)\ndf %>% \n  ggplot(aes(x=as.character(weekday), y=log(meter_reading + 1))) +\n  geom_boxplot(alpha=0.5, fill=col_pal[1], color=col_pal[6]) +\n  ggtitle(\"Places of Industry Highest Readings on Weekdays\", subtitle = \"Religious worship places have higher readings on weekends\") +\n  labs(y= \"log(Meter Reading)\") +\n  coord_flip() +\n  theme(axis.title.y = element_blank()) +\n  facet_wrap(~ primary_use, scales = \"free\") +\n    theme(axis.text.x=element_text(size=6),\n        axis.text.y=element_text(size=6))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.4.5 Outliers","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Then we take a look at the mean value of meter_reading. As mean is easier to be affected by extreme values than the median, we want to distinguish some anomaly using the mean value.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=6)\ndf %>% group_by(date) %>% summarise(mean_reading=mean(meter_reading, na.rm=T)) %>% \n  ggplot(aes(x=date, y=mean_reading)) +\n  geom_line(color=col_pal[1], size=1) +\n  ggtitle(\"Mean Readings by Day\") +\n  labs(x=\"Reading Date\", y=\"Mean Reading\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The mean meter reading explode from March to June and then drop hugely, with another peak at November. This huge difference does not seem normal. We want to distinguish the anomalous points, so first we plot the mean meter_reading by site.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=15, repr.plot.height=12)\ndf %>% \n  group_by(date, site_id) %>% \n  summarise(mean_reading=mean(meter_reading, na.rm=T)) %>% \n  ggplot(aes(x=date, y=mean_reading)) +\n  geom_line(color=col_pal[1]) +\n  ggtitle(\"Mean Readings by Site\") +\n  labs(x=\"Reading Date\", y=\"Mean Reading\") +\n  facet_wrap(~ site_id, scales=\"free\") +\n  theme(axis.text.x=element_text(size=6, angle = 45),\n        axis.text.y=element_text(size=6))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It is clear that the anomalous trend comes from site_id=13. Let's take a look at this site.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=15, repr.plot.height=12)\ndf %>% \n  filter(site_id == 13) %>% \n  group_by(date, primary_use) %>% \n  summarise(mean_reading=mean(meter_reading, na.rm=T)) %>% \n  ggplot(aes(x=date, y=mean_reading)) +\n  geom_line(color=col_pal[1]) +\n  ggtitle(\"Mean Readings by Primary Use of site_id=13\") +\n  labs(x=\"Reading Date\", y=\"Mean Reading\") +\n  facet_wrap(~ primary_use, scales=\"free\") +\n  theme(axis.text.x=element_text(size=6, angle = 45),\n        axis.text.y=element_text(size=6))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The Educational use buildings mainly contribute to the overall trend. This enable us narrow the scope further.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=10)\ndf %>% \n  filter(site_id == 13, primary_use == \"Education\") %>% \n  group_by(date, meter) %>% \n  summarise(mean_reading=mean(meter_reading, na.rm=T)) %>% \n  ggplot(aes(x=date, y=mean_reading)) +\n  geom_line(color=col_pal[1], size=1) +\n  scale_y_continuous(labels=comma) +\n  ggtitle(\"Mean Readings by Meter Type for site_id=13, primary_use=Education\") +\n  labs(x=\"Reading Date\", y=\"Mean Reading\") +\n  facet_grid(meter ~., scales=\"free\") +\n  theme(axis.text.x=element_text(size=6, angle = 45),\n        axis.text.y=element_text(size=6))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can distinguish the anomaly comes from meter=2, which is the steam meter.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=15, repr.plot.height=12)\ndf %>% \n  filter(site_id == 13, primary_use == \"Education\", meter == 2) %>% \n  group_by(date, building_id) %>% \n  summarise(mean_reading=mean(meter_reading, na.rm=T)) %>% \n  ggplot(aes(x=date, y=mean_reading)) +\n  geom_line(color=col_pal[1], size=0.5) +\n  scale_y_continuous(labels=comma) +\n  ggtitle(\"Mean Readings by Building ID for site_id=13, primary_use=Education, meter=2\") +\n  labs(x=\"Reading Date\", y=\"Mean Reading\") +\n  facet_wrap(~ building_id, scales=\"free\") +\n  theme(axis.text.x=element_text(size=6, angle = 45),\n        axis.text.y=element_text(size=6))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"By looking into the steam meter reading of each educational building in site_id 13, we found this anomaly is building_id=1099. So we remove the related meter_readings with this building_id. Now let's take a look at the remained data.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=8)\n\ndf.flt <- df %>% \n  filter(building_id!=1099) %>%\n  group_by(date) %>% summarise(mean_reading=mean(meter_reading, na.rm=T))\nggplot() + \n  geom_line(data=df.flt, aes(x=date, y=mean_reading), color=col_pal[2], size=1) +\n  scale_y_continuous(labels=comma) +\n  ggtitle(\"Mean Readings of All Data After Removing building_id=1099\") +\n  labs(x=\"Reading Date\", y=\"Mean Reading\") +\n  theme(axis.text.x=element_text(size=6, angle = 45),\n        axis.text.y=element_text(size=6))\n\nrm(df.flt)\ngc()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"outlier_bd = c(1099)\ndf <- df[-which(df$building_id %in% outlier_bd), ]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Then we have a look at the meter_reading value distribution itself.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"summary(df$meter_reading)\n\noptions(repr.plot.width=8, repr.plot.height=6)\ndf %>% \n  ggplot(aes(x=meter_reading)) +\n  geom_histogram(alpha=0.5, fill=col_pal[1], colour=col_pal[6]) +\n  scale_x_continuous(labels=comma) +\n  scale_y_continuous(labels=comma) +\n  ggtitle(\"All Meter Reading Distribution\", subtitle=\"Meter reading variables is heavily skewed by some outliers\") +\n  labs(x= \"Meter Reading\", y=\"Count\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see the meter_reading values is heavily skewed. Then we check for each meter type.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=9)\ndf %>% \n  ggplot(aes(x=meter_reading)) +\n  geom_histogram(alpha=0.5, fill=col_pal[1], colour=col_pal[6]) +\n  scale_x_continuous(labels=comma) +\n  scale_y_continuous(labels=comma) +\n  ggtitle(\"Meter Reading Distributions for Different Meters\", subtitle=\"Meter reading variables is heavily skewed by some outliers\") +\n  labs(x= \"Meter Reading\", y=\"Count\") +\n  facet_wrap(~ meter, scales=\"free\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The problem exists for each type of meter. So we split the outliers by 1.5*IQR.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"outlier <- round(1.5 * IQR(df$meter_reading),0)\n\noptions(repr.plot.width=14, repr.plot.height=6)\ndf %>% \n  mutate(outlier=ifelse(meter_reading > outlier, \"Outlier\", \"Not Outlier\")) %>% \n  ggplot(aes(x=meter_reading)) +\n  geom_histogram(alpha=0.5, fill=col_pal[1], color=col_pal[6]) +\n  scale_x_continuous(labels=comma) +\n  scale_y_continuous(labels=comma) +\n  ggtitle(\"Meter Reading wrt Outliers and Non-outliers\", subtitle=\"Meter reading variable is heavily skewed, even for non-outliers\") +\n  labs(x=\"Meter Reading\", y=\"Count\") +\n  facet_wrap(~ outlier, scales=\"free\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After removing the outliers the meter_reading distribution is more processible than before, but it is still skewed. This makes the model more sensible to the derivative (or to say, the change) of each input variables. So we will perform log transformation for modeling purpose, to make the variable more workable.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=8, repr.plot.height=6)\n\ndf %>% \n  ggplot(aes(x=log(meter_reading+1))) +\n  geom_histogram(alpha=0.5, fill=col_pal[1], color=col_pal[6], binwidth=0.5) +\n  scale_y_continuous(labels=comma) +\n  ggtitle(\"Meter Reading after Log Transforming\", subtitle=\"Log transforming for meter reading variable is neccessary\") +\n  labs(x=\"log(meter_reading+1)\", y=\"Count\")\n\ndf %>% \n  ggplot(aes(x=log(meter_reading+1), fill = as.character(meter))) +\n  geom_density(alpha=0.5, adjust=2) +\n  scale_fill_manual(values=col_pal) +\n  ggtitle(\"Meter Reading for Each Meter Types after Log Transforming\") +\n  labs(x=\"log(meter_reading+1)\") +\n  theme(axis.title.y=element_blank(), legend.position = \"top\")\n\nsummary(df$meter_reading)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Filter the dataframe.\noutlier_mr <- round(1.5 * IQR(df$meter_reading),0)\ndf <- mutate(df, outlier=ifelse(meter_reading > outlier_mr, 1, 0))\ndf <- filter(df, outlier==0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df$meter_reading = log(df$meter_reading + 1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.5 Other Variables\nThen we look at other variables that may be valuable to visualize and understand the data.\n### 3.5.1 square_feet\nFor square_feet, most of them are small buildings lie in the range of 32k to 139k, with a median value of about 107K.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=8, repr.plot.height=6)\ndf %>% \n  ggplot(aes(x=square_feet)) +\n  geom_density(color=col_pal[6], fill=col_pal[4], adjust=2, alpha=0.7) +\n  scale_x_continuous(labels=comma, name=\"Building Square Feet\") +\n  ggtitle(\"Building Size Is Positively Skewed\", subtitle = paste0(\"The median size is \", comma(median(df$square_feet)), \"square feet\")) +\n  theme(axis.title.y=element_blank(), axis.text.y=element_blank())\n\noptions(repr.plot.width=8, repr.plot.height=3)\ndf %>% \n  ggplot(aes(x=\"\", y=square_feet)) +\n  geom_boxplot(color=col_pal[6], fill=col_pal[4], alpha=0.7) +\n  labs(x=\"Building Square Feet\") +\n  ggtitle(\"Building Size (square_meter)\", subtitle=paste0(\"Building size is positively skewed\")) +\n  coord_flip() +\n  theme(axis.title.y=element_blank(), axis.text.y=element_blank())\n\nquantile(df$square_feet, c(0.25, 0.75))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Out of the same reason with meter_reading, we perform log transformation on square_feet variable.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df$square_feet = log(df$square_feet + 1) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=8, repr.plot.height=6)\ndf %>% \n  ggplot(aes(x=square_feet)) +\n  geom_density(color=col_pal[6], fill=col_pal[4], adjust=2, alpha=0.7) +\n  scale_x_continuous(labels=comma, name=\"Building Square Feet\") +\n  ggtitle(\"Building Size after Log Transformation\", subtitle = paste0(\"The median size is \", comma(median(log(df$square_feet))), \"log square feet\")) +\n  theme(axis.title.y=element_blank(), axis.text.y=element_blank())\n\noptions(repr.plot.width=8, repr.plot.height=3)\ndf %>% \n  ggplot(aes(x=\"\", y=square_feet)) +\n  geom_boxplot(color=col_pal[6], fill=col_pal[4], alpha=0.7) +\n  labs(x=\"Log Building Square Feet\") +\n  ggtitle(\"Building Size (log square_meter)\", subtitle=paste0(\"Building size density after log transformation\")) +\n  coord_flip() +\n  theme(axis.title.y=element_blank(), axis.text.y=element_blank())\n\nquantile(df$square_feet, c(0.25, 0.75))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3.5.2 Temperature\nThere are 96658 missing air_temperature and 100140 missing dew_temperature records in the training dataset. As the temperature are continuous and has its seasonality, we can interpolate the missing values using nearby days' records. We follow two rules:\n* 1. Take the air_temperature value for the previous timestamp as the current timestamp's value\n* 2. If the previous timestamp is also NA, then take the average of that day's temperature of that particular site_id","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fill_missing_temp <- function(dataframe, temp_var) {\n  \n  temp_var <- enquo(temp_var)\n  \n  dataframe %>% \n    arrange(site_id, date) %>% \n    group_by(site_id) %>% \n    mutate(last_temp = lag(!! temp_var)) %>% \n    ungroup() %>% \n    mutate(month_day = paste0(month, \"_\", weekday)) %>% \n    group_by(site_id, month_day) %>% \n    mutate(avg_daily_temp = mean(!! temp_var, na.rm = T)) %>%  \n    ungroup() %>% \n    mutate(temp_clean = ifelse(is.na(!! temp_var), last_temp, !! temp_var),\n           temp_clean = ifelse(is.na(temp_clean), avg_daily_temp, temp_clean)) %>% \n    select(-last_temp, -avg_daily_temp) %>% \n    pull(temp_clean)\n} ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Apply the function to ari_temperature and dew_temperature:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df$air_temp_clean <- fill_missing_temp(df, air_temperature)\ndf$dew_temp_clean <- fill_missing_temp(df, dew_temperature)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Air temperature with imputation:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=8, repr.plot.height=6)\ndf %>% \n  ggplot(aes(x=air_temp_clean)) +\n  geom_density(adjust=2, color=col_pal[6], fill=col_pal[4], alpha=0.7) +\n  labs(x=\"Air Temperature\", y= \"Density\") +\n  ggtitle(\"Air Temperature Distribution\", subtitle=paste0(\"Median air temperature is \", median(df$air_temp_clean, na.rm=T), \" degrees\")) +\n  theme(axis.title.y=element_blank())\n\noptions(repr.plot.width=8, repr.plot.height=3)\ndf %>% \n  ggplot(aes(x=\"\", y=air_temp_clean)) +\n  geom_boxplot(color=col_pal[6], fill=col_pal[4], alpha=0.7) +\n  labs(x=\"Air Temperature\", y= \"Density\") +\n  ggtitle(\"Most Air Temperatures Range between 8-24 Degrees\", subtitle=paste0(\"Median air temperature is \", median(df$air_temp_clean, na.rm=T), \" degrees\")) +\n  coord_flip()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Dew temperature with imputation:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=8, repr.plot.height=6)\ndf %>% \n  ggplot(aes(x=dew_temp_clean)) +\n  geom_density(adjust=2, color=col_pal[6], fill=col_pal[4], alpha=0.7) +\n  coord_cartesian(xlim=c(-30, 40)) +\n  labs(x=\"Dew Temperature\", y=\"Density\") +\n  ggtitle(\"Dew Temperatures Distribution\", subtitle=paste0(\"Median dew temperature is \", median(df$dew_temp_clean, na.rm=T), \" degrees\")) +\n  theme(axis.title.y=element_blank())\n\noptions(repr.plot.width=8, repr.plot.height=3)\ndf %>% \n  ggplot(aes(x=\"\", y=dew_temp_clean)) +\n  geom_boxplot(color=col_pal[6], fill=col_pal[4], alpha=0.7) +\n  labs(x=\"Air Temperature\", y= \"Density\") +\n  ggtitle(\"Most Dew Temperatures Range between 0-17 Degrees\", subtitle=paste0(\"Median air temperature is \", median(df$dew_temp_clean, na.rm=T), \" degrees\")) +\n  coord_flip()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3.5.3 meter type\nAmong the meters, electricity has the most records and hotwater the least","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=8, repr.plot.height=6)\ndf %>% \n  count(meter) %>% \n  ggplot(aes(x=meter, y=n)) +\n  geom_col(color=col_pal[6], fill=col_pal[4]) +\n  scale_y_continuous(labels=comma, name=\"Count of Readings\") +\n  ggtitle(\"Electricity the Most Frequent Meter Type Measured\", subtitle=\"Big drop off to Chilledwater, while Hotwater distant last\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"### 3.5.4 year_built\nFor year built, the number of buildings distribute evenly, except that in 1976, when 55 new buildings were recorded. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=8, repr.plot.height=6)\ndf %>% \n  distinct(building_id, .keep_all=T) %>% \n  group_by(year_built) %>% \n  summarise(n_buildings=n()) %>%\n  ggplot(aes(x=year_built, y=n_buildings)) +\n  geom_line(color=col_pal[4], size=1) +\n  ggtitle(\"Building Numbers wrt Year of Built\") +\n  scale_y_continuous(name=\"Number Built\", limits=c(0,100)) +\n  labs(x=\"Year Built\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3.5.5 floor counts","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=8, repr.plot.height=6)\ndf %>% \n  distinct(building_id, .keep_all = T) %>% \n  ggplot(aes(x=floor_count)) +\n  geom_density(adjust=2, color=col_pal[6], fill=col_pal[4], alpha=0.7) +\n  labs(x=\"Floor Count\", y= \"Density\") +\n  ggtitle(\"Not Many Buildings with More Than 10 Floors\", subtitle = paste0(\"Median floors is \", median(df$floor_count, na.rm = T), \" floors\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.6 Correlation between meter_reading and other variables","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"We plot the correlation between the target variable meter_reading with other input variables. Among them we can see there is a high positive correlation between square_feet and meter_reading, and a moderate positive correlation between floor_count and meter_reading. But floor_count has a lot of missing values (82.65%). Other variables have weak or no correlation with meter_reading.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\n\noptions(repr.plot.width=12, repr.plot.height=8)\ndf %>% \n  select(meter_reading, square_feet, year_built, floor_count, air_temp_clean, cloud_coverage, \n         dew_temp_clean, precip_depth_1_hr, sea_level_pressure, wind_direction, wind_speed,\n         month, day_number, weekday, hour, is_holiday) %>% na.omit() %>%\n  as.matrix() %>% \n  cor() %>% as.data.frame() %>% \n  select(meter_reading) %>% mutate(vars=row.names(.)) %>% \n  filter(vars !=\"meter_reading\") %>% \n  ggplot(aes(x=reorder(vars, meter_reading), y= meter_reading)) +\n  geom_col(fill=col_pal[1], colour = col_pal[6]) +\n  labs(x=\"Variables\", y=\"Correlation with meter_reading\") +\n  ggtitle(\"Variables Correlations with Meter Reading\", subtitle=\"Square_feet and floor_count show high positive correlation,\\nand floor_count shows moderate positive correlation\") +\n  coord_flip() +\n  theme(panel.grid.major.y=element_blank())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The correlations between variables are also plotted. Categorical variables (building_id, site_id, meter, primary_use, is_holiday) are not analyzed.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\n\noptions(repr.plot.width=16, repr.plot.height=16)\ndf %>% \n  select(meter_reading, square_feet, year_built, floor_count, air_temperature, cloud_coverage, \n         dew_temperature, precip_depth_1_hr, sea_level_pressure, wind_direction, wind_speed,\n         month, day_number, weekday, hour) %>% na.omit() %>%\n  corrgram(order=TRUE, lower.panel=panel.ellipse,\n           upper.panel=panel.pts, text.panel=panel.txt,\n           diag.panel=panel.minmax,\n           main=\"Correlation between input variables\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From above we can have some conclusions:\n1. we can hardly obeserve a correlation between the variables;\n2. among them, only air_temperature and dew_temperature tends to have an obvious positive correlation (about 0.7);\n3. the target variable meter_reading tends to have positive correlation with square_feet and floor_count(weak).\n\nWe will later analyze if the correlation affects the model.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## 3.7 Data Downsampling","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"As processing the whole dataset in the following, especially when training the model, will take up large RAM, which is not able to run locally and even with the kaggle kernel of 16G RAM. So we downsampled the data with a time interval of 4h. The original data was recorded every one hour.\n\nAnd to reduce the RAM consumption of former processing, we write the cleaned and downsampled dataframe to a new csv file. And when building and training the model, we will not run the former data cleaning procedures, but only read from the disk with the cleaned 4h-interval dataset.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\ndf <- df[which(df$hour%%4==0),]\nwrite.csv(df, \"ds4.csv\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.8 Encoding for Categorical Features\nBecause the model we use (decision tree and xgboost) does not support categorical features, so we need to encode those feature with specific methods. The concerned features are:\n* **building_id**: the unique id of each building, int, 0~1448\n* **site_id**: the unique id of each site, int, 0~15\n* **meter**: the id of meter type, int, 0~3\n* **primary_use**: the type of the building, factor, 16 categories\n\nAs the variables are all **nominal variables** but not ordinal variables, which means the value of its id does not reflect an internal order or quantity of that feature, we need to use **one-hot encoding** for these features. If we use ordinal encoding, the algorithm will try to find a splitting point of the consecutive values of the feature and split the data into two parts, and calculates the rmse with respect to each part. This will limit the combination of the elements in each part, that is, only consecutive elements can be splitted into one part. This is not what we want.\n\nFor the last three features, it is ok to encode them as one-hot vectors; but for building_id, if we also use one-hot encoding the cardinality of input features will rise stupendously. So we have to find a compromise, and we decided to use the **binary encoding** for building_id. \n\n### 3.8.1 Binary encoding for building_id\nFor convenience and speed the binary encoding is performed using python. And the encoded file is uploaded for follow-up process. Here we just paste the python code here.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"    file = r'E:\\Ms-UNICE\\DataValorization\\kaggle1\\ds4\\ds4.csv'\n    df = pd.read_csv(file)\n    \n    def addIdDigit(df, digit, idSource):\n        return df[idSource][digit]\n\n    idMap = {oriId: '{:011b}'.format(oriId) for oriId in set(df['building_id'])}\n\n    df['bdid_coded'] = df['building_id'].map(idMap)\n    \n    for i in range(11):\n        df['id_%d' % i] = df.apply(addIdDigit, axis=1, digit=i, idSource='bdid_coded')\n\n    df.to_csv(r'E:\\Ms-UNICE\\DataValorization\\kaggle1\\ds4\\train_clean_ds4_bin.csv')","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### 3.8.2 One-hot encoding for primary_use, meter and site_id\nThe data with binary-encoded building_id is read to perform the follow-up encoding.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\n\ntrain <- read.csv(\"../input/ds4-bin/train_clean_ds4_bin.csv\")\nfeat.onehot = c(\"primary_use\", \"meter\", \"site_id\")\n\ntrain$meter = as.factor(train$meter)\ntrain$site_id = as.factor(train$site_id)\ndummies <- dummyVars(~primary_use+meter+site_id, data=train)\ndf.onehot <- as.data.frame(predict(dummies, newdata=train))\ndf <- cbind(train[,-c(which(colnames(train) %in% feat.onehot))], df.onehot)\n\n# remove the useless columns\ndf <- df[, -which(colnames(df) %in% c('X', 'timestamp', 'building_id', 'air_temperature', 'dew_temperature', 'date', 'outlier', 'bdid_coded'))]\n\nhead(df)\n\nwrite.csv(df, \"train_clean_ds4_bin_oh.csv\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 4 Regression using XGBoost\n## 4.1 Load Data\nAgain, because of the RAM issue we choose to load the cleaned, downsampled, and encoded data in this part.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"rm(list=ls())\ngc()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = read.csv(\"../input/ds4-bin-oh/train_clean_ds4_bin_oh.csv\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 4.2 Split Training, Validation and Testing Datasets\nAs the data we use is only the training part in the original dataset (as kaggle does not provide the label for testing data), so we need to split the dataset into training, validation and testing parts.\n\nWe have tried to use the cross-validation in xgboost training, but again, there is no enough RAM for this even with the downsampled dataset. And we think there is no point in repeating the data cleaning and downsampling process again and again (we have done this time-consuming process for a lot of times) just to find a proper downsample rate to perform the follow-up modeling under the RAM limitation. So to examine the overfitting problem, we choose to split the dataset into training and testing dataset, and train the model with 5-fold cross-validation manuallly.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"test_ind <- sample(nrow(train_df), 0.1*nrow(train_df))\ndf.test <- train_df[test_ind, ]\ndf.train <- train_df[-test_ind, ]\n\nwrite.csv(df.train[, -1], \"split_train.csv\")\nwrite.csv(df.test[, -1], \"split_test.csv\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Split data for experiment with ordinal encoded features in experiment 4.7.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = read.csv(\"../input/ds4-ord/train_clean_ds4_ord.csv\")\n\nx.ord <- c(colnames(train_df))\ny.ord <- c(\"meter_reading\")\nx.ord <- x.ord[-c(1, 4)]\n\ndf.test.ord <- train_df[test_ind, ]\ndf.train.ord <- train_df[-test_ind, ]\nwrite.csv(df.test.ord[, -1], \"split_test_ord.csv\")\nwrite.csv(df.train.ord[, -1], \"split_train_ord.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"rm(list=ls())\ngc()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 4.3 Compose Input Datasets\nSelect corresponding columns and transform it into the type for training.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\n\nproc <- function(data, x) {\n  data %>% \n    select_at(x) %>% \n    data.matrix()\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 4.4 Training, Cross-validation and Testing\nBasically we are repeating the same training and testing process on each fold of data. Each time we use one fold as the validation dataset, and the other 4 folds as the training dataset, and after finish training we trained the model on the whole training set and test the model with the test set.\n\nAs we did many experiments, so for each model we just trained for 100 iterations (with an early stop iteration of 10).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df.train <- read.csv(\"../input/split-data/split_train.csv\")\ninds <- caret::createFolds(df.train$X, 5)\n\n# set the column names for input (x) and target (y) variables\nx <- c(colnames(df.train))\ny <- c(\"meter_reading\")\nx <- x[-c(1, 2)]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Fold 1:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\ntrain <- df.train[-inds[[1]], ]\nval <- df.train[inds[[1]], ]\ndtrain <- xgb.DMatrix(data = proc(train, x), label=proc(train, y))\ndval <- xgb.DMatrix(data = proc(val, x), label=proc(val, y))\n\ngc()\nwatchlist <- list(train=dtrain, test=dval)\nxgboost_model <- xgb.train(data=dtrain, max.depth=10, eta=0.3, nthread=3, nround=100, early_stopping_rounds=10,\n                           subsample=0.5, colsample_bytree=0.5, seed=1, objective=\"reg:squarederror\",watchlist=watchlist)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Fold 2:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\ntrain <- df.train[-inds[[2]], ]\nval <- df.train[inds[[2]], ]\ndtrain <- xgb.DMatrix(data = proc(train, x), label=proc(train, y))\ndval <- xgb.DMatrix(data = proc(val, x), label=proc(val, y))\n\ngc()\nwatchlist <- list(train=dtrain, test=dval)\nxgboost_model <- xgb.train(data=dtrain, max.depth=10, eta=0.3, nthread=3, nround=100, early_stopping_rounds=10,\n                           subsample=0.5, colsample_bytree=0.5, seed=1, objective=\"reg:squarederror\",watchlist=watchlist)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Fold 3:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\ntrain <- df.train[-inds[[3]], ]\nval <- df.train[inds[[3]], ]\ndtrain <- xgb.DMatrix(data = proc(train, x), label=proc(train, y))\ndval <- xgb.DMatrix(data = proc(val, x), label=proc(val, y))\n\ngc()\nwatchlist <- list(train=dtrain, test=dval)\nxgboost_model <- xgb.train(data=dtrain, max.depth=10, eta=0.3, nthread=3, nround=100, early_stopping_rounds=10,\n                           subsample=0.5, colsample_bytree=0.5, seed=1, objective=\"reg:squarederror\",watchlist=watchlist)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Fold 4:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\ntrain <- df.train[-inds[[4]], ]\nval <- df.train[inds[[4]], ]\ndtrain <- xgb.DMatrix(data = proc(train, x), label=proc(train, y))\ndval <- xgb.DMatrix(data = proc(val, x), label=proc(val, y))\n\ngc()\nwatchlist <- list(train=dtrain, test=dval)\nxgboost_model <- xgb.train(data=dtrain, max.depth=10, eta=0.3, nthread=3, nround=100, early_stopping_rounds=10,\n                           subsample=0.5, colsample_bytree=0.5, seed=1, objective=\"reg:squarederror\",watchlist=watchlist)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Fold 5:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\ntrain <- df.train[-inds[[4]], ]\nval <- df.train[inds[[4]], ]\ndtrain <- xgb.DMatrix(data = proc(train, x), label=proc(train, y))\ndval <- xgb.DMatrix(data = proc(val, x), label=proc(val, y))\n\ngc()\nwatchlist <- list(train=dtrain, test=dval)\nxgboost_model <- xgb.train(data=dtrain, max.depth=10, eta=0.3, nthread=3, nround=100, early_stopping_rounds=10,\n                           subsample=0.5, colsample_bytree=0.5, seed=1, objective=\"reg:squarederror\",watchlist=watchlist)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The validation loss and training loss are decreasing, there is no tendency of overfitting. And we can obtain the average relative loss reduced to around 0.28.\n### Final model\nThe final model is trained on the whole training dataset and tested on the test dataset.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"rm(list=c(\"train\", \"val\", \"dtrain\", \"dval\"))\ngc()\n\ndf.test <- read.csv(\"../input/split-data/split_test.csv\")\ndtrain <- xgb.DMatrix(data = proc(df.train, x), label=proc(df.train, y))\ndtest <- xgb.DMatrix(data = proc(df.test, x), label=proc(df.test, y))\n\nwatchlist <- list(train=dtrain, test=dtest)\nxgboost_model <- xgb.train(data=dtrain, max.depth=10, eta=0.3, nthread=3, nround=100, early_stopping_rounds=10,\n                           subsample=0.5, colsample_bytree=0.5, seed=1, objective=\"reg:squarederror\", watchlist=watchlist)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 4.5 Feature Importance","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\nimportance_matrix <- xgb.importance(model = xgboost_model)\nxgb.plot.importance(importance_matrix)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 4.6 Colinearity Analysis (simple comparison)\nAs we know in linear models the colinearity of input variables will cause problems in regression. To investigate the effects of colinearity (as we have dew_temperature and air_temperature, the correlation of which is about 0.7) in the tree model, here we did a trivial experiment on this. \n\nWe manually added one column in the data which is linear to air_temperature, and did the regression on the dataset.\n\nAs the importance of features may also affect the comparison, we also did another experiment with one added column linear to square_feet to see the effects of colinearity.\n\nBecause it is just for a simple comparison, we didn't follow the strict 5-fold cross-validation process but directly trained on the whole training dataset and tested on the test dataset.\n### With_fake air_temperature","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"rm(list=c(\"xgboost_model\"))\ngc()\n\nx <- append(x, \"fake_air_temp\", 0)\ndf.train$fake_air_temp <- 1.5*df.train$air_temp_clean\ndf.test$fake_air_temp <- 1.5*df.test$air_temp_clean\n\ngc()\n\ndtrain <- xgb.DMatrix(data = proc(df.train, x), label=proc(df.train, y))\ndval <- xgb.DMatrix(data = proc(df.test, x), label=proc(df.test, y))\n\nwatchlist <- list(train=dtrain, test=dval)\nxgboost_model <- xgb.train(data=dtrain, max.depth=10, eta=0.3, nthread=3, nround=100, early_stopping_rounds=10,\n                           subsample=0.5, colsample_bytree=0.5, seed=1, objective=\"reg:squarederror\", watchlist=watchlist)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### With fake_square_feet","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"rm(xgboost_model)\ngc()\n\nx <- x[-1]\nx <- append(x, \"fake_square_feet\", 0)\ndf.train <- subset(df.train, select = -fake_air_temp )\ndf.test <- subset(df.test, select = -fake_air_temp )\ndf.train$fake_square_feet <- 1.5*df.train$square_feet\ndf.test$fake_square_feet <- 1.5*df.test$square_feet\n\ngc()\ndtrain <- xgb.DMatrix(data = proc(df.train, x), label=proc(df.train, y))\ndval <- xgb.DMatrix(data = proc(df.test, x), label=proc(df.test, y))\n\nwatchlist <- list(train=dtrain, test=dval)\nxgboost_model <- xgb.train(data=dtrain, max.depth=10, eta=0.3, nthread=3, nround=100, early_stopping_rounds=10,\n                           subsample=0.5, colsample_bytree=0.5, seed=1, objective=\"reg:squarederror\", watchlist=watchlist)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There is no big difference when training with collinear variables. We think it is because the tree models choose one feature to split at a time using the purity score, and when two variables have strong colinearity with each other, the changes in purity of choosing them will be similar, which means the purity measure are not sensible to their collinearity. So they are actually playing similar functions as one variable.\n## 4.7 Modeling with Ordinal Encoded Features\nWe have also tried to use the ordinal encoding method for the four categorical features (building_id, site_id, meter, primary_use). And the result is only slightly worse than the one-hot encoding in terms of the relative training and validation losses. We think this may because these features are not the most critical features to meter_reading, as can be seen from the feature importance matrix above. \n### Fold 1","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"rm(list=c(\"xgboost_model\", \"df.train\", \"df.test\", \"dtrain\", \"dtest\", \"dval\"))\ngc()\n\ndf.train <- read.csv(\"../input/split-data-ord/split_train_ord.csv\")\ninds <- caret::createFolds(df.train$X, 5)\n\nx <- c(colnames(df.train))\ny <- c(\"meter_reading\")\nx <- x[-c(1, 4)]\n\nval <- df.train[inds[[1]], ]\ntrain <- df.train[-inds[[1]], ]\n              \ndtrain <- xgb.DMatrix(data = proc(train, x), label=proc(train, y))\ndval <- xgb.DMatrix(data = proc(val, x), label=proc(val, y))\n\ngc()\nwatchlist <- list(train=dtrain, test=dval)\nxgb_model_ord <- xgb.train(data=dtrain, max.depth=10, eta=0.3, nthread=3, nround=100, \n                           subsample=0.5, colsample_bytree=0.5, seed=1, early_stopping_rounds=10, objective=\"reg:squarederror\",watchlist=watchlist)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Fold 2","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\nval <- df.train[inds[[2]], ]\ntrain <- df.train[-inds[[2]], ]\ndtrain <- xgb.DMatrix(data = proc(train, x), label=proc(train, y))\ndval <- xgb.DMatrix(data = proc(val, x), label=proc(val, y))\n\ngc()\nwatchlist <- list(train=dtrain, test=dval)\nxgb_model_ord <- xgb.train(data=dtrain, max.depth=10, eta=0.3, nthread=3, nround=100, \n                           subsample=0.5, colsample_bytree=0.5, seed=1, early_stopping_rounds=10, objective=\"reg:squarederror\",watchlist=watchlist)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Fold 3","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\nval <- df.train[inds[[3]], ]\ntrain <- df.train[-inds[[3]], ]\ndtrain <- xgb.DMatrix(data = proc(train, x), label=proc(train, y))\ndval <- xgb.DMatrix(data = proc(val, x), label=proc(val, y))\n\ngc()\nwatchlist <- list(train=dtrain, test=dval)\nxgb_model_ord <- xgb.train(data=dtrain, max.depth=10, eta=0.3, nthread=3, nround=100, \n                           subsample=0.5, colsample_bytree=0.5, seed=1, early_stopping_rounds=10, objective=\"reg:squarederror\",watchlist=watchlist)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Fold 4","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\nval <- df.train[inds[[4]], ]\ntrain <- df.train[-inds[[4]], ]\ndtrain <- xgb.DMatrix(data = proc(train, x), label=proc(train, y))\ndval <- xgb.DMatrix(data = proc(val, x), label=proc(val, y))\n\ngc()\nwatchlist <- list(train=dtrain, test=dval)\nxgb_model_ord <- xgb.train(data=dtrain, max.depth=10, eta=0.3, nthread=3, nround=100, \n                           subsample=0.5, colsample_bytree=0.5, seed=1, early_stopping_rounds=10, objective=\"reg:squarederror\",watchlist=watchlist)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Fold 5","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc()\nval <- df.train[inds[[5]], ]\ntrain <- df.train[-inds[[5]], ]\ndtrain <- xgb.DMatrix(data = proc(train, x), label=proc(train, y))\ndval <- xgb.DMatrix(data = proc(val, x), label=proc(val, y))\n\ngc()\nwatchlist <- list(train=dtrain, test=dval)\nxgb_model_ord <- xgb.train(data=dtrain, max.depth=10, eta=0.3, nthread=3, nround=100, \n                           subsample=0.5, colsample_bytree=0.5, seed=1, early_stopping_rounds=10, objective=\"reg:squarederror\",watchlist=watchlist)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 4.8 Regression Tree model\nHere we also include our former results from a simple regression tree model. And we also use 5-fold cross-validation to evaluate its performance.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"rm(list=c(\"df.train\", \"df.test\", \"dtrain\", \"dval\", \"dtest\", \"xgb_model_ord\"))\ngc()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"library(rpart)\nlibrary(rpart.plot)\nset.seed(15)\n\ndf.train <- read.csv(\"../input/split-data/split_train.csv\")\ndf.test <- read.csv(\"../input/split-data/split_test.csv\")\ndf.train <- df.train[, -1]\n\ntreeFitted <- rpart(meter_reading~., \n                    data = df.train, \n                    method='anova', \n                    control = rpart.control(\n                      xval = 5, # number of cross-validations\n                      minbucket = 100, # the minimum number of observations in any terminal (leaf) node\n                      cp = 0.0006)\n                    )\n\nprintcp(treeFitted)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We selected cp with the 1-se rule. The relative loss reduced to ~0.5. But we can see the training is very time-consuming and the tree has already grows very big.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=8)\nrpart.plot(treeFitted,type=3)\nplotcp(treeFitted)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.4.2","file_extension":".r","codemirror_mode":"r"}},"nbformat":4,"nbformat_minor":4}