Statistical Data Analysis for Energy Communities

Reducing and rationalizing electric consumption for residential use is an important issue in the agenda of governments and decisions makers, and studies have tried to understand the main environmental and socio-demographic variables that influence electric energy consumptions at different scales.

In 2009, Olofsson et al. [8] analysed the effect of building-specific features on both thermal and electric energy consumption using analysis of variance (ANOVA) tests. They stated that the most important variables were the geometrical characteristics and the construction time of the buildings. However, Howard et al. [9] in 2012 applied a multiple linear regression model that evidenced the higher importance of building function than the construction type or age of construction.

Two studies conducted in Rotterdam (Netherlands) [10, 11] in recent years, using a multiple linear regression model to investigate both natural gas and electric consumption, pointed out that for electricity number of occupants, floor surface and type of dwelling are the most significant variables. Furthermore, Nouvel et al. [11] compared an engineering method to a statistical one to predict electric consumption from average floor area of dwellings, average number of occupants and the share of dwellings. They showed that the deviation from real energy consumption data was lower using a statistical approach.

In 2019, Mutani et al. [12] evaluated the thermal consumption of building in the urban context of Torino (Italy). They implemented a multiple linear regression model, using building-specific features as well as environmental and socio-demographic variables, that allowed to highlight major drivers of the energy consumption as area of dwellings, number of occupants, but also age of the residents and their social characteristics.

Chen et al. [13], used a multiple regression analysis to investigate the relationship among household variables and residential energy consumption for space heating and cooling. Data were obtained from more than 1,500 surveys applied to households in the city of Hangzhou (China). The results pointed out that up to 28.8% of the energy consumption could be explained by socio-economic and behavioural variables, and floor area accounts for 44%.

Bianco et al. [14] analysed residential and non-residential annual electric consumption in Italy, during the period 1970–2007. They implemented simple and multiple regression models using historical data on energy consumption, gross domestic product (GDP) and GDP per capita and demographic data.

In 2013, Gans et al. [15] used data on real-time usage from Continuous Household Survey of Northern Ireland, which surveyed about 300 households per month, to investigate the effect on residential electric consumption. The surveys collected socio-economic data as dwelling, health, education, employment, and welfare payments. The regression model implemented accounted electricity price, household income as well as weather, house-specific features, and type of heating system.

Nie and Kemp [16] investigated trends in energy consumption in China during the 2002-2010 years, obtaining data from open-access institutional databases. They evaluated the effect of changes in demographic parameters, house-specific characteristics, technology systems and energy sources. The key drivers for energy consumption increase were appliances, floor space per capita, whereas the energy mix and population were the less significant. They tried to forecast electric use using the regression method, predicting that the electric consumption would continue to rise despite a partial saturation of the demand.

Social, demographic, and economic variables, along with environmental and building-specific features are useful to understand electric consumption for the residential sector and even to forecast future trends. On the other hand, in-field collection of data from surveys and measurement are costly and time-consuming and cannot be performed on a vast scale. Therefore, the aim of this study is to identify the most influential variables on electric energy consumption for residential sector at a municipal scale, starting from public-access data.

To investigate major drivers of residential electric consumption, a GIS-based methodology with the implementation of statistical techniques has been developed in this study. The achievement of the best accuracy depends on the reliability of data and databases. In this study, to enhance pragmatism and ensure reproducibility, variables were obtained from institutional open-access databases [17-19], geographical and environmental data were calculated using well established procedure with GIS software [20, 21], while electric consumption observations were obtained from the database provided by the Piedmont Region technical office, which collects energy data directly from the Distribution System Operators (DSOs) [22].

The software used was R version 4.0.0 with a database composed by 1,201 observations and 117 variables related to residential electric consumptions. The observations consist in the yearly consumption of 1,201 municipalities from 2010 to 2016 (6 years), excluding the year 2011 because of missing data. Variables were collected only once for each municipality, assuming that during the time span considered they did not vary significantly. The methodological framework of this study can be summarized in the following steps:

1. Description of the statistical distribution of energy consumption data.
2. Identification of the independent variables.
3. Utilization of univariate and multivariate analysis technique (principal component analysis) to improve data quality.
4. Identify a linear regression model.
5. Evaluation of the influence of independent variables on the electric consumption.                     

4.1 Statistical distribution of energy consumption data

Considering 1,201 municipalities, statistical analysis was performed to evaluate the frequency distribution of energy consumption data separately for each of the six-year span. Anomalous observations were individually evaluated, comparing the consumption of each year to the previous and the next ones. Cut-off for maximum variation was set arbitrarily to 0.3. Observation that varied more than 30% were eliminated. The final database consists of 1,186 municipalities selected. Cullen and Frey chart has been used to inspect plausible probability distribution.

