diff --git a/boxplot_with_stats.R b/boxplot_with_stats.R
index cdc34b1..7736230 100755
--- a/boxplot_with_stats.R
+++ b/boxplot_with_stats.R
@@ -9,8 +9,7 @@ getwd()
 #=============
 # Input
 #=============
-source("boxplot_linear.R")
-
+#source("boxplot_linear.R")
 #######################################################
 med_names = c("eotaxin3", "il12p70", "itac", "il13")
 lf_test = lf_fp_npa[lf_fp_npa$mediator%in%med_names,]
@@ -38,7 +37,11 @@ p2 = ggplot(single)+ geom_boxplot(aes(x = timepoint
         , axis.title.x = element_blank()
         , axis.title.y = element_blank()
         , legend.position = "none"
-        , plot.title = element_text(size = 20, hjust = 0.5))
+        , plot.subtitle = element_text(size = 20, hjust = 0.5)
+        , plot.title = element_text(size = 20, hjust = 0.5)) +
+  labs(title = i
+       #, subtitle = "test2"
+       )
 
 stat_npa2 <- single %>%
   group_by(timepoint, mediator) %>%
@@ -50,7 +53,7 @@ stat_npa2 <- stat_npa2 %>%
   add_xy_position(x = "timepoint", dodge = 0.8)
 
 p2 = p2 + stat_pvalue_manual(stat_npa2
-                             , y.position = max_y
+                             #, y.position = max_y
                              , label = "{p} {p.signif}"
                              , hide.ns=T
                              , tip.length = 0)+
@@ -59,3 +62,91 @@ plots[[i]] <- p2
 }
 
 cowplot::plot_grid(plotlist=plots, align = 'hv', ncol=2, nrow=2)
+##################################################################
+#=======
+# facet
+#=======
+#--------
+# wilcox stats
+#--------
+stat_npa3 <- lf_test %>%
+  group_by(timepoint, mediator) %>%
+  wilcox_test(value ~ obesity, paired = F) %>%
+  add_significance("p")
+stat_npa3 
+
+stat_npa3 <- stat_npa3 %>%
+  add_xy_position(x = "timepoint", dodge = 0.8)
+head(stat_npa3)
+
+#--------
+# summary stats
+#--------
+my_summary = lf_test %>%
+  group_by(timepoint, mediator) %>%
+  get_summary_stats(value)
+my_summary
+my_max_df = subset(my_summary
+                   , select = c("mediator", "timepoint"
+                                #, "obesity"
+                                , "max"))
+
+head(my_max_df); head(stat_npa3)
+
+#-----------------------------
+# merge my_max and stat_npa3
+#-----------------------------
+head(my_max_df)
+merging_cols = intersect(names(stat_npa3), names(my_max_df)); merging_cols
+
+stat_npa3_v2 = merge(stat_npa3, my_max_df, by = merging_cols, all.x = T)
+stat_npa3_v2$my_y_pos = (stat_npa3_v2$max)*1.2
+stat_npa3_v2$my_y_pos
+
+head(stat_npa3_v2); head(stat_npa3)
+
+p3 = ggplot(lf_test)+ geom_boxplot(aes(x = timepoint
+                                      , y = value
+                                      , color = obesity))+
+  theme(axis.text.x = element_text(size = 15)
+        , axis.text.y = element_text(size = 15
+                                     , angle = 0
+                                     , hjust = 1
+                                     , vjust = 0)
+        , axis.title.x = element_blank()
+        , axis.title.y = element_blank()
+        , legend.position = "none"
+        , plot.subtitle = element_text(size = 20, hjust = 0.5)
+        , plot.title = element_text(size = 20, hjust = 0.5)) +
+  labs(title ="NPA")+
+  #facet_wrap(~mediator, scales = "free")
+
+p3_build = ggbuild(p3)  
+  
+p4 = p3 + stat_pvalue_manual(stat_npa3_v2
+                             , y.position = "my_y_pos"
+                             #, y.position = 50
+                             , label = "{p} {p.signif}"
+                             , hide.ns=T
+                             , tip.length = 0)+
+  scale_y_continuous(expand = expansion(mult = c(0.05, 0.25)))
+
+p4
+#========================================
+p3_build = ggplot_build(p3)
+p3_build$layout$panel_scales_y
+
+get_facet_ymax <- function (x){
+  ret <- x$layout$panel_scales_y
+  y_max =  NULL
+  
+  for (i in 1:length(ret)){
+    print(i)
+    y_max_i <- max(x$layout$panel_scales_y[[i]]$range$range)
+    y_max = c(y_max, y_max_i)
+  }
+  
+  return(y_max)
+}
+
+y_max_l = get_facet_ymax(p3_build); y_max_l