Greybox main vignette
Ivan Svetunkov
20230915
There are three wellknown notions of “boxes” in modelling: 1. White
box  the model that is completely transparent and does not have any
randomness. One can see how the inputs are transformed into the specific
outputs. 2. Black box  the model which does not have an apparent
structure. One can only observe inputs and outputs but does not know
what happens inside. 3. Grey box  the model that is in between the
first two. We observe inputs and outputs plus have some information
about the structure of the model, but there is still a part of
unknown.
The white boxes are usually used in optimisations (e.g. linear
programming), while black boxes are popular in machine learning. As for
the grey box models, they are more often used in analysis and
forecasting. So the package greybox
contains models that
are used for these purposes.
At the moment the package contains augmented linear model function
and several basic functions that implement model selection and
combinations using information criteria (IC). You won’t find statistical
tests in this package  there’s plenty of them in the other packages.
Here we try using the modern techniques and methods that do not rely on
hypothesis testing. This is the main philosophical point of
greybox
.
Main functions
The package includes the following functions for models
construction:
 alm()  Augmented Linear Model. This is
something similar to GLM, but with a focus on forecasting and the
information criteria usage for time series. It also supports mixture
distribution models for the intermittent data and allows adding trend to
the data via the formula.
 stepwise()  select the linear model with
the lowest IC from all the possible in the provided data. Uses partial
correlations. Works fast;
 lmCombine()  combine the linear models
into one using IC weights;
 lmDynamic()  produce model with dynamic
weights and time varying parameters based on point IC weight.
See discussion of some of these functions in this vignette below.
Models evaluation functions
 ro()  produce forecasts with a specified
function using rolling origin.
measures()
 function, returning a bunch of error
measures for the provided forecast and the holdout sample.
rmcb()
 regression on ranks of forecasting methods.
This is a fast alternative to the classical nemenyi / MCB test.
Methods
The following methods can be applied to the models, produced by
alm()
, stepwise()
, lmCombine()
and lmDynamic()
:
logLik()
 extracts loglikelihood.
AIC()
, AICc()
, BIC()
,
BICc()
 calculates the respective information
criteria.
pointLik()
 extracts the point likelihood.
pAIC()
, pAICc()
, pBIC()
,
pBICc()
 calculates the respective point information
criteria, based on pointLik.
actuals()
 extracts the actual values of the response
variable.
coefbootstrap()
 produces bootstrapped values of
parameters, taking nsim
samples of the size
size
from the data and reapplying the model.
coef()
, coefficients()
 extract the
parameters of the model.
confint()
 extracts the confidence intervals for the
parameters.
vcov()
 extracts the variancecovariance matrix of the
parameters.
sigma()
 extracts the standard deviation of the
residuals.
nobs()
 the number of the insample observations of
the model.
nparam()
 the number of all the estimated parameters
in the model.
nvariate()
 the number of variates (columns /
dimensions) of the resposne variable.
summary()
 produces the summary of the model.
predict()
 produces the predictions based on the model
and the provided newdata
. If the newdata
is
not provided, then it uses the already available data in the model. Can
also produce confidence
and prediction
intervals.
forecast()
 acts similarly to predict()
with few differences. It has a parameter h
 forecast
horizon  which is NULL
by default and is set to be equal
to the number of rows in newdata
. However, if the
newdata
is not provided, then it will produce forecasts of
the explanatory variables to the horizon h
and use them as
newdata
. Finally, if h
and
newdata
are provided, then the number of rows to use will
be regulated by h
.
plot()
 produces several plots for the analysis of the
residuals. This includes: Fitted over time, Standardised residuals vs
Fitted, Absolute residuals vs Fitted, QQ plot with the specified
distribution, Squared residuals vs Fitted, ACF of the residuals and PACF
of the residuals, which is regulated by which
parameter.
See documentation for more info: ?plot.greybox
.
detectdst()
and detectleap()
 methods
that return the ids of the hour / date for the DST / Leap year
change.
extract()
method, needed in order to produce printable
regression outputs using texreg()
function from the
texreg
package.
Distribution functions
qlaplace()
, dlaplace()
,
rlaplace()
, plaplace()
 functions for Laplace
distribution.
qalaplace()
, dalaplace()
,
ralaplace()
, palaplace()
 functions for
Asymmetric Laplace distribution.
qs()
, ds()
, rs()
,
ps()
 functions for S distribution.
qgnorm()
, dgnorm()
, rgnorm()
,
pgnorm()
 functions for the Generalised normal
