# Comparing two algorithms. Paired t-test, independent test or mixed effect linear regression

I have two algorithms running on 100 benchmark instances.
I want to test if there are significant differences in their mean average error.

A significant part of those two algorithms is not deterministic.

My data with only one run looks like this:

+--------------------+-------+-----------+
| benchmark instance | error | algorithm |
+--------------------+-------+-----------+
| BENCH1             | 1.1   | ALGO1     |
+--------------------+-------+-----------+
| BENCH2             | 5.4   | ALGO1     |
+--------------------+-------+-----------+
| BENCH3             | 7.8   | ALGO1     |
+--------------------+-------+-----------+
| ...                | ...   | ...       |
+--------------------+-------+-----------+
| BENCH1             | 2.5   | ALGO2     |
+--------------------+-------+-----------+
| BENCH2             | 4.4   | ALGO2     |
+--------------------+-------+-----------+
| BENCH3             | 7.7   | ALGO2     |
+--------------------+-------+-----------+
| ...                | ...   | ...       |
+--------------------+-------+-----------+
| BENCH1             | 1.9   | ALGO1     |
+--------------------+-------+-----------+
| BENCH2             | 5.2   | ALGO1     |
+--------------------+-------+-----------+
| BENCH3             | 6.1   | ALGO1     |
+--------------------+-------+-----------+


In this case I use 2-sample t-test.

t.test(error ~ algorithm, data = combined.data)


If I run the 100 tests for a second or a third time for each algorithm, do I need to use paired t-test or not?

t.test(error ~ algorithm, data = combined.data, paired = TRUE)


I use as a rule of thumb that when measurements on same group of objects are made twice, then in order to compare the difference paired T-test is used. Is that correct in that case? Should a repeated test on the same instance with the same algorithm considered as the same object?

The fact is that I need to run each algorithm many times on the same instances in order to be more precise. Some facts like CPU lagging or another app use some portion of the CPU or the randomized nature of some commands inside the algorithm may give an error of lets say 11% and the second time the same algorithm on the same instance will give 9%. If I run it 50 times I can say for sure that the quality of this algorithm ranges from 9-11%.

The final goal is to find out which algorithm performed better by comparing the error column.

Thanks

## --------------------------------------------------------

Data with 5 repetitions * 7 instances * 2 Algorithms.

+----------+-------+---------+
| bench_id | Algo  |  Score  |
+----------+-------+---------+
| BENCH25  | ALGO1 |  4.7601 |
| BENCH14  | ALGO1 |  2.9453 |
| BENCH13  | ALGO1 |  3.0815 |
| BENCH10  | ALGO1 |  6.3272 |
| BENCH6   | ALGO1 |  8.7662 |
| BENCH101 | ALGO1 |  6.7073 |
| BENCH200 | ALGO1 |  7.7389 |
| BENCH25  | ALGO1 |  4.7404 |
| BENCH14  | ALGO1 |  2.8845 |
| BENCH13  | ALGO1 |  4.0013 |
| BENCH10  | ALGO1 |  6.3065 |
| BENCH6   | ALGO1 |  9.6564 |
| BENCH101 | ALGO1 |  7.1119 |
| BENCH200 | ALGO1 |   8.539 |
| BENCH25  | ALGO1 |  4.4898 |
| BENCH14  | ALGO1 |  2.8129 |
| BENCH13  | ALGO1 |  3.7915 |
| BENCH10  | ALGO1 |  7.1798 |
| BENCH6   | ALGO1 |  8.1818 |
| BENCH101 | ALGO1 |  5.8875 |
| BENCH200 | ALGO1 |     8.6 |
| BENCH25  | ALGO1 |  3.0398 |
| BENCH14  | ALGO1 |  2.8487 |
| BENCH13  | ALGO1 |  4.4973 |
| BENCH10  | ALGO1 |  7.6505 |
| BENCH6   | ALGO1 |  7.9472 |
| BENCH101 | ALGO1 |  6.4629 |
| BENCH200 | ALGO1 |  9.0849 |
| BENCH25  | ALGO1 |  4.3706 |
| BENCH14  | ALGO1 |  2.7197 |
| BENCH13  | ALGO1 |  4.1304 |
| BENCH10  | ALGO1 |  6.4852 |
| BENCH6   | ALGO1 |  8.6536 |
| BENCH101 | ALGO1 |  6.7073 |
| BENCH200 | ALGO1 |  9.2053 |
| BENCH25  | ALGO2 |  5.1921 |
| BENCH14  | ALGO2 |  4.4924 |
| BENCH13  | ALGO2 |  6.1262 |
| BENCH10  | ALGO2 |  5.9123 |
| BENCH6   | ALGO2 |  9.0649 |
| BENCH101 | ALGO2 |  6.9505 |
| BENCH200 | ALGO2 | 10.0394 |
| BENCH25  | ALGO2 |   6.054 |
| BENCH14  | ALGO2 |  4.4786 |
| BENCH13  | ALGO2 |   5.641 |
| BENCH10  | ALGO2 |  9.8324 |
| BENCH6   | ALGO2 |  9.9433 |
| BENCH101 | ALGO2 |  7.6724 |
| BENCH200 | ALGO2 |  8.8431 |
| BENCH25  | ALGO2 |  6.3283 |
| BENCH14  | ALGO2 |  3.9496 |
| BENCH13  | ALGO2 |  6.8156 |
| BENCH10  | ALGO2 |  6.8474 |
| BENCH6   | ALGO2 |  7.9981 |
| BENCH101 | ALGO2 |  8.3048 |
| BENCH200 | ALGO2 | 10.5675 |
| BENCH25  | ALGO2 |  5.2637 |
| BENCH14  | ALGO2 |  4.4197 |
| BENCH13  | ALGO2 |  5.5159 |
| BENCH10  | ALGO2 |   8.963 |
| BENCH6   | ALGO2 |    8.39 |
| BENCH101 | ALGO2 |  7.5928 |
| BENCH200 | ALGO2 |  9.9212 |
| BENCH25  | ALGO2 |  5.4682 |
| BENCH14  | ALGO2 |  5.0907 |
| BENCH13  | ALGO2 |  5.5636 |
| BENCH10  | ALGO2 |  8.8328 |
| BENCH6   | ALGO2 |  8.2008 |
| BENCH101 | ALGO2 |  7.2727 |
| BENCH200 | ALGO2 |  9.5646 |
+----------+-------+---------+


