R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("AlgDesign-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('AlgDesign') > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "Federov" > > ### * Federov > > flush(stderr()); flush(stdout()) > > ### Name: optFederov > ### Title: Optimal design > ### Aliases: optFederov > ### Keywords: design > > ### ** Examples > > > # EXAMPLE 1 > # A quadratic polynomial in three variables. The resulting D will be about 0.46. > # This may be compared with a standard central composite design obtained from > # rows 1,3,5,7,9,11,13,15,17,19,21,23,25,27 of dat, which has a D value of 0.46. > # The central composite design seems to be the optimal design for all three criteria. > > dat<-gen.factorial(levels=3,nVars=3,varNames=c("A","B","C")) > > desD<-optFederov(~quad(A,B,C),dat,nTrials=14,eval=TRUE) > desA<-optFederov(~quad(.),dat,nTrials=14,eval=TRUE,crit="A") > desI<-optFederov(~quad(.),dat,nTrials=14,eval=TRUE,crit="I") > > rows<-c(1,3,5,7,9,11,13,15,17,19,21,23,25,27) > desO<-optFederov(~quad(.),dat,nTrials=14,eval=TRUE,rows=rows) > > # The I criterion may be seen to decrease as the space is expanded. > > levels<-seq(-1,1,by=.1) > dat<-expand.grid(list(A=levels,B=levels,C=levels)) > > desL<-optFederov(~quad(.),dat,nTrials=14,eval=TRUE) > > # This is not the case for A or D. For A and D, the support points are the points > # of the grid with the three levels above. Points not on this grid move > # the criteria in a non-optimal direction; hence, the enlarging space has no effect. > > # EXAMPLES 2 > # Standard designs are usually optimal designs. If nTrials is set to that for > # a standard design, and if nRepeats is large enough, the standard design will > # often be found For example, a half replicate of a 2^4 will be obtained by the > # following. > > dat<-gen.factorial(levels=2,nVars=3,varNames=c("A","B","C")) > desH<-optFederov(~.,dat,8) > > # A third replicate of a 3^3 will be obtained by the following: > > dat<-gen.factorial(levels=3,nVars=3,factor=1:3) > desT<-optFederov(~.,dat,9) > > # An orthogonal design similar to a 12 run Plackett-Burman design can be > # created by the following. > > dat<-gen.factorial(levels=2,nVars=11,varNames=c("A","B","C","D","E","F","G","H","J","K","L")) > desPB<-optFederov(~.,dat,12,nRepeats=20) > > # The above calculation is numerically difficult for the A and I criteria, > # and nRepeats=100 or more may be needed. > > # It is instructive to examine a case in which the standard design is not found. > # The following is an attempt to create a Latin square design. It is not always successful. > > lv<-factor(1:5) > dat<-expand.grid(A=lv,B=lv,C=lv) > desL<-optFederov(~.,dat,nTrials=25,nRep=100) > > # It may be summarized as follows. > > cs<-xtabs(~.,desL$design) > {xx<-matrix(0,5,5); for (i in 1:5) xx=xx+cs[1:5,1:5,i]*i;xx} B A 1 2 3 4 5 1 1 5 4 2 3 2 5 0 1 4 5 3 4 2 5 3 1 4 5 4 3 1 2 5 0 4 2 5 4 attr(,"call") xtabs(formula = ~., data = desL$design) > > > # EXAMPLE 3 > # Mixture variables have a constant sum, usually 1. This causes a linear dependency > # among terms of polynomial models. In particular the constant term is dependent. > # Squared terms in a quadratic model are confounded with interaction terms, so that > # a quadratic model for three mixture variables is ~0+(A+B+C)^2. The following > # calculation generates a set of candidate varibles using gen.mixture() with > # four values on each axis, and then creates a 15 run design. The design is optimal. > # Indeed, the candidate set produced by gen.mixture(2,5) is optimal. Note: > # nullify=TRUE is used to ensure that this example will run withough error. The > # default value of 5 for nRepeats is sometimes not enought to find a starting > # design with a mixture problem. > > dat<-gen.mixture(4,5) > desM<-optFederov(~(X1+X2+X3+X4+X5)^2-1,dat,15,nullify=TRUE) > > # EXAMPLES 4 > # Design augmenation can be obtained by setting augment=TRUE, and placing the row numbers > # of the design to be agmented in rows. Augmentation is often used to (1) add a new variable > # to an existing design or (2) to increase the complexity of the model. The following illustrates > # adding a variable to an existing design using desD above. It is assumed that all runs of the > # existing design have been made at the -1 level of the new variable: > > dat<-gen.factorial(levels=3,nVars=3,varNames=c("A","B","C")) > desA<-optFederov(~quad(.),dat,nTrials=25,augment=TRUE,rows=desD$rows) > > # The half fraction in desH, can be augmented to support an additional term: > > dat<-gen.factorial(levels=2,nVars=4,varNames=c("A","B","C","D")) > desH<-optFederov(~.,dat,8) > desH2<-optFederov(~A+B+C+D+I(A*B),dat,10,aug=TRUE,rows=desH$rows) > > # EXAMPLES 5 > # Optimal approximate theory designs have non-zero probabilities only on support points. > # For the first example above the approximate theory design is as follows. It shows > # that all points in the cubic lattice are support points. The D for this > # design is 0.474 which may be compared with the D of 0.463 of the first example to > # indicate that that exact design had a D-efficiency of 97%. The lower bound Dea > # was 82%. > > dat<-gen.factorial(levels=3,nVars=3,varNames=c("A","B","C")) > desDA<-optFederov(~quad(A,B,C),dat,eval=TRUE,approx=TRUE) > > # The largest proportions will be rounded if nTrials is specified. > > desDAN<-optFederov(~quad(A,B,C),dat,eval=TRUE,approx=TRUE,nTrials=15) > > > > cleanEx(); ..nameEx <- "MonteOpt" > > ### * MonteOpt > > flush(stderr()); flush(stdout()) > > ### Name: optMonteCarlo > ### Title: Optimal design via Monte Carlo > ### Aliases: optMonteCarlo > ### Keywords: design > > ### ** Examples > > > # EXAMPLE 1 > # The data.frame in data might look like the following: > data<-data.frame(var=paste("X",1:6,sep=""),low=c(1,1,1,0,0,0), + high=c(3,3,3,1,1,1),center=c(2,2,2,0,0,0),nLevels=3, + round=1,factor=0,mix=c(FALSE,FALSE,FALSE,TRUE,TRUE,TRUE)) > data var low high center nLevels round factor mix 1 X1 1 3 2 3 1 0 FALSE 2 X2 1 3 2 3 1 0 FALSE 3 X3 1 3 2 3 1 0 FALSE 4 X4 0 1 0 3 1 0 TRUE 5 X5 0 1 0 3 1 0 TRUE 6 X6 0 1 0 3 1 0 TRUE > > # and the design: > > optMonteCarlo(~(X1+X2+X3)^2+X4+X5+X6,data) $D [1] 0.4690664 $A [1] 32.82362 $Ge [1] 0.738 $Dea [1] 0.701 $design X1 X2 X3 X4 X5 X6 1 3 3 3 0.79 0.13 0.08 2 3 1 2 0.71 0.26 0.03 3 3 1 3 0.22 0.09 0.69 4 1 1 1 0.63 0.11 0.26 5 3 3 1 0.08 0.31 0.61 6 1 1 2 0.48 0.46 0.06 7 1 3 1 0.20 0.78 0.02 8 2 1 3 0.28 0.44 0.28 9 1 2 1 0.04 0.01 0.95 10 3 3 3 0.11 0.50 0.39 11 1 3 3 0.00 0.89 0.11 12 1 3 3 0.05 0.02 0.93 13 3 1 1 0.08 0.56 0.36 14 3 3 1 0.96 0.03 0.01 > > # Example 2 > # Standard designs will often be produced, just as > # they will with optFederov(). For example, > # a half fraction of a 2^4: > data<-data.frame(paste("X",1:4,sep=""),-1,1,0,2,0,0) > data paste..X...1.4..sep...... X.1 X1 X0 X2 X0.1 X0.2 1 X1 -1 1 0 2 0 0 2 X2 -1 1 0 2 0 0 3 X3 -1 1 0 2 0 0 4 X4 -1 1 0 2 0 0 > optMonteCarlo(~.,data,nTrials=8) $D [1] 1 $A [1] 1 $Ge [1] 1 $Dea [1] 1 $design X1 X2 X3 X4 1 -1 -1 -1 -1 2 -1 -1 1 -1 3 -1 1 1 1 4 -1 1 -1 1 5 1 -1 -1 1 6 1 1 -1 -1 7 1 -1 1 1 8 1 1 1 -1 > > # Example 3 > # optMonteCarlo() can treat much larger problems than can > # optFederov(). For example, optFederov() > # requires a candidate list of 3^20 points for > # a 20 variable, 3 level candidate list -- about > # 25 gigabytes. If the model is quadratic, this must > # be multiplied by about 12. There are other storage > # requirements internal to optFederov() which easily > # double this value. optMonteCarlo() since it only samples > # from the putative candidate list, has no difficulty > # with a problem of this size. The criterion values > # appearing in the output of optMonteCarlo() are based on > # these samples, but their values seem to be reasonable > # correct, as the following shows: (These are commented > # out for those who have a slow machine.) > > dat<-gen.factorial(levels=3,nVar=8) > #desF<-optFederov(~quad(.),dat,eval=TRUE) > #desF[1:5] > > data<-data.frame(paste("X",1:8,sep=""),-1,1,0,3,0,0) > #desH<-optMonteCarlo(~quad(.),data,Rand=FALSE,eval=TRUE) > #desH[1:5] > > # The following is a 20 variable quadratic. Uncomment > # and wait a while, even if you have a fast machine. > # Note: nRepeats has been changed from its default. > # Note: criterion values for exact designs are often > # far from approximate theory optima; hence, Ge and De > # will be small. > > data<-data.frame(paste("X",1:20,sep=""),-1,1,0,3,0,0) > #desBig<-optMonteCarlo(~quad(.),data,nRepeats=1) > > # The following will produce improved criterion values > > #desNBig<-optMonteCarlo(~quad(.),data,Rand=FALSE,nRepeats=1) > > # EXAMPLE 4 > # Practically infeasible combinations of variable are > # common. Designs may be produced which avoid such > # combinations by using a constraint function. Suppose, > # for example that one corner of a cubic box is not > # feasible, then the following will produce a design > # that makes no use of this corner. > > Constraints<-function(x){!(x[1]>0.75 && x[2]>0.75)} > data<-data.frame(paste("X",1:4,sep=""),-1,1,0,3,0,0) > desC<-optMonteCarlo(~.,data,con=Constraints) > > # The above just removes a corner. Increasing the > # number of levels will remove points along the > # boundary. > > data<-data.frame(paste("X",1:4,sep=""),-1,1,0,11,3,0) > desC2<-optMonteCarlo(~.,data,con=Constraints) > > > > > cleanEx(); ..nameEx <- "blockOpt" > > ### * blockOpt > > flush(stderr()); flush(stdout()) > > ### Name: optBlock > ### Title: Optimal design blocking > ### Aliases: optBlock > ### Keywords: design > > ### ** Examples > > # Blocking the design for a quadratic polynomial in three variables into two > # seven trial blocks: > > dat<-gen.factorial(3,3,varNames=c("A","B","C")) > desD<-optFederov(~quad(.),dat,nTrials=14,eval=TRUE) # Choose an optimum 14 trail design. > optBlock(~quad(.),desD$design,c(7,7)) $D [1] 0.4178207 $diagonality [1] 0.969 $Blocks $Blocks$B1 A B C 1 -1 -1 -1 7 -1 1 -1 9 1 1 -1 12 1 -1 0 20 0 -1 1 24 1 0 1 25 -1 1 1 $Blocks$B2 A B C 3 1 -1 -1 5 0 0 -1 13 -1 0 0 17 0 1 0 19 -1 -1 1 21 1 -1 1 27 1 1 1 $design A B C 1 -1 -1 -1 7 -1 1 -1 9 1 1 -1 12 1 -1 0 20 0 -1 1 24 1 0 1 25 -1 1 1 3 1 -1 -1 5 0 0 -1 13 -1 0 0 17 0 1 0 19 -1 -1 1 21 1 -1 1 27 1 1 1 $rows [1] 1 7 9 12 20 24 25 3 5 13 17 19 21 27 > > # Letting optBlock() search the dat candidate list instead of first choosing a > # 14 trial design. > optBlock(~quad(.),dat,c(7,7)) $D [1] 0.42645 $diagonality [1] 0.927 $Blocks $Blocks$B1 A B C 1 -1 -1 -1 3 1 -1 -1 6 1 0 -1 7 -1 1 -1 9 1 1 -1 21 1 -1 1 25 -1 1 1 $Blocks$B2 A B C 2 0 -1 -1 4 -1 0 -1 8 0 1 -1 16 -1 1 0 18 1 1 0 19 -1 -1 1 27 1 1 1 $design A B C 1 -1 -1 -1 3 1 -1 -1 6 1 0 -1 7 -1 1 -1 9 1 1 -1 21 1 -1 1 25 -1 1 1 2 0 -1 -1 4 -1 0 -1 8 0 1 -1 16 -1 1 0 18 1 1 0 19 -1 -1 1 27 1 1 1 $rows [1] 1 3 6 7 9 21 25 2 4 8 16 18 19 27 > > # A block design for 7 treatments in 7 blocks of size 3. Note how withinData > # is recycled to fill out the blocksize requirements. > > BIB<-optBlock(~.,withinData=factor(1:7),blocksizes=rep(3,7)) > > # This is a balanced incomplete block design as may be seen from: > > crossprod(table(BIB$rows,c(rep(1:7, rep(3,7))))) 1 2 3 4 5 6 7 1 3 1 1 1 1 1 1 2 1 3 1 1 1 1 1 3 1 1 3 1 1 1 1 4 1 1 1 3 1 1 1 5 1 1 1 1 3 1 1 6 1 1 1 1 1 3 1 7 1 1 1 1 1 1 3 > > # A partially balanced incomplete block design with two associate classes: > > tr<-factor(1:9) > PBIB<-optBlock(~.,withinData=tr,blocksizes=rep(3,9)) > > crossprod(table(PBIB$rows,c(rep(1:9, rep(3,9))))) 1 2 3 4 5 6 7 8 9 1 3 1 0 1 1 0 1 1 1 2 1 3 1 0 1 1 0 1 1 3 0 1 3 1 1 0 1 1 1 4 1 0 1 3 1 1 0 1 1 5 1 1 1 1 3 1 1 0 0 6 0 1 0 1 1 3 1 1 1 7 1 0 1 0 1 1 3 1 1 8 1 1 1 1 0 1 1 3 0 9 1 1 1 1 0 1 1 0 3 > > # Two fractions of a 2^(4-1). > > dat<-gen.factorial(2,4) > od<-optBlock(~.,dat,c(8,8)) > > # The blocks are not themselves orthogonal even though the entire design is optimal. > > bk<-data.matrix(od$Blocks$B1) > t(bk)%*%bk X1 X2 X3 X4 X1 8 0 -4 0 X2 0 8 4 0 X3 -4 4 8 -4 X4 0 0 -4 8 > > # Better blocks may be obtained as follows, but note that they are not generally > # the fractions that would be obtained by confounding the third order interaction. > > od<-optBlock(~.,dat,c(8,8),criterion="Dpc",nR=10) > bk<-data.matrix(od$Blocks$B1) > t(bk)%*%bk X1 X2 X3 X4 X1 8 0 0 0 X2 0 8 0 0 X3 0 0 8 0 X4 0 0 0 8 > > # Blocking with whole plot factors. Note that the 27 rows of within are recycled > # to make the 54 trial blocked design. > > within<-expand.grid(A=c(-1,0,1),B=c(-1,0,1),C=c(-1,0,1)) > whole<-expand.grid(D=factor(1:3),E=factor(1:3)) > od<-optBlock(~D+E*(quad(A,B,C)),withinData=within,blocksizes=rep(6,9),wholeBlockData=whole) > > # Either withinData, or wholeBlockData may be an approximate theory optimial design > # produced by optFederov() for nTrials. The first column in the optFederov() output > # design, named "Rep..", is used to replicate the trials. > > within<-optFederov(~quad(A,B,C),within,nT=54,approx=TRUE) > od<-optBlock(~D+E*(quad(A,B,C)),withinData=within$design,blocksizes=rep(6,9),wholeBlockData=whole) > > > > > cleanEx(); ..nameEx <- "genFact" > > ### * genFact > > flush(stderr()); flush(stdout()) > > ### Name: gen.factorial > ### Title: Generates a full factorial design > ### Aliases: gen.factorial > ### Keywords: design > > ### ** Examples > > dat<-gen.factorial(3,3) > dat<-gen.factorial(c(3,2,3)) > dat<-gen.factorial(3,3,factors="all") > dat<-gen.factorial(3,3,varNames=c("A","B","C")) > > > > > ### *