#' Stuttgart on 26 March 2017
#' Version 1.1
#'
#' create a function to calculate A-efficency factors
#' AOpt <- function(Design,nTreatments,nRows,nCols,repli)
#' input : flag = 0: default
#' 1: A++
#' 2: A22 for unreplicated treatments
#' ncheck: this is a number: 1,2,..ncheck,...,ntreatments
#' 3: A11 for checks
#' output:
#'
#' A function to calculate an average efficent factor for a block design including one blocking factor,
#' a block is identified by columns
#'
#' @param Design
#' @param nTreatments
#' @param nRows
#' @param nCols
#' @param rowgroup
#' @param colgroup
#' @param subrgroup
#' @param subcgroup
#' @param repli
#' @param ncheck
#' @param flag {0, 1, 2, 3}
#' @param cRandom=0
#' @param reps=1
#' @return an efficiency of a given design
#' @export
A3BlkDesign<- function(pDesign, nTreatments, nRows, nCols, rowgroup, colgroup, subrgroup, subcgroup, repli, ncheck, flag, cRandom=0, reps=1){
# get a model matrix
modelD <- model3Matrix(nTreatments, nRows, nCols, rowgroup, colgroup, subrgroup, subcgroup, pDesign)
# get an information matrix for a given design
Ar <- modelD$X
onesM <- matrix(1,nRows*nCols,1)
Br <- cbind(onesM,modelD$R,modelD$C,modelD$modelB)
tAr <- t(Ar)
tBr <- t(Br)
rdelta <- tAr %*% Ar
r <- tAr %*% onesM
r_a <- sum(r)/nTreatments
r_square <- r_a*r_a
#rdeltapower <- rdelta^(1/2)
#rdeltapowerinv <- solve(rdeltapower)
Cr <- tAr %*% Ar -tAr %*% Br %*%ginv(tBr %*%Br) %*% tBr %*% Ar
## need to check
#Cr_1 <- rdeltapowerinv %*% Cr %*% rdeltapowerinv
raCr <- rankMatrix(Cr)[1]
print("ranking system")
#print(Cr)
if (flag=="0"){
if(raCr==nTreatments-1){
invCr <- ginv(Cr)
#Aopt <- 1000 # this is the default value, which proofs the connected designs
temp <- 0
for (ii in c(1:(nTreatments-1))){
for(jj in c((ii+1):nTreatments)){
v <- (r[ii]*r[jj]/r_square)*(invCr[ii,ii]+invCr[jj,jj]-2*invCr[ii,jj])
temp <- temp + v
}
}
v_bar <- temp*2/(nTreatments*(nTreatments-1))
Aopt <- 2/(r_a*v_bar)
} else{
Aopt <- 0
}
}else if(flag=="1"){
invCr <- ginv(Cr)
#Aopt <- 1000 # this is the default value, which proofs the connected designs
temp <- 0
for (ii in c(1:(nTreatments-1))){
for(jj in c((ii+1):nTreatments)){
v <- (r[ii]*r[jj]/r_square)*(invCr[ii,ii]+invCr[jj,jj]-2*invCr[ii,jj])
temp <- temp + v
}
}
v_bar <- temp*2/(nTreatments*(nTreatments-1))
Aopt <- 2/(r_a*v_bar)
#a_plus <- 2*(sum(diag(invCr))-(sum(invCr)/nTreatments))/(nTreatments-1)
#Aopt <- a_plus
}else if (flag=="u"){
# insert your code here
invCr <- ginv(Cr)
# get an information matrix for unreplicated treatments
if (cRandom==0){
inforUT <- invCr[(ncheck+1):nTreatments,(ncheck+1):nTreatments]
}
else{
ftrs <- fChecks(nTreatments,Design,reps)
untrts <-ftrs$untrts
inforUT <- invCr[untrts,untrts]
}
v_bar <- 2*(sum(diag(inforUT))-(sum(inforUT))/(nTreatments-ncheck))/(nTreatments-ncheck-1)
Aopt <- 2/v_bar
#a_plus <- 2*(sum(diag(inforUT))-(sum(inforUT))/(nTreatments-ncheck))/(nTreatments-ncheck-1)
#Aopt <- a_plus
}else {
invCr <- ginv(Cr)
# get an information matrix for check treatments
checkRep <- rep(0,ncheck)
if (cRandom==0){
inforCT <- invCr[1:ncheck,1:ncheck]
checkRep <- r[c(1:ncheck)]
}
else{
ftrs <- fChecks(nTreatments,Design,reps)
checks <- ftrs$checks
Freqs <- ftrs$Freqs
checkRep <- Freqs[checks]
inforCT <- invCr[checks,checks]
}
temp <- 0
r_c <- (sum(r)-(nTreatments-ncheck))/ncheck
rc_square <- r_c*r_c
for (ii in c(1:(ncheck-1))){
for(jj in c((ii+1):ncheck)){
v <- (checkRep[ii]*checkRep[jj]/rc_square)*(inforCT[ii,ii]+inforCT[jj,jj]-2*inforCT[ii,jj])
temp <- temp + v
}
}
v_bar <- temp*2/(ncheck*(ncheck-1))
Aopt <- 2/(r_c*v_bar)
#a_plus <- 2*(sum(diag(inforCT))-(sum(inforCT))/(ncheck))/(ncheck-1)
#Aopt <- a_plus
}
}
#' create a model matrix
#' adds features on 14 March 2017
#' @author Nha Vo-Thanh
#' @description Augmented block designs
#' @param nTreatments
#' @param nRows
#' @param nCols
#' @param rowgroup
#' @param colgroup
#' @param subrgroup=0
#' @param subcgroup=0
#' @param A a design matrix
#' @return a model matrix
#' @export
#require(Matrix)
#library(MASS)
model3Matrix <- function(nTreatments, nRows, nCols, rowgroup, colgroup, subrgroup=0, subcgroup=0 ,A , flag=1)
{
#create an matrix for an additional blocking factor - 14 March 2017
modelB <- matrix(0, nRows * nCols, rowgroup* colgroup)
rInGroup <- nRows/rowgroup
cInGroup <- nCols/colgroup
trtModel <- matrix(0, nRows * nCols, nTreatments)
for (i in 1:nRows) {
for (j in 1:nCols) {
# set level for treatments in a model matrix
for (trt in 1:nTreatments) {
if (A[i, j] == trt) {
trtModel[(i - 1) * nCols + j, trt] = 1
}
}
}
}
## a model for rows
init <- rep(1,nCols)
mlist <- rep(list(init),nRows)
rModel <- as.matrix(bdiag(mlist))
## a model for columns
initCol <- rep(1,nCols)
initRow <- rep(1,nRows)
initblock <- diag(initCol)
cModel <- kronecker(initRow, initblock)
# Building a model for augmented designs, including 3 blocking factors
print("checking here")
for(i in 1:rowgroup){
for (j in 1:colgroup){
for (k in ((i-1)*rInGroup+1):(i*rInGroup)){
for (l in ((j-1)*cInGroup+1):(j*cInGroup)){
modelB[(k-1)*nCols+l,(i-1)*colgroup+j] <- 1
}
}
}
}
# in case when unequal subrows and subcolumns occur
if(flag==1){
temp1 <- 0
for(i in 1:rowgroup){
temp2 <- list()
for(j in 1:colgroup){
init <- rep(1,subcgroup[j])
temp2 <- append(temp2,list(init))
}
z <- rep(1,subrgroup[i])
block1 <- as.matrix(bdiag(temp2))
final <- kronecker(z,block1)
temp1 <- bdiag(temp1,final)
}
# a final matrix including zero column and row
modelC <- as.matrix(temp1)
modelC <- modelC[2:dim(modelC)[1],2:dim(modelC)[2]]
}else{
modelC <- 0;
}
model <- list(trtModel,rModel,cModel,modelB, modelC)
names(model) <- c("X","R","C","modelB","modelC")
return(model)
}
##################################
# Searching optimal row-column designs using simulated annealing algorithm
# input :
# output:
# Version 1.1
# Stuttgart 11 April 2017
##################################
#' @param nTreatments the number of treatments
#' @param nCheks the number of checks
#' @param nRows the number of rows
#' @param nCols the number of columns
#' @return return an optimal design
#' @export
start3DesignV1 <- function(nTreatments, nChecks, nRows, nCols, rowgroup, colgroup,repli, subrgroup=0, subcgroup=0 ){
# Construct an incidence matrix N = X' %*% R
# Starting with a design
startD <- matrix(0, nRows, nCols)
# Allocate the different treatments to one block at a time
nObs <- nRows*nCols
newObs <- nObs
N <- matrix(0,nTreatments,nRows)
newrepli <-repli
listtrt <- 1:nTreatments
for(ib in 1:nRows){
kn <- nCols/newObs
for (i in 1:nTreatments){
N[i,ib] <- floor(newrepli[i]*kn)
}
rfill <- nCols-sum(N[,ib])
apha <- newrepli-N[,ib]
ixx <- listtrt[apha>0]
selTrts <- sample(ixx, rfill, replace=F)
N[selTrts,ib] <- N[selTrts,ib]+1
newObs <- newObs-nCols
newrepli <- newrepli-N[,ib]
}
##--- Constructing a design from a given incidence matrix
# N = t(X) %*% Z
startD <- matrix(0, nRows, nCols)
for (ir in 1:nRows){
startD[ir,] <- rep(listtrt,N[,ir])
}
return(startD)
}
#' The function helps to generate a starting connected design
#' @param nTreatments the number of treatments
#' @param nCheks the number of checks
#' @param nRows the number of rows
#' @param nCols the number of columns
#' @return return an optimal design
#' @export
start3DesignV2<- function(nTreatments, nChecks, nRows, nCols, repli, nchecks, vrepChecks, rowgroup, colgroup, subrgroup=0, subcgroup=0){
# Construct an incidence matrix N = X' %*% R
# Starting with a design
#--------------------------------------------------------#
nTrts <- nTreatments-nChecks+1
repliz <- rep(1, nTrts)
repliz[1] <- sum(vrepChecks)
totalchecks <- sum(vrepChecks)
fullchecks <- rep(c(1:nChecks),vrepChecks)
optDesign <- 0
checksM <- 0
while((optDesign<=0)&&(checksM<=101)){
checksM <- 0
while(optDesign<=0){
pDesign <- startDesignV1(nTrts, nChecks, nRows, nCols, repliz)
for(i in 1:nRows){
indexx <- sample(1:nCols, nCols, replace = F)
pDesign[i,] <-pDesign[i,indexx]
}
optDesign <- aRCDesign(pDesign, nTrts, nRows, nCols, repliz, ncheck=1, flag=0, cRandom = 1)
}
nDesign <- pDesign
nDesign[which(nDesign==1)] <- rep(-3,sum(vrepChecks))
nDesign <- nDesign + (nChecks-1)
xxxx <- (nChecks-1) + -3
repchecks <- fullchecks
indexx <- sample(1:totalchecks, totalchecks, replace = F)
repchecks <- repchecks[indexx]
nDesign[which(nDesign==xxxx)] <- repchecks
optDesign <- A3BlkDesign( pDesign, nTrts, nRows, nCols, rowgroup, colgroup, subrgroup, subcgroup, repliz, nreCheck, flag=0, cRandom=0, reps=1)
# randomly assign checks to qseudo check
while((optDesign<=0)&&(checksM<=100)){
checksM <- checksM +1
repchecks <- fullchecks
indexx <- sample(1:totalchecks, totalchecks, replace = F)
repchecks <- repchecks[indexx]
nDesign[which(nDesign==-1)] <- repchecks
optDesign <- A3BlkDesign( pDesign, nTrts, nRows, nCols, rowgroup, colgroup, subrgroup, subcgroup, repliz, nreCheck, flag=0, cRandom=0, reps=1)
}
#print(optDesign)
#print(checksM)
}
mlist <- list(nDesign,optDesign)
names(mlist) <- c("nDesign","optDesign")
return(mlist)
}
##----------------------------------##
## A place to test your codes
##----------------------------------##
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.