As cities keep growing worldwide, so does the demand for key resources such as energy (electricity and gas) and water that residents consume . Meeting the demand for these resources can be challenging and requires an understanding of the ir consumptions patterns . In this work, we apply XGBoost (Extreme Gradient Boosting) to predict and analyze water and energy consumption in large -scale buildings in New York City. For this, the New York City’s local law 84 extensive dataset was merged with the Primary Land Use Tax Lot Output (PLUTO) dataset as well as with other socio -economic databases. Specifically, we developed three models: electricity, gas, and water consumption. Seven major lessons were learnt in terms of interrelationships between electricity, gas, and water consumption. In particular, water and gas consumption are highly interrelated with one another (often because gas is used for water heating). Furthermore, electricity consumption is affected by building type, and electricity and water consumption are particularly interrelated in nonresidential buildings. Overall, the knowledge gained from the models and from the SHAP analysis can help planners, engineers, and policymakers develop mor e effective strategies and help them manage the demand for energy and water in large -scale buildings. Keywords : urban sustainability ; water consumption; electricity consumption; gas consumption; machine learning _________________ *Corresponding author. Email Addresses: amovah2@uic.edu and derrible@uic.edu 2 1. Introduction As the world population keeps growing, so does the demand for key energy and resources. Moreover, most of this growth is happening in cities that already account for over 67% of energy use worldwide (Hobley, 2019) . As a result, providing the necessary reso urces and energy for a growing population has become one of the most significant challenges for engineers, planners, and policymakers, while needing to acc ount for local specificities (Mohareb et al., 2016) . Specifically, building energy consumption represents a significant portion of total energy consumption and carbon emissions worldwide (Dalia, 2014) . Buildings in the United States (US), for instance, contribute to 69% of electricity use, 39% of energy use, and more than 29% of total U.S. CO 2e emissions (Becerik -gerber et al. 2013; USEPA 2017) . Similarly, despite being challenging (Ahmad and Derrible, 2015) , the importance of studying water consumption trends is essential, especially considering the fact that 36% of the world’s population experiences water scarcity (Von Grebmer et al., 2012) . In this context, understanding building energy and water consumption is fundamental to help prepare cities meet future demand in energy and resources. In this respect, t he literature on building energy and water consumption is relatively rich. Nonetheless, existing studies tend to study energy and water separately , without looking at the interrelationships between the two . Building upon previous works , the main contribution of this article is to gain a fundamental understanding of the interrelationships between electricity, gas, and water consumption patterns . In the literature , an important knowledge gap exists to better understand how interconnected and interdependent the demand for and supply of urban infrastructure systems (Derrible, 201 7; 2018) . This article further leverages advances in Machine Learning (ML) to uncover meaningful patterns between the consumption of electricity, gas, and water in buildings. While ML techniques are generally less interpretable than traditional approaches that are b ased on strong predetermined assumptions (e.g., Kashani et al., 2019; Parsa et al., 2019a ), to partly alleviate this issue, model- agnostic methods can be adopted to help interpret what ML models have essentially “learned.” Various model -agnostic measures exist to interpret ML models such as feature importance (FI), partial dependence (PD), and individual conditional expectation (ICE) (Goldstein et al., 2013; Hastie et al., 2009) . This article adopts the SHAP (SHapley Additive exPlanation) package (Lundberg and Lee, 2017) , built in Python and previously used to study building energy consumption ( Padopoulos and Kontokosta 2019) . More specifically, this article focuses on investigating the interrelationships between electricity, gas, and water use patterns in large -scale buildings using ML . Toward this goal, we use XG Boost to model resource consumptions —that is, electricity, gas, and water —and SHAP to interpret the learned ML models. The specific objectives of this article are to : • Describe and analyze the performance of XG Boost to model electricity, gas, and water consumption, • Investigate the role of specific variables for each model separately using SHAP , • Study the interrelationships between electricity, gas, and water consumption in large -scale buildings using SHAP . In the next section, we first briefly review the ML literature on building energy and water consumption. Then, we describe the datasets used in this study and we explain the methodological approach adopted. In each model, we predict either electricity, gas, or water based on t he other 3 two variables; for example, we use both water and gas consumption to predict electricity consumption. Furthermore, we include variables related to land use and socioeconomic variables that are commonly included in energy and resource consumption m odeling. Subsequently, we use SHAP to interpret the results, identify the impact of individual variables, and study the interrelationships between electricity, gas, and water consumption. 2. Literature Review The literature on building energy and water consumption is relatively abundant . The body of work on building energy and resource use is significantly aided by energy disclosure laws that enable the provision of open- source datasets (Athanassiadis, 2020) such as the New York local law 84 dataset that is used in this study and that notably contains information on water, electricity, and gas consumption in buildings with a total floor area greater than 50,000 square feet. Kontokosta and Tull (2017) applied machine learning to predict buildi ng energy use in New York City. Kontokosta and Jain (2015) used weighted robust regression and geographically weighted regression models to investigate and understand the determinants of wat er consumption in urban multi- family housing in New York City. Hsu (2015) demonstrated the appl ication of penalized regression on commercial office and multifamily buildings energy consumption data set in New York City. Hsu (2015b) compared integrated clustering methods in terms of predicting building energy consumption and found that clusterwise regression gives more accurate predictions while K-means gives more stable clusters. Ma and Cheng (2016) analyzed the impact of 171 different kinds of features on regional energy use intensity (EUI) of residential buildings in New York City and identified the top 20 influential features. Meng et la. ( 2017) evaluated the impact of benchmarking policies on saving energy in the buildings in New York City. Using support vector machine, Jain et al. ( 2014) developed a model to forecast building energy consumption and applied the model to data from multi -famil y residential buildings in New York City. Finally, leveraging extreme gradient boosting and shapely additive explanation (used in this study as well), Papadopoulos and Kontokosta (2019) developed an interpretable and robust method to grade energy performance of residential buildings in New York City. In general, m any different factors are found to affect building energy consumption, including building physical properties and occupants behavior (Ashouri et al., 2019; Li et al., 2019) , which makes energy consumption prediction challenging (Kwok and Lee, 2011) . In particular, socioeconomic attributes can greatly affect energy and water consumption behavior. Several studies show that household income, household members age, education, and the percentage of female occupants affect building energy consumption (Pérez- Lombard et al., 2008; Tso and Yau, 2007) . Other studies show that similar variables impact water consumption (Christenson et al., 2006; Kenney et al., 2008; Kontokosta and Jain, 2015; Wentz and Gober, 2007) . Besides socioeconomic att ributes, land use has also shown to have a significant impact on energy and water consumption (Pérez- Lombard et al., 2008) . For example, energy and water consumptions patterns differ significantly depending on their type, from official to residential and commercial buildings, and (Pérez- Lombard et al., 2008) found that the average energy use intensity in offices is twice that of residential buildings. In this work, we look both at the effects of socioeconomic and land use factors on electricity, gas, and wa ter consumption. From a modeling perspective, various types of approaches have been adopted to model energy and water consumption. Beyond some of the approaches already mentioned, Suga nthi and Samuel (2012) reviewed several modeling algorithms used to mode l energy consumption, 4 including regression models, econometric models, and artificial systems. Donkor et al. (2014) reviewed the performance of econometric models, Artificial Neural Networks (ANN), and simulation or scenario -based modeling on water consumption. Recently, machine learning (ML) approaches ha ve gained significant interest (Bourdeau et al., 2019; Seyedzadeh et al., 2019; Yu et al., 2016; Froemelt et al., 2019) and tend to perform better than traditional statistical regression modeling technique s thanks to their machine -based algorithmic computation. In particular, Decision Tree (DT) ( Badhrudeen et al., 2020), ANN, ensemble methods (e.g., bagging and boosting methods), Support Vector Machine (SVM) (Parsa et al. , 2019b) , and regularized linear regression (e.g., ridge and lass o regression) tend to be the most popular ML approaches applied in the previous studies (Wei et al., 2018) . For instance, Ahmad et al. (2014) forecasted building electricity consumption by using ANN and SVM , and Li et al. (2015) predicted building electricity consumption using ANNs and principle component analysis (PCA). Deng et al. (2018) compared several ML techniques to estimate commercial building energy consumption. Furthermore, various studies (including some already mentioned) compare the performance of different methods. Papadopoulos et al. (2018a) evaluated the performance of random forest (RF), extremely randomized tree (RT), and gradient boosted regression tree (GBRT) to predict heating and cooling loads of buildings and found out GBRT is outperforming RF and RT. Robinson et al. (2017) also compared several techniques to estimate commercial bui lding energy consumption and found that Gradient Boosting Machine (GBM) regression performed well. Lee and Derrible (2019) reviewed the performance 12 different techniques (four traditional and eight machine learning techniques) to predict household water consumption under two different scenarios (one with only a selected few variables generally available to utilities and one with an extensive number of variables) and found that gradient boosting performed best in both scenarios. Finally, the application of extreme gradient boosting and SHAP in modeling energy consumption in buildings has been utilized by Papadopoulos and Kontokosta (2019) as mentioned previously. 3. Data and Methodologies 3.1 Data Description This article mainly use s three datasets: (1) the New York City’s Local Law 84 Water and Energy benchmarking dataset 2017 (NYC’s local law 84), (2) the Primary Land Use Tax Lot Output (PLUTO) dataset , and (3) the American Community Survey (ACS) from the US Census Burea u. NYC’s local law 84 applies to buildings with a gross floor area larger than 50,000 square feet and it mainly contains information about the consumption of water, electricity, and gas in buildings , as well as several other building characteristics such as building age, gross floor area, and building usage type. PLUTO dataset files contain three types of data: tax lot characteristics, building characteristics, and geographic/political/administrative districts. Notably, the PLUTO dataset provides detailed land- use information at the tax lot level in New York City. In addition, the ACS dataset is used to incorporate socioeconomic information available at the census tract level including gender proportion, age distribution, education level , and median income . 5 The original dataset consist s of a total of 11,746 buildings . Nonetheless, it contain s many erroneous entries since the data are reported by individual building owners. In terms of cleaning, first, all rows that were missing values for either electricity, gas, or water consumption were removed since these values are essential for this work . Second, t he dataset includes 937 building s that use steam from the New York City large district steam system . All were removed since energy consumption from steam can be significantly large (e.g., for building heating) but 937 is too low to consider steam as a fourth dependent variable and most cities do not possess a district steam system . Finally , a traditional cleaning process was applied to remove all errors, outliers, duplicated data points , and other buildings with missing data . The final cleaned dataset included 6,527 buildings and fifteen variables. A list of the variables collected and a summary of their statistics are shown in Table 1. Table 1 Descriptive statistics of the merged dataset Variable Definitions Mean Std Max Min Electricity Electricity consumption per square meter per year [𝑘𝑘𝑘𝑘ℎ/𝑚𝑚2/𝑦𝑦]. 49 66.95 769 1.5 Gas Gas consumption per square meter per year [𝑘𝑘𝑘𝑘ℎ/𝑚𝑚2/𝑦𝑦]. 298 307.75 2736.56 0 Water Water consumption per square meter per year [𝐿𝐿/𝑚𝑚2/𝑦𝑦] 4239 3070.68 17670.23 87.02 Residential, Office, Commercial Dummy variables to identify the type of the major part of a building. ___ ___ ___ ___ Age Age of a building [ y]. 71 31.48 191 2 More 60 The percentage of people older than 60 years old in a census tract. 18.64 8.09 58.7 1.7 Under 9 The percentage of people younger than 9 years old in a census tract. 10.91 4.84 41.7 0.8 Female Percentage The percentage of female in a census tract. 52.7 4 72 31 Median Income Median income of people in a census tract. 71253 40715 233529 9740 Building Density The total gross area of buildings in a census tract divided by census tract area [𝑚𝑚2/𝑚𝑚2]. 1.65 1.35 10.47 0.05 Residential Area The total gross area of residential buildings in a census tract divided by census tract area [𝑚𝑚2/𝑚𝑚2]. 0.98 0.58 4.46 0.004 Office Area The total gross area of offices in a census tract divide by census tract area [𝑚𝑚2/𝑚𝑚2]. 0.3 0.74 6.27 0 Factory Area The total gross area of factories in a census tract divide by census tract area [𝑚𝑚2/𝑚𝑚2]. 0.014 0.085 1.59 0 No High School Diploma The percentage of people who did not get high school degree. 13.6 12.06 57.9 0 College Education The percentage of people who got bachelor degree or more in a census tract. 53 26.92 93.9 5 3.2 Methodologies To model and understand building energy consumption, this article mainly uses extreme gradient boosting (XG Boost) . Boosting—as an ensemble method —is developed to enhance the modeling performances of decision tree (DT) , executed through a repetitive machine -based algorithm and effectively controlling the bias-variance trade -off (Friedman et all ., 2001) . Specifically, a DT has 6 a structure similar to a tree with a root node (topmost node), internal node s, and leaf node s (end node s). DT algorithms use relatively simple rules to start from a root node and split up each node into two or more branches , going through internal nodes, and ending with the leaves. As “weak learners,” DT tend to have a low performance and be unstable such that a small change in the data can make a large change in the tree structure. To overcome the weakness of DT, XG Boost models uses an ensemble of weak consecutive decision trees wh ere each tree learns from the previous tree and tries to improve the model. XG Boost is detailed in the next section. 3.2.1 Extreme Gradient Boosting - XGBoost XGBoost is an advanced version of gradient boosting machine. In this section, we briefly explain the main concepts and equations behind it, but the reader is referred to Chen and Guestrin ( 2016) for more details. In a dataset with n samples, each independent variable 𝑥𝑥 𝑖𝑖 can have m features, 𝑥𝑥𝑖𝑖∈ℝ𝑚𝑚. Corresponding to each of these variables, we define a dependent variable 𝑦𝑦𝑖𝑖, 𝑦𝑦𝑖𝑖∈ℝ. By using independent variables 𝑥𝑥𝑖𝑖, a tree ensemble model f uses K additive functions to predict the dependent variable such that : ŷ𝑖𝑖= ∑𝑓𝑓𝑘𝑘(𝑥𝑥𝑖𝑖)𝐾𝐾 𝑘𝑘=1 ,𝑓𝑓𝑘𝑘∈𝐹𝐹 (1) where ŷ 𝑖𝑖 is the predicted value of 𝑦𝑦𝑖𝑖. Moreover, here, 𝐹𝐹 is the space of regression trees and 𝑓𝑓𝑘𝑘 corresponds to an independent tree structure with leaf scores. Our goal is to minimize the following function: ℒ(𝜙𝜙)=∑𝑙𝑙(y�𝑖𝑖,𝑖𝑖 𝑦𝑦𝑖𝑖)+∑Ω(𝑓𝑓𝑘𝑘) 𝑘𝑘 (2) where Ω(𝑓𝑓)=𝛾𝛾Τ+1 2𝜆𝜆∥𝜔𝜔∥2 (3) Here 𝑙𝑙 is a loss function that measur es the difference between the predicted and the target values , Ω is a term that penaliz es the complexity of the model , T is the number of leaves , and 𝜔𝜔𝑖𝑖 is the score on i -th leaf. For a fixed structure 𝑞𝑞 (𝑥𝑥), the optimal value for 𝜔𝜔𝑗𝑗∗ and its corresponding values is : 𝜔𝜔𝑗𝑗∗=−∑𝜕𝜕y�𝑡𝑡−1𝑙𝑙�𝑦𝑦𝑖𝑖,𝑦𝑦�𝑖𝑖𝑡𝑡−1� 𝑖𝑖∈𝐼𝐼𝑗𝑗 ∑𝜕𝜕2 y�𝑡𝑡−1𝑙𝑙�𝑦𝑦𝑖𝑖,𝑦𝑦�𝑖𝑖𝑡𝑡−1�+𝜆𝜆 𝑖𝑖∈𝐼𝐼𝑗𝑗 (4) ℒ̃𝑡𝑡(𝑞𝑞)= −12∑(∑𝜕𝜕y�𝑡𝑡−1𝑙𝑙�𝑦𝑦𝑖𝑖,𝑦𝑦�𝑖𝑖𝑡𝑡−1� 𝑖𝑖∈𝐼𝐼𝑗𝑗)2 ∑𝜕𝜕2 y�𝑡𝑡−1𝑙𝑙�𝑦𝑦𝑖𝑖,𝑦𝑦�𝑖𝑖𝑡𝑡−1�+𝜆𝜆 𝑖𝑖∈𝐼𝐼𝑗𝑗𝑇𝑇 𝑗𝑗=1 +𝛾𝛾Τ (5) Here 𝑦𝑦�𝑖𝑖𝑡𝑡is the predicted value of i-th instance at the t-th iteration and define 𝐼𝐼𝑗𝑗= {𝑖𝑖|𝑞𝑞(𝑥𝑥𝑖𝑖)=𝑗𝑗} where q represents the structure of a tree, thus I j is the instance set of j . In practice , it is impossible to calculate this value for all possible tree structures q to determine which split should be chosen. Instead , the following formula is used to evaluate the split candidates: 7 ℒ𝑠𝑠𝑠𝑠𝑙𝑙𝑖𝑖𝑡𝑡= 1 2�(∑𝜕𝜕y�𝑡𝑡−1𝑙𝑙�𝑦𝑦𝑖𝑖,𝑦𝑦�𝑖𝑖𝑡𝑡−1� 𝑖𝑖∈𝐼𝐼𝐿𝐿)2 ∑𝜕𝜕2y�𝑡𝑡−1𝑙𝑙�𝑦𝑦𝑖𝑖,𝑦𝑦�𝑖𝑖𝑡𝑡−1�+𝜆𝜆 𝑖𝑖∈𝐼𝐼𝐿𝐿+(∑𝜕𝜕y�𝑡𝑡−1𝑙𝑙�𝑦𝑦𝑖𝑖,𝑦𝑦�𝑖𝑖𝑡𝑡−1�) 𝑖𝑖∈𝐼𝐼𝑅𝑅2 ∑𝜕𝜕2y�𝑡𝑡−1𝑙𝑙�𝑦𝑦𝑖𝑖,𝑦𝑦�𝑖𝑖𝑡𝑡−1�+𝜆𝜆 𝑖𝑖∈𝐼𝐼𝑅𝑅−(∑𝜕𝜕y�𝑡𝑡−1𝑙𝑙�𝑦𝑦𝑖𝑖,𝑦𝑦�𝑖𝑖𝑡𝑡−1� 𝑖𝑖∈𝐼𝐼 )2 ∑𝜕𝜕2y�𝑡𝑡−1𝑙𝑙�𝑦𝑦𝑖𝑖,𝑦𝑦�𝑖𝑖𝑡𝑡−1�+𝜆𝜆 𝑖𝑖∈𝐼𝐼�−𝛾𝛾 (6) where 𝐼𝐼𝐿𝐿 and 𝐼𝐼𝑅𝑅 are the instances sets of left and right nodes after the split and 𝐼𝐼=𝐼𝐼𝐿𝐿∪𝐼𝐼𝑅𝑅 (Chen and Guestrin, 2016) . In the context of building electricity, gas, and water consumption, 𝑥𝑥𝑖𝑖 represent s the independent variables of a building (see Table 1) and 𝑦𝑦 𝑖𝑖 represents one of the three consumption values. For example for electricity consumption , 𝑦𝑦𝑖𝑖 is electricity consumption and 𝑥𝑥𝑖𝑖 includes all other m variables that are used to predict electricity consumption . 3.2.2 Model Parameters and Selection XGBoost model s can be defined by selecting specific values for several parameters, some of which are brief ly explained in this section . The parameters should be selected to maximize the model performance on a given dataset and in such a way to prevent overfitting and unnecessary complexity . Overfitting happens when a model starts to learn the noises and random fluctuations present in the data and sees them as meaningful facts or concepts. XGBoost uses multiple parameters, which makes parameter tuning relatively difficult. Such parameters include number of boosting iterations, maximum depth, random subsampling, learning rate , lambda, and alpha . The n umber of boosting iterations in XG boost refers to the number of consecutive trees that are fitted . The m aximum depth (max_depth) of the tree indicates the maximum number of splits ; too high a maximum depth can sometimes cause overfitting . Randomly subsampling corresponds to a specific ratio of the training dataset in each iteration prior to growing trees ; the same technique of subsam pling can be applied to columns and in each iteration a ratio of columns is used to grow the trees. Learning rate (eta) is used to shrink the weights at each step and make the model more robust ; essentially , the learning rate reduces the impact of each individual tree and let fut ure trees improve the model (Chen and Guestrin, 2016) . The parameters “lambda” and “alpha” respectively ar e L2 and L1 regularization term on weights , and by increasing them the model becomes more conservative. The following values represent the final parameters we u sed for the models after the cross- validati on process (explained in the next section) : number of boosting iterations : 500; maximum tree depth: 8; subsample ratio of training dataset in each split : 0.8; subsample ratio of columns in each split: 0.5; lambda : 1.5; a lpha: 0.2; learning rate : 0.015. 3.2.3 Model Execution and Validation In this work , the hold- out method is used to split the data into two groups: 80% for training and 20% for testing. In addition, we evaluate the model through 10-fold cross validation (CV). Specifically, we randomly split the training dataset only into ten folds. Then, in each step, we set a unique fold as our test data set and use the other nine folds to train a model and assess its performance on the test fold. We repeat the process so that all ten folds are using for testing once . The performance of each of the ten models during the CV process is calculated . Specifically, we look for a stable performance across all ten models to ensure that the technique used (i.e., XG 8 Boost) is meaningful . The combined average performance measures are calculated as well. If the CV process is satisfactory, the model is retrained on the full 80% of the data used as the training set and tested on the test set. In terms of performance measure, in this work, we use goodness of fit R2, defined as: 𝑅𝑅2=1−∑(𝑦𝑦𝑖𝑖−𝑦𝑦�𝑖𝑖)2𝑖𝑖 ∑(𝑦𝑦𝑖𝑖−𝑦𝑦�)2𝑖𝑖 (7) where 𝑦𝑦𝑖𝑖 is the actual value of a data point , 𝑦𝑦�𝑖𝑖 is the predicted value, and 𝑦𝑦� is the mean of all the actual values . 3.2.4 Model Interpretation In this article , we further use SHAP (Shap ley Additive exPlanations) values to interpret the outputs. SHAP is based on game theory (Strumbelj and Kononenko, 2014) and local explanations (Ribeiro et al., 2016) , and it unifies ideas from them in order to represent a locally accurate and consistent feature attribution (Lundberg and Lee, 2017) . Let us assume a game 𝑣𝑣. In the game, there is a group N (with n members ) with a specific output and the value of the output is 𝑣𝑣(𝑁𝑁). The idea behind shaple y values is the value of the production 𝑣𝑣(𝑁𝑁) should be divided among the n members based on their marginal contributions (Shapley, 1953) . In XGBoost , n represents the features (i.e., the independent variables) and the value of the production is the target value (i.e., the dependent variable) . Therefore, we want to know how much each feature contributes (𝜙𝜙𝑖𝑖: attribution of feature i) in the final prediction of our model . There are four axioms t hat help us to fairly identify their contributions (Shapley, 1953) : • Efficiency axiom: the sum of the feature contributions is equal to the final prediction , such that ∑𝜙𝜙𝑖𝑖=𝑣𝑣(𝑁𝑁) 𝑖𝑖𝑖𝑖𝑖𝑖 • Symmetry axiom: if two features i and j have the same marginal contribution when we add them to any feature combinations for prediction , they should receive the same contribution from the final predictio n, such that if 𝑣𝑣(𝑆𝑆∪{𝑖𝑖})=𝑣𝑣(𝑆𝑆∪{𝑗𝑗}) where 𝑆𝑆⊂𝑁𝑁 𝑎𝑎𝑎𝑎𝑎𝑎 𝑖𝑖,𝑗𝑗 ∉𝑆𝑆, then 𝜙𝜙𝑖𝑖(𝑣𝑣)=𝜙𝜙𝑗𝑗(𝑣𝑣) • Dummy axiom: if by adding a feature to any combinations of features, the prediction will not be changed, the feature will not receive any contribution from the final prediction , such that if 𝑣𝑣 (𝑆𝑆∪{𝑖𝑖})=𝑣𝑣(𝑆𝑆) where 𝑆𝑆 ⊂𝑁𝑁 and 𝑖𝑖 ∉𝑆𝑆, then 𝜙𝜙 𝑖𝑖(𝑣𝑣)=0 • Linearity axiom: if the model for prediction divided into two parts, then the total contribution of a feature is equal to the sum of the feature contribution in each of them , such that if 𝑣𝑣 (𝑆𝑆)=𝑣𝑣1(𝑆𝑆)+𝑣𝑣2(𝑆𝑆) for all S, then 𝜙𝜙 (𝑣𝑣)=𝜙𝜙(𝑣𝑣1)+𝜙𝜙(𝑣𝑣2) Based on these axioms, shapl ey values are determined by the following formula: 𝜙𝜙𝑖𝑖=∑|𝑆𝑆|!(𝑛𝑛−|𝑆𝑆|−1)! 𝑛𝑛!𝑆𝑆⊆𝑖𝑖{𝑖𝑖} [𝑣𝑣(𝑆𝑆∪{𝑖𝑖})−𝑣𝑣(𝑆𝑆)] (8) In the additive feature attribution method, there is a linear function of binary variables 𝑔𝑔 defined as: 𝑔𝑔�𝑧𝑧′�=𝜙𝜙𝑖𝑖+∑𝜙𝜙𝑖𝑖𝑧𝑧𝑖𝑖′ 𝑀𝑀 𝑖𝑖=1 (9) 9 where 𝑧𝑧′=1 when a feature ob served, otherwise it is 0, and M is the number of input features (Lundberg and Lee, 2017) . 4. Results Figure 1 shows the results of the 10-fold cross validation process for the electricity, gas, and water consumption model s. The overall low variation of CV-prediction accuracy for all models indicates that the performance of the model s is consistent . Figure 2 shows the actual versus the predicted electricity, gas, and water consumption values using XG Boost . Moreover, the goodness of fit R2 of the three models are 0.80, 0.60, and 0.64 for electricity, gas, and water consumption respectively. The higher accuracy of the electricity consumption model is expected considering the independent variables selected . Indeed, gas consumption depends hea vily on whether gas is used for building heating or not ; in New York City, a significant number of large -scale buildings use steam for heating (not studied here) instead of gas . As a result, o ur model predicts the same gas consumption value for t wo buildings with the exact same characteristics, even if one uses gas for heating and the other uses electricity , hence the lower accuracy. As for water consumption model , water consumption also depends on other type s of variables that could not be included in this study, such as the frequency and duration of taking a shower / bath, whether a unit has a washing machine and dishwater or not, and so on, hence the lower accuracy. As a comparison, Papadopoulos and Kontokosta (2019) get a R2 of 0.31 using XGBoost to predict energy use intensity (EUI) of New York City large -scale buildings . Our results are significantly higher specifically because of we consider electricity, gas, and water consumption together, suggesting that these three vari ables are strong interrelated; more on interrelationships later. 10 Figure 1. Ten -fold cross -validation results for electricity, gas, and water consumption models . 11 Figure 2. Actual vs p redicted values for the electricity, gas, and water consumption models 12 Figure 3 shows the SHAP summary plots for predicting electricity, gas, and water consumption. The x-axis shows the impact of features on the outcomes based on the SHAP values. Features are sorted based on their impact , thus the variable at the top has the highest impact . The color represents the feature value, with red for high values and blue for low values. The vertical dispersion for each feature corresponds to the data point s with the same SHAP values. If red points are mostly present on the positive side of SHAP values, it means by increasing the value of the independent variable , the dependent variable increases as well. Figure 3.a shows that the most important variable when it comes to electricity consumption is whether a building is a residential building or not , and that by changing the land use type of a building from residential to something else , electricity consumption would increase. The three next important features are water consumption, gas consumption, and building age. Water and gas consumption have an overall a positive effect on electricity consumption , which means that a building that consumes more water and gas is likely to consume more electricity as well. Here, it is likely that water consumption acts as a proxy for building occupancy (which was not available as an independent variable) —an office building, for example, will have more occupants consuming both more electricity and water. A positive and significant relationship between electricity and water is therefore expected. In contrast, building age shows to have a negative relationship with electricity consumption. In other words, newer buildings tend to consume more electricity , which might reflect a higher usage of high energy consuming appliance s such as air conditioner s (EIA 2005) . For gas consumption (Figure 3.b), water and electricity consumption have the highest impact s, and by increasing water and electricity consumption, gas consumption generally increases . The third most important feature is building density, and by increasing building density , gas consumption is likely to decrease. This result is probably linked to the fact that higher building density suggests the presence of larger buildings that might have lower U-values and thus require less energy per unit area for heating. Figure 3.c shows the results for water consumption. We can see that gas and electricity consumption ha ve the highest impact in the positive direction , which means here again that buildings wit h higher electricity and gas consumption also generally consume more water . Here again, electricity likely acts as a proxy for building occupancy, and the positive interrelationship between electricity and water is expected. Furthermore, using the same dataset, Kontokosta and Jain (2015) also found a strongly relationship between electricity and gas consumption (in the form of energy use intensity) and water consumption. The third feature is college education (i.e., percentage of people with a college educat ion) and it is in the opposite direction , which means that buildings with more educated people consume proportionally less water. This feature likely relates to the presence of median income as well in the seventh place, similarly captured in Kontokosta and Jain (2015) here as well . Interestingly, we also see that building age has a negative impact on water consumption; i.e., newer buildings tend to consume comparatively more water . Finally, Figure 3.b and 3.c shows the most important variable for gas consumption model is water consumption, and the most important parameter for water consumption is gas consumption. Water and gas consumption are highly interrelated, which can be partly explained by the fact that gas is often used for water heating ; thus, a higher hot water consumption translates directly to a higher gas consumption. 13 Figure 3. Contribution of the most important features in each consumption. To further study the interrelationships between electricity, gas, and water consumption, Figure 4 shows plots with the value of a feature on the x-axis and its SHAP value on the y -axis by 14 changing a specific feature in the model . Considering the fact that SHAP values indicate a feature’s impact on a change in the model output, this type of figure s shows the impact of one feature on the model output, while the change in the second feature is displayed by changing the color of the data points from blue to red. While these plots are not trivial to analyze, they offer a n effective means to look at the interrelationships between three features. Figure 4.a shows the impact of water and gas consumption on electricity consumption. In Figure 4.a .2, using the electricity consumption model, positive SHAP values can be seen for buildings with average to high water consumption, and data points tend to be blue , thus have a low gas consumption. This means that the interrelationship between electricity and water is particularly strong in buildings with low gas consumption but with average to high water and electricity consumption. This can be seen in Figure 4.c.1 as well that applies the water consumption model , and we can see that in buildings with average to high electricity consumption, SHAP values are higher in buildings with low gas consum ption values . Buildings with average to high electricity consumption and low gas consumption likely use electricity for space heating and water heating . Conside ring the relation ship between water consumption and water heating, buildings with average to high electricity consumption probably use electricity for water heating, which can be why electricity and water are more significantly interrelated . Figure 4.b illustrates the impact of electricity and water consumption on gas consumption. Specifically, f igure 4.b.1 shows a peak SHAP value for an electricity consumption value of about 300 kWh/m2/y, which suggests that electricity consumption is increasingly interrelated with gas consumption up to a point, after which the interrelationship decreases . Moreover , the changing SHAP values between buildings with high and low water consumption suggests that electricity consumption has a higher impact on gas consumption for buildings with high water consumption. Figure 4.b.2 not only shows a direct relationship between water and gas consumption, it shows that this relationship becomes more important as water consumption increases since SHAP values increase sim ultaneously. Figure 4.c.2 echoes this observation as after a gas consumption of about 200 kWh/m2/y, an increase in gas consumption almost directly lea ds to an increasing SHAP value. 15 Figure 4. Interrelationships between electricity, g as, and w ater consumption 16 Figure 5. Impact of two features on the model output 17 We can then further study the impact of other features on electricity gas, and water consumption. Figure 5.a displays the results of water consumption as a feature to determine its impact on electricity consumption when the building type is changed from residential to non- residential. In the figure, t he red points are the residential buildings and the blue points are the nonresidential buildings. For nonresidential buildings (bl ue points), by changing water consumption, SHAP values change more significantly than for residential buildings (red points). Essentially, the figure shows that water consumption has more of an impact on electricity consumption in nonresidential buildings (e.g., in office buildings ). This is interesting because it suggests that water consumption not only capture s elements of building occupancy (as mentioned before) , it also captures some elements of consumption behavior (as building type has an impact on electricity consumption) . Furthermore, Figure 5.b shows the impact of college education and water consumption on electricity consumption. Here, the figure shows that by changing water consumption, the SHAP values change more for the red points than the blue points, where the red points represent buildings with a higher percentage of people that have a college education. Therefore, in buildings with more educated people, higher water consumption patterns tend to translate into higher electricity consumption as well. Finally, f igure 5.c shows the impact of building density and electricity consumption on gas consumption. Despite some significant noise, we can see that SHAP values are mostly positive and do not vary significantly until a building density of 2 m 2/m2. We can then see a significant decrease in SHAP value from approximately 0 to - 75 for building densities ranging from 2 to 4 m2/m2, after which SHAP values become stable again . As mentioned, gas is often us ed as the primary source of heating, and higher building densities corresponds to larger buildings that tend to have lower U-values. We note that the rate of SHAP values reduction mostly stop after a building density of 4 m2/m2. Moreover, the noise suggests that building density is one out of many important variables that impact gas consumption (others include building insulation material, resident behavior, and so on (Haas, Auer, and Biermayr 1998; Guerra Santin, Itard, and Visscher 2009) . In contra st, we do not observe any pattern between electricity consumption and building density using the gas consumption model . To summarize the results, in this work , we found that : • water and gas consumption tend to be strongly related to electricity consumption, but electricity consumption is slightly less sensitive to water and gas consumption . While this observation is interesting, we must acknowledge that it could be an artifice of the models since many important variables were not available, which also explai ns the lower accuracy of the water and gas consumption models . For example, data on building and water heating technology (i.e., electricity, gas, or steam) w ere not provided. • water and gas consumption tend to be significantly interrelated , which is likely because gas is often used for water heating. • electricity and water consumption are particularly interrelated in buildings with low gas consumption, again likely electricity is likely used for water heating in buildings with low gas consumption. • electricit y and gas consumption are interrelated up to a point after which they are not as interrelated , likely linked with the technology used for heating 18 • in buildings with low electricity consumption, water and gas consumption are negatively interrelated ; these buildings tend to be warehouses and other building types where an increase in gas consumption per square meter has no impact on water consumption. • water consumption has more of an impact on electricity consumption in nonresidential buildings ; this suggests that water consumption not only act s as a proxy to building occupancy, but it also captures some elements of consumption behavior. • after an initial threshold, increasing building density generally lowers gas consumption significantly , but only until up to a point after which the relationship plateaus; this suggests that building intensification is effective only after a certain threshold and only up to a specific building density. 5. Conclusion Understanding the consumption patterns and interrelationships between energy and water consumption in buildings is paramount to be able to meet future demand in a world with a growing population. In this work , we looked at building energy and water consumption data in New York city and merged this data with the Primary Land Use Tax Lot Output (PLUTO) dataset as well as with other socio -economic data from the American Community Survey. The final dataset allowed us to develop three models: (1) electrici ty, (2) gas, and (3) water consumption. The results show that the all three models perform relatively well, especially the electricity consumption model ; gas is used heavily for heating and building heating technology was not available, and water consumption depends on other variables that were not available. Next, we determined the importance of each variable in each model, and we calculated how the variation in these features could affect energy and water consumption by using SHAP . Subsequently, we used SHAP to analyze how changing two features at the same time could affect the energy and water consumption. In the end, we listed seven major lessons learnt from this study such that electricity and water consumption are particularly interrelated in buildings with low gas consumption. Arguably, the largest limitations of this study are the fact that only large -scale buildings are included and the specific results are not generalizable beyond the context of New York Cities. Ideally, future work would include small -scale buildings (ideall y single detached homes) located in different climate zones to fully investigate interrelationships between electricity, gas, and water consumption. Nonetheless, the knowledge gained from the models and from the SHAP analysis can provide some insights to help develop data-driven polic ies, especially in cities with large -scale buildings, and evaluate different policy alternatives to find out how effectively energy and water consumption can be managed in a sustainable way. 6. Acknowledgments This research i s partly supported by the National Science Foundation (NSF) CAREER award #1551731. 19