<jupyter_start><jupyter_text>Regression OutlineWe have a set of manually categorized FERC plants with which to test the sklearn Classifier that we're making. We'll also use them to test the regression analysis here. What does this dataset look like, how was it created?* Pull FERC Plants table.* Using a variety of search methods, to ensure that we get a wide variety of plants, identify sets of records in the FERC Form 1 Plants table that are comparable inter-year records.For regressions, both OLS and sklearn approaches are explored.  Plant generation is normalized by using capacity factor.<jupyter_code>%load_ext autoreload
%autoreload 2
import sys
import os
import numpy as np

import scipy
import pandas as pd
import statsmodels.api as sm
import sklearn
from sklearn import linear_model
from pudl import init, analysis, settings, output
import pudl.extract.ferc1
import pudl.transform.ferc1
import pudl.constants as pc
import pudl.output.pudltabl
import pudl.analysis.analysis
import matplotlib.pyplot as plt
import matplotlib as mpl
import timezonefinder

%matplotlib inline
# For some reason these things don't stick if they're in the same cell as the
# %matplotlib inline call above, but if they're separate, they work fine.
plt.style.use('ggplot')
mpl.rcParams['figure.figsize'] = (10,6)
mpl.rcParams['figure.dpi'] = 150
pd.options.display.max_columns = 56
pudl_out = pudl.output.pudltabl.PudlTabl()
ferc1_fuel = pudl_out.fuel_ferc1()
ferc1_steam = pudl_out.plants_steam_ferc1()
# eia860_plants =  pudl.output.eia860.plants_eia860()
# eia923_plants =  pudl.output.eia923.plants_eia923()<jupyter_output><empty_output><jupyter_text>Aggregate data from FERC 1 Fuel<jupyter_code># For each plant record in the plants table, we need to be able to assign
# a fraction of heat content for that plant that came from a particlar fuel
# for filtering purposes. There are different ways that we can aggregate that
# information -- e.g. plant_id_pudl, plant_id_ferc1, (utility_id_ferc1, plant_name)
def coal_fraction(fuel_df, gb_cols=['plant_id_pudl']):
    """
    Calculate the fraction of overall heat content from coal
    per year and per gb_cols, which can be:
      - [plant_id_pudl]
      - [plant_id_ferc1]
      - [utility_id_ferc1, plant_name]
    
    Assumes fuel_df is a fuel_ferc1 dataframe.
    Cannot currently work for arbitrary fuel_type becuase not all plants
    have records for all fuel types.
    """
    gb_cols1 = ['report_year'] + gb_cols
    
    # Total heat content across all fuel types
    total_mmbtu = (
        fuel_df.groupby(gb_cols1).
        agg({'fuel_consumed_mmbtu': sum}).
        reset_index().
        rename(columns={'fuel_consumed_mmbtu': 'total_mmbtu'})
    )
    
    # Same grouping... but also broken down by fuel type
    gb_cols2 = gb_cols1 + ['fuel_type_code_pudl']
    out_df = (
        fuel_df[fuel_df.fuel_type_code_pudl=='coal'].
        groupby(gb_cols2).
        agg({'fuel_consumed_mmbtu': sum}).
        reset_index().
        merge(total_mmbtu, how='left', on=gb_cols1).
        assign(coal_fraction=lambda x: x['fuel_consumed_mmbtu'] / x['total_mmbtu']).
        rename(columns={'fuel_consumed_mmbtu': 'coal_mmbtu'})
    )
    return out_df
ferc_cols_to_keep = ['report_year', 'utility_id_ferc1', 'utility_id_pudl',
                     'utility_name_ferc1', 'plant_id_pudl', 'plant_id_ferc1', 'plant_name', 'plant_type','installation_year',
                     'opex_per_mwh']
# Calculate total heat content of fuel by plant & year
fuel_consumed_mmbtu = (
    ferc1_fuel.loc[ferc1_fuel.report_year>=2012,:].
    groupby(['report_year', 'plant_id_pudl']).
    agg({'fuel_consumed_mmbtu': sum}).
    reset_index()
)

coal_fraction_by_plant_name = coal_fraction(ferc1_fuel.loc[ferc1_fuel.report_year>=2004,:], ['utility_id_ferc1', 'plant_name'])

