I created an R package called comparer that makes code comparisons simple. It is available on CRAN CRAN, or on GitHub. The vignette is shown below.

———–

Introduction to the comparer R package

Collin Erickson

2017-10-30

When coding, especially for data science, there are multiple ways to solve each problem. When presented with two options, you want to pick the one that is faster and/or more accurate. Comparing different code chunks on the same task can be tedious. It often requires creating data, writing a for loop (or using sapply), then comparing.

The comparer package makes this comparison quick and simple:

The same data can be given in to each model.
Various metrics can be used to judge results, including using the predicted errors from the code.
The results are displayed in a table that allows you to quickly judge the results.

This document introduces the main function of the comparer package, mbc.

Motivation from `microbenchmark`

The R package microbenchmark provides the fantastic eponymous function. It makes it simple to run different segments of code and see which is faster. Borrowing an example from http://adv-r.had.co.nz/Performance.html, the following shows how it gives a summary of how fast each ran.

library(microbenchmark)
x <- runif(100)
microbenchmark(sqrt(x), x ^ .5)

## Unit: nanoseconds
##     expr  min   lq    mean median   uq   max neval
##  sqrt(x)  380 1140 1418.03 1141.0 1521 18246   100
##    x^0.5 4181 6272 6880.24 7222.5 7603  9123   100

However it gives no summary of the output. For this example it is fine since the output is deterministic, but when working with randomness or model predictions we want to have some sort of summary or evaluation metric to see which has better accuracy, or to just see how the outputs differ.

`mbc` to the rescue

The function mbc in the comparer package was created to solve this problem, where a comparison of the output is desired in addition to the run time.

For example, we may wish to see how the sample size affects an estimate of the mean of a random sample. The following shows the results of finding the mean of 10 and 100 samples from a normal distribution.

library(comparer)
mbc(mean(rnorm(10)), mean(rnorm(100)))

## Run times (sec)
##           Function Sort1 Sort2 Sort3 Sort4 Sort5 mean sd neval
## 1  mean(rnorm(10))     0     0     0     0     0    0  0     5
## 2 mean(rnorm(100))     0     0     0     0     0    0  0     5
## 
## Output summary
##               Func Stat       Sort1       Sort2       Sort3       Sort4
## 1  mean(rnorm(10))    1 -0.25552091 -0.03541135 -0.03160087 -0.01748239
## 2 mean(rnorm(100))    1 -0.09895678 -0.01811174 -0.01659186  0.04176382
##        Sort5         mean       sd
## 1 0.23999128 -0.020004846 0.175673
## 2 0.07945149 -0.002489014 0.067863

By default it only runs 5 trials, but this can be changed with the times parameter. The first part of the output gives the run times. For 5 or fewer, it shows all the values in sorted order, for more than 5 it shows summary statistics. Unfortunately, the timing is only accurate up to 0.01 seconds, so these all show as 0.

The second section of the output gives the summary of the output. This also will show summary stats for more than 5 trials, but for this small sample size it shows all the values in sorted order with the mean and standard deviation given. The first column shows the name of each, and the second column shows which output statistic is given. Since there is only one output for this code it is called “1”.

Setting times changes the number of trials run. Below the same example as above is run but for 100 trials.

mbc(mean(rnorm(10)), mean(rnorm(100)), times=100)

## Run times (sec)
##           Function Min. 1st Qu. Median Mean 3rd Qu. Max. sd neval
## 1  mean(rnorm(10))    0       0      0    0       0    0  0   100
## 2 mean(rnorm(100))    0       0      0    0       0    0  0   100
## 
## Output summary
##               Func Stat       Min.     1st Qu.      Median        Mean
## 1  mean(rnorm(10))    1 -0.8175108 -0.15620118 -0.02166650 0.001056646
## 2 mean(rnorm(100))    1 -0.2170817 -0.06626509 -0.00111176 0.005715118
##      3rd Qu.      Max.        sd
## 1 0.19105787 0.6266277 0.2954722
## 2 0.06441406 0.2827851 0.1021612

