Update analysis-cassava.R

dev/analysis-cassava.R

@@ -2,18 +2,25 @@ library("gosset")
 library("ClimMobTools")
 library("PlackettLuce")
 library("ggplot2")
-library(multcompView)
+library("patchwork")
+library("multcompView")
 
 data("cassava", package = "gosset")
 
-dat = cassava[cassava$country == "Nigeria", ]
+#dat = cassava[cassava$country == "Nigeria", ]
+
+dat = cassava
 
 keep = unlist(lapply(dat[1:ncol(dat)], function(x) sum(is.na(x))))
 
-keep = keep != nrow(dat)
+keep = keep == 0
+
+keep
 
 dat = dat[, keep]
 
+names(dat)
+
 pack = c("option_a", "option_b", "option_c")
 
 items = sort(unique(unlist(dat[pack])))
@@ -28,6 +35,12 @@ traits = gsub("(^|[[:space:]])([[:alpha:]])", "\\1\\U\\2", traits, perl = TRUE)
 
 ov = which(traits %in% "Overall")
 
+table(unlist(dat[pack]), rep(dat$gender, 3))
+
+table(unlist(dat[pack]), rep(dat$country, 3))
+
+country = sort(unique(dat$country))
+
 R = lapply(trait_list, function(x) {
   rank_tricot(dat, 
               items = pack, 
@@ -39,26 +52,52 @@ lapply(R, function(x){
   likelihood_ratio(x, split = dat$gender)
 })
 
+lapply(R, function(x){
+  likelihood_ratio(x, split = dat$country)
+})
 
-mod = lapply(R, PlackettLuce)
 
-worth_map(mod, 
-          labels = traits) +
-  labs(x = "Genotype", y = "Trait") +
-  scale_fill_distiller(palette = "BrBG", 
-                       direction = 1, 
-                       na.value = "white", 
-                       name = "")
+G = lapply(trait_list, function(x) {
+  rank_tricot(dat, 
+              items = pack, 
+              input = x$string,
+              validate.rankings = TRUE,
+              group = TRUE)
+})
 
 
-plot_logworth(mod[[ov]], 
-              ref = "Akpu",
-              levels = rev(items),
-              multcomp = T)
+trait_plot = list()
+
+for (i in seq_along(country)) {
+
+  mod_i = lapply(R, function(x) {
+    PlackettLuce(x[dat$country == country[i], ])
+  })
+
+  trait_plot[[i]] = worth_map(mod_i, 
+            labels = traits) +
+    labs(x = "", 
+         y = "",
+         title = country[i]) +
+    scale_fill_distiller(palette = "BrBG", 
+                         direction = 1, 
+                         na.value = "white", 
+                         name = "")
+
+
+}
 
 
+trait_plot[[1]] + trait_plot[[2]]
 
 
+
+plot_logworth(mod[[ov]], 
+              ref = "Akpu",
+              levels = rev(items),
+              multcomp = T,
+              ci.level = 0.05)
+
 rel = reliability(mod[[ov]], ref = "Akpu")
 
 rel