Nothing
R/seqblock.R
In blockTools: Block, Assign, and Diagnose Potential Interference in Randomized Experiments
Defines functions
seqblock
Documented in
seqblock
#' @title Sequential assignment of unit(s) into experimental conditions using covariates
#'
#' @description
#' Sequentially assign units into experimental conditions. Blocking begins by creating a measure of multivariate distance between a \emph{current} unit and one or multiple \emph{prior}, already-assigned unit(s). Then, average distance between current unit and each treatment condition is calculated, and random assignment is biased toward conditions more dissimilar to current unit. Argument values can be specified either as argument to the function, or via a query. The query directly asks the user to identify the blocking variables and to input, one-by-one, each unit's variable values.
#'
#' @details
#' The \code{seqblock} function's code is primarily divided into two parts: the first half deals with instances, in which the unit being assigned is the first unit in a given study to receive an assignment; the second half addresses subsequent units that are assigned after at least one first assignment has already been made. If the \code{object} argument is left as \code{NULL}, the function will run the first half; if the \code{object} argument is specified, the second part will be executed. When \code{object = NULL}, the researcher has no past file from which to pull variable names and past data; this corresponds to the case when the unit being assigned is the first one. If the researcher does specify \code{object}, it implies the user is drawing data from a past file, which means this is not the first unit in the study to be assigned to a treatment.
#'
#' However, the function can be called for subsequent units even when \code{object} is not specified. By setting \code{query = TRUE}, the console will ask the researcher whether this is the first unit to be assigned in the study. Based on the researcher's response, it will decide which part of the code to run.
#'
#' If the \code{object} and \code{file.name} arguments are set to the same value, then \code{seqblock} overwrites the specified file with a new file, which now contains both the previously-assigned units and the newly-assigned unit. To create a new file when a new unit is assigned, use a new \code{file.name}.
#'
#' The \code{single} algorithm (see \code{exact.alg} in the Arguments section above) creates a variable that has a unique level for every possible combination of the exact variables. As an example, say that there were 3 exact blocking variables: \emph{party} (Democrat, Republican); \emph{region} (North, South); and \emph{education} (HS, NHS). The \code{single} algorithm creates one level for units with the following values: Democrat-North-HS. It would create another level for Democrat-North-NHS; a third level for Republican-North-HS; and so forth, until every possible combination of these 3 variables has its own level. Thus if there are \eqn{k} exact blocking variables and each exact blocking variable has \eqn{q_{i}} values it can take on, then there are a total of \eqn{\prod_{1}^{k} q_{i}} levels created.
#'
#' The \code{distance = "mcd"} and \code{distance = "mve"} options call \code{cov.rob} to calculate measures of multivariate spread robust to outliers. The \code{distance = "mcd"} option calculates the Minimum Covariance Determinant estimate (Rousseeuw 1985); the \code{distance = "mve"} option calculates the Minimum Volume Ellipsoid estimate (Rousseeuw and van Zomeren 1990). When \code{distance = "mcd"}, the interquartile range on blocking variables should not be zero. The \code{distance = "euclidean"} option calculates the Euclidean distance between the new unit and the previously-assigned units. The default \code{distance = "mahalanobis"} option calculates the Mahalanobis distance.
#'
#' @param object a character string giving the file name of a \code{.RData} file containing a list output from the \code{seqblock} function which contains at least one previously assigned unit.
#' @param id.vars a string or vector of strings specifying the name of the identifying variable(s); if \code{query = FALSE} and the object argument is not given, then the \code{id.vars} argument is required.
#' @param id.vals a vector of ID values for every unit being assigned to a treatment group; those are corresponding to the \code{id.vars}.
#' @param exact.vars a string or vector of strings containing the names of each of the exact blocking variables.
#' @param exact.vals a vector containing the unit's values on each of the exact blocking variables.
#' @param exact.restr a list object containing the restricted values that the exact blocking variables can take on. Thus the first element of \code{exact.restr} is a vector containing all of the possible values that the first exact blocking variable (see \code{exact.vars} above) can take on; the second element is a vector containing all of the possible values for the second exact blocking variable; and so on.
#' @param exact.alg a string specifying the blocking algorithm. Currently the only acceptable value is \code{"single"}. This algorithm creates a variable with a unique level for every possible combination of the values in all of the exact variables. See Details section below.
#' @param covar.vars a string or vector of strings containing the names of each of the non-exact blocking variables.
#' @param covar.vals a vector containing the unit's values on each of the non-exact blocking variables.
#' @param covar.restr a list object containing the restricted values that the non-exact blocking variables can take on. Thus the first element of \code{covar.restr} is a vector containing all of the possible values that the first non-exact blocking variable (see \code{covar.vars} above) can take on; the second element is a vector containing all of the possible values for the second non-exact blocking variable; and so on.
#' @param covars.ord a string or vector of strings containing the name of the non-exact blocking variables ordered so that the highest priority covariate comes first, followed by the second highest priority covariate, then the third, etc.
#' @param n.tr the number of treatment groups. If not specified, this defaults to \code{n.tr = 2}.
#' @param tr.names a string or vector of strings containing the names of the different treatment groups.
#' @param assg.prob a numeric vector containing the probabilities that a unit will be assigned to the treatment groups; this vector should sum to 1.
#' @param seed an optional integer value for the random seed, which is used when assigning units to treatment groups.
#' @param seed.dist an optional integer value for the random seed set in \code{cov.rob}, used to calculate measures of the variance-covariance matrix robust to outliers.
#' @param assg.prob.stat a string specifying which assignment probability summary statistic to use; valid values are \code{mean}, \code{median}, and \code{trimmean}. If not specified, this defaults to \code{assg.prob.stat = "mean"}.
#' @param trim a numeric value specifying what proportion of the observations are to be dropped from each tail when the assignment probability summary statistic (\code{assg.prob.stat}) is set equal to \code{trimmean}. Blocks on each tail of the distribution are dropped before the mean is calculated. If not specified, this defaults to \code{trim = 0.1}.
#' @param assg.prob.method a string specifying which algorithm should be used when assigning treatment probabilities. Acceptable values are \code{ktimes}, \code{fixed}, \code{prop}, \code{prop2}, and \code{wprop}. If not specified, this defaults to \code{assg.prob.method = "ktimes"}.
#' @param assg.prob.kfac a numeric value for \code{k}, the factor by which the most likely experimental condition will be multiplied relative to the other conditions. If not specified, this defaults to \code{assg.prob.kfac = 2}.
#' @param distance a string specifying how the multivariate distance used for blocking covariates are calculated. If not specified, this defaults to \code{distance = "mahalanobis"}.
#' @param file.name a string containing the name of the file that one would like the output to be written to. Ideally this file name should have the extension .RData.
#' @param query a logical stating whether the console should ask the user questions to input the data and assign a treatment condition. If not specified, this defaults to \code{query = FALSE}.
#' @param verbose a logical stating whether the function should print the name of the output file, the current working directory, the treatment group that the most recent unit was assigned to, and the dataframe \code{x} returned by the function as part of the \code{bdata} list. If not specified, this defaults to \code{verbose = TRUE}.
#' @param \dots additional arguments.
#'
#' @return
#' A list (called \code{bdata}) with elements
#' \itemize{
#' \item \strong{x}: a dataframe containing the names and values for the different ID and blocking variables, as well as each unit's initial treatment assignment.
#' \item \strong{nid}: a string or vector of strings containing the name(s) of the ID variable(s).
#' \item \strong{nex}: a string or vector of strings containing the name(s) of the exact blocking variable(s).
#' \item \strong{ncv}: a string or vector of strings containing the name(s) of the non-exact blocking variable(s).
#' \item \strong{rex}: a list of the restricted values of the exact blocking variables.
#' \item \strong{rcv}: a list of the restricted values of the non-exact blocking variables.
#' \item \strong{ocv}: a vector of the order of the non-exact blocking variables.
#' \item \strong{trn}: a string or vector of strings containing the name(s) of the different treatment groups.
#' \item \strong{apstat}: a string specifying the assignment probability summary statistic that was used.
#' \item \strong{mtrim}: a numeric value specifying the proportion of observations to be dropped when the assignment probability statistic takes on the value \code{"trimmean"}.
#' \item \strong{apmeth}: a string specifying the assignment probability algorithm that was used.
#' \item \strong{kfac}: the assignment probability \emph{kfactor}; see \emph{assg.prob.kfac} in the Arguments section above.
#' \item \strong{assgpr}: a vector of assignment probabilities to each treatment group.
#' \item \strong{distance}: a string specifying how the multivariate distance used for blocking was calculated.
#' \item \strong{trd}: a list with the length equal to the number of previously assigned treatment conditions; each object in the list contains a vector of the distance between each unit in one treatment group and the new unit. This will be \code{NULL} when there are no non-exact blocking variables.
#' \item \strong{tr.sort}: a string vector of treatment conditions, sorted from the largest to the smallest. Set to \code{NULL} when there are no non-exact blocking variables.
#' \item \strong{p}: a vector of assignment probabilities to each treatment group used in assigning a treatment condition to the new unit.
#' \item \strong{distance}: a string specifying how the multivariate distance used for blocking is calculated
#' \item \strong{trcount}: a table containing the counts for each experimental/treatment conditions.
#' \item \strong{datetime}: the date and time at which each unit was assigned their treatment group.
#' \item \strong{orig}: a dataframe containing the names and values for the different id and blocking variables, as well as each unit's treatment assignment.
#' }
#'
#' @examples
#' # Assign first unit (assume a 25 year old member of the Republican Party)
#' # to a treatment group. Save the results in file "sdata.RData":
#' # seqblock(query = FALSE, id.vars = "ID", id.vals = 001, exact.vars = "party",
#' # exact.vals = "Republican", covar.vars = "age", covar.vals = 25, file.name = "sdata.RData")
#'
#' # Assign next unit (age 30, Democratic Party):
#' # seqblock(query = FALSE, object = "sdata.RData", id.vals = 002, exact.vals = "Democrat",
#' # covar.vars = "age", covar.vals = 30, file.name = "sdata.RData")
#'
#' @references
#' Moore, Ryan T. and Sally A. Moore. 2013. "Blocking for Sequential Political Experiments."
#' \emph{Political Analysis} 21(4):507-523.
#'
#' Moore, Ryan T. 2012. "Multivariate Continuous Blocking to Improve Political Science Experiments." \emph{Political Analysis} 20(4):460-479.
#'
#' Rousseeuw, Peter J. 1985. "Multivariate Estimation with High Breakdown Point". \emph{Mathematical Statistics and Applications} 8:283-297.
#'
#' Rousseeuw, Peter J. and Bert C. van Zomeren. 1990. "Unmasking Multivariate Outliers and Leverage Points". \emph{Journal of the American Statistical Association} 85(411):633-639.
#'
#' @seealso \code{\link{assignment}}, \code{\link{block}}
#'
#' @author Ryan T. Moore \email{rtm@american.edu}, Tommy Carroll \email{tcarroll22@wustl.edu}, Jonathan Homola \email{homola@wustl.edu} and Jeong Hyun Kim \email{jeonghyun.kim@wustl.edu}
#'
#' @keywords design multivariate
#'
#' @export
seqblock <- function(object = NULL, id.vars, id.vals, exact.vars = NULL, exact.vals = NULL, exact.restr = NULL, exact.alg = "single", covar.vars = NULL, covar.vals = NULL, covar.restr = NULL, covars.ord = NULL, n.tr = 2, tr.names = NULL, assg.prob = NULL, seed = NULL, seed.dist, assg.prob.stat = NULL, trim = NULL, assg.prob.method = NULL, assg.prob.kfac = NULL, distance = NULL, file.name = NULL, query = FALSE, verbose = TRUE, ...){
if(is.null(object)){
if(query == TRUE){
init <- readline("Is this the first unit that is being assigned in this experiment? [y/n] ")
if(!(substr(init, 1, 1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
init <- "y"
}else(stop("The seqblock function requires the ''object'' argument to be set to a valid file name for subsequent unit assignment."))
}
if(substr(init, 1, 1) == "n"){
object <- readline("Enter the name of the input file without quotation marks. [E.g., sbout1.RData] ")
}
}
}
if(is.null(object)){ ## If assigning first unit:
if(query == TRUE){
## ID
lid <- readline("How many identification variables are there? ")
id.vars <- id.vals <- NULL
for(ii in 1:lid){ # For each identification variables,
## tmp stores the name of ID variable:
tmp <- readline(cat("Enter the name of ID variable ", ii, " without quotation marks.", sep=""))
## Add the ii-th name of ID variable to a vector id.vars:
id.vars <- append(id.vars, tmp)
## tmp2 stores the value of each ID variable:
tmp2 <- readline(cat("Enter the value of '", tmp, "'. ", sep=""))
## Add the ii-th value of ID variable to a vector id.vals
id.vals <- append(id.vals, tmp2)
} ## End loop over ID variables
nnn <- readline("Should the ID values be numeric? [y/n] ")
if(!(substr(nnn, 1, 1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline("Should the ID values be numeric? [y/n] ")
}
}
if(substr(nnn,1,1) != "n"){ # If the respondent answers that the ID values should be numeric,
id.vals <- as.numeric(id.vals) # change id.vals to numerics.
}
if(sum(is.na(id.vals)) > 0){ # If there is at least one missing value in id.vals,
tmp <- readline(cat("Warning: ID value(s)", which(is.na(id.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
if(!(substr(tmp, 1, 1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp <- readline(cat("Warning: ID value(s)", which(is.na(id.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
}
}
if(substr(tmp, 1, 1) == "n"){
stop()
}
}
## Specifying EXACT blocking variables:
lex <- readline("`Exact blocking variables' are those on which you require covariate values to be identical. How many exact blocking variables are there? ")
exact.vars <- exact.vals <- exact.restr <- NULL
if(lex > 0){ ## If there are exact blocking variables
for(ii in 1:lex){ ## Then, for each exact blocking variable
## tmp stores the name of exact blocking variable:
tmp <- readline(cat("Enter the name of exact block variable ", ii, " without quotation marks.", sep=""))
exact.vars <- append(exact.vars, tmp) ## then add to the vector exact.vars
tmp2 <- readline(cat("Enter the value of '", tmp, "'. ", sep=""))
tmp3 <- readline(cat("Should '", tmp, "' be restricted to certain values? [n/y] ", sep=""))
if(!(substr(tmp3,1,1) %in% c("y", "n"))){ ## if the answer is neither yes nor no,
ddd <- readline("The default is 'no'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp3 <- readline(cat("Should '", tmp, "' be restricted to certain values? [n/y] "))
}
} ## Close 'if'
if(substr(tmp3, 1, 1) == "y"){ ## if the ii-th exact variable should be restricted to certain values, tmp4 stores the values that variable can take that the respondent inputs.
tmp4 <- readline(cat("How many values can '", tmp, "' take? ", sep=""))
## if the value is less than one, execution is halted.
if(tmp4 < 1){
stop("Restricted variables must have at least one valid value. ")
}
exact.restr <- as.list(exact.restr)
exact.restr[[tmp]] <- NA
for(rrr in 1:tmp4){ ## for each value that the variable tmp can take,
tmp5 <- readline(cat("Enter possible value for '", tmp, "' number ", rrr, sep=""))
exact.restr[[tmp]] <- append(exact.restr[[tmp]], tmp5) ## the object named tmp in the list exact.restr will be a vector of the values that the variable can take.
} ## End loop over values 'tmp' can take
} ## Close 'if'
## if the ii-th exact variable should NOT be restricted to certain values, exact.vals simply contains the values of exact blocking variables.
exact.vals <- append(exact.vals, tmp2)
} ## Close loop over exact blocking variables 'ii in 1:lex'
nnn <- readline("Should the exact block values be numeric? [y/n] ")
if(!(substr(nnn, 1, 1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline("Should the exact block values be numeric? [y/n] ")
}
}
## If the answer is yes (the exact block values should be numeric),
if(substr(nnn,1,1) != "n"){
## change values in exact.vals into numerics:
exact.vals <- as.numeric(exact.vals)
## If there is a restriction on the values of exact variable,
if(length(exact.restr) > 0){
for(rrr in 1:length(exact.restr)){
## Each element of those values are changed to numerics.
exact.restr[[rrr]] <- as.numeric(exact.restr[[rrr]])
}
}
}
if(sum(is.na(exact.vals)) > 0){
tmp <- readline(cat("Warning: Exact blocking value(s)", which(is.na(exact.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
if(!(substr(tmp, 1, 1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp <- readline(cat("Warning: Exact blocking value(s)", which(is.na(exact.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
}
}
if(substr(tmp,1,1) == "n"){
stop()
}
} ## End confirm NA values for exact blocking vars
## CHECK that each exact.val is in exact.restr here.
} ## End 'if(lex > 0)'
## COVARIATES
lcv <- readline("How many other blocking variables are there? ")
covar.vars <- covar.vals <- covar.restr <- NULL
if(lcv > 0){ ## If there is at least one non-exact blocking variable,
for(ii in 1:lcv){ ## For each of those variables,
tmp <- readline(cat("Enter the name of blocking variable ", ii, " without quotation marks.", sep=""))
## Save its name to a vector 'covar.vars'
covar.vars <- append(covar.vars, tmp)
tmp2 <- readline(cat("Enter the value of '", tmp, "'. ", sep=""))
tmp3 <- readline(cat("Should '", tmp, "' be restricted to certain values? [n/y] ", sep=""))
if(!(substr(tmp3, 1, 1) %in% c("y", "n"))){
ddd <- readline("The default is 'no'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp3 <- readline(cat("Should '", tmp, "' be restricted to certain values? [n/y] ", sep=""))
}
} ## close if
## If the value of the blocking variable should be restricted to certain values,
if(substr(tmp3,1,1) == "y"){
tmp4 <- readline(cat("How many values can '", tmp, "' take? ", sep=""))
if(tmp4 < 1){
stop("Restricted variables must have at least one valid value. ")
}
covar.restr <- as.list(covar.restr)
covar.restr[[tmp]] <- NA
for(rrr in 1:tmp4){
tmp5 <- readline(cat("Enter possible value for '", tmp, "' number ", rrr, sep=""))
# Make a list covar.restr in which each object represents restricted values of blocking variable that should be restricted to certain values.
covar.restr[[tmp]] <- append(covar.restr[[tmp]], tmp5)
}
}
## If the values should NOT be restricted to certain values, covar.vals contains the values of the blocking variable:
covar.vals <- append(covar.vals, tmp2)
} ## End loop over blocking vars 'ii in 1:lcv'
nnn <- readline("Should the blocking values be numeric? [y/n] ")
if(!(substr(nnn,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline("Should the blocking values be numeric? [y/n] ")
}
}
if(substr(nnn,1,1) != "n"){
covar.vals <- as.numeric(covar.vals)
}else{
stop("The blocking variables are stored as character strings.\nCurrent implementation does not support distance calculations for character strings.")
} ## (Change to warning() when categorical distances incorporated)
if(sum(is.na(covar.vals)) > 0){
tmp <- readline(cat("Warning: Blocking value(s)", which(is.na(covar.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
if(!(substr(tmp,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp <- readline(cat("Warning: Blocking value(s)", which(is.na(covar.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
}
}
if(substr(tmp,1,1) == "n"){
stop()
}
}
## Covariates' order
nnn <- readline("Would you like to specify the order in which blocking variables will be used while (number units considered) < (number covariates + 2) ? (If not, covariates will be added as possible in the order in which they were entered.) [n/y] ")
if(!(substr(nnn,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `no'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline(cat("Would you like to specify the order in which blocking variables will be used while (number units considered) < (number covariates +2) ? [n/y] "))
}
}
if(substr(nnn,1,1) == "y"){
for(ii in 1:lcv){
tmp <- readline(cat("Enter the name of blocking variable prioritized number ", ii, ".", sep=""))
# if(!(tmp %in% names(x)[(lid+lex+1):(lid+lex+lcv)])){
if(!(tmp %in% covar.vars)){
tmp <- readline(cat("`", tmp, "' is not a blocking variable. Please try again. Enter the name of blocking variable prioritized number ", ii, ".", sep=""))
}
if(tmp %in% covars.ord){
tmp <- readline(cat("`", tmp, "' has already been prioritized. Please try again. Enter the name of blocking variable prioritized number ", ii, ".", sep=""))
}
covars.ord <- append(covars.ord, tmp)
}
}
} ## End 'if(lcv > 0)'
## Number of treatment conditions:
n.tr <- readline("How many experimental/treatment conditions are there? ")
n.tr <- as.numeric(n.tr)
## Names of treatment conditions:
nnn <- readline(cat("Would you like to specify the names of the experimental conditions? [y/n] [If not, `Treatment 1', ..., `Treatment ", n.tr, "', will be used] ", sep=""))
if(substr(nnn,1,1) != "n"){
tr.names <- NULL
for(ii in 1:n.tr){
tmp <- readline(cat("Enter the name of experimental condition ", ii, " without quotation marks.", sep=""))
tr.names <- append(tr.names, tmp)
}
}
## Initial condition assignment probabilities:
nnn <- readline("Would you like to specify the initial assignment probabilities? [y/n] [If not, all treatment assignment probabilities for the first unit will be equal.] ")
if(!(substr(nnn,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline("Would you like to specify the initial assignment probabilities? [y/n] [If not, all treatment assignment probabilities for the first unit will be equal.] ")
}
}
if(substr(nnn,1,1) != "n"){
assg.prob <- NULL
for(ii in 1:n.tr){
tmp <- readline(cat("Enter the assignment probability of experimental condition ", ii, ".", sep=""))
assg.prob <- append(assg.prob, tmp)
}
assg.prob <- as.numeric(assg.prob)
}
## Set random seed:
nnn <- readline("Would you like to specify a random seed for the initial experimental condition assignment? [y/n] ")
if(!(substr(nnn,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline("Would you like to specify a random seed for the initial experimental condition assignment? [y/n] ")
}
}
if(substr(nnn,1,1) != "n"){
seed <- as.numeric(readline("Enter the random seed. "))
}
## ASSIGNMENT PROBABILITY STATISTIC:
nnn <- readline("Would you like to specify the assignment probability summary statistic? [n/y] [If not, `mean' will be used.] ")
if(substr(nnn,1,1) == "y"){
assg.prob.stat <- readline("Enter the summary statistic name. [mean, median, trimmean] ")
if(assg.prob.stat == "trimmean"){
ooo <- readline("Specify the proportion of blocks to be dropped from each end of the distribution before the mean is calculated. [Defaults to 0.1] ")
if(substr(ooo,1,1) == ""){
trim <- 0.1
}else{
trim <- ooo
}
}
}else{
assg.prob.stat <- "mean"
}
## ASSIGNMENT PROBABILITY METHOD:
nnn <- readline("Would you like to specify the assignment probability algorithm? [n/y] [If not, `ktimes' with k=2 will be used.] ")
if(substr(nnn,1,1) == "y"){
assg.prob.method <- readline("Enter the algorithm name. [ktimes, fixed, prop, prop2, wprop] ")
if(assg.prob.method == "ktimes"){
assg.prob.kfac <- as.numeric(readline("Enter the value of k, the factor by which the most likely experimental condition will be multiplied, relative to the other conditions. "))
}
}else{
assg.prob.method <- "ktimes"
}
if(assg.prob.method == "fixed"){
assgpr <- readline("Enter 1st probability. ")
for(pc in 2:n.tr){
assg.prob <- append(assg.prob, readline(cat("Enter probablity number ", pc, ". ", sep="")))
}
assg.prob <- as.numeric(assg.prob)
}
## MULTIVARIATE DISTANCE CALCULATION:
nnn <- readline("Would you like to specify how the multivariate distance used for blocking is calculated? [n/y] [If not, `mahalanobis' will be used.] ")
if(substr(nnn,1,1) == "y"){
distance <- readline("Enter the calculation type. [mahalanobis, mcd, mve, euclidean] ")
}else{
distance <- "mahalanobis"
}
## OUTPUT FILE NAME:
nnn <- readline("Would you like to specify the output data file name? [y/n] [If not, `sbout1.RData' will be used.] ")
if(!(substr(nnn,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline("Would you like to specify the output data file name? [y/n] [If not, `sbout1.RData' will be used.] ")
}
}
if(substr(nnn,1,1) != "n"){
file.name <- readline("Enter the file name. Note: seqblock2k() expects a file of type `.RData'. ")
}
} ## END (query == T)
## create variables
lid <- length(id.vars)
lex <- length(exact.vars)
lcv <- length(covar.vars)
## create vector of treatment names
if(is.null(tr.names)){
tr.names <- paste("Treatment", 1:n.tr)
}
## create kfac, if null
if(!(is.null(assg.prob.kfac))){
kfac <- assg.prob.kfac
} else {kfac <- assg.prob.kfac <- 2}
## create assg.prob.stat, if null
if(!(is.null(assg.prob.stat))){
assg.prob.stat <- assg.prob.stat
} else {assg.prob.stat <- "mean"}
## create trim, if null
if(!(is.null(trim))){
trim <- trim
} else {trim <- 0.1}
## create assg.prob.method, if null
if(!(is.null(assg.prob.method))){
assg.prob.method <- assg.prob.method
} else {assg.prob.method <- "ktimes"}
## create distance, if null
if(is.null(distance)){
distance <- "mahalanobis"
}
## create vector of (equal) initial treatment probabilities
if(is.null(assg.prob)){
assg.prob <- rep(1/n.tr, n.tr)
}
## check that assg.prob sum to 1
if(sum(assg.prob)!=1){
stop("Initial assignment probabilities do not sum to 1.\nRespecify and try again.")
}
## check that if exact named, then supplied
if((!is.null(exact.vars)) && is.null(exact.vals)){
stop("Exact (grouped) blocking requested, but no values given.\nSpecify exact block values and try again.")
}
## check that exact vars and vals are same length
if(lex != length(exact.vals)){
stop(cat("Number of exact blocking variables (", lex, ") does not equal number of values (", length(exact.vals), ").\nSpecify exact block values and try again.\n", sep=""))
}
## check that if covars named, then supplied
if((!is.null(covar.vars)) && is.null(covar.vals)){
stop("Covariate blocking requested, but no values given.\nSpecify covariate values and try again.")
}
## check that covar vars and vals are same length
if(length(covar.vars) != length(covar.vals)){
stop("Number of covariate blocking variables does not equal number of values.\nSpecify covariate block values and try again.")
}
## Add NA to each existing restriction set,
## then check if any exact vals violate exact.restr
if((lex > 0) && !(is.null(exact.restr))){
for(i in names(exact.restr)){
if(!(NA %in% exact.restr[[i]])){
exact.restr[[i]] <- append(exact.restr[[i]], NA)
}
wexrstr <- which(exact.vars == i)
if(!(exact.vals[wexrstr] %in% exact.restr[[exact.vars[wexrstr]]])){
stop("Exact blocking value number ", wexrstr, " (", exact.vals[wexrstr], ") is not in the set of values to which exact blocking variable `", i, "' is restricted. Respecify either the value, the restriction set, or both as needed." )
}
}
}
## check if any covar vals violate covar.restr
if((lcv > 0) && (!is.null(covar.restr))){
for(i in names(covar.restr)){
if(!(NA %in% covar.restr[[i]])){
covar.restr[[i]] <- append(covar.restr[[i]], NA)
}
wcvrstr <- which(covar.vars == i)
if(!(covar.vals[wcvrstr] %in% covar.restr[[covar.vars[wcvrstr]]])){
stop("Covariate blocking value number ", wcvrstr, " (", covar.vals[wcvrstr], ") is not in the set of values to which covariate blocking variable `", i, "' is restricted. Respecify either the value, the restriction set, or both as needed." )
}
}
}
## set user-assigned seed
if(!is.null(seed)){
set.seed(seed)
}
## Draw initial treatment assignment
init.tr <- sample(tr.names, 1, replace = FALSE, prob = assg.prob)
## Create data frame of unit 1 data
## create id placeholders
id.tmp <- rep(NA, lid)
names(id.tmp) <- id.vars
## create group placeholders
ex.tmp <- rep(NA, lex)
names(ex.tmp) <- exact.vars
## create data frame
x <- unlist(c(id.tmp, ex.tmp, unname(covar.vals), NA))
x <- data.frame(t(x))
colnames(x) <- c(id.vars, exact.vars, covar.vars, "Tr")
## add id vals. NOTE: will be all char or all num
x[1, 1:lid] <- id.vals
## add group vals. NOTE: will be all char or all num
if(!is.null(exact.vars)){
x[1, (lid+1):(lid+lex)] <- exact.vals
}
## add initial treatment assignment
x[1,ncol(x)] <- init.tr
if(is.null(file.name)){
file.name <- "sbout1.RData"
}
if(substr(file.name, start=nchar(file.name)-5, stop=nchar(file.name)) !=".RData"){
warning("Destination file name does not end in .RData.")
}
## Create storage list 'bdata'
bdata <- list()
bdata$x <- x
bdata$nid <- id.vars
bdata$nex <- exact.vars
bdata$ncv <- covar.vars
bdata$rex <- exact.restr ## EXACT restricted values list
if(!is.null(bdata$rex)){
names(bdata$rex) <- exact.vars
}
bdata$rcv <- covar.restr ## BLOCK restricted values list
if(!is.null(bdata$rcv)){
names(bdata$rcv) <- covar.vars
}
bdata$ocv <- covars.ord ## block covariates order
bdata$trn <- tr.names
bdata$apstat <- assg.prob.stat
bdata$mtrim <- trim
bdata$apmeth <- assg.prob.method
bdata$kfac <- kfac ## multiple for method 'ktimes'
bdata$assgpr <- assg.prob
bdata$distance <- distance
bdata$datetime <- date()
bdata$orig <- x
## Write bdata to file
save(bdata, file = file.name)
if(verbose == TRUE){
cat("Unit 1 data stored as file ", file.name, ".\nThe current working directory is ", getwd(), "\n", sep="")
cat("Unit ", nrow(x), " assigned to ", init.tr, ".\n", sep="")
cat("The new data as entered:\n")
print(x)
}
} ## End 'if(is.null(object))'. I.e., end assigning first unit.
## Next case: if object is not NULL (I.e., begin assignment for subsequent unit after the first one)
if(!is.null(object)){
load(object) ## loads, doesn't export to wkspace
## Rename object elements
x <- bdata$x
nid <- bdata$nid ## names of id.vars
nex <- bdata$nex ## names of exact.vars
ncv <- bdata$ncv ## names of covar.vars
rex <- bdata$rex ## EXACT restricted values list
rcv <- bdata$rcv ## BLOCK restricted values list
ocv <- bdata$ocv ## block covariates order
trn <- bdata$trn ## treatment condition names
apstat <- bdata$apstat ## assignment prob statistic
mtrim <- bdata$mtrim ## trim fraction for apmeth=trimmean
apmeth <- bdata$apmeth ## assignment prob method
kfac <- bdata$kfac ## assg multiple for method 'ktimes'
assgpr <- bdata$assgpr ## assg probs for fixed prob
dist <- bdata$distance ## distance calculation method
orig <- bdata$orig ## original unjittered data
n.tr <- length(trn) ## number of treatment conditions
lid <- length(nid) #number of id variables
lex <- length(nex) #number of exact variables
lcv <- length(ncv) #number of covariates
if(!(is.null(exact.restr))){
rex <- exact.restr
}
if(!(is.null(covar.restr))){
rcv <- covar.restr
}
if(!(is.null(covars.ord))){
ocv <- covars.ord
}
if(!(is.null(assg.prob.stat))){
if(all.equal(assg.prob.stat,apstat) != TRUE){
warning(paste("'assg.prob.stat' of previous assignments was ", apstat, ". 'assg.prob.stat' for this assignment is ", assg.prob.stat, ".", sep = ""))
}
apstat <- assg.prob.stat
}
if(!(is.null(trim))){
if(all.equal(trim,mtrim) != TRUE){
warning(paste("'trim' of previous assignments was ", mtrim, ". 'trim' for this assignment is ", trim, ".", sep = ""))
}
mtrim <- trim
}
if(!(is.null(assg.prob.method))){
if(all.equal(assg.prob.method, apmeth) != TRUE){
warning(paste("'assg.prob.method' of previous assignments was ", apmeth, ". 'assg.prob.method' for this assignment is ", assg.prob.method, ".", sep = ""))
}
apmeth <- assg.prob.method
}
if(!(is.null((assg.prob.kfac)))){
if(all.equal(assg.prob.kfac,kfac) != TRUE){
warning(paste("'kfac' of previous assignments was ", kfac, ". 'kfac' for this assignment is ", assg.prob.kfac, ".", sep = ""))
}
kfac <- assg.prob.kfac
}
if(!(is.null(assg.prob))){
assgpr <- assg.prob
}
if(!(is.null((distance)))){
if(all.equal(distance,dist) != TRUE){
warning(paste("Distance option of previous assignments was '", dist, "'. Distance option for this assignment is '", distance, "'.", sep = ""))
}
dist <- distance
}
if(query == TRUE){
## ID
id.vals <- NULL
for(ii in 1:lid){
tmp <- names(x)[ii]
tmp2 <- readline(cat("Enter the value of '", tmp, "'. ", sep=""))
id.vals <- append(id.vals, tmp2)
}
class(id.vals) <- class(x[,lid]) ## lid could be any of 1:lid
if(sum(is.na(id.vals)) > 0){ # If there is at least one missing value in id.vals,
tmp <- readline(cat("Warning: ID value(s)", which(is.na(id.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
if(!(substr(tmp, 1, 1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp <- readline(cat("Warning: ID value(s)", which(is.na(id.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
}
}
if(substr(tmp, 1, 1) == "n"){
stop()
}
}
## Check for duplicated IDs if single id.var
if((lid == 1) && (id.vals %in% x[, nid])){
tmp <- readline(cat("Warning: ID value ", id.vals, " already used in the data. Allow duplication and proceed? [y/n]", sep = ""))
if(substr(tmp, 1, 1) == "n"){
stop()
}
}
## EXACT
exact.vals <- NULL
if(lex > 0){
for(ii in 1:lex){
tmp <- nex[ii]
tmp2 <- readline(cat("Enter the value of '", tmp, "'. ", sep=""))
exact.vals <- append(exact.vals, tmp2)
}
class(exact.vals) <- class(x[,(lid+lex)]) ## could be any of (lid+1):(lid+lex) If the exact blocking variable was set to be numeric in seqblock1(), new exact values should be also numeric. Otherwise, the values of the exact blocking variable will be coerced to NAs.
if(sum(is.na(exact.vals)) > 0){
tmp <- readline(cat("Warning: Exact blocking value(s)", which(is.na(exact.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
if(!(substr(tmp,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp <- readline(cat("Warning: Exact blocking value(s)", which(is.na(exact.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
}
}
if(substr(tmp, 1, 1) == "n"){
stop()
}
}
}
## EXACT ALGORITHM
## "single": if lex > 1, creates unique categories (4 for 2x2, e.g.)
exact.alg <- "single"
## COVARIATES
covar.vals <- NULL
if(lcv > 0){
for(ii in 1:lcv){
tmp <- ncv[ii]
tmp2 <- readline(cat("Enter the value of '", tmp, "'. ", sep=""))
covar.vals <- append(covar.vals, tmp2)
}
class(covar.vals) <- class(x[,lid+lex+lcv]) ## could be any of (lid+lex+1):(lid+lex+lcv)
if(sum(is.na(covar.vals)) > 0){
tmp <- readline(cat("Warning: Blocking value(s)", which(is.na(covar.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
if(!(substr(tmp,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp <- readline(cat("Warning: Blocking value(s)", which(is.na(covar.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
}
}
if(substr(tmp,1,1) == "n"){
stop()
}
}
}
## SET RANDOM SEED:
nnn <- readline("Would you like to specify a random seed for this experimental condition assignment? [y/n] ")
if(!(substr(nnn,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline(cat("Would you like to specify a random seed for this experimental condition assignment? [y/n] "))
}
}
if(substr(nnn,1,1) != "n"){
seed <- as.numeric(readline("Enter the random seed. "))
}
## OUTPUT FILE NAME:
nnn <- readline("Would you like to specify the output data file name? [y/n] [If not, `sbout2k.RData' will be used.] ")
if(!(substr(nnn,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline(cat("Would you like to specify the output data file name? [y/n] [If not, `sbout2k.RData' will be used.] "))
}
}
if(substr(nnn,1,1) != "n"){
file.name <- readline("Enter the file name. Note: seqblock2k() expects a file of type `.RData'. ")
}else{
file.name <- NULL
}
} ## Close query==T
## Warn if duplicated ID values if single ID var
if((lid == 1) && (id.vals %in% x[, nid])){
warning("Warning: ID value ", id.vals, " already used in the data.")
}
## check # exact vars vs. # exact vals
if(length(nex) != length(exact.vals)){
stop(cat("Number of exact blocking variables (", lex, ") does not equal number of values (", length(exact.vals), ").\nSpecify exact block values and try again.\n", sep=""))
}
## check if any exact val violates exact.restr
if((lex > 0) && !(is.null(rex))){
for(i in names(rex)){
if(!(NA %in% rex[[i]])){
rex[[i]] <- append(rex[[i]], NA)
}
wexrstr <- which(nex == i)
if(!(exact.vals[wexrstr] %in% rex[[nex[wexrstr]]])){
tmp <- readline(cat("Exact blocking value number ", wexrstr, " (", exact.vals[wexrstr], ") is not in the set of values to which exact blocking variable `", i, "' is restricted. Would you like to add ", exact.vals[wexrstr], " to the set? [y/n]", sep=""))
if(substr(tmp,1,1) != "n"){
rex[[i]] <- append(rex[[i]], exact.vals[wexrstr])
}else{
stop("Exact blocking value number ", wexrstr, " (", exact.vals[wexrstr], ") is not in the set of values to which exact blocking variable `", nex[wexrstr], "' is restricted. Respecify either the value, the restriction set, or both as needed." )
}
}
}
}
## check if any block val violates covar.restr
if((lcv > 0) && (!is.null(rcv))){ for(i in names(rcv)){
if(!(NA %in% rcv[[i]])){
rcv[[i]] <- append(rcv[[i]], NA)
}
wcvrstr <- which(ncv == i)
if(!(covar.vals[wcvrstr] %in% rcv[[ncv[wcvrstr]]])){
tmp <- readline(cat("Covariate blocking value number ", wcvrstr, " (", covar.vals[wcvrstr], ") is not in the set of values to which covariate blocking variable `", i, "' is restricted. Would you like to add ", covar.vals[wcvrstr], " to the set? [y/n]", sep=""))
if(substr(tmp,1,1) != "n"){
rcv[[i]] <- append(rcv[[i]], covar.vals[wcvrstr])
}else{
stop("Covariate blocking value number ", wcvrstr, " (", covar.vals[wcvrstr], ") is not in the set of values to which covariate blocking variable `", ncv[wcvrstr], "' is restricted. Respecify either the value, the restriction set, or both as needed.")
}
}
}
}
if(!(length(ocv) == length(unique(ocv)))){
stop("Prioritization of covariates for (number units considered) < (number covariates + 2) includes duplicate variable names. Respecify 'covars.ord'.")
}
if(lex > 0){
if(exact.alg == "single"){
x.ex <- x[apply(t(t(x[,(lid+1):(lid+lex)])==exact.vals), 1, sum) == lex,]
}
}
## When no exact matches, use all previous data:
if(is.null(exact.vals) || (nrow(x.ex)==0)){
x.ex <- x
}
prev <- nrow(x.ex)
## Create data frame of unit 1 data
## create id placeholders
id.tmp <- rep(NA, lid)
names(id.tmp) <- nid
## create exact block placeholders
ex.tmp <- rep(NA, lex)
names(ex.tmp) <- nex
qqq <- rbind(x.ex[,1:(ncol(x.ex)-1)], c(id.tmp, ex.tmp, covar.vals))
qqq[nrow(qqq), 1:lid] <- id.vals
if(!is.null(exact.vals)){
qqq[nrow(qqq), (lid+1):(lid+lex)] <- exact.vals
}
orig <- rbind(orig, c(qqq[nrow(qqq), ], Tr=NA))
## When non-exact blocking covariates are used:
## (Default: covar order = order in dataframe)
if(lcv > 0){
if(is.null(ocv)){
ocv <- ncv
}
## Select covars: if prev = 1, pick 1.
if(prev == 1){
cov.tmp <- ocv[1]
}
## if prev>=2, pick prev-1 covars.
## to define MD, need num covars, lcv >= prev.
## for distinct MD, need num covars, lcv >= prev+2
if(prev >= 2){
if(prev-2 >= lcv){
cov.tmp <- ocv ##changed to ocv from covars.ord
}else{
cov.tmp <- ocv[1:(prev-1)]
}
}
qqq.c <- data.frame(qqq[ , cov.tmp])
names(qqq.c) <- cov.tmp
## If there is a variable with no variation, remove it from the vcov matrix for this assignment:
wh.cut <- NULL
for(jj in 1:(ncol(qqq.c))){
if(isTRUE(var(qqq.c[, jj]) == 0)){
wh.cut <- append(wh.cut, jj)
}
}
if(length(wh.cut) > 0){
n.qqq.c <- names(qqq.c)[-wh.cut]
qqq.c <- data.frame(qqq.c[, -wh.cut])
names(qqq.c) <- n.qqq.c
}
## If qqq.c has zero columns, then, go recreate qqq.c.
## (This occurs if second unit is identical to first on first nonexact blocking variable, e.g.)
if(ncol(qqq.c) == 0){
qqq.c <- data.frame(qqq[ , cov.tmp])
names(qqq.c) <- cov.tmp
}
## If new unit identical to old unit, jitter new unit's covar vals
for(jj in 1:(nrow(qqq.c)-1)){
if(isTRUE(sum((qqq.c[jj,] - qqq.c[nrow(qqq.c),])==0) == length(qqq.c[1,]))){
jit.val <- array()
for(kk in 1:ncol(qqq.c)){
jit.val[kk] <- jitter(qqq.c[nrow(qqq.c), kk])
}
qqq[nrow(qqq), cov.tmp] <- qqq.c[nrow(qqq.c), ] <- jit.val
}
} ## To implement: jitter for identical rows w/ identical NAs.
## replace NA's with variable means
for(jj in 1:ncol(qqq.c)){
if(sum(is.na(qqq.c[, jj])) > 0 ){
qqq.c[which(is.na(qqq.c[, jj])), jj] <- mean(qqq.c[, jj], na.rm = TRUE)
}
}
## Different distance measures
if(is.character(dist)){
## Since mve and mcd require conditions on IQR, check, warn, use nonresistant:
if(dist %in% c("mve", "mcd")){
quan <- floor((nrow(qqq.c)+ncol(qqq.c)+1)/2)
if(quan < (ncol(qqq.c)+1)){
warning(paste("'Quantile' must be at least ", ncol(qqq.c)+1, " when using option '", dist, "'. Blocking will proceed using nonresistant Mahalanobis distance scaling matrix.", sep = ""))
dist <- "mahalanobis"
} else if(quan > (nrow(qqq.c)-1)){
warning(paste("'Quantile' must be at most ", nrow(qqq.c)-1, " when using option '", dist, "'. Blocking will proceed using nonresistant Mahalanobis distance scaling matrix.", sep = ""))
dist <- "mahalanobis"
}
}
if(dist %in% c("mve", "mcd")){
iqr.idx <- 0
while(iqr.idx < ncol(qqq.c)){
iqr.idx <- iqr.idx + 1
iqr.tmp <- unname(quantile(qqq.c[, iqr.idx], c(.25, .75)))
if(isTRUE(all.equal(iqr.tmp[1], iqr.tmp[2]))){
warning(paste("Variable ", colnames(qqq.c)[iqr.idx], " has IQR 0; blocking will proceed using nonresistant Mahalanobis distance scaling matrix.", sep = ""))
dist <- "mahalanobis"
iqr.idx <- ncol(qqq.c)
}
}
}
if(dist == "mahalanobis"){
vc.all <- var(qqq.c)
} else if(dist == "mcd"){
vc.all <- cov.rob(qqq.c, method="mcd", seed = seed.dist, ...)$cov
} else if(dist == "mve"){
vc.all <- cov.rob(qqq.c, method="mve", seed = seed.dist, ...)$cov
} else if(dist == "euclidean"){
vc.all <- diag(ncol(qqq.c))
} else{
dist <- "mahalanobis"
vc.all <- var(qqq.c)
}
} ## End 'if(is.character(dist))'
mahmat <- mahal(qqq.c, vc.all)
tr.dist <- list()
lrow <- mahmat[nrow(mahmat), 1:(ncol(mahmat)-1)]
## for(ii in levels(data1[, ncol(data1)])){} ## good for factor TR
for(ii in unique(x.ex[, ncol(x.ex)])){ ## good for character TR
tr.dist[[ii]] <- lrow[which(x.ex[, ncol(x.ex)]==ii)]
}
## Calculate summary for comparison
if(apstat == "mean"){
ms <- lapply(tr.dist, mean) ## only prev assg'd tr's
}
if(apstat == "median"){
ms <- lapply(tr.dist, median)
}
if(apstat == "trimmean"){
ms <- lapply(tr.dist, mean, trim = mtrim)
}
## Sort previously assigned treatments, largest first
tr.sort <- names(sort(unlist(ms), decreasing = TRUE))
## Add treatments not yet assigned, in random order
if(length(tr.sort) != n.tr){
trn.unassg <- sample(trn[!(trn %in% names(ms))])
nunassg <- length(trn.unassg)
tr.sort <- c(trn.unassg, tr.sort)
}
if((length(ms) != length(tr.sort)) && (apmeth %in% c("prop", "prop2", "wprop"))){
ms <- append(ms, rep(NA, nunassg), after=0)
names(ms)[1:nunassg] <- trn.unassg
## Give unassigned tr's double distance of max condition dist
## (Quasi-minimization)
ms[1:nunassg] <- 2*max(unlist(ms), na.rm=T)
}
if((length(tr.dist) != length(tr.sort)) && (apmeth == "wprop")){
tr.dist <- append(tr.dist, rep(NA, nunassg), after=0)
names(tr.dist)[1:nunassg] <- trn.unassg
## Give unassigned tr's double distance of max condition dist,
## with weight of one unit.
## (Quasi-minimization)
tr.dist[1:nunassg] <- max(unlist(ms), na.rm=T)
}
} ## End 'if(lcv > 0)'
if(lcv == 0){ # if there are ONLY exact covariates used:
tr.dist <- list()
tr.sort <- character()
tr.counts <- table(x.ex$Tr)
trn.unassg <- trn[!(trn %in% c(x.ex$Tr))]
if(length(trn.unassg) > 0){
tr.counts <- c(tr.counts, rep(0, length(trn.unassg)))
names(tr.counts)[(n.tr-length(trn.unassg)+1):n.tr] <- trn.unassg
}
p.lcv0 <- (1-tr.counts/sum(tr.counts))/sum(1-tr.counts/sum(tr.counts))
} # End 'if(lcv == 0)'
# Set assignment probability method
if(apmeth == "ktimes"){
p <- c(kfac/(kfac + n.tr - 1),
rep(1 / (kfac + n.tr - 1), n.tr - 1))
}
if(apmeth=="fixed"){
if(length(assgpr) != n.tr){
stop("assg.prob.method is `fixed', but assg.prob (vector of assignment probabilities) is length ", length(assgpr), ", while n.tr (number of treatment conditions) is ", n.tr, ". Respecify and try again.")
}
if(sum(assgpr) != 1){
stop("assg.prob.method is `fixed', but assg.prob (vector of assignment probabilities) sums to ", sum(assgpr), " instead of 1. Respecify and try again.")
}
p <- assgpr
}
if(apmeth == "prop"){
mssum <- sapply(ms, sum)
p <- c(unname((mssum/sum(mssum))[tr.sort]))
}
if(apmeth == "prop2"){
sumsq <- sapply(ms, sum)^2
p <- c(unname((sumsq/sum(sumsq))[tr.sort]))
}
# Biases toward tr's w/fewer units
if(apmeth == "wprop"){
nnnj <- lapply(tr.dist, length)
nnn <- sum(unlist(nnnj))
div.func <- function(x){return(nnn/x)}
nnnw <- lapply(nnnj, div.func)
wdist <- unlist(ms)*unlist(nnnw)
p <- c(unname((wdist/sum(wdist))[tr.sort]))
}
# set user-assigned seed
if(!is.null(seed)){
set.seed(seed)
}
# Assign new TR, given at least 1 "treated as continuous" covars
if(lcv > 0){
tr.new <- sample(tr.sort, 1, prob=p)
}
# Assign new TR using correct order for probs and Tr condition names, EXACT covars only
if(lcv == 0){
tr.new <- sample(names(tr.counts), 1, prob=p.lcv0)
}
# Assign new TR using same order for fixed probs and Tr condition names
if(apmeth == "fixed"){
tr.new <- sample(trn, 1, prob=p)
}
# Append to x
x <- rbind(x, c(qqq[nrow(qqq),], Tr=tr.new))
orig[nrow(orig),"Tr"] <- tr.new
tr.counts <- table(x$Tr)
if(is.null(file.name)){
file.name <- "sbout2k.RData"
}
if(substr(file.name, start=nchar(file.name)-5, stop=nchar(file.name)) !=".RData"){
warning("Destination file name does not end in .RData.")
}
prev.dates <- bdata$datetime
# Create storage list
bdata <- list()
bdata$x <- x
bdata$nid <- nid
bdata$nex <- nex
bdata$ncv <- ncv
bdata$rex <- rex # EXACT restricted values list
if(!is.null(bdata$rex)){
names(bdata$rex) <- nex
}
bdata$rcv <- rcv # BLOCK restricted values list
if(!is.null(bdata$rcv)){
names(bdata$rcv) <- covar.vars
}
bdata$ocv <- ocv # block covariates order
bdata$trn <- trn
bdata$apstat <- apstat # assignment prob statistic
bdata$mtrim <- mtrim # trim fraction for apmeth=trimmean
bdata$apmeth <- apmeth # assignment prob method
bdata$kfac <- kfac # assg probs for method=ktimes
bdata$assgpr <- assgpr # assg probs for fixed prob
bdata$distance <- dist ## distance calculation method
bdata$trd <- tr.dist
bdata$tr.sort <- tr.sort
if (lcv > 0){
bdata$p <- p
}else{
bdata$p <- p.lcv0
}
bdata$trcount <- tr.counts
bdata$datetime <- append(prev.dates, date())
bdata$orig <- orig
save(bdata, file = file.name)
if(verbose == TRUE){
cat("Unit 1:", nrow(x), " data stored as file ", file.name, ".\nThe current working directory is ", getwd(), "\n", sep="")
cat("Unit ", nrow(x), " assigned to ", tr.new, ".\n", sep="")
cat("The new data as entered:\n")
print(orig[nrow(orig), ])
}
} # End 'if(!is.null(object))'
return(bdata)
}
Try the blockTools package in your browser
Any scripts or data that you put into this service are public.
blockTools documentation built on April 4, 2025, 3:42 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 seqblock
Documented in seqblock
#' @title Sequential assignment of unit(s) into experimental conditions using covariates
#'
#' @description
#' Sequentially assign units into experimental conditions. Blocking begins by creating a measure of multivariate distance between a \emph{current} unit and one or multiple \emph{prior}, already-assigned unit(s). Then, average distance between current unit and each treatment condition is calculated, and random assignment is biased toward conditions more dissimilar to current unit. Argument values can be specified either as argument to the function, or via a query. The query directly asks the user to identify the blocking variables and to input, one-by-one, each unit's variable values.
#'
#' @details
#' The \code{seqblock} function's code is primarily divided into two parts: the first half deals with instances, in which the unit being assigned is the first unit in a given study to receive an assignment; the second half addresses subsequent units that are assigned after at least one first assignment has already been made. If the \code{object} argument is left as \code{NULL}, the function will run the first half; if the \code{object} argument is specified, the second part will be executed. When \code{object = NULL}, the researcher has no past file from which to pull variable names and past data; this corresponds to the case when the unit being assigned is the first one. If the researcher does specify \code{object}, it implies the user is drawing data from a past file, which means this is not the first unit in the study to be assigned to a treatment.
#'
#' However, the function can be called for subsequent units even when \code{object} is not specified. By setting \code{query = TRUE}, the console will ask the researcher whether this is the first unit to be assigned in the study. Based on the researcher's response, it will decide which part of the code to run.
#'
#' If the \code{object} and \code{file.name} arguments are set to the same value, then \code{seqblock} overwrites the specified file with a new file, which now contains both the previously-assigned units and the newly-assigned unit. To create a new file when a new unit is assigned, use a new \code{file.name}.
#'
#' The \code{single} algorithm (see \code{exact.alg} in the Arguments section above) creates a variable that has a unique level for every possible combination of the exact variables. As an example, say that there were 3 exact blocking variables: \emph{party} (Democrat, Republican); \emph{region} (North, South); and \emph{education} (HS, NHS). The \code{single} algorithm creates one level for units with the following values: Democrat-North-HS. It would create another level for Democrat-North-NHS; a third level for Republican-North-HS; and so forth, until every possible combination of these 3 variables has its own level. Thus if there are \eqn{k} exact blocking variables and each exact blocking variable has \eqn{q_{i}} values it can take on, then there are a total of \eqn{\prod_{1}^{k} q_{i}} levels created.
#'
#' The \code{distance = "mcd"} and \code{distance = "mve"} options call \code{cov.rob} to calculate measures of multivariate spread robust to outliers. The \code{distance = "mcd"} option calculates the Minimum Covariance Determinant estimate (Rousseeuw 1985); the \code{distance = "mve"} option calculates the Minimum Volume Ellipsoid estimate (Rousseeuw and van Zomeren 1990). When \code{distance = "mcd"}, the interquartile range on blocking variables should not be zero. The \code{distance = "euclidean"} option calculates the Euclidean distance between the new unit and the previously-assigned units. The default \code{distance = "mahalanobis"} option calculates the Mahalanobis distance.
#'
#' @param object a character string giving the file name of a \code{.RData} file containing a list output from the \code{seqblock} function which contains at least one previously assigned unit.
#' @param id.vars a string or vector of strings specifying the name of the identifying variable(s); if \code{query = FALSE} and the object argument is not given, then the \code{id.vars} argument is required.
#' @param id.vals a vector of ID values for every unit being assigned to a treatment group; those are corresponding to the \code{id.vars}.
#' @param exact.vars a string or vector of strings containing the names of each of the exact blocking variables.
#' @param exact.vals a vector containing the unit's values on each of the exact blocking variables.
#' @param exact.restr a list object containing the restricted values that the exact blocking variables can take on. Thus the first element of \code{exact.restr} is a vector containing all of the possible values that the first exact blocking variable (see \code{exact.vars} above) can take on; the second element is a vector containing all of the possible values for the second exact blocking variable; and so on.
#' @param exact.alg a string specifying the blocking algorithm. Currently the only acceptable value is \code{"single"}. This algorithm creates a variable with a unique level for every possible combination of the values in all of the exact variables. See Details section below.
#' @param covar.vars a string or vector of strings containing the names of each of the non-exact blocking variables.
#' @param covar.vals a vector containing the unit's values on each of the non-exact blocking variables.
#' @param covar.restr a list object containing the restricted values that the non-exact blocking variables can take on. Thus the first element of \code{covar.restr} is a vector containing all of the possible values that the first non-exact blocking variable (see \code{covar.vars} above) can take on; the second element is a vector containing all of the possible values for the second non-exact blocking variable; and so on.
#' @param covars.ord a string or vector of strings containing the name of the non-exact blocking variables ordered so that the highest priority covariate comes first, followed by the second highest priority covariate, then the third, etc.
#' @param n.tr the number of treatment groups. If not specified, this defaults to \code{n.tr = 2}.
#' @param tr.names a string or vector of strings containing the names of the different treatment groups.
#' @param assg.prob a numeric vector containing the probabilities that a unit will be assigned to the treatment groups; this vector should sum to 1.
#' @param seed an optional integer value for the random seed, which is used when assigning units to treatment groups.
#' @param seed.dist an optional integer value for the random seed set in \code{cov.rob}, used to calculate measures of the variance-covariance matrix robust to outliers.
#' @param assg.prob.stat a string specifying which assignment probability summary statistic to use; valid values are \code{mean}, \code{median}, and \code{trimmean}. If not specified, this defaults to \code{assg.prob.stat = "mean"}.
#' @param trim a numeric value specifying what proportion of the observations are to be dropped from each tail when the assignment probability summary statistic (\code{assg.prob.stat}) is set equal to \code{trimmean}. Blocks on each tail of the distribution are dropped before the mean is calculated. If not specified, this defaults to \code{trim = 0.1}.
#' @param assg.prob.method a string specifying which algorithm should be used when assigning treatment probabilities. Acceptable values are \code{ktimes}, \code{fixed}, \code{prop}, \code{prop2}, and \code{wprop}. If not specified, this defaults to \code{assg.prob.method = "ktimes"}.
#' @param assg.prob.kfac a numeric value for \code{k}, the factor by which the most likely experimental condition will be multiplied relative to the other conditions. If not specified, this defaults to \code{assg.prob.kfac = 2}.
#' @param distance a string specifying how the multivariate distance used for blocking covariates are calculated. If not specified, this defaults to \code{distance = "mahalanobis"}.
#' @param file.name a string containing the name of the file that one would like the output to be written to. Ideally this file name should have the extension .RData.
#' @param query a logical stating whether the console should ask the user questions to input the data and assign a treatment condition. If not specified, this defaults to \code{query = FALSE}.
#' @param verbose a logical stating whether the function should print the name of the output file, the current working directory, the treatment group that the most recent unit was assigned to, and the dataframe \code{x} returned by the function as part of the \code{bdata} list. If not specified, this defaults to \code{verbose = TRUE}.
#' @param \dots additional arguments.
#'
#' @return
#' A list (called \code{bdata}) with elements
#' \itemize{
#' \item \strong{x}: a dataframe containing the names and values for the different ID and blocking variables, as well as each unit's initial treatment assignment.
#' \item \strong{nid}: a string or vector of strings containing the name(s) of the ID variable(s).
#' \item \strong{nex}: a string or vector of strings containing the name(s) of the exact blocking variable(s).
#' \item \strong{ncv}: a string or vector of strings containing the name(s) of the non-exact blocking variable(s).
#' \item \strong{rex}: a list of the restricted values of the exact blocking variables.
#' \item \strong{rcv}: a list of the restricted values of the non-exact blocking variables.
#' \item \strong{ocv}: a vector of the order of the non-exact blocking variables.
#' \item \strong{trn}: a string or vector of strings containing the name(s) of the different treatment groups.
#' \item \strong{apstat}: a string specifying the assignment probability summary statistic that was used.
#' \item \strong{mtrim}: a numeric value specifying the proportion of observations to be dropped when the assignment probability statistic takes on the value \code{"trimmean"}.
#' \item \strong{apmeth}: a string specifying the assignment probability algorithm that was used.
#' \item \strong{kfac}: the assignment probability \emph{kfactor}; see \emph{assg.prob.kfac} in the Arguments section above.
#' \item \strong{assgpr}: a vector of assignment probabilities to each treatment group.
#' \item \strong{distance}: a string specifying how the multivariate distance used for blocking was calculated.
#' \item \strong{trd}: a list with the length equal to the number of previously assigned treatment conditions; each object in the list contains a vector of the distance between each unit in one treatment group and the new unit. This will be \code{NULL} when there are no non-exact blocking variables.
#' \item \strong{tr.sort}: a string vector of treatment conditions, sorted from the largest to the smallest. Set to \code{NULL} when there are no non-exact blocking variables.
#' \item \strong{p}: a vector of assignment probabilities to each treatment group used in assigning a treatment condition to the new unit.
#' \item \strong{distance}: a string specifying how the multivariate distance used for blocking is calculated
#' \item \strong{trcount}: a table containing the counts for each experimental/treatment conditions.
#' \item \strong{datetime}: the date and time at which each unit was assigned their treatment group.
#' \item \strong{orig}: a dataframe containing the names and values for the different id and blocking variables, as well as each unit's treatment assignment.
#' }
#'
#' @examples
#' # Assign first unit (assume a 25 year old member of the Republican Party)
#' # to a treatment group. Save the results in file "sdata.RData":
#' # seqblock(query = FALSE, id.vars = "ID", id.vals = 001, exact.vars = "party",
#' # exact.vals = "Republican", covar.vars = "age", covar.vals = 25, file.name = "sdata.RData")
#'
#' # Assign next unit (age 30, Democratic Party):
#' # seqblock(query = FALSE, object = "sdata.RData", id.vals = 002, exact.vals = "Democrat",
#' # covar.vars = "age", covar.vals = 30, file.name = "sdata.RData")
#'
#' @references
#' Moore, Ryan T. and Sally A. Moore. 2013. "Blocking for Sequential Political Experiments."
#' \emph{Political Analysis} 21(4):507-523.
#'
#' Moore, Ryan T. 2012. "Multivariate Continuous Blocking to Improve Political Science Experiments." \emph{Political Analysis} 20(4):460-479.
#'
#' Rousseeuw, Peter J. 1985. "Multivariate Estimation with High Breakdown Point". \emph{Mathematical Statistics and Applications} 8:283-297.
#'
#' Rousseeuw, Peter J. and Bert C. van Zomeren. 1990. "Unmasking Multivariate Outliers and Leverage Points". \emph{Journal of the American Statistical Association} 85(411):633-639.
#'
#' @seealso \code{\link{assignment}}, \code{\link{block}}
#'
#' @author Ryan T. Moore \email{rtm@american.edu}, Tommy Carroll \email{tcarroll22@wustl.edu}, Jonathan Homola \email{homola@wustl.edu} and Jeong Hyun Kim \email{jeonghyun.kim@wustl.edu}
#'
#' @keywords design multivariate
#'
#' @export
seqblock <- function(object = NULL, id.vars, id.vals, exact.vars = NULL, exact.vals = NULL, exact.restr = NULL, exact.alg = "single", covar.vars = NULL, covar.vals = NULL, covar.restr = NULL, covars.ord = NULL, n.tr = 2, tr.names = NULL, assg.prob = NULL, seed = NULL, seed.dist, assg.prob.stat = NULL, trim = NULL, assg.prob.method = NULL, assg.prob.kfac = NULL, distance = NULL, file.name = NULL, query = FALSE, verbose = TRUE, ...){
if(is.null(object)){
if(query == TRUE){
init <- readline("Is this the first unit that is being assigned in this experiment? [y/n] ")
if(!(substr(init, 1, 1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
init <- "y"
}else(stop("The seqblock function requires the ''object'' argument to be set to a valid file name for subsequent unit assignment."))
}
if(substr(init, 1, 1) == "n"){
object <- readline("Enter the name of the input file without quotation marks. [E.g., sbout1.RData] ")
}
}
}
if(is.null(object)){ ## If assigning first unit:
if(query == TRUE){
## ID
lid <- readline("How many identification variables are there? ")
id.vars <- id.vals <- NULL
for(ii in 1:lid){ # For each identification variables,
## tmp stores the name of ID variable:
tmp <- readline(cat("Enter the name of ID variable ", ii, " without quotation marks.", sep=""))
## Add the ii-th name of ID variable to a vector id.vars:
id.vars <- append(id.vars, tmp)
## tmp2 stores the value of each ID variable:
tmp2 <- readline(cat("Enter the value of '", tmp, "'. ", sep=""))
## Add the ii-th value of ID variable to a vector id.vals
id.vals <- append(id.vals, tmp2)
} ## End loop over ID variables
nnn <- readline("Should the ID values be numeric? [y/n] ")
if(!(substr(nnn, 1, 1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline("Should the ID values be numeric? [y/n] ")
}
}
if(substr(nnn,1,1) != "n"){ # If the respondent answers that the ID values should be numeric,
id.vals <- as.numeric(id.vals) # change id.vals to numerics.
}
if(sum(is.na(id.vals)) > 0){ # If there is at least one missing value in id.vals,
tmp <- readline(cat("Warning: ID value(s)", which(is.na(id.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
if(!(substr(tmp, 1, 1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp <- readline(cat("Warning: ID value(s)", which(is.na(id.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
}
}
if(substr(tmp, 1, 1) == "n"){
stop()
}
}
## Specifying EXACT blocking variables:
lex <- readline("`Exact blocking variables' are those on which you require covariate values to be identical. How many exact blocking variables are there? ")
exact.vars <- exact.vals <- exact.restr <- NULL
if(lex > 0){ ## If there are exact blocking variables
for(ii in 1:lex){ ## Then, for each exact blocking variable
## tmp stores the name of exact blocking variable:
tmp <- readline(cat("Enter the name of exact block variable ", ii, " without quotation marks.", sep=""))
exact.vars <- append(exact.vars, tmp) ## then add to the vector exact.vars
tmp2 <- readline(cat("Enter the value of '", tmp, "'. ", sep=""))
tmp3 <- readline(cat("Should '", tmp, "' be restricted to certain values? [n/y] ", sep=""))
if(!(substr(tmp3,1,1) %in% c("y", "n"))){ ## if the answer is neither yes nor no,
ddd <- readline("The default is 'no'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp3 <- readline(cat("Should '", tmp, "' be restricted to certain values? [n/y] "))
}
} ## Close 'if'
if(substr(tmp3, 1, 1) == "y"){ ## if the ii-th exact variable should be restricted to certain values, tmp4 stores the values that variable can take that the respondent inputs.
tmp4 <- readline(cat("How many values can '", tmp, "' take? ", sep=""))
## if the value is less than one, execution is halted.
if(tmp4 < 1){
stop("Restricted variables must have at least one valid value. ")
}
exact.restr <- as.list(exact.restr)
exact.restr[[tmp]] <- NA
for(rrr in 1:tmp4){ ## for each value that the variable tmp can take,
tmp5 <- readline(cat("Enter possible value for '", tmp, "' number ", rrr, sep=""))
exact.restr[[tmp]] <- append(exact.restr[[tmp]], tmp5) ## the object named tmp in the list exact.restr will be a vector of the values that the variable can take.
} ## End loop over values 'tmp' can take
} ## Close 'if'
## if the ii-th exact variable should NOT be restricted to certain values, exact.vals simply contains the values of exact blocking variables.
exact.vals <- append(exact.vals, tmp2)
} ## Close loop over exact blocking variables 'ii in 1:lex'
nnn <- readline("Should the exact block values be numeric? [y/n] ")
if(!(substr(nnn, 1, 1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline("Should the exact block values be numeric? [y/n] ")
}
}
## If the answer is yes (the exact block values should be numeric),
if(substr(nnn,1,1) != "n"){
## change values in exact.vals into numerics:
exact.vals <- as.numeric(exact.vals)
## If there is a restriction on the values of exact variable,
if(length(exact.restr) > 0){
for(rrr in 1:length(exact.restr)){
## Each element of those values are changed to numerics.
exact.restr[[rrr]] <- as.numeric(exact.restr[[rrr]])
}
}
}
if(sum(is.na(exact.vals)) > 0){
tmp <- readline(cat("Warning: Exact blocking value(s)", which(is.na(exact.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
if(!(substr(tmp, 1, 1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp <- readline(cat("Warning: Exact blocking value(s)", which(is.na(exact.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
}
}
if(substr(tmp,1,1) == "n"){
stop()
}
} ## End confirm NA values for exact blocking vars
## CHECK that each exact.val is in exact.restr here.
} ## End 'if(lex > 0)'
## COVARIATES
lcv <- readline("How many other blocking variables are there? ")
covar.vars <- covar.vals <- covar.restr <- NULL
if(lcv > 0){ ## If there is at least one non-exact blocking variable,
for(ii in 1:lcv){ ## For each of those variables,
tmp <- readline(cat("Enter the name of blocking variable ", ii, " without quotation marks.", sep=""))
## Save its name to a vector 'covar.vars'
covar.vars <- append(covar.vars, tmp)
tmp2 <- readline(cat("Enter the value of '", tmp, "'. ", sep=""))
tmp3 <- readline(cat("Should '", tmp, "' be restricted to certain values? [n/y] ", sep=""))
if(!(substr(tmp3, 1, 1) %in% c("y", "n"))){
ddd <- readline("The default is 'no'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp3 <- readline(cat("Should '", tmp, "' be restricted to certain values? [n/y] ", sep=""))
}
} ## close if
## If the value of the blocking variable should be restricted to certain values,
if(substr(tmp3,1,1) == "y"){
tmp4 <- readline(cat("How many values can '", tmp, "' take? ", sep=""))
if(tmp4 < 1){
stop("Restricted variables must have at least one valid value. ")
}
covar.restr <- as.list(covar.restr)
covar.restr[[tmp]] <- NA
for(rrr in 1:tmp4){
tmp5 <- readline(cat("Enter possible value for '", tmp, "' number ", rrr, sep=""))
# Make a list covar.restr in which each object represents restricted values of blocking variable that should be restricted to certain values.
covar.restr[[tmp]] <- append(covar.restr[[tmp]], tmp5)
}
}
## If the values should NOT be restricted to certain values, covar.vals contains the values of the blocking variable:
covar.vals <- append(covar.vals, tmp2)
} ## End loop over blocking vars 'ii in 1:lcv'
nnn <- readline("Should the blocking values be numeric? [y/n] ")
if(!(substr(nnn,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline("Should the blocking values be numeric? [y/n] ")
}
}
if(substr(nnn,1,1) != "n"){
covar.vals <- as.numeric(covar.vals)
}else{
stop("The blocking variables are stored as character strings.\nCurrent implementation does not support distance calculations for character strings.")
} ## (Change to warning() when categorical distances incorporated)
if(sum(is.na(covar.vals)) > 0){
tmp <- readline(cat("Warning: Blocking value(s)", which(is.na(covar.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
if(!(substr(tmp,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp <- readline(cat("Warning: Blocking value(s)", which(is.na(covar.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
}
}
if(substr(tmp,1,1) == "n"){
stop()
}
}
## Covariates' order
nnn <- readline("Would you like to specify the order in which blocking variables will be used while (number units considered) < (number covariates + 2) ? (If not, covariates will be added as possible in the order in which they were entered.) [n/y] ")
if(!(substr(nnn,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `no'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline(cat("Would you like to specify the order in which blocking variables will be used while (number units considered) < (number covariates +2) ? [n/y] "))
}
}
if(substr(nnn,1,1) == "y"){
for(ii in 1:lcv){
tmp <- readline(cat("Enter the name of blocking variable prioritized number ", ii, ".", sep=""))
# if(!(tmp %in% names(x)[(lid+lex+1):(lid+lex+lcv)])){
if(!(tmp %in% covar.vars)){
tmp <- readline(cat("`", tmp, "' is not a blocking variable. Please try again. Enter the name of blocking variable prioritized number ", ii, ".", sep=""))
}
if(tmp %in% covars.ord){
tmp <- readline(cat("`", tmp, "' has already been prioritized. Please try again. Enter the name of blocking variable prioritized number ", ii, ".", sep=""))
}
covars.ord <- append(covars.ord, tmp)
}
}
} ## End 'if(lcv > 0)'
## Number of treatment conditions:
n.tr <- readline("How many experimental/treatment conditions are there? ")
n.tr <- as.numeric(n.tr)
## Names of treatment conditions:
nnn <- readline(cat("Would you like to specify the names of the experimental conditions? [y/n] [If not, `Treatment 1', ..., `Treatment ", n.tr, "', will be used] ", sep=""))
if(substr(nnn,1,1) != "n"){
tr.names <- NULL
for(ii in 1:n.tr){
tmp <- readline(cat("Enter the name of experimental condition ", ii, " without quotation marks.", sep=""))
tr.names <- append(tr.names, tmp)
}
}
## Initial condition assignment probabilities:
nnn <- readline("Would you like to specify the initial assignment probabilities? [y/n] [If not, all treatment assignment probabilities for the first unit will be equal.] ")
if(!(substr(nnn,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline("Would you like to specify the initial assignment probabilities? [y/n] [If not, all treatment assignment probabilities for the first unit will be equal.] ")
}
}
if(substr(nnn,1,1) != "n"){
assg.prob <- NULL
for(ii in 1:n.tr){
tmp <- readline(cat("Enter the assignment probability of experimental condition ", ii, ".", sep=""))
assg.prob <- append(assg.prob, tmp)
}
assg.prob <- as.numeric(assg.prob)
}
## Set random seed:
nnn <- readline("Would you like to specify a random seed for the initial experimental condition assignment? [y/n] ")
if(!(substr(nnn,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline("Would you like to specify a random seed for the initial experimental condition assignment? [y/n] ")
}
}
if(substr(nnn,1,1) != "n"){
seed <- as.numeric(readline("Enter the random seed. "))
}
## ASSIGNMENT PROBABILITY STATISTIC:
nnn <- readline("Would you like to specify the assignment probability summary statistic? [n/y] [If not, `mean' will be used.] ")
if(substr(nnn,1,1) == "y"){
assg.prob.stat <- readline("Enter the summary statistic name. [mean, median, trimmean] ")
if(assg.prob.stat == "trimmean"){
ooo <- readline("Specify the proportion of blocks to be dropped from each end of the distribution before the mean is calculated. [Defaults to 0.1] ")
if(substr(ooo,1,1) == ""){
trim <- 0.1
}else{
trim <- ooo
}
}
}else{
assg.prob.stat <- "mean"
}
## ASSIGNMENT PROBABILITY METHOD:
nnn <- readline("Would you like to specify the assignment probability algorithm? [n/y] [If not, `ktimes' with k=2 will be used.] ")
if(substr(nnn,1,1) == "y"){
assg.prob.method <- readline("Enter the algorithm name. [ktimes, fixed, prop, prop2, wprop] ")
if(assg.prob.method == "ktimes"){
assg.prob.kfac <- as.numeric(readline("Enter the value of k, the factor by which the most likely experimental condition will be multiplied, relative to the other conditions. "))
}
}else{
assg.prob.method <- "ktimes"
}
if(assg.prob.method == "fixed"){
assgpr <- readline("Enter 1st probability. ")
for(pc in 2:n.tr){
assg.prob <- append(assg.prob, readline(cat("Enter probablity number ", pc, ". ", sep="")))
}
assg.prob <- as.numeric(assg.prob)
}
## MULTIVARIATE DISTANCE CALCULATION:
nnn <- readline("Would you like to specify how the multivariate distance used for blocking is calculated? [n/y] [If not, `mahalanobis' will be used.] ")
if(substr(nnn,1,1) == "y"){
distance <- readline("Enter the calculation type. [mahalanobis, mcd, mve, euclidean] ")
}else{
distance <- "mahalanobis"
}
## OUTPUT FILE NAME:
nnn <- readline("Would you like to specify the output data file name? [y/n] [If not, `sbout1.RData' will be used.] ")
if(!(substr(nnn,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline("Would you like to specify the output data file name? [y/n] [If not, `sbout1.RData' will be used.] ")
}
}
if(substr(nnn,1,1) != "n"){
file.name <- readline("Enter the file name. Note: seqblock2k() expects a file of type `.RData'. ")
}
} ## END (query == T)
## create variables
lid <- length(id.vars)
lex <- length(exact.vars)
lcv <- length(covar.vars)
## create vector of treatment names
if(is.null(tr.names)){
tr.names <- paste("Treatment", 1:n.tr)
}
## create kfac, if null
if(!(is.null(assg.prob.kfac))){
kfac <- assg.prob.kfac
} else {kfac <- assg.prob.kfac <- 2}
## create assg.prob.stat, if null
if(!(is.null(assg.prob.stat))){
assg.prob.stat <- assg.prob.stat
} else {assg.prob.stat <- "mean"}
## create trim, if null
if(!(is.null(trim))){
trim <- trim
} else {trim <- 0.1}
## create assg.prob.method, if null
if(!(is.null(assg.prob.method))){
assg.prob.method <- assg.prob.method
} else {assg.prob.method <- "ktimes"}
## create distance, if null
if(is.null(distance)){
distance <- "mahalanobis"
}
## create vector of (equal) initial treatment probabilities
if(is.null(assg.prob)){
assg.prob <- rep(1/n.tr, n.tr)
}
## check that assg.prob sum to 1
if(sum(assg.prob)!=1){
stop("Initial assignment probabilities do not sum to 1.\nRespecify and try again.")
}
## check that if exact named, then supplied
if((!is.null(exact.vars)) && is.null(exact.vals)){
stop("Exact (grouped) blocking requested, but no values given.\nSpecify exact block values and try again.")
}
## check that exact vars and vals are same length
if(lex != length(exact.vals)){
stop(cat("Number of exact blocking variables (", lex, ") does not equal number of values (", length(exact.vals), ").\nSpecify exact block values and try again.\n", sep=""))
}
## check that if covars named, then supplied
if((!is.null(covar.vars)) && is.null(covar.vals)){
stop("Covariate blocking requested, but no values given.\nSpecify covariate values and try again.")
}
## check that covar vars and vals are same length
if(length(covar.vars) != length(covar.vals)){
stop("Number of covariate blocking variables does not equal number of values.\nSpecify covariate block values and try again.")
}
## Add NA to each existing restriction set,
## then check if any exact vals violate exact.restr
if((lex > 0) && !(is.null(exact.restr))){
for(i in names(exact.restr)){
if(!(NA %in% exact.restr[[i]])){
exact.restr[[i]] <- append(exact.restr[[i]], NA)
}
wexrstr <- which(exact.vars == i)
if(!(exact.vals[wexrstr] %in% exact.restr[[exact.vars[wexrstr]]])){
stop("Exact blocking value number ", wexrstr, " (", exact.vals[wexrstr], ") is not in the set of values to which exact blocking variable `", i, "' is restricted. Respecify either the value, the restriction set, or both as needed." )
}
}
}
## check if any covar vals violate covar.restr
if((lcv > 0) && (!is.null(covar.restr))){
for(i in names(covar.restr)){
if(!(NA %in% covar.restr[[i]])){
covar.restr[[i]] <- append(covar.restr[[i]], NA)
}
wcvrstr <- which(covar.vars == i)
if(!(covar.vals[wcvrstr] %in% covar.restr[[covar.vars[wcvrstr]]])){
stop("Covariate blocking value number ", wcvrstr, " (", covar.vals[wcvrstr], ") is not in the set of values to which covariate blocking variable `", i, "' is restricted. Respecify either the value, the restriction set, or both as needed." )
}
}
}
## set user-assigned seed
if(!is.null(seed)){
set.seed(seed)
}
## Draw initial treatment assignment
init.tr <- sample(tr.names, 1, replace = FALSE, prob = assg.prob)
## Create data frame of unit 1 data
## create id placeholders
id.tmp <- rep(NA, lid)
names(id.tmp) <- id.vars
## create group placeholders
ex.tmp <- rep(NA, lex)
names(ex.tmp) <- exact.vars
## create data frame
x <- unlist(c(id.tmp, ex.tmp, unname(covar.vals), NA))
x <- data.frame(t(x))
colnames(x) <- c(id.vars, exact.vars, covar.vars, "Tr")
## add id vals. NOTE: will be all char or all num
x[1, 1:lid] <- id.vals
## add group vals. NOTE: will be all char or all num
if(!is.null(exact.vars)){
x[1, (lid+1):(lid+lex)] <- exact.vals
}
## add initial treatment assignment
x[1,ncol(x)] <- init.tr
if(is.null(file.name)){
file.name <- "sbout1.RData"
}
if(substr(file.name, start=nchar(file.name)-5, stop=nchar(file.name)) !=".RData"){
warning("Destination file name does not end in .RData.")
}
## Create storage list 'bdata'
bdata <- list()
bdata$x <- x
bdata$nid <- id.vars
bdata$nex <- exact.vars
bdata$ncv <- covar.vars
bdata$rex <- exact.restr ## EXACT restricted values list
if(!is.null(bdata$rex)){
names(bdata$rex) <- exact.vars
}
bdata$rcv <- covar.restr ## BLOCK restricted values list
if(!is.null(bdata$rcv)){
names(bdata$rcv) <- covar.vars
}
bdata$ocv <- covars.ord ## block covariates order
bdata$trn <- tr.names
bdata$apstat <- assg.prob.stat
bdata$mtrim <- trim
bdata$apmeth <- assg.prob.method
bdata$kfac <- kfac ## multiple for method 'ktimes'
bdata$assgpr <- assg.prob
bdata$distance <- distance
bdata$datetime <- date()
bdata$orig <- x
## Write bdata to file
save(bdata, file = file.name)
if(verbose == TRUE){
cat("Unit 1 data stored as file ", file.name, ".\nThe current working directory is ", getwd(), "\n", sep="")
cat("Unit ", nrow(x), " assigned to ", init.tr, ".\n", sep="")
cat("The new data as entered:\n")
print(x)
}
} ## End 'if(is.null(object))'. I.e., end assigning first unit.
## Next case: if object is not NULL (I.e., begin assignment for subsequent unit after the first one)
if(!is.null(object)){
load(object) ## loads, doesn't export to wkspace
## Rename object elements
x <- bdata$x
nid <- bdata$nid ## names of id.vars
nex <- bdata$nex ## names of exact.vars
ncv <- bdata$ncv ## names of covar.vars
rex <- bdata$rex ## EXACT restricted values list
rcv <- bdata$rcv ## BLOCK restricted values list
ocv <- bdata$ocv ## block covariates order
trn <- bdata$trn ## treatment condition names
apstat <- bdata$apstat ## assignment prob statistic
mtrim <- bdata$mtrim ## trim fraction for apmeth=trimmean
apmeth <- bdata$apmeth ## assignment prob method
kfac <- bdata$kfac ## assg multiple for method 'ktimes'
assgpr <- bdata$assgpr ## assg probs for fixed prob
dist <- bdata$distance ## distance calculation method
orig <- bdata$orig ## original unjittered data
n.tr <- length(trn) ## number of treatment conditions
lid <- length(nid) #number of id variables
lex <- length(nex) #number of exact variables
lcv <- length(ncv) #number of covariates
if(!(is.null(exact.restr))){
rex <- exact.restr
}
if(!(is.null(covar.restr))){
rcv <- covar.restr
}
if(!(is.null(covars.ord))){
ocv <- covars.ord
}
if(!(is.null(assg.prob.stat))){
if(all.equal(assg.prob.stat,apstat) != TRUE){
warning(paste("'assg.prob.stat' of previous assignments was ", apstat, ". 'assg.prob.stat' for this assignment is ", assg.prob.stat, ".", sep = ""))
}
apstat <- assg.prob.stat
}
if(!(is.null(trim))){
if(all.equal(trim,mtrim) != TRUE){
warning(paste("'trim' of previous assignments was ", mtrim, ". 'trim' for this assignment is ", trim, ".", sep = ""))
}
mtrim <- trim
}
if(!(is.null(assg.prob.method))){
if(all.equal(assg.prob.method, apmeth) != TRUE){
warning(paste("'assg.prob.method' of previous assignments was ", apmeth, ". 'assg.prob.method' for this assignment is ", assg.prob.method, ".", sep = ""))
}
apmeth <- assg.prob.method
}
if(!(is.null((assg.prob.kfac)))){
if(all.equal(assg.prob.kfac,kfac) != TRUE){
warning(paste("'kfac' of previous assignments was ", kfac, ". 'kfac' for this assignment is ", assg.prob.kfac, ".", sep = ""))
}
kfac <- assg.prob.kfac
}
if(!(is.null(assg.prob))){
assgpr <- assg.prob
}
if(!(is.null((distance)))){
if(all.equal(distance,dist) != TRUE){
warning(paste("Distance option of previous assignments was '", dist, "'. Distance option for this assignment is '", distance, "'.", sep = ""))
}
dist <- distance
}
if(query == TRUE){
## ID
id.vals <- NULL
for(ii in 1:lid){
tmp <- names(x)[ii]
tmp2 <- readline(cat("Enter the value of '", tmp, "'. ", sep=""))
id.vals <- append(id.vals, tmp2)
}
class(id.vals) <- class(x[,lid]) ## lid could be any of 1:lid
if(sum(is.na(id.vals)) > 0){ # If there is at least one missing value in id.vals,
tmp <- readline(cat("Warning: ID value(s)", which(is.na(id.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
if(!(substr(tmp, 1, 1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp <- readline(cat("Warning: ID value(s)", which(is.na(id.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
}
}
if(substr(tmp, 1, 1) == "n"){
stop()
}
}
## Check for duplicated IDs if single id.var
if((lid == 1) && (id.vals %in% x[, nid])){
tmp <- readline(cat("Warning: ID value ", id.vals, " already used in the data. Allow duplication and proceed? [y/n]", sep = ""))
if(substr(tmp, 1, 1) == "n"){
stop()
}
}
## EXACT
exact.vals <- NULL
if(lex > 0){
for(ii in 1:lex){
tmp <- nex[ii]
tmp2 <- readline(cat("Enter the value of '", tmp, "'. ", sep=""))
exact.vals <- append(exact.vals, tmp2)
}
class(exact.vals) <- class(x[,(lid+lex)]) ## could be any of (lid+1):(lid+lex) If the exact blocking variable was set to be numeric in seqblock1(), new exact values should be also numeric. Otherwise, the values of the exact blocking variable will be coerced to NAs.
if(sum(is.na(exact.vals)) > 0){
tmp <- readline(cat("Warning: Exact blocking value(s)", which(is.na(exact.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
if(!(substr(tmp,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp <- readline(cat("Warning: Exact blocking value(s)", which(is.na(exact.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
}
}
if(substr(tmp, 1, 1) == "n"){
stop()
}
}
}
## EXACT ALGORITHM
## "single": if lex > 1, creates unique categories (4 for 2x2, e.g.)
exact.alg <- "single"
## COVARIATES
covar.vals <- NULL
if(lcv > 0){
for(ii in 1:lcv){
tmp <- ncv[ii]
tmp2 <- readline(cat("Enter the value of '", tmp, "'. ", sep=""))
covar.vals <- append(covar.vals, tmp2)
}
class(covar.vals) <- class(x[,lid+lex+lcv]) ## could be any of (lid+lex+1):(lid+lex+lcv)
if(sum(is.na(covar.vals)) > 0){
tmp <- readline(cat("Warning: Blocking value(s)", which(is.na(covar.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
if(!(substr(tmp,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
tmp <- readline(cat("Warning: Blocking value(s)", which(is.na(covar.vals)), "was/were coerced to and stored as `NA'. Proceed? [y/n] ", sep=" "))
}
}
if(substr(tmp,1,1) == "n"){
stop()
}
}
}
## SET RANDOM SEED:
nnn <- readline("Would you like to specify a random seed for this experimental condition assignment? [y/n] ")
if(!(substr(nnn,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline(cat("Would you like to specify a random seed for this experimental condition assignment? [y/n] "))
}
}
if(substr(nnn,1,1) != "n"){
seed <- as.numeric(readline("Enter the random seed. "))
}
## OUTPUT FILE NAME:
nnn <- readline("Would you like to specify the output data file name? [y/n] [If not, `sbout2k.RData' will be used.] ")
if(!(substr(nnn,1,1) %in% c("y", "n"))){
ddd <- readline("The default is `yes'. Continue? [y/n] ")
if(substr(ddd, 1, 1) != "n"){
}else{
nnn <- readline(cat("Would you like to specify the output data file name? [y/n] [If not, `sbout2k.RData' will be used.] "))
}
}
if(substr(nnn,1,1) != "n"){
file.name <- readline("Enter the file name. Note: seqblock2k() expects a file of type `.RData'. ")
}else{
file.name <- NULL
}
} ## Close query==T
## Warn if duplicated ID values if single ID var
if((lid == 1) && (id.vals %in% x[, nid])){
warning("Warning: ID value ", id.vals, " already used in the data.")
}
## check # exact vars vs. # exact vals
if(length(nex) != length(exact.vals)){
stop(cat("Number of exact blocking variables (", lex, ") does not equal number of values (", length(exact.vals), ").\nSpecify exact block values and try again.\n", sep=""))
}
## check if any exact val violates exact.restr
if((lex > 0) && !(is.null(rex))){
for(i in names(rex)){
if(!(NA %in% rex[[i]])){
rex[[i]] <- append(rex[[i]], NA)
}
wexrstr <- which(nex == i)
if(!(exact.vals[wexrstr] %in% rex[[nex[wexrstr]]])){
tmp <- readline(cat("Exact blocking value number ", wexrstr, " (", exact.vals[wexrstr], ") is not in the set of values to which exact blocking variable `", i, "' is restricted. Would you like to add ", exact.vals[wexrstr], " to the set? [y/n]", sep=""))
if(substr(tmp,1,1) != "n"){
rex[[i]] <- append(rex[[i]], exact.vals[wexrstr])
}else{
stop("Exact blocking value number ", wexrstr, " (", exact.vals[wexrstr], ") is not in the set of values to which exact blocking variable `", nex[wexrstr], "' is restricted. Respecify either the value, the restriction set, or both as needed." )
}
}
}
}
## check if any block val violates covar.restr
if((lcv > 0) && (!is.null(rcv))){ for(i in names(rcv)){
if(!(NA %in% rcv[[i]])){
rcv[[i]] <- append(rcv[[i]], NA)
}
wcvrstr <- which(ncv == i)
if(!(covar.vals[wcvrstr] %in% rcv[[ncv[wcvrstr]]])){
tmp <- readline(cat("Covariate blocking value number ", wcvrstr, " (", covar.vals[wcvrstr], ") is not in the set of values to which covariate blocking variable `", i, "' is restricted. Would you like to add ", covar.vals[wcvrstr], " to the set? [y/n]", sep=""))
if(substr(tmp,1,1) != "n"){
rcv[[i]] <- append(rcv[[i]], covar.vals[wcvrstr])
}else{
stop("Covariate blocking value number ", wcvrstr, " (", covar.vals[wcvrstr], ") is not in the set of values to which covariate blocking variable `", ncv[wcvrstr], "' is restricted. Respecify either the value, the restriction set, or both as needed.")
}
}
}
}
if(!(length(ocv) == length(unique(ocv)))){
stop("Prioritization of covariates for (number units considered) < (number covariates + 2) includes duplicate variable names. Respecify 'covars.ord'.")
}
if(lex > 0){
if(exact.alg == "single"){
x.ex <- x[apply(t(t(x[,(lid+1):(lid+lex)])==exact.vals), 1, sum) == lex,]
}
}
## When no exact matches, use all previous data:
if(is.null(exact.vals) || (nrow(x.ex)==0)){
x.ex <- x
}
prev <- nrow(x.ex)
## Create data frame of unit 1 data
## create id placeholders
id.tmp <- rep(NA, lid)
names(id.tmp) <- nid
## create exact block placeholders
ex.tmp <- rep(NA, lex)
names(ex.tmp) <- nex
qqq <- rbind(x.ex[,1:(ncol(x.ex)-1)], c(id.tmp, ex.tmp, covar.vals))
qqq[nrow(qqq), 1:lid] <- id.vals
if(!is.null(exact.vals)){
qqq[nrow(qqq), (lid+1):(lid+lex)] <- exact.vals
}
orig <- rbind(orig, c(qqq[nrow(qqq), ], Tr=NA))
## When non-exact blocking covariates are used:
## (Default: covar order = order in dataframe)
if(lcv > 0){
if(is.null(ocv)){
ocv <- ncv
}
## Select covars: if prev = 1, pick 1.
if(prev == 1){
cov.tmp <- ocv[1]
}
## if prev>=2, pick prev-1 covars.
## to define MD, need num covars, lcv >= prev.
## for distinct MD, need num covars, lcv >= prev+2
if(prev >= 2){
if(prev-2 >= lcv){
cov.tmp <- ocv ##changed to ocv from covars.ord
}else{
cov.tmp <- ocv[1:(prev-1)]
}
}
qqq.c <- data.frame(qqq[ , cov.tmp])
names(qqq.c) <- cov.tmp
## If there is a variable with no variation, remove it from the vcov matrix for this assignment:
wh.cut <- NULL
for(jj in 1:(ncol(qqq.c))){
if(isTRUE(var(qqq.c[, jj]) == 0)){
wh.cut <- append(wh.cut, jj)
}
}
if(length(wh.cut) > 0){
n.qqq.c <- names(qqq.c)[-wh.cut]
qqq.c <- data.frame(qqq.c[, -wh.cut])
names(qqq.c) <- n.qqq.c
}
## If qqq.c has zero columns, then, go recreate qqq.c.
## (This occurs if second unit is identical to first on first nonexact blocking variable, e.g.)
if(ncol(qqq.c) == 0){
qqq.c <- data.frame(qqq[ , cov.tmp])
names(qqq.c) <- cov.tmp
}
## If new unit identical to old unit, jitter new unit's covar vals
for(jj in 1:(nrow(qqq.c)-1)){
if(isTRUE(sum((qqq.c[jj,] - qqq.c[nrow(qqq.c),])==0) == length(qqq.c[1,]))){
jit.val <- array()
for(kk in 1:ncol(qqq.c)){
jit.val[kk] <- jitter(qqq.c[nrow(qqq.c), kk])
}
qqq[nrow(qqq), cov.tmp] <- qqq.c[nrow(qqq.c), ] <- jit.val
}
} ## To implement: jitter for identical rows w/ identical NAs.
## replace NA's with variable means
for(jj in 1:ncol(qqq.c)){
if(sum(is.na(qqq.c[, jj])) > 0 ){
qqq.c[which(is.na(qqq.c[, jj])), jj] <- mean(qqq.c[, jj], na.rm = TRUE)
}
}
## Different distance measures
if(is.character(dist)){
## Since mve and mcd require conditions on IQR, check, warn, use nonresistant:
if(dist %in% c("mve", "mcd")){
quan <- floor((nrow(qqq.c)+ncol(qqq.c)+1)/2)
if(quan < (ncol(qqq.c)+1)){
warning(paste("'Quantile' must be at least ", ncol(qqq.c)+1, " when using option '", dist, "'. Blocking will proceed using nonresistant Mahalanobis distance scaling matrix.", sep = ""))
dist <- "mahalanobis"
} else if(quan > (nrow(qqq.c)-1)){
warning(paste("'Quantile' must be at most ", nrow(qqq.c)-1, " when using option '", dist, "'. Blocking will proceed using nonresistant Mahalanobis distance scaling matrix.", sep = ""))
dist <- "mahalanobis"
}
}
if(dist %in% c("mve", "mcd")){
iqr.idx <- 0
while(iqr.idx < ncol(qqq.c)){
iqr.idx <- iqr.idx + 1
iqr.tmp <- unname(quantile(qqq.c[, iqr.idx], c(.25, .75)))
if(isTRUE(all.equal(iqr.tmp[1], iqr.tmp[2]))){
warning(paste("Variable ", colnames(qqq.c)[iqr.idx], " has IQR 0; blocking will proceed using nonresistant Mahalanobis distance scaling matrix.", sep = ""))
dist <- "mahalanobis"
iqr.idx <- ncol(qqq.c)
}
}
}
if(dist == "mahalanobis"){
vc.all <- var(qqq.c)
} else if(dist == "mcd"){
vc.all <- cov.rob(qqq.c, method="mcd", seed = seed.dist, ...)$cov
} else if(dist == "mve"){
vc.all <- cov.rob(qqq.c, method="mve", seed = seed.dist, ...)$cov
} else if(dist == "euclidean"){
vc.all <- diag(ncol(qqq.c))
} else{
dist <- "mahalanobis"
vc.all <- var(qqq.c)
}
} ## End 'if(is.character(dist))'
mahmat <- mahal(qqq.c, vc.all)
tr.dist <- list()
lrow <- mahmat[nrow(mahmat), 1:(ncol(mahmat)-1)]
## for(ii in levels(data1[, ncol(data1)])){} ## good for factor TR
for(ii in unique(x.ex[, ncol(x.ex)])){ ## good for character TR
tr.dist[[ii]] <- lrow[which(x.ex[, ncol(x.ex)]==ii)]
}
## Calculate summary for comparison
if(apstat == "mean"){
ms <- lapply(tr.dist, mean) ## only prev assg'd tr's
}
if(apstat == "median"){
ms <- lapply(tr.dist, median)
}
if(apstat == "trimmean"){
ms <- lapply(tr.dist, mean, trim = mtrim)
}
## Sort previously assigned treatments, largest first
tr.sort <- names(sort(unlist(ms), decreasing = TRUE))
## Add treatments not yet assigned, in random order
if(length(tr.sort) != n.tr){
trn.unassg <- sample(trn[!(trn %in% names(ms))])
nunassg <- length(trn.unassg)
tr.sort <- c(trn.unassg, tr.sort)
}
if((length(ms) != length(tr.sort)) && (apmeth %in% c("prop", "prop2", "wprop"))){
ms <- append(ms, rep(NA, nunassg), after=0)
names(ms)[1:nunassg] <- trn.unassg
## Give unassigned tr's double distance of max condition dist
## (Quasi-minimization)
ms[1:nunassg] <- 2*max(unlist(ms), na.rm=T)
}
if((length(tr.dist) != length(tr.sort)) && (apmeth == "wprop")){
tr.dist <- append(tr.dist, rep(NA, nunassg), after=0)
names(tr.dist)[1:nunassg] <- trn.unassg
## Give unassigned tr's double distance of max condition dist,
## with weight of one unit.
## (Quasi-minimization)
tr.dist[1:nunassg] <- max(unlist(ms), na.rm=T)
}
} ## End 'if(lcv > 0)'
if(lcv == 0){ # if there are ONLY exact covariates used:
tr.dist <- list()
tr.sort <- character()
tr.counts <- table(x.ex$Tr)
trn.unassg <- trn[!(trn %in% c(x.ex$Tr))]
if(length(trn.unassg) > 0){
tr.counts <- c(tr.counts, rep(0, length(trn.unassg)))
names(tr.counts)[(n.tr-length(trn.unassg)+1):n.tr] <- trn.unassg
}
p.lcv0 <- (1-tr.counts/sum(tr.counts))/sum(1-tr.counts/sum(tr.counts))
} # End 'if(lcv == 0)'
# Set assignment probability method
if(apmeth == "ktimes"){
p <- c(kfac/(kfac + n.tr - 1),
rep(1 / (kfac + n.tr - 1), n.tr - 1))
}
if(apmeth=="fixed"){
if(length(assgpr) != n.tr){
stop("assg.prob.method is `fixed', but assg.prob (vector of assignment probabilities) is length ", length(assgpr), ", while n.tr (number of treatment conditions) is ", n.tr, ". Respecify and try again.")
}
if(sum(assgpr) != 1){
stop("assg.prob.method is `fixed', but assg.prob (vector of assignment probabilities) sums to ", sum(assgpr), " instead of 1. Respecify and try again.")
}
p <- assgpr
}
if(apmeth == "prop"){
mssum <- sapply(ms, sum)
p <- c(unname((mssum/sum(mssum))[tr.sort]))
}
if(apmeth == "prop2"){
sumsq <- sapply(ms, sum)^2
p <- c(unname((sumsq/sum(sumsq))[tr.sort]))
}
# Biases toward tr's w/fewer units
if(apmeth == "wprop"){
nnnj <- lapply(tr.dist, length)
nnn <- sum(unlist(nnnj))
div.func <- function(x){return(nnn/x)}
nnnw <- lapply(nnnj, div.func)
wdist <- unlist(ms)*unlist(nnnw)
p <- c(unname((wdist/sum(wdist))[tr.sort]))
}
# set user-assigned seed
if(!is.null(seed)){
set.seed(seed)
}
# Assign new TR, given at least 1 "treated as continuous" covars
if(lcv > 0){
tr.new <- sample(tr.sort, 1, prob=p)
}
# Assign new TR using correct order for probs and Tr condition names, EXACT covars only
if(lcv == 0){
tr.new <- sample(names(tr.counts), 1, prob=p.lcv0)
}
# Assign new TR using same order for fixed probs and Tr condition names
if(apmeth == "fixed"){
tr.new <- sample(trn, 1, prob=p)
}
# Append to x
x <- rbind(x, c(qqq[nrow(qqq),], Tr=tr.new))
orig[nrow(orig),"Tr"] <- tr.new
tr.counts <- table(x$Tr)
if(is.null(file.name)){
file.name <- "sbout2k.RData"
}
if(substr(file.name, start=nchar(file.name)-5, stop=nchar(file.name)) !=".RData"){
warning("Destination file name does not end in .RData.")
}
prev.dates <- bdata$datetime
# Create storage list
bdata <- list()
bdata$x <- x
bdata$nid <- nid
bdata$nex <- nex
bdata$ncv <- ncv
bdata$rex <- rex # EXACT restricted values list
if(!is.null(bdata$rex)){
names(bdata$rex) <- nex
}
bdata$rcv <- rcv # BLOCK restricted values list
if(!is.null(bdata$rcv)){
names(bdata$rcv) <- covar.vars
}
bdata$ocv <- ocv # block covariates order
bdata$trn <- trn
bdata$apstat <- apstat # assignment prob statistic
bdata$mtrim <- mtrim # trim fraction for apmeth=trimmean
bdata$apmeth <- apmeth # assignment prob method
bdata$kfac <- kfac # assg probs for method=ktimes
bdata$assgpr <- assgpr # assg probs for fixed prob
bdata$distance <- dist ## distance calculation method
bdata$trd <- tr.dist
bdata$tr.sort <- tr.sort
if (lcv > 0){
bdata$p <- p
}else{
bdata$p <- p.lcv0
}
bdata$trcount <- tr.counts
bdata$datetime <- append(prev.dates, date())
bdata$orig <- orig
save(bdata, file = file.name)
if(verbose == TRUE){
cat("Unit 1:", nrow(x), " data stored as file ", file.name, ".\nThe current working directory is ", getwd(), "\n", sep="")
cat("Unit ", nrow(x), " assigned to ", tr.new, ".\n", sep="")
cat("The new data as entered:\n")
print(orig[nrow(orig), ])
}
} # End 'if(!is.null(object))'
return(bdata)
}
Try the blockTools package in your browser
Any scripts or data that you put into this service are public.
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.