Using ConsReg Package

Josep Puig

2020-04-03

Introduction

One of my main problems in my day-to-day work developing predictive models is in the logic of the model. When I talk about the model logic I refer that: the coefficients of the predictive variables make sense from an economic reasoning, the weight of each variable introduced is reasonable, the variables requested by business department have been incorporated,… and of course the error of the model is low!

Imagine the situation, you build a model with Decision Trees for the prediction of credit default: you get low train and test errors, you detect the variables that have more influence and these are from an economic logic point of view (income, savings, loan information, educational level,…), you prepare a shiny app so that the client can “play” with the model and can introduce some random values to the predictive variables,… everything is perfect until after a few days he calls you saying that he doesn’t understand why the model returns a much higher probability of default for a person with a doctorate degree than a person with a lower level of studies.

It could also happen, for example, that by slightly modifying the income level, the probability changes considerably, while by modifying the savings, it changes practically nothing.

These facts can be influenced by how the data are generated (risk policies followed by the company), for not having enought observations or by the structure of the data itself.

The methods that are usually used to correct are, for example, elimination directly of variables from the model that do not have the desired effect, processing of variables (combination of varialbes), regularization, Bayesian methods incorporating priors in the parameters (complicated, if you are not a Bayesian expert), … However, some of these methods do not solve your problem or you spend lots of time try to solve that.

Even if the model error is low, for a person not used to using predictive models, it is really difficult to trust a model that he does not understand and that the results it proposes are illogical. Sometimes, it can be convenient to sacrifice some error if you win in logic. And if the model is logical, it is more likely to generalize well in out-of-sample.

It is in the world of time series, in my experience, where this problem is most accentuated.

ConsReg is a collection of functions, that I have been writting and that I’ve put them together in this package, in order to estimate models with a priori constraints.

It has two main functions:

For the estimation of the parameters you can use functions from several specialized optimization packages. Even with MCMC simulation method optimization. I think the package is quite simple to use and intuitive. It is my first package I write in R and it is the first version of the package, so there will be many improvements! Please contact to me.

Start

This first vignette is a first presentation of the package.

For the first example, I will use the fake_data dataset. This is a fake dataset only to show the functionalities of this package

require(ConsReg)
require(ggplot2)
require(data.table)

ConsReg

data("fake_data")

Let’s start with a simple linear regression:

fit1 = ConsReg(formula = y~x1+x2+x3+ I(x3^2) + x4, family = 'gaussian',
                     optimizer = 'solnp',
                     data = fake_data)
## 
## Iter: 1 fn: 3.6337    Pars:  -1.651009  0.129054 -0.004316 -0.127660  0.020008  0.654738
## solnp--> Completed in 1 iterations
fit1$coefficients
##  (Intercept)           x1           x2           x3      I(x3^2) 
## -1.651008854  0.129054342 -0.004316285 -0.127659853  0.020007889 
##           x4 
##  0.654737554
coef(lm(y~x1+x2+x3+I(x3^2) + x4, data = fake_data))
##  (Intercept)           x1           x2           x3      I(x3^2) 
## -1.651008854  0.129054342 -0.004316285 -0.127659853  0.020007889 
##           x4 
##  0.654737554

The use of the function ConsReg is very simple and similar to glm/lm function:

Possible optimizers are:

Additional arguments of each function can be passed to the function ConsReg.

The object fit1 has the following information:

Error metrics:

fit1$metrics
LogLik RMSE MAE MAPE MSE SMAPE
-2410.638 2.695813 2.052224 4.255439 7.267408 1.086447

Residual analysis

forecast::gghistogram(fit1$residuals, add.normal = T, add.rug = T) + 
  theme_minimal()
## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Registered S3 methods overwritten by 'forecast':
##   method             from    
##   fitted.fracdiff    fracdiff
##   residuals.fracdiff fracdiff

Let’s put some constraints to the model:

It can be easily incoporate in the function:

fit2 = ConsReg(formula = y~x1+x2+x3+ I(x3^2) + x4, data = fake_data,
            family = 'gaussian',
            constraints = '(x3 + `I(x3^2)`) > .01, x4 < .2',
            optimizer = 'mcmc',
            LOWER = -1, UPPER = 1,
            ini.pars.coef = c(-.4, .12, -.004, 0.1, 0.1, .15))
