<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/doing-magic-and-analyzing-seasonal-time-series-with-gam-generalized-additive-model-in-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://petolau.github.io/Analyzing-double-seasonal-time-series-with-GAM-in-R/"> Peter Laurinec</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>As I wrote in the previous post, I will continue in describing <strong>regression</strong> methods, which are suitable for double seasonal (or multi-seasonal) <strong>time series</strong>. In the previous post about <a target="_blank" rel="nofollow" href="https://petolau.github.io/Forecast-double-seasonal-time-series-with-multiple-linear-regression-in-R/">Multiple Linear Regression</a>, I showed how to use “simple” <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Ordinary_least_squares"><strong>OLS</strong></a> regression method to model <strong>double seasonal time series</strong> of electricity consumption and use it for accurate forecasting. Interactions between two seasonal variables were successfully used to achieve this goal. The issue of <a target="_blank" rel="nofollow" href="https://petolau.github.io/Forecast-electricity-consumption-with-similar-day-approach-in-R/">forecasting time series from smart meters</a> was discussed in my first post.</p>
<p>In this post (tutorial), I will fully introduce a “magical” <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Generalized_additive_model"><strong>Generalized Additive Model</strong></a> (GAM) to model time series of electricity consumption. You may wonder, why words “<em>magic</em>” or “<em>magical</em>” are mentioned? You will get your answer later in this post <img alt="😉" class="wp-smiley" src="https://s.w.org/images/core/emoji/2/72x72/1f609.png"></p>
<p>Now I have to warn you that this blog post is very long. <strong>GAM</strong> is a very complex method and I wanted to go little bit deeper in order to show you this interesting regression method. I think it will be worth it and in the end of the post, I hope you will understand the theory behind <strong>R</strong> functions. In my opinion, <strong>R</strong> blogging (or any data mining or machine learning blogging) without theory interpretation of used methods has zero informative value for a reader.</p>
<p>So, let’s look under the roof of the <strong>GAM</strong> regression method. I think you need some motivation to it – here it is!</p>
<p><img alt="" src="https://i1.wp.com/petolau.github.io/images/post_3/motivate.png?w=456" data-recalc-dims="1"></p>
<h3>Theory behind Generalized Additive Model (GAM)</h3>
<p>What do these three words (or letters) in the name of this method mean and where does it come from? I am sure that you know something about <strong>Linear Model</strong> (maybe because you had read my <a target="_blank" rel="nofollow" href="https://petolau.github.io/Forecast-double-seasonal-time-series-with-multiple-linear-regression-in-R/">previous post about MLR</a> :relaxed:). It is a simple regression method which models the response (dependent) variable by independent variable(s). It’s solved by the OLS method. So now we know, what the <strong>M</strong> in the name means. When we want to linearly model a response variable which is not from <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Normal_distribution">normal Gaussian distribution</a>, for example, it can be binary (logistic regression) or discrete (Poisson) variable, we can use a generalization of the linear model – <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Generalized_linear_model">Generalized Linear Model</a> (GLM). It’s solved by the <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares">iteratively reweighted least squares</a> (IRLS) method. The <strong>G</strong> in the name is also clear. The <strong>A</strong> letter <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Additive_model">Additive Model</a>, means that the response variable depends linearly on unknown smooth functions. In the other words, the goal is to model the response variable by independent variables, which are in the form of some smooth functions. Voilà, <strong>GAM</strong> is created.</p>
<p>The GAM can be formally written as:</p>
<p></p>
<p></p>
<p>where \( i = 1, \dots, N \), \( g \) is a link function (identical, logarithmic or inverse), \( y \) is a response variable, \( x_1, \dots, x_p \) are independent variables, \( \beta_0 \) is an intercept, \( f_1, \dots, f_p \) are unknown smooth functions and \( \varepsilon \) is an i.i.d. random error.</p>
<p>The smooth function \( f \) is composed by sum of basis functions \( b \) and their corresponding regression coefficients \( \beta \), formally written:</p>
<p></p>
<p></p>
<p>Smooth functions are also called <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Spline_(mathematics)">splines</a>. <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Smoothing_spline">Smoothing splines</a> are real functions that are piecewise-defined by polynomial functions (basis functions). The places, where the polynomial pieces connect are called knots. In <strong>GAMs</strong>, penalized regression splines are used in order to regularize the smoothness of a spline.</p>
<p>Therefore, the model can be written in a linear way like this:</p>
<p></p>
<p>where \( \mathbf{X} \) is a model matrix and \( \beta \) is a vector of regression coefficients.</p>
<p>Then, the objective function to be minimized is:</p>
<p></p>
<p>where \( \lambda \) is a smoothing parameter and the integral of squares of second derivatives can be written as:</p>
<p></p>
<p>where \( \mathbf{S} \) is matrix of known coefficients.</p>
<p>All this mathematical madness then implies that regression coefficients can be obtained (estimated) by the equation:</p>
<p></p>
<p>\( \hat{\beta} \) is called penalized least squares estimator in this case.<br>
The method of obtaining the estimate of the \( \beta \) is called Penalized Iteratively Re-weighted Least Squares (P-IRLS).</p>
<p>There are some more important theoretical questions to answer. For example, what kind of smoothing splines exists? Or, how to choose an optimal smoothing parameter \( \lambda \)? How basis dimensions (or a number of knots) are set?</p>
<p>There are several smoothing bases \( b \) (splines) which are suitable for regression:</p>
<ul>
<li>thin plate regression splines</li>
<li>cubic regression spline</li>
<li>cyclic cubic regression spline</li>
<li>P-splines</li>
</ul>
<p>An example of the cubic basis function for dimension \( q = 3 \):</p>
<p></p>
<p>You can read more about them in references, which I will write up. In this post, for daily seasonality cubic regression spline, and for weekly seasonality P-splines will be used. Both types of splines are knot-based, so choosing a right number of knots will be important.</p>
<p>Next, an important procedure is to choose (estimate) an optimal smoothing parameter \( \lambda \) and the number of basis dimensions (i.e. degrees of freedom). This can be done in <strong>GAM</strong> by Generalized Cross Validation score (GCV). It minimizes an equation:</p>
<p></p>
<p>where \( \mathbf{A} \) is the <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Projection_matrix">projection matrix</a> (i.e. influence or hat). It’s obvious that when lambda is near 1 then spline will be over-smoothed, in opposite side when lambda is near zero than spline isn’t penalized, so the method behaves like a classical OLS. With a number of basis dimensions (estimated degrees of freedom), it is opposite. Higher dimension implies that fit will be less smoothed (overfit), on the other side lower dimensions implies more smoothed behavior of fitted values.</p>
<h4>Interactions</h4>
<p>We are almost at the end of the explanation of <strong>GAMs</strong> theory. One more thing. As I showed in the previous blog post, <strong>interactions</strong> are a very important part of the <strong>regression</strong> model for <strong>double seasonal time series</strong>. With <strong>GAMs</strong> there are four (!) main possibilities, how to include them to the model. First is the most basic, like in <strong>MLR</strong>, the multiplication of two independent variables: \( x_1\times x_2 \). Second one is possibility to use smoothed function to one variable: \( f_1(x_1)\times x_2 \). Third one comes to use same smoothed function for both variables: \( f_1(x_1)\times f_1(x_2) \mbox{ (often denoted } f(x_1, x_2)) \). Fourth one is the most complex, with GAM it is possible to use <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Tensor_product">tensor product</a> interactions. So it is possible to use different smoothing bases for variables and penalize it in two (when we do interactions of two independent variables) different ways: \( f_1(x_1)\otimes f_2(x_2) \). More nicely, tensor product interactions can be written as:</p>
<p></p>
<p>where \( b_1 \) and \( b_2 \) are basis functions, \( I \) and \( J \) are corresponding basis dimensions and \( \delta \) is vector of unknown coefficients.</p>
<p>This allows for an overall <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Anisotropy">anisotropic</a> (different in each direction) penalty, so the overall shape of a tensor product smooth is invariant to a rescaling of its independent variables. This is a huge advantage in the comparison to usage of one smoothing function. Simply said, we have theoretically supported that it’s allowed to use different metrics of variables in the interactions term.</p>
<h4>Resources (references)</h4>
<p>Phew…that’s almost all of the difficult theory behind <strong>GAM</strong> method. I will add some more things in analytical part of this post. Honestly, it was difficult for me to find this information on the web. So if you want to read more about <strong>GAM</strong>, here is the list of useful links which were used in my post:</p>
<ul>
<li>Amazing book: Simon N. Wood: <a target="_blank" rel="nofollow" href="https://www.crcpress.com/Generalized-Additive-Models-An-Introduction-with-R/Wood/p/book/9781584884743">Generalized Additive Models: an introduction with R</a>, CRC Press, 2006.</li>
<li>Great blog post by Gavin Simpson: <a target="_blank" rel="nofollow" href="http://www.fromthebottomoftheheap.net/2014/05/09/modelling-seasonal-data-with-gam/">Modelling seasonal data with GAMs</a></li>
<li>Very helpful StackExchange (Cross Validated) <a target="_blank" rel="nofollow" href="http://stats.stackexchange.com/questions/45446/intuition-behind-tensor-product-interactions-in-gams-mgcv-package-in-r">question about the intuition behind tensor product interactions</a></li>
<li>This <a target="_blank" rel="nofollow" href="http://www.sfs.uni-tuebingen.de/~jvanrij/Tutorial/GAMM.html">tutorial about mixed GAMs</a> can be helpful too.</li>
</ul>
<p>Go ahead to modeling and analyzing time series with <strong>GAMs</strong>.</p>
<h3>Doing “magic” with GAMs for modeling time series</h3>
<p>I have prepared a file with four aggregated time series of electricity consumption for an analysis. It can be found on <a target="_blank" rel="nofollow" href="https://github.com/PetoLau/petolau.github.io/tree/master/_rmd">my GitHub repo</a>, the name of the file is <em>DT_4_ind</em>. The file was created easily by the package <code>feather</code> (<a target="_blank" rel="nofollow" href="https://cran.r-project.org/package=feather">CRAN link</a>). Data manipulations will be done then by <code>data.table</code> package.</p>
<p>GAM methods are implemented in <strong>R</strong> in the awesome package <code>mgcv</code> by Simon Wood (<a target="_blank" rel="nofollow" href="https://cran.r-project.org/package=mgcv">link to <code>mgcv</code></a>).</p>
<p>For visualizations packages <code>ggplot2</code>, <code>grid</code> and <code>animation</code> will be used. One useful function from package <code>car</code> will be used too.</p>
<p>Let’s scan all of the needed packages.</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>feather</span><span>)</span><span>
</span><span>library</span><span>(</span><span>data.table</span><span>)</span><span>
</span><span>library</span><span>(</span><span>mgcv</span><span>)</span><span>
</span><span>library</span><span>(</span><span>car</span><span>)</span><span>
</span><span>library</span><span>(</span><span>ggplot2</span><span>)</span><span>
</span><span>library</span><span>(</span><span>grid</span><span>)</span><span>
</span><span>library</span><span>(</span><span>animation</span><span>)</span></code></pre>
</figure>
<p>Read the mentioned smart meter data by <code>read_feather</code> to one <code>data.table</code>.</p>
<figure>
<pre><code data-lang="r"><span>DT</span><span> </span><span>&lt;-</span><span> </span><span>as.data.table</span><span>(</span><span>read_feather</span><span>(</span><span>&quot;DT_4_ind&quot;</span><span>))</span></code></pre>
</figure>
<p>Prepare <code>DT</code> to work with a <strong>GAM</strong> regression model. Transform the characters of weekdays to integers and use function <code>recode</code> from package <code>car</code> to recode weekdays as there are coming in the week: 1. Monday, …, 7. Sunday.</p>
<figure>
<pre><code data-lang="r"><span>DT</span><span>[,</span><span> </span><span>week_num</span><span> </span><span>:=</span><span> </span><span>as.integer</span><span>(</span><span>car</span><span>::</span><span>recode</span><span>(</span><span>week</span><span>,</span><span>
</span><span>&quot;'Monday'='1';'Tuesday'='2';'Wednesday'='3';'Thursday'='4';
'Friday'='5';'Saturday'='6';'Sunday'='7'&quot;</span><span>))]</span></code></pre>
</figure>
<p>Store informations in variables of the type of industry, date, weekday and period for simpler working.</p>
<figure>
<pre><code data-lang="r"><span>n_type</span><span> </span><span>&lt;-</span><span> </span><span>unique</span><span>(</span><span>DT</span><span>[,</span><span> </span><span>type</span><span>])</span><span>
</span><span>n_date</span><span> </span><span>&lt;-</span><span> </span><span>unique</span><span>(</span><span>DT</span><span>[,</span><span> </span><span>date</span><span>])</span><span>
</span><span>n_weekdays</span><span> </span><span>&lt;-</span><span> </span><span>unique</span><span>(</span><span>DT</span><span>[,</span><span> </span><span>week</span><span>])</span><span>
</span><span>period</span><span> </span><span>&lt;-</span><span> </span><span>48</span></code></pre>
</figure>
<p>Let’s look at some data chunk of electricity consumption and do analysis on it. I have picked aggregate consumption of commercial properties for two weeks. Store it in variable <code>data_r</code> and plot it.</p>
<figure>
<pre><code data-lang="r"><span>data_r</span><span> </span><span>&lt;-</span><span> </span><span>DT</span><span>[(</span><span>type</span><span> </span><span>==</span><span> </span><span>n_type</span><span>[</span><span>1</span><span>]</span><span> </span><span>&amp;</span><span> </span><span>date</span><span> </span><span>%in%</span><span> </span><span>n_date</span><span>[</span><span>57</span><span>:</span><span>70</span><span>])]</span><span>
</span><span>ggplot</span><span>(</span><span>data_r</span><span>,</span><span> </span><span>aes</span><span>(</span><span>date_time</span><span>,</span><span> </span><span>value</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>()</span><span> </span><span>+</span><span>
</span><span>theme</span><span>(</span><span>panel.border</span><span> </span><span>=</span><span> </span><span>element_blank</span><span>(),</span><span>
</span><span>panel.background</span><span> </span><span>=</span><span> </span><span>element_blank</span><span>(),</span><span>
</span><span>panel.grid.minor</span><span> </span><span>=</span><span> </span><span>element_line</span><span>(</span><span>colour</span><span> </span><span>=</span><span> </span><span>&quot;grey90&quot;</span><span>),</span><span>
</span><span>panel.grid.major</span><span> </span><span>=</span><span> </span><span>element_line</span><span>(</span><span>colour</span><span> </span><span>=</span><span> </span><span>&quot;grey90&quot;</span><span>),</span><span>
</span><span>panel.grid.major.x</span><span> </span><span>=</span><span> </span><span>element_line</span><span>(</span><span>colour</span><span> </span><span>=</span><span> </span><span>&quot;grey90&quot;</span><span>),</span><span>
</span><span>axis.text</span><span> </span><span>=</span><span> </span><span>element_text</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>10</span><span>),</span><span>
</span><span>axis.title</span><span> </span><span>=</span><span> </span><span>element_text</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>12</span><span>,</span><span> </span><span>face</span><span> </span><span>=</span><span> </span><span>&quot;bold&quot;</span><span>))</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;Date&quot;</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Load (kW)&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-6" src="https://i2.wp.com/petolau.github.io/images/post_3/unnamed-chunk-6-1.png?w=456" data-recalc-dims="1"></p>
<p>There is possible to see two main seasonalities in plotted time series: daily and weekly. We have 48 measurements during the day and 7 days during the week so that will be our independent variables to model response variable – electricity load. Let’s construct it:</p>
<figure>
<pre><code data-lang="r"><span>N</span><span> </span><span>&lt;-</span><span> </span><span>nrow</span><span>(</span><span>data_r</span><span>)</span><span> </span><span># number of observations in the train set
</span><span>window</span><span> </span><span>&lt;-</span><span> </span><span>N</span><span> </span><span>/</span><span> </span><span>period</span><span> </span><span># number of days in the train set
</span><span>matrix_gam</span><span> </span><span>&lt;-</span><span> </span><span>data.table</span><span>(</span><span>Load</span><span> </span><span>=</span><span> </span><span>data_r</span><span>[,</span><span> </span><span>value</span><span>],</span><span>
</span><span>Daily</span><span> </span><span>=</span><span> </span><span>rep</span><span>(</span><span>1</span><span>:</span><span>period</span><span>,</span><span> </span><span>window</span><span>),</span><span>
</span><span>Weekly</span><span> </span><span>=</span><span> </span><span>data_r</span><span>[,</span><span> </span><span>week_num</span><span>])</span></code></pre>
</figure>
<p>Here we are! Train our first <strong>GAM</strong> with function <code>gam</code>. Independent variables are modeled by smoothing function <code>s</code>, for daily seasonality cubic regression spline is used, for weekly seasonality, P-splines is used, a number of knots are logically set to the number of unique values. Let’s do it.</p>
<figure>
<pre><code data-lang="r"><span>gam_1</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>Load</span><span> </span><span>~</span><span> </span><span>s</span><span>(</span><span>Daily</span><span>,</span><span> </span><span>bs</span><span> </span><span>=</span><span> </span><span>&quot;cr&quot;</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>period</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>Weekly</span><span>,</span><span> </span><span>bs</span><span> </span><span>=</span><span> </span><span>&quot;ps&quot;</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>7</span><span>),</span><span>
</span><span>data</span><span> </span><span>=</span><span> </span><span>matrix_gam</span><span>,</span><span>
</span><span>family</span><span> </span><span>=</span><span> </span><span>gaussian</span><span>)</span></code></pre>
</figure>
<p>Package <code>mgcv</code> have many advantages and nice features. First is its visualization capabilities. Let’s try it:</p>
<figure>
<pre><code data-lang="r"><span>layout</span><span>(</span><span>matrix</span><span>(</span><span>1</span><span>:</span><span>2</span><span>,</span><span> </span><span>nrow</span><span> </span><span>=</span><span> </span><span>1</span><span>))</span><span>
</span><span>plot</span><span>(</span><span>gam_1</span><span>,</span><span> </span><span>shade</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-9" src="https://i2.wp.com/petolau.github.io/images/post_3/unnamed-chunk-9-1.png?w=456" data-recalc-dims="1"></p>
<p>That looks nice, right? We can see here the influence of variables to electricity load. In the left plot, the peak of the load is around 3 p.m. during the day. In the right plot, we can see that during weekends the consumption logically decreases.</p>
<p>Let’s use <code>summary</code> function to do the diagnostic of our first model.</p>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>gam_1</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Load ~ s(Daily, bs = &quot;cr&quot;, k = period) + s(Weekly, bs = &quot;ps&quot;, 
##     k = 7)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2731.67      18.88   144.7   &lt;2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df     F p-value    
## s(Daily)  10.159 12.688 119.8  &lt;2e-16 ***
## s(Weekly)  5.311  5.758 130.3  &lt;2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.772   Deviance explained = 77.7%
## GCV = 2.4554e+05  Scale est. = 2.3953e+05  n = 672</code></pre>
</figure>
<p>What to look at here? EDF: estimated degrees of freedom – can be interpreted like how much given variable is smoothed (higher EDF value implies more complex splines). P-values: statistical significance of given variable to response variable, tested by F-test (lower is better). \( R^2 \) – adjusted R-squared (higher is better). GCV: GCV score (mentioned above). In the summary, we can see that the R-sq. (adj) value is a little bit low…</p>
<p>Let’s plot fitted values:</p>
<figure>
<pre><code data-lang="r"><span>datas</span><span> </span><span>&lt;-</span><span> </span><span>rbindlist</span><span>(</span><span>list</span><span>(</span><span>data_r</span><span>[,</span><span> </span><span>.</span><span>(</span><span>value</span><span>,</span><span> </span><span>date_time</span><span>)],</span><span>
</span><span>data.table</span><span>(</span><span>value</span><span> </span><span>=</span><span> </span><span>gam_1</span><span>$</span><span>fitted.values</span><span>,</span><span>
</span><span>data_time</span><span> </span><span>=</span><span> </span><span>data_r</span><span>[,</span><span> </span><span>date_time</span><span>])))</span><span>
</span><span>datas</span><span>[,</span><span> </span><span>type</span><span> </span><span>:=</span><span> </span><span>c</span><span>(</span><span>rep</span><span>(</span><span>&quot;Real&quot;</span><span>,</span><span> </span><span>nrow</span><span>(</span><span>data_r</span><span>)),</span><span> </span><span>rep</span><span>(</span><span>&quot;Fitted&quot;</span><span>,</span><span> </span><span>nrow</span><span>(</span><span>data_r</span><span>)))]</span><span>
</span><span>ggplot</span><span>(</span><span>data</span><span> </span><span>=</span><span> </span><span>datas</span><span>,</span><span> </span><span>aes</span><span>(</span><span>date_time</span><span>,</span><span> </span><span>value</span><span>,</span><span> </span><span>group</span><span> </span><span>=</span><span> </span><span>type</span><span>,</span><span> </span><span>colour</span><span> </span><span>=</span><span> </span><span>type</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>0.8</span><span>)</span><span> </span><span>+</span><span>
</span><span>theme_bw</span><span>()</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;Time&quot;</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Load (kW)&quot;</span><span>,</span><span>
</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;Fit from GAM n.1&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-11" src="https://i0.wp.com/petolau.github.io/images/post_3/unnamed-chunk-11-1.png?w=456" data-recalc-dims="1"></p>
<p>That is, of course, a horrible result. We need to include interactions of two independent variables to the model.</p>
<p>One more thing before that… What is hiding behind the stored <code>gam_1</code> object and used optimizer?</p>
<figure>
<pre><code data-lang="r"><span>gam_1</span><span>$</span><span>optimizer</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] &quot;magic&quot;</code></pre>
</figure>
<p>Here is that “magic”! <img alt="🙂" class="wp-smiley" src="https://s.w.org/images/core/emoji/2/72x72/1f642.png"></p>
<p>The first type of interactions will be the third one mentioned in the section about interactions, so one smoothing function to both variables is used. In the term <code>s(Daily, Weekly)</code> simple smoothing terms <code>s(Daily)</code> and <code>s(Weekly)</code> are also automatically included. Let’s do it.</p>
<figure>
<pre><code data-lang="r"><span>gam_2</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>Load</span><span> </span><span>~</span><span> </span><span>s</span><span>(</span><span>Daily</span><span>,</span><span> </span><span>Weekly</span><span>),</span><span>
</span><span>data</span><span> </span><span>=</span><span> </span><span>matrix_gam</span><span>,</span><span>
</span><span>family</span><span> </span><span>=</span><span> </span><span>gaussian</span><span>)</span><span>
</span><span>summary</span><span>(</span><span>gam_2</span><span>)</span><span>$</span><span>r.sq</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] 0.9352108</code></pre>
</figure>
<p>R-squared value suggests that result is much better. Look at the smooth term:</p>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>gam_2</span><span>)</span><span>$</span><span>s.table</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##                     edf   Ref.df        F p-value
## s(Daily,Weekly) 28.7008 28.99423 334.4754       0</code></pre>
</figure>
<p>Seems good too, the p-value is 0, which means that independent variable is significant. Plot of fitted values:</p>
<figure>
<pre><code data-lang="r"><span>datas</span><span> </span><span>&lt;-</span><span> </span><span>rbindlist</span><span>(</span><span>list</span><span>(</span><span>data_r</span><span>[,</span><span> </span><span>.</span><span>(</span><span>value</span><span>,</span><span> </span><span>date_time</span><span>)],</span><span>
</span><span>data.table</span><span>(</span><span>value</span><span> </span><span>=</span><span> </span><span>gam_2</span><span>$</span><span>fitted.values</span><span>,</span><span>
</span><span>data_time</span><span> </span><span>=</span><span> </span><span>data_r</span><span>[,</span><span> </span><span>date_time</span><span>])))</span><span>
</span><span>datas</span><span>[,</span><span> </span><span>type</span><span> </span><span>:=</span><span> </span><span>c</span><span>(</span><span>rep</span><span>(</span><span>&quot;Real&quot;</span><span>,</span><span> </span><span>nrow</span><span>(</span><span>data_r</span><span>)),</span><span> </span><span>rep</span><span>(</span><span>&quot;Fitted&quot;</span><span>,</span><span> </span><span>nrow</span><span>(</span><span>data_r</span><span>)))]</span><span>
</span><span>ggplot</span><span>(</span><span>data</span><span> </span><span>=</span><span> </span><span>datas</span><span>,</span><span> </span><span>aes</span><span>(</span><span>date_time</span><span>,</span><span> </span><span>value</span><span>,</span><span> </span><span>group</span><span> </span><span>=</span><span> </span><span>type</span><span>,</span><span> </span><span>colour</span><span> </span><span>=</span><span> </span><span>type</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>0.8</span><span>)</span><span> </span><span>+</span><span>
</span><span>theme_bw</span><span>()</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;Time&quot;</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Load (kW)&quot;</span><span>,</span><span>
</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;Fit from GAM n.2&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-15" src="https://i0.wp.com/petolau.github.io/images/post_3/unnamed-chunk-15-1.png?w=456" data-recalc-dims="1"></p>
<p>It’s not bad, but it can be better. Now, let’s try the above mentioned tensor product interactions. That can be done by function <code>te</code>, basis functions can be defined as well.</p>
<figure>
<pre><code data-lang="r"><span>gam_3</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>Load</span><span> </span><span>~</span><span> </span><span>te</span><span>(</span><span>Daily</span><span>,</span><span> </span><span>Weekly</span><span>,</span><span>
</span><span>bs</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;cr&quot;</span><span>,</span><span> </span><span>&quot;ps&quot;</span><span>)),</span><span>
</span><span>data</span><span> </span><span>=</span><span> </span><span>matrix_gam</span><span>,</span><span>
</span><span>family</span><span> </span><span>=</span><span> </span><span>gaussian</span><span>)</span><span>
</span><span>summary</span><span>(</span><span>gam_3</span><span>)</span><span>$</span><span>r.sq</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] 0.9268452</code></pre>
</figure>
<p>Similar to the previous model <code>gam_2</code>. Look at the smooth term:</p>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>gam_3</span><span>)</span><span>$</span><span>s.table</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##                       edf   Ref.df        F p-value
## te(Daily,Weekly) 23.65709 23.98741 354.5856       0</code></pre>
</figure>
<p>Again very similar results. Let’s look at fitted values:</p>
<figure>
<pre><code data-lang="r"><span>datas</span><span> </span><span>&lt;-</span><span> </span><span>rbindlist</span><span>(</span><span>list</span><span>(</span><span>data_r</span><span>[,</span><span> </span><span>.</span><span>(</span><span>value</span><span>,</span><span> </span><span>date_time</span><span>)],</span><span>
</span><span>data.table</span><span>(</span><span>value</span><span> </span><span>=</span><span> </span><span>gam_3</span><span>$</span><span>fitted.values</span><span>,</span><span>
</span><span>data_time</span><span> </span><span>=</span><span> </span><span>data_r</span><span>[,</span><span> </span><span>date_time</span><span>])))</span><span>
</span><span>datas</span><span>[,</span><span> </span><span>type</span><span> </span><span>:=</span><span> </span><span>c</span><span>(</span><span>rep</span><span>(</span><span>&quot;Real&quot;</span><span>,</span><span> </span><span>nrow</span><span>(</span><span>data_r</span><span>)),</span><span> </span><span>rep</span><span>(</span><span>&quot;Fitted&quot;</span><span>,</span><span> </span><span>nrow</span><span>(</span><span>data_r</span><span>)))]</span><span>
</span><span>ggplot</span><span>(</span><span>data</span><span> </span><span>=</span><span> </span><span>datas</span><span>,</span><span> </span><span>aes</span><span>(</span><span>date_time</span><span>,</span><span> </span><span>value</span><span>,</span><span> </span><span>group</span><span> </span><span>=</span><span> </span><span>type</span><span>,</span><span> </span><span>colour</span><span> </span><span>=</span><span> </span><span>type</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>0.8</span><span>)</span><span> </span><span>+</span><span>
</span><span>theme_bw</span><span>()</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;Time&quot;</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Load (kW)&quot;</span><span>,</span><span>
</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;Fit from GAM n.3&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-18" src="https://i1.wp.com/petolau.github.io/images/post_3/unnamed-chunk-18-1.png?w=456" data-recalc-dims="1"></p>
<p>Just a little differences in comparison to the <code>gam_2</code> model, looks like with <code>te</code> fit is more smoothed. Are we missing something? Of course! Function <code>te</code> has a default for a number of knots \( 5^d \), where d is a number of dimensions (variables), which is in our case too small. Now, I will set a number of knots to the maximal possible value <code>k = c(period, 7)</code>, what means that upper boundary for EDF (Estimated Degrees of Freedom) will be 48*7 – 1 = 335.</p>
<figure>
<pre><code data-lang="r"><span>gam_4</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>Load</span><span> </span><span>~</span><span> </span><span>te</span><span>(</span><span>Daily</span><span>,</span><span> </span><span>Weekly</span><span>,</span><span>
</span><span>k</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>period</span><span>,</span><span> </span><span>7</span><span>),</span><span>
</span><span>bs</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;cr&quot;</span><span>,</span><span> </span><span>&quot;ps&quot;</span><span>)),</span><span>
</span><span>data</span><span> </span><span>=</span><span> </span><span>matrix_gam</span><span>,</span><span>
</span><span>family</span><span> </span><span>=</span><span> </span><span>gaussian</span><span>)</span><span>
</span><span>summary</span><span>(</span><span>gam_4</span><span>)</span><span>$</span><span>r.sq</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] 0.9727604</code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>gam_4</span><span>)</span><span>$</span><span>sp.criterion</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##   GCV.Cp 
## 34839.46</code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>gam_4</span><span>)</span><span>$</span><span>s.table</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##                       edf   Ref.df        F p-value
## te(Daily,Weekly) 119.4117 149.6528 160.2065       0</code></pre>
</figure>
<p>We can see here that R-squared jumped little bit up and edf value has increased five times (!) in comparison to the model <code>gam_3</code>.<br>
Let’s plot fitted values:</p>
<figure>
<pre><code data-lang="r"><span>datas</span><span> </span><span>&lt;-</span><span> </span><span>rbindlist</span><span>(</span><span>list</span><span>(</span><span>data_r</span><span>[,</span><span> </span><span>.</span><span>(</span><span>value</span><span>,</span><span> </span><span>date_time</span><span>)],</span><span>
</span><span>data.table</span><span>(</span><span>value</span><span> </span><span>=</span><span> </span><span>gam_4</span><span>$</span><span>fitted.values</span><span>,</span><span>
</span><span>data_time</span><span> </span><span>=</span><span> </span><span>data_r</span><span>[,</span><span> </span><span>date_time</span><span>])))</span><span>
</span><span>datas</span><span>[,</span><span> </span><span>type</span><span> </span><span>:=</span><span> </span><span>c</span><span>(</span><span>rep</span><span>(</span><span>&quot;Real&quot;</span><span>,</span><span> </span><span>nrow</span><span>(</span><span>data_r</span><span>)),</span><span> </span><span>rep</span><span>(</span><span>&quot;Fitted&quot;</span><span>,</span><span> </span><span>nrow</span><span>(</span><span>data_r</span><span>)))]</span><span>
</span><span>ggplot</span><span>(</span><span>data</span><span> </span><span>=</span><span> </span><span>datas</span><span>,</span><span> </span><span>aes</span><span>(</span><span>date_time</span><span>,</span><span> </span><span>value</span><span>,</span><span> </span><span>group</span><span> </span><span>=</span><span> </span><span>type</span><span>,</span><span> </span><span>colour</span><span> </span><span>=</span><span> </span><span>type</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>0.8</span><span>)</span><span> </span><span>+</span><span>
</span><span>theme_bw</span><span>()</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;Time&quot;</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Load (kW)&quot;</span><span>,</span><span>
</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;Fit from GAM n.4&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-20" src="https://i2.wp.com/petolau.github.io/images/post_3/unnamed-chunk-20-1.png?w=456" data-recalc-dims="1"></p>
<p>This seems much better than it was with model <code>gam_3</code>.</p>
<p>Now I will prove my statement about the upper boundary for EDF. There is a possibility to fix number of smooth basis (i.e. EDF) in smooth terms. Just use an argument <code>fx = TRUE</code>.</p>
<figure>
<pre><code data-lang="r"><span>gam_4_fx</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>Load</span><span> </span><span>~</span><span> </span><span>te</span><span>(</span><span>Daily</span><span>,</span><span> </span><span>Weekly</span><span>,</span><span>
</span><span>k</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>period</span><span>,</span><span> </span><span>7</span><span>),</span><span>
</span><span>bs</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;cr&quot;</span><span>,</span><span> </span><span>&quot;ps&quot;</span><span>),</span><span>
</span><span>fx</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>),</span><span>
</span><span>data</span><span> </span><span>=</span><span> </span><span>matrix_gam</span><span>,</span><span>
</span><span>family</span><span> </span><span>=</span><span> </span><span>gaussian</span><span>)</span><span>
</span><span>summary</span><span>(</span><span>gam_4_fx</span><span>)</span><span>$</span><span>r.sq</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] 0.965618</code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>gam_4_fx</span><span>)</span><span>$</span><span>s.table</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##                  edf Ref.df        F       p-value
## te(Daily,Weekly) 335    335 57.25389 5.289648e-199</code></pre>
</figure>
<p>EDF = 335, here we are. We can see that R-squared is lower than with the model <code>gam_4</code>, it is due to that 335 is too high and we <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Overfitting">overfitted</a> the model. It is prove to that GCV procedure (estimation of lambda and EDF) is working properly. You can read more about the optimal setting (choosing) of <code>k</code> argument at <a target="_blank" rel="nofollow" href="https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/choose.k.html"><code>?choose.k</code></a> and of course in <a target="_blank" rel="nofollow" href="https://www.crcpress.com/Generalized-Additive-Models-An-Introduction-with-R/Wood/p/book/9781584884743">the book from Simon Wood</a>.</p>
<p>With the <code>mgcv</code> package we have two more opportunities, methods, how to include tensor product interactions term – with functions <code>ti</code> and <code>t2</code>. <code>ti</code> produces a tensor product interaction, appropriate when the main effects (and any lower interactions) are also present, while <code>te</code> produces a full tensor product smooth. <code>t2</code> is an alternative function to <code>te</code> and uses different penalization method. You can read more about it in a documentation of the package <code>mgcv</code>: use <code>?ti</code> and <code>?t2</code> to see more.</p>
<p>So, let’s try both methods in the our case (model). First, let’s use <code>ti</code>:</p>
<figure>
<pre><code data-lang="r"><span>gam_5</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>Load</span><span> </span><span>~</span><span> </span><span>s</span><span>(</span><span>Daily</span><span>,</span><span> </span><span>bs</span><span> </span><span>=</span><span> </span><span>&quot;cr&quot;</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>period</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>Weekly</span><span>,</span><span> </span><span>bs</span><span> </span><span>=</span><span> </span><span>&quot;ps&quot;</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>7</span><span>)</span><span> </span><span>+</span><span>
</span><span>ti</span><span>(</span><span>Daily</span><span>,</span><span> </span><span>Weekly</span><span>,</span><span>
</span><span>k</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>period</span><span>,</span><span> </span><span>7</span><span>),</span><span>
</span><span>bs</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;cr&quot;</span><span>,</span><span> </span><span>&quot;ps&quot;</span><span>)),</span><span>
</span><span>data</span><span> </span><span>=</span><span> </span><span>matrix_gam</span><span>,</span><span>
</span><span>family</span><span> </span><span>=</span><span> </span><span>gaussian</span><span>)</span><span>
</span><span>summary</span><span>(</span><span>gam_5</span><span>)</span><span>$</span><span>r.sq</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] 0.9717469</code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>gam_5</span><span>)</span><span>$</span><span>sp.criterion</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##   GCV.Cp 
## 35772.35</code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>gam_5</span><span>)</span><span>$</span><span>s.table</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##                        edf     Ref.df          F p-value
## s(Daily)         22.583649  27.964970  444.19962       0
## s(Weekly)         5.914531   5.995934 1014.72482       0
## ti(Daily,Weekly) 85.310314 110.828814   41.22288       0</code></pre>
</figure>
<p>Then let’s use <code>t2</code>. I set argument <code>full = TRUE</code> because it gives strict invariance to penalties.</p>
<figure>
<pre><code data-lang="r"><span>gam_6</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>Load</span><span> </span><span>~</span><span> </span><span>t</span><span>2</span><span>(</span><span>Daily</span><span>,</span><span> </span><span>Weekly</span><span>,</span><span>
</span><span>k</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>period</span><span>,</span><span> </span><span>7</span><span>),</span><span>
</span><span>bs</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;cr&quot;</span><span>,</span><span> </span><span>&quot;ps&quot;</span><span>),</span><span>
</span><span>full</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>),</span><span>
</span><span>data</span><span> </span><span>=</span><span> </span><span>matrix_gam</span><span>,</span><span>
</span><span>family</span><span> </span><span>=</span><span> </span><span>gaussian</span><span>)</span><span>
</span><span>summary</span><span>(</span><span>gam_6</span><span>)</span><span>$</span><span>r.sq</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] 0.9738273</code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>gam_6</span><span>)</span><span>$</span><span>sp.criterion</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##   GCV.Cp 
## 32230.68</code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>gam_6</span><span>)</span><span>$</span><span>s.table</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##                       edf   Ref.df        F p-value
## t2(Daily,Weekly) 98.12005 120.2345 86.70754       0</code></pre>
</figure>
<p>I have also printed the GCV score value for the last three models, which is also a good criterion to choose optimal model among a set of fitted models. We can see that with <code>t2</code> term and corresponding model <code>gam_6</code> the GCV value is the lowest. It may be indicating that <code>gam_6</code> is our best model so far.</p>
<p>Other model selection criterion, which is widely used in statistics, is AIC (<a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Akaike_information_criterion">Akaike Information Criterion</a>). Criterion has equation \( AIC = 2k – 2\ln (\hat{L}) \), where k is the number of parameters to be estimated and \(\hat{L}\) is the maximized value of the likelihood function of the model. So lower values are better for AIC. Let’s look at them for our three models:</p>
<figure>
<pre><code data-lang="r"><span>AIC</span><span>(</span><span>gam_4</span><span>,</span><span> </span><span>gam_5</span><span>,</span><span> </span><span>gam_6</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##             df      AIC
## gam_4 121.4117 8912.611
## gam_5 115.8085 8932.746
## gam_6 100.1200 8868.628</code></pre>
</figure>
<p>The lowest value is in the <code>gam_6</code> model, so again it’s the winning model. Let’s again look at fitted values to be sure that everything is OK.</p>
<figure>
<pre><code data-lang="r"><span>datas</span><span> </span><span>&lt;-</span><span> </span><span>rbindlist</span><span>(</span><span>list</span><span>(</span><span>data_r</span><span>[,</span><span> </span><span>.</span><span>(</span><span>value</span><span>,</span><span> </span><span>date_time</span><span>)],</span><span>
</span><span>data.table</span><span>(</span><span>value</span><span> </span><span>=</span><span> </span><span>gam_6</span><span>$</span><span>fitted.values</span><span>,</span><span>
</span><span>data_time</span><span> </span><span>=</span><span> </span><span>data_r</span><span>[,</span><span> </span><span>date_time</span><span>])))</span><span>
</span><span>datas</span><span>[,</span><span> </span><span>type</span><span> </span><span>:=</span><span> </span><span>c</span><span>(</span><span>rep</span><span>(</span><span>&quot;Real&quot;</span><span>,</span><span> </span><span>nrow</span><span>(</span><span>data_r</span><span>)),</span><span> </span><span>rep</span><span>(</span><span>&quot;Fitted&quot;</span><span>,</span><span> </span><span>nrow</span><span>(</span><span>data_r</span><span>)))]</span><span>
</span><span>ggplot</span><span>(</span><span>data</span><span> </span><span>=</span><span> </span><span>datas</span><span>,</span><span> </span><span>aes</span><span>(</span><span>date_time</span><span>,</span><span> </span><span>value</span><span>,</span><span> </span><span>group</span><span> </span><span>=</span><span> </span><span>type</span><span>,</span><span> </span><span>colour</span><span> </span><span>=</span><span> </span><span>type</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>0.8</span><span>)</span><span> </span><span>+</span><span>
</span><span>theme_bw</span><span>()</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;Time&quot;</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Load (kW)&quot;</span><span>,</span><span>
</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;Fit from GAM n.6&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-25" src="https://i1.wp.com/petolau.github.io/images/post_3/unnamed-chunk-25-1.png?w=456" data-recalc-dims="1"></p>
<p>Seems OK. We can see that fitted values for models <code>gam_4</code> and <code>gam_6</code> are very similar. It’s very difficult to choose which model is better. It’s possible to use more visualization and model diagnostic capabilities of package <code>mgcv</code> to compare these two models.</p>
<p>The first one is function <code>gam.check</code>, which makes four plots: QQ-plot of residuals, linear predictor vs. residuals, the histogram of residuals and the plot of fitted values vs. response. Let’s make them for models <code>gam_4</code> and <code>gam_6</code>.</p>
<figure>
<pre><code data-lang="r"><span>gam.check</span><span>(</span><span>gam_4</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-26" src="https://i0.wp.com/petolau.github.io/images/post_3/unnamed-chunk-26-1.png?w=456" data-recalc-dims="1"></p>
<figure>
<pre><code data-lang="text">## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 7 iterations.
## The RMS GCV score gradiant at convergence was 0.2833304 .
## The Hessian was positive definite.
## The estimated model rank was 336 (maximum possible: 336)
## Model rank =  336 / 336 
## 
## Basis dimension (k) checking results. Low p-value (k-index&lt;1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                      k'    edf k-index p-value
## te(Daily,Weekly) 335.00 119.41    1.22       1</code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>gam.check</span><span>(</span><span>gam_6</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-27" src="https://i1.wp.com/petolau.github.io/images/post_3/unnamed-chunk-27-1.png?w=456" data-recalc-dims="1"></p>
<figure>
<pre><code data-lang="text">## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 9 iterations.
## The RMS GCV score gradiant at convergence was 0.05208856 .
## The Hessian was positive definite.
## The estimated model rank was 336 (maximum possible: 336)
## Model rank =  336 / 336 
## 
## Basis dimension (k) checking results. Low p-value (k-index&lt;1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                      k'    edf k-index p-value
## t2(Daily,Weekly) 335.00  98.12    1.18       1</code></pre>
</figure>
<p>The function <code>gam.check</code> makes also output to the console of more useful information. We can see again that models are very similar, just in histograms some differences can be seen, but it’s also insignificant.</p>
<p>We didn’t use so far the default <code>plot</code> function to <code>gam</code> object for models with tensor product interactions. Let’s use it and explore what it brings us.</p>
<figure>
<pre><code data-lang="r"><span>layout</span><span>(</span><span>matrix</span><span>(</span><span>1</span><span>:</span><span>2</span><span>,</span><span> </span><span>nrow</span><span> </span><span>=</span><span> </span><span>1</span><span>))</span><span>
</span><span>plot</span><span>(</span><span>gam_4</span><span>,</span><span> </span><span>rug</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>,</span><span> </span><span>se</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>,</span><span> </span><span>n</span><span>2</span><span> </span><span>=</span><span> </span><span>80</span><span>,</span><span> </span><span>main</span><span> </span><span>=</span><span> </span><span>&quot;gam n.4 with te()&quot;</span><span>)</span><span>
</span><span>plot</span><span>(</span><span>gam_6</span><span>,</span><span> </span><span>rug</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>,</span><span> </span><span>se</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>,</span><span> </span><span>n</span><span>2</span><span> </span><span>=</span><span> </span><span>80</span><span>,</span><span> </span><span>main</span><span> </span><span>=</span><span> </span><span>&quot;gam n.6 with t2()&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-28" src="https://i2.wp.com/petolau.github.io/images/post_3/unnamed-chunk-28-1.png?w=456" data-recalc-dims="1"></p>
<p>It’s a nice contour line plot. On axes are our independent variables, contour lines and corresponding numbers on them represents effects of ind. variables to the response variable. Now it is possible to see some little differences. The model <code>gam_6</code> with <code>t2</code> have more “wavy” contours. So it implies that it adapts more to response variable and smoothing factor is lower. In the next parts of the post, the model <code>gam_6</code> will be used for the analysis.</p>
<p>Another great feature of the package <code>mgcv</code> is the plotting function <code>vis.gam</code>. It makes 3D view (or 2D) of surface of fitted values according to independent variables. Let’s look at it.</p>
<figure>
<pre><code data-lang="r"><span>vis.gam</span><span>(</span><span>gam_6</span><span>,</span><span> </span><span>n.grid</span><span> </span><span>=</span><span> </span><span>50</span><span>,</span><span> </span><span>theta</span><span> </span><span>=</span><span> </span><span>35</span><span>,</span><span> </span><span>phi</span><span> </span><span>=</span><span> </span><span>32</span><span>,</span><span> </span><span>zlab</span><span> </span><span>=</span><span> </span><span>&quot;&quot;</span><span>,</span><span>
</span><span>ticktype</span><span> </span><span>=</span><span> </span><span>&quot;detailed&quot;</span><span>,</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;topo&quot;</span><span>,</span><span> </span><span>main</span><span> </span><span>=</span><span> </span><span>&quot;t2(D, W)&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-29" src="https://i0.wp.com/petolau.github.io/images/post_3/unnamed-chunk-29-1.png?w=456" data-recalc-dims="1"></p>
<p>We can see that the highest peak is when the Daily variable has values near 30 (3 p.m.) and the Weekly variable has value 1 (it is a Monday).</p>
<p>2D view, contour lines type of plot, can be visualized too. Just set argument <code>plot.type = &quot;contour&quot;</code>.</p>
<figure>
<pre><code data-lang="r"><span>vis.gam</span><span>(</span><span>gam_6</span><span>,</span><span> </span><span>main</span><span> </span><span>=</span><span> </span><span>&quot;t2(D, W)&quot;</span><span>,</span><span> </span><span>plot.type</span><span> </span><span>=</span><span> </span><span>&quot;contour&quot;</span><span>,</span><span>
</span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;terrain&quot;</span><span>,</span><span> </span><span>contour.col</span><span> </span><span>=</span><span> </span><span>&quot;black&quot;</span><span>,</span><span> </span><span>lwd</span><span> </span><span>=</span><span> </span><span>2</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-30" src="https://i2.wp.com/petolau.github.io/images/post_3/unnamed-chunk-30-1.png?w=456" data-recalc-dims="1"></p>
<p>Again we can see that the highest value of electricity load is on Monday at 3:00 pm, it is very similar till Thursday, then the load is decreasing (weekends).</p>
<p>Another interesting contribution to the <strong>time series regression model</strong> can be an inclusion of <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Autoregressive_model"><strong>autoregressive model</strong></a> for serially correlated errors. Let’s look at it.</p>
<h3>Autoregressive models and GAMS</h3>
<p><strong>GAM</strong> or <strong>MLR</strong> have assumptions in a model that errors (residuals) are identically and independently distributed (i.i.d.). In the case of the time series regression, it is very strong assumption, which is here, logically, not fulfilled. Present time series values are highly correlated with past values, so errors of the model will be correlated too. This phenomenon is called an <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Autocorrelation">autocorrelation</a>. This implies that estimated regression coefficients and residuals of a model might be negatively biased, which also implies that previously computed p-values of statistical tests or confidence intervals are wrong.</p>
<p>How can we handle this situation? By inclusion of autoregressive model (AR) for errors in our model. So we have the model with the term for errors like this:</p>
<p></p>
<p>where the second equation is a classical AR(1) process and \( \phi \) is an unknown autoregressive coefficient to be estimated. Errors can be also nested within a week, which is in our case more appropriate, because of the double seasonal character of our time series. You can add also higher orders of AR process and also MA (<a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Moving-average_model">moving average</a>) model. You can read more about estimation and modeling of this kind of models in an excellent book by <a target="_blank" rel="nofollow" href="http://onlinelibrary.wiley.com/book/10.1002/9781118619193">Box, Jenkins, and Reinsel: Time Series Analysis</a>.</p>
<p>It’s possible to add correlation term for errors with function <code>gamm</code>, which stands for <strong>GAM</strong> mixture models. It calls the <code>lme</code> function from package <code>nlme</code>. Now, train basic model with function <code>gamm</code> and an another model with AR(1) process nested within a week – just add this argument to <code>gamm</code> function: <code>correlation = corARMA(form = ~ 1|Weekly, p = 1)</code>.</p>
<figure>
<pre><code data-lang="r"><span>gam_6_ar0</span><span> </span><span>&lt;-</span><span> </span><span>gamm</span><span>(</span><span>Load</span><span> </span><span>~</span><span> </span><span>t</span><span>2</span><span>(</span><span>Daily</span><span>,</span><span> </span><span>Weekly</span><span>,</span><span>
</span><span>k</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>period</span><span>,</span><span> </span><span>7</span><span>),</span><span>
</span><span>bs</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;cr&quot;</span><span>,</span><span> </span><span>&quot;ps&quot;</span><span>),</span><span>
</span><span>full</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>),</span><span>
</span><span>data</span><span> </span><span>=</span><span> </span><span>matrix_gam</span><span>,</span><span>
</span><span>family</span><span> </span><span>=</span><span> </span><span>gaussian</span><span>,</span><span>
</span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;REML&quot;</span><span>)</span><span>
</span><span>gam_6_ar1</span><span> </span><span>&lt;-</span><span> </span><span>gamm</span><span>(</span><span>Load</span><span> </span><span>~</span><span> </span><span>t</span><span>2</span><span>(</span><span>Daily</span><span>,</span><span> </span><span>Weekly</span><span>,</span><span>
</span><span>k</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>period</span><span>,</span><span> </span><span>7</span><span>),</span><span>
</span><span>bs</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;cr&quot;</span><span>,</span><span> </span><span>&quot;ps&quot;</span><span>),</span><span>
</span><span>full</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>),</span><span>
</span><span>data</span><span> </span><span>=</span><span> </span><span>matrix_gam</span><span>,</span><span>
</span><span>family</span><span> </span><span>=</span><span> </span><span>gaussian</span><span>,</span><span>
</span><span>correlation</span><span> </span><span>=</span><span> </span><span>corARMA</span><span>(</span><span>form</span><span> </span><span>=</span><span> </span><span>~</span><span> </span><span>1</span><span>|</span><span>Weekly</span><span>,</span><span> </span><span>p</span><span> </span><span>=</span><span> </span><span>1</span><span>),</span><span>
</span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;REML&quot;</span><span>)</span></code></pre>
</figure>
<p>We can use an <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Analysis_of_variance">ANOVA</a> to compare these two models and pick the right one. It does generalized likelihood ratio test to our nested models.</p>
<figure>
<pre><code data-lang="r"><span>anova</span><span>(</span><span>gam_6_ar0</span><span>$</span><span>lme</span><span>,</span><span> </span><span>gam_6_ar1</span><span>$</span><span>lme</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##               Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## gam_6_ar0$lme     1 10 9070.525 9115.568 -4525.263                        
## gam_6_ar1$lme     2 11 8751.361 8800.908 -4364.680 1 vs 2 321.1644  &lt;.0001</code></pre>
</figure>
<p>AIC value of the model with AR(1) term is lower, so it seems better, also p-value is very low which indicates that second model is better to choose.</p>
<p>Estimated \( \phi \) coefficient of AR(1) process can be seen here:</p>
<figure>
<pre><code data-lang="r"><span>intervals</span><span>(</span><span>gam_6_ar1</span><span>$</span><span>lme</span><span>,</span><span> </span><span>which</span><span> </span><span>=</span><span> </span><span>&quot;var-cov&quot;</span><span>)</span><span>$</span><span>corStruct</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##         lower      est.     upper
## Phi 0.6363101 0.7107914 0.7721434
## attr(,&quot;label&quot;)
## [1] &quot;Correlation structure:&quot;</code></pre>
</figure>
<p>Value 0.71 is pretty high, which indicates a strong dependency on previous values of errors (lags).</p>
<p>Now let’s look another diagnostic for these two models. Plot of values of a <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Partial_autocorrelation_function">partial autocorrelation function</a> applied to normalized residuals. It gives the partial correlation of residuals with its own lagged values, controlling for the values of the residuals at all shorter lags.</p>
<figure>
<pre><code data-lang="r"><span>layout</span><span>(</span><span>matrix</span><span>(</span><span>1</span><span>:</span><span>2</span><span>,</span><span> </span><span>ncol</span><span> </span><span>=</span><span> </span><span>2</span><span>))</span><span>
</span><span>pacf</span><span>(</span><span>resid</span><span>(</span><span>gam_6_ar0</span><span>$</span><span>lme</span><span>,</span><span> </span><span>type</span><span> </span><span>=</span><span> </span><span>&quot;normalized&quot;</span><span>),</span><span> </span><span>lag.max</span><span> </span><span>=</span><span> </span><span>48</span><span>,</span><span> </span><span>main</span><span> </span><span>=</span><span> </span><span>&quot;pACF of gam n.6&quot;</span><span>)</span><span>
</span><span>pacf</span><span>(</span><span>resid</span><span>(</span><span>gam_6_ar1</span><span>$</span><span>lme</span><span>,</span><span> </span><span>type</span><span> </span><span>=</span><span> </span><span>&quot;normalized&quot;</span><span>),</span><span> </span><span>lag.max</span><span> </span><span>=</span><span> </span><span>48</span><span>,</span><span> </span><span>main</span><span> </span><span>=</span><span> </span><span>&quot;pACF of gam n.6 with AR(1)&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-34" src="https://i2.wp.com/petolau.github.io/images/post_3/unnamed-chunk-34-1.png?w=456" data-recalc-dims="1"></p>
<p>We can see a significant difference between these two plots and corresponding values of pACF. Optimal values of pACF should be under dashed blue lines, which isn’t completely this scenario on the right plot (model with the AR(1)). It can be better…</p>
<p>I will do little hack now. For optimal choose of AR(p) and MA(q) orders <code>auto.arima</code> function, from the library <code>forecast</code>, will be used on residuals from the model <code>gam_6_ar0</code>. It automatically chooses optimal orders of ARMA (in our case) based on AIC criterion. As we can use just ARMA models in <code>gamm</code>, so nonstationarity isn’t allowed, set an argument <code>stationary = TRUE</code>.</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>forecast</span><span>)</span><span>
</span><span>arma_res</span><span> </span><span>&lt;-</span><span> </span><span>auto.arima</span><span>(</span><span>resid</span><span>(</span><span>gam_6_ar0</span><span>$</span><span>lme</span><span>,</span><span> </span><span>type</span><span> </span><span>=</span><span> </span><span>&quot;normalized&quot;</span><span>),</span><span>
</span><span>stationary</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>,</span><span> </span><span>seasonal</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>)</span><span>
</span><span>arma_res</span><span>$</span><span>coef</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##        ar1        ma1        ma2 
##  0.9102517 -0.3763428 -0.2862549</code></pre>
</figure>
<p>It’s interesting. <code>auto.arima</code> chooses also two orders of MA model. Let’s fit a new model and again do ANOVA.</p>
<figure>
<pre><code data-lang="r"><span>gam_6_arma</span><span>&lt;-</span><span> </span><span>gamm</span><span>(</span><span>Load</span><span> </span><span>~</span><span> </span><span>t</span><span>2</span><span>(</span><span>Daily</span><span>,</span><span> </span><span>Weekly</span><span>,</span><span>
</span><span>k</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>period</span><span>,</span><span> </span><span>7</span><span>),</span><span>
</span><span>bs</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;cr&quot;</span><span>,</span><span> </span><span>&quot;ps&quot;</span><span>),</span><span>
</span><span>full</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>),</span><span>
</span><span>data</span><span> </span><span>=</span><span> </span><span>matrix_gam</span><span>,</span><span>
</span><span>family</span><span> </span><span>=</span><span> </span><span>gaussian</span><span>,</span><span>
</span><span>correlation</span><span> </span><span>=</span><span> </span><span>corARMA</span><span>(</span><span>form</span><span> </span><span>=</span><span> </span><span>~</span><span> </span><span>1</span><span>|</span><span>Weekly</span><span>,</span><span> </span><span>p</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>q</span><span> </span><span>=</span><span> </span><span>2</span><span>),</span><span>
</span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;REML&quot;</span><span>)</span><span>
</span><span>anova</span><span>(</span><span>gam_6_ar0</span><span>$</span><span>lme</span><span>,</span><span> </span><span>gam_6_ar1</span><span>$</span><span>lme</span><span>,</span><span> </span><span>gam_6_arma</span><span>$</span><span>lme</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##                Model df      AIC      BIC    logLik   Test  L.Ratio
## gam_6_ar0$lme      1 10 9070.525 9115.568 -4525.263                
## gam_6_ar1$lme      2 11 8751.361 8800.908 -4364.680 1 vs 2 321.1644
## gam_6_arma$lme     3 13 8731.750 8790.306 -4352.875 2 vs 3  23.6107
##                p-value
## gam_6_ar0$lme         
## gam_6_ar1$lme   &lt;.0001
## gam_6_arma$lme  &lt;.0001</code></pre>
</figure>
<p>The new model <code>gam_6_arma</code> is slightly better in the meaning of AIC. P-value is again very small, which means that <code>gam_6_arma</code> is significantly better than <code>gam_6_ar1</code>.</p>
<p>Now do a comparison of these two models on theirs pACF values.</p>
<figure>
<pre><code data-lang="r"><span>layout</span><span>(</span><span>matrix</span><span>(</span><span>1</span><span>:</span><span>2</span><span>,</span><span> </span><span>ncol</span><span> </span><span>=</span><span> </span><span>2</span><span>))</span><span>
</span><span>pacf</span><span>(</span><span>resid</span><span>(</span><span>gam_6_ar1</span><span>$</span><span>lme</span><span>,</span><span> </span><span>type</span><span> </span><span>=</span><span> </span><span>&quot;normalized&quot;</span><span>),</span><span> </span><span>lag.max</span><span> </span><span>=</span><span> </span><span>48</span><span>,</span><span> </span><span>main</span><span> </span><span>=</span><span> </span><span>&quot;pACF of gam n.6 with AR(1)&quot;</span><span>)</span><span>
</span><span>pacf</span><span>(</span><span>resid</span><span>(</span><span>gam_6_arma</span><span>$</span><span>lme</span><span>,</span><span> </span><span>type</span><span> </span><span>=</span><span> </span><span>&quot;normalized&quot;</span><span>),</span><span> </span><span>lag.max</span><span> </span><span>=</span><span> </span><span>48</span><span>,</span><span> </span><span>main</span><span> </span><span>=</span><span> </span><span>&quot;pACF of gam n.6 with ARMA(1,2)&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-37" src="https://i0.wp.com/petolau.github.io/images/post_3/unnamed-chunk-37-1.png?w=456" data-recalc-dims="1"></p>
<p>It’s better now. Can we be happy? Up to now, everything seems fine.</p>
<p>Let’s make one more visualization of residuals of two models to confirm our little uncertainty about the correctness of the model with ARMA(1,2). Compare it with baseline model <code>gam_6_ar0</code>.</p>
<figure>
<pre><code data-lang="r"><span>datas</span><span> </span><span>&lt;-</span><span> </span><span>data.table</span><span>(</span><span>Fitted_values</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>gam_6_ar0</span><span>$</span><span>gam</span><span>$</span><span>fitted.values</span><span>,</span><span>
</span><span>gam_6_arma</span><span>$</span><span>gam</span><span>$</span><span>fitted.values</span><span>),</span><span>
</span><span>Residuals</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>gam_6_ar0</span><span>$</span><span>gam</span><span>$</span><span>residuals</span><span>,</span><span>
</span><span>gam_6_arma</span><span>$</span><span>gam</span><span>$</span><span>residuals</span><span>),</span><span>
</span><span>Model</span><span> </span><span>=</span><span> </span><span>rep</span><span>(</span><span>c</span><span>(</span><span>&quot;Gam n.6&quot;</span><span>,</span><span> </span><span>&quot;Gam n.6 with ARMA(1,2)&quot;</span><span>),</span><span> </span><span>each</span><span> </span><span>=</span><span> </span><span>nrow</span><span>(</span><span>data_r</span><span>)))</span><span>
</span><span>ggplot</span><span>(</span><span>data</span><span> </span><span>=</span><span> </span><span>datas</span><span>,</span><span>
</span><span>aes</span><span>(</span><span>Fitted_values</span><span>,</span><span> </span><span>Residuals</span><span>))</span><span> </span><span>+</span><span>
</span><span>facet_grid</span><span>(</span><span>Model</span><span>~</span><span>.</span><span>,</span><span> </span><span>switch</span><span> </span><span>=</span><span> </span><span>&quot;y&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>geom_point</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>1.7</span><span>)</span><span> </span><span>+</span><span>
</span><span>geom_smooth</span><span>(</span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;loess&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>geom_hline</span><span>(</span><span>yintercept</span><span> </span><span>=</span><span> </span><span>0</span><span>,</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;red&quot;</span><span>,</span><span> </span><span>size</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span> </span><span>+</span><span>
</span><span>theme</span><span>(</span><span>panel.border</span><span> </span><span>=</span><span> </span><span>element_rect</span><span>(</span><span>fill</span><span> </span><span>=</span><span> </span><span>NA</span><span>,</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;black&quot;</span><span>),</span><span>
</span><span>plot.title</span><span> </span><span>=</span><span> </span><span>element_text</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>14</span><span>,</span><span> </span><span>hjust</span><span> </span><span>=</span><span> </span><span>0.5</span><span>,</span><span> </span><span>face</span><span> </span><span>=</span><span> </span><span>&quot;bold&quot;</span><span>),</span><span>
</span><span>axis.text</span><span> </span><span>=</span><span> </span><span>element_text</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>11</span><span>),</span><span>
</span><span>axis.title</span><span> </span><span>=</span><span> </span><span>element_text</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>12</span><span>,</span><span> </span><span>face</span><span> </span><span>=</span><span> </span><span>&quot;bold&quot;</span><span>),</span><span>
</span><span>strip.text</span><span> </span><span>=</span><span> </span><span>element_text</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>13</span><span>,</span><span> </span><span>face</span><span> </span><span>=</span><span> </span><span>&quot;bold&quot;</span><span>),</span><span>
</span><span>strip.background</span><span> </span><span>=</span><span> </span><span>element_rect</span><span>(</span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;black&quot;</span><span>))</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;Fitted values vs Residuals of two models&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="plot of chunk unnamed-chunk-38" src="https://i1.wp.com/petolau.github.io/images/post_3/unnamed-chunk-38-1.png?w=456" data-recalc-dims="1"></p>
<p>Here we are. The model with an ARMA(1,2) process has residuals at low fitted values somehow correlated, which is not good (heteroscedasticity). What now? I have to be honest, I can’t explain why this interesting thing happened. So, I would be very happy when someone guides me to the right direction and writes some notes to the discussion, it would really help. One more note about the inclusion of ARMA models to <strong>GAM</strong> or <strong>MLR</strong>. It didn’t help in the meaning of forecast accuracy, I made a lot of experiments, but MAPE was higher than with models without ARMAs. The big disadvantage is also higher time complexity than with the simple model without ARMA.</p>
<h3>Animations and conclusions</h3>
<p>I will end this post with a visualization analysis of the model <code>gam_6</code>. I will try to do some animated dashboards of main characteristics of the model thru time. It should be very interesting how electricity load changes and how the model can adapt to these changes. In the first two animations you can see three plots on the one figure:</p>
<ol>
<li>Fitted values by GAM and real response values of electricity consumption.</li>
<li>Fitted values vs. residuals of GAM.</li>
<li>Forecast for one week ahead vs. real consumption plus forecast accuracy in MAPE.</li>
</ol>
<p>Animations were created by functions from package <code>grid</code> to layout ggplots to one figure, and function <code>saveGIF</code> from package <code>animation</code> to create the final GIF. Here there are:</p>
<p><img alt="" src="https://i2.wp.com/petolau.github.io/images/post_3/industry_1_dashboard_v2.gif?w=456" data-recalc-dims="1"></p>
<p><img alt="" src="https://i0.wp.com/petolau.github.io/images/post_3/industry_4_dashboard_v2.gif?w=456" data-recalc-dims="1"></p>
<p>We can see the behavior of <strong>GAMs</strong> on two times series from two different types of industries: commercial property and light industrial. Notice (look), that fitted values and forecasts are smoother than where were with <a target="_blank" rel="nofollow" href="https://petolau.github.io/Forecast-double-seasonal-time-series-with-multiple-linear-regression-in-R/">Multiple linear regression from the previous post</a>. It’s good information to know, it implies that our model isn’t overfitted. On the other hand, I also compared forecast performance of <strong>GAM</strong> and <strong>MLR</strong>, but GAM was slightly worse than MLR, which can be a little bit surprising for somebody. Seems that usage of penalized least squares (GAM) for forecasting time series doesn’t bring improvement against classical OLS (MLR). But there are other advantages, which I will point out in the end.</p>
<p>Another animation, which I created, is a 3D visualization of fitted values like a surface. As I showed you before, it can be done very easily by function <code>vis.gam</code>.</p>
<p><img alt="" src="https://i0.wp.com/petolau.github.io/images/post_3/industry_1_v2.gif?w=456" data-recalc-dims="1"></p>
<p>We can compare it with the first dashboard because again it’s time series from commercial properties. EDF is also printed to the title of GIF for reasons to see how it changes depending on the behavior of electricity consumption.</p>
<p>With this, I would like to end the main part of this tutorial and conclude it with some remarks.</p>
<p>In the beginning of the post, motivation and the theory behind GAM method was introduced. Next, the analytical (“magical”) part continued with explaining various features of the package <code>mgcv</code>. Different types of <strong>interactions</strong> were showed and analyzed – best chose was tensor product interactions (<code>t2</code>). Next, consideration about the correctness of an autoregressive model (or ARMA) for errors was discussed – ended with an open question. In the end, I tried to make a big view of a behavior of <strong>GAMs</strong> by animations of its performance on <strong>double seasonal time series</strong>.</p>
<p>Everything that was said implies these advantages and disadvantages of using <strong>GAM</strong> for your problem:</p>
<p>Advantages:</p>
<ul>
<li>Many possibilities how to model your independent variables (many types of smooth functions).</li>
<li>Model response variable by all known distribution families (Normal, Gamma, Poisson etc.).</li>
<li>Define interactions by tensor product, which means usage of different basis functions for interacted variables.</li>
<li>Test statistical significance of a non-linear relationship of the independent variable to dependent.</li>
</ul>
<p>Disadvantages:</p>
<ul>
<li>Complex method, difficult to learn.</li>
<li>It’s sometimes difficult to interpret results.</li>
<li>Several parameters to tune, which have to be set carefully.</li>
</ul>
<p>In the future post, I want to focus on some machine learning method like Support Vector Regression or Regression Trees with alongside of grid search for tuning parameters.</p>
<p>The script (<em>3post-GAM.R</em>) for a creation of the whole tutorial is located on <a target="_blank" rel="nofollow" href="https://github.com/PetoLau/petolau.github.io/tree/master/_Rscripts">my GitHub repository</a>. I would like to encourage you to add any feedback to the discussion below. Thank you for reading this long post.</p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://petolau.github.io/Analyzing-double-seasonal-time-series-with-GAM-in-R/"> Peter Laurinec</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/vwexLvQXb0U"></a></div>