<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/greybox-0-3-0-whats-new/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://forecasting.svetunkov.ru/en/2018/08/07/greybox-0-3-0-whats-new/"> R – Modern Forecasting</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Three months have passed since the initial release of </p>
<pre>greybox</pre>
<p> on CRAN. I would not say that the package develops like crazy, but there have been some changes since May. Let’s have a look. We start by loading both </p>
<pre>greybox</pre>
<p> and </p>
<pre>smooth</pre>
<p>:</p>
<pre>library(greybox)
library(smooth)</pre>

<h3>Rolling Origin</h3>
<p>First of all, </p>
<pre>ro()</pre>
<p> function now has its own class and works with </p>
<pre>plot()</pre>
<p> function, so that you can have a visual representation of the results. Here’s an example:</p>
<pre>x &lt;- rnorm(100,100,10)
ourCall &lt;- &quot;es(data, h=h, intervals=TRUE)&quot;
ourValue &lt;- c(&quot;forecast&quot;, &quot;lower&quot;, &quot;upper&quot;)
ourRO &lt;- ro(x,h=20,origins=5,ourCall,ourValue,co=TRUE)
plot(ourRO)</pre>

<div class="wp-caption"><a rel="nofollow" href="https://i2.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/roPlotExample.png?ssl=1" target="_blank"><img sizes="(max-width: 300px) 100vw, 300px" src="https://i2.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/roPlotExample-300x175.png?resize=300%2C175&ssl=1" alt="" width="300" class="wp-image-1781" data-recalc-dims="1" height="175"></a>
<p class="wp-caption-text">Example of the plot of rolling origin function</p>
</div>
<p>Each point on the produced graph corresponds to an origin and straight lines correspond to the forecasts. Given that we asked for point forecasts and for lower and upper bounds of prediction interval, we have three respective lines. By plotting the results of rolling origin experiment, we can see if the model is stable or not. Just compare the previous graph with the one produced from the call to Holt’s model:</p>
<pre>ourCall &lt;- &quot;es(data, model='AAN', h=h, intervals=TRUE)&quot;
ourRO &lt;- ro(x,h=20,origins=5,ourCall,ourValue,co=TRUE)
plot(ourRO)</pre>
<div class="wp-caption"><a rel="nofollow" href="https://i2.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/roPlotExampleAAN.png?ssl=1" target="_blank"><img sizes="(max-width: 300px) 100vw, 300px" src="https://i1.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/roPlotExampleAAN-300x175.png?resize=300%2C175&ssl=1" alt="" width="300" class="wp-image-1782" data-recalc-dims="1" height="175"></a>
<p class="wp-caption-text">Example of the plot of rolling origin function with ETS(A,A,N)</p>
</div>
<p>Holt’s model is not suitable for this time series, so it’s forecasts are less stable than the forecasts of the automatically selected model in the previous case (which is ETS(A,N,N)).</p>
<p>Once again, there is a vignette with examples for the </p>
<pre>ro()</pre>
<p> function, <a rel="nofollow" href="https://cran.r-project.org/web/packages/greybox/vignettes/ro.html" target="_blank">have a look</a> if you want to know more.</p>
<h3>ALM – Advanced Linear Model</h3>
<p>Yes, there is “Generalised Linear Model” in R, which implements Poisson, Gamma, Binomial and other regressions. Yes, there are smaller packages, implementing models with more exotic distributions. But I needed several regression models with: Laplace distribution, Folded normal distribution, Chi-squared distribution and one new mysterious distribution, which is currently called “S distribution”. I needed them in one place and in one format: properly estimated using likelihoods, returning confidence intervals, information criteria and being able to produce forecasts. I also wanted them to work similar to </p>
<pre>lm()</pre>
<p>, so that the learning curve would not be too steep. So, here it is, the function </p>
<pre>alm()</pre>
<p>. It works quite similar to </p>
<pre>lm()</pre>
<p>:</p>
<pre>xreg &lt;- cbind(rfnorm(100,1,10),rnorm(100,50,5))
xreg &lt;- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rlaplace(100,0,3),xreg,rnorm(100,300,10))
colnames(xreg) &lt;- c(&quot;y&quot;,&quot;x1&quot;,&quot;x2&quot;,&quot;Noise&quot;)
inSample &lt;- xreg[1:80,]
outSample &lt;- xreg[-c(1:80),]

