Chapter 4 Genotypic fitness sharing

Results for the genotypic fitness sharing parameter sweep on the diagnostics with valleys.

4.1 Data setup

over_time_df <- read.csv(paste(DATA_DIR,'OVER-TIME-MVC/gfs.csv', sep = "", collapse = NULL), header = TRUE, stringsAsFactors = FALSE)
over_time_df$Sigma <- factor(over_time_df$Sigma, levels = FS_LIST)

best_df <- read.csv(paste(DATA_DIR,'BEST-MVC/gfs.csv', sep = "", collapse = NULL), header = TRUE, stringsAsFactors = FALSE)
best_df$Sigma <- factor(best_df$Sigma, levels = FS_LIST)

4.2 Exploitation rate results

Here we present the results for best performances found by each selection scheme parameter on the exploitation rate diagnostic with valleys. 50 replicates are conducted for each scheme explored.

4.2.1 Performance over time

Best performance in a population over time. Data points on the graph is the average performance across 50 replicates every 2000 generations. Shading comes from the best and worse performance across 50 replicates.

lines = filter(over_time_df, acro == 'exp') %>%
  group_by(Sigma, gen) %>%
  dplyr::summarise(
    min = min(pop_fit_max) / DIMENSIONALITY,
    mean = mean(pop_fit_max) / DIMENSIONALITY,
    max = max(pop_fit_max) / DIMENSIONALITY
  )

## `summarise()` has grouped output by 'Sigma'. You can override using the
## `.groups` argument.

over_time_plot = ggplot(lines, aes(x=gen, y=mean, group = Sigma, fill = Sigma, color = Sigma, shape = Sigma)) +
  geom_ribbon(aes(ymin = min, ymax = max), alpha = 0.1) +
  geom_line(size = 0.5) +
  geom_point(data = filter(lines, gen %% 2000 == 0 & gen != 0), size = 1.5, stroke = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Average trait score",
    limits = c(0,50)
  ) +
  scale_x_continuous(
    name="Generations",
    limits=c(0, 50000),
    breaks=c(0, 10000, 20000, 30000, 40000, 50000),
    labels=c("0e+4", "1e+4", "2e+4", "3e+4", "4e+4", "5e+4")

  ) +
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Performance over time')+
  p_theme +
  guides(
    shape=guide_legend(nrow=1, title.position = "left", title = 'Sigma'),
    color=guide_legend(nrow=1, title.position = "left", title = 'Sigma'),
    fill=guide_legend(nrow=1, title.position = "left", title = 'Sigma')
  )

over_time_plot

4.2.2 Best performance throughout

Best performance reached throughout 50,000 generations in a population.

plot = filter(best_df, acro == 'exp' &  var == 'pop_fit_max') %>%
  ggplot(., aes(x = Sigma, y = val / DIMENSIONALITY, color = Sigma, fill = Sigma, shape = Sigma)) +
  geom_flat_violin(position = position_nudge(x = .1, y = 0), scale = 'width', alpha = 0.2, width = 1.5) +
  geom_boxplot(color = 'black', width = .07, outlier.shape = NA, alpha = 0.0, size = 1.0, position = position_nudge(x = .16, y = 0)) +
  geom_point(position = position_jitter(width = 0.03, height = 0.02), size = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Average trait score",
    limits=c(30,50)
  ) +
  scale_x_discrete(
    name="Sigma"
  )+
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette, ) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Best performance throughout')+
  p_theme

plot_grid(
  plot +
    theme(legend.position="none"),
  legend,
  nrow=2,
  rel_heights = c(3,1)
)

4.2.2.1 Stats

Summary statistics for the best performance.

performance = filter(best_df, acro == 'exp' &  var == 'pop_fit_max')
performance$Sigma = factor(performance$Sigma, levels = FS_LIST)
performance %>%
  group_by(Sigma) %>%
  dplyr::summarise(
    count = n(),
    na_cnt = sum(is.na(val)),
    min = min(val / DIMENSIONALITY, na.rm = TRUE),
    median = median(val / DIMENSIONALITY, na.rm = TRUE),
    mean = mean(val / DIMENSIONALITY, na.rm = TRUE),
    max = max(val / DIMENSIONALITY, na.rm = TRUE),
    IQR = IQR(val / DIMENSIONALITY, na.rm = TRUE)
  )

## # A tibble: 7 x 8
##   Sigma count na_cnt   min median  mean   max   IQR
##   <fct> <int>  <int> <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 0        50      0  42.0   43.9  43.9  46.7 1.58 
## 2 0.1      50      0  41.8   43.3  43.5  45.7 1.35 
## 3 0.3      50      0  40.8   43.1  43.2  45.5 1.11 
## 4 0.6      50      0  41.0   42.4  42.4  44.4 0.963
## 5 1.2      50      0  39.0   40.6  40.5  41.9 1.01 
## 6 2.5      50      0  31.3   33.3  33.3  35.4 1.07 
## 7 5        50      0  33.4   35.4  35.3  37.1 1.25

Kruskal–Wallis test illustrates evidence of statistical differences.

kruskal.test(val ~ Sigma, data = performance)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  val by Sigma
## Kruskal-Wallis chi-squared = 294.25, df = 6, p-value < 2.2e-16

Results for post-hoc Wilcoxon rank-sum test with a Bonferroni correction.

pairwise.wilcox.test(x = performance$val, g = performance$Sigma, p.adjust.method = "bonferroni",
                     paired = FALSE, conf.int = FALSE, alternative = 't')

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  performance$val and performance$Sigma 
## 
##     0       0.1     0.3     0.6     1.2     2.5    
## 0.1 0.6669  -       -       -       -       -      
## 0.3 0.0125  1.0000  -       -       -       -      
## 0.6 1.5e-09 1.5e-05 0.0024  -       -       -      
## 1.2 < 2e-16 < 2e-16 1.1e-15 6.0e-15 -       -      
## 2.5 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 -      
## 5   < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 1.7e-13
## 
## P value adjustment method: bonferroni

