---
title: 'A Thorough Analysis of ASHRAE'
author: "Ethan Schacht"
date: "10/10/2019"
output: 
  html_document:
    toc: yes
    number_sections: yes
    code_folding: hide
    theme: simplex
    highlight: espresso
---

```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(out.width = "75%", 
               fig.align = "center")
```

<style type="text/css">

body{
  font-size: 14pt;
}

</style>

# **Introduction**

"Founded in 1894, ASHRAE  (American Society of Heating, Refrigerating and Air-Conditioning Engineers) serves to advance the arts and sciences of heating, ventilation, air conditioning refrigeration and their allied fields. ASHRAE members represent building system design and industrial process professionals around the world. With over 54,000 members serving in 132 countries, ASHRAE supports research, standards writing, publishing and continuing education - shaping tomorrow’s built environment today."

The goal of this competition is to develop accurate models from ASHRAE's 2016 data in order to better understand metered building energy usage in the following areas: chilled water, electric, hot water, and steam meters.  As a result, large scale investors and financial institutions will be more inclined to invest in this area to enable progress in building efficiencies.

For this analysis, I will provide a complete EDA of the ASHRAE data by joining the train, weather_train, and building datasets.  Feel free to leave a comment and/or an upvote if you enjoy the analysis or have any feedback.  I'm a young college student who is relatively new to data science and I'm always aiming to improve :)

