R: General block and treatment designs.

design {blocksdesign}

R Documentation

General block and treatment designs.

Description

Constructs D-optimal block and treatment designs for any feasible linear treatment model and any feasible combination of block factors.

Usage

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

Arguments

`treatments`	a candidate set of treatments for any combination of qualitative or quantitative level factors.
`blocks`	a data frame of block factors.
`treatments_model`	a list containing one or more nested treatment model formula.
`weighting`	a weighting factor between 0 and 1 for weighting the 2-factor interaction effects of factorial blocks.
`searches`	the maximum number of local optima searched at each stage of an optimization.
`seed`	an integer initializing the random number generator.
`jumps`	the number of pairwise random treatment swaps used to escape a local maxima.

Details

treatments a factor or data frame or list of generators for the candidate set of treatments. The candidate treatment set can include both quantitative and qualitative level factors and the required design is selected from the candidate treatment set without replacement. The number of times a treatment can be included in a design depends on the replication of that treatment in the candidate set therefore increasing the treatment replication in the candidate set allows increased replication for individual treatments.

blocks a data frame of nested or crossed block factors. The length of the block factors defines the required size of the treatment design. the treatment design is optimized for the required blocks design by maximising the determinant of the block-adjusted treatment information matrix (D-optimality. For crossed block designs, the information matrix is based on a weighted combination of block main effects and two-factor interaction effects where the information on the interaction effects is down-weighted by a coefficient 0 < w < 1 where w = 0 gives blocks main effects only whereas w = 1 gives the full two-factor blocks interaction model. This process ensures that factorial main effects can be given more importance than the interaction effects in the optimization of a crossed blocks design.

treatments_model a character vector containing one or more treatments model formula. The treatment models are optimized sequentially for each model formula in turn assuming the treatments of any previously fitted models are held constant. Sequential fitting provides improved flexibility for fitting treatment factors or variables of different status or importance (see examples below). The design criterion for each treatment model is maximization of the determinant of the information matrix of that treatment model conditional on the the information matrices of all previously fitted treatment models.

The D-optimal treatment design efficiency is the ratio of the generalized variance of the full candidate treatment model relative to the generalized variance of the optimized design. The efficiency will be less than or equal to 1 for factorial models but may exceed 1 for polynomial models.

The design outputs include a table of efficiency factors for first and second-order factorial block effects and comparison of the efficiency factors for different choices of w can be used to find a good compromise design for fitting both main block and interaction block effects.

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

Value

`Replication`	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 in standardized block order.
`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.

References

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

Gupta, S. C., and B. Jones. Equireplicate Balanced Block Designs with Unequal Block Sizes. Biometrika, vol. 70, no. 2, 1983, pp. 433–440.

Examples

## For best results, the number of searches may need to be increased.

## 4 replicates of 12 treatments with 16 nested blocks of size 3
## giving a rectangular lattice: see Plan 10.10 Cochran and Cox 1957.

blocks = data.frame(Main = gl(4,12), Sub = gl(16,3))
treatments = data.frame(treatments = gl(12,1,48))
Z=design(treatments, blocks)
print(Z)

## 3 replicates of 15 treatments in 3 main blocks with two sets of
## nested blocks and one control treatment

blocks=data.frame( Main = gl(3,18,54),Sub1 = gl(9,6,54),Sub2 = gl(27,2,54))
treatments=factor(rep(c(1:15,rep("control",3)),3))
Z=design(treatments,blocks)
print(Z)
incid1=table(interaction(Z$Design$Main,Z$Design$Sub1,lex.order = TRUE),Z$Design$treatments)
crossprod(incid1) # print pairwise concurrences within Sub1 blocks
incid2=table(interaction(Z$Design$Main,Z$Design$Sub2,lex.order = TRUE),Z$Design$treatments)
crossprod(incid2) # print pairwise concurrences within Sub2 blocks

 
## 4 x 12 design for 4 replicates of 12 treatments with 3 plots in each intersection block
## and  3 sub-columns nested within each main column. The optimal row-and column design
## is Trojan with A-efficiency = 22/31 for the intersection blocks
 
blocks = data.frame(Rows = gl(4,12), Cols = gl(4,3,48), subCols = gl(12,1,48))
treatments = data.frame(treatments =gl(12,1,48))
Z=design(treatments,blocks)
print(Z)
incid=table(interaction(Z$Design$Rows,Z$Design$Cols,lex.order = TRUE),Z$Design$treatments)
crossprod(incid) # print pairwise concurrences within blocks


## 4 x 13 Row-and-column design for 4 replicates of 13 treatments 
## Youden design Plan 13.5 Cochran and Cox (1957).

blocks = data.frame(Rows = gl(4,13), Cols = gl(13,1,52))
treatments = data.frame(treatments= gl(13,1,52))
Z=design(treatments,blocks,searches = 700)
print(Z)
incid=table(Z$Design$Cols,Z$Design$treatments)
crossprod(incid) # print pairwise concurrences of treatments within column blocks (BIB's)


## 48 treatments in 2 replicate blocks with 2 nested rows in each replicate and 3 main columns

blocks = data.frame(Reps = gl(2,48), Rows = gl(4,24,96), Cols = gl(3,8,96))
treatments = data.frame(treatments=gl(48,1,96))
Z=design(treatments,blocks,searches=5)
print(Z)

## 48 treatments in 2 replicate blocks with 2 main columns
## The default weighting gives non-estimable Reps:Cols effects due to inherent aliasing
## Increased weighting gives estimable Reps:Cols effects but non-orthogonal main effects

blocks = data.frame(Rows = gl(2,48), Cols = gl(2,24,96))
Z1=design(treatments=gl(48,1,96),blocks,searches=5)
print(Z1)
Z2=design(treatments=gl(48,1,96),blocks,searches=5,weighting=.9)
print(Z2)


## Equireplicate balanced designs with unequal block sizes. Gupta & Jones (1983). 
## Check for equality of D and A-efficiency in optimized design

t=factor(c(rep(1:12,each=7)))
b=factor(c(rep(1:12,each=6),rep(13:18,each=2)))
Z=design(t,b,searches=100)$Blocks_model # max efficiency = 6/7
print(Z)


t=factor(c(rep(1:14,each=8)))
b=factor(c(rep(1:14,each=4),rep(15:21,each=8)))
Z=design(t,b,searches=500)$Blocks_model # max efficiency = 7/8
print(Z)


t=factor(c(rep(1:16,each=7)))
b=factor(c(rep(1:16,each=4),rep(17:22,each=8)))
Z=design(t,b,searches=1000)$Blocks_model # max efficiency = 6/7
print(Z)


t=factor(c(rep(1:18,each=7)))
b=factor(c(rep(1:18,each=6),rep(19:24,each=3)))
Z=design(t,b,searches=500)$Blocks_model # max efficiency = 6/7
print(Z)


## Factorial treatment designs defined by a factorial treatment model

## Main effects of five 2-level factors in a half-fraction in 2/2/2 nested blocks design 
## The algorithmic search method is likely to require a very long search time to reach the
## known orthogonal solution for this regular fractional factorial design 

treatments = list(F1 = gl(2,1), F2 = gl(2,1),F3 = gl(2,1), F4 = gl(2,1), F5 = gl(2,1))
blocks = data.frame(b1 = gl(2,8),b2 = gl(4,4),b3 = gl(8,2))
model= ~ F1 + F2 + F3 + F4 + F5
repeat {Z = design(treatments,blocks,treatments_model=model,searches=50)
if ( isTRUE(all.equal(Z$Blocks_model[3,3],1) ) ) break }
print(Z)

 
# Second-order model for five qualitative 2-level factors in 4 randomized blocks with two
# nested sub-blocks in each main plot
treatments = list(F1 = gl(2,1), F2 = gl(2,1),F3 = gl(2,1), F4 = gl(2,1), F5 = gl(2,1))
blocks = data.frame(main = gl(4,8),sub=gl(8,4))
model = ~ (F1 + F2 + F3 + F4 + F5)^2
Z=design(treatments,blocks,treatments_model=model,searches = 10)
print(Z)

# Main effects of five 2-level factors in a half-fraction of a 4 x 4 row-and column design

treatments = list(F1 = gl(2,1), F2 = gl(2,1),F3 = gl(2,1), F4 = gl(2,1), F5 = gl(2,1))
blocks = data.frame(rows = gl(4,4), cols = gl(4,1,16))
model = ~ F1 + F2 + F3 + F4 + F5
repeat {Z = design(treatments,blocks,treatments_model=model,searches=50)
if ( isTRUE(all.equal(Z$Blocks_model[2,3],1) ) ) break}
print(Z)


# Quadratic regression for three 3-level numeric factor assuming a 10/27 fraction
treatments = list(A = 1:3, B = 1:3, C = 1:3)
blocks=factor(rep(1,10))
model = ~ ( A + B + C)^2 + I(A^2) + I(B^2) + I(C^2)
Z=design(treatments,blocks,treatments_model=model,searches=10) 
print(Z)

## response surface designs with doubled treatment candidate sets 
## allowing duplicated design points

## two treatment factor design showing double replication on the 
## four corner points and double replication on the centre point
treatments = list( V1 = 1:3, V2 = rep(1:3,2))
blocks=factor(rep(1,14))
model =  ~ (V1 + V2)^2 + I(V1^2) + I(V2^2)
design(treatments,blocks,treatments_model=model,searches=10)

## three treatment factor design with eight corner points, 12 edge points and 1 centre 
## point. All 21 design points are accounted for leaving nothing for extra replication.
treatments = list( V1 = 1:3, V2 = 1:3, V3 = rep(1:3,2))
blocks=factor(rep(1,21))
model =  ~ (V1 + V2 + V3)^2 + I(V1^2) + I(V2^2) + I(V3^2)
design(treatments,blocks,treatments_model=model,searches=20)

## as above but with extra replication on each corner point and on the centre point
treatments = list( V1 = 1:3, V2 = 1:3, V3 = rep(1:3,2))
blocks=factor(rep(1,30))
model =  ~ (V1 + V2 + V3)^2 + I(V1^2) + I(V2^2) + I(V3^2)
design(treatments,blocks,treatments_model=model,searches=20)

## Factorial treatment designs defined by sequentially fitted factorial treatment models

## 4 varieties by 3 levels of N by 3 levels of K assuming degree-2 treatment
## interaction effects and two blocks of 12 plots
## the single stage model gives an unequal split for the replication of the four varieties
## which may be undesirable whereas the two stage model forces an equal split of 6 plots
## per variety. The single stage model is slightly more efficient than the two stage model
## (1.052045 versus 0.9761135 standardized relative to the full factorial design)  but
## in this example global D-optimality does not give the most suitable design structure. 

treatments = list(Variety = factor(1:4), N = 1:3, K = 1:3)
blocks = data.frame(main=gl(2,12))
treatments_model = ~ (Variety + N + K)^2  + I(N^2) + I(K^2)
Z1 = design(treatments,blocks,treatments_model=treatments_model,searches=10)
print(Z1)
treatments_model = list( ~ Variety , ~ Variety + (Variety + N + K)^2 + I(N^2) + I(K^2))
Z2 = design(treatments,blocks,treatments_model=treatments_model,searches=10)
print(Z2)

[Package blocksdesign version 4.9 Index]