Remembrance Day

Lest we forget

1:00

Plan for today

  • Setup check
  • Why spatial optimisation?
  • How speed works
  • Spatially optimised common designs
  • Evaluating designs
  • Designs with more complex structure
  • MET designs
  • Extensions (relationship matrices, custom objectives)
  • Gotchas, getting help, Q&A

Setup check

Getting set up

# install.packages("pak")
pak::pak("biometryhub/speed")



I will assume a working knowledge of experimental design, but please ask if something isn’t clear.

Does it work for you?

library(speed)
d <- initialise_design_df(4, 4, 4)
autoplot(speed(d, 
      swap = "treatment", 
      seed = 1))
Optimising level: single treatment within whole design 
Optimal score reached at iteration 787 for level single treatment within whole design

Why spatial optimisation?

Field heterogeneity and neighbours

  • Spatial trends are well known to exist in field trials1
  • Neighbour effects impact treatments
  • Spatially optimised designs guard against a sub-optimal randomisation
  • Alternatives are not FOSS and/or lack documentation and/or features

Restricted randomisation

  • Impose some restrictions \(\rightarrow\) improve spatial balance of treatments.
  • Optimisation problem…

Systematic

Fully random

Restricted randomisation

How speed works

The workflow

  1. Set up the layout (initialise_design_df() or otherwise)
  2. speed(): optimise the layout
  3. autoplot(): visualise it
  4. summary(): evaluate it

Four questions to ask of every design

Question Argument
What moves? swap, linked_cols*
Where is it allowed to move to? swap_within
What should it be balanced across? spatial_factors
What counts as a ‘neighbour’? grid_factors + ...

*: optional

The default objective

  • Adjacency: penalises the same treatment in neighbouring plots (what counts as a neighbour?)
  • Balance: penalises a treatment repeating in the same rows or columns (what should it be balanced across?)
  • Minimise the (optionally weighted) sum of these two components

Simulated annealing in one slide

%%{init: {'theme': 'base', 'themeVariables': {'fontSize': '22px', 'fontFamily': 'Arial', 'lineColor': '#414042'}}}%%
flowchart LR
  A(["<span style='display:inline-block;padding:14px 30px;font-size:30px;font-weight:bold'>Design</span>"]) --> B["Swap two plots,<br/>score the new design"]
  B --> C{"Better?"}
  C -- "yes" --> D["Accept"]
  C -- "no" --> E["Sometimes accept<br/>(less often as it cools)"]
  D --> A
  E --> A
  A -- "optimal score or<br/>no improvement" --> H(["Stop"])

  classDef hub fill:#00808b,stroke:#004d54,stroke-width:5px,color:#ffffff
  classDef design fill:#00808b,stroke:#00808b,color:#ffffff
  classDef step fill:#f2f2f2,stroke:#00808b,color:#414042
  classDef good fill:#b6d438,stroke:#54921e,color:#414042
  classDef maybe fill:#ffbc42,stroke:#ec8525,color:#414042
  class A hub
  class H design
  class B,C step
  class D good
  class E maybe

Spatially optimised common designs

Spatially balanced design1: set up the layout

sbd_design <- initialise_design_df(items = 8, nrows = 8, ncols = 4)
head(sbd_design)
  row col treatment
1   1   1        T1
2   2   1        T2
3   3   1        T3
4   4   1        T4
5   5   1        T5
6   6   1        T6

SBD: the four questions

Question Argument Input
What moves? swap treatment
Where is it allowed to move to? swap_within default
What should it be balanced across? spatial_factors default
What counts as a ‘neighbour’? grid_factors default

SBD: optimise

sbd_result <- speed(sbd_design,
                    swap = "treatment",
                    seed = 42)
Optimising level: single treatment within whole design 
Optimal score reached at iteration 730 for level single treatment within whole design 
sbd_result
Optimised Experimental Design
----------------------------
Score: 2.285714 
Iterations Run: 730 
Stopped Early: TRUE 
Treatments: T1, T2, T3, T4, T5, T6, T7, T8 
Seed: 42

SBD: plot

autoplot(sbd_result)

Spatially balanced block design: set up the layout

sbbd_design <- initialise_design_df(
  items = paste0("V", 1:6), 
  nrows = 4, ncols = 6,
  block_nrows = 1, block_ncols = 6
)

SBBD: the four questions

Question Argument Input
What moves? swap treatment
Where is it allowed to move to? swap_within block
What should it be balanced across? spatial_factors default
What counts as a ‘neighbour’? grid_factors default

SBBD: optimise

sbbd_result <- speed(sbbd_design,
                     swap = "treatment",
                     swap_within = "block",
                     seed = 42)
Optimising level: single treatment within block 
Optimal score reached at iteration 329 for level single treatment within block

SBBD layout

autoplot(sbbd_result)

Your turn

Exercise 1 (about 10 minutes, on the Exercises page)

  • 8 varieties × 4 replicates in a 4 × 8 field
  • Each block is a 2 × 4 patch, not a whole row
  • Answer the four questions, optimise, and check each block has every variety once

Evaluating designs

summary(): structure and optimisation