4.2.3 Largest valley reached throughout

The largest valley reached in a single trait throughout an entire evolutionary run. To collect this value, we look through all the best-performing solutions each generation and find the largest valley reached.

plot = filter(best_df, acro == 'exp' & var == 'ele_big_peak') %>%
  ggplot(., aes(x = Sigma, y = val, color = Sigma, fill = Sigma, shape = Sigma)) +
  geom_flat_violin(position = position_nudge(x = .1, y = 0), scale = 'width', alpha = 0.2, width = 1.5) +
  geom_boxplot(color = 'black', width = .07, outlier.shape = NA, alpha = 0.0, size = 1.0, position = position_nudge(x = .16, y = 0)) +
  geom_point(position = position_jitter(width = 0.03, height = 0.02), size = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Valley reached",
    limits=c(11.9,14.1),
    breaks = c(11,12,13,14)
  ) +
  scale_x_discrete(
    name="Sigma"
  )+
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette, ) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Largest valley reached')+
  p_theme

plot_grid(
  plot +
    theme(legend.position="none"),
  legend,
  nrow=2,
  rel_heights = c(3,1)
)

4.2.3.1 Stats

Summary statistics for the largest valley crossed.

valleys = filter(best_df, acro == 'exp' & var == 'ele_big_peak')
valleys$Sigma = factor(valleys$Sigma, levels = FS_LIST)
valleys %>%
  group_by(Sigma) %>%
  dplyr::summarise(
    count = n(),
    na_cnt = sum(is.na(val)),
    min = min(val, na.rm = TRUE),
    median = median(val, na.rm = TRUE),
    mean = mean(val, na.rm = TRUE),
    max = max(val, na.rm = TRUE),
    IQR = IQR(val, na.rm = TRUE)
  )

## # A tibble: 7 x 8
##   Sigma count na_cnt   min median  mean   max   IQR
##   <fct> <int>  <int> <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 0        50      0    12     13  13.1    14  0   
## 2 0.1      50      0    12     13  13.2    14  1   
## 3 0.3      50      0    12     13  13.3    14  1   
## 4 0.6      50      0    12     13  13.3    14  1   
## 5 1.2      50      0    13     13  13.5    14  1   
## 6 2.5      50      0    12     13  13.0    14  0.75
## 7 5        50      0    12     13  13.3    14  1

Kruskal–Wallis test illustrates evidence of statistical differences.

kruskal.test(val ~ Sigma, data = valleys)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  val by Sigma
## Kruskal-Wallis chi-squared = 16.366, df = 6, p-value = 0.01192

Results for post-hoc Wilcoxon rank-sum test with a Bonferroni correction.

pairwise.wilcox.test(x = valleys$val, g = valleys$Sigma, p.adjust.method = "bonferroni",
                     paired = FALSE, conf.int = FALSE, alternative = 't')

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  valleys$val and valleys$Sigma 
## 
##     0     0.1   0.3   0.6   1.2   2.5  
## 0.1 1.000 -     -     -     -     -    
## 0.3 0.869 1.000 -     -     -     -    
## 0.6 1.000 1.000 1.000 -     -     -    
## 1.2 0.042 1.000 1.000 1.000 -     -    
## 2.5 1.000 1.000 0.589 1.000 0.046 -    
## 5   0.869 1.000 1.000 1.000 1.000 0.589
## 
## P value adjustment method: bonferroni

4.3 Ordered exploitation results

Here we present the results for best performances found by each selection scheme parameter on the ordered exploitation diagnostic with valleys. 50 replicates are conducted for each scheme explored.

4.3.1 Performance over time

Best performance in a population over time. Data points on the graph is the average performance across 50 replicates every 2000 generations. Shading comes from the best and worse performance across 50 replicates.

lines = filter(over_time_df, acro == 'ord') %>%
  group_by(Sigma, gen) %>%
  dplyr::summarise(
    min = min(pop_fit_max) / DIMENSIONALITY,
    mean = mean(pop_fit_max) / DIMENSIONALITY,
    max = max(pop_fit_max) / DIMENSIONALITY
  )

## `summarise()` has grouped output by 'Sigma'. You can override using the
## `.groups` argument.

over_time_plot = ggplot(lines, aes(x=gen, y=mean, group = Sigma, fill = Sigma, color = Sigma, shape = Sigma)) +
  geom_ribbon(aes(ymin = min, ymax = max), alpha = 0.1) +
  geom_line(size = 0.5) +
  geom_point(data = filter(lines, gen %% 2000 == 0 & gen != 0), size = 1.5, stroke = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Average trait score",
    limits = c(0,15)
  ) +
  scale_x_continuous(
    name="Generations",
    limits=c(0, 50000),
    breaks=c(0, 10000, 20000, 30000, 40000, 50000),
    labels=c("0e+4", "1e+4", "2e+4", "3e+4", "4e+4", "5e+4")

  ) +
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Performance over time')+
  p_theme +
  guides(
    shape=guide_legend(nrow=1, title.position = "left", title = 'Sigma'),
    color=guide_legend(nrow=1, title.position = "left", title = 'Sigma'),
    fill=guide_legend(nrow=1, title.position = "left", title = 'Sigma')
  )

over_time_plot

4.3.2 Best performance throughout

Best performance reached throughout 50,000 generations in a population.