# New coal fraction estimation using this function - still need to fix column name changes this induces
# coal_fraction_by_plant_name = pudl.transform.ferc1.fuel_by_plant_ferc1(ferc1_fuel)
coal_fraction_by_plant_name.columns
steam_agg_ferc1 = (
    # Only interested in 2004 and later
    ferc1_steam.loc[ferc1_steam.report_year>=2004,:].
    # This ensures we only include plants that use *some* coal
    merge(coal_fraction_by_plant_name, how='inner', on=['report_year', 'utility_id_ferc1', 'plant_name']).
    # Break out the non-fuel OpEx
    assign(opex_nonfuel=lambda x: x['opex_production_total'] - x['opex_fuel']).
    # Aggregate by plant_id_pudl
    groupby(['report_year','plant_id_pudl']).
    agg({'capacity_mw': sum,
         'net_generation_mwh': sum,
         'opex_fuel': sum,
         'opex_nonfuel': sum,
        'opex_allowances': sum,
        'opex_boiler': sum,
        'avg_num_employees': sum,
         'opex_electric': sum,
         'opex_engineering': sum,
         'opex_misc_power': sum,
         'opex_misc_steam': sum,
         'opex_operations': sum,
         'opex_plants': sum,
         'opex_production_total': sum,
         'opex_rents': sum,
         'opex_steam': sum,
         'opex_steam_other': sum,
         'opex_structures': sum,
         'opex_transfer': sum       
        }).
    reset_index().
    # Bring in total heat content per pudl_plant
    merge(fuel_consumed_mmbtu, how='left', on=['report_year', 'plant_id_pudl']).
    # Calculate some other interesting metrics:
    assign(opex_fuel_per_mwh=lambda x: x['opex_fuel'] / x['net_generation_mwh'],
           opex_nonfuel_per_mwh=lambda x: x['opex_nonfuel'] / x['net_generation_mwh'],
           heat_rate_mmbtu_mwh=lambda x: x['fuel_consumed_mmbtu'] / x['net_generation_mwh'],
           capacity_factor=lambda x: x['net_generation_mwh'] / (8760*x['capacity_mw']),
           opex_total_per_mwh=lambda x: x['opex_nonfuel_per_mwh'] + x['opex_fuel_per_mwh']).
    merge(ferc1_steam[ferc_cols_to_keep],how='left',on=['report_year','plant_id_pudl',])
)
coal_plants = coal_fraction_by_plant_name[coal_fraction_by_plant_name.coal_fraction >= 0.85] # returns ?% of records
coal_plants[coal_plants.plant_name=='dolet hills']
## This function call replaced by coal fraction calculation
# fuel_df_plant = ferc1_aggregate_by_fuel_type(ferc1_fuel, fuels=['coal'])
new = pd.merge(steam_agg_ferc1, coal_plants, on=['plant_name', 'utility_id_ferc1', 'report_year'], how='right')
new.head(2)
# new.groupby(['report_year','plant_id_pudl', 'capacity_mw', 'net_generation_mwh', 'opex_fuel',
#             'opex_nonfuel','avg_num_employees','opex_production_total','fuel_consumed_mmbtu', 'opex_fuel_per_mwh',
#             'opex_nonfuel_per_mwh','heat_rate_mmbtu_mwh','capacity_factor','opex_total_per_mwh', 'installation_year'], as_index=False).agg({'coal_fraction_mmbtu':sum,
#          'fuel_mmbtu':sum}, as_index=False)

# With previous version of coal_fraction function
new.groupby(['report_year','plant_id_pudl', 'capacity_mw', 'net_generation_mwh', 'opex_fuel',
            'opex_nonfuel','avg_num_employees','opex_production_total','fuel_consumed_mmbtu', 'opex_fuel_per_mwh',
            'opex_nonfuel_per_mwh','heat_rate_mmbtu_mwh','capacity_factor','opex_total_per_mwh', 'installation_year'], as_index=False).agg({'coal_mmbtu': sum,
         'total_mmbtu':sum}, as_index=False);<jupyter_output><empty_output><jupyter_text>pd.options.mode.chained_assignment = None   default='warn'<jupyter_code>new.installation_year = pd.to_numeric(new.installation_year, errors='coerce')