We see that the mean of both is around zero, but that the larger sample size (mean(rnorm(100))) has a tighter distribution and a standard deviation a third as large as the other, which is about what we expect for a sample that is 10 times larger (it should be $\sqrt{10} \approx 3.16$ times smaller on average).

In this example each function had its own input, but many times we want to compare the functions on the same input for better comparison.

Shared input

Input can be passed in to the input argument as a list, and then the code will be evaluated in an environment with that data. In this example we compare the functions mean and median on random data from an exponential distribution. The mean should be about 1, while the median should be about $\ln(2)=0.693$.

mbc(mean(x), median(x), input=list(x=rexp(30)))

## Run times (sec)
##    Function Sort1 Sort2 Sort3 Sort4 Sort5 mean sd neval
## 1   mean(x)     0     0     0     0     0    0  0     5
## 2 median(x)     0     0     0     0     0    0  0     5
## 
## Output summary
##        Func Stat     Sort1     Sort2     Sort3     Sort4     Sort5
## 1   mean(x)    1 0.7360094 0.7360094 0.7360094 0.7360094 0.7360094
## 2 median(x)    1 0.5679456 0.5679456 0.5679456 0.5679456 0.5679456
##        mean sd
## 1 0.7360094  0
## 2 0.5679456  0

In this case each evaluation is identical since the input is not random. The data passed to input is kept as is, so there is no randomness from the data. If we want randomness in the data, we can use inputi, which evaluates its argument as an expression, meaning that each time it will be different.

Below is the same code as above except inputi is used with x set in brackets instead of a list. We see there is randomness and we can get an idea of the distribution of the median and mean.

mbc(mean(x), median(x), inputi={x=rexp(30)})

## Run times (sec)
##    Function Sort1 Sort2 Sort3 Sort4 Sort5 mean sd neval
## 1   mean(x)     0     0     0     0     0    0  0     5
## 2 median(x)     0     0     0     0     0    0  0     5
## 
## Output summary
##        Func Stat        V1        V2        V3       V4        V5
## 1   mean(x)    1 0.9041031 1.1820833 0.9131991 1.483039 0.9399961
## 2 median(x)    1 0.7290917 0.8445266 0.5957285 1.163204 0.6192530
##        mean        sd
## 1 1.0844841 0.2505615
## 2 0.7903606 0.2306833

When the code chunks to evaluate are simple functions of a single variable, this can be simplified. Look how simple it is to run a test on these!

mbc(mean, median, inputi=rexp(30))

## Run times (sec)
##   Function Sort1 Sort2 Sort3 Sort4 Sort5 mean sd neval
## 1     mean     0     0     0     0     0    0  0     5
## 2   median     0     0     0     0     0    0  0     5
## 
## Output summary
##     Func Stat        V1        V2        V3        V4        V5      mean
## 1   mean    1 0.8241886 0.8783592 0.7454609 0.8459507 1.1176282 0.8823175
## 2 median    1 0.6591866 0.4701343 0.5269150 0.4204864 0.7466348 0.5646714
##          sd
## 1 0.1403727
## 2 0.1352737

Comparing to expected values

The previous comparisons showed a summary of the outputs, but many times we want to compare output values to true values, then calculate a summary statistic, such as an average error. The argument target specifies the values the code chunks should give, then summary statistics can be calculated by specifying metrics, which defaults to calculating the rmse.

For example, suppose we have data from a linear function, and want to see how accurate the model is when the output values are corrupted with noise. Below we compare two linear models: the first with an intercept term, and the second without. The model with the intercept term should be much better since the data has an intercept of $-0.6$.

We see that the output is different in a few ways now. The Stat column tells what the row is showing. These all say rmse, meaning they are giving the root mean squared error of the predicted values compared to the true y. There’s also a new section at the bottom title Compare. This compares the rmse values from the two methods, and does a t-test to see if the difference is significant. However, since there is no randomness, it fails to perform the t-test.