ourModel &lt;- alm(y~x1+x2, inSample, distribution=&quot;laplace&quot;)
summary(ourModel)</pre>
<p>Here’s the output of the summary:</p>
<pre>Distribution used in the estimation: Laplace
Coefficients:
            Estimate Std. Error Lower 2.5% Upper 97.5%
(Intercept) 95.85207    0.36746   95.12022    96.58392
x1           0.59618    0.02479    0.54681     0.64554
x2          -0.67865    0.00622   -0.69103    -0.66626
ICs:
     AIC     AICc      BIC     BICc 
474.2453 474.7786 483.7734 484.9419</pre>

<p>And here’s the respective plot of the forecast:</p>
<pre>plot(forecast(ourModel,outSample))</pre>
<p><div class="wp-caption"><a rel="nofollow" href="https://i2.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/almLaplaceExample.png?ssl=1" target="_blank"><img sizes="(max-width: 300px) 100vw, 300px" src="https://i2.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/almLaplaceExample-300x175.png?resize=300%2C175&ssl=1" alt="" width="300" class="wp-image-1787" data-recalc-dims="1" height="175"></a>
<p class="wp-caption-text">Forecast from lm with Laplace distribution</p>
</div>
<p>The thing that is currently missing in the function is prediction intervals, but this will be added in the upcoming releases.</p>
<p>Having the likelihood approach, allows comparing different models with different distributions using information criteria. Here’s, for example, what model we get if we assume S-distribution (which has fatter tails than Laplace):</p>
<pre>summary(alm(y~x1+x2, inSample, distribution=&quot;s&quot;))</pre>
<pre>Distribution used in the estimation: S
Coefficients:
            Estimate Std. Error Lower 2.5% Upper 97.5%
(Intercept) 95.61244    0.23386   95.14666    96.07821
x1           0.56144    0.00721    0.54708     0.57581
x2          -0.66867    0.00302   -0.67470    -0.66265
ICs:
     AIC     AICc      BIC     BICc 
482.9358 483.4692 492.4639 493.6325</pre>
<p>As you see, the information criteria for S distribution are higher than for Laplace, so we can conclude that the previous model was better than the second in terms of ICs.</p>
<p><strong>Note</strong> that at this moment the AICc and BICc are not correct for non-normal models (at least the derivation of them needs to be double checked, which I haven’t done yet), so don’t rely on them too much.</p>
<p>I intent to add several other distributions that either are not available in R or are implemented unsatisfactory (from my point of view) – the function is written in a quite flexible way, so this should not be difficult to do. If you have any preferences, please add them on github, <a rel="nofollow" href="https://github.com/config-i1/greybox/issues/13" target="_blank">here</a>.</p>
<p>I also want to implement the mixture distributions, so that things discussed in <a rel="nofollow" href="https://forecasting.svetunkov.ru/en/2017/11/07/multiplicative-state-space-models-for-intermittent-time-series/" target="_blank">the paper on intermittent state-space model</a> can also be implemented using pure regression.</p>
<p>Finally, now that I have alm, we can select between the regression models with different distributions (with </p>
<pre>stepwise()</pre>
<p> function) or even combine them using AIC weights (hello, </p>
<pre>lmCombine()</pre>
<p>!). Yes, I know that it sounds crazy (think of the pool of models in this case), but this should be fun!</p>
<h3>Regression for Multiple Comparison</h3>
<p>One of the reasons why I needed </p>
<pre>alm()</pre>
<p> is because I wanted to implement a parametric analogue of </p>
<pre>nemenyi()</pre>
<p> test. This should be handy, when you have large samples (because the asymptotic properties start kicking in) and it should be more powerful than non-parametric tests. So, long story short, I implemented the test in </p>
<pre>rmc()</pre>
<p> function. It allows, for example, comparing different forecasting methods based on some known error measures. For example, if we can assume that the forecast error is distributed normally, then we can use this test in order to see if the forecasts are biased or not and how they compare with each other. On large samples it should be enough at least to have symmetric distribution.</p>
<p>The test is based on the simple linear model with dummy variables for each provided method (1 if the error corresponds to the method and 0 otherwise). Here’s an example of how this thing works:</p>
<pre>ourData &lt;- cbind(rnorm(100,0,10), rnorm(100,-2,5), rnorm(100,2,6), rlaplace(100,1,5))
colnames(ourData) &lt;- c(&quot;Method A&quot;,&quot;Method B&quot;,&quot;Method C&quot;,&quot;Method D&quot;)

