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

Question

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

--------------------------------------------------------

Addition,

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

Answer 1

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.