Running an experiment using CAISEr
With the definitions above it is possible now to run an experiment using the iterative sample size determination implemented in CAISEr. For that, all we have to do is define the desired experimental parameters and use run_experiment():
library(CAISEr)
my.results <- run_experiment(Instance.list = Instance.list,
Algorithm.list = Algorithm.list,
power = 0.8, # Desired power: 80%
d = 0.5, # to detect differences greater
# than 0.5 standard deviations
sig.level = 0.05, # at a 95% confidence level.
se.max = 0.05, # Measurement error: 5%
dif = "perc", # on the paired percent
# differences of means,
method = "boot", # calculated using bootstrap.
nstart = 15, # Start with 20 runs/algo/inst
nmax = 200, # and do no more than 200 runs/inst
seed = 1234) # PRNG seed (for reproducibility)After that we can interrogate the results and perform inference, if we are so inclined. For instance, we can check if our sample of paired differences in performance is (at least approximately) Normal, so that we can assume a Normal sampling distribution of the means and use a t test with a clean conscience:
# Take a look at the data summary:
# summary(my.results)
# print(my.results)
suppressPackageStartupMessages(library(car))
car::qqPlot(my.results$data.summary$phi.j,
pch = 20, las = 1,
ylab = "observed results", xlab = "theoretical quantiles")
## [1] 3 18The normal QQ plot indicates that no expressive deviations of normality are present, which gives us confidence in using the t test as our inferential procedure of choice (as the sampling distribution of the means will be even more “well-behaved” than the data distribution).
It is also interesting to observe a few things from the summary table. First, we observed negative values of phi.j in the majority of instances tested, which suggests an advantage of the MOEA/D-DE over the original MOEA/D (remember, smaller = better for the quality indicator used). The MOEA/D-DE seems to require smaller sample sizes in most instances, suggesting a smaller variance of performance, which is also a desirable feature. Also, in three instances (UF_5_28, UF_3_29 and UF__7) the maximum number of runs/instance (nmax = 200 in the run_experiment() call) was not enough to reduce the standard error (our “measurement error” on the values of phi.j) below the predefined threshold of \(0.05\). There is no reason to worry in this particular case, however, since the resulting standard errors were not particularly high, and therefore their effect on the test power (resulting from the increased uncertainty in the estimation of these particular phi.j values) will be insignificant.
Since our observations phi.j already express paired differences per instance, we can compare the two algorithms using a simple, one-sample t.test:
t.test(my.results$data.summary$phi.j)##
## One Sample t-test
##
## data: my.results$data.summary$phi.j
## t = -5.627, df = 33, p-value = 2.897e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.5165943 -0.2422327
## sample estimates:
## mean of x
## -0.3794135which indicates a statistically significant advantage of the MOEA/D-DE over the original MOEA/D (\(p = 2.90\times10^{-6}, df=33\)), with estimated mean IGD gains of \(37.94\%\) (\(CI_{0.95} = [24.22\%,51.66\%]\)) over the original MOEA/D for instances belonging to the problem class of interest.
These results could also be used to motivate further analyses. For instance, we can observe in the summary table that the only two cases for which the MOEA/D was substantially better than the MOEA/D-DE were for different dimensions of problem UF_5, which could suggest that some specific feature of this problem jeopardizes the latter algorithm’s search mechanism. This could motivate research on which particular aspect of this problem results in this loss of performance, and on how to improve the MOEA/D-DE.
Finally, the full data of the experiment is contained in other fields of the output list my.results, and the user is encouraged to explore these. As an example, we can generate box plots and confidence intervals on the mean performance of each algorithm on each sampled instance, which could inspire new questions for the researcher.
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(ggridges))
# Adjust instance names for plotting
mydata <- my.results$data.raw
mydata$Instance <- gsub(pattern = "UF\\_", replacement = "UF", mydata$Instance)
mydata$Instance <- gsub(pattern = "\\_", replacement = " (", mydata$Instance)
mydata$Instance <- sapply(mydata$Instance, FUN = function(x){paste0(x, ")")})
ggplot2::ggplot(mydata,
aes(x = Observation, y = Instance, fill = Algorithm)) +
ggridges::geom_density_ridges(alpha = 0.7) +
ggplot2::ggtitle("Estimated IGD distribution",
subtitle = "for each algorithm on each instance") +
ggplot2::theme(legend.position = "bottom")## Picking joint bandwidth of 0.0132
# Calculate confidence intervals for each instance
algos <- unique(mydata$Algorithm)
ninstances <- length(my.results$instances.sampled)
CIs <- data.frame(instance = rep(unique(mydata$Instance, times = 2)),
algorithm = rep(algos, each = ninstances),
x.est = 0, CI.l = 0, CI.u = 0)
for (i in 1:ninstances){
tmpdata <- mydata %>%
filter(Instance == unique(mydata$Instance)[i])
myt1 <- t.test(tmpdata$Observation[tmpdata$Algorithm == algos[1]])
myt2 <- t.test(tmpdata$Observation[tmpdata$Algorithm == algos[2]])
CIs[i,3:5] <- c(myt1$estimate, as.numeric(myt1$conf.int))
CIs[i + ninstances,3:5] <- c(myt2$estimate, as.numeric(myt2$conf.int))
}
# Plot individual confidence intervals for each instance
myplot <- ggplot2::ggplot(CIs, aes(x = instance,
y = x.est, ymin = CI.l, ymax = CI.u,
group = algorithm, colour = algorithm,
fill = algorithm))
myplot +
ggplot2::geom_pointrange(position = position_dodge(width = 0.5), alpha = 0.7) +
ggplot2::xlab("Instance") + ggplot2::ylab("IGD") +
ggplot2::ggtitle("Estimated mean IGD",
subtitle = "for each algorithm on each instance") +
ggplot2::theme(legend.position = "bottom",
axis.text.x = element_text(angle = 55, hjust = 1, size = 6))