rmc(ourData, distribution=&quot;norm&quot;, level=0.95)</pre>
<p><p>By default the function produces graph in the MCB (Multiple Comparison with the Best) style:</p>
<div class="wp-caption"><a rel="nofollow" href="https://i2.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/rmcExampleNorm.png?ssl=1" target="_blank"><img sizes="(max-width: 300px) 100vw, 300px" src="https://i1.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/rmcExampleNorm-300x175.png?resize=300%2C175&ssl=1" alt="" width="300" class="wp-image-1791" data-recalc-dims="1" height="175"></a>
<p class="wp-caption-text">RMC example with normal distribution</p>
</div>
<p>Don’t forget that we deal here with errors, so closer to zero is better (Method A). Judging by the graph, we can conclude that the least biased method is A, and that the method B is negatively biased (the bounds lie below zero).</p>
<p>If we compare the results of the test with the mean values, we will see similar ranking:</p>
<pre>apply(ourData,2,mean)</pre>

<pre>Method A   Method B   Method C   Method D 
 -0.5523353 -2.4700249  1.8401501  1.1383028</pre>
<p><p>This also reflects how the data was generated. Notice that Method D was generated from Laplace distribution with mean 1, but the test managed to give the correct answer in this situation, because Laplace distribution is symmetric and the sample size is large enough. But the main point of the test is that we can get the confidence intervals for each parameter, so we can see if the differences between the methods are significant: if the intervals intersect, then they are not.</p>
<p>The function also calculates the relative importance of the model (based on AIC), comparing it with the simple model of the type y = const + e. In my example, the constructed model explains the variance quite well, so the importance is 0.999</p>
<p>Similarly, we can compare absolute errors, but we would need to use Folded Normal distribution in this case:</p>
<pre>rmc(abs(ourData), distribution=&quot;fnorm&quot;, level=0.95)</pre>
<p><div class="wp-caption"><a rel="nofollow" href="https://i2.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/rmcExampleFNorm.png?ssl=1" target="_blank"><img sizes="(max-width: 300px) 100vw, 300px" src="https://i1.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/rmcExampleFNorm-300x175.png?resize=300%2C175&ssl=1" alt="" width="300" class="wp-image-1794" data-recalc-dims="1" height="175"></a>
<p class="wp-caption-text">RMC example with folded normal distribution</p>
</div>
<p>Once again, just to compare and make sense of the plot, here are the mean absolute values for the methods:</p>
<pre>apply(abs(ourData),2,mean)</pre>

<pre>Method A   Method B   Method C   Method D 
  7.703774   4.569836   4.840496   5.569483</pre>