![](https://www.ashrae.org/Image%20Library/Main%20Nav/About/New%20HQ/ASHRAE-mission.png)

# **Understanding the Data**

First, we will load the appropriate libraries and import / join the data and see what we are working with: 

```{r message = FALSE, warning = FALSE}

library(tidyverse)
library(lubridate)
library(extrafont)
library(scales)
library(ggridges)
library(corrplot)
library(Hmisc)
library(reshape2)

theme_set(theme_dark() + theme(text = element_text(family = "Impact")))

train <- read_csv("../input/ashrae-energy-prediction/train.csv")
building <- read_csv("../input/ashrae-energy-prediction/building_metadata.csv")
weather <- read_csv("../input/ashrae-energy-prediction/weather_train.csv")

ashrae <- train %>%
  left_join(building, by = "building_id") %>%
  left_join(weather, by = c("site_id", "timestamp"))

rm(train, building, weather)

str(ashrae)
head(ashrae, 10)

```

After joining the data, we have over 20 million rows with 16 variables.  The following is a key for all of our variables by original dataset: 

**train.csv**

* building_id - Foreign key for the building metadata.
* meter - The meter id code. Read as {0: electricity, 1: chilledwater, 2: steam, 3: hotwater}. Not every building has all meter types.
* timestamp - When the measurement was taken
* meter_reading - The target variable. Energy consumption in kWh (or equivalent). Note that this is real data with measurement error, which we expect will impose a baseline level of modeling error.

**building_meta.csv**

* site_id - Foreign key for the weather files.
* building_id - Foreign key for training.csv
* primary_use - Indicator of the primary category of activities for the building based on EnergyStar property type definitions
* square_feet - Gross floor area of the building
* year_built - Year building was opened
* floor_count - Number of floors of the building

**weather_train.csv**

Weather data from a meteorological station as close as possible to the site.

* site_id
* air_temperature - Degrees Celsius
* cloud_coverage - Portion of the sky covered in clouds, in oktas
* dew_temperature - Degrees Celsius
* precip_depth_1_hr - Millimeters
* sea_level_pressure - Millibar/hectopascals
* wind_direction - Compass direction (0-360)
* wind_speed - Meters per second

The meter_reading variable will be the primary target of the analysis for each type of meter.

## Missing Values

Now, lets count the missing values by column.

```{r message = FALSE, warning = FALSE}

colSums(is.na(ashrae))

```

Most of the NA's fall under the weather related variables.  

## Feature Engineering

### Meter

We can refer to the key in the variable descriptions to see what 0, 1, and 2 represent for the type of meter.  However, it would be easier to show these values in the dataset.

```{r message = FALSE, warning = FALSE}

ashrae2 <- ashrae %>%
  mutate(meter = if_else(meter == 0, "Electricity", if_else(meter == 1, "Chilled Water", if_else(meter == 2, "Steam", "Hot Water"))))

ashrae2 %>%
  select(meter, meter_reading) %>%
  sample_n(10)

```

### Time Components

From the timestamp variable, we will extract various time components, month of year, day of week, and hour of day.  This will be required when creating various time series.  All data takes place throughout the entire year of 2016.  It also appears that a reading is taken at each building and each meter once an hour.

```{r message = FALSE, warning = FALSE}

ashrae3 <- ashrae2 %>%
  mutate(month = month(timestamp, label = TRUE), day = wday(timestamp, label = TRUE), hour = hour(timestamp))

ashrae3 %>%
  select(month, day, hour, meter_reading) %>%
  sample_n(10)

```

# **Exploratory Data Analysis**

## Basic Distributions

Before we analyze the effects of certain independant variables on meter reading, lets visualize their distributions so we can better understand various aspects of ASHRAE.

### Site and Building ID

How many sites and buildings within those sites are involved with ASHRAE?

```{r message = FALSE, warning = FALSE}

ashrae_sample <- ashrae3

sites <- ashrae_sample %>%
  select(site_id, building_id) %>%
  arrange(site_id, building_id) %>%
  distinct() %>%
  group_by(site_id) %>%
  summarise(count = n())

sites %>%
  ggplot(aes(site_id, count)) +
  geom_bar(stat = "identity", fill = "green") +
  geom_label(aes(label = count)) +
  xlab("Site ID") +
  ylab("Number of Buildings") +
  ggtitle("Number of Buildings by ASHRAE Sight")

median(sites$count)
mean(sites$count)

```

There are 15 building sites involved with ASHRAE and an average of 90 buildings per site.

### Building Use

```{r message = FALSE, warning = FALSE}

ashrae3 %>%
  select(building_id, primary_use) %>%
  distinct() %>%
  group_by(primary_use) %>%
  summarise(count = n()) %>%
  ggplot(aes(fct_reorder(primary_use, count), count)) +
  geom_col(fill = "green") +
  geom_label(aes(label = count)) +
  coord_flip() +
  xlab(NULL) +
  ylab("Count") +
  ggtitle("Building Use Distribution") +
  theme(legend.position = "none") 

```

The majority of ASHRAE's target buildings are for the use of education.

### Square Feet

```{r message = FALSE, warning = FALSE}

sqfeet <- ashrae3 %>%
  select(building_id, square_feet) %>%
  distinct()

sqfeet %>%
  ggplot(aes(square_feet)) +
  geom_density(fill = "green") +
  geom_vline(xintercept = median(sqfeet$square_feet)) +
  geom_text(aes(250000, 6e-6, label = paste0("Median: ", median(sqfeet$square_feet), " Sq Feet")), color = "blue") +
  xlab("Square Feet") +
  ylab(NULL) +
  ggtitle("Square Footage Distribution") +
  theme(legend.position = "none", axis.text.y = element_blank())

```

The average square footage for ASHRAE buildings is over 57,000 Sq Feet.

### Year Built

We will have to remove missing values here.

```{r message = FALSE, warning = FALSE}

yearbuilt <- ashrae3 %>%
  select(building_id, year_built) %>%
  filter(!is.na(year_built)) %>%
  distinct()

yearbuilt %>%
  ggplot(aes(year_built)) +
  geom_density(fill = "green") +
  geom_vline(xintercept = median(yearbuilt$year_built)) +
  geom_text(aes(1925, 0.01, label = paste0("Median Year: ", median(yearbuilt$year_built))), color = "blue") +
  xlab("Year Built") +
  ylab(NULL) +
  ggtitle("Year Built Distribution") +
  theme(legend.position = "none", axis.text.y = element_blank())

```

The majority of ASHRAE structures were built in the late 60's and early 70's.

### Floor Count

```{r message = FALSE, warning = FALSE}

floors <- ashrae3 %>%
  select(building_id, floor_count) %>%
  filter(!is.na(floor_count)) %>%
  distinct()

floors %>%
  ggplot(aes(floor_count)) +
  geom_density(fill = "green") +
  geom_vline(xintercept = median(floors$floor_count)) +
  geom_text(aes(10, 0.15, label = paste0("Median Floor Count: ", median(floors$floor_count))), color = "blue") +
  xlab("Floor Count") +
  ylab(NULL) +
  ggtitle("Floor Count Distribution") +
  theme(legend.position = "none", axis.text.y = element_blank())

```

Most ASHRAE building are on the shorter side, with a median height of 3 stories.

### Meter

When we analyze the effects of our independant variables on meter reading, we want to do this separately for all the different meters.  Therefore, it's important to know which meter types are most common and uncommon.

```{r message = FALSE, warning = FALSE}

ashrae3 %>%
  group_by(meter) %>%
  summarise(count = n()) %>%
  ggplot(aes(fct_reorder(meter, count), count)) +
  geom_bar(stat = "identity", fill = "green") +
  geom_label(aes(label = count)) +
  xlab("Meter Type") +
  ylab("Count")
  ggtitle("Number of Readings for Each Meter Type") +
  theme(legend.position = "none")

```

Electricity by far receives the most readings.

### Meter Reading (Dependant Variable)

```{r message = FALSE, warning = FALSE}

ashrae3 %>%
  sample_n(500000) %>%
  ggplot(aes(meter_reading)) +
  geom_histogram(fill = "green") +
  xlab("Meter Reading") +
  ylab("Count") +
  ggtitle("Distribution of Meter Readings")

```

Nearly all of the meter readings are 0 or somewhere close, yet there are some that read as high as the 10's of millions.  We will put this on a logarithmic scale to better visualize the distribution.

```{r message = FALSE, warning = FALSE}

ashrae_sample %>%
  ggplot(aes(meter_reading)) +
  geom_histogram(fill = "green") +
  xlab("Meter Reading") +
  ylab("Count") +
  ggtitle("Logarithmic Distribution of Meter Readings") +
  scale_x_log10() +
  geom_vline(xintercept = median(ashrae_sample$meter_reading))
  
median(ashrae_sample$meter_reading)
```

Meter readings normalize logarithmically. The median meter reading is 79 kWh.  **We will use a logarithmic scale throughout this entire analysis.**

## Effects on Meter Reading  

Now it's time to analyze effects of our independant variables on our variable of interest.  We will always group by meter type for a sufficiently detailed analysis.

### Building Effects

First, we will measure the effect of our building-related variables on the four different meter readings.

#### Which Sites Use the Most Energy?

```{r message = FALSE, warning = FALSE}

sites2 <- ashrae_sample %>%
  select(site_id, meter_reading, meter) %>%
  arrange(site_id) %>%
  mutate(site_id = paste(site_id)) 
  
sites2 %>%
  mutate(side_id = factor(site_id, levels = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15"))) %>%
  ggplot(aes(x = site_id, y = meter_reading)) +
  geom_boxplot(fill = "green") +
  scale_y_log10() +
  facet_wrap(~meter) +
  xlab("Site ID") +
  ylab("Meter Reading (kWh)") +
  ggtitle("Average Site Meter Readings for Each Meter Type", subtitle = "Missing plots represent lack of meter type")

```

Site 7 appears to use the most energy on average.  Note that this plot is also on a logarithmic scale.

#### Do Larger Buildings Use More Energy?

I would expect larger buildings to take up more energy.  Let's take a look...

```{r message = FALSE, warning = FALSE}

sqfeet2 <- ashrae_sample %>%
  select(meter, meter_reading, square_feet, building_id) %>%
  group_by(building_id, square_feet, meter) %>%
  summarise(avg_meter = median(meter_reading))

sqfeet2 %>%
  ggplot(aes(square_feet, avg_meter)) +
  geom_point() +
  geom_smooth(color = "green", method = "lm") +
  facet_wrap(~meter) +
  xlab("Building Square Footage") +
  ylab("Average Meter Reading (kWh)") +
  ggtitle("Building Square Footage vs. Average Meter Reading for Each Meter Type") 
  
feetcor1 <- sqfeet2 %>%
  filter(meter == "Electricity")
paste0("Electricity Meter Correlation Value: ", round(cor(feetcor1$square_feet, feetcor1$avg_meter), 2))

feetcor2 <- sqfeet2 %>%
  filter(meter == "Chilled Water")
paste0("Chilled Water Meter Correlation Value: ", round(cor(feetcor2$square_feet, feetcor2$avg_meter), 2))

feetcor3 <- sqfeet2 %>%
  filter(meter == "Hot Water")
paste0("Hot Water Meter Correlation Value: ", round(cor(feetcor3$square_feet, feetcor3$avg_meter), 2))

feetcor4 <- sqfeet2 %>%
  filter(meter == "Steam")
paste0("Steam Meter Correlation Value: ", round(cor(feetcor4$square_feet, feetcor4$avg_meter), 2))
```

Square Footage has the strongest positive correlation with meter reading for Electricity Meters, though there is significant correlation for all meter types.  Steam meters have the capability to show signicifant higher readings than the other meter types.

#### Which Building Uses Take up the Most Energy?

```{r message = FALSE, warning = FALSE}

build_use <- ashrae_sample %>%
  select(meter, meter_reading, primary_use, building_id) %>%
  group_by(building_id, primary_use, meter) %>%
  summarise(avg_meter = median(meter_reading))

build_use %>%
  ggplot(aes(primary_use, log(avg_meter) + 1)) +
  geom_bar(stat = "identity", fill = "green") +
  facet_wrap(~meter) +
  xlab("Building Use") +
  ylab("Average Meter Reading (kWh)") +
  ggtitle("Building Usage vs. Meter Reading for Each Meter Type") +
  coord_flip()

```

Buildings used for education take up the most energy by far on average.  However, is this correlation or causation? It's possible that educational buildings are larger...

```{r message = FALSE, warning = FALSE}

bigbuild <- ashrae_sample %>%
  select(building_id, primary_use, square_feet) %>%
  distinct() %>%
  group_by(primary_use) %>%
  summarise(avg_sqfeet = mean(square_feet), n = n())

bigbuild %>%
  ggplot(aes(fct_reorder(primary_use, avg_sqfeet), avg_sqfeet, fill = n)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  ylab("Average Square Feet") +
  xlab(NULL) +
  ggtitle("Average Square Footage by Building Usage") +
  labs(fill = "Building Count")

```

Education is ranked only 4th on the list, implying that square footage has little to do with the relationship between primary use and average meter reading.  However, it is noteable that education is by far the most common building usage that ASHRAE services.

### Weather Effects {.tabset .tabset-fade .tabset-pills}

Now, for the effects of the weather on average meter reading.  As there are many different variables and they are all continuous, a correlation plot would be perfect for each meter type.  Note that omitting all the rows with missing values eliminates over half of the data, though there still plenty of statistical significance here.

#### Overall

```{r message = FALSE, warning = FALSE}

weather <- ashrae_sample %>%
  select(10:16, meter, meter_reading) %>%
  na.omit() 

overall <- weather %>%
  select(-meter)

corrmat1 <- rcorr(as.matrix(overall))

corrplot(corrmat1$r, type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)

corrmat1$r

```

Overall, there exists almost no correlation between any of the weather related variables and meter reading.  However, some of the variables seem to affect hot water more than they affect the other meters.  The most significant correlation is the negative relationship between air temperature and meter reading on the hot water meter.  This could be due to lower temperatures causing more difficulty in extracting hot water from the pipes.

#### Chilled Water

```{r message = FALSE, warning = FALSE}

chilled <- weather %>%
  filter(meter == "Chilled Water") %>%
  select(-meter)

corrmat2 <- rcorr(as.matrix(chilled))

corrplot(corrmat2$r, type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)

corrmat2$r

```

Overall, there exists almost no correlation between any of the weather related variables and meter reading.  However, some of the variables seem to affect hot water more than they affect the other meters.  The most significant correlation is the negative relationship between air temperature and meter reading on the hot water meter.  This could be due to lower temperatures causing more difficulty in extracting hot water from the pipes.

#### Hot Water

```{r message = FALSE, warning = FALSE}

hot <- weather %>%
  filter(meter == "Hot Water") %>%
  select(-meter)

corrmat3 <- rcorr(as.matrix(hot))

corrplot(corrmat3$r, type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)

corrmat3$r

```

Overall, there exists almost no correlation between any of the weather related variables and meter reading.  However, some of the variables seem to affect hot water more than they affect the other meters.  The most significant correlation is the negative relationship between air temperature and meter reading on the hot water meter.  This could be due to lower temperatures causing more difficulty in extracting hot water from the pipes.

#### Electricity

```{r message = FALSE, warning = FALSE}

electric <- weather %>%
  filter(meter == "Electricity") %>%
  select(-meter)

corrmat4 <- rcorr(as.matrix(electric))

corrplot(corrmat4$r, type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)

corrmat4$r

```

Overall, there exists almost no correlation between any of the weather related variables and meter reading.  However, some of the variables seem to affect hot water more than they affect the other meters.  The most significant correlation is the negative relationship between air temperature and meter reading on the hot water meter.  This could be due to lower temperatures causing more difficulty in extracting hot water from the pipes.

#### Steam

```{r message = FALSE, warning = FALSE}

steam <- weather %>%
  filter(meter == "Steam") %>%
  select(-meter)

corrmat5 <- rcorr(as.matrix(steam))

corrplot(corrmat5$r, type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)

corrmat5$r

```

Overall, there exists almost no correlation between any of the weather related variables and meter reading.  However, some of the variables seem to affect hot water more than they affect the other meters.  The most significant correlation is the negative relationship between air temperature and meter reading on the hot water meter.  This could be due to lower temperatures causing more difficulty in extracting hot water from the pipes.

### Time Series

In general, how did different meter readings behave over periods of time in 2016? Let's take a look...

#### Meter Readings Throughout the Year

First, how did meter readings differ depending on the month?

```{r message = FALSE, warning = FALSE}

months <- ashrae_sample %>%
  select(month, meter, meter_reading) %>%
  group_by(month, meter) %>%
  summarise(avg_meter = median(meter_reading)) %>%
  arrange(meter)

months %>%
  ggplot(aes(month, avg_meter, group = 1, color = meter)) +
  geom_point(stat='summary', fun.y=sum) +
  stat_summary(fun.y=sum, geom="line") +
  xlab("Month") +
  ylab("Average Meter Reading (kWh)") +
  ggtitle("Average Meter Readings Throughout the Year") +
  facet_wrap(~meter) +
  theme(legend.position = "none")
  

```

All of these time series make sense as people will want to use more cold forms of energy in the hotter American months of June through August and more hot forms of energy in the colder Amercian months of November through February.  The magnitude to which output of steam energy changes throughout the year, however, is greater than the others.

#### Meter Readings Throughout the Week

Now, for the seven weekdays.

```{r message = FALSE, warning = FALSE}

months <- ashrae_sample %>%
  select(day, meter, meter_reading) %>%
  group_by(day, meter) %>%
  summarise(avg_meter = median(meter_reading)) %>%
  arrange(meter)

months %>%
  ggplot(aes(day, avg_meter, group = 1, color = meter)) +
  geom_point(stat='summary', fun.y=sum) +
  stat_summary(fun.y=sum, geom="line") +
  xlab("Weekday") +
  ylab("Average Meter Reading (kWh)") +
  ggtitle("Average Meter Readings Throughout the Week") +
  facet_wrap(~meter) +
  theme(legend.position = "none")
  

```

In general, the least amount of energy is used on the weekends (Saturday and Sunday), which is due to the fact that Americans usually don't work on these days.  The magnitude of energy output change, however, is a bit smaller than I expected.

#### Meter Readings Throughout the Day

Lastly, throughout the 24 hours of the day.

```{r message = FALSE, warning = FALSE}

months <- ashrae_sample %>%
  select(hour, meter, meter_reading) %>%
  group_by(hour, meter) %>%
  summarise(avg_meter = median(meter_reading)) %>%
  arrange(meter)

months %>%
  ggplot(aes(hour, avg_meter, group = 1, color = meter)) +
  geom_point(stat='summary', fun.y=sum) +
  stat_summary(fun.y=sum, geom="line") +
  xlab("Hour") +
  ylab("Average Meter Reading (kWh)") +
  ggtitle("Average Meter Readings Throughout the Day") +
  scale_x_continuous(breaks = c(0,9,17,23)) +
  facet_wrap(~meter) +
  theme(legend.position = "none")
  

```

The 0th hour is midnight and the 12th hour is noon.  The average American work day takes place from 9 AM to 5 PM, or hour 9 through hour 17 on this chart.  This is unsruprisingly the general timeframe of greatest energy output for all of the meters, though steam output interestingly peaks before the average work hours.

# **Conclusions** 

From this analysis, we can draw the following conclusions on the ASHRAE dataset in regards to effects on various meter readings:

* Site 7 uses the greatest amount of energy, especially when it comes to electricity.
* Square Footage of buildings positively effects electricity meter readings more than any other meter type.
* Educational Buildings exert the most energy in all meter categories by far in comparison to the other building uses, and ASHRAE servives these types of buildings the most.
* The weather has little to no effect on meter readings, though the strongest of these effects is the negative correlation between air temperature and hot water meter reading.
* Electricity Meter Readings don't change much throughout the year.  Hot water and steam meters peak in the winter months and trough in the summer months, whereas cold water meters do the opposite.
* The weekends (Saturday and Sunday) generally experience the lowest energy output in comparison to all the weekdays.
* Throughout the day, energy output tends to be at its highest levels throughout the common work hours of 9 AM - 5 PM, though steam output peaks before the work hours.

** Thank you for reading my analysis and remember to leave an upvote if you liked it! :) I am a young college student and consider myself very new to data science so any feedback would be appriciated, as well, to help me improve.  I will likely update this kernel with more analysis as there is a lot more information that can be extracted. **