Nothing
# file Utility.R
# copyright (C) 2002-2004 by Robert E. Wheeler
#
"expand.formula"<-
function (frml, varNames, const = TRUE, numerics = NULL)
{
env <- environment(frml)
noNumerics <- missing(numerics)
nameargs <- function(...) {
dots <- as.list(substitute(list(...)))[-1]
nm <- names(dots)
fixup <- if (is.null(nm))
seq(along = dots)
else nm == ""
dep <- sapply(dots[fixup], function(x) deparse(x, width.cutoff = 500)[1])
if (is.null(nm))
nm <- dep
else {
nm[fixup] <- dep
}
if (!noNumerics) {
if ((nm[1] == "." && !all(numerics)) || !all(numerics[is.element(varNames,
nm)]))
stop("All arguments to special functions such as quad() must be numeric.")
}
nm
}
find.dot<-function (strng) {
# Returns location of dot in strng or -1 if not found
loc<-regexpr("([^0-9a-zA-Z_]+ *\\.)|(^ *\\.)",strng) #not variable followed by one or more spaces
locp<-loc[1]
if (locp!=-1) { # dot found
leng<-attr(loc,"match.length")
subs<-substring(strng,locp,locp+leng-1) # substring containing dot
locsub<-regexpr("\\.",subs) # loc of dot in substring
return(locp+locsub[1]-1)
}
locp
}
quad <- function(...) {
nms<-nameargs(...)
if (-1!=find.dot(nms)) {
nms <- varNames
nVars <- length(nms)
}
else nVars <- nargs()
strg <- paste(paste("(", paste("X", 1:nVars, sep = "",
collapse = "+"), ")^2", sep = ""), "+", paste("I(X",
1:nVars, "^2)", sep = "", collapse = "+"))
for (i in 1:nVars) {
ag <- paste("X", i, sep = "")
strg <- gsub(ag, nms[i], strg)
}
strg
}
cubic <- function(...) {
nms<-nameargs(...)
if (-1!=find.dot(nms)) {
nms <- varNames
nVars <- length(nms)
}
else nVars <- nargs()
strg <- paste(paste("(", paste("X", 1:nVars, sep = "",
collapse = "+"), ")^3", sep = ""), "+", paste("I(X",
1:nVars, "^2)", sep = "", collapse = "+"), "+", paste("I(X",
1:nVars, "^3)", sep = "", collapse = "+"))
for (i in 1:nVars) {
ag <- paste("X", i, sep = "")
strg <- gsub(ag, nms[i], strg)
}
strg
}
cubicS <- function(...) {
nms<-nameargs(...)
if (-1!=find.dot(nms)) {
nms <- varNames
nVars <- length(nms)
}
else nVars <- nargs()
strg <- paste("(", paste("X", 1:nVars, sep = "", collapse = "+"),
")^3", sep = "")
for (i in 1:(nVars - 1)) {
var <- paste("X", i, sep = "")
strg <- paste(strg, "+", paste(paste("I(", var, paste("*X",
(i + 1):nVars, sep = ""), sep = ""), paste("*(",
var, paste("-X", (i + 1):nVars, "))", sep = "")),
collapse = "+"), sep = "", collapse = "+")
}
for (i in 1:nVars) {
ag <- paste("X", i, sep = "")
strg <- gsub(ag, nms[i], strg)
}
strg
}
findFunction <- function(name, string) {
if (-1 == (strt <- regexpr(name, string)))
return(c(0, 0))
head <- substr(string, 1, strt - 1)
tail <- substr(string, strt, nchar(string))
if (-1 == (fin <- regexpr(")", tail)))
return(c(0, 0))
c(strt, strt + fin - 1)
}
frml <- deparse(frml, width.cutoff = 500)
while ((0 != (pos <- findFunction("quad", frml))[1]) || (0 !=
(pos <- findFunction("cubicS", frml))[1]) || (0 != (pos <- findFunction("cubic",
frml))[1])) {
prog <- substr(frml, pos[1], pos[2])
strHead <- substr(frml, 1, pos[1] - 1)
strTail <- substr(frml, pos[2] + 1, nchar(frml))
prog <- eval(parse(text = prog))
frml <- paste(strHead, prog, strTail, sep = "")
}
if (!const)
frml <- paste(frml, "+0", sep = "")
frml <- as.formula(frml)
environment(frml) <- env
frml
}
"actuals" <-
function(frmals){
# returns the actual argument list of a function
nargs<-length(frmals)
fnames<-names(frmals)
for (i in 1:nargs) {
if (i==1)
actuls=list(eval(parse(text=fnames[i]),sys.frame(-1)))
else
actuls<-c(actuls,list(eval(parse(text=fnames[i]),sys.frame(-1))))
}
names(actuls)<-fnames
actuls
}
"efficient.rounding" <-
function(proportions,n,random=TRUE){
# Efficient rounding, following Pulkesheim and Rieder, (1992). Efficient rounding of
# approximate designs. Biometrika. 79, 763-770
# proportions: a vector of proportions
# n: the integer sum of the rounded proportions
# random: if TRUE will randomly select the weight to change, otherwise the first available will be changed.
l<-length(proportions)
n<-as.integer(n)
if (n<l) # in case input is too small
n<-l
if (any(proportions<0))
stop("Negative weights are not allowed")
if (abs(sum(proportions)-1)>1e-12)
stop("proportions must sum to unity")
ni<-ceiling(proportions*(n-l/2))
N<-sum(ni)
while (N!=n) {
if (N<n) {
m<-ni/proportions
ee<-min(m)
}
else {
m<-(ni-1)/proportions
ee<-max(m)
}
v<-(abs(ee-m)<1e-8)
i<-(1:l)[v]
if (random) {
i<-sample(i)
}
if (N<n) {
ni[i[1]]<-ni[i[1]]+1
}
else {
found=FALSE
for (j in i) {
if (ni[j]>1) {
ni[j]<-ni[j]-1
found=TRUE
break
}
}
if (!found) # this will result in a zero frequency
ni[i[1]]<-ni[i[1]]-1
}
N<-sum(ni)
}
ni
}
"eval.blockdesign" <-
function(frml,design,blocksizes,rho=1,confounding=FALSE,center=FALSE){
# evaluates a blocked design.
# a formula
# a design
# a vector of blocksizes
# rho is a vector of variance component ratios.
# confounding=TRUE, returns the confounding matrix
# center=TRUE, the numeric columns will be centered before model expansion
rownames(design)<-1:nrow(design) # prevents warning messages in case of duplicate rownames
frml<-expand.formula(frml,colnames(design))
design<-data.frame(design)
numericColumns<-sapply(design,is.numeric)
if (center) {
means<-apply(design[,numericColumns,drop=FALSE],2,mean)
design[,numericColumns]<-sweep(design[,numericColumns,drop=FALSE],2,means)
}
des<-model.matrix(frml,design)
desBlkCtr<-des
if (missing(blocksizes))
stop("blocksizes must be specified")
if (any(rho<1e-10 || rho>1e10))
stop("Some rho's are too large or too small")
nB<-length(blocksizes)
fn<-0
S<-matrix(0,nB,ncol(des))
for (i in 1:nB) {
st<-fn+1
fn<-fn+blocksizes[i]
ad<-desBlkCtr[st:fn,]
S[i,]<-apply(ad,2,sum)
ad<-scale(ad,scale=FALSE)
desBlkCtr[st:fn,]<-ad
}
index<-apply((abs(desBlkCtr)<1e-10),2,all) # TRUE corresponds to constant or whole block vars
desUnblocked<-scale(des,scale=FALSE)
# the following combines uncenterd whole terms with centered within terms
desUnblocked<-cbind(des[,index,drop=FALSE],desUnblocked[,!index,drop=FALSE]) # rearrange columns
des<-cbind(des[,index,drop=FALSE],desBlkCtr[,!index,drop=FALSE]) # to conform to these
desBlkCtr<-cbind(desBlkCtr[,index,drop=FALSE],desBlkCtr[,!index,drop=FALSE])
S<-cbind(S[,index,drop=FALSE],S[,!index,drop=FALSE])
constantPosition<-match("(Intercept)",dimnames(des)[[2]],nomatch=0)
k.whole<-sum(index)
k.within<-sum(!index)
nB<-length(blocksizes)
wholeIndex<-c(rep(TRUE,k.whole),rep(FALSE,k.within))
N<-nrow(desUnblocked)
k<-ncol(desUnblocked)
output<-NULL
Mu<-(t(desUnblocked)%*%desUnblocked)/N
Mw<-(t(desBlkCtr[,!wholeIndex,drop=FALSE])%*%desBlkCtr[,!wholeIndex,drop=FALSE])/N
determinant<-det(Mu)
if (determinant<=0)
stop("Singular design")
determinant<-(determinant)^(1/k)
MI<-solve(Mu)
Mbc<-solve((t(des)%*%des)/N)
determinant.within<-det(Mw)
if (determinant.within<=0)
stop("Singular within design")
# If the model contains nesting, the determinant.within will differ
# from that reported by optBlock(), which uses multiplication instead
# of nesting. This is a constant shift of no importance to the optimization
determinant.within<-(determinant.within)^(1/k.within)
if (confounding == TRUE) {
Mbc<-solve((t(des)%*%des)/N)
dg<-diag(-1/diag(Mbc))
output <- list(confounding = round(Mbc%*%dg,4))
}
# calcualate the blocking efficiencies
T<-S[,wholeIndex,drop=FALSE]
S<-S[,!wholeIndex,drop=FALSE];
S<-scale(S,scale=FALSE)
Xa<-desUnblocked[,!wholeIndex,drop=TRUE]
Xa<-scale(Xa,scale=FALSE)
Xb<-desUnblocked[,wholeIndex,drop=TRUE]
M<-t(Xa)%*%Xa
Nb<-t(Xb)%*%Xb
H<-t(Xa)%*%Xb
D<-blocksizes
dM<-det(M)
trM<-sum(diag(solve(M)))
lambda<-matrix(0,3,length(rho))
lambda[1,]<-rho
colnames(lambda)<-rep(" ",length(rho))
rownames(lambda)<-c("rho","lambda.det","lambda.trace")
lambdaD<-lambdaT<-NULL
for (r in rho) {
G<-diag(1/((1/r)+D))
MM<-M-t(S)%*%G%*%S
A1<-H-t(S)%*%G%*%T
A2<-solve(Nb-t(T)%*%G%*%T)
MM<-MM-A1%*%A2%*%t(A1)
dMM<-det(MM)
lambdaD<-c(lambdaD,round((dMM/dM)^(1/k.within),3))
if (dMM>dM*1e-16) {
lambdaT<-c(lambdaT,round(trM/sum(diag(solve(MM))),3))
}
else {
lambdaT=c(lambdaT,0)
}
}
lambda[2,]<-lambdaD
lambda[3,]<-lambdaT
outs<-data.frame(constant=I(rep(" ",4)),whole=I(rep(" ",4)),within=I(rep(" ",4)))
rownames(outs)<-c("df","determinant","gmean.variance","gmean.efficiencies")
whIndex<-wholeIndex
if (constantPosition && k.whole>0) {
df<-k.whole-1
whIndex[constantPosition]<-FALSE
}
else
df<-k.whole
outs[1,]<-c(sprintf("1"),sprintf("%.0f",df),sprintf("%.0f",k.within))
outs[2,3]<-sprintf("%f",determinant.within)
output<-c(output,list(determinant.all.terms.within.terms.centered=determinant))
output<-c(output,list(within.block.efficiencies=lambda))
if (k>nB) {
output<-c(output,list(comment="Too few blocks to recover interblock information."))
}
# variances of blockcentered design
dg<-diag(Mbc)
if (k.whole>0) {
if (constantPosition) {
s<-dg[constantPosition]
outs[3,1]<-sprintf("%f",s)
}
if (sum(whIndex)>0) {
s<-dg[whIndex]
s<-exp(mean(log(s)))
outs[3,2]<-sprintf("%f",s)
}
}
t<-dg[!wholeIndex]
t<-exp(mean(log(t)))
outs[3,3]<-sprintf("%f",t)
# efficiencies of un-blockcentered to blockcentered designs
dg<-(diag(MI)/diag(Mbc))
if (k.whole>0) {
if (constantPosition) {
s<-dg[constantPosition]
outs[4,1]<-sprintf("%f",s)
}
if (sum(whIndex)>0) {
s<-dg[whIndex]
s<-exp(mean(log(s)))
outs[4,2]<-sprintf("%.3f",round(s,3))
}
}
t<-dg[!wholeIndex]
t<-exp(mean(log(t)))
outs[4,3]<-sprintf("%.3f",round(t,3))
output<-c(output,list(block.centered.properties=outs))
output
}
"eval.design" <-
function (frml, design, confounding = FALSE, variances = TRUE,
center = FALSE, X=NULL)
{
rownames(design)<-1:nrow(design) # prevents warning messages in case of duplicated rownames
proportions<-NULL
if (colnames(design)[1]=="Rep..") { # Expand data according to Rep.. from optFederov() output.
ND<-nrow(design)
reps<-design[,1]
design<-design[,-1,drop=FALSE]
design<-design[rep(1:ND,reps),]
rownames(design)<-1:nrow(design)
}
else
if (colnames(design)[1]=="Proportion") {
proportions<-design[,1]
design<-design[,-1,drop=FALSE]
}
frml <- expand.formula(frml, colnames(design))
design <- data.frame(design)
numericColumns <- sapply(design, is.numeric)
if (center) {
means <- apply(design[, numericColumns, drop = FALSE],
2, mean)
design[, numericColumns] <- sweep(design[, numericColumns,
drop = FALSE], 2, means)
}
des <- model.matrix(frml, design)
if (!missing(X)){
X<-data.frame(X)
if (center) {
means <- apply(X[, numericColumns, drop = FALSE],
2, mean)
X[, numericColumns] <- sweep(X[, numericColumns,
drop = FALSE], 2, means)
}
Xd <- model.matrix(frml, X)
}
output <- NULL
if (is.vector(proportions)) {
M<-t(des)%*%diag(proportions)%*%des
}
else {
M <- (t(des) %*% des)/nrow(des)
}
k<-ncol(M)
determinant<-det(M)
if (determinant<=0)
stop("Singular design")
determinant<-(determinant)^(1/k)
trac<-0
Ival<-0
Geff<-0
Deff<-0
orthog<-0
if (determinant>0) {
MI<-solve(M)
if (confounding == TRUE) {
dg<-diag(-1/diag(MI))
output <- list(confounding = round(MI%*%dg,4))
}
trac<-sum(diag(MI))/k
if (!missing(X)) {
CC<-t(Xd)%*%Xd/nrow(Xd)
Ival<-sum(diag(qr.solve(M,CC)))
Z<-t(qr.solve(M,t(Xd)))
Geff<-ncol(Xd)/max(apply(Z*Xd,1,sum))
Deff<-exp(1-1/Geff)
}
cp<-match("(Intercept)",dimnames(M)[[2]],nomatch=0)
if (cp>0) {
orthog<-exp((log(det(M[-cp,-cp]))-sum(log(diag(M)[-cp])))/(k-1))
}
else {
orthog<-exp(log(determinant)-sum(log(diag(M)))/k)
}
}
output <- c(output, list(determinant = determinant,A = trac))
if (!missing(X)) {
output<-c(output,list(I = Ival,Geff = round(Geff,3), Deffbound = round(Deff,3)))
}
output<-c(output,list(diagonality = round(orthog,3)))
if (variances) {
if (determinant > 0) {
t <- diag(solve(M))
if (cp>0)
t<-t[-cp]
t <- exp(mean(log(t)))
output <- c(output, list(gmean.variances = t))
}
}
output
}
"gen.factorial" <-
function(levels,nVars=0,center=TRUE,factors="none",varNames=NULL){
# Generates a full factorial design
# nVars: number of factors. May be missing if length(levels)>1
# levels: a vector of levels
# center: centers non-factorial variables
# factors: variable numbers of the factors or "all" if all are factors
# varNames: the names of the variables
if (length(levels)==1){
levels<-rep(levels,nVars)
if (missing(nVars))
stop("nVars is needed when levels is a scalar")
}
else {
nVars<-length(levels)
}
if (any(levels<1))
stop("All levels must be greater than 1")
N<-prod(levels)
factorVec<-rep(0,nVars)
if (missing(factors)) {
factors<-NULL
}
else
if (is.character(factors)) {
if (factors=="all") {
factors<-1:nVars
factorVec<-rep(1,nVars)
}
else
stop("factors error")
}
else {
nFactors<-length(factors);
for (i in factors)
factorVec[i]<-1
}
design<-matrix(0,N,nVars)
.Call("GetFactorial",design,as.integer(levels),as.integer(center),as.integer(factorVec),PACKAGE="AlgDesign")
design<-data.frame(design)
if (!missing(varNames) && length(varNames==nVars))
colnames(design)<-varNames
else
colnames(design)<-paste("X",1:nVars,sep="")
if (!missing(factors)) {
for (i in factors) {
design[,i]<-factor(design[,i])
}
}
design
}
"model.matrix.formula" <-
function(frml,data=sys.frame(sys.parent()),...){
if (!missing(data)) {
if (!inherits(data,"data.frame"))
stop("data must be a data.frame")
if (!inherits(frml,"formula"))
stop("frml must be a formula")
frml<-expand.formula(frml,colnames(data))
}
model.matrix.default(frml,data,...)
}
"model.matrix.terms"<-
function(...){
model.matrix.default(...)
}
"gen.mixture" <-
function(levels,vars) {
if (is.numeric(vars)) {
nvars<-vars
vars<-paste("X",1:nvars,sep="")
}
else {
nvars<-length(vars)
}
levels<-levels[1]-1
if (levels<1)
stop("levels error")
if (nvars<2)
stop("vars error")
N<-choose(levels+nvars-1,levels);
X<-matrix(0,N,nvars)
.Call("GetMixture", X,as.integer(levels),PACKAGE="AlgDesign")
X<-data.frame(X)
colnames(X)<-vars
X
}
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.