林木子代测定林田间试验较常采用随机完全区组设计(RCB)和不完全区组设计(IB)。在R
中利用agricolae
包能够实现RCB设计,对于后者则只能进行平衡不完全区组设计(BIB)。BIB要求bkvrλ全为整数,b是区组数、k是区组大小即小区数、v是处理数比如家系个数、r是重复次数,λ=r(k-1)/(v-1),对于测定处理较多的材料不能满足λ为整数这个限制。ibd
这个包开发者说能够实现不平衡IB设计,但同样在处理较多时给不出设计结果。最终在网页找到方法,函数的代码如下
incomplete.block.design <- function(number.alternatives, number.blocks, alternatives.per.block, n.repeats = 1000){
# Check that the parameters are appropriate
# Sawtooth recommends that number.blocks >= 3 * number.alternatives / alternatives.per.block
if (number.blocks < 3 * number.alternatives / alternatives.per.block)
warning("It is recomended that number.blocks >= 3 * number.alternatives / alternatives.per.block");
library(AlgDesign)
best.result = NULL
best.D = -Inf
for (i in 1:n.repeats){
alg.results <- optBlock(~.,withinData=factor(1:number.alternatives),blocksizes=rep(alternatives.per.block,number.blocks), nRepeats=5000) #BIB
if (alg.results$D > best.D){
best.result = alg.results
best.D = alg.results$D
}
}
design <- matrix(NA,number.blocks,alternatives.per.block, dimnames= list(block = 1:number.blocks, Alternative = 1:alternatives.per.block))
binary.design <- matrix(0,number.blocks,number.alternatives, dimnames= list(block = 1:number.blocks, alternative = 1:number.alternatives))
counter <- 0
for (block in best.result$Blocks){
counter <- counter + 1
blck <- unlist(block)
design[counter,] <- blck
for (a in blck)
binary.design[counter,a] <- 1
}
n.appearances.per.alternative <- table(as.numeric(design))
combinations.of.alternatives <- crossprod(table(c(rep(1:number.blocks, rep(alternatives.per.block,number.blocks))), best.result$design[,1]))
list(binary.design = t(binary.design), design = t(design), frequencies = n.appearances.per.alternative, pairwise.frequencies=combinations.of.alternatives, binary.correlations = round(cor(binary.design),2))
}
套入White T L, Adams W T, Neale D B. 2007. Forest genetics. 1st ed. CABI, Cambridge (MA).384页的例子
incomplete.block.design(number.alternatives = 24, number.blocks = 18, alternatives.per.block = 4, n.repeats = 100)
重复设定为100次后得到结果
与书中所使用的软件给出的结果不一致,但应该都能用于实践(现实情况是我们过去的测定林设计多少存在些问题(包括我自己之前的一个设计),或者设计问题不大但造林时没有严格按照设计进行,比较可惜)。
林木周期长,测定很关键,设计要专业