n <- 20
x <- seq(0, 1, length.out = n)
y <- 1.8 * x - .6
ynoise <- y + rnorm(n, 0, .2)

mbc(predict(lm(ynoise ~ x), data.frame(x)),
    predict(lm(ynoise ~ x - 1), data.frame(x)),
    target = y)

## Run times (sec)
##                                     Function Sort1 Sort2 Sort3 Sort4 Sort5
## 1     predict(lm(ynoise ~ x), data.frame(x))     0     0     0  0.01  0.02
## 2 predict(lm(ynoise ~ x - 1), data.frame(x))     0     0     0  0.00  0.00
##    mean          sd neval
## 1 0.006 0.008944272     5
## 2 0.000 0.000000000     5
## 
## Output summary
##                                         Func Stat         V1         V2
## 1     predict(lm(ynoise ~ x), data.frame(x)) rmse 0.03031246 0.03031246
## 2 predict(lm(ynoise ~ x - 1), data.frame(x)) rmse 0.31245490 0.31245490
##           V3         V4         V5       mean sd
## 1 0.03031246 0.03031246 0.03031246 0.03031246  0
## 2 0.31245490 0.31245490 0.31245490 0.31245490  0
## 
## Compare
##                                                                                   Func
## 1 predict(lm(ynoise ~ x), data.frame(x)) vs predict(lm(ynoise ~ x - 1), data.frame(x))
##   Stat conf.low conf.up  t  p
## 1 rmse       NA      NA NA NA

To add randomness we can simply define ynoise in the inputi argument, as shown below. Now there is randomness in the data, so a paired t-test can be computed. It is paired since the same ynoise is given to each model. We see that even with only a sample size of 5, the p-value is highly significant.

mbc(predict(lm(ynoise ~ x), data.frame(x)),
    predict(lm(ynoise ~ x - 1), data.frame(x)),
    inputi={ynoise <- y + rnorm(n, 0, .2)},
    target = y)

## Run times (sec)
##                                     Function Sort1 Sort2 Sort3 Sort4 Sort5
## 1     predict(lm(ynoise ~ x), data.frame(x))     0     0     0  0.00  0.00
## 2 predict(lm(ynoise ~ x - 1), data.frame(x))     0     0     0  0.01  0.02
##    mean          sd neval
## 1 0.000 0.000000000     5
## 2 0.006 0.008944272     5
## 
## Output summary
##                                         Func Stat        V1         V2
## 1     predict(lm(ynoise ~ x), data.frame(x)) rmse 0.1209678 0.01310184
## 2 predict(lm(ynoise ~ x - 1), data.frame(x)) rmse 0.3258303 0.31153273
##            V3         V4         V5       mean          sd
## 1 0.009311599 0.01142981 0.04681271 0.04032476 0.047652966
## 2 0.311442443 0.31147247 0.31202430 0.31446045 0.006360358
## 
## Compare
##                                                                                Func
## t predict(lm(ynoise ~ x), data.frame(x))-predict(lm(ynoise ~ x - 1), data.frame(x))
##   Stat         V1         V2         V3         V4         V5       mean
## t rmse -0.2048625 -0.2984309 -0.3021308 -0.3000427 -0.2652116 -0.2741357
##           sd         t            p
## t 0.04160401 -14.73382 0.0001235005

Simplifying with `evaluator`

Many times the code chunks we want to compare only differ by a small amount, such as asingle argument. In the example above, the only difference is the formula in the lm command. With mbc, the evaluator can be set to make these cases easier. The argument for evaluator should be an expression including ., which will be replaced with the code chunks provided. The example below rewrites the above comparison using evaluator.

mbc(ynoise ~ x,
    ynoise ~ x - 1,
    evaluator=predict(lm(.), data.frame(x)),
    inputi={ynoise <- y + rnorm(n, 0, .2)},
    target = y)