plot = filter(best_df, acro == 'ord' &  var == 'pop_fit_max') %>%
  ggplot(., aes(x = Sigma, y = val / DIMENSIONALITY, color = Sigma, fill = Sigma, shape = Sigma)) +
  geom_flat_violin(position = position_nudge(x = .1, y = 0), scale = 'width', alpha = 0.2, width = 1.5) +
  geom_boxplot(color = 'black', width = .07, outlier.shape = NA, alpha = 0.0, size = 1.0, position = position_nudge(x = .16, y = 0)) +
  geom_point(position = position_jitter(width = 0.03, height = 0.02), size = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Average trait score",
    limits=c(8,14)
  ) +
  scale_x_discrete(
    name="Sigma"
  )+
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette, ) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Best performance throughout')+
  p_theme

plot_grid(
  plot +
    theme(legend.position="none"),
  legend,
  nrow=2,
  rel_heights = c(3,1)
)

4.3.2.1 Stats

Summary statistics for the best performance.

performance = filter(best_df, acro == 'ord' &  var == 'pop_fit_max')
performance$Sigma = factor(performance$Sigma, levels = FS_LIST)
performance %>%
  group_by(Sigma) %>%
  dplyr::summarise(
    count = n(),
    na_cnt = sum(is.na(val)),
    min = min(val / DIMENSIONALITY, na.rm = TRUE),
    median = median(val / DIMENSIONALITY, na.rm = TRUE),
    mean = mean(val / DIMENSIONALITY, na.rm = TRUE),
    max = max(val / DIMENSIONALITY, na.rm = TRUE),
    IQR = IQR(val / DIMENSIONALITY, na.rm = TRUE)
  )

## # A tibble: 7 x 8
##   Sigma count na_cnt   min median  mean   max   IQR
##   <fct> <int>  <int> <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 0        50      0  8.62   10.4  10.5  12.5 0.761
## 2 0.1      50      0  9.96   11.6  11.5  12.8 0.898
## 3 0.3      50      0 10.1    11.6  11.5  13.3 0.840
## 4 0.6      50      0  9.93   12.0  12.0  13.9 0.908
## 5 1.2      50      0 11.2    12.2  12.2  13.9 0.821
## 6 2.5      50      0 10.1    11.5  11.4  12.4 0.579
## 7 5        50      0 10.2    11.4  11.4  12.7 0.629

Kruskal–Wallis test illustrates evidence of statistical differences.

kruskal.test(val ~ Sigma, data = performance)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  val by Sigma
## Kruskal-Wallis chi-squared = 128.96, df = 6, p-value < 2.2e-16

Results for post-hoc Wilcoxon rank-sum test with a Bonferroni correction.

pairwise.wilcox.test(x = performance$val, g = performance$Sigma, p.adjust.method = "bonferroni",
                     paired = FALSE, conf.int = FALSE, alternative = 't')

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  performance$val and performance$Sigma 
## 
##     0       0.1     0.3     0.6     1.2     2.5    
## 0.1 5.3e-09 -       -       -       -       -      
## 0.3 6.1e-09 1.00000 -       -       -       -      
## 0.6 4.5e-12 0.04909 0.06661 -       -       -      
## 1.2 8.0e-15 3.3e-05 0.00010 1.00000 -       -      
## 2.5 1.9e-09 1.00000 1.00000 0.00035 2.6e-08 -      
## 5   2.0e-09 1.00000 1.00000 0.00056 8.8e-08 1.00000
## 
## P value adjustment method: bonferroni

4.3.3 Largest valley reached throughout

The largest valley reached in a single trait throughout an entire evolutionary run. To collect this value, we look through all the best-performing solutions each generation and find the largest valley reached.

plot = filter(best_df, acro == 'ord' & var == 'ele_big_peak') %>%
  ggplot(., aes(x = Sigma, y = val, color = Sigma, fill = Sigma, shape = Sigma)) +
  geom_flat_violin(position = position_nudge(x = .1, y = 0), scale = 'width', alpha = 0.2, width = 1.5) +
  geom_boxplot(color = 'black', width = .07, outlier.shape = NA, alpha = 0.0, size = 1.0, position = position_nudge(x = .16, y = 0)) +
  geom_point(position = position_jitter(width = 0.03, height = 0.02), size = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Valley reached",
    limits=c(10.9,14.1),
    breaks = c(11,12,13,14)
  ) +
  scale_x_discrete(
    name="Sigma"
  )+
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette, ) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Largest valley reached')+
  p_theme

plot_grid(
  plot +
    theme(legend.position="none"),
  legend,
  nrow=2,
  rel_heights = c(3,1)
)

4.3.3.1 Stats

Summary statistics for the largest valley crossed.

valleys = filter(best_df, acro == 'ord' & var == 'ele_big_peak')
valleys$Sigma = factor(valleys$Sigma, levels = FS_LIST)
valleys %>%
  group_by(Sigma) %>%
  dplyr::summarise(
    count = n(),
    na_cnt = sum(is.na(val)),
    min = min(val, na.rm = TRUE),
    median = median(val, na.rm = TRUE),
    mean = mean(val, na.rm = TRUE),
    max = max(val, na.rm = TRUE),
    IQR = IQR(val, na.rm = TRUE)
  )

## # A tibble: 7 x 8
##   Sigma count na_cnt   min median  mean   max   IQR
##   <fct> <int>  <int> <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 0        50      0    11     12  11.9    13     0
## 2 0.1      50      0    12     13  12.7    14     1
## 3 0.3      50      0    12     13  12.8    14     1
## 4 0.6      50      0    11     13  13.2    14     1
## 5 1.2      50      0    13     14  13.6    14     1
## 6 2.5      50      0    13     14  13.9    14     0
## 7 5        50      0    14     14  14      14     0

Kruskal–Wallis test illustrates evidence of statistical differences.