summary(sbbd_result), part 1 of 2

Design Summary
==============

Structure
---------
Layout:       4 rows x 6 cols (24 plots)
Treatments:   6
Replication:  4 each
Spatial:      row (4), col (6)

Optimisation
------------
Seed:         42
Objective:    objective_function
Score:        1.6  (initial 34 -> final 1.6)
              adjacency  0
              balance    1.6
Optimal:      1.6  (reached)
Iterations:   329 / 10,000 (optimal reached)
Temperature:  start 100, cooling 0.99

summary(): evaluation

summary(sbbd_result), part 2 of 2

Design Summary
==============

Evaluation
----------
Connected:    connected - treatment estimable given row + col + block [model (row + col + block)]
Concurrence:  complete blocks - not informative
Blk. spread:  each treatment in 4 of 4 blocks [block: block]
Repl. span:   worst-case 2 (row), 2 (col) across 6 replicated treatment(s)
Efficiency:   not requested (set efficiency = TRUE)
Self-adj.:    none
Neighbour:    min 1, max 4 over 15 pairs (variance 0.6952)

What the sections tell you

  • Structure: layout, replication, spatial factors
  • Optimisation: final score split into adjacency and balance
  • Evaluation: connectedness, concurrence, replicate span
  • Flags: iteration cap reached, unequal replication, disconnected design

Designs with more complex structure

More complex design types

  • Designs with hierarchy
  • Designs with multiple treatments
  • Designs with multiple sites

Split-plot: describe the structure

splits <- list(
  subplot   = list(items = letters[1:4], nrows = 1, ncols = 1),
  wholeplot = list(items = LETTERS[1:3], nrows = 1, ncols = 4),
  block     = list(nrows = 3, ncols = 4)
)
sp_design <- initialise_split_design_df(splits = splits, rep_dim = c(4, 1))

Split-plot: the four questions

One answer per level (wp, then sp):

Question Argument Whole plots Subplots
What moves? swap wholeplot_treatment subplot_treatment
Where? swap_within block wholeplot
Balance spatial_factors default default
Neighbour grid_factors default default

Split-plot: optimise level by level

optimise <- list(
  wp = list(
    swap        = "wholeplot_treatment",
    swap_within = "block",
    swap_all    = TRUE
  ),
  sp = list(
    swap        = "subplot_treatment",
    swap_within = "wholeplot"
  )
)

sp_result <- speed(sp_design, 
  optimise = optimise, 
  seed = 42)
Optimising level: wp 
Level: wp Iteration: 1000 Score: 100 Best: 100 Since Improvement: 1000 
Level: wp Iteration: 2000 Score: 100 Best: 100 Since Improvement: 2000 
Early stopping at iteration 2000 for level wp 
Optimising level: sp 
Optimal score reached at iteration 570 for level sp

Split-plot layout

p-rep: any data frame will do

prep_design <- data.frame(
  row   = rep(1:14, each = 10),
  col   = rep(1:10, times = 14),
  block = rep(rep(1:2, each = 5), 14),
  entry = c(rep("check", 10),           # 1 check, 10 replicates
            paste0("V", rep(1:30, 2)),  # 30 varieties, 2 replicates each
            paste0("V", 31:100))        # 70 varieties, 1 replicate each
)

p-rep: the four questions

Question Argument Input
What moves? swap entry
Where? swap_within default
Balance spatial_factors ~ block + row + col
Neighbour grid_factors default

p-rep: optimise

prep_result <- speed(prep_design,
                     swap = "entry",
                     spatial_factors = ~ block + row + col,
                     seed = 42)
Optimising level: single entry within whole design 
Level: single entry within whole design Iteration: 1000 Score: 3.41703 Best: 3.41703 Since Improvement: 0 
Level: single entry within whole design Iteration: 2000 Score: 3.29703 Best: 3.29703 Since Improvement: 27 
Level: single entry within whole design Iteration: 3000 Score: 3.29703 Best: 3.29703 Since Improvement: 1027 
Early stopping at iteration 3973 for level single entry within whole design

p-rep layout

autoplot(prep_result, treatments = "entry")

p-rep layout

Also covered in the vignettes

  • Split-split plot and strip-plot
  • Incomplete block designs
  • Row–column designs with blocking
  • linked_cols: keep variety names with their codes

Your turn

Exercise 2 (about 10 minutes, on the Exercises page)

  • 4 tillage methods (whole plots) × 3 varieties (subplots), in 3 blocks
  • Answer the four questions for each level, then optimise
  • Check every block has every tillage method once, and every whole plot every variety once

MET designs

Several sites, one treatment set

all_treatments <- c(rep(1:50, 7), rep(51:57, 8))
met_design <- initialise_design_df(
  items = all_treatments,
  designs = list(
    a = list(nrows = 28, ncols = 5, block_nrows = 7,  block_ncols = 5),
    b = list(nrows = 20, ncols = 3, block_nrows = 10, block_ncols = 3),
    c = list(nrows = 20, ncols = 4, block_nrows = 10, block_ncols = 4),
    d = list(nrows = 15, ncols = 4, block_nrows = 5,  block_ncols = 4),
    e = list(nrows = 22, ncols = 3, block_nrows = 11, block_ncols = 3)
  )
)
met_design$site_col   <- paste(met_design$site, met_design$col,   sep = "_")
met_design$site_block <- paste(met_design$site, met_design$block, sep = "_")