distribution.
qfnorm()
, dfnorm()
, rfnorm()
,
pfnorm()
 functions for folded normal distribution.
qtplnorm()
, dtplnorm()
,
rtplnorm()
, ptplnorm()
 functions for three
parameter log normal distribution.
qbcnorm()
, dbcnorm()
,
rbcnorm()
, pbcnorm()
 functions for the
BoxCox normal distribution.
qlogitnorm()
, dlogitnorm()
,
rlogitnorm()
, plogitnorm()
 functions for the
Logitnormal distribution.
Additional functions
graphmaker()
 produces linear plots for the variable,
its forecasts and fitted values.
xregExpander
The function xregExpander()
is useful in cases when the
exogenous variable may influence the response variable either via some
lags or leads. As an example, consider BJsales.lead
series
from the datasets
package. Let’s assume that the
BJsales
variable is driven by the today’s value of the
indicator, the value five and 10 days ago. This means that we need to
produce lags of BJsales.lead
. This can be done using
xregExpander()
:
BJxreg < xregExpander(BJsales.lead,lags=c(5,10))
The BJxreg
is a matrix, which contains the original
data, the data with the lag 5 and the data with the lag 10. However, if
we just move the original data several observations ahead or backwards,
we will have missing values in the beginning / end of series, so
xregExpander()
fills in those values with the forecasts
using es()
and iss()
functions from
smooth
package (depending on the type of variable we are
dealing with). This also means that in cases of binary variables you may
have weird averaged values as forecasts (e.g. 0.7812), so beware and
look at the produced matrix. Maybe in your case it makes sense to just
substitute these weird numbers with zeroes…
You may also need leads instead of lags. This is regulated with the
same lags
parameter but with positive values:
BJxreg < xregExpander(BJsales.lead,lags=c(7,5,10))
Once again, the values are shifted, and now the first 7 values are
backcasted. In order to simplify things we can produce all the values
from 10 lags till 10 leads, which returns the matrix with 21
variables:
BJxreg < xregExpander(BJsales.lead,lags=c(10:10))
stepwise
The function stepwise() does the selection based on an information
criterion (specified by user) and partial correlations. In order to run
this function the response variable needs to be in the first column of
the provided matrix. The idea of the function is simple, it works
iteratively the following way:
 The basic model of the first variable and the constant is
constructed (this corresponds to simple mean). An information criterion
is calculated;
 The correlations of the residuals of the model with all the original
exogenous variables are calculated;
 The regression model of the response variable and all the variables
in the previous model plus the new most correlated variable from (2) is
constructed using
lm()
function;
 An information criterion is calculated and is compared with the one
from the previous model. If it is greater or equal to the previous one,
then we stop and use the previous model. Otherwise we go to step 2.
This way we do not do a blind search, going forward or backwards, but
we follow some sort of “trace” of a good model: if the residuals contain
a significant part of variance that can be explained by one of the
exogenous variables, then that variable is included in the model.
Following partial correlations makes sure that we include only
meaningful (from technical point of view) variables in the model. In
general the function guarantees that you will have the model with the
lowest information criterion. However this does not guarantee that you
will end up with a meaningful model or with a model that produces the
most accurate forecasts. So analyse what you get as a result.
Let’s see how the function works with the BoxJenkins data. First we
expand the data and form the matrix with all the variables:
BJxreg < as.data.frame(xregExpander(BJsales.lead,lags=c(10:10)))
BJxreg < cbind(as.matrix(BJsales),BJxreg)
colnames(BJxreg)[1] < "y"
ourModel < stepwise(BJxreg)
This way we have a nice data frame with nice names, not something
weird with strange long names. It is important to note that the response
variable should be in the first column of the resulting matrix. After
that we use stepwise function:
ourModel < stepwise(BJxreg)
And here’s what it returns (the object of class lm
):
ourModel
#> Time elapsed: 0.27 seconds
#> Call:
#> alm(formula = y ~ xLag4 + xLag9 + xLag3 + xLag10 + xLag5 + xLag6 +
#> xLead9 + xLag7 + xLag8, data = data, distribution = "dnorm")
#>
#> Coefficients:
#> (Intercept) xLag4 xLag9 xLag3 xLag10 xLag5
#> 17.6668913 3.3703773 1.3725674 4.6768101 1.5415708 2.3209845
#> xLag6 xLead9 xLag7 xLag8
#> 1.7074084 0.3765307 1.4027852 1.3371727
The values in the function are listed in the order of most correlated
with the response variable to the least correlated ones. The function
works very fast because it does not need to go through all the variables
and their combinations in the dataset.
All the basic methods can be used together with the final model
(e.g. predict()
, forecast()
,
summary()
etc).
Furthermore, the greybox
package implements
extract()
method from texreg
package for the
production of printable outputs from the regression, here is an
example:
texreg::htmlreg(ourModel)
Statistical models

