Repeatability for stability parameters

repstab = function(data, boot = 3, nenv = 16, initial = 3, strata = FALSE){

  if(nenv < 6){
    stop("O número de ambientes deve ser maior que seis, para que estatisticas como o desvio da regressão possam ser calculadas para cada subconjunto de ambientes")
  }

sets = nenv/2
nenv = (nenv/2-initial)+1


Ncorr=(10*(10-1))/2

names = c("mean", "R2", "RMSE", "sdij", "CV", "Ecovalencia",
          "Shukla", "WAASB", "PC1", "WAASY")
combnam = combn(names,2, paste, collapse = " x ")

temp_cor_1 = data.frame(matrix(nrow = boot, ncol = Ncorr + 2))
names(temp_cor_1) = c("NENV", "BOOT", names(sapply(combnam, names)))
temp_cor_2 = data.frame(matrix(nrow = boot, ncol = Ncorr + 2))
names(temp_cor_2) = c("NENV", "BOOT", names(sapply(combnam, names)))

temp = data.frame(matrix(nrow = boot, ncol = 12))
names(temp) = c("NENV", "BOOT", names)

results_rand_cor1 = list()
results_rand_cor2 = list()
results_random = list()


trigoup = data %>% group_by(ENV) %>%
  dplyr::summarise(mean = mean(RG)) %>%
  arrange(desc(mean)) %>%
  slice(1:sets)
trigoudown = data %>% group_by(ENV) %>%
  dplyr::summarise(mean = mean(RG)) %>%
  arrange(mean) %>%
  slice(1:sets)

for (k in 1:nenv){
  pb <- winProgressBar(title = "Repeatability by resampling approach",
                       label = "processing data",
                       min = 0,
                       width = 400,
                       max = boot,
                       initial = 0)
for (i in 1:boot){

  if (strata == TRUE){

    AMB = sample(unique(trigoup$ENV), initial, replace = FALSE)
    AMB2 = sample(unique(trigoudown$ENV), initial, replace = FALSE)

  } else{
    AMB = sample(unique(data$ENV), initial, replace = FALSE)
    AMB2 = unique(data$ENV)[-c(AMB)]
    AMB2 = sample(AMB2, initial, replace = FALSE)
  }

data1 = data %>% filter(ENV %in% AMB) %>% as.data.frame()
data1$ENV = factor(data1$ENV, levels = AMB)
data2 = data %>% filter(ENV %in% AMB2) %>% as.data.frame()
data2$ENV = factor(data2$ENV, levels = AMB2)

data1_stat = ge_stats(data1, env = "ENV", gen = "CULTIVAR", rep = "REP", resp = "RG")
data2_stat = ge_stats(data2, env = "ENV", gen = "CULTIVAR", rep = "REP", resp = "RG")

stab_1 = cbind(data1_stat[["ER"]][["Regression"]] %>%
         select(Gen, Mean, R.Sqr, RMSE, sdij),
         data1_stat[["stab_meas"]] %>%
           select(CV, Ecov, ShuklaVar))

stab_2 = cbind(data2_stat[["ER"]][["Regression"]] %>%
                 select(Gen, Mean, R.Sqr, RMSE, sdij),
               data2_stat[["stab_meas"]] %>%
                 select(CV, Ecov, ShuklaVar))

stab_1_rank = cbind(
  stab_1 %>% mutate(mean_r = rank(-Mean),
                    R.Sqr_r = rank(-R.Sqr),
                    RMSE_r = rank(RMSE),
                    sdij_r = rank(abs(sdij)),
                    CV_r = rank(CV),
                    Ecov_r = rank(Ecov),
                    ShuklaVar_r = rank(ShuklaVar)) %>%
    select(c(mean_r, R.Sqr_r, RMSE_r, sdij_r, CV_r, Ecov_r, ShuklaVar_r)),

  WAASB(data1, resp = "RG")[["model"]] %>%
    subset(type == "GEN", select = c(OrWAASB, OrPC1, OrWAASY))
)

stab_2_rank = cbind(
             stab_2 %>% mutate(mean_r = rank(-Mean),
                                R.Sqr_r = rank(-R.Sqr),
                                RMSE_r = rank(RMSE),
                                sdij_r = rank(abs(sdij)),
                                CV_r = rank(CV),
                                Ecov_r = rank(Ecov),
                                ShuklaVar_r = rank(ShuklaVar)) %>%
  select(c(mean_r, R.Sqr_r, RMSE_r, sdij_r, CV_r, Ecov_r, ShuklaVar_r)),

  WAASB(data2, resp = "RG")[["model"]] %>%
    subset(type == "GEN", select = c(OrWAASB, OrPC1, OrWAASY))
)

setWinProgressBar(pb, i, title=paste("Bootstrapping", i," of ", boot, "replicates ", initial, "envs - ",
                                        round(i/boot*100,2),"% Concluded -"))


temp[i,1] = initial
temp[i,2] = i
temp[i,3:12] = diag(cor(stab_1_rank, stab_2_rank, method = "spearman"))

temp_cor_1[i,1] = initial
temp_cor_1[i,2] = i
temp_cor_1[i,3:ncol(temp_cor_1)] = as.vector(t(cor(stab_1_rank, method = "spearman"))
                                 [lower.tri(cor(stab_1_rank, method = "spearman"),diag=F)])

temp_cor_2[i,1] = initial
temp_cor_2[i,2] = i
temp_cor_2[i,3:ncol(temp_cor_2)] = as.vector(t(cor(stab_2_rank, method = "spearman"))
                                 [lower.tri(cor(stab_2_rank, method = "spearman"),diag=F)])
}

results_random[[paste(initial, "ENV")]] = temp
results_rand_cor1[[paste(initial, "ENV")]] = temp_cor_1
results_rand_cor2[[paste(initial, "ENV")]] = temp_cor_2

initial = initial + 1
close(pb)
}

return(list(results_random = ldply(results_random),
            results_rand_cor1 = ldply(results_rand_cor1),
            results_rand_cor2 = ldply(results_rand_cor2)))
}