kruskal.test(val ~ Sigma, data = valleys)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  val by Sigma
## Kruskal-Wallis chi-squared = 231.15, df = 6, p-value < 2.2e-16

Results for post-hoc Wilcoxon rank-sum test with a Bonferroni correction.

pairwise.wilcox.test(x = valleys$val, g = valleys$Sigma, p.adjust.method = "bonferroni",
                     paired = FALSE, conf.int = FALSE, alternative = 't')

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  valleys$val and valleys$Sigma 
## 
##     0       0.1     0.3     0.6     1.2     2.5   
## 0.1 7.9e-06 -       -       -       -       -     
## 0.3 1.1e-07 1.0000  -       -       -       -     
## 0.6 6.4e-11 0.0013  0.1424  -       -       -     
## 1.2 1.0e-15 1.2e-10 9.6e-08 0.0093  -       -     
## 2.5 < 2e-16 3.3e-16 9.9e-14 1.6e-08 0.0165  -     
## 5   < 2e-16 < 2e-16 3.9e-16 3.8e-11 6.8e-05 0.9092
## 
## P value adjustment method: bonferroni

4.4 Contradictory objectives results

Here we present the results for activation gene coverage and satisfactory trait coverage found by each selection scheme parameter on the contradictory objectives diagnostic with valleys. 50 replicates are conducted for each scheme parameters explored.

4.4.1 Activation gene coverage over time

Activation gene coverage in a population over time. Data points on the graph is the average activation gene coverage across 50 replicates every 2000 generations. Shading comes from the best and worse coverage across 50 replicates.

lines = filter(over_time_df, acro == 'con')  %>%
  group_by(Sigma, gen) %>%
  dplyr::summarise(
    min = min(uni_str_pos),
    mean = mean(uni_str_pos),
    max = max(uni_str_pos)
  )

## `summarise()` has grouped output by 'Sigma'. You can override using the
## `.groups` argument.

ggplot(lines, aes(x=gen, y=mean, group = Sigma, fill = Sigma, color = Sigma, shape = Sigma)) +
  geom_ribbon(aes(ymin = min, ymax = max), alpha = 0.1) +
  geom_line(size = 0.5) +
  geom_point(data = filter(lines, gen %% 2000 == 0 & gen != 0), size = 1.5, stroke = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Coverage",
    limits=c(0, 100),
    breaks=seq(0,100, 20),
    labels=c("0", "20", "40", "60", "80", "100")
  ) +
  scale_x_continuous(
    name="Generations",
    limits=c(0, 50000),
    breaks=c(0, 10000, 20000, 30000, 40000, 50000),
    labels=c("0e+4", "1e+4", "2e+4", "3e+4", "4e+4", "5e+4")

  ) +
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Activation gene coverage over time')+
  p_theme +
  guides(
    shape=guide_legend(nrow=1, title.position = "left", title = 'Sigma'),
    color=guide_legend(nrow=1, title.position = "left", title = 'Sigma'),
    fill=guide_legend(nrow=1, title.position = "left", title = 'Sigma')
  )

4.4.2 Final activation gene coverage

Activation gene coverage found in the final population at 50,000 generations.

plot = filter(over_time_df, gen == 50000 & acro == 'con') %>%
  ggplot(., aes(x = Sigma, y = uni_str_pos, color = Sigma, fill = Sigma, shape = Sigma)) +
  geom_flat_violin(position = position_nudge(x = .1, y = 0), scale = 'width', alpha = 0.2, width = 1.5) +
  geom_boxplot(color = 'black', width = .07, outlier.shape = NA, alpha = 0.0, size = 1.0, position = position_nudge(x = .16, y = 0)) +
  geom_point(position = position_jitter(width = 0.03, height = 0.02), size = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Coverage",
    limits=c(0, 50.1)
  ) +
  scale_x_discrete(
    name="Sigma"
  )+
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette, ) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Final activation gene coverage')+
  p_theme

plot_grid(
  plot +
    theme(legend.position="none"),
  legend,
  nrow=2,
  rel_heights = c(3,1)
)

4.4.2.1 Stats

Summary statistics for the generation a satisfactory solution is found.

act_coverage = filter(over_time_df, gen == 50000 & acro == 'con')
act_coverage$Sigma = factor(act_coverage$Sigma, levels = FS_LIST)
act_coverage %>%
  group_by(Sigma) %>%
  dplyr::summarise(
    count = n(),
    na_cnt = sum(is.na(uni_str_pos)),
    min = min(uni_str_pos, na.rm = TRUE),
    median = median(uni_str_pos, na.rm = TRUE),
    mean = mean(uni_str_pos, na.rm = TRUE),
    max = max(uni_str_pos, na.rm = TRUE),
    IQR = IQR(uni_str_pos, na.rm = TRUE)
  )

## # A tibble: 7 x 8
##   Sigma count na_cnt   min median  mean   max   IQR
##   <fct> <int>  <int> <int>  <dbl> <dbl> <int> <dbl>
## 1 0        50      0     5   30    30.2    45 10.5 
## 2 0.1      50      0    10   31    30.3    43  7.5 
## 3 0.3      50      0     1   32.5  30      40 10.8 
## 4 0.6      50      0     1   31    27.9    46 15.5 
## 5 1.2      50      0     2   32    30.4    50 16   
## 6 2.5      50      0     1   14    13.6    34  9.75
## 7 5        50      0     2   13    12.1    27  9.5

Kruskal–Wallis test illustrates evidence of statistical differences.

kruskal.test(uni_str_pos ~ Sigma, data = act_coverage)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  uni_str_pos by Sigma
## Kruskal-Wallis chi-squared = 157.28, df = 6, p-value < 2.2e-16

