add(admin game init): implement pheno params in game stups script

timflutre · Jul 23, 2024 · d75fdbb · d75fdbb
1 parent 45035c8
commit d75fdbb
Show file tree

Hide file tree

Showing 6 changed files with 408 additions and 37 deletions.
diff --git a/initialise_data.R b/initialise_data.R
@@ -2,6 +2,7 @@
 
 params = list(
   rng_seed = 1993,
+
   cost.pheno.field = 50,
   cost.pheno.patho = 0.1,
   cost.allof = 0.1,
@@ -11,7 +12,20 @@ params = list(
   cost.geno.ld = 0.5,
   cost.geno.single = 0.02,
   cost.register = 4,
-  initialBudget = 3900
+  initialBudget = 3900,
+
+  t1_mu = 100,
+  t1_min = 20,
+  t1_cv_g = 0.1,
+  t1_h2 = 0.3,
+
+  t2_mu = 15,
+  t2_min = 5,
+  t2_cv_g = 0.06,
+  t2_h2 = 0.6,
+
+  prop_pleio = 0.4,
+  cor_pleio = -0.7
 )
 
 

diff --git a/src/fun/func_gameInit_validation.R b/src/fun/func_gameInit_validation.R
@@ -100,15 +100,39 @@ valid_number <- function(x, accept_null = TRUE, raise_error = FALSE) {
     error("Mandatory and should be a number")
   }
 
+  if (!is.numeric(x)) {
+    error("Mandatory and should be a number")
+  }
+
   return(NULL)
 }
 
 valid_range <- function(x, min, max, incl_min = TRUE, incl_max = TRUE, accept_null = TRUE, raise_error = FALSE) {
 
+  error <- return
   error <- return
   if (raise_error) {
     error <- stop
   }
+  error <- return
+  if (raise_error) {
+    error <- stop
+  }
+  error <- return
+  if (raise_error) {
+    error <- stop
+  }
+  error <- return
+  if (raise_error) {
+    error <- stop
+  }
+  error <- return
+  if (raise_error) {
+    error <- stop
+  }
+  if (raise_error) {
+    error <- stop
+  }
 
   if (is.null(x)) {
     if (accept_null) {
@@ -121,6 +145,10 @@ valid_range <- function(x, min, max, incl_min = TRUE, incl_max = TRUE, accept_nu
     error("Mandatory and should be a number")
   }
 
+  if (!is.numeric(x)) {
+    error("Mandatory and should be a number")
+  }
+
   min_msg <- paste(min, "(excluded)")
   if (incl_min) {
     min_msg <- paste(min, "(included)")
@@ -151,7 +179,31 @@ valid_range <- function(x, min, max, incl_min = TRUE, incl_max = TRUE, accept_nu
 }
 
 valid_mu <- function(x, accept_null = TRUE, raise_error = FALSE) {
-  valid_number(x, accept_null, raise_error)
+  error <- return
+  if (raise_error) {
+    error <- stop
+  }
+
+  if (is.null(x)) {
+    if (accept_null) {
+      return(NULL)
+    }
+    error("Must not be NULL")
+  }
+
+  if (is.na(x)) {
+    error("Mandatory and should be a number")
+  }
+
+  if (!is.numeric(x)) {
+    error("Mandatory and should be a number")
+  }
+
+  if (x == 0) {
+    error("μ cannot be null")
+  }
+
+  return(NULL)
 }
 
 valid_Tmin <- function(x, mu, accept_null = TRUE, raise_error = FALSE) {
@@ -171,6 +223,10 @@ valid_Tmin <- function(x, mu, accept_null = TRUE, raise_error = FALSE) {
     error("Mandatory and should be a number")
   }
 
+  if (!is.numeric(x)) {
+    error("Mandatory and should be a number")
+  }
+
   if (is.na(mu)) {
     return(NULL)
   }
@@ -228,6 +284,10 @@ valid_variance <- function(x, name = NULL, accept_na = FALSE, accept_null = TRUE
     error(paste0("variance is NA", var_name))
   }
 
+  if (!is.numeric(x)) {
+    error(paste0("variance is not a number", var_name))
+  }
+
   if (x <= 0) {
     error(paste0("This value leads to negative variance", var_name))
   }

diff --git a/src/fun/module_gameInit_params.R b/src/fun/module_gameInit_params.R
@@ -663,10 +663,10 @@ gameInit_traits_server <- function(id, iv) {
       }
       input$t1_mu
     })
-    output$T1_sig_p2 <- renderText({round(T1_sig_p2(), 4)})
-    output$T1_sig_a2 <- renderText({round(T1_sig_a2(), 4)})
-    output$T1_sig2 <- renderText({round(T1_sig2(), 4)})
-    output$T1_sig_y2 <- renderText({round(T1_sig_y2(), 4)})
+    output$T1_sig_p2 <- renderText({signif(T1_sig_p2(), 6)})
+    output$T1_sig_a2 <- renderText({signif(T1_sig_a2(), 6)})
+    output$T1_sig2 <- renderText({signif(T1_sig2(), 6)})
+    output$T1_sig_y2 <- renderText({signif(T1_sig_y2(), 6)})
 
     output$T2_mu <- renderText({
       id = session$ns("prev_t2_mu")
@@ -677,10 +677,10 @@ gameInit_traits_server <- function(id, iv) {
       }
       input$t2_mu
     })
-    output$T2_sig_p2 <- renderText({round(T2_sig_p2(), 4)})
-    output$T2_sig_a2 <- renderText({round(T2_sig_a2(), 4)})
-    output$T2_sig2 <- renderText({round(T2_sig2(), 4)})
-    output$T2_sig_y2 <- renderText({round(T2_sig_y2(), 4)})
+    output$T2_sig_p2 <- renderText({signif(T2_sig_p2(), 6)})
+    output$T2_sig_a2 <- renderText({signif(T2_sig_a2(), 6)})
+    output$T2_sig2 <- renderText({signif(T2_sig2(), 6)})
+    output$T2_sig_y2 <- renderText({signif(T2_sig_y2(), 6)})
 
 
     output$pheno_plot <- renderPlotly({
@@ -766,11 +766,32 @@ gameInit_traits_server <- function(id, iv) {
 
       return(fig)
     })
-  })
-}
-
-
-
 
+    return(
+      list(
+        value = reactive({
+          if (pheno_params_validator$is_valid()) {
+            return(list(
+              t1_mu = input$t1_mu,
+              t1_min = input$t1_min,
+              t1_cv_g = input$t1_cv_g,
+              t1_h2 = input$t1_h2,
+
+              t2_mu = input$t2_mu,
+              t2_min = input$t2_min,
+              t2_cv_g = input$t2_cv_g,
+              t2_h2 = input$t2_h2,
+
+              prop_pleio = input$prop_pleio,
+              cor_pleio = input$cor_pleio
+            ))
+          }
+          return(NA)
+        }),
+        iv = iv
+      )
+    )
 
+  })
+}
 
