R/Block3Designs.R

In nhavt/Block3AugDesigns:

#' 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
##----------------------------------##