Results for post-hoc Wilcoxon rank-sum test with a Bonferroni correction.

pairwise.wilcox.test(x = act_coverage$uni_str_pos, g = act_coverage$Sigma, p.adjust.method = "bonferroni",
                     paired = FALSE, conf.int = FALSE, alternative = 't')

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  act_coverage$uni_str_pos and act_coverage$Sigma 
## 
##     0       0.1     0.3     0.6     1.2     2.5
## 0.1 1       -       -       -       -       -  
## 0.3 1       1       -       -       -       -  
## 0.6 1       1       1       -       -       -  
## 1.2 1       1       1       1       -       -  
## 2.5 2.7e-12 7.4e-13 6.8e-12 2.1e-09 9.8e-11 -  
## 5   1.5e-13 2.0e-14 4.5e-13 1.1e-10 1.8e-12 1  
## 
## P value adjustment method: bonferroni

4.4.3 Satisfactory trait coverage over time

Satisfactory trait coverage in a population over time. Data points on the graph is the average activation gene coverage across 50 replicates every 2000 generations. Shading comes from the best and worse coverage across 50 replicates.

lines = filter(over_time_df, acro == 'con')  %>%
  group_by(Sigma, gen) %>%
  dplyr::summarise(
    min = min(pop_uni_obj),
    mean = mean(pop_uni_obj),
    max = max(pop_uni_obj)
  )

## `summarise()` has grouped output by 'Sigma'. You can override using the
## `.groups` argument.

ggplot(lines, aes(x=gen, y=mean, group = Sigma, fill = Sigma, color = Sigma, shape = Sigma)) +
  geom_ribbon(aes(ymin = min, ymax = max), alpha = 0.1) +
  geom_line(size = 0.5) +
  geom_point(data = filter(lines, gen %% 2000 == 0 & gen != 0), size = 1.5, stroke = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Coverage",
    limits=c(0, 2),
    breaks=c(0,1,2)
  ) +
  scale_x_continuous(
    name="Generations",
    limits=c(0, 50000),
    breaks=c(0, 10000, 20000, 30000, 40000, 50000),
    labels=c("0e+4", "1e+4", "2e+4", "3e+4", "4e+4", "5e+4")

  ) +
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Satisfactory trait coverage over time')+
  p_theme +
  guides(
    shape=guide_legend(nrow=1, title.position = "left", title = 'Sigma'),
    color=guide_legend(nrow=1, title.position = "left", title = 'Sigma'),
    fill=guide_legend(nrow=1, title.position = "left", title = 'Sigma')
  )

4.4.4 Final satisfactory trait coverage

Satisfactory trait coverage found in the final population at 50,000 generations.

plot = filter(over_time_df, gen == 50000 & acro == 'con') %>%
  ggplot(., aes(x = Sigma, y = pop_uni_obj, color = Sigma, fill = Sigma, shape = Sigma)) +
  geom_flat_violin(position = position_nudge(x = .1, y = 0), scale = 'width', alpha = 0.2, width = 1.5) +
  geom_boxplot(color = 'black', width = .07, outlier.shape = NA, alpha = 0.0, size = 1.0, position = position_nudge(x = .16, y = 0)) +
  geom_point(position = position_jitter(width = 0.03, height = 0.02), size = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Coverage",
    limits=c(0, 2),
    breaks=c(0,1,2)
  ) +
  scale_x_discrete(
    name="Sigma"
  )+
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette, ) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Final satisfactory trait coverage')+
  p_theme

plot_grid(
  plot +
    theme(legend.position="none"),
  legend,
  nrow=2,
  rel_heights = c(3,1)
)

## Warning: Removed 181 rows containing missing values (`geom_point()`).

4.4.4.1 Stats

Summary statistics for the generation a satisfactory solution is found.

sat_coverage = filter(over_time_df, gen == 50000 & acro == 'con')
sat_coverage$Sigma = factor(sat_coverage$Sigma, levels = FS_LIST)
sat_coverage %>%
  group_by(Sigma) %>%
  dplyr::summarise(
    count = n(),
    na_cnt = sum(is.na(pop_uni_obj)),
    min = min(pop_uni_obj, na.rm = TRUE),
    median = median(pop_uni_obj, na.rm = TRUE),
    mean = mean(pop_uni_obj, na.rm = TRUE),
    max = max(pop_uni_obj, na.rm = TRUE),
    IQR = IQR(pop_uni_obj, na.rm = TRUE)
  )

## # A tibble: 7 x 8
##   Sigma count na_cnt   min median  mean   max   IQR
##   <fct> <int>  <int> <int>  <dbl> <dbl> <int> <dbl>
## 1 0        50      0     0      0     0     0     0
## 2 0.1      50      0     0      0     0     0     0
## 3 0.3      50      0     0      0     0     0     0
## 4 0.6      50      0     0      0     0     0     0
## 5 1.2      50      0     0      0     0     0     0
## 6 2.5      50      0     0      0     0     0     0
## 7 5        50      0     0      0     0     0     0

4.4.5 Largest valley reached throughout

The largest valley reached in a single trait throughout an entire evolutionary run. To collect this value, we look through all the best-performing solutions each generation and find the largest valley reached.

plot = filter(best_df, acro == 'con' & var == 'ele_big_peak') %>%
  ggplot(., aes(x = Sigma, y = val, color = Sigma, fill = Sigma, shape = Sigma)) +
  geom_flat_violin(position = position_nudge(x = .1, y = 0), scale = 'width', alpha = 0.2, width = 1.5) +
  geom_boxplot(color = 'black', width = .07, outlier.shape = NA, alpha = 0.0, size = 1.0, position = position_nudge(x = .16, y = 0)) +
  geom_point(position = position_jitter(width = 0.03, height = 0.02), size = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Valley reached",
    limits=c(4.9,10.1)
  ) +
  scale_x_discrete(
    name="Sigma"
  )+
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette, ) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Largest valley reached')+
  p_theme