The Normal and Log-Normal distributions were checked by graphically evaluating the goodness of fit, drawing probability density functions and histogram together [23]. The ANOVA test has been used for evaluating whether the electric consumption differs significantly among years.

4.2 Identification of independent variables

Residential electric energy consumption of 1,186 municipalities were georeferenced using a GIS tool, to combine the energy consumption data with the socio-demographic and environmental characteristics of the municipalities. The socio-economic and built environment characteristics refer spatially to the census sections, which represent the territorial units commonly used to collect data. In this work, the 2011 ISTAT census database has been managed using GIS tool referring all available information to the municipal area. Some variables have been transformed: the values of some observations (absolute values) have been calculated in percent of the total value, new variables were created from the former variables.

4.3 Univariate and multivariate analysis techniques

First, univariate statistical methods were used to describe and investigate variables and their relationship with the outcome. Pearson’s rho [r] was used to assess the level of correlation among variables and between each variable and the electric consumption data. A correlation matrix was used to identify collinear variables and remove them from subsequent analysis to avoid singularity of the correlation matrix itself.

The principal component analysis (PCA) was used to identify variables that are able to explain the majority of data variability and to divide observations into homogeneous groups. The Kaiser-Meyer-Olkin factor has been referred as a measure of adequacy of PCA. Hierarchical clustering has been performed using on the selected principal components the criterion of Ward that is based on the multidimensional variance. To obtain homogeneous groups the cluster dendrogram, obtained from the PCA, has been cut off at the height of 1.0. The contribution of variables to main principal components was graphically evaluated to identify uninfluential variables.

Multiple linear regression models

The linear regression model is the most used statistical technique for investigating and modelling the relationship between a dependent and two or more independent variables. The electric energy consumption was estimated using a linear regression, which is expressed by Eq. (1):

$Y_{i}=\beta_{0}+\beta_{1} \cdot x_{i 1}+\cdots+\beta_{p} \cdot x_{i p}+\varepsilon_{i}$            (1)

where, Y_i is the annual electric consumption (dependent variable), x_ij are the independent variables, $\beta_{\mathrm{i}}$ are the estimated coefficients and ε_i correspond to the random errors of each observation i, i=1, …, N. The standard errors must be independent and normally distributed with mean 0 and constant variance: ε_i ~IIND (0, σ²).

The Ordinary Least Squares (OLS) method was used to identify the model express by Eq. (1). The observed values y_i (i=1, …, N) can be written as Eq. (2):

$y_{i}=b_{0}+b_{1} \cdot x_{i 1}+\cdots+b_{\mathrm{p}} \cdot X_{\mathrm{ip}}+e_{i}$            (2)

where, b_j are the least squares estimates of β_j (j=0,1, …, p) and e_i (i=1, …, N) are the residuals. Predicted values y_i are computed as b₀ + b₁ ∙ x_i1 + ⋯ + b_p ∙ x_ip (i=1, …, N).

Linear dependency, no multicollinearity among variables, normality of ε_i, and homoscedasticity among ε_i are the conditions on a multiple linear regression model.

The stepwise method in both directions (backward and forward) was performed to identify the linear model minimizing the Akaike information criterion (AIC), that is asymptotically equivalent to cross validation, and selecting the most influential variables. The stepwise method is an automatic selection procedure which combines backward elimination steps (computing the t-ratio for each regressor in a subset and eliminating the ones that have absolute value of the t-ratios smaller than a prespecified value) and forward selection steps (adding a new variable if the corresponding t-ratio is the largest and its value is greater than a prespecified value).

The models resulted from the stepwise must be checked, evaluating the residuals for their normality and homoscedasticity, and the variables for multicollinearity by their Variance Inflation Factor (VIF).

Homoskedasticity (homogeneous variance for residuals) is a crucial assumption in regression analysis. It was tested through the scatter plot of residuals predicted values, normality, and probability graphs (Q-Q-plot). The outlier and leverage diagnostic graphs have been used to identify observations that are influential points. Observations, whose residuals (Student residuals) and leverage values were higher than 2 standard deviations, have been selected and removed from the dataset and the analysis was performed again.

The appropriate diagnostics (e.g., Breusch–Pagan test and White test) were performed to check the accuracy of the model prediction. To evaluating homoscedasticity of residuals the White test has been conducted, considering the following hypotheses:

$H_{0}: \sigma_{i}^{2}=\sigma^{2}$         (3)

$H_{1}: \exists_{i, j}$ such that $\sigma_{i}^{2} \neq \sigma_{j}^{2}$         (4)