diff --git a/src/plantbreedgame_setup.Rmd b/src/plantbreedgame_setup.Rmd
@@ -15,6 +15,16 @@ params:
   cost.geno.single: 0.02
   cost.register: 4
   initialBudget: 3900
+  t1_mu: 20
+  t1_min: 100
+  t1_cv_g: 0.1
+  t1_h2: 0.3
+  t2_mu: 15
+  t2_min: 5
+  t2_cv_g: 0.06
+  t2_h2: 0.6
+  prop_pleio: 0.4
+  cor_pleio: -0.7
 colorlinks: true
 output:
   html_document:
@@ -140,6 +150,38 @@ invisible({
   valid_positive_number(params$cost.geno.single, accept_null = FALSE, raise_error = TRUE)
   valid_positive_number(params$cost.register, accept_null = FALSE, raise_error = TRUE)
   valid_positive_number(params$initialBudget, accept_null = FALSE, raise_error = TRUE)
+
+  valid_mu(params$t1_mu, accept_null = FALSE, raise_error = TRUE)
+  valid_Tmin(params$t1_min, params$t1_mu, accept_null = FALSE, raise_error = TRUE)
+  valid_cv_g(params$t1_cv_g, accept_null = FALSE, raise_error = TRUE)
+  valid_h2(params$t1_h2, accept_null = FALSE, raise_error = TRUE)
+
+  valid_mu(params$t2_mu, accept_null = FALSE, raise_error = TRUE)
+  valid_Tmin(params$t2_min, params$t2_mu, accept_null = FALSE, raise_error = TRUE)
+  valid_cv_g(params$t2_cv_g, accept_null = FALSE, raise_error = TRUE)
+  valid_h2(params$t2_h2, accept_null = FALSE, raise_error = TRUE)
+
+  t1_sig_p2 <- calc_sigma_p2(params$t1_mu, params$t1_min)
+  valid_variance(t1_sig_p2, "phenotypes T1", accept_na = FALSE, accept_null = FALSE, raise_error = TRUE)
+  t1_sig_a2 <- calc_sigma_a2(params$t1_cv_g, params$t1_mu)
+  valid_variance(t1_sig_a2, "genetic values T1", accept_na = FALSE, accept_null = FALSE, raise_error = TRUE)
+  t1_sig2 <- calc_sigma2(params$t1_h2, t1_sig_a2)
+  valid_variance(t1_sig2, "noise T1", accept_na = FALSE, accept_null = FALSE, raise_error = TRUE)
+  t1_sig_y2 <- calc_sigma_y2(t1_sig_p2, t1_sig_a2, t1_sig2 )
+  valid_variance(t1_sig_y2, "year effects T1", accept_na = FALSE, accept_null = FALSE, raise_error = TRUE)
+
+  t2_sig_p2 <- calc_sigma_p2(params$t2_mu, params$t2_min)
+  valid_variance(t2_sig_p2, "phenotypes T2", accept_na = FALSE, accept_null = FALSE, raise_error = TRUE)
+  t2_sig_a2 <- calc_sigma_a2(params$t2_cv_g, params$t2_mu)
+  valid_variance(t2_sig_a2, "genetic values T2", accept_na = FALSE, accept_null = FALSE, raise_error = TRUE)
+  t2_sig2 <- calc_sigma2(params$t2_h2, t2_sig_a2)
+  valid_variance(t2_sig2, "noise T2", accept_na = FALSE, accept_null = FALSE, raise_error = TRUE)
+  t2_sig_y2 <- calc_sigma_y2(t2_sig_p2, t2_sig_a2, t2_sig2 )
+  valid_variance(t2_sig_y2, "year effects T2", accept_na = FALSE, accept_null = FALSE, raise_error = TRUE)
+
+  valid_prop_pleio(params$prop_pleio, accept_null = FALSE, raise_error = TRUE)
+  valid_cor_pleio(params$cor_pleio, accept_null = FALSE, raise_error = TRUE)
+
 })
 
 print(params[names(params) != "progressBar"])
