diff --git a/DESCRIPTION b/DESCRIPTION index 089d05a..fc2150d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: FielDHub Title: A Shiny App for Design of Experiments in Life Sciences -Version: 1.4.0 +Version: 1.4.1 Authors@R: c(person(given = "Didier", family = "Murillo", diff --git a/NAMESPACE b/NAMESPACE index 4584522..2602546 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -12,12 +12,14 @@ export(alpha_lattice) export(diagonal_arrangement) export(do_optim) export(full_factorial) +export(improve_efficiency) export(incomplete_blocks) export(latin_square) export(multi_location_prep) export(optimized_arrangement) export(partially_replicated) export(rectangular_lattice) +export(report_efficiency) export(row_column) export(run_app) export(sparse_allocation) @@ -27,6 +29,7 @@ export(split_split_plot) export(square_lattice) export(strip_plot) export(swap_pairs) +export(swap_treatments) import(shiny) importFrom(config,get) importFrom(dplyr,glimpse) diff --git a/R/app_ui.R b/R/app_ui.R index 9db2040..89aaf3e 100644 --- a/R/app_ui.R +++ b/R/app_ui.R @@ -17,7 +17,7 @@ app_ui <- function(request) { tagList( golem_add_external_resources(), fluidPage(theme = shinythemes::shinytheme("flatly"), - navbarPage(title = "FielDHub v1.4.0", + navbarPage(title = "FielDHub v1.4.1", tabPanel( " Welcome!", icon = icon("home", lib = "glyphicon"), suppressWarnings( diff --git a/R/fct_incomplete_blocks.R b/R/fct_incomplete_blocks.R index 6530366..f60d43f 100644 --- a/R/fct_incomplete_blocks.R +++ b/R/fct_incomplete_blocks.R @@ -57,8 +57,8 @@ #' head(ibd2$fieldBook) #' #' @export -incomplete_blocks <- function(t = NULL, k = NULL, r = NULL, l = 1, plotNumber = 101, locationNames = NULL, - seed = NULL, data = NULL) { +incomplete_blocks <- function(t = NULL, k = NULL, r = NULL, l = 1, plotNumber = 101, + locationNames = NULL, seed = NULL, data = NULL) { if (is.null(seed) || !is.numeric(seed)) seed <- runif(1, min = -50000, max = 50000) set.seed(seed) @@ -137,17 +137,12 @@ incomplete_blocks <- function(t = NULL, k = NULL, r = NULL, l = 1, plotNumber = blocks_model <- list() for (i in 1:l) { mydes <- blocksdesign::blocks(treatments = nt, replicates = r, blocks = list(r, b), seed = NULL) - # print("---Blocks Model original design:---") - # print(mydes$Blocks_model) mydes <- rerandomize_ibd(ibd_design = mydes) - # print("---Blocks Model re-randomized design:---") - # print(mydes$Blocks_model_new) matdf <- base::data.frame(list(LOCATION = rep(locationNames[i], each = N))) matdf$PLOT <- as.numeric(unlist(ibd_plots[[i]])) matdf$BLOCK <- rep(c(1:r), each = nt) matdf$iBLOCK <- rep(c(1:b), each = k) matdf$UNIT <- rep(c(1:k), nincblock) - # matdf$TREATMENT <- mydes$Design[,4] matdf$TREATMENT <- mydes$Design_new[,4] colnames(matdf) <- c("LOCATION","PLOT", "REP", "IBLOCK", "UNIT", "ENTRY") outIBD_loc[[i]] <- matdf diff --git a/R/fct_row_column.R b/R/fct_row_column.R index 2c62ede..d1e1477 100644 --- a/R/fct_row_column.R +++ b/R/fct_row_column.R @@ -12,6 +12,7 @@ #' @param plotNumber Numeric vector with the starting plot number for each location. By default \code{plotNumber = 101}. #' @param seed (optional) Real number that specifies the starting seed to obtain reproducible designs. #' @param locationNames (optional) Names for each location. +#' @param iterations Number of iterations for design optimization. By default \code{iterations = 1000}. #' @param data (optional) Data frame with label list of treatments #' #' @author Didier Murillo [aut], @@ -67,14 +68,15 @@ #' #' #' @export -row_column <- function(t = NULL, nrows = NULL, r = NULL, l = 1, plotNumber= 101, locationNames = NULL, - seed = NULL, data = NULL) { +row_column <- function(t = NULL, nrows = NULL, r = NULL, l = 1, plotNumber= 101, + locationNames = NULL, seed = NULL, iterations = 1000, + data = NULL) { if (is.null(seed) || !is.numeric(seed)) seed <- runif(1, min = -50000, max = 50000) - set.seed(seed) + # set.seed(seed) k <- nrows lookup <- FALSE - if(is.null(data)) { + if (is.null(data)) { if (is.null(t) || is.null(k) || is.null(r) || is.null(l)) { shiny::validate('Some of the basic design parameters are missing (t, k, r or l).') } @@ -92,18 +94,20 @@ row_column <- function(t = NULL, nrows = NULL, r = NULL, l = 1, plotNumber= 101, nt <- length(t) TRT <- t } - }else if (is.character(t) || is.factor(t)) { + } else if (is.character(t) || is.factor(t)) { if (length(t) == 1) { shiny::validate('incomplete_blocks() requires more than one treatment.') } nt <- length(t) - }else if ((length(t) > 1)) { + } else if ((length(t) > 1)) { nt <- length(t) } - df <- data.frame(list(ENTRY = 1:nt, TREATMENT = paste0("G-", 1:nt))) - data_RowCol <- df + data_up <- data.frame(list(ENTRY = 1:nt, TREATMENT = paste0("G-", 1:nt))) + colnames(data_up) <- c("ENTRY", "TREATMENT") lookup <- TRUE - }else if (!is.null(data)) { + df <- data.frame(list(ENTRY = 1:nt, LABEL_TREATMENT = paste0("G-", 1:nt))) + dataLookUp <- df + } else if (!is.null(data)) { if (is.null(t) || is.null(r) || is.null(k) || is.null(l)) { shiny::validate('Some of the basic design parameters are missing (t, r, k or l).') } @@ -116,43 +120,95 @@ row_column <- function(t = NULL, nrows = NULL, r = NULL, l = 1, plotNumber= 101, if (t != new_t) base::stop("Number of treatments do not match with data input.") TRT <- data_up$TREATMENT nt <- length(TRT) - data_RowCol <- data_up lookup <- TRUE + dataLookUp <- data.frame(list(ENTRY = 1:nt, LABEL_TREATMENT = TRT)) } if (k >= nt) shiny::validate('incomplete_blocks() requires k < t.') + if (nt %% k != 0) { + shiny::validate('Number of treatments can not be fully distributed over the specified incomplete block specification.') + } if(is.null(locationNames) || length(locationNames) != l) locationNames <- 1:l nunits <- k - matdf <- incomplete_blocks(t = nt, k = nunits, r = r, l = l, plotNumber = plotNumber, - seed = seed, locationNames = locationNames, - data = data_RowCol) - matdf <- matdf$fieldBook - matdf <- as.data.frame(matdf) - colnames(matdf)[5] <- "COLUMN" - matdf$ROW <- matdf$UNIT - OutRowCol <- matdf[,-6] - OutRowCol$LOCATION <- factor(OutRowCol$LOCATION, levels = locationNames) - OutRowCol <- OutRowCol[order(OutRowCol$LOCATION, OutRowCol$REP, OutRowCol$ROW),] - RowCol_plots <- ibd_plot_numbers(nt = nt, plot.number = plotNumber, r = r, l = l) - OutRowCol$PLOT <- as.vector(unlist(RowCol_plots)) + + + ## New code + N <- nt * r + out_row_col_loc <- vector(mode = "list", length = l) + blocks_model <- list() + for (i in 1:l) { + reps <- r + ncols <- nt / nunits + mydes <- blocksdesign::blocks( + treatments = nt, + replicates = reps, + blocks = list(reps, ncols), + seed = seed + i + ) + mydes <- rerandomize_ibd(ibd_design = mydes) + # Create row and column design + row_col_design <- mydes$Design_new |> + dplyr::mutate(Level_3 = rep(rep(paste0("B", 1:nrows), times = ncols), times = reps)) |> + dplyr::mutate(Level_3 = paste(Level_1, Level_3, sep = ".")) |> + dplyr::mutate(Level_3 = factor(Level_3, levels = unique(Level_3))) |> + dplyr::select(Level_1, Level_2, Level_3, plots, treatments) + + improved_design <- improve_efficiency(row_col_design, iterations, seed = seed + i) + field_book_best_design <- improved_design$best_design + row_column_efficiency <- report_efficiency(improved_design$best_design) + blocks_model[[i]] <- row_column_efficiency + row_col_fieldbook <- field_book_best_design |> + dplyr::rename( + REP = Level_1, + COLUMN = Level_2, + ROW = Level_3, + PLOT = plots, + ENTRY = treatments) |> + dplyr::mutate( + REP = as.numeric(factor(REP, levels = unique(REP))), + COLUMN = as.numeric(factor(COLUMN, levels = unique(COLUMN))), + ROW = as.numeric(factor(ROW, levels = unique(ROW))) + ) |> + dplyr::select(PLOT, REP, COLUMN, ROW, ENTRY) + + locations_df <- data.frame(list(LOCATION = rep(locationNames[i], each = N))) + row_col_fieldbook <- dplyr::bind_cols(locations_df, row_col_fieldbook) + out_row_col_loc[[i]] <- row_col_fieldbook + + } + + out_row_col <- dplyr::bind_rows(out_row_col_loc) + out_row_col$ENTRY <- as.numeric(out_row_col$ENTRY) + if(lookup) { - OutRowCol <- OutRowCol[,c(2,3,4,8,5,6,7)] - }else OutRowCol <- OutRowCol[,c(2,3,4,7,5,6)] - ID <- 1:nrow(OutRowCol) - OutRowCol <- cbind(ID, OutRowCol) - rownames(OutRowCol) <- 1:nrow(OutRowCol) - loc <- levels(OutRowCol$LOCATION) + out_row_col <- dplyr::inner_join(out_row_col, dataLookUp, by = "ENTRY") + out_row_col <- out_row_col |> + dplyr::rename(TREATMENT = LABEL_TREATMENT) |> + dplyr::select(-ENTRY) + out_row_col <- dplyr::inner_join(out_row_col, data_up, by = "TREATMENT") |> + dplyr::select(LOCATION, PLOT, REP, ROW, COLUMN, ENTRY, TREATMENT) + } + + ID <- 1:nrow(out_row_col) + out_row_col_id <- cbind(ID, out_row_col) + + row_col_plots <- ibd_plot_numbers(nt = nt, plot.number = plotNumber, r = r, l = l) + out_row_col_id$PLOT <- as.vector(unlist(row_col_plots)) + + # return(list(fieldBook = out_row_col_id, blocks_model = blocks_model)) + + loc <- levels(out_row_col_id$LOCATION) ib <- nt/k Resolvable_rc_reps <- vector(mode = "list", length = r*l) w <- 1 for (sites in 1:l) { for (j in 1:r) { - z <- OutRowCol + z <- out_row_col_id z <- subset(z, z$LOCATION == loc[sites] & z$REP == j) if (is.null(data)){ - Resolvable_rc_reps[[w]] <- matrix(data = as.vector(z$ENTRY), nrow = nunits, + Resolvable_rc_reps[[w]] <- matrix(data = as.vector(z$ENTRY), nrow = nunits, ncol = ib, byrow = TRUE) }else { - Resolvable_rc_reps[[w]] <- matrix(data = as.vector(z$TREATMENT), nrow = nunits, + Resolvable_rc_reps[[w]] <- matrix(data = as.vector(z$TREATMENT), nrow = nunits, ncol = ib, byrow = TRUE) } w <- w + 1 @@ -164,31 +220,35 @@ row_column <- function(t = NULL, nrows = NULL, r = NULL, l = 1, plotNumber= 101, y <- seq(r, r * l, r) z <- 1 for (loc in 1:l) { - NEW_Resolvable[[loc]] <- setNames(Resolvable_rc_reps[x[z]:y[z]], + NEW_Resolvable[[loc]] <- setNames(Resolvable_rc_reps[x[z]:y[z]], paste0(rep("rep", r), 1:r)) z <- z + 1 } - - df <- OutRowCol - trt <- "ENTRY" + + df <- out_row_col_id + trt <- "ENTRY" c1 <- concurrence_matrix(df=df, trt=trt, target='REP') c2 <- concurrence_matrix (df=df, trt=trt, target='ROW') c3 <- concurrence_matrix (df=df, trt=trt, target='COLUMN') summ <- merge(c1, c2, by="Concurrence", all=TRUE) new_summ <- merge(summ, c3, by='Concurrence', all=TRUE) infoDesign <- list( - rows = nrows, - columns = ib, - reps = r, - treatments = nt, - locations = l, - location_names = locationNames, + rows = nrows, + columns = ib, + reps = r, + treatments = nt, + locations = l, + location_names = locationNames, seed = seed, id_design = 9 ) - output <- list(infoDesign = infoDesign, resolvableBlocks = NEW_Resolvable, - concurrence = new_summ, - fieldBook = OutRowCol) + output <- list( + infoDesign = infoDesign, + blocksModel = blocks_model, + resolvableBlocks = NEW_Resolvable, + concurrence = new_summ, + fieldBook = out_row_col_id + ) class(output) <- "FielDHub" return(invisible(output)) } \ No newline at end of file diff --git a/R/mod_RowCol.R b/R/mod_RowCol.R index 5f9e126..d82bcc0 100644 --- a/R/mod_RowCol.R +++ b/R/mod_RowCol.R @@ -46,7 +46,7 @@ mod_RowCol_ui <- function(id){ ns = ns, numericInput(ns("t.rcd"), label = "Input # of Treatments:", - value = 36, + value = 24, min = 2), ), fluidRow( @@ -110,6 +110,16 @@ mod_RowCol_ui <- function(id){ width = 8, fluidRow( tabsetPanel( + tabPanel( + "Summary Design", + br(), + shinycssloaders::withSpinner( + verbatimTextOutput(outputId = ns("summary_row_column"), + placeholder = FALSE), + type = 4 + ), + style = "padding-right: 40px;" + ), tabPanel("Field Layout", shinyjs::useShinyjs(), shinyjs::hidden(downloadButton(ns("downloadCsv.rcd"), @@ -333,6 +343,12 @@ mod_RowCol_server <- function(id){ }) |> bindEvent(input$RUN.rcd) + output$summary_row_column <- renderPrint({ + req(RowCol_reactive()) + cat("Randomization was successful!", "\n", "\n") + print(RowCol_reactive(), n = 6) + }) + upDateSites <- reactive({ req(input$l.rcd) locs <- as.numeric(input$l.rcd) diff --git a/R/utils_S3_methods.R b/R/utils_S3_methods.R index 25affa9..b1b4224 100644 --- a/R/utils_S3_methods.R +++ b/R/utils_S3_methods.R @@ -46,7 +46,7 @@ print.FielDHub <- function(x, n=10, ...){ "First observations of the data frame with the CRD field book:", "\n") print(head(x$fieldBook, n=nhead_print, ...)) - }else if (x$infoDesign$id_design == 2){ + } else if (x$infoDesign$id_design == 2){ cat("Randomized Complete Block Design (RCBD):", "\n\n") #----------------------------------------------------------- cat("Information on the design parameters:", "\n") @@ -62,7 +62,7 @@ print.FielDHub <- function(x, n=10, ...){ "First observations of the data frame with the RCBD field book:", "\n") print(head(x$fieldBook, n=nhead_print, ...)) - }else if (x$infoDesign$id_design == 3){ + } else if (x$infoDesign$id_design == 3){ cat("Latin Square Design:", "\n\n") cat("Information on the design parameters:", "\n") len <- length(x$infoDesign) @@ -92,7 +92,7 @@ print.FielDHub <- function(x, n=10, ...){ "First observations of the data frame with the full_factorial field book:", "\n") print(head(x$fieldBook, n=nhead_print, ...)) - }else if (x$infoDesign$id_design == 5) { + } else if (x$infoDesign$id_design == 5) { cat("Split Plot Design", "\n\n") cat("Information on the design parameters:", "\n") len <- length(x$infoDesign) @@ -107,7 +107,7 @@ print.FielDHub <- function(x, n=10, ...){ "First observations of the data frame with the split_plot field book:", "\n") print(head(x$fieldBook, n=nhead_print, ...)) - }else if (x$infoDesign$id_design == 6) { + } else if (x$infoDesign$id_design == 6) { cat("Split-Split Plot Design", "\n\n") cat("Information on the design parameters:", "\n") len <- length(x$infoDesign) @@ -122,7 +122,7 @@ print.FielDHub <- function(x, n=10, ...){ "First observations of the data frame with the split_split_plot field book:", "\n") print(head(x$fieldBook, n=nhead_print, ...)) - }else if (x$infoDesign$id_design == 7) { + } else if (x$infoDesign$id_design == 7) { cat("Strip Plot Design", "\n\n") cat("Information on the design parameters:", "\n") len <- length(x$infoDesign) @@ -137,7 +137,7 @@ print.FielDHub <- function(x, n=10, ...){ "First observations of the data frame with the strip_plot field book:", "\n") print(head(x$fieldBook, n=nhead_print, ...)) - }else if (x$infoDesign$id_design == 8) { + } else if (x$infoDesign$id_design == 8) { cat("Incomplete Blocks Design", "\n\n") #--------------------------------------------------------------------- cat("Efficiency of design:", "\n") @@ -156,8 +156,12 @@ print.FielDHub <- function(x, n=10, ...){ "First observations of the data frame with the incomplete_blocks field book:", "\n") print(head(x$fieldBook, n=nhead_print, ...)) - }else if (x$infoDesign$id_design == 9) { + } else if (x$infoDesign$id_design == 9) { cat("Row Column Design", "\n\n") + #--------------------------------------------------------------------- + cat("First location efficiency of design:", "\n") + print(x$blocksModel[[1]]) + cat("\n") cat("Information on the design parameters:", "\n") len <- length(x$infoDesign) str(x$infoDesign[1:(len-1)]) @@ -171,7 +175,7 @@ print.FielDHub <- function(x, n=10, ...){ "First observations of the data frame with the row_column field book:", "\n") print(head(x$fieldBook, n=nhead_print, ...)) - }else if (x$infoDesign$id_design == 10) { + } else if (x$infoDesign$id_design == 10) { cat("Square Lattice Design", "\n\n") #--------------------------------------------------------------------- cat("Efficiency of design:", "\n") @@ -190,7 +194,7 @@ print.FielDHub <- function(x, n=10, ...){ "First observations of the data frame with the square_lattice field book:", "\n") print(head(x$fieldBook, n=nhead_print, ...)) - }else if (x$infoDesign$id_design == 11) { + } else if (x$infoDesign$id_design == 11) { cat("Rectangular Lattice Design", "\n\n") #--------------------------------------------------------------------- cat("Efficiency of design:", "\n") @@ -277,7 +281,7 @@ print.FielDHub <- function(x, n=10, ...){ "First observations of the data frame with the RCBD_augmented field book:", "\n") print(head(x$fieldBook, n=nhead_print, ...)) - }else if (x$infoDesign$id_design == 15) { + } else if (x$infoDesign$id_design == 15) { cat("Un-replicated Diagonal Arrangement Design", "\n\n") cat("Information on the design parameters:", "\n") len <- length(x$infoDesign) @@ -322,7 +326,7 @@ print.FielDHub <- function(x, n=10, ...){ "First observations of the data frame with the optimized_arrangement field book:", "\n") print(head(x$fieldBook, n=nhead_print, ...)) - }else if (x$infoDesign$id_design == 17) { + } else if (x$infoDesign$id_design == 17) { cat("Split families:", "\n\n") #--------------------------------------------------------------------- # Head diff --git a/R/utils_row_col_optimization.R b/R/utils_row_col_optimization.R new file mode 100644 index 0000000..2bdbf97 --- /dev/null +++ b/R/utils_row_col_optimization.R @@ -0,0 +1,106 @@ +# Function to randomly swap a pair of treatments within a random +# level of Level_2 for all levels of Level_1 +#' @export +swap_treatments <- function(df) { + # Split the dataframe by Level_1 + df_split <- split(df, df$Level_1) + + # Initialize an empty dataframe to store the results + result <- data.frame() + + # Loop through each level of Level_1 + for (level1 in names(df_split)) { + df_level1 <- df_split[[level1]] + + # Get the unique values of Level_2 + unique_levels2 <- unique(df_level1$Level_2) + + # Randomly select one level of Level_2 + random_level2 <- sample(unique_levels2, 1) + df_level2 <- df_level1[df_level1$Level_2 == random_level2, ] + + if (nrow(df_level2) >= 2) { + # Randomly select two different rows to swap treatments + rows_to_swap <- sample(1:nrow(df_level2), 2) + temp <- df_level2$treatments[rows_to_swap[1]] + df_level2$treatments[rows_to_swap[1]] <- df_level2$treatments[rows_to_swap[2]] + df_level2$treatments[rows_to_swap[2]] <- temp + } + + # Combine the modified dataframe with the rest + df_level1[df_level1$Level_2 == random_level2, ] <- df_level2 + result <- rbind(result, df_level1) + } + + return(result) +} + +# Function to improve A-Efficiency for Level 2 +#' @export +improve_efficiency <- function(design, iterations, seed) { + set.seed(seed) + # Initial design + best_design <- design + + # Calculate initial efficiencies + row_blocks <- best_design |> + dplyr::select(Level_1, Level_3, plots, treatments) + efficiencies <- BlockEfficiencies(row_blocks) + best_a_efficiency <- efficiencies$`A-Efficiency`[efficiencies$Level == 2] + + # Run iterations to improve A-Efficiency + for (i in 1:iterations) { + # Generate a new design by swapping treatments + new_design <- swap_treatments(best_design) + + # Calculate efficiencies for the new design + new_row_blocks <- new_design |> + dplyr::select(Level_1, Level_3, plots, treatments) + new_efficiencies <- BlockEfficiencies(new_row_blocks) + new_a_efficiency <- new_efficiencies$`A-Efficiency`[new_efficiencies$Level == 2] + + # Update the best design if the new A-Efficiency is higher + if (new_a_efficiency > best_a_efficiency) { + best_design <- new_design + best_a_efficiency <- new_a_efficiency + best_efficiencies <- new_efficiencies + } + } + return( + list( + best_design = best_design, + best_efficiencies = best_efficiencies, + best_a_efficiency = best_a_efficiency + ) + ) +} + +# Function to calculate and return combined BlockEfficiencies +#' @export +report_efficiency <- function(design) { + # Calculate row block efficiencies + row_blocks <- design |> + dplyr::select(Level_1, Level_3, plots, treatments) + row_efficiencies <- BlockEfficiencies(row_blocks) + row_efficiencies <- row_efficiencies |> + dplyr::filter(Level == 2) |> + dplyr::mutate(Level = "Row") + + # Calculate column block efficiencies + col_blocks <- design |> + dplyr::select(Level_1, Level_2, plots, treatments) + col_efficiencies <- BlockEfficiencies(col_blocks) + col_efficiencies <- col_efficiencies |> + dplyr::filter(Level == 2) |> + dplyr::mutate(Level = "Column") + + # Get replication efficiencies + rep_efficiencies <- BlockEfficiencies(row_blocks) |> + dplyr::filter(Level == 1) |> + dplyr::mutate(Level = "Rep") + + # Combine the results + combined_efficiencies <- dplyr::bind_rows(rep_efficiencies, row_efficiencies, col_efficiencies) + + return(combined_efficiencies) +} diff --git a/man/row_column.Rd b/man/row_column.Rd index c3f540b..66fc530 100644 --- a/man/row_column.Rd +++ b/man/row_column.Rd @@ -12,6 +12,7 @@ row_column( plotNumber = 101, locationNames = NULL, seed = NULL, + iterations = 1000, data = NULL ) } @@ -30,6 +31,8 @@ row_column( \item{seed}{(optional) Real number that specifies the starting seed to obtain reproducible designs.} +\item{iterations}{Number of iterations for design optimization. By default \code{iterations = 1000}.} + \item{data}{(optional) Data frame with label list of treatments} } \value{