<p><p>In order to compare the different methods using squared errors, we may use Chi-Squared distribution, but some preliminary transformations are in order. Specifically, there is only the basic Chi-Squared distribution implemented in RMC, so we need to centre our data (noncentral Chi-Squared distribution would be the next thing to implement):</p>
<pre>newData &lt;- ourData - matrix(apply(ourData,2,mean),nrow(ourData),ncol(ourData),byrow=TRUE)</pre>
<p><p>This way we remove bias and in a way do the estimation of accuracy only. In order to get to the conventional Chi-Squared distribution, we would also need to standardise the data (divide by standard deviation), but this will make everything similar (in our case – several Chi-Squared distributions with h=1 degree of freedom). Instead I would suggest not to do that, and use the new values as is, implying that each method has a Chi-Squared distribution with \(h \times \sigma^2 \) degrees of freedom, where \(h\) is the forecasting horizon and \(\sigma^2\) is the variance of the error of the method. Then the differences in \(\sigma^2\) between the methods will be captured by the test:</p>
<pre>rmc(newData^2, distribution=&quot;chisq&quot;, level=0.95)</pre>
<p><div class="wp-caption"><a rel="nofollow" href="https://i1.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/rmcExampleChisq.png?ssl=1" target="_blank"><img sizes="(max-width: 300px) 100vw, 300px" src="https://i2.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/rmcExampleChisq-300x175.png?resize=300%2C175&ssl=1" alt="" width="300" class="wp-image-1813" data-recalc-dims="1" height="175"></a>
<p class="wp-caption-text">RMC example with Chi-Squared distribution</p>
</div>
<p>You may notice that we get the results very similar to the ones with the Folded Normal distribution. This is because the data is well behaved and does not violate the assumptions of normality.</p>
<p>Overall, in this example we see that although Method A is the least biased of all, it is also considered as the least accurate one due to the higher variance (the parameter for it is significantly higher than for the other methods).</p>
<p>We can also produce plots with vertical lines, that connect the models that are in the same group (no statistical difference, intersection of respective intervals). Here’s the example for the normal distribution:</p>
<pre>rmc(ourData, distribution=&quot;norm&quot;, level=0.95, style=&quot;lines&quot;)</pre>
<p><div class="wp-caption"><a rel="nofollow" href="https://i1.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/rmcExampleNormLines.png?ssl=1" target="_blank"><img sizes="(max-width: 300px) 100vw, 300px" src="https://i0.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/rmcExampleNormLines-300x175.png?resize=300%2C175&ssl=1" alt="" width="300" class="wp-image-1797" data-recalc-dims="1" height="175"></a>
<p class="wp-caption-text">RMC example with normal distribution, lines plot</p>
</div>
<p>If you want to tune the plot, you can always save the model and then replot it with the desired parameters:</p>
<pre>ourModel &lt;- rmc(ourData, distribution=&quot;norm&quot;, level=0.95)
plot(ourModel, xlab=&quot;Models&quot;, ylab=&quot;Errors&quot;)</pre>
<p><p>Finally, it is worth pointing out that the regression on ranks on large samples will give the results similar to the nemenyi test:</p>
<pre>rmc(t(apply(abs(ourData),1,rank)), distribution=&quot;norm&quot;, level=0.95)</pre>
<p><div class="wp-caption"><a rel="nofollow" href="https://i2.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/rmcExampleNormRanks.png?ssl=1" target="_blank"><img sizes="(max-width: 300px) 100vw, 300px" src="https://i1.wp.com/forecasting.svetunkov.ru/wp-content/uploads/2018/08/rmcExampleNormRanks-300x175.png?resize=300%2C175&ssl=1" alt="" width="300" class="wp-image-1799" data-recalc-dims="1" height="175"></a>
<p class="wp-caption-text">RMC example with normal distribution, ranked data</p>
</div>
<p>This function is currently under development, and I decided to show it here just to give an impression of the ongoing research. I still don’t have good solutions for many popular error measures, but the idea right now is to develop the respective elements of the function.</p>
<p>Also, keep in mind that at the moment the function can only be applied to the errors for a specific horizon (for example, for 5 steps ahead). If you feed it the averaged multiple steps ahead forecast errors (e.g. MAE or MSE), you might get unreliable results, because the respective distribution might be complicated and poorly approximated by the existing ones. Having said that, the larger the dataset that you use is, the more consistent results you will get even with Normal distribution (say hi to the Central Limit Theorem).</p>
<h3>What else?</h3>
<p>Several methods have been moved from </p>
<pre>smooth</pre>
<p> to </p>
<pre>greybox</pre>
<p>. These include:</p>
<ul>
<li>pointLik() – returns point Likelihoods, discussed in <a rel="nofollow" href="http://kourentzes.com/forecasting/2018/06/20/isf2018-presentation-beyond-summary-performance-metrics-for-forecast-selection-and-combination/" target="_blank">our research with Nikos</a>;</li>
<li>pAIC, pBIC, pAICc, pBICc – point values of respective information criteria, from <a rel="nofollow" href="http://kourentzes.com/forecasting/2018/06/20/isf2018-presentation-beyond-summary-performance-metrics-for-forecast-selection-and-combination/" target="_blank">the same research</a>;</li>
<li>nParam() – returns number of the estimated parameters in the model (+ variance);</li>
<li>errorType() – returns the type of error used in the model (Additive / Multiplicative);</li>
</ul>
<p>Furthermore, as you might have already noticed, I’ve implemented several distribution functions:</p>
<ul>
<li>Folded normal distribution;</li>
<li>Laplace distribution;</li>
<li>S distribution.</li>
</ul>
<p>Finally, there is also a function, called </p>
<pre>lmDynamic()</pre>
<p>, which uses pAIC in order to produce dynamic linear regression models. But this should be discussed separately in a separate post.</p>
<p>That’s it for now. See you in greybox 0.4.0!</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="https://forecasting.svetunkov.ru/en/2018/08/07/greybox-0-3-0-whats-new/"> R – Modern Forecasting</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 rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/BZv26jdmq38" height="1"></a></div>