# What types of fuel are included in the training data set?
new.plant_type.unique()
# Look at the anomolous nuclear plant - it's actually coal
new[new.plant_type=='nuclear']
# Unique plants
# steam_plants = new.sort_values('installation_year').drop_duplicates(subset=['capacity_mw', 'plant_name'], keep='last') # 1007 'plants'
# ct = out_df_ct.sort_values('installation_year').drop_duplicates(subset=['plant_id_ferc1'], keep='last') # 584 'plants'
# cc = out_df_cc.sort_values('installation_year').drop_duplicates(subset=['plant_id_ferc1'], keep='last') # 188 'plants'
# ic = out_df_ic.sort_values('installation_year').drop_duplicates(subset=['plant_id_ferc1'], keep='last') # 147 'plants'
# nuc = out_df_n.sort_values('installation_year').drop_duplicates(subset=['plant_id_ferc1'], keep='last') # 5 'plants'
# pv = out_df_pv.sort_values('installation_year').drop_duplicates(subset=['plant_id_ferc1'], keep='last') # 3 'plants'

# # Remove null values for variables of interest, for purposes of histogram plotting
# steam_plants = steam_plants[steam_plants.construction_year.notnull()]
new = new[new.installation_year.notnull()]
new = new.loc[(new['installation_year'] > 1900)]
new = new.loc[(new['capacity_factor']> 0)]

# steam_plants = steam_plants[steam_plants.capacity_mw.notnull()]
# steam_plants = steam_plants[steam_plants.heat_rate.notnull()]
# Look at distribution of construction years vs. installation years.
fig = plt.figure()

plt.xlabel('Year')
plt.ylabel('Number plants')
plt.title('FERC 1 Steam Plants')

# mpl.pyplot.hist(new.construction_year, bins = 117, range = [1900,2016], alpha=0.5, label='Construction year')
mpl.pyplot.hist(new.installation_year, bins = 117, range = [1900,2016], alpha=0.5, label='Installation year')
mpl.pyplot.legend(loc='upper right');
# Look at distribution of construction years vs. installation years.
fig = plt.figure()

plt.xlabel('Non-fuel costs')
plt.ylabel('Number of plants?')
plt.title('FERC 1 Steam Plants')

# mpl.pyplot.hist(new.construction_year, bins = 117, range = [1900,2016], alpha=0.5, label='Construction year')
mpl.pyplot.hist(new.capacity_mw, bins = 100, alpha=0.5, label='Installation year')
mpl.pyplot.legend(loc='upper right');
plt.xlabel('Installation year')
plt.ylabel('Capacity (MW)')
plt.title('FERC 1 Steam Plants')
# plt.legend()

mpl.pyplot.hist2d(new.installation_year, new.capacity_mw, bins=117);

# fig.colorbar(im);
#1940 and forward, up to 3K megawatts, just steam - natural groupings?
# install, constru, capacity, and HR then we should find groupings to work wtih.
new.plant_id_pudl.nunique()<jupyter_output><empty_output><jupyter_text>Manual binning of construction year, installation year - regressions with statsmodel<jupyter_code>oldest = new.loc[(new['installation_year'] < 1956)]
# older = new.loc[(new['installation_year'] >= 1946) & (new['installation_year']< 1956)]
mid = new.loc[(new['installation_year'] >= 1956) & (new['installation_year']< 1966)]
newer = new.loc[(new['installation_year'] >= 1966) & (new['installation_year']< 1986)]
newest = new.loc[new['installation_year'] >= 1986]
res = sm.OLS.from_formula(formula='opex_nonfuel ~ net_generation_mwh + capacity_mw', data=new).fit()

# res = sm.OLS.from_formula(formula='opex_nonfuel ~ net_generation_mwh', data=oldest).fit()
# res = sm.OLS.from_formula(formula='opex_nonfuel ~ capacity_mw', data=oldest).fit()


# print(res.params)
print("R-squared = ")
print(res.rsquared_adj)
print("")
print(res.params)
print(res.summary())
res = sm.OLS.from_formula(formula='opex_nonfuel ~ net_generation_mwh + capacity_mw', data=mid).fit()