MET: the four questions

Two levels: first connect the sites, then balance within each site.

Question Argument Connectivity Balance
What moves? swap treatment treatment
Where? swap_within default site
Balance spatial_factors ~ site ~ site_col + site_block
Neighbour grid_factors by = "site" by = "site"

Optimise connectivity, then balance

MET Part 1

optimise <- list(
  connectivity = list(
    spatial_factors = ~ site
  ),
  balance = list(
    swap_within     = "site",
    spatial_factors = ~ site_col + site_block
  )
)

Optimise connectivity, then balance

MET Part 2

met_result <- speed(
  data = met_design,
  swap = "treatment",
  optimise = optimise,
  grid_factors = list(dim1 = "row", dim2 = "col", by = "site"),
  early_stop_iterations = 5000,
  optimise_params = optim_params(random_initialisation = TRUE, adj_weight = 0),
  seed = 112, quiet = TRUE
)

MET layout

Extensions

Relationship matrices

Treat related entries as partial repeats when scoring adjacency.

items  <- paste0("T", 1:16)
groups <- setNames(rep(c("A", "B", "C", "D"), each = 4), items)
rel <- matrix(0, 16, 16, dimnames = list(items, items))
within <- c(A = 0.8, B = 0.4, C = 0.1, D = 0.0)
for (g in names(within)) rel[items[groups == g], 
items[groups == g]] <- within[[g]]
diag(rel) <- 1
head(rel)
    T1  T2  T3  T4  T5  T6  T7  T8 T9 T10 T11 T12 T13 T14 T15 T16
T1 1.0 0.8 0.8 0.8 0.0 0.0 0.0 0.0  0   0   0   0   0   0   0   0
T2 0.8 1.0 0.8 0.8 0.0 0.0 0.0 0.0  0   0   0   0   0   0   0   0
T3 0.8 0.8 1.0 0.8 0.0 0.0 0.0 0.0  0   0   0   0   0   0   0   0
T4 0.8 0.8 0.8 1.0 0.0 0.0 0.0 0.0  0   0   0   0   0   0   0   0
T5 0.0 0.0 0.0 0.0 1.0 0.4 0.4 0.4  0   0   0   0   0   0   0   0
T6 0.0 0.0 0.0 0.0 0.4 1.0 0.4 0.4  0   0   0   0   0   0   0   0

Relationship matrices: optimise

rel_design <- initialise_design_df(items = rep(items, 4), nrows = 8, ncols = 8,
                                   block_nrows = 4, block_ncols = 4)
rel_result <- speed(rel_design, swap = "treatment", swap_within = "block",
                    relationship = rel, seed = 42, quiet = TRUE)

rel_result               
Optimised Experimental Design
----------------------------
Score: 4.266667 
Iterations Run: 3511 
Stopped Early: TRUE 
Treatments: T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11, T12, T13, T14, T15, T16 
Seed: 42

Relationship matrices: the effect

Tiles show each treatment’s relationship group (A: highly related, D: unrelated) not treatment.

Custom objectives

An objective returns a list with a score; lower is better.

get_cooccurrence <- function(design_df, swap) {
  pairs <- unlist(lapply(unique(design_df$block), function(block) {
    combn(sort(as.character(design_df[design_df$block == block, swap])), 2,
          FUN = function(x) paste(x, collapse = ","))
  }))
  return(as.list(table(pairs)))
}

bibd_objective <- function(layout_df, swap, spatial_cols,
                           row_column = "row", col_column = "col", ...) {
  adj <- calculate_adjacency_score(layout_df, swap, row_column, col_column)
  bal <- calculate_balance_score(layout_df, swap, spatial_cols)
  pair_var <- var(as.numeric(get_cooccurrence(layout_df, swap)))
  return(list(score = round(adj + bal + 10 * pair_var, 10)))
}

Custom objectives: a BIBD

bibd_df <- initialise_design_df(items = 5, nrows = 3, ncols = 10,
                                block_nrows = 3, block_ncols = 1)
bibd_result <- speed(bibd_df, swap = "treatment",
                     spatial_factors = ~ block + row + col,
                     obj_function = bibd_objective,
                     optimise_params = optim_params(random_initialisation = TRUE),
                     seed = 42, quiet = TRUE)
unique(unlist(get_cooccurrence(bibd_result$design_df, "treatment")))
[1] 3

Custom objectives: every pair together 3 times

T1 and T2 share a block exactly 3 times, like every other pair.

Gotchas, getting help, Q&A

Gotchas

  • Treatment column not called treatment? Pass treatments = to autoplot()
  • Companion columns (variety names) don’t move unless named in linked_cols
  • MET sites reuse coordinates: set grid_factors = list(..., by = "site")
  • Check summary() flags: a run that hit the iteration cap may need more iterations

Getting help

Questions?