{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# **Generalized Additive Model (GAM) for Predicting Smart Homes Indoor Temperature**\n## ** Contents:**\n\n* [1. Import the libraries](#Import)\n* [2. Dataset Information](#Dataset)\n* [3. Model Identification](#Identification)\n > * [3.1 Adversarial Validation](#Adversarial)\n > * [3.2 Normality tests & Distribution plots](#Distribution)\n* [4. Feature Engineering & Data preprocessing](#Engineering)\n > * [4.1 Feature Extraction](#Extraction)\n > * [4.2 Train Validation split](#Split)\n* [5. Build & Train Generalized additive model (GAM) with pyGAM](#GAM)\n > * [Partial dependency plots](#dependency-plot)\n* [6. GAM Model Evaluation & Diagnostics](#GAM-Evaluation)\n > * [6.1 Evaluation Metrics](#Metrics)\n > * [6.2 Residual Plots](#Residual)\n > * [6.3 Normality Test for Residuals](#Normality)\n > * [6.4 Heteroscedasticity Analysis of Residuals: Breusch-Pagan, Goldfeld-Quandt and White Tests in Python](#Heteroscedasticity)\n* [7. Model Prediction & Kaggle Submission](#Prediction)\n* [8. Concluding remarks](#Conclusion)\n* [9. References](#References)\n\n## **1. Import the required libraries**","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"# Import data handling & numerical libraries\nimport pandas as pd\nimport numpy as np\nfrom copy import copy\nimport datetime\n\n# Import Data Visualization libraries\nimport seaborn as sb\nimport matplotlib.pyplot as plt\n\n#import libraries for muting unnecessary warnings if needed\nimport warnings\nwarnings.filterwarnings('ignore')","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:08.137007Z","iopub.execute_input":"2022-07-22T12:41:08.138071Z","iopub.status.idle":"2022-07-22T12:41:09.318867Z","shell.execute_reply.started":"2022-07-22T12:41:08.137947Z","shell.execute_reply":"2022-07-22T12:41:09.317813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **2. Dataset information**\n[ Back to Table of contents
](#Table)\n\nThe dataset collected from the monitor system mounted in a solar house corresponds to approximately 40 days of monitoring data. The Goal is to predict indoor temperature of a room (the Bedroom), in order to choose whether or not to activate the HVAC (Heating, Ventilation, and Air Conditioning) system. The data was sampled every minute, computing and uploading it smoothed with 15 minute means. The dataset includes dates, other sensor measurements, weather measurements and other information. It is a multivariate time-series dataset. It has been established that the power consumption attributed to HVAC accounts for 53.9 % of total consumption, and the energy required to maintain the temperature is less than that required to drop or raise it. As a result, a predictive model capable of predicting a room's indoor temperature (a short-term forecast of indoor temperature) would help in lowering overall energy consumption, by deciding whether or not to activate the HVAC system, at the appropriate time. Below displayed map provides an idea of locations of solar house sensors and actuators.","metadata":{}},{"cell_type":"code","source":"from IPython.display import Image\nurl = '../input/smart-homes-temperature-time-series-forecasting/Solar house sensors and actuators map.png'\nImage(url,width=700, height=700)","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:09.320819Z","iopub.execute_input":"2022-07-22T12:41:09.321343Z","iopub.status.idle":"2022-07-22T12:41:09.346757Z","shell.execute_reply.started":"2022-07-22T12:41:09.321270Z","shell.execute_reply":"2022-07-22T12:41:09.345704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Reading datasets\ndf_train=pd.read_csv(\"../input/smart-homes-temperature-time-series-forecasting/train.csv\")\ndf_test=pd.read_csv(\"../input/smart-homes-temperature-time-series-forecasting/test.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:09.347988Z","iopub.execute_input":"2022-07-22T12:41:09.348412Z","iopub.status.idle":"2022-07-22T12:41:09.395903Z","shell.execute_reply.started":"2022-07-22T12:41:09.348369Z","shell.execute_reply":"2022-07-22T12:41:09.394780Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train.info()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:09.400016Z","iopub.execute_input":"2022-07-22T12:41:09.400935Z","iopub.status.idle":"2022-07-22T12:41:09.434591Z","shell.execute_reply.started":"2022-07-22T12:41:09.400887Z","shell.execute_reply":"2022-07-22T12:41:09.433378Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train.describe()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:09.435893Z","iopub.execute_input":"2022-07-22T12:41:09.436373Z","iopub.status.idle":"2022-07-22T12:41:09.508883Z","shell.execute_reply.started":"2022-07-22T12:41:09.436344Z","shell.execute_reply":"2022-07-22T12:41:09.507702Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To make sure that the dataset does not contain any missing values and/or duplicates create a complete list of DaTime with 15 minutes interval from the starting to end point and check whether it matches with the index list of dataset.","metadata":{}},{"cell_type":"code","source":"# sort by dates\ndf_train.sort_index(inplace = True)\n\n#creating datetime list with boundaries of raw data series, hourly frequency\ndatelist = pd.date_range(datetime.datetime(2012,3,13,11,45,0), datetime.datetime(2012,4,11,6,30,0), freq='15min').tolist()\n\n#extracting raw data series indices\nidx_list = df_train.index.to_list()\n\n#checking for anomalies by comparing the two\nprint(idx_list == datelist)\n#searching for anomalies\nprint(\"\\n No. of elements in full list:\", len(datelist), \"\\n No. of indices:\", len(idx_list), \"\\n No. of elements in set of indices:\", len(set(idx_list)))","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:09.510253Z","iopub.execute_input":"2022-07-22T12:41:09.510571Z","iopub.status.idle":"2022-07-22T12:41:09.536714Z","shell.execute_reply.started":"2022-07-22T12:41:09.510544Z","shell.execute_reply":"2022-07-22T12:41:09.535157Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.express as px\nfig = px.line(df_train['Indoor_temperature_room'])\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:09.538081Z","iopub.execute_input":"2022-07-22T12:41:09.539360Z","iopub.status.idle":"2022-07-22T12:41:12.050275Z","shell.execute_reply.started":"2022-07-22T12:41:09.539297Z","shell.execute_reply":"2022-07-22T12:41:12.048581Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **3. Model Identification**\n[ Back to Table of contents
](#Table)\n\n### **3.1 Adversarial Validation**\nAdversarial Validation is to check whether both the training and the test datasets come from the independent and identically distributed (IID) samples. This is done by classifying and encoding them into a separate feature classes by giving labels such as 0 for the training data and 1 for the test data, mix them up, and then do a classification analysis to check if we are able to correctly re-identify them using a binary classifier. More on Adversaial Validation can be read from [here](https://www.kaggle.com/code/carlmcbrideellis/what-is-adversarial-validation/notebook)","metadata":{}},{"cell_type":"code","source":"import xgboost as xgb\nfrom xgboost import XGBClassifier\nfrom xgboost import plot_importance\nfrom xgboost import cv\n# drop the target column from the training data\ntrain = df_train.drop(['Id', 'Indoor_temperature_room', 'Date', 'Time'], axis=1)\ntest = df_test.drop(['Id', 'Date', 'Time'], axis = 1)\n# select only the numerical features\n# add the train/test labels\ntrain[\"AV_class\"] = 0\ntest[\"AV_class\"] = 1\n\n# make one big dataset\ntrain_test = pd.concat([train, test], axis=0, ignore_index=True)\n\n# shuffle\ntrain_test_shuffled = train_test.sample(frac=1)\n\n# create our DMatrix (the XGBoost data structure)\nX = train_test_shuffled.drop(['AV_class'], axis=1)\ny = train_test_shuffled['AV_class']\nXGBdata = xgb.DMatrix(data=X,label=y)\n\n# our XGBoost parameters\nparams = {\"objective\":\"binary:logistic\",\n \"eval_metric\":\"logloss\",\n 'learning_rate': 0.05,\n 'max_depth': 5, }\n\n# perform cross validation with XGBoost\ncross_val_results = cv(dtrain=XGBdata, params=params, \n nfold=5, metrics=\"auc\", \n num_boost_round=200,early_stopping_rounds=20,\n as_pandas=True)\n\n# print out the final result\nprint((cross_val_results[\"test-auc-mean\"]).tail(1))","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:12.051551Z","iopub.execute_input":"2022-07-22T12:41:12.052329Z","iopub.status.idle":"2022-07-22T12:41:24.490485Z","shell.execute_reply.started":"2022-07-22T12:41:12.052289Z","shell.execute_reply":"2022-07-22T12:41:24.489477Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"AUC of 0.999735 indicates that the classifier is able to perfectly distinguish between the original training and test data. This means that both training and testing data are clearly distinguishable from each other. Feature importance plot gives an idea of that particular feature that contributes the most to the separation of training and testing data. ","metadata":{}},{"cell_type":"code","source":"classifier = XGBClassifier(eval_metric='logloss',use_label_encoder=False)\nclassifier.fit(X, y)\nfig, ax = plt.subplots(figsize=(12,4),dpi=100)\nplt.title('Feature Importances', size = 18, color='purple', loc='center', backgroundcolor='lavender', pad='10.0')\nplot_importance(classifier, ax=ax, color='#087E8B', edgecolor= 'cyan', linewidth=3)\nplt.show();","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:24.492364Z","iopub.execute_input":"2022-07-22T12:41:24.492722Z","iopub.status.idle":"2022-07-22T12:41:25.449781Z","shell.execute_reply.started":"2022-07-22T12:41:24.492683Z","shell.execute_reply":"2022-07-22T12:41:25.448711Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Since all features are equally important in predicting indoor tempearture it is not possible to drop any of the above features. \n### **3.2 Normality tests & Distribution plots**\n[ Back to Table of contents
](#Table)\nA more detailed analysis using normality tests determine how likely a data sample is from a normally distributed population using p-values. The null hypothesis for each test is that “the sample came from a normally distributed population”.","metadata":{}},{"cell_type":"code","source":"from scipy import stats\nfeatures_list = test.columns.values.tolist()\nfor feature in features_list:\n statistic, pvalue = stats.kstest(train[feature], test[feature]) #Kolmogorov-Smirnov test \n print(\"p-value %.2f\" %pvalue, \"for the feature\",feature)","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:25.453765Z","iopub.execute_input":"2022-07-22T12:41:25.454116Z","iopub.status.idle":"2022-07-22T12:41:26.187727Z","shell.execute_reply.started":"2022-07-22T12:41:25.454085Z","shell.execute_reply":"2022-07-22T12:41:26.186527Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"All the above p-values for Kolmogorov-Smirnov normality tests returned an alpha value less than 0.05, which means that the null hypothesis is to be rejected with 95% confidence level and it is likely that the data points do not come from a normal distribution. These features have completely different distributions between the training and the test datasets. The below displayed plot gives a visual experience of difference of distribution between training and testing data sets. ","metadata":{}},{"cell_type":"code","source":"!pip install -q ptitprince\nimport ptitprince as pt\nnum_feats=[col for col in df_test.columns if df_test[col].dtypes != 'object' and col !='Day_of_the_week' and col != 'Id']\nfig=plt.figure(figsize=(30,80))\nfor i, col in enumerate(num_feats):\n plt.subplot(len(num_feats),2,1*i+1)\n pt.RainCloud(data=df_train, y=df_train[col], bw=0.1, cut=0, hue=['dodgerblue'], orient='h', label=\"train data\", palette=['dodgerblue'], alpha = .65)\n pt.RainCloud(data=df_test, y=df_test[col], bw=0.1, cut=0, orient='h',label=\"test data\", hue=['crimson'], palette=['crimson'], alpha = .65)\n legend = plt.legend()\n plt.title(f'{col} distribution')\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:26.189555Z","iopub.execute_input":"2022-07-22T12:41:26.189915Z","iopub.status.idle":"2022-07-22T12:41:45.936012Z","shell.execute_reply.started":"2022-07-22T12:41:26.189884Z","shell.execute_reply":"2022-07-22T12:41:45.934538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"When the training data is completely different from testing data, it becomes impossible to predit a test data target value from a train data feature and target values. In this case The target values of test data lie outside of the range of the target values it was originally trained on.","metadata":{}},{"cell_type":"code","source":"df_train.Date.unique(), df_test.Date.unique()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:45.938482Z","iopub.execute_input":"2022-07-22T12:41:45.939073Z","iopub.status.idle":"2022-07-22T12:41:45.949643Z","shell.execute_reply.started":"2022-07-22T12:41:45.939023Z","shell.execute_reply":"2022-07-22T12:41:45.948746Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As can be seen above the train data contains features captured on end of winter and start of spring season and test data contain features from end of spring and beginning of summer months. Training on this type of data will never return a prediction value outside of what the model has seen in training data while using ensemble tree based regression models. This can result in erroneous predictions on unseen data. Here comes the importance of linear regression models which can extrapolate the prediction values beyond the prediction intervals of training data. ","metadata":{}},{"cell_type":"markdown","source":"## **Feature Engineering & Data preprocessing **\n[ Back to Table of contents
](#Table)\n\n### **Feature Extraction**\n\n1. Having a feature giving special importance to seasons is required because of the significance in difference that may arise in data collected in different months. \n\n2. From Solar house sensors and actuators map displayed above, it is clearly visible that the room is exposed to sunlight from east, west and south facades. Hence the heating effect from sunlight would be the combined effect of all three. The distribution plots also confirms the fact that they have similar distribution pattern. Instead of taking them as separate, it would be better to take the average of these three.\n\n3. CO2 dining and CO2 room also shows a similarity between each other. Hence those features are also averaged and combined into one. \n\n3. Extracting a minute feature is also important, its importance will be demonstrated further in the analysis after generating this feature. ","metadata":{}},{"cell_type":"code","source":"# To extract seasons\ndef month2seasons(x):\n if x in [12, 1, 2, 3]:\n season = 0\n elif x in [6, 7, 8]:\n season = 1\n elif x in [4, 5]:\n season = 2\n elif x in [9, 10, 11]:\n season = 3\n return season\n\nfor df in (df_train, df_test):\n df['Date'] = pd.to_datetime(df['Date'], dayfirst=True)\n df['Date'] = df['Date'].dt.strftime('%Y-%m-%d') \n df['DateTime'] = df['Date'] + ' ' + df['Time']\n df['DateTime'] = pd.to_datetime(df['DateTime'])\n df['hour'] = df['DateTime'].apply(lambda x : x.hour)\n df['Month'] = df['DateTime'].apply(lambda x : x.month)\n df['Minutes']=df['DateTime'].apply(lambda x: x.hour *60 + x.minute).astype(float)\n df[\"CO2_avg\"] = (df[\"CO2_(dinning-room)\"] + df[\"CO2_room\"])/2 # To extract average of CO@ features\n # To extract average of Meteo Sun light features\n df['Meteo_Sun_light_AVG'] = (df['Meteo_Sun_light_in_west_facade'] + df['Meteo_Sun_light_in_east_facade'] + df['Meteo_Sun_light_in_south_facade'])/3\n df['Season'] = df['Month'].apply(month2seasons)\n df.drop(['Day_of_the_week', 'Date','Time', \"CO2_(dinning-room)\", \"CO2_room\", 'Month','Meteo_Sun_light_in_west_facade',\n 'Meteo_Sun_light_in_east_facade', 'Meteo_Sun_light_in_south_facade'], axis = 1, inplace = True)\ndf_train.shape, df_train.columns","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:45.950866Z","iopub.execute_input":"2022-07-22T12:41:45.951612Z","iopub.status.idle":"2022-07-22T12:41:46.100673Z","shell.execute_reply.started":"2022-07-22T12:41:45.951579Z","shell.execute_reply":"2022-07-22T12:41:46.099370Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The following steps are done to demonstrate the importance of minutes columns. For this, first the dataset is grouped to get mean temperature for each 15 minute interval.","metadata":{}},{"cell_type":"code","source":"minutes=df_train[['Indoor_temperature_room','Minutes']].groupby(['Minutes']).mean().reset_index()\nminutes.columns = ['Minutes', 'mean_temp']\n#minutes['Minutes']= minutes['Minutes'].astype(int)\nscaled_temp = minutes['mean_temp'] - minutes['mean_temp'].mean()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:46.102132Z","iopub.execute_input":"2022-07-22T12:41:46.102500Z","iopub.status.idle":"2022-07-22T12:41:46.116815Z","shell.execute_reply.started":"2022-07-22T12:41:46.102467Z","shell.execute_reply":"2022-07-22T12:41:46.115876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Plotting mean temperature versus minutes.","metadata":{}},{"cell_type":"code","source":"# Define a function to plot the entire dataframe to performs data visualization\ndef display_plot(x, y, fig_title):\n plt.figure(figsize = (20,7))\n plt.title(fig_title, loc='center', fontsize=20)\n sb.barplot(x = x, y = y, palette = 'cool') \n plt.tight_layout();\ndisplay_plot(minutes['Minutes'],minutes['mean_temp'], \"Mean indoor room temperature\")\ndisplay_plot(minutes['Minutes'],scaled_temp, \"Scaled Mean indoor room temperature\")\n\nfig,ax=plt.subplots(figsize = (12,5),dpi=100)\n#ax.plot(minutes['Minutes'],minutes['mean_temp'])\nimport matplotlib.patheffects as pe\nax.plot(minutes['Minutes'],minutes['mean_temp'], lw = 3, color='#087E8B', \n path_effects=[pe.SimpleLineShadow(shadow_color='cyan'), pe.Stroke(linewidth=5, foreground='cyan'),pe.Normal()])\nplt.title('Sine wave pattern recognized', loc = 'Center', fontsize=18, color='purple', style='normal', backgroundcolor='lavender', pad='10.0')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:46.118428Z","iopub.execute_input":"2022-07-22T12:41:46.118834Z","iopub.status.idle":"2022-07-22T12:41:49.555165Z","shell.execute_reply.started":"2022-07-22T12:41:46.118790Z","shell.execute_reply":"2022-07-22T12:41:49.553904Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Obviously a sinosoidal wave pattern with respect to daily period, , which is quite similar to the pattern generated from a sine function as in the following","metadata":{}},{"cell_type":"code","source":"a1=np.arange(0,1400)\na2=-np.sin(2*np.pi*a1/1400) # sin (2pift)\ndisplay_plot(a1,a2, \"Sinosoidal wave pattern\\n $\\mathcal{A}\\mathrm{sin}(2 \\omega t)$\")\nfig,ax=plt.subplots(figsize = (12,5),dpi=100)\nax.plot(a1,a2, lw = 3, color='#087E8B', \n path_effects=[pe.SimpleLineShadow(shadow_color='cyan'), pe.Stroke(linewidth=5, foreground='cyan'),pe.Normal()])\nmono_font = {'fontname':'monospace'}\nplt.title('Sine wave with period 1400', loc = 'Center', fontsize=18, **mono_font,\n color='purple', style='oblique', backgroundcolor='lavender', pad='10.0')\nplt.title(\"$\\mathcal{A}\\mathrm{sin}(2 \\pi f t)$\", fontsize=16, loc = 'right', pad='10.0')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:41:49.556576Z","iopub.execute_input":"2022-07-22T12:41:49.556951Z","iopub.status.idle":"2022-07-22T12:42:07.852860Z","shell.execute_reply.started":"2022-07-22T12:41:49.556910Z","shell.execute_reply":"2022-07-22T12:42:07.851588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **Train Validation split**\n[ Back to Table of contents
](#Table)","metadata":{}},{"cell_type":"code","source":"df_train1 = df_train.drop(['Id', 'DateTime','hour'], axis = 1)\nX = df_train1.drop('Indoor_temperature_room', axis = 1)\ny = df_train1['Indoor_temperature_room']\nfrom sklearn.model_selection import train_test_split\nX_train, X_val, y_train, y_val = train_test_split(X, y, shuffle = True, test_size = 0.2)","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:07.854698Z","iopub.execute_input":"2022-07-22T12:42:07.855169Z","iopub.status.idle":"2022-07-22T12:42:07.866528Z","shell.execute_reply.started":"2022-07-22T12:42:07.855124Z","shell.execute_reply":"2022-07-22T12:42:07.865175Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **5. Build & Train Generalized additive model (GAM) with pyGAM**\n[ Back to Table of contents
](#Table)\nGAM allows us to easily examine the partial relationships between the response variable and the predictors. As the name implies its addictive nature ensures that the partial impact of each variable does not depend on the others in the model. This additive nature is used to explore and interpret individual features by holding others at their mean.","metadata":{}},{"cell_type":"code","source":"!pip install pygam --quiet\nfrom pygam import LinearGAM\ngam = LinearGAM(n_splines=10).fit(X, y)\ngam.summary()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:07.868356Z","iopub.execute_input":"2022-07-22T12:42:07.868763Z","iopub.status.idle":"2022-07-22T12:42:20.413517Z","shell.execute_reply.started":"2022-07-22T12:42:07.868732Z","shell.execute_reply":"2022-07-22T12:42:20.411877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **Partial dependency plots**\n[ Back to Table of contents
](#Table)\n\nPartial dependence functions can be used to inspect the contribution of each feature and demonstrate partial relationships. ","metadata":{}},{"cell_type":"code","source":"## plotting\nplt.figure();\nfig, axes = plt.subplots(3,ncols =int(len(X.columns)/3), figsize = (30, 6));\n\ntitles = X.columns\nfor i, (col,ax) in enumerate(zip(X.columns, axes.flatten())):\n XX = gam.generate_X_grid(term=i)\n ax.plot(XX[:, i], gam.partial_dependence(term=i, X=XX))\n ax.plot(XX[:, i], gam.partial_dependence(term=i, X=XX, width=.95)[1], c='r', ls='--')\n #if i == 0:\n #ax.set_ylim(-30,30)\n ax.set_title(titles[i]);\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:20.416143Z","iopub.execute_input":"2022-07-22T12:42:20.417164Z","iopub.status.idle":"2022-07-22T12:42:21.958380Z","shell.execute_reply.started":"2022-07-22T12:42:20.417106Z","shell.execute_reply":"2022-07-22T12:42:21.957601Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **6. GAM Model Evaluation & Diagnostics**\n[ Back to Table of contents
](#Table)\n### **6.1 Evaluation Metrics**","metadata":{}},{"cell_type":"code","source":"def model_train_evaluation(y, ypred, model_name): \n \n # Model Evaluation metrics\n from sklearn.metrics import mean_squared_error,mean_absolute_error,explained_variance_score, r2_score, mean_absolute_percentage_error\n print(\"\\n \\n Model Evaluation Report: \")\n print('Mean Absolute Error(MAE) of', model_name,':', mean_absolute_error(y, ypred))\n print('Mean Squared Error(MSE) of', model_name,':', mean_squared_error(y, ypred))\n print('Root Mean Squared Error (RMSE) of', model_name,':', mean_squared_error(y, ypred, squared = False))\n print('Mean absolute percentage error (MAPE) of', model_name,':', mean_absolute_percentage_error(y, ypred))\n print('Explained Variance Score (EVS) of', model_name,':', explained_variance_score(y, ypred))\n print('R2 of', model_name,':', (r2_score(y, ypred)).round(2))\n print('\\n \\n')\n \n # Actual vs Predicted Plot\n f, ax = plt.subplots(figsize=(12,6),dpi=100);\n plt.scatter(y, ypred, label=\"Actual vs Predicted\")\n # Perfect predictions\n plt.xlabel('Indoor Temperature in celsius')\n plt.ylabel('Indoor Temperature in celsius')\n plt.title('Expection vs Prediction')\n plt.plot(y,y,'r', label=\"Perfect Expected Prediction\")\n plt.legend()\n f.text(0.95, 0.06, 'AUTHOR: RINI CHRISTY',\n fontsize=12, color='green',\n ha='left', va='bottom', alpha=0.5);\n \n print('\\n \\n \\n \\n')\n fig,ax=plt.subplots(figsize=(15,8))\n plt.plot(y.values, lw = 4, label='Actual values', color = 'blue')\n plt.plot(ypred, label='Predicted values', color = 'red')\n plt.legend(loc='best')\n plt.title(f'Actual vs Predicted for {model_name}')\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:21.959622Z","iopub.execute_input":"2022-07-22T12:42:21.960453Z","iopub.status.idle":"2022-07-22T12:42:21.974173Z","shell.execute_reply.started":"2022-07-22T12:42:21.960421Z","shell.execute_reply":"2022-07-22T12:42:21.972742Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### **6.1.1 Evaluation metrics for n_splines = 10**\n[ Back to Table of contents
](#Table)","metadata":{}},{"cell_type":"code","source":"yhat = gam.predict(X)\nmodel_train_evaluation(y, yhat, 'GAM with n_splines=10')","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:21.976487Z","iopub.execute_input":"2022-07-22T12:42:21.977015Z","iopub.status.idle":"2022-07-22T12:42:22.610248Z","shell.execute_reply.started":"2022-07-22T12:42:21.976967Z","shell.execute_reply":"2022-07-22T12:42:22.608862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### **6.1.2 Evaluation metrics for n_splines = 40**\n[ Back to Table of contents
](#Table)","metadata":{}},{"cell_type":"code","source":"gam = LinearGAM(n_splines=40).fit(X, y)\n## plotting\nplt.figure();\nfig, axes = plt.subplots(3,ncols =int(len(X.columns)/3), figsize = (30, 6));\n\ntitles = X.columns\nfor i, (col,ax) in enumerate(zip(X.columns, axes.flatten())):\n XX = gam.generate_X_grid(term=i)\n ax.plot(XX[:, i], gam.partial_dependence(term=i, X=XX))\n ax.plot(XX[:, i], gam.partial_dependence(term=i, X=XX, width=.95)[1], c='r', ls='--')\n #if i == 0:\n #ax.set_ylim(-30,30)\n ax.set_title(titles[i]);\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:22.611994Z","iopub.execute_input":"2022-07-22T12:42:22.612442Z","iopub.status.idle":"2022-07-22T12:42:25.311736Z","shell.execute_reply.started":"2022-07-22T12:42:22.612406Z","shell.execute_reply":"2022-07-22T12:42:25.310332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"yhat = gam.predict(X)\nmodel_train_evaluation(y, yhat, 'GAM (n_splines=40) prediction using training set ')","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:25.313684Z","iopub.execute_input":"2022-07-22T12:42:25.314142Z","iopub.status.idle":"2022-07-22T12:42:25.995619Z","shell.execute_reply.started":"2022-07-22T12:42:25.314095Z","shell.execute_reply":"2022-07-22T12:42:25.994481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ypred = gam.predict(X_val)\nmodel_train_evaluation(y_val, ypred, 'GAM prediction using validation set')","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:25.997139Z","iopub.execute_input":"2022-07-22T12:42:25.999406Z","iopub.status.idle":"2022-07-22T12:42:26.551891Z","shell.execute_reply.started":"2022-07-22T12:42:25.999354Z","shell.execute_reply":"2022-07-22T12:42:26.550585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gam.summary()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:26.553616Z","iopub.execute_input":"2022-07-22T12:42:26.554103Z","iopub.status.idle":"2022-07-22T12:42:26.566123Z","shell.execute_reply.started":"2022-07-22T12:42:26.554059Z","shell.execute_reply":"2022-07-22T12:42:26.564845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **6.2 Residual Plots**\n[ Back to Table of contents
](#Table)\nMany of the assumptions that are necessary to have a valid model can be checked by identifying patterns in the residuals of that model. We can make a quick visual check by looking at the residual plot of a given model.\nWith a residual plot, we look at the predicted values of the model versus the residuals themselves. What we expect to see is just a cloud of unrelated points****","metadata":{}},{"cell_type":"code","source":"residuals = y_val.values-ypred\nplt.scatter(ypred, residuals);\nplt.axhline(0, color='red')\nplt.xlabel('Predicted Values');\nplt.ylabel('Residuals');","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:26.567999Z","iopub.execute_input":"2022-07-22T12:42:26.568467Z","iopub.status.idle":"2022-07-22T12:42:26.777208Z","shell.execute_reply.started":"2022-07-22T12:42:26.568423Z","shell.execute_reply":"2022-07-22T12:42:26.776034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# line plot of residuals\nresiduals = pd.DataFrame(residuals)\nresiduals.plot()\n#plt.ylim(-0.5, 0.5)\nplt.show()\n# density plot of residuals\nresiduals.plot(kind='kde')\nplt.xlim(-10, 10)\nplt.show()\n# summary stats of residuals\nprint(residuals.describe())","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:26.778860Z","iopub.execute_input":"2022-07-22T12:42:26.779408Z","iopub.status.idle":"2022-07-22T12:42:27.198142Z","shell.execute_reply.started":"2022-07-22T12:42:26.779377Z","shell.execute_reply":"2022-07-22T12:42:27.197069Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **6.3 Normality Test for Residuals**\n[ Back to Table of contents
](#Table)\n\nThe above residual plots visually confirms that the model satisfies independent and identically distributed (i.i.d) data requirement for residuals. Mathematical confirmation of i.i.d requirement can be executed using normality tests like: Shapiro Wilk test, Kolmogorov-Smirnov test, Jarque-Bera test, Anderson-Darling test (from both scipy and statsmodels), Kolmogorov-Smirnov, and D’Agostino K-squared tests.\n\nThe normality tests determine how likely a data sample is from a normally distributed population using p-values. The null hypothesis for each test is that “the sample came from a normally distributed population”. This means that if the resulting p-values are below a chosen alpha value, then the null hypothesis is rejected. Thus there is evidence to suggest that the data comes from a non-normal distribution. For this project, use an Alpha value of 0.01.\nThe below executed normality tests determine how likely a data sample is from a normally distributed population using p-values. The null hypothesis for each test is that “the sample came from a normally distributed population”.","metadata":{}},{"cell_type":"code","source":"from scipy import stats\nfrom statsmodels.stats.diagnostic import normal_ad\nimport statsmodels\nsw_result = stats.shapiro(residuals)#Shapiro Wilk test\nKS_result = stats.kstest(residuals, 'norm')#Kolmogorov-Smirnov test \nJB_result = stats.jarque_bera(residuals) #Jarque-Bera test\n#AD_result = stats.anderson(residuals)#Anderson-Darling test\nad_result = normal_ad(np.array(residuals), axis=0)#Anderson-Darling test\ndag_result = stats.normaltest(residuals, axis=0, nan_policy='propagate')#D’Agostino’s K-squared test\nlf_result = statsmodels.stats.diagnostic.lilliefors(residuals)#Lilliefors test\n#print(f'\\n \\n Residulas for Indoor_temperature_room')\n#print(residuals)\nprint(f'\\n \\n Shapiro Wilk test results for Indoor_temperature_room')\nprint(sw_result)\nprint(f'\\n \\n Kolmogorov-Smirnov test results for Indoor_temperature_room')\nprint(KS_result)\nprint(f'\\n \\n Jarque-Bera test results for Indoor_temperature_room')\nprint(JB_result)\n#print(f'\\n \\n Anderson-Darling test results for Indoor_temperature_room')\n#print(AD_result)\nprint(f'\\n \\n Anderson-Darling test results for Indoor_temperature_room')\nprint(ad_result)\nprint(f'\\n \\n D’Agostino’s K-squared test results for Indoor_temperature_room')\nprint(dag_result)\nprint(f'\\n \\n Lilliefors test results for Indoor_temperature_room')\nprint(lf_result)","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:27.203888Z","iopub.execute_input":"2022-07-22T12:42:27.204237Z","iopub.status.idle":"2022-07-22T12:42:27.332037Z","shell.execute_reply.started":"2022-07-22T12:42:27.204206Z","shell.execute_reply":"2022-07-22T12:42:27.330721Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"p values of >0.01 indicates that the null hypothesis cannot be rejected and hence all residuals come from normally distributed sample population. \n\n### **6.4 Heteroscedasticity Analysis of Residuals: Breusch-Pagan, Goldfeld-Quandt and White Tests in Python**\n[ Back to Table of contents
](#Table)\n\nTest for heteroskedasticity of standardized residuals\n\nTests whether the sum-of-squares in the first third of the sample is significantly different than the sum-of-squares in the last third of the sample. Analogous to a Goldfeld-Quandt test. The null hypothesis is of no heteroskedasticity.\n\n#### **6.4.1 Breush-Pagan test**\nA Breusch-Pagan test uses the following null and alternative hypotheses:\n\nThe null hypothesis (H0): Homoscedasticity is present.\n\nThe alternative hypothesis: (Ha): Homoscedasticity is not present (i.e. heteroscedasticity exists).\n\nIf the p-value of the test results is smaller than alpha (significance level) then we can reject H0 and conclude that the data is heteroscedastic. We will use 0.05 as the signifcance level parameter.\n\nWe need to convert y_val into a 2d array since this is required for as one of the input in Breusch-Pagan test","metadata":{}},{"cell_type":"code","source":"def test_model(col):\n s = []\n for i in col:\n a = [1,i]\n s.append(a)\n return (np.array(s))\nHet = df_train.iloc[8:,]\n\nexog = test_model(y_val)\n# Method 1\nfrom statsmodels.stats.diagnostic import het_breuschpagan\nbreusch_pagan_test = het_breuschpagan(residuals, exog)\nprint(breusch_pagan_test)\nprint ('\\n Het_breuschpagan-test p_value:', breusch_pagan_test[1])\nif breusch_pagan_test[1] > 0.05:\n print(\"The residuals are not heteroscedastic.\")\nif breusch_pagan_test[1] < 0.05:\n print(\"The residuals are heteroscedastic.\")","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:27.333393Z","iopub.execute_input":"2022-07-22T12:42:27.334218Z","iopub.status.idle":"2022-07-22T12:42:27.348060Z","shell.execute_reply.started":"2022-07-22T12:42:27.334182Z","shell.execute_reply":"2022-07-22T12:42:27.346807Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Method 2\nimport statsmodels.stats.api as sms\nfrom statsmodels.compat import lzip\nimport statsmodels.tools.tools as smt\nname = [\"Lagrange multiplier statistic\", \"p-value\", \"f-value\", \"f p-value\"]\ntest = sms.het_breuschpagan(residuals, exog_het= exog )\nlzip(name, test)","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:27.349536Z","iopub.execute_input":"2022-07-22T12:42:27.350572Z","iopub.status.idle":"2022-07-22T12:42:27.407022Z","shell.execute_reply.started":"2022-07-22T12:42:27.350532Z","shell.execute_reply":"2022-07-22T12:42:27.406094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### **6.4.2 Goldfeld-Quandt test**\n[ Back to Table of contents
](#Table)","metadata":{}},{"cell_type":"code","source":"name = [\"F statistic\", \"p-value\"]\nGQ_test = sms.het_goldfeldquandt(residuals, exog)\nprint(lzip(name, test))\nif GQ_test[1] > 0.05:\n print(\"The residuals are not heteroscedastic.\")\nif GQ_test[1] < 0.05:\n print(\"The residuals are heteroscedastic.\")","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:27.408412Z","iopub.execute_input":"2022-07-22T12:42:27.408769Z","iopub.status.idle":"2022-07-22T12:42:27.419019Z","shell.execute_reply.started":"2022-07-22T12:42:27.408738Z","shell.execute_reply":"2022-07-22T12:42:27.417775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### **6.4.3 White test**\n[ Back to Table of contents
](#Table)\nWhite’s test uses the following null and alternative hypotheses:\n\nNull (H0): Homoscedasticity is present (residuals are equally scattered)\n\nAlternative (HA): Heteroscedasticity is present (residuals are not equally scattered)","metadata":{}},{"cell_type":"code","source":"from statsmodels.stats.diagnostic import het_white\n#define labels to use for output of White's test\nlabels = ['Test Statistic', 'Test Statistic p-value', 'F-Statistic', 'F-Test p-value']\nwhite_test = het_white(residuals.values,exog) \n#print results of White's test\nprint(dict(zip(labels, white_test)))\nif white_test[1] > 0.05:\n print(\"The residuals are not heteroscedastic.\")\nif white_test[1] < 0.05:\n print(\"The residuals are heteroscedastic.\")","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:27.420767Z","iopub.execute_input":"2022-07-22T12:42:27.421867Z","iopub.status.idle":"2022-07-22T12:42:27.432882Z","shell.execute_reply.started":"2022-07-22T12:42:27.421821Z","shell.execute_reply":"2022-07-22T12:42:27.431815Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The residuals satisfy the i.i.d (independent and identically distributed) requirements and hence the model seems quite alright. Let’s proceed with the prediction.\n\n## **7. Model Prediction & Kaggle Submission**\n[ Back to Table of contents
](#Table)","metadata":{}},{"cell_type":"code","source":"X_test = df_test.drop(['Id','DateTime','hour'], axis = 1)\nX_test.shape, X_test.columns, X_train.shape, X_train.columns ","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:27.434189Z","iopub.execute_input":"2022-07-22T12:42:27.434692Z","iopub.status.idle":"2022-07-22T12:42:27.445132Z","shell.execute_reply.started":"2022-07-22T12:42:27.434626Z","shell.execute_reply":"2022-07-22T12:42:27.444268Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test['forecast'] = gam.predict(X_test)\n#test['forecast'] = forecast\nFinal = df_test[['Id', 'forecast']]\nFinal.columns = ['Id', 'Indoor_temperature_room']\nFinal.to_csv('submission.csv', index = False)\nFinal","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:27.446521Z","iopub.execute_input":"2022-07-22T12:42:27.446883Z","iopub.status.idle":"2022-07-22T12:42:27.520254Z","shell.execute_reply.started":"2022-07-22T12:42:27.446851Z","shell.execute_reply":"2022-07-22T12:42:27.519093Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize= (12,8))\nplt.plot(Final['Indoor_temperature_room']);","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:27.521511Z","iopub.execute_input":"2022-07-22T12:42:27.521829Z","iopub.status.idle":"2022-07-22T12:42:27.738724Z","shell.execute_reply.started":"2022-07-22T12:42:27.521802Z","shell.execute_reply":"2022-07-22T12:42:27.737889Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **8. Concluding remarks**\n[ Back to Table of contents
](#Table)\nThe training dataset and testing dataset come from different distributions due to data collection done in different months and since there exits a possibility that test prediction values lies outside the range of trained values, a linear regression seemed a right model for this data. The GAM model used seems a perfect fit with residues satisfying a i.i.d requirements and a mean error in the range of 0.7. \n\n## **9. References:**\n[ Back to Table of contents
](#Table)\n\n1. [What is Adversarial Validation?](https://www.kaggle.com/code/carlmcbrideellis/what-is-adversarial-validation/notebook)\n2. [Multivariate time series forecasting: Linear-tree](https://www.kaggle.com/code/carlmcbrideellis/multivariate-time-series-forecasting-linear-tree)\n3. [Building interpretable models with Generalized additive models in Python](https://medium.com/just-another-data-scientist/building-interpretable-models-with-generalized-additive-models-in-python-c4404eaf5515)\n4. [pyGAM : Getting Started with Generalized Additive Models in Python](https://codeburst.io/pygam-getting-started-with-generalized-additive-models-in-python-457df5b4705f)","metadata":{}},{"cell_type":"code","source":"%%html\n","metadata":{"execution":{"iopub.status.busy":"2022-07-22T12:42:27.739773Z","iopub.execute_input":"2022-07-22T12:42:27.740664Z","iopub.status.idle":"2022-07-22T12:42:27.747139Z","shell.execute_reply.started":"2022-07-22T12:42:27.740621Z","shell.execute_reply":"2022-07-22T12:42:27.746372Z"},"trusted":true},"execution_count":null,"outputs":[]}]}