diff --git a/.Rbuildignore b/.Rbuildignore index d0aba75..7e1a4c2 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -21,3 +21,4 @@ codemeta.json ^Meta$ ^CODE_OF_CONDUCT\.md$ ^doc$ +^\.dir-locals\.el$ diff --git a/DESCRIPTION b/DESCRIPTION index df08a3a..4038760 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Type: Package Package: speed Title: Generate Spatially Efficient Experimental Designs Version: 0.0.6 -Authors@R: +Authors@R: c(person(given = "Sam", family = "Rogers", role = c("aut", "cre"), @@ -19,31 +19,39 @@ Authors@R: family = "Pipattungsakul", role = "aut", email = "wasin.pipattungsakul@adelaide.edu.au"), + person(given = "Kai", + family = "Bagley", + email = "kai.bagley@curtin.edu.au", + role = "aut", + comment = c(ORCID = "0009-0004-6579-6959")), person(given = "University of Adelaide", role = c("cph", "fnd"), comment = "https://adelaide.edu.au/"), person(given = "Grains Research and Development Corporation", role = c("cph", "fnd"), - comment = "https://grdc.com.au/")) -Description: The Speed package optimises spatial experimental designs by rearranging - treatments to improve statistical efficiency while maintaining statistical validity. - It employs customisable optimisation metrics, with a default approach that balances - minimising treatment adjacency and maintaining spatial balance. Users can supply their - own optimisation metrics to tailor designs to specific experimental requirements. - The package also includes visualisation tools that allow users to inspect + comment = "https://grdc.com.au/"), + person(given = "Curtin University", + role = c("cph"), + comment = c("http://www.curtin.edu.au/", ROR = "02n415q13"))) +Description: The Speed package optimises spatial experimental designs by rearranging + treatments to improve statistical efficiency while maintaining statistical validity. + It employs customisable optimisation metrics, with a default approach that balances + minimising treatment adjacency and maintaining spatial balance. Users can supply their + own optimisation metrics to tailor designs to specific experimental requirements. + The package also includes visualisation tools that allow users to inspect optimised designs. License: MIT + file LICENSE URL: https://biometryhub.github.io/speed/ BugReports: https://github.com/biometryhub/speed/issues Depends: R (>= 4.1) -Imports: +Imports: farver, ggplot2, matrixStats, rlang, scales, stringi -Suggests: +Suggests: igraph, patchwork, quarto, diff --git a/NAMESPACE b/NAMESPACE index 889d8fa..b1c9732 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,11 +4,18 @@ S3method(autoplot,design) S3method(print,design) export(add_buffers) export(autoplot) +export(calc_concurrence_matrix) +export(calc_incidence_matrix) +export(calc_info_matrix) export(calculate_adjacency_score) export(calculate_balance_score) export(calculate_ed) export(calculate_efficiency_factor) +export(calculate_efficiency_factors) export(calculate_nb) +export(compute_L_projection) +export(cor_ar1) +export(cor_ar1_ar1) export(create_pair_mapping) export(get_edges) export(get_vertices) @@ -16,6 +23,7 @@ export(initialise_design_df) export(initialize_design_df) export(objective_function) export(objective_function_factorial) +export(objective_function_info) export(objective_function_piepho) export(objective_function_signature) export(optim_params) diff --git a/R/objectives.R b/R/objectives.R new file mode 100644 index 0000000..9fd72cb --- /dev/null +++ b/R/objectives.R @@ -0,0 +1,434 @@ +#' Objective function using the treatment information matrix +#' +#' Creates an objective function that optimises experimental designs using the +#' Fisher information for treatment contrasts after adjusting for nuisance +#' effects, and takes an optional spatial covariance structure for spatial +#' optimisation. +#' +#' @param layout_df Data frame representing the spatial layout of the +#' experiment. +#' @param swap Column name to swap, usually the treatment. +#' @param spatial_cols Column name of the spatial factors. +#' @param criterion Either \code{"A"} or \code{"D"}, representing A or D +#' optimality. +#' - A-optimality: Minimises \eqn{\mathrm{tr} \left( \mathcal{I}^- +#' \right)}{tr(I⁻)}. +#' - D-optimality: Minimises \eqn{-\log \left| \mathcal{I} \right|}{- +#' log(|I|)} +#' @param L_matrix Precomputed projection matrix. Use +#' \code{precompute_projection} to generate it. If \code{NULL}, then the +#' identity covariance structure will be assumed, and the projection will be +#' computed using the structure of the design's blocks. +#' @param block_column Column name of the design's block factor. Used when +#' \code{L_matrix} isn't supplied. +#' @param ... Extra parameters passed from \code{speed}. +#' +#' @details +#' This function computes the treatment information matrix: +#' \deqn{I = X_1^\intercal L X_1}{I = X₁ᵀ L X₁} +#' Where \eqn{X_1}{X₁} is the treatment design matrix, and \eqn{L} is the +#' projection that removes nuisance fixed effects from the GLS's inverse +#' covariance: +#' \deqn{L = \Sigma^{-1} - \Sigma^{-1} X_2 \left( X_2^\intercal \Sigma^{-1} X_2 +#' \right)^{-1}}{L = Σ⁻¹ - Σ⁻¹ X₂ (X₂ᵀ Σ⁻¹ X₂)⁻¹} +#' The user specifies the spatial covariance \eqn{\Sigma}{Σ} they intend to use +#' at the time of analysis, e.g. lag-1 autoregressive \eqn{\text{AR}_1 \otimes +#' \text{AR}_1}{AR₁ ⊗ AR₁}, and this objective will optimise under that +#' structure. +#' +#' @examples +#' # Small RCBD layout: 6 treatments in 4 blocks of 6 plots +#' df <- initialise_design_df( +#' items = 6, nrows = 4, ncols = 6, +#' block_nrows = 1, block_ncols = 6 +#' ) +#' +#' # Non-spatial: identity covariance +#' result <- speed( +#' df, +#' swap = "treatment", +#' swap_within = "block", +#' obj_function = objective_function_info, +#' criterion = "A", +#' seed = 42, +#' quiet = TRUE +#' ) +#' +#' # Spatial: AR(1) x AR(1) +#' Sigma <- cor_ar1_ar1( +#' n_rows = 4, n_cols = 6, +#' rho_row = 0.6, rho_col = 0.3 +#' ) +#' L <- compute_L_projection(df, Sigma, block_column = "block") +#' result_spatial <- speed( +#' df, +#' swap = "treatment", +#' swap_within = "block", +#' obj_function = objective_function_info, +#' L_matrix = L, +#' criterion = "A", +#' optimise_params = optim_params(random_initialisation = TRUE), +#' seed = 42, +#' quiet = TRUE +#' ) +#' +#' # Inspect the optimised design +#' calc_info_matrix(result_spatial$design_df, L_matrix = L) +#' +#' @return A named list with \code{score}, \code{info_matrix}, +#' \code{eigenvalues}, \code{criterion}. +#' +#' @export +objective_function_info <- function( + layout_df, + swap, + spatial_cols, + criterion = c("A", "D"), + L_matrix = NULL, + block_column = "block", + ... +) { + criterion <- match.arg(criterion) + + ci <- .compute_info(layout_df, swap, L_matrix, block_column) + info <- ci$info + v <- ci$v + + ## get positive eigenvals from info matrix + eig <- eigen(info, symmetric = TRUE, only.values = TRUE)$values + max_eig <- max(eig) + if (max_eig <= 0) { + score <- 1e10 + pos_eig <- numeric(0) + } else { + pos_eig <- eig[eig > max_eig * 1e-10] + ## if design is disconnected, penalise greatly + ## otherwise A or D + if (length(pos_eig) < v - 1) { + score <- 1e10 + } else { + score <- switch(criterion, + A = sum(1 / pos_eig), + D = -sum(log(pos_eig)) + ) + if (!is.finite(score)) score <- 1e10 + } + } + + return(list( + score = score, + info_matrix = info, + eigenvalues = sort(pos_eig, decreasing = TRUE), + criterion = criterion + )) +} + +#' Build the treatment indicator matrix X1 +#' @noRd +.build_treatment_matrix <- function(treatments, trt_levels, n, v) { + X1 <- matrix(0, n, v) + trt_idx <- match(treatments, trt_levels) + for (i in seq_len(n)) { + X1[i, trt_idx[i]] <- 1 + } + return(X1) +} + +#' Function to precompute \eqn{L} +#' +#' \eqn{L} is the projection matrix that removes nuisance effects under a +#' spatial correlation structure, given an \eqn{n \times n}{n * n} covariance +#' \eqn{\Sigma}{Sigma} and nuisance effects design \eqn{X_2}. +#' +#' This only depends on lambda and the block structure, so precompute this. +#' @noRd +.build_L_from_df <- function(layout_df, block_column, n) { + if (!block_column %in% names(layout_df)) { + ## L = I - (1/n) ones_matrix + return(diag(n) - matrix(1 / n, n, n)) + } + + blocks <- factor(layout_df[[block_column]]) + b <- nlevels(blocks) + + X2 <- matrix(0, n, b) + for (j in seq_len(b)) { + X2[blocks == levels(blocks)[j], j] <- 1 + } + + ## L = I - X2 (X2t X2)^{-1} X2t + X2tX2_inv <- diag(1 / colSums(X2), nrow = b) + diag(n) - X2 %*% X2tX2_inv %*% t(X2) +} + +#' Function for computing the information matrix +#' @noRd +.compute_info <- function( + layout_df, + swap, + L_matrix, + block_column +) { + ## a lot of the args come from the design now + treatments <- layout_df[[swap]] + trt_levels <- levels(factor(treatments)) + v <- length(trt_levels) + n <- nrow(layout_df) + + ## L + ### projection matrix + if (is.null(L_matrix)) { + L_matrix <- .build_L_from_df(layout_df, block_column, n) + } + + ## X1 + ### treatment design matrix. + X1 <- .build_treatment_matrix(treatments, trt_levels, n, v) + ## information matrix I = X1t L X1 + info <- t(X1) %*% L_matrix %*% X1 + + list(info = info, v = v, trt_levels = trt_levels) +} + +# Helper functions + +#' Compute L projection +#' +#' This is the projection matrix that removes nuisance effects from the GLS's +#' covariance +#' +#' @param layout_df Data frame representing the spatial layout of the +#' experiment. +#' @param Sigma Covariance structure to use. +#' @param block_column Column name of the design's block factor in +#' \code{layout_df}. +#' +#' @return \eqn{(n \times n)} numeric matrix. +#' +#' @export +compute_L_projection <- function( + layout_df, + Sigma, + block_column = "block" +) { + n <- nrow(layout_df) + + ## arg checks + if (!is.matrix(Sigma)) stop("Sigma must be a matrix") + if (nrow(Sigma) != n) { + stop( + "Sigma matrix must have dimension equal ", + "to dimension of layout_df dataframe") + } + if (ncol(Sigma) != n) { + stop( + "Sigma matrix must have dimension equal ", + "to dimension of layout_df dataframe" + ) + } + if (!isSymmetric(Sigma, tol = 1e-8)) stop("Sigma must be symmetric") + if (!(block_column %in% names(layout_df))) { + stop("Block column name not in layout_df dataframe") + } + + blocks <- factor(layout_df[[block_column]]) + b <- nlevels(blocks) + + ## X2 + ### block indicator matrix with the intercept removed + X2 <- matrix(0, n, b) + for (j in seq_len(b)) { + X2[blocks == levels(blocks)[j], j] <- 1 + } + + Sigma_inv <- solve(Sigma) + SiX2 <- Sigma_inv %*% X2 + return(Sigma_inv - SiX2 %*% solve(t(X2) %*% SiX2) %*% t(SiX2)) +} + +#' Calculate info matrix +#' +#' Calculate the Fisher information in the experimental design +#' +#' @param layout_df Data frame representing the spatial layout of the +#' experiment. +#' @param treatment_column Column name containing the design's treatments in +#' \code{layout_df}. +#' @param L_matrix Precomputed projection matrix from +#' \code{compute_L_projection}. +#' @param block_column Column name containing the designs block in +#' \code{layout_df}. +#' +#' @export +calc_info_matrix <- function( + layout_df, + treatment_column = "treatment", + L_matrix = NULL, + block_column = "block" +) { + ci <- .compute_info(layout_df, treatment_column, L_matrix, block_column) + info <- ci$info + v <- ci$v + rownames(info) <- ci$trt_levels + colnames(info) <- ci$trt_levels + + ## get positive eigenvals from info matrix + eig <- eigen(info, symmetric = TRUE, only.values = TRUE)$values + max_eig <- max(eig) + if (max_eig <= 0) { + pos_eig <- numeric(0) + } else { + pos_eig <- eig[eig > max_eig * 1e-10] + } + + if (length(pos_eig) < v - 1) { + A_val <- Inf + D_val <- Inf + } else { + A_val = sum(1 / pos_eig) + D_val = -sum(log(pos_eig)) + if (!is.finite(A_val)) A_val <- Inf + if (!is.finite(D_val)) D_val <- Inf + } + + return(list( + info_matrix = info, + eigenvalues = sort(pos_eig, decreasing = TRUE), + rank = length(pos_eig), + A_value = A_val, + D_value = D_val + )) +} + +#' Calculate incidence matrix +#' +#' @param layout_df Data frame representing the spatial layout of the +#' experiment. +#' @param treatment_column Column name containing the design's treatments in +#' \code{layout_df}. +#' @param block_column Column name containing the designs block in +#' \code{layout_df}. +#' +#' @return \eqn{v \times b} matrix of integers, with treatment levels as row +#' names, and block levels as column names. +#' +#' @export +calc_incidence_matrix <- function( + layout_df, + treatment_column = "treatment", + block_column = "block" +) { + treatments <- factor(layout_df[[treatment_column]]) + blocks <- factor(layout_df[[block_column]]) + + ## make a matrix of integers + N <- table(treatments, blocks) + class(N) <- "matrix" + storage.mode(N) <- "integer" + + return(N) +} + +#' Calculate concurrence matrix +#' +#' @param layout_df Data frame representing the spatial layout of the +#' experiment. +#' @param treatment_column Column name containing the design's treatments in +#' \code{layout_df}. +#' @param block_column Column name containing the designs block in +#' \code{layout_df}. +#' +#' @return Symmetric \eqn{N N^\intercal}{N Nᵀ} matrix. +#' +#' @export +calc_concurrence_matrix <- function( + layout_df, + treatment_column = "treatment", + block_column = "block" +) { + N <- calc_incidence_matrix(layout_df, treatment_column, block_column) + C <- N %*% t(N) + storage.mode(C) <- "integer" + return(C) +} + +#' Calculate canonical efficiency factors +#' +#' @param layout_df Data frame representing the spatial layout of the +#' experiment. +#' @param treatment_column Column name containing the design's treatments in +#' \code{layout_df}. +#' @param L_matrix Precomputed projection matrix from +#' \code{compute_L_projection}. +#' @param block_column Column name containing the designs block in +#' \code{layout_df}. +#' +#' @return A list with \code{efficiency_factors}, \code{E}, \code{replication}. +#' +#' @export +calculate_efficiency_factors <- function( + layout_df, + treatment_column = "treatment", + L_matrix = NULL, + block_column = "block" +) { + info <- calc_info_matrix( + layout_df, + treatment_column, + L_matrix, + block_column + ) + + treatments <- factor(layout_df[[treatment_column]]) + reps <- table(treatments) + r_bar <- mean(reps) + + eff_factors <- info$eigenvalues / r_bar + + v <- nlevels(treatments) + E <- if (info$rank >= v - 1) { + (v - 1) / sum(1 / eff_factors) + } else { + NA_real_ + } + + list( + efficiency_factors = eff_factors, + E = E, + replication = setNames(as.integer(reps), levels(treatments)) + ) +} + + +# Convenience correlation structure constructors + +#' Construct a 1 dimensional lag-1 autoregressive covariance +#' +#' @param n Size of the covariance vector to generate +#' @param rho Correlation parameter +#' +#' @return A symmetric \eqn{(n \times n)} correlation matrix. +#' +#' @export +cor_ar1 <- function(n, rho) { + rho^abs(outer(seq_len(n), seq_len(n), "-")) +} + + +#' Construct a 2 dimensional lag-1 autoregressive covariance +#' +#' @param n_rows Number of rows in the design. +#' @param n_cols Number of columns in the design. +#' @param rho_row Correlation parameter along the row direction. +#' @param rho_col Correlation parameter along the column direction. +#' +#' @return A symmetric \eqn{(\text{n_rows} \cdot \text{n_cols}) \times +#' (\text{n_rows} \cdot \text{n_cols})} correlation matrix. +#' +#' @export +cor_ar1_ar1 <- function(n_rows, n_cols, rho_row, rho_col) { + kronecker( + cor_ar1(n_rows, rho_row), + cor_ar1(n_cols, rho_col) + ) +} diff --git a/man/calc_concurrence_matrix.Rd b/man/calc_concurrence_matrix.Rd new file mode 100644 index 0000000..463ee5c --- /dev/null +++ b/man/calc_concurrence_matrix.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/objectives.R +\name{calc_concurrence_matrix} +\alias{calc_concurrence_matrix} +\title{Calculate concurrence matrix} +\usage{ +calc_concurrence_matrix( + layout_df, + treatment_column = "treatment", + block_column = "block" +) +} +\arguments{ +\item{layout_df}{Data frame representing the spatial layout of the +experiment.} + +\item{treatment_column}{Column name containing the design's treatments in +\code{layout_df}.} + +\item{block_column}{Column name containing the designs block in +\code{layout_df}.} +} +\value{ +Symmetric \eqn{N N^\intercal}{N Nᵀ} matrix. +} +\description{ +Calculate concurrence matrix +} diff --git a/man/calc_incidence_matrix.Rd b/man/calc_incidence_matrix.Rd new file mode 100644 index 0000000..8f113bb --- /dev/null +++ b/man/calc_incidence_matrix.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/objectives.R +\name{calc_incidence_matrix} +\alias{calc_incidence_matrix} +\title{Calculate incidence matrix} +\usage{ +calc_incidence_matrix( + layout_df, + treatment_column = "treatment", + block_column = "block" +) +} +\arguments{ +\item{layout_df}{Data frame representing the spatial layout of the experiment.} + +\item{treatment_column}{Column name containing the design's treatments in +\code{layout_df}.} + +\item{block_column}{Column name containing the designs block in +\code{layout_df}.} +} +\value{ +\eqn{v \times b} matrix of integers, with treatment levels as row +names, and block levels as column names. +} +\description{ +Calculate incidence matrix +} diff --git a/man/calc_info_matrix.Rd b/man/calc_info_matrix.Rd new file mode 100644 index 0000000..5887f42 --- /dev/null +++ b/man/calc_info_matrix.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/objectives.R +\name{calc_info_matrix} +\alias{calc_info_matrix} +\title{Calculate info matrix} +\usage{ +calc_info_matrix( + layout_df, + treatment_column = "treatment", + L_matrix = NULL, + block_column = "block" +) +} +\arguments{ +\item{layout_df}{Data frame representing the spatial layout of the experiment.} + +\item{treatment_column}{Column name containing the design's treatments in +\code{layout_df}.} + +\item{L_matrix}{Precomputed projection matrix from +\code{compute_L_projection}.} + +\item{block_column}{Column name containing the designs block in +\code{layout_df}.} +} +\description{ +Calculate the Fisher information in the experimental design +} diff --git a/man/calculate_efficiency_factors.Rd b/man/calculate_efficiency_factors.Rd new file mode 100644 index 0000000..cadc96f --- /dev/null +++ b/man/calculate_efficiency_factors.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/objectives.R +\name{calculate_efficiency_factors} +\alias{calculate_efficiency_factors} +\title{Calculate canonical efficiency factors} +\usage{ +calculate_efficiency_factors( + layout_df, + treatment_column = "treatment", + L_matrix = NULL, + block_column = "block" +) +} +\arguments{ +\item{layout_df}{Data frame representing the spatial layout of the experiment.} + +\item{treatment_column}{Column name containing the design's treatments in +\code{layout_df}.} + +\item{L_matrix}{Precomputed projection matrix from +\code{compute_L_projection}.} + +\item{block_column}{Column name containing the designs block in +\code{layout_df}.} +} +\value{ +A list with \code{efficiency_factors}, \code{E}, \code{replication}. +} +\description{ +Calculate canonical efficiency factors +} diff --git a/man/compute_L_projection.Rd b/man/compute_L_projection.Rd new file mode 100644 index 0000000..0d401ea --- /dev/null +++ b/man/compute_L_projection.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/objectives.R +\name{compute_L_projection} +\alias{compute_L_projection} +\title{Compute L projection} +\usage{ +compute_L_projection(layout_df, Sigma, block_column = "block") +} +\arguments{ +\item{layout_df}{Data frame representing the spatial layout of the experiment.} + +\item{Sigma}{Covariance structure to use.} + +\item{block_column}{Column name of the design's block factor in \code{layout_df}.} +} +\value{ +\eqn{(n \times n)} numeric matrix. +} +\description{ +This is the projection matrix that removes nuisance effects from the GLS's +covariance +} diff --git a/man/objective_function_info.Rd b/man/objective_function_info.Rd new file mode 100644 index 0000000..0f0993d --- /dev/null +++ b/man/objective_function_info.Rd @@ -0,0 +1,101 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/objectives.R +\name{objective_function_info} +\alias{objective_function_info} +\title{Objective function using the treatment information matrix} +\usage{ +objective_function_info( + layout_df, + swap, + spatial_cols, + criterion = c("A", "D"), + L_matrix = NULL, + block_column = "block", + ... +) +} +\arguments{ +\item{layout_df}{Data frame representing the spatial layout of the experiment.} + +\item{swap}{Column name to swap, usually the treatment.} + +\item{spatial_cols}{Column name of the spatial factors.} + +\item{criterion}{Either \code{"A"} or \code{"D"}, representing A or D +optimality. +- A-optimality: Minimises \eqn{\mathrm{tr} \left( \mathcal{I}^- \right)}{tr(I⁻)}. + - D-optimality: Minimises \eqn{-\log \left| \mathcal{I} \right|}{-log(|I|)}} + +\item{L_matrix}{Precomputed projection matrix. Use +\code{precompute_projection} to generate it. If \code{NULL}, then the +identity covariance structure will be assumed, and the projection will be +computed using the structure of the design's blocks.} + +\item{block_column}{Column name of the design's block factor. Used when +\code{L_matrix} isn't supplied.} + +\item{...}{Extra parameters passed from \code{speed}.} +} +\value{ +A named list with \code{score}, \code{info_matrix}, +\code{eigenvalues}, \code{criterion}. +} +\description{ +Creates an objective function that optimises experimental designs using the +Fisher information for treatment contrasts after adjusting for nuisance +effects, and takes an optional spatial covariance structure for spatial +optimisation. +} +\details{ +This function computes the treatment information matrix: +\deqn{I = X_1^\intercal L X_1}{I = X₁ᵀ L X₁} +Where \eqn{X_1}{X₁} is the treatment design matrix, and \eqn{L} is the +projection that removes nuisance fixed effects from the GLS's inverse +covariance: +\deqn{L = \Sigma^{-1} - \Sigma^{-1} X_2 \left( X_2^\intercal \Sigma^{-1} X_2 +\right)^{-1}}{L = Σ⁻¹ - Σ⁻¹ X₂ (X₂ᵀ Σ⁻¹ X₂)⁻¹} +The user specifies the spatial covariance \eqn{\Sigma}{Σ} they intend to use +at the time of analysis, e.g. lag-1 autoregressive \eqn{\text{AR}_1 \otimes +\text{AR}_1}{AR₁ ⊗ AR₁}, and this objective will optimise under that +structure. +} +\examples{ +# Small RCBD layout: 6 treatments in 4 blocks of 6 plots +df <- initialise_design_df( + items = 6, nrows = 4, ncols = 6, + block_nrows = 1, block_ncols = 6 +) + +# Non-spatial: identity covariance +result <- speed( + df, + swap = "treatment", + swap_within = "block", + obj_function = objective_function_info, + criterion = "A", + seed = 42, + quiet = TRUE +) + +# Spatial: AR(1) x AR(1) +Sigma <- cor_ar1_ar1( + n_rows = 4, n_cols = 6, + rho_row = 0.6, rho_col = 0.3 +) +L <- compute_L_projection(df, Sigma, block_column = "block") +result_spatial <- speed( + df, + swap = "treatment", + swap_within = "block", + obj_function = objective_function_info, + L_matrix = L, + criterion = "A", + optimise_params = optim_params(random_initialisation = TRUE), + seed = 42, + quiet = TRUE +) + +# Inspect the optimised design +calc_info_matrix(result_spatial$design_df, L_matrix = L) + +} diff --git a/man/speed-package.Rd b/man/speed-package.Rd index 49e43ba..746f653 100644 --- a/man/speed-package.Rd +++ b/man/speed-package.Rd @@ -23,12 +23,14 @@ Authors: \item Julian Taylor \email{julian.taylor@adelaide.edu.au} \item Russell Edson \email{russell.edson@adelaide.edu.au} \item Wasin Pipattungsakul \email{wasin.pipattungsakul@adelaide.edu.au} + \item Kai Bagley \email{kai.bagley@curtin.edu.au} (\href{https://orcid.org/0009-0004-6579-6959}{ORCID}) } Other contributors: \itemize{ \item University of Adelaide (https://adelaide.edu.au/) [copyright holder, funder] \item Grains Research and Development Corporation (https://grdc.com.au/) [copyright holder, funder] + \item Curtin University (\href{https://ror.org/02n415q13}{ROR}) (http://www.curtin.edu.au/) [copyright holder] } } diff --git a/tests/testthat/test-objective_function_info.R b/tests/testthat/test-objective_function_info.R new file mode 100644 index 0000000..9c89379 --- /dev/null +++ b/tests/testthat/test-objective_function_info.R @@ -0,0 +1,378 @@ +# tests/testthat/test-objective_function_info.R + +# ---- objective_function_info + +# all RCBDs are equally optimal +test_that("all valid RCBD randomisations produce the same score under identity covariance", { + set.seed(1) + n_treatments <- 4 + n_blocks <- 3 + + df <- data.frame( + row = rep(1:n_blocks, each = n_treatments), + col = rep(1:n_treatments, times = n_blocks), + block = rep(1:n_blocks, each = n_treatments), + treatment = NA_character_ + ) + + scores <- numeric(20) + for (i in seq_len(20)) { + for (b in 1:n_blocks) { + idx <- df$block == b + df$treatment[idx] <- paste0("T", sample(n_treatments)) + } + result <- objective_function_info( + df, swap = "treatment", spatial_cols = c("row", "col"), + criterion = "A" + ) + scores[i] <- result$score + } + + expect_true(max(scores) - min(scores) < 1e-10) +}) + +# when there IS spatial correlation, two different arrangements of the same +# treatments must give different scores. df1 should be worse, since treatments +# are always in the same column position, but in df2 they're "randomised". +test_that("spatial designs differentiate between arrangements", { + n_treatments <- 4 + n_blocks <- 3 + + df <- data.frame( + row = rep(1:n_blocks, each = n_treatments), + col = rep(1:n_treatments, times = n_blocks), + block = rep(1:n_blocks, each = n_treatments) + ) + + Sigma <- cor_ar1_ar1(n_blocks, n_treatments, rho_row = 0.7, rho_col = 0.5) + L <- compute_L_projection(df, Sigma, block_col = "block") + + # Same treatment ordering in every block + df1 <- df + df1$treatment <- rep(paste0("T", 1:n_treatments), n_blocks) + + # Alternated ordering + df2 <- df + df2$treatment <- c( + paste0("T", c(1, 2, 3, 4)), + paste0("T", c(3, 4, 1, 2)), + paste0("T", c(4, 3, 2, 1)) + ) + + s1 <- objective_function_info( + df1, "treatment", c("row", "col"), + criterion = "A", L_matrix = L + )$score + + s2 <- objective_function_info( + df2, "treatment", c("row", "col"), + criterion = "A", L_matrix = L + )$score + + expect_false(abs(s1 - s2) < 1e-10) +}) + +# just make sure this works lel +test_that("D-optimality criterion works", { + df <- data.frame( + row = rep(1:3, each = 4), + col = rep(1:4, times = 3), + block = rep(1:3, each = 4), + treatment = rep(paste0("T", 1:4), 3) + ) + + result <- objective_function_info( + df, "treatment", c("row", "col"), + criterion = "D" + ) + + expect_true(is.finite(result$score)) + expect_equal(result$criterion, "D") + expect_equal(length(result$eigenvalues), 3) # v - 1 +}) + +# if treatments never exist in the same block together, you cant compare them, +# so the design is fully broken +test_that("disconnected designs have massive penalty", { + df <- data.frame( + row = 1:6, + col = c(1, 2, 1, 2, 1, 2), + block = c(1, 1, 2, 2, 3, 3), + treatment = c("A", "A", "B", "B", "C", "C") + ) + + result <- objective_function_info( + df, "treatment", c("row", "col"), + criterion = "A" + ) + + expect_gte(result$score, 1e9) + expect_true(is.finite(result$score)) +}) + + +# ---- compute_L_projection + +# input validation test +test_that("compute_L_projection validates inputs", { + df <- data.frame( + row = 1:4, col = 1:4, block = c(1, 1, 2, 2), + treatment = c("A", "B", "A", "B") + ) + + # Wrong dimensions + expect_error(compute_L_projection(df, diag(3))) + + # Not symmetric + expect_error(compute_L_projection(df, matrix(1:16, 4, 4))) + + # Missing block column + expect_error(compute_L_projection(df, diag(4), block_col = "nonexistent")) +}) + +# identity sigma should be the default L projection. i.e. spatial optimisation +# with identity correlation is the same as non-spatial +test_that("L with identity Sigma matches .build_L_from_df", { + df <- data.frame( + row = rep(1:3, each = 4), + col = rep(1:4, times = 3), + block = rep(1:3, each = 4), + treatment = rep(paste0("T", 1:4), 3) + ) + + L_spatial <- compute_L_projection(df, diag(12), block_col = "block") + L_direct <- .build_L_from_df(df, "block", 12) + + expect_true(max(abs(L_spatial - L_direct)) < 1e-10) +}) + + +# ---- calc_info_matrix + +# balanced RCBD should have rank = v - 1 +test_that("info matrix for RCBD has correct properties", { + df <- data.frame( + row = rep(1:3, each = 5), + col = rep(1:5, times = 3), + block = rep(1:3, each = 5), + treatment = rep(paste0("T", 1:5), 3) + ) + + info <- calc_info_matrix(df) + + # v = 5, so rank = v - 1 = 4 + expect_equal(info$rank, 4) + + # For balanced RCBD all non-zero eigenvalues = b = 3 + expect_true(all(abs(info$eigenvalues - 3) < 1e-10)) + + # 5x5 matrix + expect_equal(dim(info$info_matrix), c(5, 5)) + + # Named + expect_equal(rownames(info$info_matrix), paste0("T", 1:5)) + expect_equal(colnames(info$info_matrix), paste0("T", 1:5)) +}) + + +# spatial should always be different to non-spatial +test_that("spatial info matrix differs from non-spatial", { + df <- data.frame( + row = rep(1:3, each = 4), + col = rep(1:4, times = 3), + block = rep(1:3, each = 4), + treatment = rep(paste0("T", 1:4), 3) + ) + + Sigma <- cor_ar1_ar1(3, 4, 0.5, 0.3) + L <- compute_L_projection(df, Sigma) + + info_plain <- calc_info_matrix(df) + info_spatial <- calc_info_matrix(df, L_matrix = L) + + expect_false( + all(abs(info_plain$eigenvalues - info_spatial$eigenvalues) < 1e-10) + ) +}) + +# disconnected design should give invalid criteria +test_that("disconnected design gives Inf for A and D values", { + df <- data.frame( + row = 1:6, + col = c(1, 2, 1, 2, 1, 2), + block = c(1, 1, 2, 2, 3, 3), + treatment = c("A", "A", "B", "B", "C", "C") + ) + + info <- calc_info_matrix(df) + + expect_equal(info$A_value, Inf) + expect_equal(info$D_value, Inf) +}) + + +# ---- calc_incidence_matrix + +# complete design should have all incidence matrix entries = 1 +# should also encode replication (rows) and block size (cols) +test_that("incidence matrix has correct dimensions and values for RCBD", { + df <- data.frame( + row = rep(1:3, each = 4), + col = rep(1:4, times = 3), + block = rep(1:3, each = 4), + treatment = rep(paste0("T", 1:4), 3) + ) + + N <- calc_incidence_matrix(df) + + expect_equal(dim(N), c(4, 3)) # v x b + expect_true(all(N == 1)) # complete blocks + expect_true(all(rowSums(N) == 3)) # replication = b + expect_true(all(colSums(N) == 4)) # block size = v +}) + +# ensure incidence works fine for imbalanced designs +test_that("incidence matrix handles unequal replication", { + df <- data.frame( + row = 1:5, + col = c(1, 2, 1, 2, 1), + block = c(1, 1, 2, 2, 2), + treatment = c("A", "B", "A", "B", "A") + ) + + N <- calc_incidence_matrix(df) + + expect_equal(dim(N), c(2, 2)) + expect_equal(sum(N["A", ]), 3) # A appears 3 times + expect_equal(sum(N["B", ]), 2) # B appears 2 times +}) + +# ---- calc_concurrence_matrix + +# balanced RCBD with N all = 1, N Nt should have entries equal to block number +test_that("concurrence matrix is correct for RCBD", { + df <- data.frame( + row = rep(1:3, each = 4), + col = rep(1:4, times = 3), + block = rep(1:3, each = 4), + treatment = rep(paste0("T", 1:4), 3) + ) + + C <- calc_concurrence_matrix(df) + + expect_equal(dim(C), c(4, 4)) + expect_true(all(C == 3)) +}) + + +# ---- calculate_efficiency_factors + +# balanced RCBD has all efficiency factors = 1 +test_that("efficiency factors are 1 for complete blocks", { + df <- data.frame( + row = rep(1:4, each = 5), + col = rep(1:5, times = 4), + block = rep(1:4, each = 5), + treatment = rep(paste0("T", 1:5), 4) + ) + + info <- calc_info_matrix(df) + treatments <- factor(df$treatment) + r_bar <- mean(table(treatments)) + eff <- info$eigenvalues / r_bar + + expect_true(all(abs(eff - 1) < 1e-10)) +}) + + +# ---- correlation constructors + +# make sure ar1 builder works fine +test_that("cor_ar1 produces valid correlation matrix", { + R <- cor_ar1(5, 0.7) + + expect_equal(dim(R), c(5, 5)) + expect_true(isSymmetric(R)) + expect_true(all(diag(R) == 1)) + expect_true(all(eigen(R, only.values = TRUE)$values > 0)) + + # Check specific entries + expect_equal(R[1, 2], 0.7) + expect_equal(R[1, 3], 0.7^2) + expect_equal(R[1, 5], 0.7^4) +}) + +# yes +test_that("cor_ar1 with rho = 0 gives identity", { + expect_equal(cor_ar1(4, 0), diag(4)) +}) + +# make sure ar1 x ar1 builder works fine +test_that("cor_ar1_ar1 has correct dimensions and is valid", { + Sigma <- cor_ar1_ar1(3, 5, 0.6, 0.4) + + expect_equal(dim(Sigma), c(15, 15)) + expect_true(isSymmetric(Sigma)) + expect_true(all(diag(Sigma) == 1)) + expect_true(all(eigen(Sigma, only.values = TRUE)$values > 0)) +}) + + +# ---- Test with speed function + +# Actually test by passing to speed +test_that("speed() optimises an RCBD with the info objective under AR1xAR1", { + skip_if_not_installed("speed") + + df <- initialise_design_df( + items = 6, nrows = 4, ncols = 6, + block_nrows = 1, block_ncols = 6 + ) + + Sigma <- cor_ar1_ar1(4, 6, rho_row = 0.6, rho_col = 0.3) + L <- compute_L_projection(df, Sigma, block_column = "block") + + result <- speed( + df, + swap = "treatment", + swap_within = "block", + spatial_factors = ~ block + row + col, + obj_function = objective_function_info, + L_matrix = L, + criterion = "A", + optimise_params = optim_params(random_initialisation = TRUE), + seed = 42, + quiet = TRUE + ) + + expect_s3_class(result, "design") + expect_lt(result$score, Inf) + + info <- calc_info_matrix(result$design_df, L_matrix = L) + expect_equal(info$rank, 5) + expect_true(is.finite(info$A_value)) +}) + +test_that("speed() optimises a BIBD with the info objective", { + skip_if_not_installed("speed") + + df <- initialise_design_df( + items = 5, nrows = 3, ncols = 10, + block_nrows = 3, block_ncols = 1 + ) + + result <- speed( + df, + swap = "treatment", + spatial_factors = ~ block + row + col, + obj_function = objective_function_info, + criterion = "A", + optimise_params = optim_params(random_initialisation = TRUE), + seed = 42, + quiet = TRUE + ) + + info <- calc_info_matrix(result$design_df) + expect_equal(info$rank, 4) + expect_true(is.finite(info$A_value)) +})