# res = sm.OLS.from_formula(formula='opex_nonfuel ~ net_generation_mwh', data=mid).fit()
# res = sm.OLS.from_formula(formula='opex_nonfuel ~ capacity_mw', data=mid).fit()

# print(res.params)
print("R-squared = ")
print(res.rsquared_adj)
print("")
print(res.params)
print(res.summary())

res = sm.OLS.from_formula(formula='opex_nonfuel ~ net_generation_mwh + capacity_mw', data=newer).fit()

# res = sm.OLS.from_formula(formula='opex_nonfuel ~ net_generation_mwh', data=newer).fit()
# res = sm.OLS.from_formula(formula='opex_nonfuel ~ capacity_mw', data=newer).fit()

# print(res.params)
print("R-squared = ")
print(res.rsquared_adj)
print("")
print(res.params)
# res = sm.OLS.from_formula(formula='opex_nonfuel ~ net_generation_mwh + capacity_mw', data=newest).fit()

# res = sm.OLS.from_formula(formula='opex_nonfuel ~ net_generation_mwh', data=newest).fit()
res = sm.OLS.from_formula(formula='opex_nonfuel ~ capacity_mw', data=newest).fit()


# print(res.params)
print("R-squared = ")
print(res.rsquared_adj)
print("")
print(res.params)<jupyter_output>R-squared = 
0.8478313804037859

Intercept     -2.156691e+07
capacity_mw    7.756056e+04
dtype: float64<jupyter_text>sklearn regressions<jupyter_code>from sklearn import preprocessing
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
train, test = train_test_split(new,test_size=0.2)
# train.fillna(train.mean(), inplace=True)
# train=train.astype('float16')

# test.fillna(test.mean(), inplace=True)
# test=test.astype('float16')

features = ['installation_year', 'heat_rate_mmbtu_mwh', 'capacity_mw', 'capacity_factor', 'opex_fuel_per_mwh','net_generation_mwh']
#features = ['installation_year', 'capacity_mw', 'capacity_factor', 'net_generation_mwh']

y_var='opex_nonfuel'
x_train = train[features]
x_train.fillna(x_train.mean(), inplace=True)
x_train=x_train.astype('float32')
# divided the dollar amount by a million
y_train = train[y_var]/1000000
y_train.fillna(y_train.mean(), inplace=True)
y_train=y_train.astype('float32')

# for testing
x_test = test[features]
x_test.fillna(x_test.mean(), inplace=True)
x_test= x_test.astype('float32')

y_test = test[y_var]/1000000
y_test.fillna(y_test.mean(), inplace=True)
y_test = y_test.astype('float32')<jupyter_output>/home/alana/anaconda3/lib/python3.6/site-packages/pandas/core/generic.py:6120: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  self._update_inplace(new_data)<jupyter_text>scaling the variables<jupyter_code>min_max_scaler = preprocessing.MinMaxScaler()
# make sure that the scaler is applied to each col separately
x_train_scaled = min_max_scaler.fit_transform(x_train)
x_test_scaled = min_max_scaler.transform(x_test)
(x_train_scaled[:,1]).mean()
reg = linear_model.LinearRegression()
reg.fit(x_train_scaled,y_train)

f=reg.coef_
#f=regressorAirTemp.feature_importances_

# Returns the coefficient of determination R^2 of the prediction.
r2_model2=reg.score(x_test_scaled, y_test)

# printing the results of the regression
print ('R^2 of the prediction is {0}'.format(r2_model2))
for i in range(0,len(f)):
    print ('The importance of ', features[i],'-------',round(f[i],3))
## Extra random forest
# Overfits the data?

from sklearn.ensemble import ExtraTreesRegressor
reg = ExtraTreesRegressor(n_estimators=200,min_samples_split=5,max_depth=200,random_state=0)
reg.fit(x_train_scaled, y_train)

f=reg.feature_importances_

# Returns the coefficient of determination R^2 of the prediction.
r2_model2=reg.score(x_test_scaled, y_test)

# printing the results of the regression
print ('R^2 of the prediction is {0}'.format(r2_model2))
for i in range(0,len(f)):
    print ('The importance of ', features[i],'-------',round(f[i],3))