where, the null hypothesis (Eq. 3) implies that residuals have constant variance (σ²) and it is similar across all the values without showing any pattern, while the alternative hypothesis (Eq. 4) signifies a different variance among them. In this latter case a variance-stabilizing transformation is required to obtain more accurate variables estimators of the model.

The presence of heteroscedasticity in the final model was reduced through Heteroskedasticity-Consistent Standard Errors (HCSE or White correction), that allowed the fitting of a linear regression model that contains heteroscedastic residuals [24].

The multicollinearity of a multiple regression model makes difficult to understand how much the dependent variable is affected by the independent variables since they are all influencing each other. The validity of the multiple regression was assessed by the Variance Inflation Factors (VIF). It is calculated as the reciprocal of the inverse of R_j² (where R² is the coefficient of determination) of an independent variable x_j as it is expressed by Eq. (5):

$V I F_{j}=\frac{1}{1-R_{j}^{2}}$        (5)

The VIF gives the proportional increase in the variance of $\beta_{i}$ with respect to what it would have been if the explanatory variables were completely uncorrelated. The evaluation criterion refers to the threshold calculated as Eq. (6)

$V I F<\max \left(10, \frac{1}{1-R_{\text {model }}^{2}}\right)$               (6)

where, $R_{\text {model}}^{2}$ is the usual R-squared of the regression model. Variables with high VIF have been transformed or redefined and analysis was performed again.

The coefficient of determination R² correspond to the proportion of the variance in the dependent variable that is predictable from the independent variables and it was used to evaluate the adequacy of the model. The adjusted R² (Adj R²) has been adjusted for the number of predictors in the model. R² increases with the addition of variables in the model, while Adj R² increases only if that addition improves the model more than would be expected by chance.

In this study, the transformation of variables in order to improve the accuracy of the model consisted in three types of intervention.

1. Some absolute variables have been normalized according to Eq. (7):

$z_{i j}=\frac{x_{i j}-\bar{x}_{j}}{\sigma_{j}}$         (7)

where, $\bar{x}_{j}$ is the mean of the j^th variable and σ_j its standard deviation {x_ij: i=1, …, N};

2. Some variables were squared;
3. Some variables have been transformed and redefined.

Particularly, this latter is the case for the variables related to three characteristics of the population (Education degree, Family members and Age) and one characteristic of residential building (Building Construction Year).

Education degree variables. The 4 initial variables selected for the linear model consisted of the percentage of the population with different education levels: primary school diploma (P), junior high school diploma (JH), senior high school diploma (SH) or university degree (U); these variables were calculated as % values relative to the total population Ptot of each municipality. Each of them has been associated with the total number of years of study envisaged by the Italian school system, respectively: 5, 8, 13 and 18 years. Then the 4 variables have been redefined in a single variable (EDU) that expresses the average number of years of education of the population of the municipality, calculating the weighted average as Eq. (8):

$E D U=\frac{\operatorname{Ptot} *[(P * 5)+(J H * 8)+(S H * 13)+(U * 18)]}{P t o t}$               (8)

Family members variables. The 5 initial variables selected consisted of the percentage of the families with 1, 2, 3, 4, or 5 members per family relative to the total number of families (Ftot) of each municipality. Then the 4 variables have been redefined in a single variable (FAM) that expresses the average number of members per family for each municipality, calculating the weighted average.

Age variables. The 3 initial variables selected for the initial model consisted of the percentage on the total population (Ptot) that had 0-19, 20-69 or 70-99 years old. Then the 4 variables have been redefined in a single variable (AGE) that expresses the average age of the population in each municipality, considering the median value for each of the 3 starting age categories and calculating the weighted average.

Building Construction Year variables. The 9 initial variables selected for the initial model consisted of the percentage on the total number of residential building in each municipality that have been built respectively before 1918, 1919-1945, 1946-1960, 1961-1970, 1971-1980, 1981-1990, 1991-2000, 2001-2005, after 2005. Then the 9 variables have been redefined in a single variable (BCY) that expresses the average age of residential buildings in each municipality. The oldness of building has been identified as the difference between the corresponding year of energy data (2015) and the median value for each of the 9 starting categories, then the weighted average has been calculated.

4.4 Influence of independent variables

A further verification of the validity of the model took place by checking whether the observed values fall within the interval of the predicted values. Firstly, the predicted values for the year 2015 and the respective 95% confidence interval were calculated. The prediction interval was then calculated. This information was graphically evaluated in comparison with the observed values of the annual electricity consumption for the year 2010, 2012, 2013,2014,2016. For each of the 5 years, the observations that fell outside the prediction interval were identified. The information relating to them (e.g., belonging to the clusters and other variables) was obtained and their distribution was assessed.