@@ -742,13 +784,13 @@ Here is the procedure for trait 1 of setting the variances required by the gener
 We set the initial phenotypic mean, $\mu_p$, to a global intercept which will be held constant during the whole game: $\mu_p = \mu$.
 
 ```{r simul_proc_t1_mu}
-(mu <- 100)
+(mu <- params$t1_mu)
 ```
 
 We set the individual-plant phenotypic variance so that approximately 99% of all the phenotypes are above a minimal value: $\mu - 3 \sigma_p \ge y_{\text{min}}$.
 
 ```{r simul_proc_t1_sigma_p2}
-y.min <- 20
+y.min <- params$t1_min
 (sigma.p2 <- ((mu - y.min) / 3)^2)
 sqrt(sigma.p2)
 ```
@@ -777,7 +819,7 @@ We deduce the genotypic coef of variation:
 
 However, because the previous computations don't account for the negative correlation between traits, we still fix the CVg by hand at this point:
 ```{r}
-CV.g <- 0.10
+CV.g <- params$t1_cv_g
 ```
 
 We deduce the (additive) genotypic variance:
@@ -791,7 +833,7 @@ stopifnot(sigma.a2 < sigma.p2)
 We deduce the error variance given that the intra-annual, individual-plant heritability is set at a particular value: $h^2 = \frac{\sigma_a^2}{\sigma_a^2 + \sigma^2}$.
 
 ```{r simul_proc_t1_sigma2}
-h2 <- 0.3
+h2 <- params$t1_h2
 (sigma2 <- ((1 - h2) / h2) * sigma.a2)
 sqrt(sigma2)
 stopifnot(sigma2 + sigma.a2 < sigma.p2)
@@ -839,13 +881,13 @@ Here is the procedure for trait 2 of setting the variances required by the gener
 We fix the initial phenotypic mean, $\mu_p$, to a global intercept which will be held constant during the whole game: $\mu_p = \mu$.
 
 ```{r simul_proc_t2_mu}