## number of accepted runs: 893 out of 1000 (89.3%)

To put in the function is just:

Observe that now, all coefficient fulfill our constraints:

rbind(coef(fit1), 
      coef(fit2))
##      (Intercept)        x1           x2         x3     I(x3^2)        x4
## [1,]  -1.6510089 0.1290543 -0.004316285 -0.1276599  0.02000789 0.6547376
## [2,]  -0.9731923 0.3265623 -0.011224336  0.3108502 -0.09635921 0.1972126

Also we can compare the errors to see that there is no much difference:

rbind(fit1$metrics, 
      fit2$metrics)
LogLik RMSE MAE MAPE MSE SMAPE
-2410.638 2.695813 2.052224 4.255439 7.267408 1.086447
-2494.864 2.932706 2.152670 3.303897 8.600763 1.228690

For predictions, it follows the same system as a glm or lm object:

pred = data.frame(
  fit1 = predict(fit1, newdata = fake_data[2:3,]), 
  fit2 = predict(fit2, newdata = fake_data[2:3,])
  )
pred
fit1 fit2
2 -1.676538 -0.7073757
3 -2.263643 -1.6110852

Setting the parameter component = T, returns a matrix for the weight of each variable to the predictions

pr = predict(fit2, components = T, newdata = fake_data[5,])
pr
##        Total (Intercept)        x1         x2        x3    I(x3^2)
## 5 -0.5119405  -0.9731923 0.4005724 0.01181268 0.9325506 -0.8672329
##            x4
## 5 -0.01645103

ConsRegArima

As I said, it is in the time series models where the problem mentioned above arises.

In this case, a function has been implemented that estimates a regression with the Arima errors.

This functions is quite similar to stats::arima function or in the forecast package, but the restrictions and constraints have been introduced. Also it can be write more friendly by using formula class. Let me show you.

For this example, I will use another fake dataset

data('series')

The objective function has the following trend:

plot(series$y, t='l')

And the data set has 4 predictive variables:

head(series)
y x1 x2 x3 x4
1.4240582 -0.6051647 -0.2025728 0 -1.7669537
2.9434282 0.0531372 -1.9044371 0 -0.3468814
0.0227284 1.1121740 0.7198898 0 1.8846833
2.2191330 1.2894990 -1.7833160 0 0.1728917
2.9136307 -2.1313672 -2.1982599 0 -2.5840346
0.6459775 0.3245861 0.3002453 0 -0.7843327

We will estimate a first arma model (1, 1) with no regressors and no intercept

fit_ts1 = ConsRegArima(y ~ -1, order = c(1, 1), data = series[1:60, ])
## 
## Iter: 1 fn: 0.1970    Pars:   0.97980 -0.72928
## Iter: 2 fn: 0.1970    Pars:   0.97980 -0.72928
## solnp--> Completed in 2 iterations
fit_ts1$coefficients
##        ar1        ma1 
##  0.9798030 -0.7292797
coef(arima(series$y[1:60], order = c(1, 0, 1), include.mean = F, method = 'CSS'))
##        ar1        ma1 
##  0.9798008 -0.7292652

Next I will add some regressors to the model:

fit_ts2 = ConsRegArima(y ~ x1+x2+x3+x4, order = c(1, 1), data = series[1:60,])
## 
## Iter: 1 fn: -0.7311   Pars:   0.82849  0.39639 -0.80844  1.35773 -0.40255 -0.33677  1.16707
## Iter: 2 fn: -0.7311   Pars:   0.82849  0.39639 -0.80844  1.35773 -0.40255 -0.33677  1.16707
## solnp--> Completed in 2 iterations

Next I will add some constraints to the model:

fit_ts3 = ConsRegArima(y ~ x1+x2+x3+x4, order = c(1, 1), data = series[1:60,], 
                       LOWER = -1, UPPER = 1, 
                       constraints = "x4 < x2", 
                       ini.pars.coef = c(.9, .3, -.1, .3, -.3), 
                       control = list(trace = 0) #  not show the trace of the optimizer 
                       )