plot_grid(
  plot +
    theme(legend.position="none"),
  legend,
  nrow=2,
  rel_heights = c(3,1)
)

4.4.5.1 Stats

Summary statistics for the largest valley crossed.

valleys = filter(best_df, acro == 'con' & var == 'ele_big_peak')
valleys$Sigma = factor(valleys$Sigma, levels = FS_LIST)
valleys %>%
  group_by(Sigma) %>%
  dplyr::summarise(
    count = n(),
    na_cnt = sum(is.na(val)),
    min = min(val, na.rm = TRUE),
    median = median(val, na.rm = TRUE),
    mean = mean(val, na.rm = TRUE),
    max = max(val, na.rm = TRUE),
    IQR = IQR(val, na.rm = TRUE)
  )

## # A tibble: 7 x 8
##   Sigma count na_cnt   min median  mean   max   IQR
##   <fct> <int>  <int> <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 0        50      0     6      6  6.2      7  0   
## 2 0.1      50      0     6      6  6.28     7  1   
## 3 0.3      50      0     6      6  6.26     7  0.75
## 4 0.6      50      0     5      6  6.26     7  1   
## 5 1.2      50      0     6      6  6.36     8  1   
## 6 2.5      50      0     7      8  8.36    10  1   
## 7 5        50      0     8      9  8.6     10  1

Kruskal–Wallis test illustrates evidence of statistical differences.

kruskal.test(val ~ Sigma, data = valleys)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  val by Sigma
## Kruskal-Wallis chi-squared = 248.66, df = 6, p-value < 2.2e-16

Results for post-hoc Wilcoxon rank-sum test with a Bonferroni correction.

pairwise.wilcox.test(x = valleys$val, g = valleys$Sigma, p.adjust.method = "bonferroni",
                     paired = FALSE, conf.int = FALSE, alternative = 't')

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  valleys$val and valleys$Sigma 
## 
##     0      0.1    0.3    0.6    1.2    2.5
## 0.1 1      -      -      -      -      -  
## 0.3 1      1      -      -      -      -  
## 0.6 1      1      1      -      -      -  
## 1.2 1      1      1      1      -      -  
## 2.5 <2e-16 <2e-16 <2e-16 <2e-16 <2e-16 -  
## 5   <2e-16 <2e-16 <2e-16 <2e-16 <2e-16 1  
## 
## P value adjustment method: bonferroni

4.5 Multi-path exploration results

Here we present the results for best performances and activation gene coverage found by each selection scheme parameter on the multi-path exploration diagnostic with valleys. 50 replicates are conducted for each scheme parameter explored.

4.5.1 Activation gene coverage over time

Activation gene coverage in a population over time. Data points on the graph is the average activation gene coverage across 50 replicates every 2000 generations. Shading comes from the best and worse coverage across 50 replicates.

lines = filter(over_time_df, acro == 'mpe')  %>%
  group_by(Sigma, gen) %>%
  dplyr::summarise(
    min = min(uni_str_pos),
    mean = mean(uni_str_pos),
    max = max(uni_str_pos)
  )

## `summarise()` has grouped output by 'Sigma'. You can override using the
## `.groups` argument.

ggplot(lines, aes(x=gen, y=mean, group = Sigma, fill = Sigma, color = Sigma, shape = Sigma)) +
  geom_ribbon(aes(ymin = min, ymax = max), alpha = 0.1) +
  geom_line(size = 0.5) +
  geom_point(data = filter(lines, gen %% 2000 == 0 & gen != 0), size = 1.5, stroke = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Coverage",
    limits=c(0, 100),
    breaks=seq(0,100, 20),
    labels=c("0", "20", "40", "60", "80", "100")
  ) +
  scale_x_continuous(
    name="Generations",
    limits=c(0, 50000),
    breaks=c(0, 10000, 20000, 30000, 40000, 50000),
    labels=c("0e+4", "1e+4", "2e+4", "3e+4", "4e+4", "5e+4")

  ) +
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Activation gene coverage over time')+
  p_theme +
  guides(
    shape=guide_legend(nrow=1, title.position = "left", title = 'Sigma'),
    color=guide_legend(nrow=1, title.position = "left", title = 'Sigma'),
    fill=guide_legend(nrow=1, title.position = "left", title = 'Sigma')
  )

4.5.2 Final activation gene coverage

Activation gene coverage found in the final population at 50,000 generations.

plot = filter(over_time_df, gen == 50000 & acro == 'mpe') %>%
  ggplot(., aes(x = Sigma, y = uni_str_pos, color = Sigma, fill = Sigma, shape = Sigma)) +
  geom_flat_violin(position = position_nudge(x = .1, y = 0), scale = 'width', alpha = 0.2, width = 1.5) +
  geom_boxplot(color = 'black', width = .07, outlier.shape = NA, alpha = 0.0, size = 1.0, position = position_nudge(x = .16, y = 0)) +
  geom_point(position = position_jitter(width = 0.03, height = 0.02), size = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Coverage",
    limits=c(0, 6.1),
    breaks=seq(0,6,2)
  ) +
  scale_x_discrete(
    name="Sigma"
  )+
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette, ) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Final activation gene coverage')+
  p_theme

plot_grid(
  plot +
    theme(legend.position="none"),
  legend,
  nrow=2,
  rel_heights = c(3,1)
)

4.5.2.1 Stats

Summary statistics for the generation a satisfactory solution is found.

