Nothing
#' Search for a Cross-Over Design
#'
#' Search for a Cross-Over Design
#'
#' See the vignette of this package for further details.
#'
#' @param s Number of sequences.
#' @param p Number of periods.
#' @param v Number of treatments.
#' @param model Model - one of the following: "Standard additive model" (2),
#' "Second-order carry-over effects" (3), "Full set of interactions" (3),
#' "Self-adjacency model" (3), "Placebo model" (2), "No carry-over into self
#' model" (2), "Treatment decay model" (2), "Proportionality model" (1), "No carry-over effects" (0).
#' The number in parentheses is the number of different efficiency factors that can
#' be specified.
#' @param eff.factor Weights for different efficiency factors. (Not used in the
#' moment.)
#' @param v.rep Integer vector specifying how often each treatment should be
#' assigned (sum must equal s*p).
#' @param balance.s Boolean specifying whether to allocate the treatments as
#' equally as possible to each sequence (can result in loss of efficiency).
#' @param balance.p Boolean specifying whether to allocate the treatments as
#' equally as possible to each period (can result in loss of efficiency).
#' @param verbose Level of verbosity, a number between 0 and 10. The default
#' \code{verbose=0} does not print any output, while \code{verbose=10} prints
#' any available notes.
#' @param model.param List of additional model specific parameters. In the
#' moment these are \code{ppp}, the proportionality parameter for the
#' proportionality model, and \code{placebos}, the number of placebo treatments
#' in the placebo model.
#' @param n \code{n=c(n1,n2)} with \var{n1} the number of hill climbing steps
#' per trial and \var{n2} the number of searches from random start matrices.
#' @param jumps To reduze the possibility of the hill-climbing algorithm to get
#' stuck in local extrema long jumps of distance \var{d} can be performed all
#' \var{k} steps. This can be specified as \code{long.jumps=c(d,k)}. If
#' \var{long.jumps} has only length 1 the default for \var{k} is 50. If after
#' \var{k/2} hill-climbing steps the old design criterion is not enhanced (or
#' at least reached), the algorithm returns to the design from before the jump.
#' @param start.designs A single design or a list of start designs. If missing or to few start
#' designs are specified (with regard to parameter \var{n} which specifies a
#' number of 20 start designs as default) the start designs are generated
#' randomly with the sample function. Alternatively
#' \code{start.designs="catalog"} can be used to take start designs from the
#' catalog to which random designs are added till \var{n2} start designs are at
#' hand.
#' @param contrast Contrast matrix to be optimised. TODO: Example and better
#' explanation for contrast.
#' @param random.subject Should the subject effects be random (\code{random.subject=TRUE})
#' or fixed effects (\code{random.subject=FALSE}).
#' @param correlation Either a correlation matrix for the random subject effects or one
#' of the following character strings: "equicorrelated", "autoregressive"
#' @param rho Parameter for the correlation if parameter \code{correlation} is a character string.
#' @return Returns the design as an integer matrix.
#' @author Kornelius Rohmeyer \email{rohmeyer@@small-projects.de}
#' @references John, J. A., Russell, K. G., & Whitaker, D. (2004). CrossOver:
#' an algorithm for the construction of efficient cross-over designs.
#' Statistics in medicine, 23(17), 2645-2658.
#' @keywords misc
#' @examples
#'
#' \dontrun{
#' x <- searchCrossOverDesign(s=9, p=5, v=4, model=4)
#'
#' jumps <- c(10000, 200) # Do a long jump (10000 changes) every 200 steps
#' n <- c(1000, 5) # Do 5 trials with 1000 steps in each trial
#' result <- searchCrossOverDesign(s=9, p=5, v=4, model=4, jumps=jumps, n=n)
#' plot(result)
#' }
#'
#'
#' @export searchCrossOverDesign
searchCrossOverDesign <- function(s, p, v, model="Standard additive model", eff.factor=1,
v.rep, balance.s=FALSE, balance.p=FALSE, verbose=0, model.param=list(),
n=c(5000, 20), jumps=c(5, 50), start.designs, random.subject=FALSE, contrast, correlation=NULL, rho=0) {
if (!exists(".summary_table", envir=Crossover.env)) {
assign(".summary_table", buildSummaryTable(), envir=Crossover.env)
}
#seed <<- .Random.seed #TODO Do not forget to remove this after testing! :)
start.time <- proc.time()
if (length(n)==1) {
if (missing(start.designs)) { n <- c(n, 20) } else { n <- c(n, length(start.designs)) }
}
if (length(jumps)==1) jumps <- c(jumps, 50)
if (jumps[2]==0) stop("The second component of 'jumps' must be a positive integer.")
model <- getModelNr(model)
H <- do.call( linkMatrix, c(list(model=model, v=v), model.param) )
interchange <- TRUE
if (missing(v.rep)) {
v.rep <- rep((s*p) %/% v, v) + c(rep(1, (s*p) %% v), rep(0, v-((s*p) %% v)))
interchange <- FALSE
} else if (sum(v.rep)!=s*p) { # TODO Feature: Allow NA or sum(v.rep)<s*p
stop("The sum of argument v.rep must equal s times p.")
}
if (balance.s || balance.p) interchange <- TRUE
if (balance.s && balance.p) stop("Balancing sequences AND periods simultaneously is a heavy restriction and not supported (yet?).")
if (missing(contrast)) {
Csub <- contrMat(n=rep(1, v), type="Tukey")
class(Csub) <- "matrix" #TODO Package matrix can be improved here (IMO)!
C <- appendZeroColumns(Csub, model, v)
} else {
if (is.matrix(contrast)) {
C <- contrast
} else {
C <- contrMat2(type=contrast, v, model, eff.factor)
}
}
if (!is.null(correlation) && !is.matrix(correlation)) correlation <- corMat(correlation, s=s, p=p, rho=rho)
if (missing(start.designs)) { start.designs <- list() } # In this list we save n[2] random start designs.
if (isTRUE(start.designs %in% c("catalog","catalogue"))) {
st <- get(".summary_table", envir=Crossover.env)
start.designs <- lapply(st[st$t==v & st$p==p & st$s==s,]$dataset, get, envir=Crossover.env)
x <- getStartDesigns(s=s, p=p, v=v)
if (length(x)!=0) start.designs <- c(start.designs, x)
}
if (!is.list(start.designs)) {
start.designs <- list(start.designs)
}
i <- length(start.designs) + 1
while (i <= n[2]) {
start.designs[[i]] <- randomDesign(s, p, v, v.rep, balance.s=balance.s, balance.p=balance.p, model=model, C=C)
i <- i + 1
}
if (length(start.designs)!=n[2]) { warning(paste("More start designs specified than number of search runs. Search runs are increased to ", length(start.designs), ".", sep="")) }
CC <- t(C) %*% C
if (model==8) { # Second-order carry-over effects
r <- c(rep(s/v, v), rep((p-1)*s/v^2, v^2), rep((p-2)*s/v^3, v^3))
} else if (model==9) {
r <- c(rep(s*p/v, v))
} else {
r <- c(rep(s/v, v), rep((p-1)*s/v^2, v^2))
}
S2 <- sum(diag(ginv(t(H) %*% diag(r) %*% H) %*% CC))
result <- .Call( "searchCOD", as.integer(s), as.integer(p), as.integer(v),
start.designs, H, C, model,
v.rep, balance.s, balance.p, verbose,
as.integer(n), as.integer(jumps), S2, TRUE, random.subject, correlation, interchange, PACKAGE = "Crossover")
design <- result$design
if (!estimable(design, v=v, model=model, C=C)) {
stop("Something went wrong. Specified contrasts are not estimable with this design.")
}
time <- proc.time()-start.time
class(time) <- NULL
#varTrtPair <- paste(capture.output(print(general.carryover(t(design), model=model))), collapse = "\n")
#return(list(design=design, varTrtPair=varTrtPair, eff=eff, search=, time=))
return(new("CrossoverSearchResult", design=new("CrossoverDesign", design, v=v, model=model, description="Found by search algorithm"), startDesigns=start.designs, eff=result$eff,
search=list(n=n, jumps=jumps), model=model, time=time, misc=list(designs=result$designs)))
}
getS1 <- function(design, v, model, C, randomS=FALSE, verbose=0) {
if(missing(C)) {
Csub <- contrMat(n=rep(1, v), type="Tukey")
class(Csub) <- "matrix" #TODO Package matrix can be improved here (IMO)!
C <- appendZeroColumns(Csub, model, v)
}
linkM <- linkMatrix(model, v)
return(.Call( "getS12R", design, v, model, linkM, C, randomS))
}
appendZeroColumns <- function(Csub, model, v) {
model <- getModelNr(model)
if (model %in% c(2,8)) {
C <- as.matrix(cbind(Csub, matrix(0, dim(Csub)[1], 2*v)))
} else if (model %in% c(3,9)) {
C <- Csub
} else if (model == 7) {
C <- as.matrix(cbind(Csub,matrix(0,dim(Csub)[1], v+v*v)))
} else if (model %in% c(1,4,5,6) ) {
C <- as.matrix(cbind(Csub, matrix(0,dim(Csub)[1], v)))
}
return(C)
}
getTrtPair <- function(design, model=1) {
gco <- general.carryover(t(design), model=model)
return(triu(gco$Var.trt.pair))
}
getValues <- function(design, model=1, C, v, ppp=0.5, placebos=1) {
model <- getModelNr(model)
if (missing(C)) {
Csub <- contrMat(n=rep(1, v), type="Tukey")
class(Csub) <- "matrix" #TODO Package matrix can be improved here (IMO)!
C <- appendZeroColumns(Csub, model, v)
CC <- t(C) %*% C
}
rcDesign <- rcd(design, v, model=model)
Ar <- infMatrix(rcDesign, v, model=model)
H <- linkMatrix(model, v, ppp=ppp, placebos=placebos)
return(diag(C %*% ginv(t(H) %*% Ar %*% H) %*% t(C)))
}
getDesignMatrix <- function(design, v) {
H <- linkMatrix(model="Standard additive model", v)
rcDesign <- rcd(design, v=v, model=1)
Xr <- rcdMatrix(rcDesign, v, model=1)
return(Xr %*% H)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.