We want to plot the heterogeneity properties of the two correlation methods depending on their parameters. Let’s compute the values.

heterogeneity <- NULL
CV <- 1
M <- 1e-3
breaks <- pnorm(seq(qnorm(M), -qnorm(M), length.out = 100))
for (rhoR in breaks)
  for (rhoC in breaks) {
    N1 <- -sqrt((rhoR - rhoC) ^ 2*CV ^ 4 +
                  2*(rhoR*(rhoC - 1) ^ 2 + rhoC*(rhoR - 1) ^ 2)*CV ^ 2 + (rhoC*rhoR - 1) ^ 2)
    N2 <- -(2*rhoR*rhoC - rhoR - rhoC)*CV ^ 2 - rhoC*rhoR + 1
    D <- 2*rhoR*rhoC*(CV ^ 2 + 1)
    Vnoise <- sqrt((N1 + N2)/D)
    Vtask <- 1 / sqrt(1 / Vnoise ^ 2 * (1 / rhoC - 1) - 1)
    Vmach <- 1 / sqrt(1 / Vnoise ^ 2 * (1 / rhoR - 1) - 1)
    Vmutask <- Vtask
    Vmumach <- Vmach
    muVtask <- sqrt(Vtask ^ 2*Vnoise ^ 2 + Vtask ^ 2 + Vnoise ^ 2)
    muVmach <- sqrt(Vmach ^ 2*Vnoise ^ 2 + Vmach ^ 2 + Vnoise ^ 2)
    heterogeneity <- rbind(heterogeneity,
                           data.frame(rhoR = rhoR, rhoC = rhoC,
                                      method = "CNB",
                                      Vmutask = Vmutask, Vmumach = Vmumach,
                                      muVtask = muVtask, muVmach = muVmach))
    Psi <- CV/(sqrt(rhoR*(1 - rhoC)) + sqrt(1 - rhoR)*(sqrt(rhoC) + sqrt(1 - rhoC)))
    Vmutask <- sqrt((1 - rhoR)*rhoC)*Psi
    Vmumach <- sqrt((1 - rhoC)*rhoR)*Psi
    muVtask <- sqrt(1 - rhoR)*Psi
    muVmach <- sqrt(1 - rhoC)*Psi
    heterogeneity <- rbind(heterogeneity,
                           data.frame(rhoR = rhoR, rhoC = rhoC,
                                      method = "CB",
                                      Vmutask = Vmutask, Vmumach = Vmumach,
                                      muVtask = muVtask, muVmach = muVmach))
  }

Let’s plot them:

heterogeneity_long <- gather(heterogeneity, measure, hetero, matches("Vmu|muV")) %>%
  filter(measure == "Vmutask" | measure == "muVtask")

levels(heterogeneity_long$method) <- c("Correlation~noise-based", "Combination-based")
levels(heterogeneity_long$measure) <- c(expression(V ~ mu[task]), expression(V ~ mu[mach]),
                                        expression(mu ~ V[task]), expression(mu ~ V[mach]))

library(scales)
probit_trans <- trans_new("probit", transform = qnorm, inverse = pnorm)

breaks_axis <- c(0.01, 0.1, 0.5, 0.9, 0.99)
breaks_fill <- c(0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1)
ggplot(heterogeneity_long, aes(x = rhoR, y = rhoC, fill = hetero, z = hetero)) +
  scale_x_continuous(name = expression(r[task]), trans = probit_trans,
                     breaks = breaks_axis) +
  scale_y_continuous(name = expression(r[mach]), trans = probit_trans,
                     breaks = breaks_axis) +
  facet_grid(measure ~ method, labeller = label_parsed) +
  geom_tile() +
  scale_fill_gradientn(name = "", colours = rainbow(100), trans = "log",
                       breaks = breaks_fill) +
  scale_color_continuous(low = "black", high = "black", guide = "none") +
  stat_contour(aes(color = ..level..), breaks = breaks_fill)