act_coverage = filter(over_time_df, gen == 50000 & acro == 'mpe')
act_coverage$Sigma = factor(act_coverage$Sigma, levels = FS_LIST)
act_coverage %>%
  group_by(Sigma) %>%
  dplyr::summarise(
    count = n(),
    na_cnt = sum(is.na(uni_str_pos)),
    min = min(uni_str_pos, na.rm = TRUE),
    median = median(uni_str_pos, na.rm = TRUE),
    mean = mean(uni_str_pos, na.rm = TRUE),
    max = max(uni_str_pos, na.rm = TRUE),
    IQR = IQR(uni_str_pos, na.rm = TRUE)
  )

## # A tibble: 7 x 8
##   Sigma count na_cnt   min median  mean   max   IQR
##   <fct> <int>  <int> <int>  <dbl> <dbl> <int> <dbl>
## 1 0        50      0     1      3  2.8      4     1
## 2 0.1      50      0     1      3  2.86     5     1
## 3 0.3      50      0     2      3  2.9      4     1
## 4 0.6      50      0     2      3  3.08     6     2
## 5 1.2      50      0     2      3  3.38     5     1
## 6 2.5      50      0     3      4  4.14     6     2
## 7 5        50      0     3      5  4.76     6     1

Kruskal–Wallis test illustrates evidence of statistical differences.

kruskal.test(uni_str_pos ~ Sigma, data = act_coverage)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  uni_str_pos by Sigma
## Kruskal-Wallis chi-squared = 138.2, df = 6, p-value < 2.2e-16

Results for post-hoc Wilcoxon rank-sum test with a Bonferroni correction.

pairwise.wilcox.test(x = act_coverage$uni_str_pos, g = act_coverage$Sigma, p.adjust.method = "bonferroni",
                     paired = FALSE, conf.int = FALSE, alternative = 't')

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  act_coverage$uni_str_pos and act_coverage$Sigma 
## 
##     0       0.1     0.3     0.6     1.2     2.5   
## 0.1 1.0000  -       -       -       -       -     
## 0.3 1.0000  1.0000  -       -       -       -     
## 0.6 1.0000  1.0000  1.0000  -       -       -     
## 1.2 0.0128  0.0347  0.0682  0.7788  -       -     
## 2.5 1.7e-08 6.2e-08 7.8e-08 1.1e-05 0.0012  -     
## 5   4.8e-14 1.6e-13 8.5e-14 4.3e-11 9.6e-11 0.0129
## 
## P value adjustment method: bonferroni

4.5.3 Performance over time

Best performance in a population over time. Data points on the graph is the average performance across 50 replicates every 2000 generations. Shading comes from the best and worse performance across 50 replicates.

lines = filter(over_time_df, acro == 'mpe') %>%
  group_by(Sigma, gen) %>%
  dplyr::summarise(
    min = min(pop_fit_max) / DIMENSIONALITY,
    mean = mean(pop_fit_max) / DIMENSIONALITY,
    max = max(pop_fit_max) / DIMENSIONALITY
  )

## `summarise()` has grouped output by 'Sigma'. You can override using the
## `.groups` argument.

over_time_plot = ggplot(lines, aes(x=gen, y=mean, group = Sigma, fill = Sigma, color = Sigma, shape = Sigma)) +
  geom_ribbon(aes(ymin = min, ymax = max), alpha = 0.1) +
  geom_line(size = 0.5) +
  geom_point(data = filter(lines, gen %% 2000 == 0 & gen != 0), size = 1.5, stroke = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Average trait score",
    limits=c(0, 15)
  ) +
  scale_x_continuous(
    name="Generations",
    limits=c(0, 50000),
    breaks=c(0, 10000, 20000, 30000, 40000, 50000),
    labels=c("0e+4", "1e+4", "2e+4", "3e+4", "4e+4", "5e+4")

  ) +
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Performance over time')+
  p_theme +
  guides(
    shape=guide_legend(nrow=1, title.position = "left", title = 'Sigma'),
    color=guide_legend(nrow=1, title.position = "left", title = 'Sigma'),
    fill=guide_legend(nrow=1, title.position = "left", title = 'Sigma')
  )

over_time_plot

4.5.4 Best performance throughout

Best performance reached throughout 50,000 generations in a population.

plot = filter(best_df, acro == 'mpe' &  var == 'pop_fit_max') %>%
  ggplot(., aes(x = Sigma, y = val / DIMENSIONALITY, color = Sigma, fill = Sigma, shape = Sigma)) +
  geom_flat_violin(position = position_nudge(x = .1, y = 0), scale = 'width', alpha = 0.2, width = 1.5) +
  geom_boxplot(color = 'black', width = .07, outlier.shape = NA, alpha = 0.0, size = 1.0, position = position_nudge(x = .16, y = 0)) +
  geom_point(position = position_jitter(width = 0.03, height = 0.02), size = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Average trait score",
    limits=c(1.9, 14.1),
    breaks=c(2,5,8,11,14)
  ) +
  scale_x_discrete(
    name="Sigma"
  )+
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette, ) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Best performance throughout')+
  p_theme

plot_grid(
  plot +
    theme(legend.position="none"),
  legend,
  nrow=2,
  rel_heights = c(3,1)
)

4.5.4.1 Stats

Summary statistics for the best performance.

performance = filter(best_df, acro == 'mpe' &  var == 'pop_fit_max')
performance$Sigma = factor(performance$Sigma, levels =FS_LIST)
performance %>%
  group_by(Sigma) %>%
  dplyr::summarise(
    count = n(),
    na_cnt = sum(is.na(val)),
    min = min(val / DIMENSIONALITY, na.rm = TRUE),
    median = median(val / DIMENSIONALITY, na.rm = TRUE),
    mean = mean(val / DIMENSIONALITY, na.rm = TRUE),
    max = max(val / DIMENSIONALITY, na.rm = TRUE),
    IQR = IQR(val / DIMENSIONALITY, na.rm = TRUE)
  )

