Skip to content

Instantly share code, notes, and snippets.

@jstults jstults/ofat.R
Created Feb 2, 2013

What would you like to do?
why one-factor-at-a-time testing is a bad idea
# why one-factor-at-a-time testing is a bad idea
# full factorial of 6 factors at 2 levels:
fullfac = gen.factorial(levels=c(2,2,2,2,2,2), center=TRUE)
# D-optimal fraction, 32 runs:
dx.1 = optFederov(~(X1+X2+X3+X4+X5+X6)^2, data=fullfac, nTrials=32, nullify=2, maxIteration=64, nRepeats=200)
# ofat design with replication, 36 runs:
dx.2 = data.frame(rbind(diag(6),-diag(6),diag(6),-diag(6),diag(6),-diag(6)))
ntrials = 1000 # number of times to fit the model to synthetic data
std.err.1 = rep(0,ntrials) # record average standard error for each design
std.err.2 = rep(0,ntrials)
for(i in seq(ntrials)){
# add a response:
dat.1 = transform(dx.1$design, Y=X1+X2+X3+X4+X5+X6+rnorm(nrow(dx.1$design),mean=0,sd=1))
dat.2 = transform(dx.2, Y=X1+X2+X3+X4+X5+X6+rnorm(nrow(dx.2),mean=0,sd=1))
# fit the models
lm.1 = lm(Y~X1+X2+X3+X4+X5+X6, data=dat.1)
lm.2 = lm(Y~X1+X2+X3+X4+X5+X6, data=dat.2)
# average the standard errors for the model coefficients
std.err.1[i] = mean(summary(lm.1)$coefficients[,2])
std.err.2[i] = mean(summary(lm.2)$coefficients[,2])
# estimate probability densities from the trials:
std.err.dens.1 = density(std.err.1, bw="sj")
std.err.dens.2 = density(std.err.2, bw="sj")
#plots and tables:
png("ofat_doe_compare.png", width=6.4, height=0.5*6.4, units="in", res=360)
plot(c(std.err.dens.1$x,std.err.dens.2$x), c(std.err.dens.1$y,std.err.dens.2$y), type="n", xlab="average coefficient standard error", ylab="probability density")
lines(std.err.dens.1$x, std.err.dens.1$y, col="blue")
lines(std.err.dens.2$x, std.err.dens.2$y, col="red")
legend(x="topright",legend=c("D-optimal, 32 runs","OFAT, 36 runs"), col=c("blue","red"), lty=1)
print(xtable(dx.1$design,caption="D-Optimal Design",label="t:dx1",digits=0),type="latex",file="dx.1.tex",table.placement="htbp",caption.placement="top")
print(xtable(dx.2,caption="OFAT Design",label="t:dx2",digits=0),type="latex",file="dx.2.tex",table.placement="htbp",caption.placement="top")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.