blocksdesign (version 3.6)

design: General block and treatment designs.

Description

Constructs D-optimal block and treatment designs for any feasible combinations of nested or crossed block factors and any feasible linear treatment model.

Usage

design(treatments, blocks = NULL, treatments_model = NULL,
  weighting = 0.5, searches = NULL, seed = NULL, jumps = 1)

Arguments

treatments

a single treatment factor or data frame containing one or more qualitative or quantitative (numeric) level treatment factors.

blocks

a single block factor or data frame containing one or more qualitative level block factors in the required order of fitting.

treatments_model

a model formula for the required treatments design.

weighting

a weighting factor between 0 and 1 for weighting the interaction effects of crossed blocks factors where the default weighting is 0.5

searches

the maximum number of local optima searched at each stage of an optimization. The default depends on the design size.

seed

an integer initializing the random number generator. The null default gives an arbitrary random initialization.

jumps

the number of pairwise random treatment swaps used to escape a local maxima. The default is a single swap.

Value

treatments

The treatments included in the design and the replication of each individual treatment taken in de-randomized standard order.

design

The design layout showing the randomized allocation of treatments to blocks and plots.

treatments_model

The fitted treatment model, the number of model parameters (DF) and the D-efficiency of each sequentially fitted treatment model

blocks_model

The blocks sub-model design and the D- and A-efficiency factors of each successively fitted sub-blocks model.

seed

Numerical seed for random number generator.

searches

Maximum number of searches in each stratum.

jumps

Number of random treatment swaps to escape a local maxima.

Details

The treatments object is a factor or a data frame containing one or more qualitative or quantitative level treatment factors defining a set of candidate treatments from which the optimized treatment design is selected.

The blocks object is a factor or a data frame containing one or more qualitative level block factors. The default blocks object is a single level factor equal in length to the treatments object. The size of the blocks object defines the required design size.

The treatments_model can be either a single formula for a single design matrix or a compound formula for two or more design matrices. A compound formula contains one or more occurrences of a splitting operator | which splits-off a partial design formula on the left hand side of each |. Assuming the left hand part of each split is a well defined model design formula and replacing all remaining | by + in each partial design formula gives a hierarchical set of design matrices for sequential model fitting. The advantage of sequential model fitting is that it provides improved flexibility for fitting factors or variables of different status or importance and allows a wider range of choices of optimized design for different situations (see examples below).

Treatments are selected from the design candidate set with replacement unless the size of the candidate set exactly equals the size of the required design and the treatments_model is null, in which case the full candidate set is used for the treatment design. This option allows any arbitrary treatment set with any arbitrary treatment replication to be input as the required treatment design. If the size of the candidate set is different from the size of the required design, or if the treatments_model is non-null, the treatment design is optimized by selection with replacement.

The treatment design criterion is the ratio of the generalized variance of the treatment design based on the full treatment candidate set relative to the generalized variance of the treatment design based on the optimized treatment set (D-optimality). If the required design is a fractional factorial and the candidate set is a full factorial, the candidate set will be orthogonal and any design selected from the candidate set will have a relative efficiency less than or equal to 1. For quantitative level treatments, however, a full factorial design may not provide an optimal design and then the relative efficiency of the optimized design may exceed 1. The efficiency factor is used to compare different optimizations of the same design with the best design having the largest efficiency.

The design algorithm fits the blocks design by sequentially adding the block factors in the column order of the blocks data frame. Each block factor is optimized conditional on all preceding block factors remaining fixed but ignoring all succeeding block factors. This method allows the blocking factors to be fitted in order of importance with the largest and most important blocks fitted first and the smaller and less important blocks fitted subsequently.

For crossed blocks designs, a differential weighting factor w is used to determine the relative importance of the block main effects versus the block interaction effects. If w = 0 the algorithm fits a simple additive main effects design whereas if w = 1 the algorithm fits a fully crossed blocks design. For intermediate 0 < w < 1, the block factor interaction effects are downweighted relative to the main effects, the smaller the value of w the greater the downweighting. The default weighting is 0.5 and provided that all block effects are estimable, this weighting gives a compromise design where all factorial block effects are included in the blocks design but where the main block effects are given greater importance than the block interaction effects.See vignette(package = "blocksdesign") for more details.

The blocks design criterion is the D-efficiency of the treatment design without block effects versus the efficiency of the treatment design with block effects. The D-efficiency factor is necessarily less than or equal to one with the most efficient design giving the largest D-efficiency factor. For unstructured treatment designs, the A-efficiency factor is also shown together with an estimated A-efficiency upper-bound, where available.

For more details see vignette(package = "blocksdesign")

References

Cochran W. G. & Cox G. M. (1957) Experimental Designs 2nd Edition John Wiley & Sons.

Examples

# NOT RUN {
## For optimum results, the number of searches may need to be increased.

