comparison of D- and I-optimal samples for a 6th order 2D response surface

idopt.R

# compare point distributions for I and D optimal designs 
# file: idopt.R
# 
rm(list=ls()) 
library(pwr) 
library(AlgDesign)
library(xtable)

dice = data.frame(var=c("x","y"), low=c(-1,-1), high=c(1,1), center=c(0,0), nLevels=111, round=4, factor=FALSE)
Ispace = gen.factorial(levels=c(111,111),nVars=2,center=TRUE,varNames=c("x","y"))/55 

dx.I = optMonteCarlo(~(x+I(x^2)+I(x^3)+I(x^4)+I(x^5)+I(x^6))*(y+I(y^2)+I(y^3)+I(y^4)+I(y^5)+I(y^6)), data=dice, nTrials=49, criterion="I", evaluateI=TRUE, space=Ispace, RandomStart=TRUE, nRepeats=10, nCand=2000)
dx.D = optMonteCarlo(~(x+I(x^2)+I(x^3)+I(x^4)+I(x^5)+I(x^6))*(y+I(y^2)+I(y^3)+I(y^4)+I(y^5)+I(y^6)), data=dice, nTrials=49, criterion="D", RandomStart=TRUE, nRepeats=10, nCand=2000)

pw = 469.75502/72.0 # inches
png("ID_points.png", width=pw, height=pw, units="in", res=150)
plot(c(-1,1),c(-1,1),type="n",xlab="x",ylab="y")
points(dx.I$design$x, dx.I$design$y, pch=1, col="darkblue")
points(dx.D$design$x, dx.D$design$y, pch=2, col="darkred")
legend(-0.4,0.4,legend=c("I","D"),pch=c(1,2),col=c("darkblue","darkred"))
title(main="Comparison of D- and I-optimal Samples", sub="supports 6th order 2D response surface")
dev.off()