Automating Optimal Design for the Normal Linear Model

# design.splus  Seek a configuration of RELATIVE sample sizes that maximizes
#               power for the normal linear model. Represent one categorical
#               independent variable with p categories (possibly representing
#               combinations of several factors in a factorial design) with
#               cell means dummy variable coding. The null hypothesis is
#               C Beta = h. You specify the d x p matrix C (called CC) and the
#               d x 1 matrix C_Beta_minus_h (called CBMH), and also an
#               arbitrary total sample size n.
#
#               The function check may be used to compare sample size
#               configurations manually. This is useful because this program
#               just finds ONE good design; there may be others just as good.
#
CC <- rbind( c(1,-1,0),
             c(1,0,-1) )        # Give rows of C matrix
CBMH <-  as.matrix(c(1,1))      # C Beta Minus H is a d x 1 matrix
n <- 150
#############################################################################
# Define functions
#
dsq <- function(xx) # Returns a quantity proportional tp Minus delta-squared.
                    # Function nlminb will MINIMIZE this (That's why minus).
                    # xx is all the sample sizes except the last
                    # Assumes n, CC, CBMH
    {
    if(n < sum(xx)) stop("Sum of xx cannot exceed n")
    sizes <- c(xx,n-sum(xx))
    dsq <- -1 *t(CBMH)%*%solve(CC%*%diag(1/sizes)%*%t(CC))%*%CBMH
    dsq
    }
check <- function(sizes) # Check a vector of sample sizes. Assumes CC & CBMH
    {
    check <- as.numeric(t(CBMH)%*%solve(CC%*%diag(1/sizes)%*%t(CC))%*%CBMH)
    check # Return function value
    }
#############################################################################
# Now do calculations
#
dd <- length(CBMH) ; pp <- dim(CC)[2]
if(dd != dim(CC)[1]) stop("CC, CBMH wrong sizes")
eql <-  numeric(pp-1) + n/pp     # Equal starting values
loww <- numeric(pp-1) + .1       # Lower limit almost zero
upp  <- numeric(pp-1) + n-.1     # Upper limit almost n
best <- nlminb( start = eql, obj=dsq, lower=loww, upper=upp)$parameters
bestd <- c(best,n-sum(best))
cat("  \n")
cat("Power is maximized for relative sample sizes of \n")
cat("   ",bestd,"  \n")
cat("  \n")
cat("Delta-squared is proportional to ",check(bestd),"  \n")
cat("Use the function check to compare other designs. \n")
cat("  \n")

Now get into S-plus and use the program above.

/dos/brunner/2201 > OSplus
Warning: Cannot open audit file
S-PLUS : Copyright (c) 1988, 1996 MathSoft, Inc.
S : Copyright AT&T.
Version 3.4 Release 1 for Silicon Graphics Iris, IRIX 5.3 : 1996
Working data will be in .Data
> source("design.splus")

Power is maximized for relative sample sizes of
    75.0000249670751 50.0000000041319 24.999975028793

Delta-squared is proportional to  37.4999999999959
Use the function check to compare other designs.

> check(c(75,50,25))
[1] 37.5
> check(c(75,25,50))
[1] 37.5
> check(c(75,74,1))
[1] 37.5
> check(c(70,30,50))
[1] 37.33333
> check(c(50,50,50))
[1] 33.33333

But we were interested in what would happen, say if the difference betw control & treatment 2 were twice as big. Go into design.splus, modify CBMH <- as.matrix(c(1,2)), save, and

> source("design.splus")

Power is maximized for relative sample sizes of
    74.9500249780114 0.1 74.9499750219886

Delta-squared is proportional to  149.899999999983
Use the function check to compare other designs.


> check(c(75,0,75))
Error in .Fortran(if(!cmplx) "dqr" else "zqr",: subroutine dqr: 3 missing
        value(s) in argument 1
Dumped

Well, of course. You can't invert x`x, can you? Modify the first 2 lines to read

CC <- rbind( c(1,-1))
CBMH <-  as.matrix(2)

Save, and

> source("design.splus")

Power is maximized for relative sample sizes of
    75 75

Delta-squared is proportional to  150
Use the function check to compare other designs.

So the moral of the story seems to be, place all your eggs in the basket with the biggest effect. Is this rational?

Yes, if you're sure.