From 260f73b136f4a68461dc80374b31778560ded1e9 Mon Sep 17 00:00:00 2001 From: kauedesousa Date: Thu, 22 Jun 2023 18:02:40 -0400 Subject: [PATCH 1/2] fix-issue-in-handling-empty-logworth-plots --- ClimMob.R | 30 +++++---------- modules/03_organize_quantitative_data.R | 9 ++++- modules/06_PlackettLuce_models.R | 50 ++++++++++++------------- 3 files changed, 40 insertions(+), 49 deletions(-) diff --git a/ClimMob.R b/ClimMob.R index fbe989b..66f7363 100755 --- a/ClimMob.R +++ b/ClimMob.R @@ -161,7 +161,8 @@ if (any_error(org_lonlat)) { org_agroclim = tryCatch({ agroclimate = get_agroclimatic_data(cmdata, - coords = trial_map$coords) + coords = trial_map$coords, + ndays = 60) }, error = function(cond) { return(cond) @@ -333,7 +334,7 @@ dir.create(chartdir, recursive = TRUE, showWarnings = FALSE) # log worth plot by trait for(m in seq_along(PL_models$logworth_plot)){ - try(ggsave(paste0(chartdir, rank_dat$trait_code[m], "_logworth.png"), + try(ggsave(paste0(chartdir, m, "-logworth.pdf"), plot = PL_models$logworth_plot[[m]], width = 21, height = 15, @@ -342,7 +343,7 @@ for(m in seq_along(PL_models$logworth_plot)){ } # plot kendall tau plot -try(ggsave(paste0(chartdir, "kendall_tau.png"), +try(ggsave(paste0(chartdir, "kendall_tau.pdf"), plot = PL_models$kendall$kendall_plot, width = 15, height = 18, @@ -350,7 +351,7 @@ try(ggsave(paste0(chartdir, "kendall_tau.png"), dpi = 200), silent = TRUE) # plot worth map -try(ggsave(paste0(chartdir, "worth_map.png"), +try(ggsave(paste0(chartdir, "worth_map.pdf"), plot = PL_models$worthmap, width = 25, height = 25, @@ -358,7 +359,7 @@ try(ggsave(paste0(chartdir, "worth_map.png"), dpi = 200), silent = TRUE) # plot worth map -try(ggsave(paste0(chartdir, "reliability.png"), +try(ggsave(paste0(chartdir, "reliability.pdf"), plot = PL_models$reliability_plot, width = 25, height = 25, @@ -368,21 +369,8 @@ try(ggsave(paste0(chartdir, "reliability.png"), try(write.csv(PL_models$reliability_data, paste0(chartdir, "reliability_data.csv"), row.names = FALSE), silent = TRUE) -if(length(unique(rank_dat$group)) > 1) { - g = unique(rank_dat$group) - # log worth plot by group - for(m in seq_along(PL_models$logworth_plot_groups)){ - try(ggsave(paste0(chartdir, "Group", m, "_", g[m], "_logworth_grouped_rank.png"), - plot = PL_models$logworth_plot_groups[[m]], - width = 21, - height = 15, - units = "cm", - dpi = 200), silent = TRUE) - } -} - if(PL_tree$isTREE){ - try(ggsave(paste0(chartdir, "PlackettLuce.png"), + try(ggsave(paste0(chartdir, "PlackettLuceTree.pdf"), plot = PL_tree$PLtree_plot, width = 18, height = 25, @@ -401,14 +389,14 @@ if (isTRUE(agroclimate$agroclimate)) { file = paste0(chartdir, "weekly_temperature_indices.csv"), row.names = FALSE) - try(ggsave(paste0(chartdir, "weekly_precipitation_indices.png"), + try(ggsave(paste0(chartdir, "weekly_precipitation_indices.pdf"), plot = agroclimate$rain_plot, width = 20, height = 20, units = "cm", dpi = 200), silent = TRUE) - try(ggsave(paste0(chartdir, "weekly_temperature_indices.png"), + try(ggsave(paste0(chartdir, "weekly_temperature_indices.pdf"), plot = agroclimate$temperature_plot, width = 20, height = 20, diff --git a/modules/03_organize_quantitative_data.R b/modules/03_organize_quantitative_data.R index 72a818a..99f73e1 100644 --- a/modules/03_organize_quantitative_data.R +++ b/modules/03_organize_quantitative_data.R @@ -36,6 +36,13 @@ organize_quantitative_data = function(cmdata, id = "id", tech_index = c("package_item_A", "package_item_B", "package_item_C")) { + + quanti_traits = pars[["linear"]] + + if (length(quanti_traits) == 0) { + return(list(quantitative = FALSE)) + } + # from json to data.frame cmdata = as.data.frame(x = cmdata, tidynames = FALSE, @@ -43,8 +50,6 @@ organize_quantitative_data = function(cmdata, ntech = length(tech_index) - quanti_traits = pars[["linear"]] - # check if a request to split the data by groups (segments) # (gender, location, etc.) is provided if (isTRUE(length(groups) > 0)) { diff --git a/modules/06_PlackettLuce_models.R b/modules/06_PlackettLuce_models.R index 4f1f6b7..73f01c8 100644 --- a/modules/06_PlackettLuce_models.R +++ b/modules/06_PlackettLuce_models.R @@ -36,9 +36,9 @@ get_PlackettLuce_models = function(cmdata, rank_dat) { # first a list with rankings R = lapply(trait_list, function(x){ rank_tricot(cmdata[, c(technologies_index, x$strings)], - items = technologies_index, - input = x$strings, - validate.rankings = TRUE) + items = technologies_index, + input = x$strings, + validate.rankings = TRUE) }) @@ -51,8 +51,8 @@ get_PlackettLuce_models = function(cmdata, rank_dat) { logworth_plot = list() for(m in seq_along(mod)) { - logworth_plot[[m]] = - plot_logworth(mod[[m]], ref = reference_tech, ci.level = 0.5) + + lwp = + try(plot_logworth(mod[[m]], ref = reference_tech, ci.level = 0.5) + labs(title = paste0(rank_dat$trait_names[m], " (n = ", length(mod[[m]]$rankings),")")) + coord_flip() + @@ -63,9 +63,25 @@ get_PlackettLuce_models = function(cmdata, rank_dat) { strip.placement = "outside", strip.text = element_text(size = 10, color = "grey20"), legend.text = element_text(size = 10, color = "grey20"), - axis.title = element_text(size = 10, color = "grey20")) + axis.title = element_text(size = 10, color = "grey20")), + silent = TRUE) + + if ("try-error" %in% class(lwp)) next + + logworth_plot[[m]] = lwp + + } + + logworth_overall = logworth_plot[[reference_trait_index]] + + if ("try-error" %in% class(logworth_overall)) { + logworth_overall = ggplot() } + rmv = unlist(lapply(logworth_plot, is.null)) + + logworth_plot = logworth_plot[!rmv] + #........................................................... # Head to head visualization of each technology performance by trait worth_map_data = data.frame() @@ -177,14 +193,6 @@ get_PlackettLuce_models = function(cmdata, rank_dat) { }) - # for(j in seq_along(trait_names)) { - # - # if (is.null(rel_i[[j]])) next - # - # rel_i[[j]]$Trait = trait_names[j] - # - # } - rel_i = do.call("rbind", rel_i) reliability_data = rbind(reliability_data, rel_i) @@ -193,12 +201,7 @@ get_PlackettLuce_models = function(cmdata, rank_dat) { # remove the checks reliability_data = reliability_data[!reliability_data$item %in% checks, ] - - # # put traits in the right order - # reliability_data$Trait = factor(reliability_data$Trait, - # levels = rev(trait_names)) - - + # plot the reliability reliability_plot = ggplot(data = reliability_data, @@ -212,7 +215,6 @@ get_PlackettLuce_models = function(cmdata, rank_dat) { linewidth = 1) + scale_fill_manual(values = "#b2df8a") + facet_wrap(~ Check, strip.position = "bottom") + - #theme_bw() + theme_classic() + theme(panel.grid.major = element_blank(), strip.background =element_rect(fill="white"), @@ -225,11 +227,7 @@ get_PlackettLuce_models = function(cmdata, rank_dat) { axis.title = element_text(size = 12, color = "grey20")) + labs(y = "Probability of outperforming", x = "") - - - reliability_plot - #..................................................... #..................................................... #..................................................... # Get the Kendall tau #### @@ -313,7 +311,7 @@ get_PlackettLuce_models = function(cmdata, rank_dat) { # export results result = list(PL_models = mod, - logworth_overall = logworth_plot[[reference_trait_index]], + logworth_overall = logworth_overall, worthmap = worthmap, ANOVA = aov_tbl, logworth_plot = logworth_plot, From d657f2b4d587d00d24cede9142f11ce54aba5a9d Mon Sep 17 00:00:00 2001 From: kauedesousa Date: Thu, 25 Jan 2024 15:37:03 +0100 Subject: [PATCH 2/2] new-version --- NEWS.md | 13 +++++ modules/01_functions.R | 84 +++++++++++++++++++++++++++----- modules/05_spatial_overview.R | 11 ++--- modules/06_PlackettLuce_models.R | 23 +++++---- modules/08_PlackettLuce_tree.R | 1 - report/mainreport.Rmd | 4 +- 6 files changed, 101 insertions(+), 35 deletions(-) diff --git a/NEWS.md b/NEWS.md index 7778acc..91ec913 100755 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,16 @@ +ClimMob-analysis v3.0 (2024-01-25) +========================= + +### Improvements + +* Adds analysis of variance for variety performance +* Adds pseudo ranking when network is poorly connected + +### Bug fixes + +* Fixes changes in reference for the log-worth plot + + ClimMob-analysis v2.1 (2022-12-15) ========================= diff --git a/modules/01_functions.R b/modules/01_functions.R index 76e8891..273d0a3 100755 --- a/modules/01_functions.R +++ b/modules/01_functions.R @@ -35,6 +35,59 @@ library("janitor") library("lubridate") library("ggchicklet") + +#' Add pseudo ranking +#' Adds pseudo values to weakly connected networks +#' @param x a PlackettLuce ranking object +force_pseudo_rank = function(x) { + + # get membership in the network + members = PlackettLuce::connectivity(x)$membership + # put the members in order + members = sort(members) + # rankings into a matrix + r = unclass(x) + + performance = coefficients(PlackettLuce::PlackettLuce(r)) + + # get the worst item per cluster + members = split(members, members) + + members = lapply(members, function(z){ + p = performance[names(z)] + worst = which.min(p) + names(p)[worst] + }) + + members = as.character(unlist(members)) + + # create a pseudo ranking for these members where they will always + # lose and win to each other + # number to start ranking + max_rank = max(r) + 1 + + # rows to add the pseudo rankings + to_input = rowSums(r) != 0 + + to_sample = c(rep(0, ceiling(length(members)/2)), + max_rank:(max_rank + ceiling(length(members)/2))) + + # rows to add the pseudo rankings + r[to_input, members] = apply(r[to_input, members], 1, function(y){ + + where = y == 0 + + y[where] = sample(to_sample, size = length(y[where])) + + y + + }) + + r = as.rankings(r) + +} + + #'Get colour pallet #' @param x an integer #' @examples @@ -577,6 +630,12 @@ decode_pars = function(x) { tr = toupper(tr) } + if (any(grepl("generalappreciation", tr))) { + i = which(grepl("generalappreciation", tr))[1] + questions$traitOrder[i] = "referenceTrait" + tr = toupper(tr) + } + if (any(grepl("yield", tr))) { i = which(grepl("yield", tr))[1] questions$traitOrder[i] = "referenceTrait" @@ -794,7 +853,8 @@ multcompPL = function(mod, items = NULL, threshold = 0.05, adjust = "none", ...) #' @param ci.level the confidence interval level #' @param multcomp logical to add group letters #' @param levels an optional vector with factor levels to plot -plot_logworth = function(x, ci.level = 0.95, ref = NULL, multcomp = TRUE, levels = NULL, ...) { +plot_logworth = function(x, ci.level = 0.95, ref = NULL, + multcomp = TRUE, levels = NULL, ...) { frame = data.frame() @@ -806,7 +866,7 @@ plot_logworth = function(x, ci.level = 0.95, ref = NULL, multcomp = TRUE, levels } if (is.null(levels)) { - levels = unique(frame$items) + levels = union(ref, sort(unique(frame$items))) } items = factor(frame$items, levels = levels) @@ -838,29 +898,29 @@ plot_logworth = function(x, ci.level = 0.95, ref = NULL, multcomp = TRUE, levels pdat$items = factor(pdat$items, levels = levels) p = ggplot(data = pdat, - aes(x = items, - y = est, - ymax = tops, - ymin = tails, + aes(y = items, + x = est, + xmax = tops, + xmin = tails, label = group)) + - geom_hline(yintercept = 0, + geom_vline(xintercept = 0, colour = "#E5E7E9", linewidth = 0.8) + geom_point() + geom_errorbar(width = 0.1) + - geom_text(vjust = 1.2, hjust = 1.2) + + geom_text(hjust = 1.2, vjust = 1.2) + theme_bw() + facet_wrap(~ ref, strip.position = "bottom") + theme(panel.grid.major = element_blank(), - axis.text.x = element_text(angle = 45, vjust = 1, hjust=1, - size = 10, color = "grey20"), + strip.background.x = element_blank(), axis.text.y = element_text(size = 10, color = "grey20"), + axis.text.x = element_text(size = 10, color = "grey20"), text = element_text(color = "grey20"), legend.position = "bottom", legend.title = element_blank(), - strip.background.x = element_blank(), + strip.background.y = element_blank(), strip.placement = "outside") + - labs(x = "", y = "Log-worth") + labs(y = "", x = "Log-worth") p diff --git a/modules/05_spatial_overview.R b/modules/05_spatial_overview.R index 61a4bdc..87bea98 100644 --- a/modules/05_spatial_overview.R +++ b/modules/05_spatial_overview.R @@ -60,20 +60,15 @@ get_testing_sites_map = function(cmdata, output_path, backward_path){ minimap = TRUE, map_provider = "OpenStreetMap.Mapnik") - tempmap = paste0(getwd(), "/tempmap/") - - dir.create(tempmap, recursive = TRUE, showWarnings = FALSE) - try(mapview::mapshot(trial_map, - url = paste0(tempmap, "/trial_map.html"), - file = paste0(tempmap, "/trial_map.png")), + url = paste0(output_path, "/trial_map.html"), + file = paste0(output_path, "/trial_map.png")), silent = TRUE) } result = list(geoTRUE = TRUE, - mapDIR = tempmap, - map_path = paste0(tempmap, "/trial_map.png"), + map_path = paste0(fullpath, "/", output_path, "/trial_map.png"), map = trial_map, coords = lonlat) diff --git a/modules/06_PlackettLuce_models.R b/modules/06_PlackettLuce_models.R index 73f01c8..23aabcf 100644 --- a/modules/06_PlackettLuce_models.R +++ b/modules/06_PlackettLuce_models.R @@ -41,7 +41,14 @@ get_PlackettLuce_models = function(cmdata, rank_dat) { validate.rankings = TRUE) }) + # # Handle rankings with poor connectivity + connection = lapply(R, function(x) connectivity(x, verbose = FALSE)$no) + connection = as.vector(unlist(connection) > 1) + + # force a pseudo rank + R[connection] = lapply(R[connection], force_pseudo_rank) + # fit the model mod = lapply(R, function(x){ PlackettLuce(x) @@ -51,19 +58,11 @@ get_PlackettLuce_models = function(cmdata, rank_dat) { logworth_plot = list() for(m in seq_along(mod)) { - lwp = - try(plot_logworth(mod[[m]], ref = reference_tech, ci.level = 0.5) + + lwp = try(plot_logworth(mod[[m]], + ref = reference_tech, + ci.level = 0.5) + labs(title = paste0(rank_dat$trait_names[m], - " (n = ", length(mod[[m]]$rankings),")")) + - coord_flip() + - theme(axis.text.x = element_text(angle = 0, - vjust = 0.5, - hjust = 0.5), - strip.background.x = element_blank(), - strip.placement = "outside", - strip.text = element_text(size = 10, color = "grey20"), - legend.text = element_text(size = 10, color = "grey20"), - axis.title = element_text(size = 10, color = "grey20")), + " (n = ", length(mod[[m]]$rankings),")")), silent = TRUE) if ("try-error" %in% class(lwp)) next diff --git a/modules/08_PlackettLuce_tree.R b/modules/08_PlackettLuce_tree.R index d20ffd8..07a61dc 100644 --- a/modules/08_PlackettLuce_tree.R +++ b/modules/08_PlackettLuce_tree.R @@ -151,7 +151,6 @@ get_PlackettLuce_tree = function(cmdata, rank_dat, agroclimate = NULL) { validate.rankings = TRUE, group = TRUE) - #.......................................................... # PlackettLuce tree #### # data frame of explanatory variables diff --git a/report/mainreport.Rmd b/report/mainreport.Rmd index 11a9379..7ad765f 100755 --- a/report/mainreport.Rmd +++ b/report/mainreport.Rmd @@ -29,9 +29,9 @@ This report presents the results from your tricot experiment entitled `r project # What you should expect from this report -This report will provide to you an overview of your tricot experiment data and results, with insights of the performance of the technologies tested in the tricot trial. This is a standard automated report that tries to accommodate the most important outputs from a tricot trial. After reading through this report, you should be able to: (i) identify the technology that outperforms the others in your trial, and under which conditions, and (ii) support a decision on advancement of a technology to the next stage of your product development program (e.g. breeding program, market assessment). +This report provides the results from your tricot experiment, with insights of the performance of the technologies tested in the trial. This is a standard automated report that accommodates the most important outputs from a tricot trial. After reading through this report, you should be able to: (i) identify the technology that outperforms the others in your trial, and under which conditions, and (ii) support a decision on advancement of a technology to the next stage of your product development program (e.g. breeding program, market assessment). -However, since this report was done in an automated process, is likely that it will not provide to you the full overview on the performance of the tested technologies as the automation may not capture all the factors that influence the performance of a technology in the target environment (under real-world conditions). If you would like to have an in-deep analysis, or to merge data from different projects or is planning a peer-review publication with your tricot data, send a proposal/request to Kauê de Sousa (k.desousa@cgiar.org) and Jacob van Etten (j.vanetten@cgiar.org). Even being implemented with high-quality statistical analysis, the automated process for the production of this report limits the ability to test specific hypothesis with the data to match the needs of international peer-review journals. Therefore, we recommend to extend the analysis beyond the results described here. +However, since this report was done in an automated process, is likely that it will not provide to you the full overview on the performance of the tested technologies as the automation may not capture all the factors that influence the performance of a technology in the target environment (under real-world conditions). If you would like to have an in-deep analysis, or to merge data from different projects or is planning a peer-review publication with your tricot data, send a proposal/request to the ClimMob Team . # Methods