fit_ts3$coefficients
## (Intercept)          x1          x2          x3          x4         ar1 
##   0.9983209   0.2870705  -0.4862388   0.4804189  -0.4862388  -0.3641870 
##         ma1 
##   0.6373165

To put in the function is just:

Next, we will change the optimizer. Let’s try with a genetic algorithm:

fit_ts4 = ConsRegArima(y ~ x1+x2+x3+x4, order = c(1, 1), data = series[1:60,],
                       LOWER = -1, UPPER = 1,
                       constraints = "x4 < x2",
                       penalty = 10000,
                       optimizer = 'GA', maxiter = 1000,
                       monitor = NULL, #  not show the trace of the optimizer
                       ini.pars.coef = c(.9, .2, 0, .3, -.6)
                       )
fit_ts4$coefficients
## (Intercept)          x1          x2          x3          x4         ar1 
##   0.8739497   0.4262889  -0.5202935   0.9974570  -0.5436331  -0.2806282 
##         ma1 
##   0.9539885

The restrictions are still fulfilled We can compare the errors of the 4 models:

data.frame(
  metrics = colnames(fit_ts1$metrics),
  fit_ts1 = as.numeric(fit_ts1$metrics), 
  fit_ts2 = as.numeric(fit_ts2$metrics), 
  fit_ts3 = as.numeric(fit_ts3$metrics),
  fit_ts4 = as.numeric(fit_ts4$metrics)
  )
metrics fit_ts1 fit_ts2 fit_ts3 fit_ts4
ME 0.1176876 -0.0022704 0.1183575 0.0523606
RMSE 1.2075971 0.4773731 0.7518715 0.5971884
MAE 0.9814655 0.3849611 0.6203518 0.4842597
MPE -179.8131183 -24.5363218 -33.0289775 -52.0671108
MAPE 322.1692219 80.1435850 133.7496066 126.4981995

For predictions you will see that is very easy:

pred = predict(fit_ts4, newdata = series[61:63, ], h=3, intervals = 90)
pred$predict
Prediction Prediction_High Prediction_Low
61 1.720644 2.711354 0.7299348
62 2.208070 3.402294 1.0138447
63 2.895303 4.104100 1.6865051

And this object, you can plot to see the predictions as well as the fitted values:

plot(pred) + theme_minimal()

In the ConsRegArima function, you can introduce the seasonal part P,Q, or in the formula, you can introduce lags in the predictor variables:

fit_ts5 = ConsRegArima(y ~ x1+x3+
                         shift(x3, 2) + # x2 from 2 periods above
                         x4, 
                       order = c(1, 1), data = series[1:60,], 
                       seasonal = list(order = c(1, 0), period = 4), # seasonal component
                       control = list(trace = 0)
                       )

If you have used lags in the predictive variables, then, in the predict function, you must add the original data:

pred = predict(fit_ts5, newdata = series[61:63,], origdata = series[1:60,])
pred$predict
Prediction Prediction_High Prediction_Low
59 2.191149 3.601062 0.7812351
60 2.503279 4.006080 1.0004775
61 2.193245 3.706991 0.6794988

Finally, I have implemented a feature that I miss in many time series packages which is the possibility of backtesting.

For a ConsRegArima object, I have implemented the “rolling” function that allows rolling-forecast with recalibration every \(n\) periods, and projections to \(h\) periods.

And of course, very easy to carry out!

ro = rolling(object = fit_ts3, used.sample = 50, 
             refit = 4, h = 4, orig.data = series)

In this case, the arguments are:

The errors of the rolling are:

And graphically:

plot(ro) + theme_minimal()

We can compare that errors of 4-step-forecasts are greater than 1-step-forecast:

ro$results
xx Real Prediction Prediction_High Prediction_Low Fitted
55 2.2588072 1.3297014 2.702927 -0.0435246 1.3470754
56 0.9998795 0.3919824 1.750749 -0.9667838 0.6849166
57 3.5387356 3.1040633 4.453682 1.7544449 3.0891555
58 2.5365840 1.9089582 3.247440 0.5704760 2.0334505
59 0.1263169 0.2512654 1.612710 -1.1101792 0.3981201
60 2.0548210 1.7922128 3.124912 0.4595141 1.6842149