Model 1

(Intercept)

17.67^{*}


[16.07; 19.26]

xLag4

3.37^{*}


[ 2.75; 3.99]

xLag9

1.37^{*}


[ 0.75; 1.99]

xLag3

4.68^{*}


[ 4.10; 5.25]

xLag10

1.54^{*}


[ 0.98; 2.11]

xLag5

2.32^{*}


[ 1.68; 2.96]

xLag6

1.71^{*}


[ 1.06; 2.35]

xLead9

0.38^{*}


[ 0.12; 0.63]

xLag7

1.40^{*}


[ 0.76; 2.05]

xLag8

1.34^{*}


[ 0.69; 1.98]

Num. obs.

150.00

Num. param.

11.00

Num. df

139.00

AIC

416.45

AICc

418.36

BIC

449.57

BICc

454.36

^{*} 0 outside the confidence interval.

Similarly, you can produce pdf tables via texreg()
function from that package. Alternatively, you can use
kable()
function from knitr
package on the
summary to get a table for LaTeX / HTML.
lmCombine
lmCombine()
function creates a pool of linear models
using lm()
, writes down the parameters, standard errors and
information criteria and then combines the models using IC weights. The
resulting model is of the class “lm.combined”. The speed of the function
deteriorates exponentially with the increase of the number of variables
\(k\) in the dataset, because the
number of combined models is equal to \(2^k\). The advanced mechanism that uses
stepwise()
and removes a large chunk of redundant models is
also implemented in the function and can be switched using
bruteforce
parameter.
Here’s an example of the reduced data with combined model and the
parameter bruteforce=TRUE
:
ourModel < lmCombine(BJxreg[,c(3:7,18:22)],bruteforce=TRUE)
summary(ourModel)
#> The AICc combined model
#> Response variable: y
#> Distribution used in the estimation: Normal
#> Coefficients:
#> Estimate Std. Error Importance Lower 2.5% Upper 97.5%
#> (Intercept) 20.9213 0.2326 1.0000 20.4616 21.3811 *
#> x 0.0432 0.0286 0.2591 0.0998 0.0134
#> xLag5 6.3977 0.0839 1.0000 6.2318 6.5636 *
#> xLag4 5.8466 0.0899 1.0000 5.6688 6.0244 *
#> xLag3 5.6854 0.0901 1.0000 5.5074 5.8635 *
#> xLag2 0.1251 0.0382 0.2876 0.0496 0.2006 *
#> xLag1 0.0844 0.0342 0.2717 0.1521 0.0167 *
#> xLead1 0.0908 0.0324 0.2782 0.1547 0.0268 *
#> xLead2 0.0356 0.0257 0.2600 0.0865 0.0152
#> xLead3 0.1200 0.0343 0.2971 0.1877 0.0523 *
#> xLead4 0.0068 0.0229 0.2585 0.0520 0.0384
#> xLead5 0.1159 0.0300 0.3030 0.0566 0.1751 *
#>
#> Error standard deviation: 2.2075
#> Sample size: 150
#> Number of estimated parameters: 7.2151
#> Number of degrees of freedom: 142.7849
#> Approximate combined information criteria:
#> AIC AICc BIC BICc
#> 670.6603 671.4964 692.3824 694.4770
summary()
function provides the table with the
parameters, their standard errors, their relative importance and the 95%
confidence intervals. Relative importance indicates in how many cases
the variable was included in the model with high weight. So, in the
example above variables xLag5, xLag4, xLag3 were included in the models
with the highest weights, while all the others were in the models with
lower ones. This may indicate that only these variables are needed for
the purposes of analysis and forecasting.
The more realistic situation is when the number of variables is high.
In the following example we use the data with 21 variables. So if we use
brute force and estimate every model in the dataset, we will end up with
\(2^{21}\) = 2^21
combinations of models, which is not possible to estimate in the
adequate time. That is why we use bruteforce=FALSE
:
ourModel < lmCombine(BJxreg,bruteforce=FALSE)
summary(ourModel)
#> The AICc combined model
#> Response variable: y
#> Distribution used in the estimation: Normal
#> Coefficients:
#> Estimate Std. Error Importance Lower 2.5% Upper 97.5%
#> (Intercept) 17.6924 0.7760 1.0000 16.1580 19.2268 *
#> xLag4 3.3747 0.3020 1.0000 2.7775 3.9719 *
#> xLag9 1.3711 0.3029 0.9998 0.7723 1.9699 *
#> xLag3 4.6846 0.2809 1.0000 4.1292 5.2400 *
#> xLag10 1.5424 0.2749 1.0000 0.9989 2.0859 *
#> xLag5 2.3222 0.3118 1.0000 1.7058 2.9386 *
#> xLag6 1.7075 0.3145 1.0000 1.0858 2.3293 *
#> xLead9 0.3638 0.1247 0.9663 0.1174 0.6103 *
#> xLag7 1.4017 0.3152 0.9997 0.7786 2.0249 *
#> xLag8 1.3363 0.3132 0.9994 0.7170 1.9556 *
#>
#> Error standard deviation: 0.936
#> Sample size: 150
#> Number of estimated parameters: 10.9652
#> Number of degrees of freedom: 139.0348
#> Approximate combined information criteria:
#> AIC AICc BIC BICc
#> 416.7030 418.6040 449.7152 454.4777
In this case first, the stepwise()
function is used,
which finds the best model in the pool. Then each variable that is not
in the model is added to the model and then removed iteratively. IC,
parameters values and standard errors are all written down for each of
these expanded models. Finally, in a similar manner each variable is
removed from the optimal model and then added back. As a result the pool
of combined models becomes much smaller than it could be in case of the
brute force, but it contains only meaningful models, that are close to
the optimal. The rationale for this is that the marginal contribution of
variables deteriorates with the increase of the number of parameters in
case of the stepwise function, and the IC weights become close to each
other around the optimal model. So, whenever the models are combined,
there is a lot of redundant models with very low weights. By using the
mechanism described above we remove those redundant models.
There are several methods for the lm.combined
class,
including:
predict.greybox()
 returns the point and interval