## 4 replicates of 12 treatments with 16 nested blocks of size 3
# rectangular lattice see Plan 10.10 Cochran and Cox 1957.
treatments = factor(1:12)
blocks = data.frame(Main = gl(4,12), Sub = gl(16,3))
# }
# NOT RUN {
design(treatments,blocks)
# }
# NOT RUN {
## 4 x 12 design for 4 replicates of 12 treatments with 3 plots in each intersection block
## The optimal design is Trojan with known A-efficiency = 22/31 for the intersection blocks
treatments = factor(1:12)
blocks = data.frame(Rows = gl(4,12), Cols = gl(4,3,48))
# }
# NOT RUN {
design(treatments,blocks)
# }
# NOT RUN {
## 4 x 12 design for 4 replicates of 12 treatments with 3 sub-column blocks nested 
## as above but showing 3 sub-columns nested within each main column
treatments = factor(1:12)
blocks = data.frame(Rows = gl(4,12), Cols = gl(4,3,48), subCols = gl(12,1,48))
# }
# NOT RUN {
design(treatments,blocks,searches=200)
# }
# NOT RUN {
## 4 x 13 Row-and-column design for 4 replicates of 13 treatments 
## Youden design Plan 13.5 Cochran and Cox (1957).
treatments = factor(1:13)
blocks = data.frame(Rows = gl(4,13), Cols = gl(13,1,52))
# }
# NOT RUN {
design(treatments,blocks,searches = 700)
# }
# NOT RUN {
## differential replication 
treatments=factor(c(rep(1:12,2),rep(13,12)))
blocks = data.frame(Main = gl(2,18),  Sub = gl(12,3,36))
design(treatments,blocks,searches = 5)

## 48 treatments in 2 replicate blocks of size 48 for a 24 x 4 array 
## with 2 main rows and 3 main columns the cols factor must precede 
## the rows factor otherwise the design will confound one treatment contrast
## in the replicates.rows x columns interactions due to inherent aliasing 
treatments=factor(1:48)
blocks = data.frame(Reps = gl(2,48),Cols = gl(3,8,96),Rows = gl(2,24,96))
design(treatments,blocks,searches=5)


## Factorial treatment designs defined by a single factorial treatment model

## Main effects of five 2-level factors in a half-fraction in 2/2/2 nested blocks design 
## (may require many 100's of repeats to find a fully orthogonal solution)
treatments = expand.grid(F1 = factor(1:2), F2 = factor(1:2),
 F3 = factor(1:2), F4 = factor(1:2), F5 = factor(1:2))
blocks = data.frame(b1 = gl(2,8),b2 = gl(4,4),b3 = gl(8,2))
model=" ~ F1 + F2 + F3 + F4 + F5"
# }
# NOT RUN {
repeat {z = design(treatments,blocks,treatments_model=model,searches=50)
if ( isTRUE(all.equal(z$blocks_model[3,3],1) ) ) break }
 print(z)
# }
# NOT RUN {
 
# Second-order model for five qualitative 2-level factors in 4 randomized blocks
treatments = expand.grid(F1 = factor(1:2), F2 = factor(1:2), F3 = factor(1:2), 
F4 = factor(1:2), F5 = factor(1:2))
blocks = data.frame(blocks = gl(4,8))
model = " ~ (F1 + F2 + F3 + F4 + F5)^2"
design(treatments,blocks,treatments_model=model,searches = 10)

# Main effects of five 2-level factors in a half-fraction of 
# a 4 x 4 row-and column design.
treatments = expand.grid(F1 = factor(1:2), F2 = factor(1:2), F3 = factor(1:2), 
F4 = factor(1:2), F5 = factor(1:2))
blocks = data.frame(rows = gl(4,4), cols = gl(4,1,16))
model = "~ F1 + F2 + F3 + F4 + F5"
# }
# NOT RUN {
repeat {z = design(treatments,blocks,treatments_model=model,searches=50)
if ( isTRUE(all.equal(z$blocks_model[2,3],1) ) ) break }
 print(z)
# }
# NOT RUN {
# Quadratic regression for three 3-level numeric factor assuming a 10/27 fraction
treatments = expand.grid(A = 1:3, B = 1:3, C = 1:3)
blocks=data.frame(main=gl(1,10))
model = " ~ ( A + B + C)^2 + I(A^2) + I(B^2) + I(C^2)"
design(treatments,blocks,treatments_model=model,searches=10) 

# Quadratic regression for three 3-level numeric factor crossed with a qualitative 2-level factor
treatments = expand.grid(F = factor(1:2), A = 1:3, B = 1:3, C = 1:3)
blocks=data.frame(main=gl(1,18))
model = " ~ F + A + B + C + F:A + F:B + F:C + A:B + A:C + B:C + I(A^2) + I(B^2) + I(C^2)"
design(treatments,blocks,treatments_model=model,searches=5) 

# 1st-order model for 1/3rd fraction of four qualitative 3-level factors in 3  blocks
treatments = expand.grid(F1 = factor(1:3), F2 = factor(1:3), F3 = factor(1:3), 
F4 = factor(1:3))
blocks = data.frame(main = gl(3,9))
model = " ~ F1 + F2 + F3 + F4"
# }
# NOT RUN {
 design(treatments,blocks,treatments_model=model,searches=25)