-(mu <- 15)
+(mu <- params$t2_mu)
 ```
 
 We deduce the individual-plant phenotypic variance so that approximately 99% of all the phenotypes are above a minimal value: $\mu - 3 \sigma_p \ge y_{\text{min}}$.
 
 ```{r simul_proc_t2_sigma_p2}
-y.min <- 5
+y.min <- params$t2_min
 (sigma.p2 <- ((mu - y.min) / 3)^2)
 sqrt(sigma.p2)
 ```
@@ -859,7 +901,7 @@ register.min.trait2 <- mu
 We deduce the (additive) genotypic variance given that the genotypic coef of variation is set at a particular value:
 
 ```{r simul_proc_t2_sigma_a2}
-CV.g <- 0.06
+CV.g <- params$t2_cv_g
 (sigma.a2 <- (CV.g * mu)^2)
 sqrt(sigma.a2)
 stopifnot(sigma.a2 < sigma.p2)
@@ -868,7 +910,7 @@ stopifnot(sigma.a2 < sigma.p2)
 We deduce the error variance given that the intra-annual, individual-plant heritability is set at a particular value: $h^2 = \frac{\sigma_a^2}{\sigma_a^2 + \sigma^2}$.
 
 ```{r simul_proc_t2_sigma2}
-h2 <- 0.6
+h2 <- params$t2_h2
 (sigma2 <- ((1 - h2) / h2) * sigma.a2)
 sqrt(sigma2)
 stopifnot(sigma2 + sigma.a2 < sigma.p2)
@@ -920,7 +962,7 @@ Simulate the SNP effects for traits 1 and 2 jointly:
 snp.effects12 <- simulSnpEffectsTraits12(
   snp.ids = thin.snps,
   sigma.beta2 = sigma.beta2,
-  prop.pleio = 0.4, cor.pleio = -0.7
+  prop.pleio = params$prop_pleio, cor.pleio = params$cor_pleio
 )
 ```
 

diff --git a/src/server/server_admin.R b/src/server/server_admin.R
@@ -776,6 +776,7 @@ observeEvent(input$initialiseGame, {
   params <- list(
     progressBar = progress_bar,
     rng_seed = gameInit_seed$value(),
+
     cost.pheno.field = gameInit_costs$value()$cost.pheno.field,
     cost.pheno.patho = gameInit_costs$value()$cost.pheno.patho,
     cost.allof = gameInit_costs$value()$cost.allof,
@@ -785,9 +786,24 @@ observeEvent(input$initialiseGame, {
     cost.geno.ld = gameInit_costs$value()$cost.geno.ld,
     cost.geno.single = gameInit_costs$value()$cost.geno.single,
     cost.register = gameInit_costs$value()$cost.register,
-    initialBudget = gameInit_costs$value()$initialBudget
+    initialBudget = gameInit_costs$value()$initialBudget,
+
+    t1_mu = gameInit_traits$value()$t1_mu,
+    t1_min = gameInit_traits$value()$t1_min,
+    t1_cv_g = gameInit_traits$value()$t1_cv_g,
+    t1_h2 = gameInit_traits$value()$t1_h2,
+
+    t2_mu = gameInit_traits$value()$t2_mu,
+    t2_min = gameInit_traits$value()$t2_min,
+    t2_cv_g = gameInit_traits$value()$t2_cv_g,
+    t2_h2 = gameInit_traits$value()$t2_h2,
+
+    prop_pleio = gameInit_traits$value()$prop_pleio,
+    cor_pleio = gameInit_traits$value()$cor_pleio
   )
 
+  print(params)
+
   out_report <- rmarkdown::render("./src/plantbreedgame_setup.Rmd",
     output_file = tempfile(),
     encoding = "UTF-8",