From 61f0ec9bf268541abcd1f8a62aa12dd9222a5a4b Mon Sep 17 00:00:00 2001 From: Kai Bagley Date: Tue, 28 Apr 2026 17:13:49 +0800 Subject: [PATCH 01/10] add everything --- NAMESPACE | 1 + R/objectives.R | 301 ++++++++++++++++++++++++++++ man/calc_concurrence_matrix.Rd | 15 ++ man/calc_incidence_matrix.Rd | 11 + man/calc_info_matrix.Rd | 16 ++ man/calculate_efficiency_factors.Rd | 16 ++ man/compute_L_projection.Rd | 11 + man/dot-build_L_from_df.Rd | 16 ++ man/dot-build_treatment_matrix.Rd | 11 + man/objective_function_info.Rd | 39 ++++ 10 files changed, 437 insertions(+) create mode 100644 R/objectives.R create mode 100644 man/calc_concurrence_matrix.Rd create mode 100644 man/calc_incidence_matrix.Rd create mode 100644 man/calc_info_matrix.Rd create mode 100644 man/calculate_efficiency_factors.Rd create mode 100644 man/compute_L_projection.Rd create mode 100644 man/dot-build_L_from_df.Rd create mode 100644 man/dot-build_treatment_matrix.Rd create mode 100644 man/objective_function_info.Rd diff --git a/NAMESPACE b/NAMESPACE index 889d8fa..c4c819e 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -16,6 +16,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..1e3701f --- /dev/null +++ b/R/objectives.R @@ -0,0 +1,301 @@ +#' 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 design 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 row_col Column name of the design's rows +#' @param range_col Column name of the design's range +#' @param criterion Either \code{"A"} or \code{"D"}, representing A or D +#' optimality. +#' - A-optimality: Minimises \eqn{\mathrm{tr} \left( A^- \right)}{Sigma^(-1)}. +#' - D-optimality: Minimises \eqn{-\log \left| I \right|} +#' @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_col Column name of the design's block factor. Used when +#' \code{L_matrix} isn't supplied. +#' +#' @details +#' This function computes the treatment information matrix: +#' +#' +#' @export +objective_function_info <- function( + design, + swap, + spatial_cols, + row_col = "row", + range_col = "range", + criterion = c("A", "D"), + L_matrix = NULL, + block_col = "block", + ... +) { + criterion <- match.arg(criterion) + + ## a lot of the args come from the design now + treatments <- design[[swap]] + trt_levels <- levels(factor(treatments)) + v <- length(trt_levels) + n <- nrow(design) + + ## L + ### projection matrix + if (is.null(L_matrix)) { + L_matrix <- .build_L_from_df(design, block_col, 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 + + ## get positive eigenvals from info matrix + eig <- eigen(info, symmetric = TRUE, only.values = TRUE)$values + pos_eig <- eig[eig > 1e-10] + + ## if design is disconnected, penalise greatly + ## otherwise A or D + if (length(pos_eig) < v - 1) { + score <- Inf + } else { + score <- switch(criterion, + A = sum(1 / pos_eig), + D = -sum(log(pos_eig)) + ) + } + + return(list( + score = score, + info_matrix = info, + eigenvalues = sort(pos_eig, decreasing = TRUE), + criterion = criterion + )) +} + +#' Build the treatment indicator matrix X1 +.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. +.build_L_from_df <- function(design, block_col, n) { + if (!block_col %in% names(design)) { + ## L = I - (1/n) ones_matrix + return(diag(n) - matrix(1 / n, n, n)) + } + + blocks <- factor(design[[block_col]]) + 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) +} + + +# Helper functions + +#' Compute L projection +compute_L_projection <- function(design, Sigma, block_col = "block") { + n <- nrow(design) + + ## 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 design dataframe") + if (ncol(Sigma) != n) stop("Sigma matrix must have dimension equal to dimension of design dataframe") + if (!isSymmetric(Sigma, tol = 1e-8)) stop("Sigma must be symmetric") + if (!(block_col %in% names(design))) stop("Block column name not in design dataframe") + + blocks <- factor(design[[block_col]]) + 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 +calc_info_matrix <- function( + design, + treatment_col = "treatment", + L_matrix = NULL, + block_col = "block" +) { + treatments <- design[[treatment_col]] + trt_levels <- levels(factor(treatments)) + v <- length(trt_levels) + n <- nrow(design) + + # TODO: Duplicated from above. Probably try and simplify later + ## L + ### projection matrix + if (is.null(L_matrix)) { + L_matrix <- .build_L_from_df(design, block_col, 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 + + rownames(info) <- trt_levels + colnames(info) <- trt_levels + + ## get positive eigenvals from info matrix + eig <- eigen(info, symmetric = TRUE, only.values = TRUE)$values + pos_eig <- eig[eig > 1e-10] + + ## if design is disconnected, penalise greatly + ## otherwise A or D + 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)) + } + + 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 +calc_incidence_matrix <- function( + design, + treatment_col = "treatment", + block_col = "block" +) { + treatments <- factor(design[[treatment_col]]) + blocks <- factor(design[[block_col]]) + + ## make a matrix of integers + N <- table(treatments, blocks) + class(N) <- "matrix" + storage.mode(N) <- "integer" + + return(N) +} + +#' Calculate concurrence matrix +calc_concurrence_matrix <- function( + design, + treatment_col = "treatment", + block_col = "block" +) { + N <- calc_incidence_matrix(design, treatment_col, block_col) + C <- N %*% t(N) + storage.mode(C) <- "integer" + return(C) +} + +#' Calculate canonical efficiency factors +calculate_efficiency_factors <- function( + design, + treatment_col = "treatment", + L_matrix = NULL, + block_col = "block" +) { + info <- calc_info_matrix( + design, + treatment_col, + L_matrix, + block_col + ) + + treatments <- factor(design[[treatment_col]]) + 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 +cor_ar1 <- function(n, rho) { + rho^abs(outer(seq_len(n), seq_len(n), "-")) +} + +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) + ) +} + + +## Testing +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 +} diff --git a/man/calc_concurrence_matrix.Rd b/man/calc_concurrence_matrix.Rd new file mode 100644 index 0000000..308cae3 --- /dev/null +++ b/man/calc_concurrence_matrix.Rd @@ -0,0 +1,15 @@ +% 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( + design, + treatment_col = "treatment", + block_col = "block" +) +} +\description{ +Calculate concurrence matrix +} diff --git a/man/calc_incidence_matrix.Rd b/man/calc_incidence_matrix.Rd new file mode 100644 index 0000000..6bf2ddb --- /dev/null +++ b/man/calc_incidence_matrix.Rd @@ -0,0 +1,11 @@ +% 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(design, treatment_col = "treatment", block_col = "block") +} +\description{ +Calculate incidence matrix +} diff --git a/man/calc_info_matrix.Rd b/man/calc_info_matrix.Rd new file mode 100644 index 0000000..120ee04 --- /dev/null +++ b/man/calc_info_matrix.Rd @@ -0,0 +1,16 @@ +% 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( + design, + treatment_col = "treatment", + L_matrix = NULL, + block_col = "block" +) +} +\description{ +Calculate info matrix +} diff --git a/man/calculate_efficiency_factors.Rd b/man/calculate_efficiency_factors.Rd new file mode 100644 index 0000000..fc3a1b9 --- /dev/null +++ b/man/calculate_efficiency_factors.Rd @@ -0,0 +1,16 @@ +% 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( + design, + treatment_col = "treatment", + L_matrix = NULL, + block_col = "block" +) +} +\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..6dc4353 --- /dev/null +++ b/man/compute_L_projection.Rd @@ -0,0 +1,11 @@ +% 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(design, Sigma, block_col = "block") +} +\description{ +Compute L projection +} diff --git a/man/dot-build_L_from_df.Rd b/man/dot-build_L_from_df.Rd new file mode 100644 index 0000000..7b0b379 --- /dev/null +++ b/man/dot-build_L_from_df.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/objectives.R +\name{.build_L_from_df} +\alias{.build_L_from_df} +\title{Function to precompute \eqn{L}} +\usage{ +.build_L_from_df(design, block_col, n) +} +\description{ +\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}. +} +\details{ +This only depends on lambda and the block structure, so precompute this. +} diff --git a/man/dot-build_treatment_matrix.Rd b/man/dot-build_treatment_matrix.Rd new file mode 100644 index 0000000..c785382 --- /dev/null +++ b/man/dot-build_treatment_matrix.Rd @@ -0,0 +1,11 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/objectives.R +\name{.build_treatment_matrix} +\alias{.build_treatment_matrix} +\title{Build the treatment indicator matrix X1} +\usage{ +.build_treatment_matrix(treatments, trt_levels, n, v) +} +\description{ +Build the treatment indicator matrix X1 +} diff --git a/man/objective_function_info.Rd b/man/objective_function_info.Rd new file mode 100644 index 0000000..0642a50 --- /dev/null +++ b/man/objective_function_info.Rd @@ -0,0 +1,39 @@ +% 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( + design, + swap, + spatial_cols, + current_score_obj = NULL, + swapped_items = NULL, + row_col = "row", + col_col = "col", + criterion = c("A", "D"), + L_matrix = NULL, + block_col = "block", + ... +) +} +\arguments{ +\item{design}{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. + +Information matrix based objectives +\itemize{ +\item A-optimality: \eqn{\mathrm{tr} \left( A^- \right)}{Sigma^(-1)} +\item D-optimality: \eqn{-\log \left| I \right|} +}} +} +\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. +} From 7efaf560bbbc7b5302b4bbc88381d41b01fc5a6c Mon Sep 17 00:00:00 2001 From: Kai Bagley Date: Thu, 30 Apr 2026 12:01:20 +0800 Subject: [PATCH 02/10] update tests and add some comments --- tests/testthat/test-objective_function_info.R | 317 ++++++++++++++++++ 1 file changed, 317 insertions(+) create mode 100644 tests/testthat/test-objective_function_info.R diff --git a/tests/testthat/test-objective_function_info.R b/tests/testthat/test-objective_function_info.R new file mode 100644 index 0000000..6a64ed6 --- /dev/null +++ b/tests/testthat/test-objective_function_info.R @@ -0,0 +1,317 @@ +# 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 return Inf", { + 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_equal(result$score, Inf) +}) + + +# ---- 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)) +}) From d827d90a29163ad834fe9d1eb4a7be079e44e748 Mon Sep 17 00:00:00 2001 From: Kai Bagley Date: Thu, 30 Apr 2026 12:03:22 +0800 Subject: [PATCH 03/10] add Kai and Curtin Uni to the description --- DESCRIPTION | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index df08a3a..fc034de 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, From d678f00239525bf2ecbb435934c34a6c25715e69 Mon Sep 17 00:00:00 2001 From: Kai Bagley Date: Thu, 30 Apr 2026 12:45:40 +0800 Subject: [PATCH 04/10] update roxygen2 stuff for exported functions --- DESCRIPTION | 2 +- R/objectives.R | 76 +++++++++++++++++++++++++++++++--- man/objective_function_info.Rd | 31 ++++++++++---- man/speed-package.Rd | 2 + 4 files changed, 95 insertions(+), 16 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index fc034de..4038760 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -29,7 +29,7 @@ Authors@R: comment = "https://adelaide.edu.au/"), person(given = "Grains Research and Development Corporation", role = c("cph", "fnd"), - comment = "https://grdc.com.au/") + comment = "https://grdc.com.au/"), person(given = "Curtin University", role = c("cph"), comment = c("http://www.curtin.edu.au/", ROR = "02n415q13"))) diff --git a/R/objectives.R b/R/objectives.R index 1e3701f..2e06dc1 100644 --- a/R/objectives.R +++ b/R/objectives.R @@ -8,12 +8,10 @@ #' @param design 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 row_col Column name of the design's rows -#' @param range_col Column name of the design's range #' @param criterion Either \code{"A"} or \code{"D"}, representing A or D #' optimality. -#' - A-optimality: Minimises \eqn{\mathrm{tr} \left( A^- \right)}{Sigma^(-1)}. -#' - D-optimality: Minimises \eqn{-\log \left| I \right|} +#' - 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 @@ -23,7 +21,14 @@ #' #' @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. autoregressive-lag-1 \eqn{\text{AR}_1 \otimes \text{AR}_1}{AR₁ ⊗ AR₁}, +#' and this objective will optimise under that structure. #' #' @export objective_function_info <- function( @@ -121,6 +126,15 @@ objective_function_info <- function( # Helper functions #' Compute L projection +#' +#' This is the projection matrix that removes nuisance effects from the GLS's +#' covariance +#' +#' @param design Data frame representing the spatial layout of the experiment. +#' @param Sigma Covariance structure to use. +#' @param block_col Column name of the design's block factor in \code{design}. +#' +#' @export compute_L_projection <- function(design, Sigma, block_col = "block") { n <- nrow(design) @@ -147,6 +161,17 @@ compute_L_projection <- function(design, Sigma, block_col = "block") { } #' Calculate info matrix +#' +#' Calculate the Fisher information in the experimental design +#' +#' @param design Data frame representing the spatial layout of the experiment. +#' @param treatment_col Column name containing the design's treatments in +#' \code{design}. +#' @param L_matrix Precomputed projection matrix from +#' \code{compute_L_projection}. +#' @param block_col Column name containing the designs block in \code{design}. +#' +#' @export calc_info_matrix <- function( design, treatment_col = "treatment", @@ -199,6 +224,13 @@ calc_info_matrix <- function( } #' Calculate incidence matrix +#' +#' @param design Data frame representing the spatial layout of the experiment. +#' @param treatment_col Column name containing the design's treatments in +#' \code{design}. +#' @param block_col Column name containing the designs block in \code{design}. +#' +#' @export calc_incidence_matrix <- function( design, treatment_col = "treatment", @@ -216,6 +248,13 @@ calc_incidence_matrix <- function( } #' Calculate concurrence matrix +#' +#' @param design Data frame representing the spatial layout of the experiment. +#' @param treatment_col Column name containing the design's treatments in +#' \code{design}. +#' @param block_col Column name containing the designs block in \code{design}. +#' +#' @export calc_concurrence_matrix <- function( design, treatment_col = "treatment", @@ -228,6 +267,15 @@ calc_concurrence_matrix <- function( } #' Calculate canonical efficiency factors +#' +#' @param design Data frame representing the spatial layout of the experiment. +#' @param treatment_col Column name containing the design's treatments in +#' \code{design}. +#' @param L_matrix Precomputed projection matrix from +#' \code{compute_L_projection}. +#' @param block_col Column name containing the designs block in \code{design}. +#' +#' @export calculate_efficiency_factors <- function( design, treatment_col = "treatment", @@ -262,11 +310,27 @@ calculate_efficiency_factors <- function( } -## Convenience correlation structure constructors +# 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 +#' +#' @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. +#' +#' @export cor_ar1_ar1 <- function(n_rows, n_cols, rho_row, rho_col) { kronecker( cor_ar1(n_rows, rho_row), diff --git a/man/objective_function_info.Rd b/man/objective_function_info.Rd index 0642a50..a884b94 100644 --- a/man/objective_function_info.Rd +++ b/man/objective_function_info.Rd @@ -8,10 +8,8 @@ objective_function_info( design, swap, spatial_cols, - current_score_obj = NULL, - swapped_items = NULL, row_col = "row", - col_col = "col", + range_col = "range", criterion = c("A", "D"), L_matrix = NULL, block_col = "block", @@ -23,13 +21,24 @@ objective_function_info( \item{swap}{Column name to swap, usually the treatment.} -\item{spatial_cols}{Column name of the spatial factors. +\item{spatial_cols}{Column name of the spatial factors.} -Information matrix based objectives -\itemize{ -\item A-optimality: \eqn{\mathrm{tr} \left( A^- \right)}{Sigma^(-1)} -\item D-optimality: \eqn{-\log \left| I \right|} -}} +\item{row_col}{Column name of the design's rows} + +\item{range_col}{Column name of the design's range} + +\item{criterion}{Either \code{"A"} or \code{"D"}, representing A or D +optimality. +- A-optimality: Minimises \eqn{\mathrm{tr} \left( A^- \right)}{tr(A^-)}. +- D-optimality: Minimises \eqn{-\log \left| I \right|}{-log(det(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_col}{Column name of the design's block factor. Used when +\code{L_matrix} isn't supplied.} } \description{ Creates an objective function that optimises experimental designs using the @@ -37,3 +46,7 @@ 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₁} +} 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] } } From f48a670f8438f312da6e5895bd650041ab80cf29 Mon Sep 17 00:00:00 2001 From: Kai Bagley Date: Fri, 1 May 2026 11:15:18 +0800 Subject: [PATCH 05/10] add emacs dir locals --- .Rbuildignore | 1 + 1 file changed, 1 insertion(+) 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$ From b31c28379810dbc18c8a6cdad79d31144fffe407 Mon Sep 17 00:00:00 2001 From: Kai Bagley Date: Fri, 1 May 2026 11:15:38 +0800 Subject: [PATCH 06/10] many fixes so examples work, and other issues lel --- R/objectives.R | 537 +++++++++++++++++++++++++++---------------------- 1 file changed, 299 insertions(+), 238 deletions(-) diff --git a/R/objectives.R b/R/objectives.R index 2e06dc1..59a01b5 100644 --- a/R/objectives.R +++ b/R/objectives.R @@ -5,7 +5,7 @@ #' effects, and takes an optional spatial covariance structure for spatial #' optimisation. #' -#' @param design Data frame representing the spatial layout of the experiment. +#' @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 @@ -16,8 +16,9 @@ #' \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_col Column name of the design's block factor. Used when +#' @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: @@ -25,75 +26,109 @@ #' 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₂)⁻¹} +#' \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. autoregressive-lag-1 \eqn{\text{AR}_1 \otimes \text{AR}_1}{AR₁ ⊗ AR₁}, -#' and this objective will optimise under that structure. +#' 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( - design, - swap, - spatial_cols, - row_col = "row", - range_col = "range", - criterion = c("A", "D"), - L_matrix = NULL, - block_col = "block", - ... + layout_df, + swap, + spatial_cols, + criterion = c("A", "D"), + L_matrix = NULL, + block_column = "block", + ... ) { - criterion <- match.arg(criterion) - - ## a lot of the args come from the design now - treatments <- design[[swap]] - trt_levels <- levels(factor(treatments)) - v <- length(trt_levels) - n <- nrow(design) - - ## L - ### projection matrix - if (is.null(L_matrix)) { - L_matrix <- .build_L_from_df(design, block_col, 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 - - ## get positive eigenvals from info matrix - eig <- eigen(info, symmetric = TRUE, only.values = TRUE)$values - pos_eig <- eig[eig > 1e-10] - + 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 <- Inf + score <- 1e10 } else { - score <- switch(criterion, - A = sum(1 / pos_eig), - D = -sum(log(pos_eig)) - ) + 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 - )) + } + + 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) + 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} @@ -103,25 +138,54 @@ objective_function_info <- function( #' \eqn{\Sigma}{Sigma} and nuisance effects design \eqn{X_2}. #' #' This only depends on lambda and the block structure, so precompute this. -.build_L_from_df <- function(design, block_col, n) { - if (!block_col %in% names(design)) { - ## L = I - (1/n) ones_matrix - return(diag(n) - matrix(1 / n, n, n)) - } - - blocks <- factor(design[[block_col]]) - 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) +#' @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 @@ -130,183 +194,201 @@ objective_function_info <- function( #' This is the projection matrix that removes nuisance effects from the GLS's #' covariance #' -#' @param design Data frame representing the spatial layout of the experiment. +#' @param layout_df Data frame representing the spatial layout of the experiment. #' @param Sigma Covariance structure to use. -#' @param block_col Column name of the design's block factor in \code{design}. +#' @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(design, Sigma, block_col = "block") { - n <- nrow(design) - - ## 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 design dataframe") - if (ncol(Sigma) != n) stop("Sigma matrix must have dimension equal to dimension of design dataframe") - if (!isSymmetric(Sigma, tol = 1e-8)) stop("Sigma must be symmetric") - if (!(block_col %in% names(design))) stop("Block column name not in design dataframe") - - blocks <- factor(design[[block_col]]) - 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)) +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 design Data frame representing the spatial layout of the experiment. -#' @param treatment_col Column name containing the design's treatments in -#' \code{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_col Column name containing the designs block in \code{design}. +#' @param block_column Column name containing the designs block in +#' \code{layout_df}. #' #' @export calc_info_matrix <- function( - design, - treatment_col = "treatment", - L_matrix = NULL, - block_col = "block" + layout_df, + treatment_column = "treatment", + L_matrix = NULL, + block_column = "block" ) { - treatments <- design[[treatment_col]] - trt_levels <- levels(factor(treatments)) - v <- length(trt_levels) - n <- nrow(design) - - # TODO: Duplicated from above. Probably try and simplify later - ## L - ### projection matrix - if (is.null(L_matrix)) { - L_matrix <- .build_L_from_df(design, block_col, 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 - - rownames(info) <- trt_levels - colnames(info) <- trt_levels - - ## get positive eigenvals from info matrix - eig <- eigen(info, symmetric = TRUE, only.values = TRUE)$values - pos_eig <- eig[eig > 1e-10] - - ## if design is disconnected, penalise greatly - ## otherwise A or D - 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)) - } - - return(list( - info_matrix = info, - eigenvalues = sort(pos_eig, decreasing = TRUE), - rank = length(pos_eig), - A_value = A_val, - D_value = D_val - )) + 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 design Data frame representing the spatial layout of the experiment. -#' @param treatment_col Column name containing the design's treatments in -#' \code{design}. -#' @param block_col Column name containing the designs block in \code{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 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( - design, - treatment_col = "treatment", - block_col = "block" + layout_df, + treatment_column = "treatment", + block_column = "block" ) { - treatments <- factor(design[[treatment_col]]) - blocks <- factor(design[[block_col]]) + 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" + ## make a matrix of integers + N <- table(treatments, blocks) + class(N) <- "matrix" + storage.mode(N) <- "integer" - return(N) + return(N) } #' Calculate concurrence matrix #' -#' @param design Data frame representing the spatial layout of the experiment. -#' @param treatment_col Column name containing the design's treatments in -#' \code{design}. -#' @param block_col Column name containing the designs block in \code{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 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( - design, - treatment_col = "treatment", - block_col = "block" + layout_df, + treatment_column = "treatment", + block_column = "block" ) { - N <- calc_incidence_matrix(design, treatment_col, block_col) - C <- N %*% t(N) - storage.mode(C) <- "integer" - return(C) + 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 design Data frame representing the spatial layout of the experiment. -#' @param treatment_col Column name containing the design's treatments in -#' \code{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_col Column name containing the designs block in \code{design}. +#' @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( - design, - treatment_col = "treatment", - L_matrix = NULL, - block_col = "block" + layout_df, + treatment_column = "treatment", + L_matrix = NULL, + block_column = "block" ) { - info <- calc_info_matrix( - design, - treatment_col, - L_matrix, - block_col - ) - - treatments <- factor(design[[treatment_col]]) - 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)) - ) + 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)) + ) } @@ -317,6 +399,8 @@ calculate_efficiency_factors <- function( #' @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), "-")) @@ -330,36 +414,13 @@ cor_ar1 <- function(n, rho) { #' @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) - ) -} - - -## Testing -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 + kronecker( + cor_ar1(n_rows, rho_row), + cor_ar1(n_cols, rho_col) + ) } From 35a6d3ff0ac9df158aad4a94f1a4f365c1f86a70 Mon Sep 17 00:00:00 2001 From: Kai Bagley Date: Fri, 1 May 2026 11:15:58 +0800 Subject: [PATCH 07/10] fix test with new penalty --- tests/testthat/test-objective_function_info.R | 479 ++++++++++-------- 1 file changed, 270 insertions(+), 209 deletions(-) diff --git a/tests/testthat/test-objective_function_info.R b/tests/testthat/test-objective_function_info.R index 6a64ed6..9c89379 100644 --- a/tests/testthat/test-objective_function_info.R +++ b/tests/testthat/test-objective_function_info.R @@ -4,109 +4,110 @@ # 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 + 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) + 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) + 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 + 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 return Inf", { - 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_equal(result$score, Inf) +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)) }) @@ -114,35 +115,35 @@ test_that("disconnected designs return Inf", { # 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") - ) + 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))) + # Wrong dimensions + expect_error(compute_L_projection(df, diag(3))) - # Not symmetric - expect_error(compute_L_projection(df, matrix(1:16, 4, 4))) + # 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")) + # 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) - ) + 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) + 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) + expect_true(max(abs(L_spatial - L_direct)) < 1e-10) }) @@ -150,63 +151,63 @@ test_that("L with identity Sigma matches .build_L_from_df", { # 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) - ) + 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) + info <- calc_info_matrix(df) - # v = 5, so rank = v - 1 = 4 - expect_equal(info$rank, 4) + # 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)) + # 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)) + # 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)) + # 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) - ) + 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") - ) + 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) + info <- calc_info_matrix(df) - expect_equal(info$A_value, Inf) - expect_equal(info$D_value, Inf) + expect_equal(info$A_value, Inf) + expect_equal(info$D_value, Inf) }) @@ -215,52 +216,52 @@ test_that("disconnected design gives Inf for A and D values", { # 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 + 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 + 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) - ) + 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) + C <- calc_concurrence_matrix(df) - expect_equal(dim(C), c(4, 4)) - expect_true(all(C == 3)) + expect_equal(dim(C), c(4, 4)) + expect_true(all(C == 3)) }) @@ -268,19 +269,19 @@ test_that("concurrence matrix is correct for RCBD", { # 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)) + 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)) }) @@ -288,30 +289,90 @@ test_that("efficiency factors are 1 for complete blocks", { # make sure ar1 builder works fine test_that("cor_ar1 produces valid correlation matrix", { - R <- cor_ar1(5, 0.7) + 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)) + 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) + # 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)) + 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) + 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)) +}) - 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_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)) }) From 9275f8d4b745b2ae986010e8533586907df981bb Mon Sep 17 00:00:00 2001 From: Kai Bagley Date: Fri, 1 May 2026 11:16:09 +0800 Subject: [PATCH 08/10] run devtools::check() --- NAMESPACE | 7 +++ man/calc_concurrence_matrix.Rd | 19 ++++++-- man/calc_incidence_matrix.Rd | 19 +++++++- man/calc_info_matrix.Rd | 20 ++++++-- man/calculate_efficiency_factors.Rd | 21 +++++++-- man/compute_L_projection.Rd | 15 +++++- man/dot-build_L_from_df.Rd | 16 ------- man/dot-build_treatment_matrix.Rd | 11 ----- man/objective_function_info.Rd | 73 ++++++++++++++++++++++++----- 9 files changed, 149 insertions(+), 52 deletions(-) delete mode 100644 man/dot-build_L_from_df.Rd delete mode 100644 man/dot-build_treatment_matrix.Rd diff --git a/NAMESPACE b/NAMESPACE index c4c819e..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) diff --git a/man/calc_concurrence_matrix.Rd b/man/calc_concurrence_matrix.Rd index 308cae3..463ee5c 100644 --- a/man/calc_concurrence_matrix.Rd +++ b/man/calc_concurrence_matrix.Rd @@ -5,11 +5,24 @@ \title{Calculate concurrence matrix} \usage{ calc_concurrence_matrix( - design, - treatment_col = "treatment", - block_col = "block" + 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 index 6bf2ddb..8f113bb 100644 --- a/man/calc_incidence_matrix.Rd +++ b/man/calc_incidence_matrix.Rd @@ -4,7 +4,24 @@ \alias{calc_incidence_matrix} \title{Calculate incidence matrix} \usage{ -calc_incidence_matrix(design, treatment_col = "treatment", block_col = "block") +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 index 120ee04..5887f42 100644 --- a/man/calc_info_matrix.Rd +++ b/man/calc_info_matrix.Rd @@ -5,12 +5,24 @@ \title{Calculate info matrix} \usage{ calc_info_matrix( - design, - treatment_col = "treatment", + layout_df, + treatment_column = "treatment", L_matrix = NULL, - block_col = "block" + 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 info matrix +Calculate the Fisher information in the experimental design } diff --git a/man/calculate_efficiency_factors.Rd b/man/calculate_efficiency_factors.Rd index fc3a1b9..cadc96f 100644 --- a/man/calculate_efficiency_factors.Rd +++ b/man/calculate_efficiency_factors.Rd @@ -5,12 +5,27 @@ \title{Calculate canonical efficiency factors} \usage{ calculate_efficiency_factors( - design, - treatment_col = "treatment", + layout_df, + treatment_column = "treatment", L_matrix = NULL, - block_col = "block" + 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 index 6dc4353..0d401ea 100644 --- a/man/compute_L_projection.Rd +++ b/man/compute_L_projection.Rd @@ -4,8 +4,19 @@ \alias{compute_L_projection} \title{Compute L projection} \usage{ -compute_L_projection(design, Sigma, block_col = "block") +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{ -Compute L projection +This is the projection matrix that removes nuisance effects from the GLS's +covariance } diff --git a/man/dot-build_L_from_df.Rd b/man/dot-build_L_from_df.Rd deleted file mode 100644 index 7b0b379..0000000 --- a/man/dot-build_L_from_df.Rd +++ /dev/null @@ -1,16 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/objectives.R -\name{.build_L_from_df} -\alias{.build_L_from_df} -\title{Function to precompute \eqn{L}} -\usage{ -.build_L_from_df(design, block_col, n) -} -\description{ -\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}. -} -\details{ -This only depends on lambda and the block structure, so precompute this. -} diff --git a/man/dot-build_treatment_matrix.Rd b/man/dot-build_treatment_matrix.Rd deleted file mode 100644 index c785382..0000000 --- a/man/dot-build_treatment_matrix.Rd +++ /dev/null @@ -1,11 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/objectives.R -\name{.build_treatment_matrix} -\alias{.build_treatment_matrix} -\title{Build the treatment indicator matrix X1} -\usage{ -.build_treatment_matrix(treatments, trt_levels, n, v) -} -\description{ -Build the treatment indicator matrix X1 -} diff --git a/man/objective_function_info.Rd b/man/objective_function_info.Rd index a884b94..0f0993d 100644 --- a/man/objective_function_info.Rd +++ b/man/objective_function_info.Rd @@ -5,40 +5,40 @@ \title{Objective function using the treatment information matrix} \usage{ objective_function_info( - design, + layout_df, swap, spatial_cols, - row_col = "row", - range_col = "range", criterion = c("A", "D"), L_matrix = NULL, - block_col = "block", + block_column = "block", ... ) } \arguments{ -\item{design}{Data frame representing the spatial layout of the experiment.} +\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{row_col}{Column name of the design's rows} - -\item{range_col}{Column name of the design's range} - \item{criterion}{Either \code{"A"} or \code{"D"}, representing A or D optimality. -- A-optimality: Minimises \eqn{\mathrm{tr} \left( A^- \right)}{tr(A^-)}. -- D-optimality: Minimises \eqn{-\log \left| I \right|}{-log(det(I))}} +- 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_col}{Column name of the design's block factor. Used when +\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 @@ -49,4 +49,53 @@ 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) + } From c292872895f9392a7a6fb439d323dd1fc540753a Mon Sep 17 00:00:00 2001 From: Kai Bagley Date: Fri, 1 May 2026 11:16:38 +0800 Subject: [PATCH 09/10] fix indenting --- R/objectives.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/objectives.R b/R/objectives.R index 59a01b5..f4bbc3d 100644 --- a/R/objectives.R +++ b/R/objectives.R @@ -403,7 +403,7 @@ calculate_efficiency_factors <- function( #' #' @export cor_ar1 <- function(n, rho) { - rho^abs(outer(seq_len(n), seq_len(n), "-")) + rho^abs(outer(seq_len(n), seq_len(n), "-")) } From 179b4847e08f83d5c2b85e9c1625cbfbf6e695c6 Mon Sep 17 00:00:00 2001 From: Kai Bagley Date: Fri, 1 May 2026 11:18:16 +0800 Subject: [PATCH 10/10] fix some lintr stuff --- R/objectives.R | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/R/objectives.R b/R/objectives.R index f4bbc3d..9fd72cb 100644 --- a/R/objectives.R +++ b/R/objectives.R @@ -5,13 +5,16 @@ #' effects, and takes an optional spatial covariance structure for spatial #' optimisation. #' -#' @param layout_df Data frame representing the spatial layout of the experiment. +#' @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|)} +#' - 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 @@ -194,9 +197,11 @@ objective_function_info <- function( #' 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 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}. +#' @param block_column Column name of the design's block factor in +#' \code{layout_df}. #' #' @return \eqn{(n \times n)} numeric matrix. #' @@ -245,7 +250,8 @@ compute_L_projection <- function( #' #' Calculate the Fisher information in the experimental design #' -#' @param layout_df Data frame representing the spatial layout of the experiment. +#' @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 @@ -296,7 +302,8 @@ calc_info_matrix <- function( #' Calculate incidence matrix #' -#' @param layout_df Data frame representing the spatial layout of the experiment. +#' @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 @@ -347,7 +354,8 @@ calc_concurrence_matrix <- function( #' Calculate canonical efficiency factors #' -#' @param layout_df Data frame representing the spatial layout of the experiment. +#' @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