# }
# NOT RUN {
# 2nd-order model for a 1/3rd fraction of five qualitative 3-level factors in 3 blocks
# (may require many repeats to find a fully orthogonal solution)
treatments = expand.grid(F1 = factor(1:3), F2 = factor(1:3), F3 = factor(1:3), 
F4 = factor(1:3), F5 = factor(1:3))
blocks=data.frame(main=gl(3,27))
model = " ~ (F1 + F2 + F3 + F4 + F5)^2"
# }
# NOT RUN {
 repeat {z = design(treatments,blocks,treatments_model=model,searches=50)
if ( isTRUE(all.equal(z$blocks_model[1,3],1) ) ) break}
 print(z) 
# }
# NOT RUN {
# 2nd-order model for two qualitative and two quantitative level factors in 2 blocks of size 18
treatments = expand.grid(F1 = factor(1:2), F2 = factor(1:3), V1 = 1:3, V2 = 1:4)
blocks = data.frame(main = gl(2,18))
model = " ~ (F1 + F2 + V1 + V2)^2 +  I(V1^2) +  I(V2^2)"
# }
# NOT RUN {
 design(treatments,blocks,treatments_model=model,searches=5)
# }
# NOT RUN {
 
# Plackett and Burman design for eleven 2-level factors in 12 runs 
GF = expand.grid(F1 = factor(1:2,labels=c("a","b")), F2 = factor(1:2,labels=c("a","b")), 
                 F3 = factor(1:2,labels=c("a","b")), F4 = factor(1:2,labels=c("a","b")),
                 F5 = factor(1:2,labels=c("a","b")), F6 = factor(1:2,labels=c("a","b")),
                 F7 = factor(1:2,labels=c("a","b")), F8 = factor(1:2,labels=c("a","b")), 
                 F9 = factor(1:2,labels=c("a","b")), F10= factor(1:2,labels=c("a","b")), 
                 F11= factor(1:2,labels=c("a","b")) )
blocks=data.frame(main=gl(1,12))
model = "~ F1 + F2 + F3 + F4 + F5 + F6 + F7 + F8 + F9 + F10 + F11"
# }
# NOT RUN {
design(GF,blocks,treatments_model=model,searches=25)
# }
# NOT RUN {
## Factorial treatment designs defined by sequentially fitted factorial treatment models

## 2 varieties x 3 levels of N x 3 levels of K assuming 1st-order interactions and 12 plots
## the single stage model gives an unequal 7 + 5 split for the two varieties
## whereas the two stage model forces an equal 6 + 6 split
## NB the two stage model is slightly less efficient than the single stage model 
treatments = expand.grid(Variety = factor(rep(1:2)), N = 1:3, K = 1:3)
blocks=data.frame(main=gl(1,12))
treatments_model = " ~  (Variety + N + K)^2  + I(N^2) + I(K^2)"
# }
# NOT RUN {
design(treatments,blocks,treatments_model=treatments_model,searches=10)
# }
# NOT RUN {
treatments_model = " ~ Variety | (Variety + N + K)^2 + I(N^2) + I(K^2)"
# }
# NOT RUN {
design(treatments,blocks,treatments_model=treatments_model,searches=10)
# }
# NOT RUN {
## 8 varieties x 4 levels of N x 4 levels of K assuming 1st-order interactions and 2 blocks
## of size 16, A two stage treatment model gives fully orthogonal blocks with
## no loss of treatment efficiency comapred with a single stage treatment design
treatments = expand.grid(variety = factor(1:8), N = 1:4, K = 1:4)
blocks=data.frame(blocks=gl(2,16,32))
treatments_model = "~ (variety + N + K)^2  + I(N^2) + I(K^2) + I(N^3) + I(K^3)"
# }
# NOT RUN {
design(treatments,blocks,treatments_model=treatments_model,searches=50)
# }
# NOT RUN {
treatments_model = "~variety | (variety + N + K)^2   + I(N^2) + I(N^3) + I(K^2)+ I(K^3)"
# }
# NOT RUN {
design(treatments,blocks,treatments_model=treatments_model,searches=50)
# }
# NOT RUN {
## A 6 x 6 row-and-column design with linear row by linear column interaction.
## Crossed blocks with interactions fitted in the treatments model and additive 
## treatments fitted inthe blocks model as a dual design
##  see vignette(package = "blocksdesign") for further discussion
## may require many separate attempts to get the best overall design efficiency 
LS_grid   = expand.grid(rows=factor(1:6), cols=factor(1:6))
blocks = data.frame(varieties=factor(rep(1:6,6)))
lin_rows = as.numeric(levels(LS_grid$rows))[LS_grid$rows]
lin_cols = as.numeric(levels(LS_grid$cols))[LS_grid$cols]
latin_sq = "~  rows | cols + lin_rows:lin_cols "
# }
# NOT RUN {
design(LS_grid,blocks,latin_sq,searches=2000)
# }
# NOT RUN {
# }