## Run times (sec)
##         Function Sort1 Sort2 Sort3 Sort4 Sort5  mean          sd neval
## 1     ynoise ~ x     0     0     0     0  0.01 0.002 0.004472136     5
## 2 ynoise ~ x - 1     0     0     0     0  0.02 0.004 0.008944272     5
## 
## Output summary
##             Func Stat         V1         V2        V3         V4
## 1     ynoise ~ x rmse 0.05038184 0.02971854 0.1055691 0.06501551
## 2 ynoise ~ x - 1 rmse 0.31377557 0.31206988 0.3214883 0.31179433
##           V5       mean          sd
## 1 0.04056214 0.05824943 0.029468982
## 2 0.31132893 0.31409141 0.004236993
## 
## Compare
##                        Func Stat         V1         V2         V3
## t ynoise ~ x-ynoise ~ x - 1 rmse -0.2633937 -0.2823513 -0.2159192
##           V4         V5      mean        sd         t            p
## t -0.2467788 -0.2707668 -0.255842 0.0257803 -22.19059 2.441287e-05

K-Fold Cross Validation

K-fold cross validation can also be done using mbc using the kfold parameter. K-fold cross validation involves splitting $N$ data points into $k$ groups. kfold should specify what this $N$ is, since it depends on the data. By default it will set the number of folds, $k$, to be times. Then each replicate will be evaluating a single fold. Note that this will not do $k$ folds $times$ times.

To make $k$ different from $times$, pass in $kfold$ as a vector whose second element in the number of folds. For example, suppose you have 100 data points, want to do 5 folds, and repeat this process twice (i.e. evaluate 10 folds). Then you should pass in $kfold=c(100,5)$ and $times=10$. The first five trials would then the the five separate folds. The sixth through tenth trials would be a new partition of the data into five folds.

Then to use this folds you must use ki as part of an expression in the code chunk or inputi. The following shows how to use k-fold cross validation fitting a linear model to the cars dataset. Setting kfold=c(nrow(cars), 5) tells it that you want to use 5 folds on the cars data set. It has 50 rows, so in each trial ki is a subset of 1:50 of 40 elements. Setting times=30 means that we are repeating the five folds six times. The code chunk fits the model, makes predictions on the hold-out data, and calculates the RMSE.

mbc({mod <- lm(dist ~ speed, data=cars[ki,])
     p <- predict(mod,cars[-ki,])
     sqrt(mean((p - cars$dist[-ki])^2))
     },
    kfold=c(nrow(cars), 5),
    times=30)

## Run times (sec)
##                                                                                                           Function
## 1 { mod <- lm(dist ~ speed, data = cars[ki, ]) p <- predict(mod, cars[-ki, ]) sqrt(mean((p - cars$dist[-ki])^2)) }
##   Min. 1st Qu. Median        Mean 3rd Qu. Max.          sd neval
## 1    0       0      0 0.003666667       0 0.02 0.007183954    30
## 
## Output summary
##                                                                                                               Func
## 1 { mod <- lm(dist ~ speed, data = cars[ki, ]) p <- predict(mod, cars[-ki, ]) sqrt(mean((p - cars$dist[-ki])^2)) }
##   Stat     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.       sd
## 1    1 8.685515 12.44126 15.03503 15.16872 16.92814 23.77031 3.815427

The following example simplifies this a little. Setting targetin tells it what the input to predict should be and setting target="dist" tells it that the target is the dist element from targetin. You cannot set target=cars$dit[-ki] since target cannot be evaluated as an expression.

mbc(lm(dist ~ speed, data=cars[ki,]),
    targetin=cars[-ki,], target="dist",
    kfold=c(nrow(cars), 5),
    times=30)

