bean_plot

inst/paper/cbean.R

@@ -2,6 +2,7 @@ library("climatrends")
 library("tidyverse")
 library("PlackettLuce")
 library("patchwork")
+library("qvcalc")
 library("ggparty")
 
 
@@ -30,96 +31,94 @@ temp <- temperature(modis,
 
 # round values to 3 decimals, this will not change the main result, but will 
 # help in visualizing the tree 
-temp[1:ncol(temp)] <- lapply(temp[1:ncol(temp)], function(x) round(x, 3))
+temp[1:ncol(temp)] <- lapply(temp[1:ncol(temp)], function(x) round(x, 2))
 
 # combine the indices with the main data
 cbean <- cbind(cbean, temp)
 
 # fit a PL tree
-plt <- pltree(G ~ maxNT, data = cbean, minsize = 200)
+plt <- pltree(G ~ maxNT + SU, data = cbean, minsize = 50)
+
+plot(plt)
 
 # ..................................
 # ..................................
 # Prepare charts ####
 A <- 
 ggplot(cbean, aes(gdd)) +
   geom_density() +
-  labs(x = "GDD", y = "Density", title = "A") +
+  labs(x = "GDD", y = "Density", title = "(a)") +
   theme_bw() + 
   theme(panel.background = element_blank(), 
         panel.grid = element_blank(),
         plot.title = element_text(size=16, 
                                   colour = "grey20", 
                                   face = "bold"),
-        axis.text.x = element_text(size = 12, angle = 0,
+        axis.text.x = element_text(size = 10, angle = 0,
                                    face="plain", colour = "grey20"),
-        axis.title.x = element_text(size=12, colour = "grey20"),
-        axis.text.y = element_text(size=12, angle = 0,
+        axis.title.x = element_text(size=10, colour = "grey20"),
+        axis.text.y = element_text(size=10, angle = 0,
                                    hjust=1, vjust=0.5,
                                    face="plain", colour = "grey20"),
-        axis.title.y = element_text(size=12, colour = "grey20"))
+        axis.title.y = element_text(size=10, colour = "grey20"))
 
 B <- 
 ggplot(cbean, aes(season, maxNT)) +
-  geom_violin(aes(fill = season), show.legend = FALSE) +
+  geom_violin(aes(fill = season), alpha = 0.5, show.legend = FALSE) +
   #geom_jitter(height = 0, width = 0.1) +
   scale_fill_manual(values = c('#d7191c','#fdae61','#ffffbf','#abd9e9','#2c7bb6')) +
-  labs(x = "", y = "maxNT", title = "B") +
+  labs(x = "", y = "maxNT", title = "(b)") +
   theme_bw() + 
   theme(panel.background = element_blank(), 
         panel.grid = element_blank(),
         plot.title = element_text(size=16, 
                                   colour = "grey20", 
                                   face = "bold"),
-        axis.text.x = element_text(size = 12, angle = 0,
+        axis.text.x = element_text(size = 10, angle = 0,
                                    face="plain", colour = "grey20"),
-        axis.title.x = element_text(size=12, colour = "grey20"),
-        axis.text.y = element_text(size=12, angle = 0,
+        axis.title.x = element_text(size=10, colour = "grey20"),
+        axis.text.y = element_text(size=10, angle = 0,
                                    hjust=1, vjust=0.5,
                                    face="plain", colour = "grey20"),
-        axis.title.y = element_text(size=12, colour = "grey20"))
-
+        axis.title.y = element_text(size=10, colour = "grey20"))
 
+C <-
+ggparty(plt) +
+  geom_edge() +
+  labs(title = "(c)") +
+  geom_edge_label() +
+  geom_node_splitvar() +
+  theme(plot.title = element_text(size=16, 
+                                  colour = "grey20", 
+                                  face = "bold"))
 
 coeff <- as.vector(coef(plt, log = FALSE))
-items <- rep(dimnames(coef(plt))[[2]], each = 2)
+items <- rep(dimnames(coef(plt))[[2]], each = dim(coef(plt))[[1]])
 nit <- length(unique(items))
-
-coeff <- data.frame(coeff = coeff, items = items, 
-                    node = rep(c("Node 2","Node 3"),
+coeff <- data.frame(estimate = coeff,
+                    items = items,
+                    node = rep(paste("Node", nodeids(plt, terminal = TRUE)),
                                times = length(items)))
 
 D <-
-ggplot(coeff, aes(x = coeff, y = items, group = node)) +
+  ggplot(coeff, aes(x = estimate, y = items, group = node)) +
   geom_vline(xintercept = 1/nit, 
              colour = "#E5E7E9", size = 0.8) +
   geom_point() +
   facet_wrap(~ node) +
+  scale_x_continuous(labels = c(".08",".10",".12",".14")) +
   labs(x = "", y="") +
   theme_bw() + 
   theme(panel.background = element_blank(), 
         strip.background = element_rect(fill = "white"),
         panel.grid = element_blank(),
-        plot.title = element_text(size=16, 
-                                  colour = "grey20", 
-                                  face = "bold"),
-        axis.text.x = element_text(size = 12, angle = 50, vjust = 0.5,
-                                   face="plain", colour = "grey20"),
-        axis.title.x = element_text(size=12, colour = "grey20"),
-        axis.text.y = element_text(size=12, angle = 0,
-                                   hjust=1, vjust=0.5,
+        axis.text.x = element_text(size = 8, vjust = 0.5,
+                                   face ="plain", colour = "grey20"),
+        axis.title.x = element_text(size = 10, colour = "grey20"),
+        axis.text.y = element_text(size = 10, angle = 0,
+                                   hjust = 1, vjust=0.5,
                                    face="plain", colour = "grey20"),
-        axis.title.y = element_text(size=12, colour = "grey20"))
-
-C <-
-ggparty(plt) +
-  geom_edge() +
-  labs(title = "C") +
-  geom_edge_label() +
-  geom_node_splitvar() +
-  theme(plot.title = element_text(size=16, 
-                                  colour = "grey20", 
-                                  face = "bold"))
+        axis.title.y = element_text(size=10, colour = "grey20"))
 
 p <- A / B | (C / D)