## # A tibble: 7 x 8
##   Sigma count na_cnt   min median  mean   max   IQR
##   <fct> <int>  <int> <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 0        50      0  1.94   9.90  9.12  11.8 2.03 
## 2 0.1      50      0  2.83  11.2  10.5   13.0 0.846
## 3 0.3      50      0  2.88  11.3  10.7   14.0 1.69 
## 4 0.6      50      0  3.25  11.6  11.5   13.5 1.29 
## 5 1.2      50      0  3.20  12.1  11.8   13.1 0.870
## 6 2.5      50      0  6.68  11.3  11.1   12.5 0.677
## 7 5        50      0  6.85  11.2  10.8   12.9 1.07

Kruskal–Wallis test illustrates evidence of statistical differences.

kruskal.test(val ~ Sigma, data = performance)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  val by Sigma
## Kruskal-Wallis chi-squared = 80.554, df = 6, p-value = 2.745e-15

Results for post-hoc Wilcoxon rank-sum test with a Bonferroni correction.

pairwise.wilcox.test(x = performance$val, g = performance$Sigma, p.adjust.method = "bonferroni",
                     paired = FALSE, conf.int = FALSE, alternative = 't')

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  performance$val and performance$Sigma 
## 
##     0       0.1     0.3     0.6     1.2     2.5    
## 0.1 0.00157 -       -       -       -       -      
## 0.3 0.00067 1.00000 -       -       -       -      
## 0.6 1.6e-08 0.16193 1.00000 -       -       -      
## 1.2 6.8e-12 8.0e-05 0.01514 1.00000 -       -      
## 2.5 2.3e-06 1.00000 1.00000 0.58011 8.3e-06 -      
## 5   0.00075 1.00000 1.00000 0.08218 3.1e-06 1.00000
## 
## P value adjustment method: bonferroni

4.5.5 Largest valley reached throughout

The largest valley reached in a single trait throughout an entire evolutionary run. To collect this value, we look through all the best-performing solutions each generation and find the largest valley reached.

plot = filter(best_df, acro == 'mpe' & var == 'ele_big_peak') %>%
  ggplot(., aes(x = Sigma, y = val, color = Sigma, fill = Sigma, shape = Sigma)) +
  geom_flat_violin(position = position_nudge(x = .1, y = 0), scale = 'width', alpha = 0.2, width = 1.5) +
  geom_boxplot(color = 'black', width = .07, outlier.shape = NA, alpha = 0.0, size = 1.0, position = position_nudge(x = .16, y = 0)) +
  geom_point(position = position_jitter(width = 0.03, height = 0.02), size = 2.0, alpha = 1.0) +
  scale_y_continuous(
    name="Valley reached",
    limits=c(7.9,14.1)
  ) +
  scale_x_discrete(
    name="Sigma"
  )+
  scale_shape_manual(values=SHAPE)+
  scale_colour_manual(values = cb_palette, ) +
  scale_fill_manual(values = cb_palette) +
  ggtitle('Largest valley reached')+
  p_theme

plot_grid(
  plot +
    theme(legend.position="none"),
  legend,
  nrow=2,
  rel_heights = c(3,1)
)

4.5.5.1 Stats

Summary statistics for the largest valley crossed.

valleys = filter(best_df, acro == 'mpe' & var == 'ele_big_peak')
valleys$Sigma = factor(valleys$Sigma, levels = FS_LIST)
valleys %>%
  group_by(Sigma) %>%
  dplyr::summarise(
    count = n(),
    na_cnt = sum(is.na(val)),
    min = min(val, na.rm = TRUE),
    median = median(val, na.rm = TRUE),
    mean = mean(val, na.rm = TRUE),
    max = max(val, na.rm = TRUE),
    IQR = IQR(val, na.rm = TRUE)
  )

## # A tibble: 7 x 8
##   Sigma count na_cnt   min median  mean   max   IQR
##   <fct> <int>  <int> <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 0        50      0     8     12  11.5    13     1
## 2 0.1      50      0    10     13  12.5    14     1
## 3 0.3      50      0    11     13  12.8    14     0
## 4 0.6      50      0    11     13  13.1    14     0
## 5 1.2      50      0    12     14  13.6    14     1
## 6 2.5      50      0    13     14  13.9    14     0
## 7 5        50      0    14     14  14      14     0

Kruskal–Wallis test illustrates evidence of statistical differences.

kruskal.test(val ~ Sigma, data = valleys)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  val by Sigma
## Kruskal-Wallis chi-squared = 238.55, df = 6, p-value < 2.2e-16

Results for post-hoc Wilcoxon rank-sum test with a Bonferroni correction.

pairwise.wilcox.test(x = valleys$val, g = valleys$Sigma, p.adjust.method = "bonferroni",
                     paired = FALSE, conf.int = FALSE, alternative = 't')

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  valleys$val and valleys$Sigma 
## 
##     0       0.1     0.3     0.6     1.2     2.5    
## 0.1 2.2e-06 -       -       -       -       -      
## 0.3 4.4e-10 0.19650 -       -       -       -      
## 0.6 2.1e-13 0.00022 0.88310 -       -       -      
## 1.2 6.9e-16 4.6e-10 1.2e-06 0.00092 -       -      
## 2.5 < 2e-16 5.9e-15 1.2e-12 9.9e-10 0.07000 -      
## 5   < 2e-16 < 2e-16 < 2e-16 2.1e-14 1.4e-05 0.13628
## 
## P value adjustment method: bonferroni