Nothing
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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.