## Run times (sec)
##                              Function Min. 1st Qu. Median        Mean
## 1 lm(dist ~ speed, data = cars[ki, ])    0       0      0 0.002666667
##   3rd Qu. Max.          sd neval
## 1       0 0.02 0.006914918    30
## 
## Output summary
##                                  Func Stat     Min.  1st Qu.   Median
## 1 lm(dist ~ speed, data = cars[ki, ]) rmse 6.076282 12.12332 14.96145
##       Mean  3rd Qu.     Max.       sd
## 1 15.34269 17.76412 25.26908 4.848105

Metrics

In the previous example, the output shows that the “Stat” is “rmse”, meaning that it calculated the root-mean-square error from the. predictions and target values. The metric, or statistic, calculated can be changed using the metric argument, which defaults to rmse. Three of the other options for metric are t, mis90, and sr27. These three all compare target values ($y$) to predicted values ($\hat{y}$) and predicted errors ($s$). These only work for models that give predicted errors, such as Gaussian process models.

`metric=t`

Using these the target value, predicted value, and predicted error, we can calculate a t-score.

\[ t = \frac{\hat{y} – y}{s} \]

The output then shows the distribution of these t-scores by showing the six number summary.

mbc(lm(dist ~ speed, data=cars[ki,]),
    targetin=cars[-ki,], target="dist",
    kfold=c(nrow(cars), 5),
    times=30,
    metric='t')

## Run times (sec)
##                              Function Min. 1st Qu. Median        Mean
## 1 lm(dist ~ speed, data = cars[ki, ])    0       0      0 0.001666667
##   3rd Qu. Max.          sd neval
## 1       0 0.02 0.005306686    30
## 
## Output summary
##                                  Func      Stat       Min.     1st Qu.
## 1 lm(dist ~ speed, data = cars[ki, ])    Min. t -20.470581 -11.4997530
## 2 lm(dist ~ speed, data = cars[ki, ]) 1st Qu. t  -9.568407  -3.0442251
## 3 lm(dist ~ speed, data = cars[ki, ])  Median t  -2.907563  -0.3339233
## 4 lm(dist ~ speed, data = cars[ki, ])    Mean t  -4.247618  -1.0958828
## 5 lm(dist ~ speed, data = cars[ki, ]) 3rd Qu. t  -1.074035   2.0308401
## 6 lm(dist ~ speed, data = cars[ki, ])    Max. t   3.122954   5.1532598
##       Median       Mean    3rd Qu.      Max.       sd
## 1 -9.4512544 -9.7160202 -4.9702389 -2.335433 5.301566
## 2 -1.0500575 -1.9403262 -0.5463969  1.446661 2.387126
## 3  0.8057603  0.7977990  2.3060365  4.065288 1.847873
## 4  0.4079040  0.1078723  1.3042845  4.196559 1.950348
## 5  2.6977232  2.9430711  4.0937671  8.182076 1.966320
## 6  6.4729256  6.8863984  8.8773532 10.434104 2.208616

`metric=mis90`

The t-score metric is not very informative because you can get the same t-scores by having a large error and large predicted error as having a small error and small predicted error. mis90 is the mean interval score for 90% coverage intervals as described by Gneiting and Raftery (2007 Equation 43).

\[ 3.28s + 20 \left( \hat{y} – y – 1.64s \right)^+ + 20 \left( y – \hat{y} – 1.64s \right)^+ \] where $()+$ denotes the positive part of what is in the parantheses. Smaller values are better. This metric penalizes having large predicted errors and having actual errors different from the predicted errors, so it is very good for judging the accuracy of a prediction interval.

`metric=sr27`

The scoring rule in Equation 27 Gneiting and Raftery (2007) is another proper scoring rule.

\[ -\left( \frac{\hat{y} – y}{s} \right)^2 – \log s^2 \] For this metric, larger values are better. A problem with this metric is that if $s=0$, which can happen from numerical issues, then it will go to infinity, which does not happen with the mean interval score.

References

Gneiting, Tilmann, and Adrian E. Raftery. 2007. “Strictly Proper Scoring Rules, Prediction, and Estimation.” Journal of the American Statistical Association. Taylor & Francis, 359–78.

comparer