Using 1st approach like this:

b_var <- 1/var(combined.data$score) res <- lmer( score ~ algo + (1 | bench_id), combined.data, weights = 1 / b_var) summary(res)  I get Error variable lengths differ (found for '(weights)') Using 2nd approach I get: res1 <- lmer( score ~ algo + (1 | bench_id), combined.data) summary(res1)  Linear mixed model fit by REML ['lmerMod'] Formula: score ~ algo + (1 | bench_id) Data: combined.data REML criterion at convergence: 183.2 Scaled residuals: Min 1Q Median 3Q Max -2.91976 -0.48884 0.00928 0.59819 2.47335 Random effects: Groups Name Variance Std.Dev. bench_id (Intercept) 4.3060 2.0751 Residual 0.5283 0.7269 Number of obs: 70, groups: bench_id, 7 Fixed effects: Estimate Std. Error t value (Intercept) 5.9518 0.7939 7.497 algoALGO2 1.2229 0.1738 7.038 Correlation of Fixed Effects: (Intr) algoALGO2 -0.109  Using 3rd approach I get: res2 <- lmer( score ~ algo + (0 + factor(algo) | bench_id), combined.data) summary(res2)  Linear mixed model fit by REML ['lmerMod'] Formula: score ~ algo + (0 + factor(algo) | bench_id) Data: combined.data REML criterion at convergence: 177.5 Scaled residuals: Min 1Q Median 3Q Max -2.91175 -0.51624 -0.00599 0.60168 2.73685 Random effects: Groups Name Variance Std.Dev. Corr bench_id factor(algo)ALGO1 5.2709 2.296 factor(algo)ALGO2 3.4470 1.857 1.00 Residual 0.4816 0.694 Number of obs: 70, groups: bench_id, 7 Fixed effects: Estimate Std. Error t value (Intercept) 5.9518 0.8756 6.797 algoALGO2 1.2229 0.2347 5.210 Correlation of Fixed Effects: (Intr) algoALGO2 -0.768 convergence code: 0 singular fit  Using t-test (paired): data: score by algo t = -7.0242, df = 34, p-value = 4.165e-08 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -1.5766716 -0.8690655 sample estimates: mean of the differences -1.222869  Using Wilcoxon test: wilcox.test(score ~ algo, data = combined.data)  Wilcoxon rank sum test with continuity correction data: score by algo W = 425, p-value = 0.02805 alternative hypothesis: true location shift is not equal to 0  • You have two samples, algo1 and algo2. But they are linked by the benches assuming bench1 is the same task for both algorithms. So you should use a paired t-test. If the benches are arbitrary, then independent samples. – Heteroskedastic Jim Nov 15 '18 at 12:56 • @HeteroskedasticJim would it have any negative impact if I run the test 2-3 times and add those rows? I mean would it contribute something or make it worse? – Antonis Nov 15 '18 at 12:58 • Why would you do this? Help explain. – Heteroskedastic Jim Nov 15 '18 at 12:59 • Because when I run only one test on one instance with Algorithm A for example. some facts like CPU lagging or another app use some portion of the CPU or the randomized nature of some commands inside the algorithm may give an error of lets say 11% and the second time the same algorithm on the same instance will give 9%. If I run it 50 times I can say for sure that the quality of this algorithm ranges from 9-11% – Antonis Nov 15 '18 at 13:04 • You should modify your question to make this part clear. The most natural approach would be a mixed effects regression. I'll write an answer. – Heteroskedastic Jim Nov 15 '18 at 13:07 ## 1 Answer The most natural approach to handle your situation would be a mixed effects or multilevel regression. You can also perform an inverse variance weighted regression. You have 100 benchmarks for both algorithms. Additionally, you plan to attempt each benchmark about 50 times to account for randomness in the benchmark results. ### Inverse variance weighting I'll begin with the inverse variance approach since it is easier. Assuming that the data for the 50 instances of a benchmark for an algorithm are reasonably normal, then the data can be summarized by its mean and variance. So for each benchmark and algorithm, calculate the mean and the variance of the 50 instances. And you will arrive at 100 average benchmarks for each algorithm with 100 variances. The averages become the outcome for a regression model where you correlate the scores by benchmark to account for the fact that each benchmark produced two averages. You also need to supply the variances as weights to the regression. Assume the benchmark averages are score, benchmark unique identifier is bench_id, algorithm indicator or dummy variable is algo, and benchmark variances are b_var, all contained in the dataset, dat. dat should be 200 rows long, then: library(lme4) res <- lmer( score ~ algo + (1 | bench_id), dat, weights = 1 / b_var) summary(res)  Should work. You'll get the difference between algorithms from its coefficients, and you'll get a bench random intercept standard deviation which I'll call $$v$$, and a residual standard deviation which I'll call $$s$$. If you perform: $$\frac{v^2}{v^2+s^2}$$ You'll get the average correlation between measurements from the same benchmark, this is known as the (residual) ICC. The syntax should be: as.numeric(VarCorr(res)) / (as.numeric(VarCorr(res)) + sigma(res) ^ 2)  ### Complete multilevel setup This time, you have all the data in long form. Each row pertains to one of 50 measurements for an algorithm for one of 100 benchmarks. So your number of rows should be $$10000\ (50\times 100\times 2)$$. The model you'll run still has a similar syntax to the 10000 long dataset: res1 <- lmer( score ~ algo + (1 | bench_id), dat.longer) summary(res1)  I would prefer to run: res2 <- lmer( score ~ algo + (0 + factor(algo) | bench_id), dat.longer) summary(res2)  This would allow you to observe the correlation between the scores from both algorithms and test whether the algorithms vary differently. Of the three models, I prefer the first and the last. The first is a simple approach and if the normality assumption for the 50 instances per benchmark per algorithm is satisfied, then it is reasonable and defensible. Alternatively, the final model uses all the data in the most natural way. However, it assumes homogeneity of variance which the first approach does not. Final model assumes the variance of the 5000 measurements within each algorithm is the same across the 100 benchmarks. Hope this helps. It might be useful to pick up a basic guide to multilevel modeling using lme4 if you are not familiar with such models. For your problem, especially if you stick with the first approach, an introductory guide is sufficient. EDIT Assuming the dataset is called data.long, then: # to create dataset of mean scores dat <- aggregate(score ~ algo + bench_id, data.long, mean) dat$$b_var <- aggregate(score ~ algo + bench_id, data.long, var)$$score  Then use dat for model 1 syntax above. This dataset will have only 14 rows since you currently have only 7 bench_ids. Try model 3 when you get more data with your combined.data. • Thank you so much for your thorough answer. I need to wait to accept it because I'm a new user. – Antonis Nov 15 '18 at 14:40 • @Antonis glad it helps. Welcome to the site. I'm used to new users not getting the formatting correct so I assumed you'd been here longer. – Heteroskedastic Jim Nov 15 '18 at 14:42 • I just wanted to give some feedback and ask some guidance. I hope I'm not bothering you. I run 7 instances for five times for both algos and I got back 70 lines. I tried all the models. Sadly only the second (the one that you do not recommend worked). I must do something wrong. At the first model I calculate weights like this: 1 / var(mydata$scores) and when I run the code it returns variable lengths differ (found for '(weights)'). When I run the 3rd model it returns singular fit. Does this imply that the algorithms do not differ? – Antonis Nov 15 '18 at 20:55
• @Antonis see the edit at the bottom of the question. – Heteroskedastic Jim Nov 16 '18 at 3:21