R/Block2Designs.R
In nhavt/Block3AugDesigns:
Defines functions
model2Matrix
model2MatrixV2
info2Matrix
delInfoMatrix
addInfoMatrix
startDesignV1
startDesignV2
genBlkDesign
updateFormulaeT
Swapping
updateFormulaeT3
layoutAugDesign
NeighboorUpdate
CoordinateExchange
CompareSwappingDesign
SwappingTwoPlots
RandomSearch
AeffRowColumnDesign
#' Develop by Nha Vo-Thanh
#' Stuttgart on 22 March 2017
#' Augmented row-column designs
#' Version 1.0
#'
##################################
#' create a model matrix
#' @param nTreatments number of treatments
#' @param nRows number of rows
#' @param nCols number of columns
#' @return a model matrix
#' @export
model2Matrix <- function(nTreatments, nRows, nCols, A)
{
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)
model <- list(trtModel,rModel,cModel)
names(model) <- c("X","R","C")
return(model)
}
##################################
#' create a model matrix
#' @param nTreatments number of treatments
#' @param nRows number of rows
#' @param nCols number of columns
#' @return a model matrix
#' @export
model2MatrixV2 <- function(nTreatments, nRows, nCols, A)
{
trtModel <- matrix(0, nRows * nCols, nTreatments)
vecA <- c(t(A))
mylist <- rep(list(vecA),nTreatments)
fxn <- function(var1,var2){
as.numeric(var1==var2)
}
var2 <- as.list(c(1:nTreatments))
trtModel <- mapply(fxn,mylist,var2)
## 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)
model <- list(trtModel,rModel,cModel)
names(model) <- c("X","R","C")
return(model)
}
##################################
# return an information matrix
# input :
# output:
# Version 1.1
##################################
#' @param Design a given design
#' @param nTreatments the number of treatments
#' @param nRows the number of rows
#' @param nCols the number of columns
#' @param flag a flag to decide whether to include inverse matrix
#' @return information matrix
#' @export
info2Matrix <- function(Design,nTreatments, nRows, nCols, flag=0){
# connectedFlag
connectedFlag <- 0
#modelA <- model2Matrix(nTreatments, nRows, nCols, Design)
modelA <- model2MatrixV2(nTreatments, nRows, nCols, Design)
# information for treatments
X <- modelA$X
# information for rows and columns
R <- modelA$R
C <- modelA$C
# tranpose matrix
tX <- t(X)
tR <- t(R)
tC <- t(C)
# incidence matrix
NR <- tX %*% R
NC <- tX %*% C
# tranpose matrix
tNR <- t(NR)
tNC <- t(NC)
# creates a vector of ones
onevec <- matrix(1, nrow = dim(tX)[2], ncol = 1)
r <- tX %*% onevec
tr <- t(r)
rdelta <- tX %*% X
rdeltapower <- rdelta^(1/2)
rdeltapowerinv <- solve(rdeltapower)
# information matrix
A <- rdelta - (1/nCols) * NR %*% tNR -(1/nRows)*NC %*% tNC + (1/(nRows*nCols)) * r %*% tr
A_star <- rdeltapowerinv %*% A %*% rdeltapowerinv
AnBn <- cbind(NR,NC,r)
YY <- NR # for the row blocking factor
ZZ <- NC # for the column blocking factor
if(flag){
InvA <- ginv(A)
InvAStar <- ginv(A_star)
}else{
InvA <- 0
InvAStar <- 0
}
if(rankMatrix(A)[1]==nTreatments-1){
connectedFlag <- 1
}
infoMat <- list(A,A_star,AnBn,InvA,InvAStar,YY,ZZ,rdeltapowerinv,connectedFlag,r)
names(infoMat) <- c("Cn","CnStar","AnBn","InvCn","InvAStar","YY","ZZ","rdeltapowerinv","connectedness","repli")
return(infoMat)
}
##################################
# return an information matrix after deleting a treatment in irow and icol
# input :
# output:
# Version 1.1
##################################
#' @param Design a given design
#' @param nTreatments the number of treatments
#' @param nRows the number of rows
#' @param nCols the number of columns
#' @param irow the element in a row will be deleted
#' @param icol the element in a column will be deleted
#' @return informatrix
#' @export
delInfoMatrix <- function(Design,nTreatments, nRows, nCols, irow, icol){
# Cn1
# get infomration for the orginal matrix
infoMat <- info2Matrix(Design,nTreatments, nRows, nCols)
Cn <- infoMat$Cn
AnBn <- infoMat$AnBn
InvCn <- infoMat$InvCn
# get new vectors
a1 <- rep(0,nTreatments)
a1[Design[irow,icol]] <- 1
b1 <- rep(0,nRows+nCols+1)
b1[irow] <- 1
b1[nRows+icol] <- 1
b1[nRows+nCols+1]<- 1
Ik <- (1/nCols) * diag(x = 1, nRows, nRows)
Ib <- (1/nRows) * diag(x = 1, nCols, nCols)
InvB <- bdiag(Ik,Ib,-1/(nCols*nRows))
#
xxx1 <- a1-AnBn %*% InvB %*% b1
xxx2 <- sqrt(1-t(b1) %*% InvB %*% b1)
c1 <- as.matrix(xxx1/xxx2[1])
Cn1 <- Cn - c1 %*% t(c1)
InvCn1 <- InvCn + (InvCn %*% c1 %*% t(c1) %*% InvCn)/(1-t(c1) %*% InvCn %*% c1)[1]
mlist <- list(Cn,Cn1,a1,b1,InvB,InvCn1)
names(mlist) <- c("Cn","Cn1","a1","b1","InvB","InvCn1")
return(mlist)
}
##################################
# return an information matrix after adding a treatment in irow and icol
# input :
# output:
# Version 1.1
##################################
#' @param Design a given design
#' @param nTreatments the number of treatments
#' @param nRows the number of rows
#' @param nCols the number of columns
#' @param irow the element in a row will be deleted
#' @param icol the element in a column will be deleted
#' @param itrt the treatment will be deleted
#' @return informatrix
#' @export
addInfoMatrix <- function(Design,nTreatments, nRows, nCols, irow, icol, itrt){
mlist <- delInfoMatrix(Design,nTreatments, nRows, nCols, irow, icol)
Cn1 <- mlist$Cn1
a1 <- mlist$a1
b1 <- mlist$b1
InvB <- as.matrix(mlist$InvB)
InvCn1 <- mlist$InvCn1
a2 <- rep(0,nTreatments)
a2[itrt]<- 1
b2 <- rep(0,nRows+nCols+1)
b2[irow] <- 1
b2[nRows+icol] <- 1
b2[nRows+nCols+1]<- 1
An_1Bn_1 <- AnBn-a1%*%t(b1)
InvB_1 <- InvB+ (InvB %*% b1 %*% t(b1) %*% InvB)/(1-t(b1) %*% InvB %*% b1)[1]
xxx1 <- a2-An_1Bn_1 %*% InvB_1 %*% b2
xxx2 <- sqrt(1+t(b2) %*% InvB_1 %*% b2)
c2 <- xxx1/xxx2[1]
Cn <- Cn1 + c2 %*% t(c2)
#InvCn1 <- ginv(Cn1)
InvCn <- InvCn1 - (InvCn1 %*% c2 %*% t(c2) %*% InvCn1)/(1+t(c2) %*% InvCn1 %*% c2)[1]
mlist <- list(Cn,Cn1,c2,a2,b2,InvB_1,InvCn)
names(mlist) <- c("Cn","Cn1","c2","a2","b2","InvB_1","InvCn")
return(mlist)
}
##################################
# 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
startDesignV1 <- function(nTreatments, nChecks, nRows, nCols, repli){
# 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
startDesignV2<- function(nTreatments, nChecks, nRows, nCols, repli, nchecks, vrepChecks, weighted=0, flag){
# 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)
#write.xlsx(pDesign, "pDesign1.xlsx",sheetName = "Step1")
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=0, reps=1,weighted)
}
#write.xlsx(pDesign, "pDesign2.xlsx",sheetName = "Step2")
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 <- aRCDesign(nDesign, nTreatments, nRows, nCols, repliz, ncheck=nChecks, flag=0, cRandom=0, reps=1,weighted)
# 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 <- aRCDesign(nDesign, nTreatments, nRows, nCols, repliz, ncheck=nChecks, flag=0, cRandom=0, reps=1,weighted)
}
#write.xlsx(nDesign, "pDesign3.xlsx",sheetName = "Step3")
}
optDesign <- aRCDesign(nDesign, nTreatments, nRows, nCols, repliz, ncheck=nChecks, flag=flag, cRandom=0, reps=1,weighted)
mlist <- list(nDesign,optDesign)
names(mlist) <- c("nDesign","optDesign")
return(mlist)
}
#' @param nTreatments
#' @param nRows
#' @param nCols
#' @param repliz
#' @param flag type of design, CRD or row-column design, or ...
#' @return a connected design
#' @export
genBlkDesign <- function(nTreatments, nRows, nCols, nChecks, vrepChecks, flag=0){
repliz <- rep(1, nTreatments)
randIndex <- sample(1:nTreatments, nChecks, replace = F)
repliz[randIndex] <- vrepChecks
pDesign <- StartDesign(nTreatments, nChecks, nRows, nCols, repliz)
## fixed flag for Aopt function-- can be changed
t <- AOpt(pDesign, nTreatments, nRows, nCols , repli, 0)
t <- 0
iter <- 0
niter <- 100
olditer <- 0
# Generating a connected design
t <- 0
while (t == 0) {
iter <- iter +1
for (ii in 1:nRows) {
indexx <- sample(1:nCols, nCols, replace = F)
pDesign[ii, ] <- pDesign[ii, indexx]
}
t <- AOpt(pDesign, nTreatments, nRows, nCols, repli=-1, 0)
if((iter==niter) && (t==0)){
repliz <- rep(1, nTreatments)
randIndex <- sample(1:nTreatments, nChecks, replace = F)
repliz[randIndex] <- vrepChecks
pDesign <- StartDesign(nTreatments, nChecks, nRows, nCols, repliz)
iter <- 0
olditer <- niter
}
}
return(pDesign)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
updateFormulaeT<- function(pDesign, infoMat, iRow1, iCol1, iRow2,iCol2, flag=0){
# get information of a given matrix
Cn1 <- infoMat$CnStar
InvCO1 <- infoMat$InvAStar
nRows <- dim(pDesign)[1]
nCols <- dim(pDesign)[2]
nTrts <- length(unique(c(pDesign)))
# a design after exchanging treatments
nDesign <- pDesign
nDesign[iRow1,iCol1] <- pDesign[iRow2,iCol2]
nDesign[iRow2,iCol2] <- pDesign[iRow1,iCol1]
infoMat2 <- info2Matrix(nDesign,nTrts, nRows, nCols)
zz <- 0
if(infoMat2$connectedness==1){
Cn2 <- infoMat2$CnStar
deltaM0 <- infoMat$Cn - infoMat2$Cn
deltaM1 <- Cn1-Cn2
if(sum(abs(deltaM1))!=0){
zz0 <- round(eigen(deltaM1)$values*10^10)/10^10
zz <- zz0[which(zz0!=0)]
if(length(zz)==2){
omega <- diag(zz)
}else{
omega <- zz
}
X <- eigen(deltaM1)$vectors[,which(zz0!=0)]
siM <- solve(omega)- t(X) %*%InvCO1 %*% X
invSiM <- solve(siM)
checkInvCO2 <- InvCO1 + InvCO1 %*% X %*% invSiM %*% t(X) %*% InvCO1
Echeck <- (nTrts-1)/sum(diag(checkInvCO2))
}else{
Echeck <- 0
checkInvCO2 <- 0
}
}else{
Echeck <- 0
checkInvCO2 <- 0
Cn2 <- 0
zz <- 0
deltaM0 <- 0
}
mlist <- list(Echeck,checkInvCO2,Cn2,zz,deltaM0)
names(mlist) <- c("Aopt","InvAStar","CnStar","Eigen","deltaM0")
return(mlist)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
Swapping<- function(pDesign,irow1,icol1,irow2,icol2){
copyDesign <- pDesign
xx <- pDesign[irow1, icol1]
yy <- pDesign[irow2, icol2]
copyDesign[irow1, icol1] <- yy
copyDesign[irow2, icol2] <- xx
return(copyDesign)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param infoMat
#' @param nTreatments
#' @param twocombn
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @param weighted = 0
#' @param flag = 0, 'u', 'c'
#' @export
updateFormulaeT3 <- function(pDesign, infoMat, nTreatments, nChecks, twoCombn, check2Combn , iRow1, iCol1, iRow2,iCol2, weighted = 0, flag=0){
#initialization
Echeck <- 0
checkInvCO2 <- 0
Cn2 <- 0
r <- 0
YY <- 0
ZZ <- 0
x0 <- 0
final <- 0
#initialization
# get information of a given matrix
r <- infoMat$repli
r_a <- sum(r)/nTreatments
r_square <- r_a*r_a
Cn1 <- infoMat$Cn
InvCO1 <- infoMat$InvCn
nRows <- dim(pDesign)[1]
nCols <- dim(pDesign)[2]
# get two different treatments will be exchanged
xx <- pDesign[iRow1, iCol1]
yy <- pDesign[iRow2, iCol2]
# for row-component design
alpha1 <- matrix(0,nTreatments,nRows)
alpha1[xx,iRow1] <- alpha1[xx,iRow1] - 1
alpha1[xx,iRow2] <- alpha1[xx,iRow2] + 1
alpha1[yy,iRow1] <- alpha1[yy,iRow1] + 1
alpha1[yy,iRow2] <- alpha1[yy,iRow2] - 1
#for column-component design
alpha2 <- matrix(0,nTreatments,nCols)
alpha2[xx,iCol1] <- alpha2[xx,iCol1] - 1
alpha2[xx,iCol2] <- alpha2[xx,iCol2] + 1
alpha2[yy,iCol1] <- alpha2[yy,iCol1] + 1
alpha2[yy,iCol2] <- alpha2[yy,iCol2] - 1
# calculate differences
#I need to check here
#print(sum(abs(alpha1)==0))
#print(sum(abs(alpha2)==0))
#if(sum(abs(alpha1)==0)){
#rtt <- matrix(0,nTreatments,nTreatments)
#}else{
rtt <- infoMat$YY %*% t(alpha1) + alpha1 %*%t(infoMat$YY) + alpha1 %*% t(alpha1)
#}
#if(sum(abs(alpha2)==0)){
#ctt <- matrix(0,nTreatments,nTreatments)
#}else{
ctt <- infoMat$ZZ %*% t(alpha2) + alpha2 %*%t(infoMat$ZZ) + alpha2 %*% t(alpha2)
#}
final <- (1/nCols)*rtt+(1/nRows)*ctt
YY <- infoMat$YY + alpha1 # for the row blocking factor
ZZ <- infoMat$ZZ + alpha2 # for the column blocking factor
Cn2 <- Cn1 - final
# check connectedness of a new design
diffdesign <- sum(abs(final))!=0
connected <- matrix.rank(Cn2)==(nTreatments-1)
if(diffdesign&&connected){
selColumn <- setdiff(c(1:nTreatments), c(xx,yy))
ordervec <- c(xx,yy,selColumn)
newdesign <- final[,ordervec]
newdesign <- newdesign[ordervec,]
newordervec <- matrix(0,1,nTreatments)
for(iii in 1:nTreatments){
newordervec[iii] <- c(which(ordervec==iii))
}
#get elements in a differences matrix
tz <- newdesign[3:nTreatments,1]
a1 <- -newdesign[1,2]
a2 <- newdesign[1,1] + newdesign[1,2]
tvalue <- t(tz)%*%tz
# calculate eigen values
lambda1 <- as.numeric(a1 + sqrt(a1*a1+a2*a2+2*tvalue))
lambda2 <- as.numeric(a1 - sqrt(a1*a1+a2*a2+2*tvalue))
# calculate eigen vectors
a <- newdesign[1,1]
b <- newdesign[1,2]
# There will be two different eigen vectors, which are denoted by p1 and p2
if(round(lambda1*10^10)/10^10!=0){
if((a+b+lambda1)==0){
x11 <- 0
x12 <- as.numeric((tvalue)/b)
x13 <- -t(tz)*x12/lambda1
}else{
rate1 <- (a+b-lambda1)/(a+b+lambda1)
if((lambda1-a-b*rate1)==0){
x11 <- as.numeric(sqrt((1-(tvalue))/(1+rate1*rate1)))
x12 <- as.numeric(rate1*x11)
x13 <- (t(tz)*x11-t(tz)*x12)/lambda1
}else{
x11 <- as.numeric((tvalue)/(lambda1-a-b*rate1))
x12 <- as.numeric(rate1*x11)
x13 <- (t(tz)*x11-t(tz)*x12)/lambda1
}
}
p1 <- c(x11,x12,x13)
}else{
p1 <- rep(0,nTreatments)
p1[1] <- 1
}
# vector p2
if(round(lambda2*10^10)/10^10!=0){
if((a+b+lambda2)==0){
x21 <- 0
x22 <- as.numeric((tvalue)/b)
x23 <- -t(tz)*x22/lambda2
}else{
rate2 <- (a+b-lambda2)/(a+b+lambda2)
if((lambda2-a-b*rate2)==0){
x21 <- as.numeric(sqrt((1-(tvalue))/(1+rate2*rate2)))
x22 <- as.numeric(rate2*x21)
x23 <- (t(tz)*x21-t(tz)*x22)/lambda2
}else{
x21 <- (tvalue)/(lambda2-a-b*rate2)
x22 <- rate2*x21[1]
x23 <- (t(tz)*x21[1]-t(tz)*x22[1])/lambda2
}
}
p2 <- c(x21,x22,x23)
}else{
p2 <- rep(0,nTreatments)
p2[1] <- 1
}
lambda <- round(c(lambda1,lambda2)*10^10)/10^10
ieigen <- which(lambda!=0)
reallambda <- lambda[ieigen]
ifelse(length(reallambda)==2,omega <- diag(reallambda),omega <- reallambda)
# normalized vectors
normp1 <- sqrt(1/(p1%*%p1))*p1
normp2 <- sqrt(1/(p2%*%p2))*p2
#get update for an inverse matrix
InvCO1 <- infoMat$InvCn
x0 <- cbind(matrix(normp1),matrix(normp2))
X <- x0[,ieigen]
ifelse(is.null(dim(X)),X <- t(t(X)), X <- X)
#get update for an inverse matrix
X <- X[newordervec,]
siM <- solve(omega)- t(X) %*%InvCO1 %*% X
invSiM <- solve(siM)
checkInvCO2 <- InvCO1 + InvCO1 %*% X %*% invSiM %*% t(X) %*% InvCO1
# add on 23 May 2017
# calculate A-efficiency factor
temp <- 0
#for (ii in c(1:(nTreatments-1))) {
# for(jj in c((ii+1):nTreatments)){
# #v <- (r[ii]*r[jj]/r_square)*(checkInvCO2[ii,ii]+checkInvCO2[jj,jj])
# v <- (r[ii]*r[jj]/r_square)*(checkInvCO2[ii,ii]+checkInvCO2[jj,jj]-2*checkInvCO2[ii,jj])
# temp2 <- temp2 + v
# }
#}
## incorporate with flag option
if ( weighted == 1 ) {
if (flag=="0") {
rarb <- r[twocombn[,1]]*r[twocombn[,2]]
diagElements <- diag(checkInvCO2)
AAA <- rarb*(diagElements[twocombn[,1]]+diagElements[twocombn[,2]])
BBB <- diag(as.vector(r))%*%checkInvCO2%*%diag(as.vector(r))
CCC <- sum(BBB)-sum(diag(BBB))
temp <- (sum(AAA) - CCC)/r_square
v_bar <- temp*2/(nTreatments*(nTreatments-1))
Aopt <- 2/(r_a*v_bar)
Echeck <- Aopt
} else if(flag=="u") {
nEntries <- nTreatments - nChecks
InvCO2Unrep <- checkInvCO2[c((nChecks+1):nTreatments),c((nChecks+1):nTreatments)]
temp <- nEntries*sum(diag(InvCO2Unrep))-sum(InvCO2Unrep)
v_bar <- temp*2/(nEntries*(nEntries-1))
Aopt <- 2/v_bar
Echeck <- Aopt
} else {
InvCO2Checks <- checkInvCO2[c(1:nChecks),c(1:nChecks)]
rCheck <- sum(r[1:nChecks])/nChecks
rSquareCheck <- rCheck*rCheck
temp <- 0
for (ii in c(1:(nChecks-1))) {
for(jj in c((ii+1):nChecks)){
v <- (r[ii]*r[jj]/rSquareCheck)*(InvCO2Checks[ii,ii]+InvCO2Checks[jj,jj]-2*InvCO2Checks[ii,jj])
temp <- temp + v
}
}
v_bar <- temp*2/(nChecks*(nChecks-1))
Aopt <- 2/(rCheck*v_bar)
Echeck <- Aopt
if(Echeck<0){
save(checkInvCO2,file = "INV.r")
print(InvCO2Checks)
print(solve(siM))
print(connected)
print(Echeck)
}else{
#print(InvCO2Checks)
}
}
} else {
temp <- nTreatments*sum(diag(checkInvCO2))-sum(checkInvCO2)
v_bar <- temp*2/(nTreatments*(nTreatments-1))
Aopt <- 2/(r_a*v_bar)
Echeck <- Aopt
}
} # diffdesign
mlist <- list(Echeck,checkInvCO2,Cn2,r,YY,ZZ,x0,final)
names(mlist) <- c("Aopt","InvCn","Cn","repli","YY","ZZ","Eigen","final")
return(mlist)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
layoutAugDesign <- function(pDesign,nTreatments,nChecks){
layoutDesign <- pDesign
for(i in (nChecks+1):nTreatments){
layoutDesign[pDesign==i] <- 0
}
return(layoutDesign)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
NeighboorUpdate <- function(optDeg,optA, inforDesign, iRow1, iCol1, iRow2, iCol2){
check <- updateFormulaeT3(optDeg, inforDesign, iRow1, iCol1, iRow2, iCol2)
optNew <- check$Aopt
ap <- exp((optNew-optA)/temper)
if ((ap > runif(1, 0, 1))&&(optNew>0)) {
optDeg <- Swapping(optDeg, iRow1, iCol1, iRow2, iCol2)
optA <- optNew
inforDesign <- check
}
mlist <- list(optDeg,optA,inforDesign)
names(mlist) <- c("optDeg","optA","inforDesign")
return(mlist)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
CoordinateExchange<- function(pDesign, nTreatments, nChecks , twocombn, check2Combn, nRows, nCols, infoDesign, cOpt, pos, weighted=0, flag){
# obtain row and column positions
irow <- pos[1]
icol <- pos[2]
tempCOpt <- cOpt
tempDesign <- pDesign
tempInfoDesign <- infoDesign
colPos <- icol+1
## ----> move to right middle
while ( colPos <= nCols ) {
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, irow, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
colPos <- colPos + 1
}
## ----> move to left middle
colPos <- icol-1
while(colPos>=1){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, irow, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
colPos <- colPos - 1
}
## ----> move to forward
rowPos <- irow-1
while(rowPos>=1){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, icol, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos - 1
}
## -----> move to backward
rowPos <- irow + 1
while(rowPos<=nRows){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, icol, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos + 1
}
# backward right
rowPos <- irow+1
colPos <- icol+1
while(((rowPos<=nRows))&&(colPos<=nCols)){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos + 1
colPos <- colPos + 1
}
# forward left
rowPos <- irow-1
colPos <- icol-1
while(((rowPos>=1))&&(colPos>=1)){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos - 1
colPos <- colPos - 1
}
# forward right
rowPos <- irow-1
colPos <- icol+1
while(((rowPos>=1))&&(colPos<=nCols)){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos - 1
colPos <- colPos + 1
}
# backward left
rowPos <- irow+1
colPos <- icol-1
while(((rowPos<=nRows))&&(colPos>=1)){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos + 1
colPos <- colPos - 1
}
mlist <- list(tempDesign,tempInfoDesign,tempCOpt)
names(mlist) <- c("pDesign","infoDesign","Aopt")
return(mlist)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
CompareSwappingDesign <- function(pDesign, infoDesign, cDesign, cOpt, cInfoDesign , nTreatments, nChecks , twocombn, check2Combn, iRow1, iCol1, iRow2, iCol2, weighted, flag){
iTrt1 <- pDesign[iRow1, iCol1]
iTrt2 <- pDesign[iRow2, iCol2]
# set variables
tempDesign <- cDesign
tempCOpt <- cOpt
tempInfoDesign <- cInfoDesign
if(iTrt1!=iTrt2){
uptDesign <- updateFormulaeT3(pDesign, infoDesign, nTreatments, nChecks , twocombn, check2Combn, iRow1, iCol1, iRow2, iCol2, weighted,flag) #print(uptDesign$Aopt)
if (uptDesign$Aopt> tempCOpt) {
tempCOpt <- uptDesign$Aopt
tempDesign <- Swapping(pDesign, iRow1, iCol1, iRow2, iCol2)
tempInfoDesign <- uptDesign
}
}
mlist <- list(tempDesign,tempInfoDesign,tempCOpt)
names(mlist) <- c("pDesign","infoDesign","Aopt")
return(mlist)
}
#' Swapping two different plots in a given design
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
SwappingTwoPlots <- function(pDesign, infoDesign, AOpt, nTreatments, nChecks , twocombn, check2Combn, iRow1, iCol1, iRow2, iCol2, weighted, flag){
iTrt1 <- pDesign[iRow1, iCol1]
iTrt2 <- pDesign[iRow2, iCol2]
# set variables
tempDesign <- pDesign
tempCOpt <- round(AOpt*10^10)/10^10
tempInfoDesign <- infoDesign
if(iTrt1!=iTrt2){
uptDesign <- updateFormulaeT3(pDesign, infoDesign, nTreatments, nChecks , twocombn, check2Combn, iRow1, iCol1, iRow2, iCol2, weighted,flag)
#print(uptDesign$Aopt)
if (round(uptDesign$Aopt*10^10)/10^10> tempCOpt) {
tempCOpt <- uptDesign$Aopt
tempDesign <- Swapping(pDesign, iRow1, iCol1, iRow2, iCol2)
tempInfoDesign <- uptDesign
}
}
mlist <- list(tempDesign,tempInfoDesign,tempCOpt)
names(mlist) <- c("pDesign","infoDesign","Aopt")
return(mlist)
}
RandomSearch <- function(pDesign,nTreatments, nRows, nCols, optDesign, infoDesign, twocombn, check2Combn, weighted, flag, niter){
copyDesign <- pDesign
iter <- 0
recordImp <- rep(0,niter)
while (iter < niter) {
cTreatments <- as.matrix(rbind(which(pDesign==1,arr.ind=T),which(pDesign==2,arr.ind=T),which(pDesign==3,arr.ind=T),which(pDesign==4,arr.ind=T)))
icheck <- sample(1:lcTreatment, 1, replace = T)
siter <- 0
while (siter < 1000) {
irow <- c(as.double(cTreatments[icheck,1]),sample(1:nRows, 1, replace = T))
icol <- c(as.double(cTreatments[icheck,2]),sample(1:nCols, 1, replace = T))
if (pDesign[irow[1], icol[1]] != pDesign[irow[2], icol[2]]) {
check <- updateFormulaeT3(pDesign, infoDesign, nTreatments, nChecks, twocombn, check2Combn, irow[1], icol[1], irow[2], icol[2], weighted,flag)
if (round(check$Aopt*10^10)/10^10 > round(optDesign*10^10)/10^10) {
#print('I am here')
optDesign <- check$Aopt
# We need to swap treatments
iPDesign <- Swapping(pDesign, irow[1], icol[1], irow[2], icol[2])
pDesign <- iPDesign
infoDesign <- check
}
}
siter <- siter + 1
}
recordImp[iter] <- optDesign
iter <- iter + 1
}
mlist <- list(optDesign,pDesign,infoDesign,recordImp)
names(mlist) <- c("optDesign","pDesign","infoDesign","recordImp")
return(mlist)
}
#' Swapping two different plots in a given design
#' @param InvInfoMatrix
#' @param TwoCombn
#' @param Repli
#' @param Weight
#' @param Flag
#' @export
AeffRowColumnDesign <- function(InvInfoMatrix, nTreatments, nChecks, r_a, ...){
# get variables
arguments <- list(...)
TwoCombn <- arguments$TwoCombn
ifelse(!is.null(arguments$TwoCombn),TwoCombn <- arguments$TwoCombn,TwoCombn <- -1)
ifelse(!is.null(arguments$Repli),Repli <- arguments$Repli, Repli <- -1)
ifelse(!is.null(arguments$Weight),Weight <- arguments$Weight, Weight <- -1)
ifelse(!is.null(arguments$Flag),Flag <- arguments$Flag,Flag <- -1)
# set variables
v_bar <- 0
Aopt <- 0
Echeck <- 0
r_square <- r_a*r_a
if ( Weight == 1 ) {
if (Flag=="0") {
rarb <- Repli[TwoCombn[,1]]*Repli[TwoCombn[,2]]
diagElements <- diag(InvInfoMatrix)
AAA <- rarb*(diagElements[TwoCombn[,1]]+diagElements[TwoCombn[,2]])
BBB <- diag(as.vector(Repli))%*%InvInfoMatrix%*%diag(as.vector(Repli))
CCC <- sum(BBB)-sum(diag(BBB))
temp <- (sum(AAA) - CCC)/r_square
v_bar <- temp*2/(nTreatments*(nTreatments-1))
Aopt <- 2/(r_a*v_bar)
Echeck <- Aopt
} else if(Flag=="u") {
nEntries <- nTreatments - nChecks
InvCOUnrep <- InvInfoMatrix[c((nChecks+1):nTreatments),c((nChecks+1):nTreatments)]
temp <- nEntries*sum(diag(InvCOUnrep))-sum(InvCOUnrep)
v_bar <- temp*2/(nEntries*(nEntries-1))
Aopt <- 2/v_bar
Echeck <- Aopt
} else if(Flag=="c") {
InvCO2Checks <- InvInfoMatrix[c(1:nChecks),c(1:nChecks)]
rCheck <- sum(r[1:nChecks])/nChecks
rSquareCheck <- rCheck*rCheck
temp <- 0
for (ii in c(1:(nChecks-1))) {
for(jj in c((ii+1):nChecks)){
v <- (r[ii]*r[jj]/rSquareCheck)*(InvCO2Checks[ii,ii]+InvCO2Checks[jj,jj]-2*InvCO2Checks[ii,jj])
temp <- temp + v
}
}
v_bar <- temp*2/(nChecks*(nChecks-1))
Aopt <- 2/(rCheck*v_bar)
Echeck <- Aopt
}
} else {
temp <- nTreatments*sum(diag(InvInfoMatrix))-sum(InvInfoMatrix)
v_bar <- temp*2/(nTreatments*(nTreatments-1))
Aopt <- 2/(r_a*v_bar)
Echeck <- Aopt
}
#mlist <- list(optDesign,pDesign,infoDesign,recordImp)
#names(mlist) <- c("optDesign","pDesign","infoDesign","recordImp")
return(Echeck)
}
##----------------------------------##
## A place to test your codes
##----------------------------------##
#setwd("D:/user/vthanh/Google Drive/Project Stuttgart/R_programing/UpdateFormular")
#source("D:/user/vthanh/Google Drive/Project Stuttgart/R_programing/A_Eff_RowColumnDesign.r")
#data <- split(read.table("./Designs/Design.txt"), gl(1, 4))
#repli <- 3
#pDesign <- data$`1`
#modelD <- modelMatrix(10, 4, 4, pDesign)
nhavt/Block3AugDesigns documentation built on May 7, 2019, 11:15 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions model2Matrix model2MatrixV2 info2Matrix delInfoMatrix addInfoMatrix startDesignV1 startDesignV2 genBlkDesign updateFormulaeT Swapping updateFormulaeT3 layoutAugDesign NeighboorUpdate CoordinateExchange CompareSwappingDesign SwappingTwoPlots RandomSearch AeffRowColumnDesign
#' Develop by Nha Vo-Thanh
#' Stuttgart on 22 March 2017
#' Augmented row-column designs
#' Version 1.0
#'
##################################
#' create a model matrix
#' @param nTreatments number of treatments
#' @param nRows number of rows
#' @param nCols number of columns
#' @return a model matrix
#' @export
model2Matrix <- function(nTreatments, nRows, nCols, A)
{
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)
model <- list(trtModel,rModel,cModel)
names(model) <- c("X","R","C")
return(model)
}
##################################
#' create a model matrix
#' @param nTreatments number of treatments
#' @param nRows number of rows
#' @param nCols number of columns
#' @return a model matrix
#' @export
model2MatrixV2 <- function(nTreatments, nRows, nCols, A)
{
trtModel <- matrix(0, nRows * nCols, nTreatments)
vecA <- c(t(A))
mylist <- rep(list(vecA),nTreatments)
fxn <- function(var1,var2){
as.numeric(var1==var2)
}
var2 <- as.list(c(1:nTreatments))
trtModel <- mapply(fxn,mylist,var2)
## 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)
model <- list(trtModel,rModel,cModel)
names(model) <- c("X","R","C")
return(model)
}
##################################
# return an information matrix
# input :
# output:
# Version 1.1
##################################
#' @param Design a given design
#' @param nTreatments the number of treatments
#' @param nRows the number of rows
#' @param nCols the number of columns
#' @param flag a flag to decide whether to include inverse matrix
#' @return information matrix
#' @export
info2Matrix <- function(Design,nTreatments, nRows, nCols, flag=0){
# connectedFlag
connectedFlag <- 0
#modelA <- model2Matrix(nTreatments, nRows, nCols, Design)
modelA <- model2MatrixV2(nTreatments, nRows, nCols, Design)
# information for treatments
X <- modelA$X
# information for rows and columns
R <- modelA$R
C <- modelA$C
# tranpose matrix
tX <- t(X)
tR <- t(R)
tC <- t(C)
# incidence matrix
NR <- tX %*% R
NC <- tX %*% C
# tranpose matrix
tNR <- t(NR)
tNC <- t(NC)
# creates a vector of ones
onevec <- matrix(1, nrow = dim(tX)[2], ncol = 1)
r <- tX %*% onevec
tr <- t(r)
rdelta <- tX %*% X
rdeltapower <- rdelta^(1/2)
rdeltapowerinv <- solve(rdeltapower)
# information matrix
A <- rdelta - (1/nCols) * NR %*% tNR -(1/nRows)*NC %*% tNC + (1/(nRows*nCols)) * r %*% tr
A_star <- rdeltapowerinv %*% A %*% rdeltapowerinv
AnBn <- cbind(NR,NC,r)
YY <- NR # for the row blocking factor
ZZ <- NC # for the column blocking factor
if(flag){
InvA <- ginv(A)
InvAStar <- ginv(A_star)
}else{
InvA <- 0
InvAStar <- 0
}
if(rankMatrix(A)[1]==nTreatments-1){
connectedFlag <- 1
}
infoMat <- list(A,A_star,AnBn,InvA,InvAStar,YY,ZZ,rdeltapowerinv,connectedFlag,r)
names(infoMat) <- c("Cn","CnStar","AnBn","InvCn","InvAStar","YY","ZZ","rdeltapowerinv","connectedness","repli")
return(infoMat)
}
##################################
# return an information matrix after deleting a treatment in irow and icol
# input :
# output:
# Version 1.1
##################################
#' @param Design a given design
#' @param nTreatments the number of treatments
#' @param nRows the number of rows
#' @param nCols the number of columns
#' @param irow the element in a row will be deleted
#' @param icol the element in a column will be deleted
#' @return informatrix
#' @export
delInfoMatrix <- function(Design,nTreatments, nRows, nCols, irow, icol){
# Cn1
# get infomration for the orginal matrix
infoMat <- info2Matrix(Design,nTreatments, nRows, nCols)
Cn <- infoMat$Cn
AnBn <- infoMat$AnBn
InvCn <- infoMat$InvCn
# get new vectors
a1 <- rep(0,nTreatments)
a1[Design[irow,icol]] <- 1
b1 <- rep(0,nRows+nCols+1)
b1[irow] <- 1
b1[nRows+icol] <- 1
b1[nRows+nCols+1]<- 1
Ik <- (1/nCols) * diag(x = 1, nRows, nRows)
Ib <- (1/nRows) * diag(x = 1, nCols, nCols)
InvB <- bdiag(Ik,Ib,-1/(nCols*nRows))
#
xxx1 <- a1-AnBn %*% InvB %*% b1
xxx2 <- sqrt(1-t(b1) %*% InvB %*% b1)
c1 <- as.matrix(xxx1/xxx2[1])
Cn1 <- Cn - c1 %*% t(c1)
InvCn1 <- InvCn + (InvCn %*% c1 %*% t(c1) %*% InvCn)/(1-t(c1) %*% InvCn %*% c1)[1]
mlist <- list(Cn,Cn1,a1,b1,InvB,InvCn1)
names(mlist) <- c("Cn","Cn1","a1","b1","InvB","InvCn1")
return(mlist)
}
##################################
# return an information matrix after adding a treatment in irow and icol
# input :
# output:
# Version 1.1
##################################
#' @param Design a given design
#' @param nTreatments the number of treatments
#' @param nRows the number of rows
#' @param nCols the number of columns
#' @param irow the element in a row will be deleted
#' @param icol the element in a column will be deleted
#' @param itrt the treatment will be deleted
#' @return informatrix
#' @export
addInfoMatrix <- function(Design,nTreatments, nRows, nCols, irow, icol, itrt){
mlist <- delInfoMatrix(Design,nTreatments, nRows, nCols, irow, icol)
Cn1 <- mlist$Cn1
a1 <- mlist$a1
b1 <- mlist$b1
InvB <- as.matrix(mlist$InvB)
InvCn1 <- mlist$InvCn1
a2 <- rep(0,nTreatments)
a2[itrt]<- 1
b2 <- rep(0,nRows+nCols+1)
b2[irow] <- 1
b2[nRows+icol] <- 1
b2[nRows+nCols+1]<- 1
An_1Bn_1 <- AnBn-a1%*%t(b1)
InvB_1 <- InvB+ (InvB %*% b1 %*% t(b1) %*% InvB)/(1-t(b1) %*% InvB %*% b1)[1]
xxx1 <- a2-An_1Bn_1 %*% InvB_1 %*% b2
xxx2 <- sqrt(1+t(b2) %*% InvB_1 %*% b2)
c2 <- xxx1/xxx2[1]
Cn <- Cn1 + c2 %*% t(c2)
#InvCn1 <- ginv(Cn1)
InvCn <- InvCn1 - (InvCn1 %*% c2 %*% t(c2) %*% InvCn1)/(1+t(c2) %*% InvCn1 %*% c2)[1]
mlist <- list(Cn,Cn1,c2,a2,b2,InvB_1,InvCn)
names(mlist) <- c("Cn","Cn1","c2","a2","b2","InvB_1","InvCn")
return(mlist)
}
##################################
# 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
startDesignV1 <- function(nTreatments, nChecks, nRows, nCols, repli){
# 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
startDesignV2<- function(nTreatments, nChecks, nRows, nCols, repli, nchecks, vrepChecks, weighted=0, flag){
# 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)
#write.xlsx(pDesign, "pDesign1.xlsx",sheetName = "Step1")
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=0, reps=1,weighted)
}
#write.xlsx(pDesign, "pDesign2.xlsx",sheetName = "Step2")
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 <- aRCDesign(nDesign, nTreatments, nRows, nCols, repliz, ncheck=nChecks, flag=0, cRandom=0, reps=1,weighted)
# 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 <- aRCDesign(nDesign, nTreatments, nRows, nCols, repliz, ncheck=nChecks, flag=0, cRandom=0, reps=1,weighted)
}
#write.xlsx(nDesign, "pDesign3.xlsx",sheetName = "Step3")
}
optDesign <- aRCDesign(nDesign, nTreatments, nRows, nCols, repliz, ncheck=nChecks, flag=flag, cRandom=0, reps=1,weighted)
mlist <- list(nDesign,optDesign)
names(mlist) <- c("nDesign","optDesign")
return(mlist)
}
#' @param nTreatments
#' @param nRows
#' @param nCols
#' @param repliz
#' @param flag type of design, CRD or row-column design, or ...
#' @return a connected design
#' @export
genBlkDesign <- function(nTreatments, nRows, nCols, nChecks, vrepChecks, flag=0){
repliz <- rep(1, nTreatments)
randIndex <- sample(1:nTreatments, nChecks, replace = F)
repliz[randIndex] <- vrepChecks
pDesign <- StartDesign(nTreatments, nChecks, nRows, nCols, repliz)
## fixed flag for Aopt function-- can be changed
t <- AOpt(pDesign, nTreatments, nRows, nCols , repli, 0)
t <- 0
iter <- 0
niter <- 100
olditer <- 0
# Generating a connected design
t <- 0
while (t == 0) {
iter <- iter +1
for (ii in 1:nRows) {
indexx <- sample(1:nCols, nCols, replace = F)
pDesign[ii, ] <- pDesign[ii, indexx]
}
t <- AOpt(pDesign, nTreatments, nRows, nCols, repli=-1, 0)
if((iter==niter) && (t==0)){
repliz <- rep(1, nTreatments)
randIndex <- sample(1:nTreatments, nChecks, replace = F)
repliz[randIndex] <- vrepChecks
pDesign <- StartDesign(nTreatments, nChecks, nRows, nCols, repliz)
iter <- 0
olditer <- niter
}
}
return(pDesign)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
updateFormulaeT<- function(pDesign, infoMat, iRow1, iCol1, iRow2,iCol2, flag=0){
# get information of a given matrix
Cn1 <- infoMat$CnStar
InvCO1 <- infoMat$InvAStar
nRows <- dim(pDesign)[1]
nCols <- dim(pDesign)[2]
nTrts <- length(unique(c(pDesign)))
# a design after exchanging treatments
nDesign <- pDesign
nDesign[iRow1,iCol1] <- pDesign[iRow2,iCol2]
nDesign[iRow2,iCol2] <- pDesign[iRow1,iCol1]
infoMat2 <- info2Matrix(nDesign,nTrts, nRows, nCols)
zz <- 0
if(infoMat2$connectedness==1){
Cn2 <- infoMat2$CnStar
deltaM0 <- infoMat$Cn - infoMat2$Cn
deltaM1 <- Cn1-Cn2
if(sum(abs(deltaM1))!=0){
zz0 <- round(eigen(deltaM1)$values*10^10)/10^10
zz <- zz0[which(zz0!=0)]
if(length(zz)==2){
omega <- diag(zz)
}else{
omega <- zz
}
X <- eigen(deltaM1)$vectors[,which(zz0!=0)]
siM <- solve(omega)- t(X) %*%InvCO1 %*% X
invSiM <- solve(siM)
checkInvCO2 <- InvCO1 + InvCO1 %*% X %*% invSiM %*% t(X) %*% InvCO1
Echeck <- (nTrts-1)/sum(diag(checkInvCO2))
}else{
Echeck <- 0
checkInvCO2 <- 0
}
}else{
Echeck <- 0
checkInvCO2 <- 0
Cn2 <- 0
zz <- 0
deltaM0 <- 0
}
mlist <- list(Echeck,checkInvCO2,Cn2,zz,deltaM0)
names(mlist) <- c("Aopt","InvAStar","CnStar","Eigen","deltaM0")
return(mlist)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
Swapping<- function(pDesign,irow1,icol1,irow2,icol2){
copyDesign <- pDesign
xx <- pDesign[irow1, icol1]
yy <- pDesign[irow2, icol2]
copyDesign[irow1, icol1] <- yy
copyDesign[irow2, icol2] <- xx
return(copyDesign)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param infoMat
#' @param nTreatments
#' @param twocombn
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @param weighted = 0
#' @param flag = 0, 'u', 'c'
#' @export
updateFormulaeT3 <- function(pDesign, infoMat, nTreatments, nChecks, twoCombn, check2Combn , iRow1, iCol1, iRow2,iCol2, weighted = 0, flag=0){
#initialization
Echeck <- 0
checkInvCO2 <- 0
Cn2 <- 0
r <- 0
YY <- 0
ZZ <- 0
x0 <- 0
final <- 0
#initialization
# get information of a given matrix
r <- infoMat$repli
r_a <- sum(r)/nTreatments
r_square <- r_a*r_a
Cn1 <- infoMat$Cn
InvCO1 <- infoMat$InvCn
nRows <- dim(pDesign)[1]
nCols <- dim(pDesign)[2]
# get two different treatments will be exchanged
xx <- pDesign[iRow1, iCol1]
yy <- pDesign[iRow2, iCol2]
# for row-component design
alpha1 <- matrix(0,nTreatments,nRows)
alpha1[xx,iRow1] <- alpha1[xx,iRow1] - 1
alpha1[xx,iRow2] <- alpha1[xx,iRow2] + 1
alpha1[yy,iRow1] <- alpha1[yy,iRow1] + 1
alpha1[yy,iRow2] <- alpha1[yy,iRow2] - 1
#for column-component design
alpha2 <- matrix(0,nTreatments,nCols)
alpha2[xx,iCol1] <- alpha2[xx,iCol1] - 1
alpha2[xx,iCol2] <- alpha2[xx,iCol2] + 1
alpha2[yy,iCol1] <- alpha2[yy,iCol1] + 1
alpha2[yy,iCol2] <- alpha2[yy,iCol2] - 1
# calculate differences
#I need to check here
#print(sum(abs(alpha1)==0))
#print(sum(abs(alpha2)==0))
#if(sum(abs(alpha1)==0)){
#rtt <- matrix(0,nTreatments,nTreatments)
#}else{
rtt <- infoMat$YY %*% t(alpha1) + alpha1 %*%t(infoMat$YY) + alpha1 %*% t(alpha1)
#}
#if(sum(abs(alpha2)==0)){
#ctt <- matrix(0,nTreatments,nTreatments)
#}else{
ctt <- infoMat$ZZ %*% t(alpha2) + alpha2 %*%t(infoMat$ZZ) + alpha2 %*% t(alpha2)
#}
final <- (1/nCols)*rtt+(1/nRows)*ctt
YY <- infoMat$YY + alpha1 # for the row blocking factor
ZZ <- infoMat$ZZ + alpha2 # for the column blocking factor
Cn2 <- Cn1 - final
# check connectedness of a new design
diffdesign <- sum(abs(final))!=0
connected <- matrix.rank(Cn2)==(nTreatments-1)
if(diffdesign&&connected){
selColumn <- setdiff(c(1:nTreatments), c(xx,yy))
ordervec <- c(xx,yy,selColumn)
newdesign <- final[,ordervec]
newdesign <- newdesign[ordervec,]
newordervec <- matrix(0,1,nTreatments)
for(iii in 1:nTreatments){
newordervec[iii] <- c(which(ordervec==iii))
}
#get elements in a differences matrix
tz <- newdesign[3:nTreatments,1]
a1 <- -newdesign[1,2]
a2 <- newdesign[1,1] + newdesign[1,2]
tvalue <- t(tz)%*%tz
# calculate eigen values
lambda1 <- as.numeric(a1 + sqrt(a1*a1+a2*a2+2*tvalue))
lambda2 <- as.numeric(a1 - sqrt(a1*a1+a2*a2+2*tvalue))
# calculate eigen vectors
a <- newdesign[1,1]
b <- newdesign[1,2]
# There will be two different eigen vectors, which are denoted by p1 and p2
if(round(lambda1*10^10)/10^10!=0){
if((a+b+lambda1)==0){
x11 <- 0
x12 <- as.numeric((tvalue)/b)
x13 <- -t(tz)*x12/lambda1
}else{
rate1 <- (a+b-lambda1)/(a+b+lambda1)
if((lambda1-a-b*rate1)==0){
x11 <- as.numeric(sqrt((1-(tvalue))/(1+rate1*rate1)))
x12 <- as.numeric(rate1*x11)
x13 <- (t(tz)*x11-t(tz)*x12)/lambda1
}else{
x11 <- as.numeric((tvalue)/(lambda1-a-b*rate1))
x12 <- as.numeric(rate1*x11)
x13 <- (t(tz)*x11-t(tz)*x12)/lambda1
}
}
p1 <- c(x11,x12,x13)
}else{
p1 <- rep(0,nTreatments)
p1[1] <- 1
}
# vector p2
if(round(lambda2*10^10)/10^10!=0){
if((a+b+lambda2)==0){
x21 <- 0
x22 <- as.numeric((tvalue)/b)
x23 <- -t(tz)*x22/lambda2
}else{
rate2 <- (a+b-lambda2)/(a+b+lambda2)
if((lambda2-a-b*rate2)==0){
x21 <- as.numeric(sqrt((1-(tvalue))/(1+rate2*rate2)))
x22 <- as.numeric(rate2*x21)
x23 <- (t(tz)*x21-t(tz)*x22)/lambda2
}else{
x21 <- (tvalue)/(lambda2-a-b*rate2)
x22 <- rate2*x21[1]
x23 <- (t(tz)*x21[1]-t(tz)*x22[1])/lambda2
}
}
p2 <- c(x21,x22,x23)
}else{
p2 <- rep(0,nTreatments)
p2[1] <- 1
}
lambda <- round(c(lambda1,lambda2)*10^10)/10^10
ieigen <- which(lambda!=0)
reallambda <- lambda[ieigen]
ifelse(length(reallambda)==2,omega <- diag(reallambda),omega <- reallambda)
# normalized vectors
normp1 <- sqrt(1/(p1%*%p1))*p1
normp2 <- sqrt(1/(p2%*%p2))*p2
#get update for an inverse matrix
InvCO1 <- infoMat$InvCn
x0 <- cbind(matrix(normp1),matrix(normp2))
X <- x0[,ieigen]
ifelse(is.null(dim(X)),X <- t(t(X)), X <- X)
#get update for an inverse matrix
X <- X[newordervec,]
siM <- solve(omega)- t(X) %*%InvCO1 %*% X
invSiM <- solve(siM)
checkInvCO2 <- InvCO1 + InvCO1 %*% X %*% invSiM %*% t(X) %*% InvCO1
# add on 23 May 2017
# calculate A-efficiency factor
temp <- 0
#for (ii in c(1:(nTreatments-1))) {
# for(jj in c((ii+1):nTreatments)){
# #v <- (r[ii]*r[jj]/r_square)*(checkInvCO2[ii,ii]+checkInvCO2[jj,jj])
# v <- (r[ii]*r[jj]/r_square)*(checkInvCO2[ii,ii]+checkInvCO2[jj,jj]-2*checkInvCO2[ii,jj])
# temp2 <- temp2 + v
# }
#}
## incorporate with flag option
if ( weighted == 1 ) {
if (flag=="0") {
rarb <- r[twocombn[,1]]*r[twocombn[,2]]
diagElements <- diag(checkInvCO2)
AAA <- rarb*(diagElements[twocombn[,1]]+diagElements[twocombn[,2]])
BBB <- diag(as.vector(r))%*%checkInvCO2%*%diag(as.vector(r))
CCC <- sum(BBB)-sum(diag(BBB))
temp <- (sum(AAA) - CCC)/r_square
v_bar <- temp*2/(nTreatments*(nTreatments-1))
Aopt <- 2/(r_a*v_bar)
Echeck <- Aopt
} else if(flag=="u") {
nEntries <- nTreatments - nChecks
InvCO2Unrep <- checkInvCO2[c((nChecks+1):nTreatments),c((nChecks+1):nTreatments)]
temp <- nEntries*sum(diag(InvCO2Unrep))-sum(InvCO2Unrep)
v_bar <- temp*2/(nEntries*(nEntries-1))
Aopt <- 2/v_bar
Echeck <- Aopt
} else {
InvCO2Checks <- checkInvCO2[c(1:nChecks),c(1:nChecks)]
rCheck <- sum(r[1:nChecks])/nChecks
rSquareCheck <- rCheck*rCheck
temp <- 0
for (ii in c(1:(nChecks-1))) {
for(jj in c((ii+1):nChecks)){
v <- (r[ii]*r[jj]/rSquareCheck)*(InvCO2Checks[ii,ii]+InvCO2Checks[jj,jj]-2*InvCO2Checks[ii,jj])
temp <- temp + v
}
}
v_bar <- temp*2/(nChecks*(nChecks-1))
Aopt <- 2/(rCheck*v_bar)
Echeck <- Aopt
if(Echeck<0){
save(checkInvCO2,file = "INV.r")
print(InvCO2Checks)
print(solve(siM))
print(connected)
print(Echeck)
}else{
#print(InvCO2Checks)
}
}
} else {
temp <- nTreatments*sum(diag(checkInvCO2))-sum(checkInvCO2)
v_bar <- temp*2/(nTreatments*(nTreatments-1))
Aopt <- 2/(r_a*v_bar)
Echeck <- Aopt
}
} # diffdesign
mlist <- list(Echeck,checkInvCO2,Cn2,r,YY,ZZ,x0,final)
names(mlist) <- c("Aopt","InvCn","Cn","repli","YY","ZZ","Eigen","final")
return(mlist)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
layoutAugDesign <- function(pDesign,nTreatments,nChecks){
layoutDesign <- pDesign
for(i in (nChecks+1):nTreatments){
layoutDesign[pDesign==i] <- 0
}
return(layoutDesign)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
NeighboorUpdate <- function(optDeg,optA, inforDesign, iRow1, iCol1, iRow2, iCol2){
check <- updateFormulaeT3(optDeg, inforDesign, iRow1, iCol1, iRow2, iCol2)
optNew <- check$Aopt
ap <- exp((optNew-optA)/temper)
if ((ap > runif(1, 0, 1))&&(optNew>0)) {
optDeg <- Swapping(optDeg, iRow1, iCol1, iRow2, iCol2)
optA <- optNew
inforDesign <- check
}
mlist <- list(optDeg,optA,inforDesign)
names(mlist) <- c("optDeg","optA","inforDesign")
return(mlist)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
CoordinateExchange<- function(pDesign, nTreatments, nChecks , twocombn, check2Combn, nRows, nCols, infoDesign, cOpt, pos, weighted=0, flag){
# obtain row and column positions
irow <- pos[1]
icol <- pos[2]
tempCOpt <- cOpt
tempDesign <- pDesign
tempInfoDesign <- infoDesign
colPos <- icol+1
## ----> move to right middle
while ( colPos <= nCols ) {
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, irow, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
colPos <- colPos + 1
}
## ----> move to left middle
colPos <- icol-1
while(colPos>=1){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, irow, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
colPos <- colPos - 1
}
## ----> move to forward
rowPos <- irow-1
while(rowPos>=1){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, icol, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos - 1
}
## -----> move to backward
rowPos <- irow + 1
while(rowPos<=nRows){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, icol, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos + 1
}
# backward right
rowPos <- irow+1
colPos <- icol+1
while(((rowPos<=nRows))&&(colPos<=nCols)){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos + 1
colPos <- colPos + 1
}
# forward left
rowPos <- irow-1
colPos <- icol-1
while(((rowPos>=1))&&(colPos>=1)){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos - 1
colPos <- colPos - 1
}
# forward right
rowPos <- irow-1
colPos <- icol+1
while(((rowPos>=1))&&(colPos<=nCols)){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos - 1
colPos <- colPos + 1
}
# backward left
rowPos <- irow+1
colPos <- icol-1
while(((rowPos<=nRows))&&(colPos>=1)){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos + 1
colPos <- colPos - 1
}
mlist <- list(tempDesign,tempInfoDesign,tempCOpt)
names(mlist) <- c("pDesign","infoDesign","Aopt")
return(mlist)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
CompareSwappingDesign <- function(pDesign, infoDesign, cDesign, cOpt, cInfoDesign , nTreatments, nChecks , twocombn, check2Combn, iRow1, iCol1, iRow2, iCol2, weighted, flag){
iTrt1 <- pDesign[iRow1, iCol1]
iTrt2 <- pDesign[iRow2, iCol2]
# set variables
tempDesign <- cDesign
tempCOpt <- cOpt
tempInfoDesign <- cInfoDesign
if(iTrt1!=iTrt2){
uptDesign <- updateFormulaeT3(pDesign, infoDesign, nTreatments, nChecks , twocombn, check2Combn, iRow1, iCol1, iRow2, iCol2, weighted,flag) #print(uptDesign$Aopt)
if (uptDesign$Aopt> tempCOpt) {
tempCOpt <- uptDesign$Aopt
tempDesign <- Swapping(pDesign, iRow1, iCol1, iRow2, iCol2)
tempInfoDesign <- uptDesign
}
}
mlist <- list(tempDesign,tempInfoDesign,tempCOpt)
names(mlist) <- c("pDesign","infoDesign","Aopt")
return(mlist)
}
#' Swapping two different plots in a given design
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
SwappingTwoPlots <- function(pDesign, infoDesign, AOpt, nTreatments, nChecks , twocombn, check2Combn, iRow1, iCol1, iRow2, iCol2, weighted, flag){
iTrt1 <- pDesign[iRow1, iCol1]
iTrt2 <- pDesign[iRow2, iCol2]
# set variables
tempDesign <- pDesign
tempCOpt <- round(AOpt*10^10)/10^10
tempInfoDesign <- infoDesign
if(iTrt1!=iTrt2){
uptDesign <- updateFormulaeT3(pDesign, infoDesign, nTreatments, nChecks , twocombn, check2Combn, iRow1, iCol1, iRow2, iCol2, weighted,flag)
#print(uptDesign$Aopt)
if (round(uptDesign$Aopt*10^10)/10^10> tempCOpt) {
tempCOpt <- uptDesign$Aopt
tempDesign <- Swapping(pDesign, iRow1, iCol1, iRow2, iCol2)
tempInfoDesign <- uptDesign
}
}
mlist <- list(tempDesign,tempInfoDesign,tempCOpt)
names(mlist) <- c("pDesign","infoDesign","Aopt")
return(mlist)
}
RandomSearch <- function(pDesign,nTreatments, nRows, nCols, optDesign, infoDesign, twocombn, check2Combn, weighted, flag, niter){
copyDesign <- pDesign
iter <- 0
recordImp <- rep(0,niter)
while (iter < niter) {
cTreatments <- as.matrix(rbind(which(pDesign==1,arr.ind=T),which(pDesign==2,arr.ind=T),which(pDesign==3,arr.ind=T),which(pDesign==4,arr.ind=T)))
icheck <- sample(1:lcTreatment, 1, replace = T)
siter <- 0
while (siter < 1000) {
irow <- c(as.double(cTreatments[icheck,1]),sample(1:nRows, 1, replace = T))
icol <- c(as.double(cTreatments[icheck,2]),sample(1:nCols, 1, replace = T))
if (pDesign[irow[1], icol[1]] != pDesign[irow[2], icol[2]]) {
check <- updateFormulaeT3(pDesign, infoDesign, nTreatments, nChecks, twocombn, check2Combn, irow[1], icol[1], irow[2], icol[2], weighted,flag)
if (round(check$Aopt*10^10)/10^10 > round(optDesign*10^10)/10^10) {
#print('I am here')
optDesign <- check$Aopt
# We need to swap treatments
iPDesign <- Swapping(pDesign, irow[1], icol[1], irow[2], icol[2])
pDesign <- iPDesign
infoDesign <- check
}
}
siter <- siter + 1
}
recordImp[iter] <- optDesign
iter <- iter + 1
}
mlist <- list(optDesign,pDesign,infoDesign,recordImp)
names(mlist) <- c("optDesign","pDesign","infoDesign","recordImp")
return(mlist)
}
#' Swapping two different plots in a given design
#' @param InvInfoMatrix
#' @param TwoCombn
#' @param Repli
#' @param Weight
#' @param Flag
#' @export
AeffRowColumnDesign <- function(InvInfoMatrix, nTreatments, nChecks, r_a, ...){
# get variables
arguments <- list(...)
TwoCombn <- arguments$TwoCombn
ifelse(!is.null(arguments$TwoCombn),TwoCombn <- arguments$TwoCombn,TwoCombn <- -1)
ifelse(!is.null(arguments$Repli),Repli <- arguments$Repli, Repli <- -1)
ifelse(!is.null(arguments$Weight),Weight <- arguments$Weight, Weight <- -1)
ifelse(!is.null(arguments$Flag),Flag <- arguments$Flag,Flag <- -1)
# set variables
v_bar <- 0
Aopt <- 0
Echeck <- 0
r_square <- r_a*r_a
if ( Weight == 1 ) {
if (Flag=="0") {
rarb <- Repli[TwoCombn[,1]]*Repli[TwoCombn[,2]]
diagElements <- diag(InvInfoMatrix)
AAA <- rarb*(diagElements[TwoCombn[,1]]+diagElements[TwoCombn[,2]])
BBB <- diag(as.vector(Repli))%*%InvInfoMatrix%*%diag(as.vector(Repli))
CCC <- sum(BBB)-sum(diag(BBB))
temp <- (sum(AAA) - CCC)/r_square
v_bar <- temp*2/(nTreatments*(nTreatments-1))
Aopt <- 2/(r_a*v_bar)
Echeck <- Aopt
} else if(Flag=="u") {
nEntries <- nTreatments - nChecks
InvCOUnrep <- InvInfoMatrix[c((nChecks+1):nTreatments),c((nChecks+1):nTreatments)]
temp <- nEntries*sum(diag(InvCOUnrep))-sum(InvCOUnrep)
v_bar <- temp*2/(nEntries*(nEntries-1))
Aopt <- 2/v_bar
Echeck <- Aopt
} else if(Flag=="c") {
InvCO2Checks <- InvInfoMatrix[c(1:nChecks),c(1:nChecks)]
rCheck <- sum(r[1:nChecks])/nChecks
rSquareCheck <- rCheck*rCheck
temp <- 0
for (ii in c(1:(nChecks-1))) {
for(jj in c((ii+1):nChecks)){
v <- (r[ii]*r[jj]/rSquareCheck)*(InvCO2Checks[ii,ii]+InvCO2Checks[jj,jj]-2*InvCO2Checks[ii,jj])
temp <- temp + v
}
}
v_bar <- temp*2/(nChecks*(nChecks-1))
Aopt <- 2/(rCheck*v_bar)
Echeck <- Aopt
}
} else {
temp <- nTreatments*sum(diag(InvInfoMatrix))-sum(InvInfoMatrix)
v_bar <- temp*2/(nTreatments*(nTreatments-1))
Aopt <- 2/(r_a*v_bar)
Echeck <- Aopt
}
#mlist <- list(optDesign,pDesign,infoDesign,recordImp)
#names(mlist) <- c("optDesign","pDesign","infoDesign","recordImp")
return(Echeck)
}
##----------------------------------##
## A place to test your codes
##----------------------------------##
#setwd("D:/user/vthanh/Google Drive/Project Stuttgart/R_programing/UpdateFormular")
#source("D:/user/vthanh/Google Drive/Project Stuttgart/R_programing/A_Eff_RowColumnDesign.r")
#data <- split(read.table("./Designs/Design.txt"), gl(1, 4))
#repli <- 3
#pDesign <- data$`1`
#modelD <- modelMatrix(10, 4, 4, pDesign)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.