{"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":"markdown","source":"References - \n* [Steering Wheel of Fortune - Porto Seguro EDA](https://www.kaggle.com/code/headsortails/steering-wheel-of-fortune-porto-seguro-eda/report)\n* [Clustered bar plot in R using ggplot2](https://stackoverflow.com/questions/17273949/clustered-bar-plot-in-r-using-ggplot2)\n* [Machine Learning with XGBoost (in R)](https://www.kaggle.com/code/rtatman/machine-learning-with-xgboost-in-r/notebook)\n* [How to one hot encode several categorical variables in R](https://stackoverflow.com/questions/48649443/how-to-one-hot-encode-several-categorical-variables-in-r)\n* [Increase plot width ggplot2](https://stackoverflow.com/questions/29587881/increase-plot-size-width-in-ggplot2)\n* [Cyclical features encoding](https://towardsdatascience.com/cyclical-features-encoding-its-about-time-ce23581845ca)\n* [xgboost parameters](https://xgboost.readthedocs.io/en/latest/parameter.html)\n* [How to Convert Multiple Columns to Numeric Using dplyr](https://www.statology.org/convert-multiple-columns-to-numeric-dplyr/)","metadata":{}},{"cell_type":"markdown","source":"Load required libraries","metadata":{}},{"cell_type":"code","source":"# general visualisation\nsuppressMessages(library('ggplot2')) # visualisation\nsuppressMessages(library('scales')) # visualisation\nsuppressMessages(library('grid')) # visualisation\nsuppressMessages(library('ggthemes')) # visualisation\nsuppressMessages(library('gridExtra')) # visualisation\nsuppressMessages(library('RColorBrewer')) # visualisation\nsuppressMessages(library('corrplot')) # visualisation\n\n# general data manipulation\nsuppressMessages(library('dplyr')) # data manipulation\nsuppressMessages(library('readr')) # input/output\nsuppressMessages(library('data.table')) # data manipulation\nsuppressMessages(library('tibble')) # data wrangling\nsuppressMessages(library('tidyr')) # data wrangling\nsuppressMessages(library('stringr')) # string manipulation\nsuppressMessages(library('forcats')) # factor manipulation\nsuppressMessages(library('rlang')) # data manipulation\nsuppressMessages(library('moments')) # data manipulation\n\n# specific visualisation\nsuppressMessages(library('alluvial')) # visualisation\n#library('ggfortify') # visualisation\nsuppressMessages(library('ggrepel')) # visualisation\nsuppressMessages(library('ggridges')) # visualisation\n#library('VIM') # NAs\nsuppressMessages(library('plotly')) # interactive\nsuppressMessages(library('ggforce')) # visualisation\n\n# modelling\nsuppressMessages(library('xgboost')) # modelling\nsuppressMessages(library('caret')) # modelling\nsuppressMessages(library('MLmetrics')) # gini metric","metadata":{"_uuid":"a2d42175-0b90-4545-806e-1fef8dcdfc76","_cell_guid":"28a1ff4c-3fa3-4293-8486-d206e1f186b5","collapsed":false,"jupyter":{"outputs_hidden":false},"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-05-25T06:47:59.377740Z","iopub.execute_input":"2022-05-25T06:47:59.380526Z","iopub.status.idle":"2022-05-25T06:48:03.444858Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# list.files('../input/energy-anomaly-detection/')\ntrain_features = fread('../input/energy-anomaly-detection/train_features.csv',stringsAsFactors = F)\ntest_features = fread('../input/energy-anomaly-detection/test_features.csv',stringsAsFactors = F)","metadata":{"_uuid":"3bd391f8-c0eb-4325-a36d-c3ea65a2931e","_cell_guid":"39d6b603-a890-4834-b74e-689360cdbdca","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2022-05-25T06:51:57.016532Z","iopub.execute_input":"2022-05-25T06:51:57.018250Z","iopub.status.idle":"2022-05-25T06:52:05.445721Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# helper functions\nget_binCI <- function(x,n) as.list(setNames(binom.test(x,n)$conf.int, c(\"lwr\", \"upr\")))\nmultiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {\n  plots <- c(list(...), plotlist)\n  numPlots = length(plots)\n  if (is.null(layout)) {\n    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),\n                    ncol = cols, nrow = ceiling(numPlots/cols))}\n if (numPlots==1) {\n    print(plots[[1]])\n  } else {\n    grid.newpage()\n    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))\n    for (i in 1:numPlots) {\n      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))\n      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,\n                                      layout.pos.col = matchidx$col))}}}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-05-25T06:48:19.705099Z","iopub.execute_input":"2022-05-25T06:48:19.707038Z","iopub.status.idle":"2022-05-25T06:48:19.723663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Univariates\n### Cat variables\n* building_id : near uniform signal count from all buildings\n* site id : skew observed in count of singals per site id, need to understand what site id means\n* primary_use : education and office appear to be generating highest signal count\n* year_built : year built is difference of actual year built and 1900, it's unclear why 255 has a spike\n* floor_count : most buildings appear to having 0 floors? Only ground floor? Wouldn't office and education imply more than one floor per bulding?\n* cloud coverage : 255 appears to be a value filled in case of NAs, 0 is mode, might imply sunny days with clear skies == more power consumption?\n* hour : even number of signals generated each hour\n* weekday : equal amount of signals on all weekdays\n* month : equal amount of singnals on all months\n* features ending with '_x' and '_y' are cyclic transformations of the timestamp field, one of the appraoches used to retain timestamp information as distinct features\n\n","metadata":{}},{"cell_type":"code","source":"# univariates\n############# cat variables #################\noptions(repr.plot.width=25, repr.plot.height=50)\n# building_id : near uniform signal count from all buildings\np1 <- train_features %>%\n    mutate(bulding_id = as.factor(building_id)) %>%\n    ggplot(aes(building_id,fill=(building_id))) + \n    geom_bar() +\n    theme(legend.position = \"none\")\n\n# site id : skew observed in count of singals per site id, need to understand what site id means\np2 <- train_features %>%\n    mutate(site_id = as.factor(site_id)) %>%\n    ggplot(aes(site_id,fill=(site_id))) + \n    geom_bar() +\n    scale_y_log10() +\n    theme(legend.position = \"none\")\n\n# primary_use : education and office appear to be generating highest signal count\np3 <- train_features %>%\n    ggplot(aes(primary_use,fill=primary_use)) +\n    geom_bar() +\n    theme(text = element_text(size = 15)) +\n    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +\n    theme(legend.position = \"none\")\n\n# year_built : year built is difference of actual year built and 1900, it's unclear why 255 has a spike\np4 <- train_features %>%\n    mutate(year_built = as.factor(year_built)) %>%\n    ggplot(aes(year_built,fill=(year_built))) +\n    geom_bar() + \n    theme(legend.position = \"none\") +\n    scale_y_log10()\n\n# floor_count : most buildings appear to having 0 floors? Only ground floor? Wouldn't office and education imply more than one floor per bulding?\np5 <- train_features %>%\n    ggplot(aes(as.factor(floor_count),fill=as.factor(floor_count))) + \n    geom_bar() + \n    theme(legend.position = \"none\") +\n    scale_y_log10()\n\n# cloud coverage : 255 appears to be a value filled in case of NAs, 0 is mode, might imply sunny days with clear skies == more power consumption?\np6 <- train_features %>%\n    mutate(cloud_coverage = as.factor(cloud_coverage)) %>%\n    ggplot(aes(cloud_coverage,fill = cloud_coverage)) +\n    theme(legend.position = \"none\") +\n    geom_bar()\n\n# hour : even number of signals generated each hour\np7 <- train_features %>%\n    mutate(hour = as.factor(hour)) %>%\n    ggplot(aes(hour,fill = hour)) +\n    theme(legend.position = \"none\") +\n    geom_bar()\n\n# weekday : equal amount of signals on all weekdays\np8 <- train_features %>%\n    mutate(weekday = as.factor(weekday)) %>%\n    ggplot(aes(weekday,fill=weekday)) +\n    theme(legend.position = \"none\") +\n    geom_bar()\n\n# month : equal amount of singnals on all months\np9 <- train_features %>%\n    mutate(month = as.factor(month)) %>%\n    ggplot(aes(month,fill=month)) +\n    theme(legend.position = \"none\") +\n    geom_bar()\n\n#! all features starting with gte are calcualted using meter field, it's unclear how meter field is calculated\n# _x and _y are cyclic transformations of the timestamp field, one of the appraoches used to retain timestamp information as distinct features\n\n# is_holiday\np10 <- train_features %>%\n    mutate(is_holiday = as.factor(is_holiday)) %>%\n    ggplot(aes(is_holiday,fill=is_holiday)) +\n    theme(legend.position = \"none\") +\n    geom_bar()\n\nlayout <- matrix(c(1,2,3,4,5,6,7,8,9),9,1,byrow=TRUE)\nmultiplot(p2,p3,p4,p5,p6,p7,p8,p9,p10,layout = layout)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-05-25T07:12:10.804967Z","iopub.execute_input":"2022-05-25T07:12:10.807271Z","iopub.status.idle":"2022-05-25T07:12:23.501933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Univariates\n### Numeric variables\n* square feet : highly right skewed, checking for relation with primary use\n* temprature : normally distributed\n* Q. is cloud coverage related with temprature?\n* **Kaggle has inbuilt dataset visualizer for univariates, will refer that instead going forward!!**","metadata":{}},{"cell_type":"code","source":"############# numeric variables #################\noptions(repr.plot.width=25, repr.plot.height=12)\n# square feet : highly right skewed, checking for relation with primary use\np1 <- train_features %>%\n    #mutate(square_feet = log(square_feet)) %>%\n    ggplot(aes(square_feet)) +\n    geom_histogram(bins = 50)\np2 <- train_features %>%\n    mutate(square_feet = log(square_feet)) %>%\n    ggplot(aes(square_feet,fill=as.factor(primary_use))) + \n    geom_boxplot()\n\n# temprature : normally distributed\np3 <- train_features %>%\n    ggplot(aes(air_temperature)) +\n    geom_histogram(bins = 20)\n# is cloud coverage related with temprature?\np4 <- train_features %>%\n    filter(cloud_coverage != 255) %>%\n    #filter(cloud_coverage == 4) %>%\n    ggplot(aes(air_temperature,fill = as.factor(cloud_coverage))) +\n    geom_density(alpha = 0.9, bw = 1)\n\n# ummm, will refer to the dataset visualizer in Kaggle for further univariates!\nlayout <- matrix(c(1,2,3,4),2,2,byrow=TRUE)\nmultiplot(p1,p2,p3,p4, layout = layout)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-05-25T06:56:43.396743Z","iopub.execute_input":"2022-05-25T06:56:43.398515Z","iopub.status.idle":"2022-05-25T06:56:50.200450Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Bivariates\n### Cat variables\n* building_id : significant seperation in anomaly% building wise | contended for multivariate deep dive\n* site_id : significant seperateion in anomaly% building wise | contended for multivariate deep dive\n* primary_use : education and office have the highest signal count, although the anomaly% for both is in similar range. Other's have a high anomaly% but less #of records\n* year_built : interestingly buildings built in 1965 have very high anomaly%, why would that be the case?\n* floor_count : floor count x primary_use could be a new variable\n* cloud_coverage : cloud coverage doesn't seem to have significant difference across levels, although a sign. portion of cloud coverage column has NULLs, candidate for filling NAs\n* hour : no significant seperation, this might need to be split by month\n* July seems to have highest anomalies while April, the least, sig. seperation in anomaly% by month, month-hr shows even stronger seperation\n* Monday has lower anomaly%, Sunday evenings have a higher anomaly% | weekday_hr shows good separetion in anomaly\n* list of new features to be created : \n    * month_hour, month_weekday_hour\n    * building_month_hour, building_month_weekday_hr","metadata":{}},{"cell_type":"code","source":"########## bivariates ###########\n# cat columns\noptions(repr.plot.width=25, repr.plot.height=50)\n# building_id : significant seperation in anomaly% building wise | contended for multivariate deep dive\np1 <- train_features %>%\n    mutate(building_id = as.factor(building_id)) %>%\n    group_by(building_id,anomaly) %>%\n    count() %>%\n    spread(anomaly, n) %>%\n    mutate(anomaly_pct = `1`/(`1` + `0`)*100,\n         lwr = get_binCI(`1`,(`1`+`0`))[[1]]*100,\n         upr = get_binCI(`1`,(`1`+`0`))[[2]]*100\n         ) %>%\n    ggplot(aes(building_id,anomaly_pct,fill=building_id)) +\n    geom_col(width = 0.8) +\n    geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.5, size = 0.7, color = \"gray30\") +\n    theme(legend.position = \"none\") +\n    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +\n    labs(y = \"Anomaly[%]\")\n# site_id : significant seperateion in anomaly% building wise | contended for multivariate deep dive\np2 <- train_features %>%\n    mutate(site_id = as.factor(site_id)) %>%\n    group_by(site_id,anomaly) %>%\n    count() %>%\n    spread(anomaly, n) %>%\n    mutate(anomaly_pct = `1`/(`1` + `0`)*100,\n         lwr = get_binCI(`1`,(`1`+`0`))[[1]]*100,\n         upr = get_binCI(`1`,(`1`+`0`))[[2]]*100\n         ) %>%\n    ggplot(aes(site_id,anomaly_pct,fill=site_id)) +\n    geom_col(width = 0.8) +\n    geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.1, size = 1, color = \"gray30\") +\n    theme(legend.position = \"none\") +\n    #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +\n    labs(y = \"Anomaly[%]\")\n# primary_use : education and office have the highest signal count, although the anomaly% for both is in similar range. Other's have a high anomaly% but less #of records\np3 <- train_features %>%\n    mutate(primary_use = as.factor(primary_use)) %>%\n    group_by(primary_use,anomaly) %>%\n    count() %>%\n    spread(anomaly, n) %>%\n    mutate(anomaly_pct = `1`/(`1` + `0`)*100,\n         lwr = get_binCI(`1`,(`1`+`0`))[[1]]*100,\n         upr = get_binCI(`1`,(`1`+`0`))[[2]]*100\n         ) %>%\n    ggplot(aes(primary_use,anomaly_pct,fill=primary_use)) +\n    geom_col(width = 0.8) +\n    geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.1, size = 1, color = \"gray30\") +\n    theme(legend.position = \"none\") +\n    #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +\n    labs(y = \"Anomaly[%]\")\n# year_built : interestingly buildings built in 1965 have very high anomaly%, why would that be the case?\np4 <- train_features %>%\n    mutate(year_built = as.factor(year_built)) %>%\n    group_by(year_built,anomaly) %>%\n    count() %>%\n    spread(anomaly, n) %>%\n    mutate(anomaly_pct = `1`/(`1` + `0`)*100,\n         lwr = get_binCI(`1`,(`1`+`0`))[[1]]*100,\n         upr = get_binCI(`1`,(`1`+`0`))[[2]]*100\n         ) %>%\n    ggplot(aes(year_built,anomaly_pct,fill=year_built)) +\n    geom_col(width = 0.8) +\n    geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.1, size = 1, color = \"gray30\") +\n    theme(legend.position = \"none\") +\n    #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +\n    labs(y = \"Anomaly[%]\") + ggtitle(label = 'Anomaly% by year_builtbuilt_year')\n# floor_count : floor count x primary_use could be a new variable\np5 <- train_features %>%\n    mutate(floor_count = as.factor(floor_count)) %>%\n    group_by(floor_count,anomaly) %>%\n    count() %>%\n    spread(anomaly, n) %>%\n    mutate(anomaly_pct = `1`/(`1` + `0`)*100,\n         lwr = get_binCI(`1`,(`1`+`0`))[[1]]*100,\n         upr = get_binCI(`1`,(`1`+`0`))[[2]]*100\n         ) %>%\n    ggplot(aes(floor_count,anomaly_pct,fill=as.factor(floor_count))) +\n    geom_col(width = 0.8,position = 'dodge') +\n    geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.1, size = 1, color = \"gray30\") +\n    theme(legend.position = \"right\") +\n    #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +\n    labs(y = \"Anomaly[%]\")\n# cloud_coverage : cloud coverage doesn't seem to have significant difference across levels, although a sign. portion of cloud coverage column has NULLs, candidate for filling NAs\np6 <- train_features %>%\n    mutate(cloud_coverage = as.factor(cloud_coverage)) %>%\n    group_by(cloud_coverage,anomaly) %>%\n    count() %>%\n    spread(anomaly, n) %>%\n    mutate(anomaly_pct = `1`/(`1` + `0`)*100,\n         lwr = get_binCI(`1`,(`1`+`0`))[[1]]*100,\n         upr = get_binCI(`1`,(`1`+`0`))[[2]]*100\n         ) %>%\n    ggplot(aes(cloud_coverage,anomaly_pct,fill=as.factor(cloud_coverage))) +\n    geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.1, size = 1, color = \"gray30\") +\n    theme(legend.position = \"none\") +\n    #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +\n    labs(y = \"Anomaly[%]\")\n# hour : no significant seperation, this might need to be split by month\np7 <- train_features %>%\n    mutate(hour = as.factor(hour)) %>%\n    group_by(hour,anomaly) %>%\n    count() %>%\n    spread(anomaly, n) %>%\n    mutate(anomaly_pct = `1`/(`1` + `0`)*100,\n         lwr = get_binCI(`1`,(`1`+`0`))[[1]]*100,\n         upr = get_binCI(`1`,(`1`+`0`))[[2]]*100\n         ) %>%\n    ggplot(aes(hour,anomaly_pct,fill=as.factor(hour))) +\n    geom_col(width = 0.8,position = 'dodge') +\n    geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.1, size = 1, color = \"gray30\") +\n    theme(legend.position = \"none\") +\n    #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +\n    labs(y = \"Anomaly[%]\")\n# July seems to have highest anomalies while April, the least, sig. seperation in anomaly% by month, month-hr shows even stronger seperation\np8 <- train_features %>%\n    mutate(month = as.factor(month),\n          month_hr = as.factor(paste(month,hour,sep = \"_\"))) %>%\n    group_by(month_hr,anomaly) %>%\n    count() %>%\n    spread(anomaly, n) %>%\n    mutate(anomaly_pct = `1`/(`1` + `0`)*100,\n         lwr = get_binCI(`1`,(`1`+`0`))[[1]]*100,\n         upr = get_binCI(`1`,(`1`+`0`))[[2]]*100\n         ) %>%\n    ggplot(aes(month_hr,anomaly_pct,fill=as.factor(month_hr))) +\n    geom_col(width = 0.8,position = 'dodge') +\n    geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.1, size = 1, color = \"gray30\") +\n    theme(legend.position = \"none\") +\n    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +\n    labs(y = \"Anomaly[%]\")\n# Monday has lower anomaly%, Sunday evenings have a higher anomaly% | weekday_hr shows good separetion in anomaly\np9 <- train_features %>%\n    mutate(weekday = as.factor(weekday),\n          weekday_hr = as.factor(paste(weekday,hour,sep = \"_\"))) %>%\n    group_by(weekday_hour,anomaly) %>%\n    count() %>%\n    spread(anomaly, n) %>%\n    mutate(anomaly_pct = `1`/(`1` + `0`)*100,\n         lwr = get_binCI(`1`,(`1`+`0`))[[1]]*100,\n         upr = get_binCI(`1`,(`1`+`0`))[[2]]*100\n         ) %>%\n    ggplot(aes(weekday_hour,anomaly_pct,fill=as.factor(weekday_hour))) +\n    geom_col(width = 0.8,position = 'dodge') +\n    geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.1, size = 1, color = \"gray30\") +\n    theme(legend.position = \"none\") +\n    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +\n    labs(y = \"Anomaly[%]\")\n# list of new features to be created : \n    # month_hour, month_weekday_hour\n    # building_month_hour, building_month_weekday_hr\nlayout <- matrix(c(1,2,3,4,5,6,7,8,9),9,1,byrow=TRUE)\nmultiplot(p1,p2,p3,p4,p5,p6,p7,p8,p9,layout = layout)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-05-25T07:13:09.905404Z","iopub.execute_input":"2022-05-25T07:13:09.907278Z","iopub.status.idle":"2022-05-25T07:13:18.258809Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Bivariates\n### Numeric variables\n* meter_reading : strong seperation observed\n* square_feet : we don't observe a strong seperation\n* air_temperature : seperation observed between 10-20 degrees - how to use this info for FE\n    * is there a relationship between temprature trend and anomaly appearance | how is a relevance established here?","metadata":{}},{"cell_type":"code","source":"# numeric columns\n# meter_reading : strong seperation observed\noptions(repr.plot.width=25, repr.plot.height=15)\np1 <- train_features %>%\n  mutate(anomaly = as.factor(anomaly)) %>%\n  ggplot(aes(log(meter_reading), fill = as.factor(anomaly))) +\n  geom_density(alpha = 0.5, bw = 0.1)+\n  theme(legend.position = \"right\")\n# square_feet : we don't observe a strong seperation\np2 <- train_features %>%\n    ggplot(aes(log(square_feet),fill = as.factor(anomaly))) + \n    geom_density(alpha = 0.5,bw = 'nrd0')\n# air_temperature : seperation observed between 10-20 degrees\np3 <- train_features %>%\n    ggplot(aes(air_temperature,fill = as.factor(anomaly))) + \n    geom_density(alpha = 0.5, bw = 'nrd0') + \n    scale_x_continuous(breaks = seq(-50, 50, 5))\np31 <- train_features %>%\n    filter(building_id == 1319) %>%\n    ggplot(aes(as.factor(month),air_temperature,fill=as.factor(anomaly))) +\n    geom_boxplot(position = 'dodge')\n# is there a relationship between temprature trend and anomaly appearance | how is a relevance established here?\nhighlight_df <- train_features %>% \n    filter(building_id == 1319 &\n             anomaly == 1) %>%\n    filter(as.IDate(timestamp) >= '2016-09-01' &\n           as.IDate(timestamp) <= '2016-10-31')\np32 <- train_features %>%\n    filter(building_id == 1319) %>%\n    filter(as.IDate(timestamp) >= '2016-09-01' &\n           as.IDate(timestamp) <= '2016-10-31') %>%\n    ggplot(aes(timestamp,air_temperature_max_lag73)) +\n    geom_line() +\n    geom_point(data = highlight_df,\n               color = 'red') +\n    geom_line(aes(y=air_temperature)) +\n    geom_line(aes(y=dew_temperature))\n# dew temprature : slight seperation observed, anomaly concentrated on lower and higher end of tempratures, similar to air temprature\np4 <- train_features %>%\n    ggplot(aes(x=dew_temperature,fill = as.factor(anomaly))) + \n    geom_density(alpha = 0.5)\n# precip_depth_1_hr : highly right skewed, most cases have 0 precip. not much seperation being observed\np5 <- train_features %>%\n    ggplot(aes(precip_depth_1_hr,fill = as.factor(anomaly))) + \n    geom_density(alpha = 0.5)\n# sea_level_pressure : no sig. seperation\np6 <- train_features %>%\n    ggplot(aes(sea_level_pressure,fill = as.factor(anomaly))) + \n    geom_density(alpha = 0.5)\n# wind_speed : no sig. seperation\np7 <- train_features %>%\n    ggplot(aes(wind_speed,fill = as.factor(anomaly))) + \n    geom_density()\nlayout <- matrix(c(1,2,3,4,5,6,5,6,7,8),4,2,byrow=TRUE)\nmultiplot(p1,p2,p3,p4,p5,p6,p7,cols = 2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-05-25T07:13:28.312372Z","iopub.execute_input":"2022-05-25T07:13:28.314638Z","iopub.status.idle":"2022-05-25T07:13:36.843767Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# multivariate\n# year_built : interestingly buildings built in 1965 have very high anomaly%, why would that be the case? | floor count x primary_use could be a new variable\np5 <- train_features %>%\n    mutate(floor_count = as.factor(floor_count)) %>%\n    group_by(floor_count,primary_use,anomaly) %>%\n    count() %>%\n    spread(anomaly, n) %>%\n    mutate(anomaly_pct = `1`/(`1` + `0`)*100,\n         lwr = get_binCI(`1`,(`1`+`0`))[[1]]*100,\n         upr = get_binCI(`1`,(`1`+`0`))[[2]]*100\n         ) %>%\n    ggplot(aes(floor_count,anomaly_pct,fill=as.factor(primary_use))) +\n    geom_col(width = 0.8,position = 'dodge') +\n    geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.1, size = 1, color = \"gray30\") +\n    theme(legend.position = \"right\") +\n    #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +\n    labs(y = \"Anomaly[%]\")","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-05-25T06:49:07.319658Z","iopub.execute_input":"2022-05-25T06:49:07.321290Z","iopub.status.idle":"2022-05-25T06:49:07.735163Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#features with unknown relevance\nunknown_relevance_columns = c('building_weekday_hour',\t'building_weekday',\t'building_month',\t'building_hour',\t'building_meter',\t'gte_weekday',\t'gte_month',\t'gte_building_id',\t'gte_primary_use',\t'gte_site_id',\t'gte_meter',\t'gte_meter_hour',\t'gte_meter_weekday',\t'gte_meter_month',\t'gte_meter_building_id',\t'gte_meter_site_id',\t'gte_meter_primary_use',\t'gte_meter_building_id_hour',\t'gte_meter_building_id_weekday',\t'gte_meter_building_id_month')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-05-25T06:49:07.737839Z","iopub.execute_input":"2022-05-25T06:49:07.739378Z","iopub.status.idle":"2022-05-25T06:49:07.751999Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train <- train_features %>%\n    select(setdiff(colnames(train_features),c(unknown_relevance_columns,'year')))\ntest <- test_features %>%\n    select(setdiff(colnames(test_features),c(unknown_relevance_columns,'year','row_id')))","metadata":{"execution":{"iopub.status.busy":"2022-05-25T06:49:07.754637Z","iopub.execute_input":"2022-05-25T06:49:07.756104Z","iopub.status.idle":"2022-05-25T06:49:07.781234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"features_train_numeric <- train %>%\n    select(-c(anomaly,timestamp,hour,weekday,month,is_holiday)) %>%# remove non-numeric and non-relevant features\n    select_if(is.numeric)# select numeric features only\nfeatures_train_cat <- train %>%\n    select(setdiff(colnames(train),colnames(features_train_numeric))) %>%\n    select(-c(timestamp,anomaly)) %>%\n    mutate_all(as.factor)\ndmy <- dummyVars(\" ~ .\", data = features_train_cat)\n#trsf <- data.frame(predict(dmy, newdata = features_train_cat))\nfeatures_train <- cbind(features_train_numeric,features_train_cat)\n#write.csv(features_train,'features_train.csv',row.names = F)\nrm(features_train_numeric,features_train_cat)\ngc()","metadata":{"execution":{"iopub.status.busy":"2022-05-25T06:49:07.783908Z","iopub.execute_input":"2022-05-25T06:49:07.785388Z","iopub.status.idle":"2022-05-25T06:49:08.754902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"features_train <- data.matrix(features_train)\n# create train and validation sets\n# get the numb 70/30 training test split\nnumberOfTrainingSamples <- round(nrow(features_train) * .7)\n\n# training data\ntrain_data <- features_train[1:numberOfTrainingSamples,]\ntrain_labels <- train$anomaly[1:numberOfTrainingSamples]\n\n# testing data\ntest_data <- features_train[-(1:numberOfTrainingSamples),]\ntest_labels <- train$anomaly[-(1:numberOfTrainingSamples)]\n\n# put our testing & training data into two seperates Dmatrixs objects\ndtrain <- xgb.DMatrix(data = train_data, label= train_labels)\ndtest <- xgb.DMatrix(data = test_data, label= test_labels)","metadata":{"execution":{"iopub.status.busy":"2022-05-25T06:49:08.757790Z","iopub.execute_input":"2022-05-25T06:49:08.759416Z","iopub.status.idle":"2022-05-25T06:49:11.994380Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train a model using our training data\nmodel <- xgboost(data = dtrain, # the data   \n                 max_depth = 10,\n                 nround = 10, # max number of boosting iterations\n                 objective = \"binary:logistic\")  # the objective function","metadata":{"execution":{"iopub.status.busy":"2022-05-25T06:49:11.997705Z","iopub.execute_input":"2022-05-25T06:49:11.999658Z","iopub.status.idle":"2022-05-25T06:49:38.104196Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# generate predictions for our held-out testing data\npred <- predict(model, dtest)\n\n# get & print the classification error\nerr <- mean(as.numeric(pred > 0.5) != test_labels)\nprint(paste(\"test-error=\", err))","metadata":{"execution":{"iopub.status.busy":"2022-05-25T06:49:38.108010Z","iopub.execute_input":"2022-05-25T06:49:38.109905Z","iopub.status.idle":"2022-05-25T06:49:38.257632Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# variable importance plot\nimportance_matrix <- xgb.importance(names(features_train), model = model)\nxgb.plot.importance(importance_matrix, )","metadata":{"execution":{"iopub.status.busy":"2022-05-25T06:49:38.260396Z","iopub.execute_input":"2022-05-25T06:49:38.261784Z","iopub.status.idle":"2022-05-25T06:49:38.438869Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# clear memory\nls()\nrm(dtest,dtrain,features_train_cat,features_train_numeric,importance_matrix,numberOfTrainingSamples,p1,test_data,test_labels,train,train_data,train_features,train_labels)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-05-25T06:49:38.441509Z","iopub.execute_input":"2022-05-25T06:49:38.443011Z","iopub.status.idle":"2022-05-25T06:49:38.838273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# generate and submit predictions \nfeatures_test_numeric <- test %>%\n    select(-c(timestamp,hour,weekday,month,is_holiday)) %>%# remove non-numeric and non-relevant features\n    select_if(is.numeric)# select numeric features only\nfeatures_test_cat <- test %>%\n    select(setdiff(colnames(test),colnames(features_test_numeric))) %>%\n    select(-c(timestamp)) %>%\n    mutate_all(as.factor)\ndmy <- dummyVars(\" ~ .\", data = features_test_cat)\ntrsf <- data.frame(predict(dmy, newdata = features_test_cat))\nfeatures_test <- cbind(features_test_numeric,features_test_cat)\nfeatures_test <- cbind(features_test_numeric,features_test_cat)\nfeatures_test <- data.matrix(features_test)","metadata":{"execution":{"iopub.status.busy":"2022-05-25T06:49:38.840846Z","iopub.execute_input":"2022-05-25T06:49:38.842293Z","iopub.status.idle":"2022-05-25T06:51:21.062555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pred <- predict(model, features_test)\npred <- as.numeric(pred > 0.5)\nsubmission = cbind(test_features$row_id,pred)\ncolnames(submission) <- c('row_id','anomaly')\nsubmission <- as.data.frame(submission) %>% mutate_all(as.integer)\nwrite.csv(submission,'submission.csv',row.names = FALSE)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Thanks for reading! Please let me know which areas need more improvement and which ones you liked :)","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}}]}