# Unscale the data
# y_new_inverse = scalery.inverse_transform(y_new)<jupyter_output><empty_output><jupyter_text>statmodel<jupyter_code># model = sm.OLS(y_train,x_train_scaled)
X2 = sm.add_constant(x_train_scaled)
model=sm.OLS(y_train, X2)
results = model.fit()
print(results.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:           opex_nonfuel   R-squared:                       0.656
Model:                            OLS   Adj. R-squared:                  0.655
Method:                 Least Squares   F-statistic:                     1091.
Date:                Thu, 06 Jun 2019   Prob (F-statistic):               0.00
Time:                        07:55:14   Log-Likelihood:                -18217.
No. Observations:                3438   AIC:                         3.645e+04
Df Residuals:                    3431   BIC:                         3.649e+04
Df Model:                           6                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|[...]<jupyter_text>Further/future regression exlorations - scratch<jupyter_code># # Set up variables for OLS regression 

# # X = new[["net_generation_mwh", "heat_rate"]]
# X = new[["capacity_mw", "net_generation_mwh"]]
# # X = new[["non_fuel_costs"]]
# # X = new[["expns_fuel"]]

# # y = new["opex_fuel"]
# # y = new[["non_fuel_costs"]]
# y = new[["opex_nonfuel"]]

# # y = new[["net_generation_mwh"]]


# X = sm.add_constant(X) ## add an intercept

# model = sm.OLS(y, X).fit()
# predictions = model.predict(X) # make the predictions by the model

# # Print out the statistics
# model.summary()
# # Use sklearn linear_model to fit the same regression
# lm = linear_model.LinearRegression()
# model = lm.fit(X,y)
# predictions = lm.predict(X)
# # print(predictions)[0:5]
# lm.score(X,y)
# scipy.stats.linregress(out_df.net_generation_mwh, out_df.non_fuel_costs)
# from sklearn.linear_model import LinearRegression, Lasso, Ridge, ElasticNet, SGDRegressor
# from sklearn.model_selection import train_test_split
# x = out_df.net_generation_mwh
# y = out_df.non_fuel_costs
# x_reshape = out_df_s.net_generation_mwh.values.reshape(-1,1)
# y_reshape = out_df_s.non_fuel_costs.values.reshape(-1,1)
# intercept, slope, corr_coeff = plot_linear_regression(x_reshape, y_reshape)
# plt.show()

# plt.plot(x,y,'.')
# plt.title('FERC 1 steam records')
# plt.xlabel('net generation (mwh)')
# plt.ylabel('non-fuel costs')
### Exploration of other approaches...
# x_train, x_test, y_train, y_test = train_test_split(x, y, test_size= 0.2, random_state=0)
# x_train= x_train.values.reshape(-1, 1)
# y_train= y_train.values.reshape(-1, 1)
# x_test = x_test.values.reshape(-1, 1)
# ridge = Ridge(fit_intercept=True, alpha=0.5)
# ridge.fit(x_reshape,y_reshape)
# sklearn.linear_model.Ridge(alpha=1.0, fit_intercept=True, normalize=False, copy_X=True, max_iter=None, tol=0.001, solver='auto', random_state=None)
# ferc1_training_steam.net_generation_mwh, ferc1_training_steam.non_fuel_costs)
# clf = sklearn.linear_model.Ridge(alpha=1)
# clf = sklearn.linear_model.Ridge(alpha=1, copy_X=True, fit_intercept=True, max_iter=None, normalize=False, random_state=None, solver='auto', tol=0.001)

# clf.fit(x_reshape, y_reshape)
# Ridge(alpha=1, copy_X=True, fit_intercept=True, max_iter=None,normalize=False, random_state=None, solver='auto', tol=0.001)
# result = clf.fit(x_reshape, y_reshape, sample_weight=None)
# clf.score(x_reshape,y_reshape)
# result.get_params()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13736677790935514, 'fraction_numerical': 0.021207880288513295, 'mean_word_length': 4.950992953235105, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '145356', 'license': 'MIT', 'n_tokens_mistral': 7115, 'n_tokens_neox': 6565, 'n_words': 1282, 'path': 'results/notebooks/cpi_ferc_plants/regressions.ipynb', 'repo_name': 'claytonpbarrows/pudl', 'size': '18578'}