Nothing
##############################################################################
# Copyright (c) 2008-2010, 2012-2016 Russell V. Lenth #
# #
# This file is part of the rsm package for R (*rsm*) #
# #
# *rsm* is free software: you can redistribute it and/or modify #
# it under the terms of the GNU General Public License as published by #
# the Free Software Foundation, either version 2 of the License, or #
# (at your option) any later version. #
# #
# *rsm* is distributed in the hope that it will be useful, #
# but WITHOUT ANY WARRANTY; without even the implied warranty of #
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
# GNU General Public License for more details. #
# #
# A copy of the GNU General Public License is available at #
# <https://www.r-project.org/Licenses/> #
##############################################################################
### Code for CCDs and related designs
# Modified by RVL starting Feb 10, 2011
### Find list of row indexes for each block in a design
.block.indices = function(design) {
rsd = attr(design, "rsdes")
if (!is.null(blknm <- rsd$block)) {
blk = design[[blknm[1]]]
if (length(blknm) > 1) for (nm in blknm[-1])
blk = paste(blk, design[[nm]], sep=".")
# Even now, it's possible to have a null blk, if that var not in data yet
if (!is.null(blk))
i.init = split(1:nrow(design), blk)
else
i.init = list(1:nrow(design))
}
else
i.init = list(1:nrow(design))
i.init
}
### Randomizes; adds columns "std.order" and "run.order" if it's a not-previously-randomized design
.randomize = function(design, randomize=TRUE) {
OneToN = 1:nrow(design)
if (is.na(match("std.order", names(design)))) {
design$run.order = design$std.order = OneToN
extcols = match(c("run.order","std.order"), names(design))
oldcols = (1:ncol(design))[-extcols]
design = `[.coded.data`(design, , c(extcols, oldcols))
}
if (randomize) {
i.init = .block.indices(design)
for (idx in i.init) {
design[sample(idx), ] = design[idx, ]
design$run.order[idx] = 1:length(idx)
}
row.names(design) = OneToN
}
design
}
### Calculate design moments needed for determining orthogonal blocks
### Computes (N, [i], [ii], [ij]) for each block (really just the [i1])
### Returns the matrix of [ii] values (col names are "N.obs" and pri vars, rows are blocks)
### For any block where any [i] != 0 or [ij] != 0, we return NA
### For orth blocks, need [ii]/N.obs constant on each block for each i
### (but not necessary for [ii] == [jj])
.orth.moments = function(design) {
blkidx = .block.indices(design)
if (is.null(rsd <- attr(design, "rsdes")))
pricols = sapply(design, is.numeric)
else
pricols = match(rsd$primary, names(design))
t(sapply(blkidx, function(idx) {
N.obs = length(idx)
X = as.matrix(design[idx, pricols])
x1 = X[,1]
moms = apply(X, 2, function(x)
c(i = sum(x), ii = sum(x^2), i1 = sum(x * x1)))
if (any(abs(moms[1, ]) > 0.01) || (any(abs(moms[3,-1]) > .01)))
moms[2, ] = NA
c(N.obs=N.obs, moms[2, ])
}))
}
### Calculate design moments needed for determining rotatability
### Computes design moments of order 4 or less and checks whether rotatability
### is possible. We need all except [iiii] and [iijj] to be zero,
### and the latter to be the same for all i <> j
### Returns iiii and iijj
### NOTE: For the design to be rotatable, we need [iiii]/[iijj] = 3
### For mixed moments, we actually check only [1j], [12k], etc.
.rot.moments = function(design) {
if (is.null(rsd <- attr(design, "rsdes")))
pricols = sapply(design, is.numeric)
else
pricols = match(rsd$primary, names(design))
k = length(pricols)
if (k < 2)
stop("Need at least 2 primary variables for a rotatable design")
X = as.matrix(design[, pricols])
x1 = X[,1];
x2 = X[,2];
x3 = if (k==2) 0 else X[,3]
moms = apply(X, 2, function(x)
c(i = sum(x), ij = sum(x*x1),
iii = sum(x^3), iij = sum(x^2*x1), ijk = sum(x*x1*x2),
iiij = sum(x^3*x1), iijk = sum(x^2*x1*x2), ijkl = sum(x*x1*x2*x3),
iiii = sum(x^4), iijj = sum(x^2 * x1^2)))
# We need to check abs values of 1st 8 rows of moms
moms[1:8,] = abs(moms[1:8,])
# Zero out mixed parts that really aren't mixed
moms["ij",1] = moms["iij",1] = moms["ijk",1] = moms["ijk",2] =
moms["iiij",1] = moms["iijk",1] = moms["iijk",2] = 0
if (k > 2) moms["ijkl",1:3] = 0
# for checking purposes...
moms["iijj", 1] = moms["iijj", 2]
nzflag = any(moms[1:8,] > 0.01)
if (nzflag || (sd(moms["iiii",]) > .1) || (sd(moms["iijj",]) > .1))
moms["iiii", ] = moms["iijj", ] = NA
c(iiii=moms["iiii",1], iijj=moms["iijj",2])
}
### Cube (1st-order) design -- this is a 2^k or 2^{k-p}, possibly replicated, plus center points
# inscribed may be boolean or anticipated value of alpha. If TRUE, we use sqrt(k) (spherical)
cube = function(basis, generators, n0=4, reps=1, coding, randomize=TRUE,
blockgen, bid=1, inscribed=FALSE)
{
if (inherits(basis, "formula"))
xvars = all.vars(basis[[length(basis)]])
else if (is.numeric(basis))
xvars = paste("x", 1:basis, sep="")
else
stop("'basis' must be an integer or a formula")
CALL = match.call()
args = lapply(xvars, function(nm) c(-1,1))
names(args) = xvars
des = do.call(expand.grid, args)
if (!missing(generators)) {
if (!is.list(generators)) generators = list(generators)
for (gen in generators) {
gen = as.character(gen)
des[[gen[[2]]]] = with(des, eval (parse (text = as.character(gen[[3]]))))
}
}
k = ncol(des)
xvars = names (des)
# within-block reps
if (reps > 1) {
des = des[rep(1:nrow(des), each=reps), ]
}
rsd = list(primary=xvars)
if (!missing(blockgen)) {
blk = NULL
if (!is.vector(blockgen)) {
if (inherits(blockgen, "formula")) {
bg = with(des, eval(blockgen[[length(blockgen)]])) ###blockgen = c(blockgen)
if (is.numeric(bg)) blk = matrix(bg, nrow(des))
}
else {
bg = with(des, eval(blockgen))
# should either give a list of formulas, or results
if (is.numeric(bg))
blk = matrix(bg, nrow(des))
else
blockgen = bg
}
}
if (is.vector(blockgen) || is.null(blk))
blk = with(des, sapply(blockgen, function(f) {
if (is.character(f))
f = as.formula(paste("~", f))
if (!inherits(f, "formula"))
stop("block generators must be formulas or strings")
if (length(f)==3) f = f[-2]
eval(f[[2]])
}))
if (is.matrix(blk)) {
fac = 2^(1:ncol(blk)) / 2 # = 1, 2, 4, ...
blk = (blk %*% fac + 2^(ncol(blk)) + 1) / 2
}
else
blk = (blk + 3) / 2
rep.rows = sapply(1:max(blk), function(b) which(blk==b)[1])
rsd$blk.signs = des[rep.rows, ]
rsd$bid = bid
des = des[blk == bid, ]
}
# center points
if (n0 > 0) {
zero = as.data.frame(matrix(0, nrow=n0, ncol=k))
names(zero) = names(des)
des = rbind(des, zero)
}
if (inscribed) {
if (is.logical(inscribed)) inscribed = sqrt(k)
rsd$divisor = inscribed
# Note: divisor is used only by star()
# Some trickery here. ccd() will give a negative value; we'll wait for
# star() to fill-in a value, and scale the design then
if (inscribed > 0)
des = des / inscribed
}
# Add vars from left-hand side, if any
if (inherits(basis, "formula") & (length(basis) > 2)) {
yvars = all.vars(basis[[2]])
for (v in yvars) des[[v]] = NA
### not used rsd$other=yvars
}
# Figure out sort order
### row.names(des) = 1:nrow(des)
des = .randomize(des, randomize=randomize)
if (missing(coding)) {
coding = sapply(xvars, function(v) as.formula(paste(v,"~",v,".as.is", sep="")))
}
des = as.coded.data (des, formulas=coding)
# rsd$n0 = c(n0,NA)
# rsd$non0 = c(nrow(des)-n0,NA)
# rsd$alpha = NA
rsd$call = CALL
# rsd$blk.info = list(rsd[c("n0","non0","alpha")])
attr(des, "rsdes") = rsd
des
}
### new version of ccd -- constructs it using cube, star, dupe, foldover, and djoin
ccd = function(basis, generators, blocks="Block", n0=4, alpha="orthogonal",
wbreps=1, bbreps=1, randomize=TRUE, inscribed=FALSE, coding,
oneblock = FALSE)
{
n0 = rep(n0,2)[1:2]
wbreps = rep(wbreps,2)[1:2]
bbreps = rep(bbreps,2)[1:2]
# Note: We play tricks and pass a negative value if inscribed is TRUE
if (inscribed && is.logical(inscribed)) inscribed = -999
if (is.character(blocks)) { # need to call cube with missing blockgen
cp = des = cube(basis, generators, n0=n0[1], reps=wbreps[1], coding=coding,
randomize=randomize, inscribed=inscribed)
blkname = blocks
}
else if (inherits(blocks, "formula")) { # we have fractional blocks
blkname = if(length(blocks) == 3) as.character(blocks[[2]]) else "Block"
### get only the rhs to make life easier for cube
blocks = blocks[[length(blocks)]]
cp = des = cube(basis, generators, blockgen=blocks, n0=n0[1], reps=wbreps[1],
coding=coding, randomize=randomize, bid=1, inscribed=inscribed)
# we need to add the other fractional blocks
nfb = nrow(attr(cp, "rsdes")$blk.signs)
if(nfb > 1) for (i in 2:nfb)
des = djoin(des, foldover(cp, bid=i, randomize=randomize), blkname=blkname)
cp = des
}
else
stop("improper 'blocks' argument; must be a string or a formula")
# Take care of replicate blocks
if(bbreps[1] > 1) for (i in 2:bbreps[1]) {
des = djoin(des, dupe(cp, randomize=randomize), blkname=blkname)
}
# Star part
sp = star(des, n0=n0[2], alpha=alpha, reps=wbreps[2], randomize=randomize)
# and add it the right no of times
for (i in 1:bbreps[2])
des = djoin(des, dupe(sp, randomize=randomize), blkname=blkname)
if (inscribed) {
div = attr(des, "rsdes")$divisor = attr(sp, "rsdes")$divisor
pri = attr(des, "rsdes")$primary
des[, pri] = des[, pri] / div
}
if (oneblock) {
des[[blkname]] = NULL
if (randomize)
des = .randomize(des, TRUE)
}
des
}
### Make a copy of a design, re-randomizing as instructed
dupe = function(design, randomize=TRUE, coding) {
des = .randomize(design, randomize=randomize)
if(!missing(coding))
codings(des) <- coding
des
}
### New star function
# slightly older version, commented out - we're replacing with version that ac=tually
# checks design moments
# star = function(basis, n0=4, alpha="orthogonal", reps=1, randomize=TRUE) {
# CALL = match.call()
# if (missing(basis)) return(CALL)
# rsd = NULL # so these are named outside of a program block
# if (inherits(basis, "data.frame")) {
# rsd = attr(basis, "rsdes")
# if (is.null(rsd))
# stop("'basis' design must have 'rsdes' attribute.\nTry updating it via 'as.coded.data'")
# }
# else
# stop("Requires a coded.data object")
# privars = rsd$primary
# k = length(privars)
#
# # internal code to test for whether we have a cube block
# .is.cube = function(bi) {
# if (is.null(bi)) FALSE
# else {
# n0 = bi$n0
# if (is.null(n0)) FALSE
# else !is.na(bi$n0[1])
# }
# }
#
# if (!is.numeric(alpha)) {
# mt = pmatch(alpha, c("orthogonal", "rotatable", "spherical", "faces"))
# if (is.na(mt)) stop("illegal 'alpha' option")
# got.it = FALSE
# F.tot = 0
# # Should be simpler than this, but this works...
# if (is.null(rsd$blk.info)) {
# if (.is.cube(rsd)) {
# f = c(rsd$non0, rsd$n0)
# F.tot = F.tot + rsd$non0[1]
# got.it = TRUE
# }
# }
# else {
# for (ann in rev(rsd$blk.info)) {
# if (.is.cube(ann)) {
# if(!got.it) f = c(ann$non0, ann$n0)
# F.tot = F.tot + rsd$non0[1]
# got.it = TRUE
# }
# }
# }
#
#
# if (!got.it) stop("No cube blocks available to use in determining alpha")
# f = f[c(1,3)]
# alpha = switch(mt,
# sqrt(f[1]*(2*k*reps + n0)/(2*reps*sum(f))),
# sqrt(sqrt(F.tot / reps)),
# sqrt(k),
# 1
# )
# if (!is.null(rsd$divisor)) {
# if (rsd$divisor > 0) alpha = alpha / rsd$divisor
# else rsd$divisor = alpha
# }
#
# }
#
# des = as.data.frame(matrix(0, nrow=2*reps*k + n0, ncol=k))
# names(des) = privars
# for (i in 1:k) {
# des[2*i-1, i] = -alpha
# des[2*i, i] = alpha
# }
# if (reps > 1) {
# rowrng = 1:(2*k)
# for (i in 1:(reps-1))
# des[rowrng + 2*k*i, ] = des[rowrng, ]
# }
#
# if (is.coded.data(basis))
# des = as.coded.data(des, formulas = codings(basis))
#
# des = .randomize(des, randomize=randomize)
#
# rsd$n0 = c(NA,n0)
# rsd$non0 = c(NA,nrow(des)-n0)
# rsd$alpha = alpha
# attr(des, "rsdes") = rsd
#
# des
# }
star = function(basis, n0=4, alpha="orthogonal", reps=1, randomize=TRUE) {
CALL = match.call()
if (missing(basis)) return(CALL)
rsd = NULL # so these are named outside of a program block
if (inherits(basis, "data.frame")) {
rsd = attr(basis, "rsdes")
if (is.null(rsd))
stop("'basis' design must have 'rsdes' attribute.\nTry updating it via 'as.coded.data'")
}
else
stop("Requires a coded.data object")
privars = rsd$primary
k = length(privars)
N = 2*k*reps + n0
if (!is.numeric(alpha)) {
mt = pmatch(alpha, c("orthogonal", "rotatable", "spherical", "faces"))
if (mt == 1) { # orthogonal -- NOTE could have different alpha for each variable
mom = .orth.moments(basis)
if (all(is.na(mom)))
stop("Orthogonal blocks are not achievable")
if (any(is.na(mom))) {
warning("Some blocks are not orthogonal first-order designs: these were ignored")
mom = mom[!is.na(mom[,2]), , drop=FALSE]
}
else if (!is.matrix(mom))
mom = matrix(mom, nrow=1)
ratios = apply(mom[, -1, drop = FALSE], 2, function(x) x / mom[, 1])
alpha = alphas = sqrt(N*ratios/(2*reps))
if (is.matrix(alphas)) {
alpha = alphas[nrow(alphas), ]
tst = apply(alphas[-nrow(alphas), , drop=FALSE], 1, function(a) {
any(abs(a-alpha) > .01)
})
if (any(tst))
warning(paste("New block is orthogonal to block", nrow(alphas),
"but not to block(s)", which(tst)))
}
}
else if (mt == 2) { # rotatable
mom = .rot.moments(basis)
if (any(is.na(mom)))
stop("Rotatable design is not achievable: inconsistent design moments")
else if (is.matrix(mom))
mom = apply(mom, 2, sum) # overall [iiii] and [iijj] values
alpha = (3*mom[2] - mom[1]) / (2 * reps)
if (alpha < .0001) # after 4th root, need alpha > .1 to be valid
stop("Rotatable design is not achievable: requires small or imaginary alpha")
alpha = alpha^.25
}
else # 3 = spherical, 4 = face
alpha = switch(mt - 2, sqrt(k), 1)
}
# possible to have several alphas -- expand so we have at least k of 'em
if(length(alpha) < k) alpha = rep(alpha, k)
# if (is.na(mt)) stop("illegal 'alpha' option")
# alpha = switch(mt,
# sqrt(f[1]*(2*k*reps + n0)/(2*reps*sum(f))),
# sqrt(sqrt(F.tot / reps)),
# sqrt(k),
# 1
# )
if (!is.null(rsd$divisor)) {
if (rsd$divisor > 0) alpha = alpha / rsd$divisor
else rsd$divisor = max(alpha)
}
des = as.data.frame(matrix(0, nrow=2*reps*k + n0, ncol=k))
names(des) = privars
for (i in 1:k) {
des[2*i-1, i] = -alpha[i]
des[2*i, i] = alpha[i]
}
if (reps > 1) {
rowrng = 1:(2*k)
for (i in 1:(reps-1))
des[rowrng + 2*k*i, ] = des[rowrng, ]
}
if (is.coded.data(basis))
des = as.coded.data(des, formulas = codings(basis))
des = .randomize(des, randomize=randomize)
# rsd$n0 = c(NA,n0)
# rsd$non0 = c(NA,nrow(des)-n0)
# rsd$alpha = alpha
rsd$call = CALL
attr(des, "rsdes") = rsd
des
}
# fold over a design. If variables is missing, use all primary ones
foldover = function(basis, variables, bid, randomize=TRUE) {
CALL = match.call()
if (missing(basis)) {
return (CALL)
}
design = basis
if (!inherits(design, "data.frame")) { # Try a switcheroo
CALL$basis = NULL
CALL$variables = basis
return(CALL)
}
rsd = attr(basis, "rsdes")
if (!missing(bid)) {
if(is.null(rsd)) stop("information required for 'bid' is unavailable")
signs = rsd$blk.signs
if (is.null(signs)) stop("'bid' requires that basis be generated as a fractional block")
s = as.numeric(signs[bid, ]) * as.numeric(signs[rsd$bid, ])
variables = names(signs)[s < 0]
rsd$bid = bid
}
else if (missing(variables)) {
if (is.null(rsd))
variables = names(design)[which(sapply(design, is.numeric))]
else
variables = rsd$primary
}
for (nm in variables)
design[[nm]] = -design[[nm]]
design = .randomize(design, randomize=randomize)
attr(design, "rsdes") = rsd
design
}
djoin = function(design1, design2, ..., blkname = "Block", blocklev) {
if (!inherits(design1, "data.frame"))
stop("design1 must inherit from 'data.frame'.")
rsd1 = attr(design1, "rsdes")
if (is.null(rsd1))
stop(paste("Requires an updated 'coded.data' object\n",
"Please update", as.character(substitute(design1)),
"via 'as.coded.data' before calling"))
# if design2 is actually a call, supply design1 as its 'basis' and call it
if (is.call(design2)) {
design2$basis = design1
design2 = eval(design2, parent.frame())
}
# info for design 2 -- before messing with columns
if (!is.coded.data(design2)) {
newcode = NULL
# if primary vars missing, make a coded dataset
mat = match(rsd1$primary, names(design2))
if (any(is.na(mat)))
design2 = coded.data(design2, formulas=codings(design1))
}
if (is.coded.data(design2)) {
newcode = codings(design2)
# are they the same as design1?
oldcode = codings(design1)
comp = sapply(1:length(oldcode), function(i) newcode[[i]]==oldcode[[i]])
if (!all(comp)) {
design1 = recode.data(design1, formulas=newcode)
}
}
rsd2 = attr(design2, "rsdes")
# preserve annotations of separate blocks -- before modifying anything ...
# if(is.null(rsd1$blk.info))
# rsd1$blk.info = list(rsd1[c("n0","non0","alpha")])
# rsd1$blk.info[[1+length(rsd1$blk.info)]] = rsd2[c("n0","non0","alpha")]
# check names2 vs names1 - any NAs here need to be discarded in design2
cols = match(names(design2), names(design1))
design2 = design2[ , !is.na(cols)]
# check names1 vs names2 - need to add NA columns in design2 for these
cols = match(names(design1), names(design2))
na.nm = names(design1)[is.na(cols)]
for (nm in na.nm) design2[[nm]] = NA
# --- at this point, design1 and design2 now have the same column names
# --- (but in possibly different orders)
# update blocks
blknm = rsd1$block[1]
if (is.null(blknm) || is.na(match(blknm, names(design1)))) { # Create a blocking factor
blknm = blkname
design1[[blknm]] = factor(1)
design2[[blknm]] = factor(2)
}
if (is.na(design2[[blknm]][1]))
design2[[blknm]] = factor(2)
blevs = levels(design1[[blknm]])
nlevs = levels(design2[[blknm]])
if (!missing(blocklev)) {
if (length(blocklev) == length(nlevs)) {
if (all(is.na(match(blocklev, blevs))))
nlevs = blocklev
else
warning("blocklev ignored -- duplicates existing block levels")
}
else
warning("blocklev ignored -- wrong number of levels")
}
# redefine design2's block variable if it matches existing levels
if (any(!is.na(match(nlevs, blevs)))) {
# try extending a numeric sequence
nlevs = as.character(seq(1+length(blevs), length=length(nlevs)))
if (any(!is.na(match(nlevs, blevs))))
nlevs = paste("A", blevs, sep="")
}
levels(design2[[blknm]]) = nlevs
### Join the designs
des = rbind(design1, design2)
# # update attributes
# if (!is.null(rsd2)) {
# for (i in 1:2) {
# if (!is.na(rsd2$n0[i])) rsd1$n0[i] = rsd2$n0[i]
# if (!is.na(rsd2$non0[i])) rsd1$non0[i] = rsd2$non0[i]
# }
# if (!is.na(rsd2$alpha)) rsd1$alpha = rsd2$alpha
# }
rsd1$block = blknm
attr(des, "rsdes") = rsd1
# handle dots
for (dzn in list(...))
des = djoin(des, dzn)
des
}
stdorder = function(design) {
if (!is.coded.data(design))
stop("Requires a coded.data object")
if (!("std.order" %in% names(design)))
stop("No \"std.order\" variable found")
if (!is.null(blknm <- attr(design, "rsdes")$block))
key = apply(design[, c(blknm, "std.order")], 1, paste, collapse="")
else
key = design$std.order
design[order(key), ]
}
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.