predictions.
forecast.greybox()
 wrapper around
predict()
The forecast horizon is defined by the length of
the provided sample of newdata
.
plot.lm.combined()
 plots actuals and fitted
values.
plot.predict.greybox()
 which uses
graphmaker()
function from smooth
in order to
produce graphs of actuals and forecasts.
As an example, let’s split the whole sample with BoxJenkins data
into insample and the holdout:
BJInsample < BJxreg[1:130,];
BJHoldout < BJxreg[(1:130),];
ourModel < lmCombine(BJInsample,bruteforce=FALSE)
A summary and a plot of the model:
summary(ourModel)
#> The AICc combined model
#> Response variable: y
#> Distribution used in the estimation: Normal
#> Coefficients:
#> Estimate Std. Error Importance Lower 2.5% Upper 97.5%
#> (Intercept) 19.4134 0.8558 1.0000 17.7189 21.1079 *
#> xLag4 3.3480 0.2966 1.0000 2.7607 3.9353 *
#> xLag9 1.3340 0.2983 0.9990 0.7434 1.9247 *
#> xLag3 4.7559 0.2787 1.0000 4.2040 5.3078 *
#> xLag10 1.5364 0.2701 1.0000 1.0016 2.0713 *
#> xLag5 2.3209 0.3063 1.0000 1.7144 2.9274 *
#> xLag6 1.6611 0.3090 1.0000 1.0493 2.2730 *
#> xLead9 0.2948 0.1260 0.8920 0.0454 0.5443 *
#> xLag8 1.3692 0.3084 0.9989 0.7585 1.9799 *
#> xLag7 1.3274 0.3093 0.9982 0.7150 1.9397 *
#>
#> Error standard deviation: 0.9541
#> Sample size: 130
#> Number of estimated parameters: 10.8881
#> Number of degrees of freedom: 119.1119
#> Approximate combined information criteria:
#> AIC AICc BIC BICc
#> 367.7981 369.9900 399.0203 404.3546
plot(ourModel)
Importance tells us how important the respective variable is in the
combination. 1 means 100% important, 0 means not important at all.
And the forecast using the holdout sample:
ourForecast < predict(ourModel,BJHoldout)
plot(ourForecast)
These are the main functions implemented in the package for now. If
you want to read more about IC model selection and combinations, I would
recommend (Burnham and Anderson 2004)
textbook.
References
Burnham, Kenneth P, and David R Anderson. 2004.
Model Selection and Multimodel Inference.
Edited by Kenneth P Burnham and David R Anderson. Springer New York.
https://doi.org/10.1007/b97636.