Nothing
R/Compare2FDS.R
In Vdgraph: Variance Dispersion Graphs and Fraction of Design Space Plots for Response Surface Designs
Defines functions
Compare2FDS
Documented in
Compare2FDS
Compare2FDS <- function(des1, des2, name1, name2, mod=2)
{
#####################################################################
# This function creates a fraction of design space plot for a design
# stored in the matrix or data frame des
# Inputs;
# des1 - a matrix or data frame containing a experimental design
# in coded or uncoded units. There should be one column
# for each column in the matrix and one row for each run
# in the design. The maximum number of columns allowed is
# 7.
# des2 - a matrix or data frame containing a experimental design
# in coded or uncoded units. There should be one column
# for each column in the matrix and one row for each run
# in the design. The maximum number of columns allowed is
# 7.
#
# mod - the model to be represented:
# 0 = linear model
# 1 = linear main effects plus linear by linear 2-factor
# interactions
# 2 = full quadratic response surface model (default)
####################################################################
#
# First check calculate the coordinates for des1
####################################################################
#First check to see if des1 is a matrix, and if so convert it to a data frame
if(is.matrix(des1)) des1 <- as.data.frame(des1)
#Next check to see if des1 is centered, and if not center it
chkm <-apply(des1, 2, f)
if (mean(chkm) !=0) des1<-scale(des1, scale=FALSE)
#Check to see if scaled, and if not scale it
chkm <-apply(des1, 2, mx)
if(mean(chkm) !=1) des1<-scale(des1, center=FALSE, scale=chkm)
rc <-dim(des1)
ndpts<-rc[1]
kvar1<-rc[2]
if (kvar1<2) {
stop("The minimum number of independent variables must be 2","\n")
}
if (kvar1>7) {
stop("The maximum number of independent variables allowed is 7","\n")
}
des1<-as.data.frame(des1)
# get the model
y<-stats::runif(nrow(des1),0,1)
### Case for two variables #############
if(kvar1==2) {
names<-list("x1","x2")
names(des1)<-names
lmod<-stats::lm(y~x1+x2,data=des1)
imod<-stats::lm(y~x1+x2+x1:x2,data=des1)
qmod<-stats::lm(y~x1+x2+I(x1^2)+I(x2^2)+x1:x2,data=des1)
if(mod==0){
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2)
}
if(mod==1){
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X1X2<-X1*X2
fX<-cbind(Int,X1,X2,X1X2)
}
if(mod==2){
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X1X2<-X1*X2
fX<-cbind(Int,X1,X2,X1sq,X2sq,X1X2)
}
}
# Case for 3 variables ########################
if(kvar1==3) {
names<-list("x1","x2","x3")
names(des1)<-names
lmod<-stats::lm(y~x1+x2+x3,data=des1)
imod<-stats::lm(y~x1+x2+x3+x1:x2+x1:x3+x2:x3,data=des1)
qmod<-stats::lm(y~x1+x2+x3+I(x1^2)+I(x2^2)+I(x3^2)+x1:x2+x1:x3+x2:x3,data=des1)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X2X3<-X2*X3
fX<-cbind(Int,X1,X2,X3,X1X2,X1X3,X2X3)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X1X2<-X1*X2
X1X3<-X1*X3
X2X3<-X2*X3
fX<-cbind(Int,X1,X2,X3,X1sq,X2sq,X3sq,X1X2,X1X3,X2X3)
}
}
# Case for 4 variables ########################
if(kvar1==4) {
names<-list("x1","x2","x3","x4")
names(des1)<-names
lmod<-stats::lm(y~x1+x2+x3+x4,data=des1)
imod<-stats::lm(y~x1+x2+x3+x4+x1:x2+x1:x3+x1:x4+x2:x3+x2:x4+x3:x4,data=des1)
qmod<-stats::lm(y~x1+x2+x3+x4+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+x1:x2+x1:x3+x1:x4+x2:x3+x2:x4+x3:x4,data=des1)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X2X3<-X2*X3
X2X4<-X2*X4
X3X4<-X3*X4
fX<-cbind(Int,X1,X2,X3,X4,X1X2,X1X3,X1X4,X2X3,X2X4,X3X4)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X2X3<-X2*X3
X2X4<-X2*X4
X3X4<-X3*X4
fX<-cbind(Int,X1,X2,X3,X4,X1sq,X2sq,X3sq,X4sq,X1X2,X1X3,X1X4,X2X3,X2X4,X3X4)
}
}
# Case for 5 variables ########################
if(kvar1==5) {
names<-list("x1","x2","x3","x4","x5")
names(des1)<-names
lmod<-stats::lm(y~x1+x2+x3+x4+x5,data=des1)
imod<-stats::lm(y~x1+x2+x3+x4+x5+x1:x2+x1:x3+x1:x4+x1:x5+x2:x3+x2:x4+x2:x5+x3:x4+x3:x5+x4:x5,data=des1)
qmod<-stats::lm(y~x1+x2+x3+x4+x5+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+I(x5^2)+x1:x2+x1:x3+x1:x4+x1:x5+x2:x3+x2:x4+x2:x5+x3:x4+x3:x5+x4:x5,data=des1)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4,X5)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X3X4<-X3*X4
X3X5<-X3*X5
X4X5<-X4*X5
fX<-cbind(Int,X1,X2,X3,X4,+X5+X1X2,X1X3,X1X4,X1X5,X2X3,X2X4,X2X5,X3X4,X3X5,X4X5)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X5sq<-X5*X5
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X3X4<-X3*X4
X3X5<-X3*X5
X4X5<-X4*X5
fX<-cbind(Int,X1,X2,X3,X4,X5,X1sq,X2sq,X3sq,X4sq,X5sq,X1X2,X1X3,X1X4,X1X5,X2X3,X2X4,X2X5,X3X4,X3X5,X4X5)
}
}
# Case for 6 variables ########################
if(kvar1==6) {
names<-list("x1","x2","x3","x4","x5","x6")
names(des1)<-names
lmod<-stats::lm(y~x1+x2+x3+x4+x5+x6,data=des1)
imod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x2:x3+x2:x4+x2:x5+x2:x6+x3:x4+x3:x5+x3:x6+x4:x5+x4:x6+x5:x6,data=des1)
qmod<-stats::lm(y~x1+x2+x3+x4+x5+x6+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+I(x5^2)+I(x6^2)+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x2:x3+x2:x4+x2:x5+x2:x6+x3:x4+x3:x5+x3:x6+x4:x5+x4:x6+x5:x6,data=des1)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4,X5,X6)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X4X5<-X4*X5
X4X6<-X4*X6
X5X6<-X5*X6
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X1X2,X1X3,X1X4,X1X5,X1X6,X2X3,X2X4,X2X5,X2X6,X3X4,X3X5,X3X6,X4X5,X4X6,X5X6)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X5sq<-X5*X5
X6sq<-X6*X6
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X4X5<-X4*X5
X4X6<-X4*X6
X5X6<-X5*X6
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X1sq,X2sq,X3sq,X4sq,X5sq,X6sq,X1X2,X1X3,X1X4,X1X5,X1X6,X2X3,X2X4,X2X5,X2X6,X3X4,X3X5,X3X6,X4X5,X4X6,X5X6)
}
}
# Case for 7 variables ########################
if(kvar1==7) {
names<-list("x1","x2","x3","x4","x5","x6","x7")
names(des1)<-names
lmod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x7,data=des1)
imod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x7+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x1:x7+x2:x3+x2:x4+x2:x5+x2:x6+x2:x7+x3:x4+x3:x5+x3:x6+x3:x7+x4:x5+x4:x6+x4:x7+x5:x6+x5:x7+x6:x7,data=des1)
qmod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x7+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+I(x5^2)+I(x6^2)+I(x7^2)+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x1:x7+x2:x3+x2:x4+x2:x5+x2:x6+x2:x7+x3:x4+x3:x5+x3:x6+x3:x7+x4:x5+x4:x6+x4:x7+x5:x6+x5:x7+x6:x7,data=des1)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X7<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X7)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X7<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X1X7<-X1*X7
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X2X7<-X2*X7
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X3X7<-X3*X7
X4X5<-X4*X5
X4X6<-X4*X6
X4X7<-X4*X7
X5X6<-X5*X6
X5X7<-X5*X7
X6X7<-X6*X7
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X7,X1X2,X1X3,X1X4,X1X5,X1X6,X1X7,X2X3,X2X4,X2X5,X2X6,X2X7,X3X4,X3X5,X3X6,X3X7,X4X5,X4X6,X4X7,X5X6,X5X7,X6X7)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X7<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X5sq<-X5*X5
X6sq<-X6*X6
X7sq<-X7*X7
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X1X7<-X1*X7
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X2X7<-X2*X7
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X3X7<-X3*X7
X4X5<-X4*X5
X4X6<-X4*X6
X4X7<-X4*X7
X5X6<-X5*X6
X5X7<-X5*X7
X6X7<-X6*X7
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X7,X1sq,X2sq,X3sq,X4sq,X5sq,X6sq,X7sq,X1X2,X1X3,X1X4,X1X5,X1X6,X1X7,X2X3,X2X4,X2X5,X2X6,X2X7,X3X4,X3X5,X3X6,X3X7,X4X5,X4X6,X4X7,X5X6,X5X7,X6X7)
}
}
# get v1 and v1i
v1<-diag(fX%*%XtXI%*%t(fX))
v1i<-sort(v1)
# calculate fraction of design space
fs1<-(1:length(v1i))/length(v1i)
####################################################################
#
# Next check calculate the coordinates for des2
####################################################################
#First check to see if des2 is a matrix, and if so convert it to a data frame
if(is.matrix(des2)) des2 <- as.data.frame(des2)
#Next check to see if des2 is centered, and if not center it
chkm <-apply(des2, 2, f)
if (mean(chkm) !=0) des2<-scale(des2, scale=FALSE)
#Check to see if scaled, and if not scale it
chkm <-apply(des2, 2, mx)
if(mean(chkm) !=1) des2<-scale(des2, center=FALSE, scale=chkm)
rc <-dim(des2)
ndpts<-rc[1]
kvar1<-rc[2]
if (kvar1<2) {
stop("The minimum number of independent variables must be 2","\n")
}
if (kvar1>7) {
stop("The maximum number of independent variables allowed is 7","\n")
}
des2<-as.data.frame(des2)
# get the model
y<-stats::runif(nrow(des2),0,1)
### Case for two variables #############
if(kvar1==2) {
names<-list("x1","x2")
names(des2)<-names
lmod<-stats::lm(y~x1+x2,data=des2)
imod<-stats::lm(y~x1+x2+x1:x2,data=des2)
qmod<-stats::lm(y~x1+x2+I(x1^2)+I(x2^2)+x1:x2,data=des2)
if(mod==0){
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2)
}
if(mod==1){
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X1X2<-X1*X2
fX<-cbind(Int,X1,X2,X1X2)
}
if(mod==2){
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X1X2<-X1*X2
fX<-cbind(Int,X1,X2,X1sq,X2sq,X1X2)
}
}
# Case for 3 variables ########################
if(kvar1==3) {
names<-list("x1","x2","x3")
names(des2)<-names
lmod<-stats::lm(y~x1+x2+x3,data=des2)
imod<-stats::lm(y~x1+x2+x3+x1:x2+x1:x3+x2:x3,data=des2)
qmod<-stats::lm(y~x1+x2+x3+I(x1^2)+I(x2^2)+I(x3^2)+x1:x2+x1:x3+x2:x3,data=des2)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X2X3<-X2*X3
fX<-cbind(Int,X1,X2,X3,X1X2,X1X3,X2X3)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X1X2<-X1*X2
X1X3<-X1*X3
X2X3<-X2*X3
fX<-cbind(Int,X1,X2,X3,X1sq,X2sq,X3sq,X1X2,X1X3,X2X3)
}
}
# Case for 4 variables ########################
if(kvar1==4) {
names<-list("x1","x2","x3","x4")
names(des2)<-names
lmod<-stats::lm(y~x1+x2+x3+x4,data=des2)
imod<-stats::lm(y~x1+x2+x3+x4+x1:x2+x1:x3+x1:x4+x2:x3+x2:x4+x3:x4,data=des2)
qmod<-stats::lm(y~x1+x2+x3+x4+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+x1:x2+x1:x3+x1:x4+x2:x3+x2:x4+x3:x4,data=des2)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X2X3<-X2*X3
X2X4<-X2*X4
X3X4<-X3*X4
fX<-cbind(Int,X1,X2,X3,X4,X1X2,X1X3,X1X4,X2X3,X2X4,X3X4)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X2X3<-X2*X3
X2X4<-X2*X4
X3X4<-X3*X4
fX<-cbind(Int,X1,X2,X3,X4,X1sq,X2sq,X3sq,X4sq,X1X2,X1X3,X1X4,X2X3,X2X4,X3X4)
}
}
# Case for 5 variables ########################
if(kvar1==5) {
names<-list("x1","x2","x3","x4","x5")
names(des2)<-names
lmod<-stats::lm(y~x1+x2+x3+x4+x5,data=des2)
imod<-stats::lm(y~x1+x2+x3+x4+x5+x1:x2+x1:x3+x1:x4+x1:x5+x2:x3+x2:x4+x2:x5+x3:x4+x3:x5+x4:x5,data=des2)
qmod<-stats::lm(y~x1+x2+x3+x4+x5+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+I(x5^2)+x1:x2+x1:x3+x1:x4+x1:x5+x2:x3+x2:x4+x2:x5+x3:x4+x3:x5+x4:x5,data=des2)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4,X5)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X3X4<-X3*X4
X3X5<-X3*X5
X4X5<-X4*X5
fX<-cbind(Int,X1,X2,X3,X4,+X5+X1X2,X1X3,X1X4,X1X5,X2X3,X2X4,X2X5,X3X4,X3X5,X4X5)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X5sq<-X5*X5
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X3X4<-X3*X4
X3X5<-X3*X5
X4X5<-X4*X5
fX<-cbind(Int,X1,X2,X3,X4,X5,X1sq,X2sq,X3sq,X4sq,X5sq,X1X2,X1X3,X1X4,X1X5,X2X3,X2X4,X2X5,X3X4,X3X5,X4X5)
}
}
# Case for 6 variables ########################
if(kvar1==6) {
names<-list("x1","x2","x3","x4","x5","x6")
names(des2)<-names
lmod<-stats::lm(y~x1+x2+x3+x4+x5+x6,data=des2)
imod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x2:x3+x2:x4+x2:x5+x2:x6+x3:x4+x3:x5+x3:x6+x4:x5+x4:x6+x5:x6,data=des2)
qmod<-stats::lm(y~x1+x2+x3+x4+x5+x6+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+I(x5^2)+I(x6^2)+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x2:x3+x2:x4+x2:x5+x2:x6+x3:x4+x3:x5+x3:x6+x4:x5+x4:x6+x5:x6,data=des2)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4,X5,X6)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X4X5<-X4*X5
X4X6<-X4*X6
X5X6<-X5*X6
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X1X2,X1X3,X1X4,X1X5,X1X6,X2X3,X2X4,X2X5,X2X6,X3X4,X3X5,X3X6,X4X5,X4X6,X5X6)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X5sq<-X5*X5
X6sq<-X6*X6
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X4X5<-X4*X5
X4X6<-X4*X6
X5X6<-X5*X6
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X1sq,X2sq,X3sq,X4sq,X5sq,X6sq,X1X2,X1X3,X1X4,X1X5,X1X6,X2X3,X2X4,X2X5,X2X6,X3X4,X3X5,X3X6,X4X5,X4X6,X5X6)
}
}
# Case for 7 variables ########################
if(kvar1==7) {
names<-list("x1","x2","x3","x4","x5","x6","x7")
names(des2)<-names
lmod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x7,data=des2)
imod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x7+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x1:x7+x2:x3+x2:x4+x2:x5+x2:x6+x2:x7+x3:x4+x3:x5+x3:x6+x3:x7+x4:x5+x4:x6+x4:x7+x5:x6+x5:x7+x6:x7,data=des2)
qmod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x7+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+I(x5^2)+I(x6^2)+I(x7^2)+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x1:x7+x2:x3+x2:x4+x2:x5+x2:x6+x2:x7+x3:x4+x3:x5+x3:x6+x3:x7+x4:x5+x4:x6+x4:x7+x5:x6+x5:x7+x6:x7,data=des2)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X7<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X7)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X7<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X1X7<-X1*X7
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X2X7<-X2*X7
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X3X7<-X3*X7
X4X5<-X4*X5
X4X6<-X4*X6
X4X7<-X4*X7
X5X6<-X5*X6
X5X7<-X5*X7
X6X7<-X6*X7
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X7,X1X2,X1X3,X1X4,X1X5,X1X6,X1X7,X2X3,X2X4,X2X5,X2X6,X2X7,X3X4,X3X5,X3X6,X3X7,X4X5,X4X6,X4X7,X5X6,X5X7,X6X7)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X7<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X5sq<-X5*X5
X6sq<-X6*X6
X7sq<-X7*X7
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X1X7<-X1*X7
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X2X7<-X2*X7
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X3X7<-X3*X7
X4X5<-X4*X5
X4X6<-X4*X6
X4X7<-X4*X7
X5X6<-X5*X6
X5X7<-X5*X7
X6X7<-X6*X7
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X7,X1sq,X2sq,X3sq,X4sq,X5sq,X6sq,X7sq,X1X2,X1X3,X1X4,X1X5,X1X6,X1X7,X2X3,X2X4,X2X5,X2X6,X2X7,X3X4,X3X5,X3X6,X3X7,X4X5,X4X6,X4X7,X5X6,X5X7,X6X7)
}
}
# get v2 and v2i
v2<-diag(fX%*%XtXI%*%t(fX))
v2i<-sort(v2)
# calculate fraction of design space
fs2<-(1:length(v2i))/length(v2i)
# make Comparison FDS Plots
maxv<-max(max(v1i),max(v2i))
plot(fs1,v1i,type='l',ylim=c(0.0,(maxv+.15)),xlab="Fraction of Space",ylab="Relative Prediction Variance")
graphics::abline(v=.5,lty=3)
graphics::abline(h=(v1i[2500]+v1i[2501])/2,lty=3)
graphics::lines(fs2,v2i,col="red",lty=2)
graphics::abline(h=(v2i[2500]+v2i[2501])/2,lty=3)
graphics::lines(c(0,.10),c(maxv,maxv),lty=1,col="black")
graphics::text(.11,maxv,name1,adj=-.1)
graphics::lines(c(0,.10),c(maxv-maxv/20,maxv-maxv/20),lty=3,col="red")
graphics::text(.11,maxv-maxv/20,name2,col="red",adj=-.1)
}
Try the Vdgraph package in your browser
Any scripts or data that you put into this service are public.
Vdgraph documentation built on Sept. 9, 2023, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions Compare2FDS
Documented in Compare2FDS
Compare2FDS <- function(des1, des2, name1, name2, mod=2)
{
#####################################################################
# This function creates a fraction of design space plot for a design
# stored in the matrix or data frame des
# Inputs;
# des1 - a matrix or data frame containing a experimental design
# in coded or uncoded units. There should be one column
# for each column in the matrix and one row for each run
# in the design. The maximum number of columns allowed is
# 7.
# des2 - a matrix or data frame containing a experimental design
# in coded or uncoded units. There should be one column
# for each column in the matrix and one row for each run
# in the design. The maximum number of columns allowed is
# 7.
#
# mod - the model to be represented:
# 0 = linear model
# 1 = linear main effects plus linear by linear 2-factor
# interactions
# 2 = full quadratic response surface model (default)
####################################################################
#
# First check calculate the coordinates for des1
####################################################################
#First check to see if des1 is a matrix, and if so convert it to a data frame
if(is.matrix(des1)) des1 <- as.data.frame(des1)
#Next check to see if des1 is centered, and if not center it
chkm <-apply(des1, 2, f)
if (mean(chkm) !=0) des1<-scale(des1, scale=FALSE)
#Check to see if scaled, and if not scale it
chkm <-apply(des1, 2, mx)
if(mean(chkm) !=1) des1<-scale(des1, center=FALSE, scale=chkm)
rc <-dim(des1)
ndpts<-rc[1]
kvar1<-rc[2]
if (kvar1<2) {
stop("The minimum number of independent variables must be 2","\n")
}
if (kvar1>7) {
stop("The maximum number of independent variables allowed is 7","\n")
}
des1<-as.data.frame(des1)
# get the model
y<-stats::runif(nrow(des1),0,1)
### Case for two variables #############
if(kvar1==2) {
names<-list("x1","x2")
names(des1)<-names
lmod<-stats::lm(y~x1+x2,data=des1)
imod<-stats::lm(y~x1+x2+x1:x2,data=des1)
qmod<-stats::lm(y~x1+x2+I(x1^2)+I(x2^2)+x1:x2,data=des1)
if(mod==0){
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2)
}
if(mod==1){
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X1X2<-X1*X2
fX<-cbind(Int,X1,X2,X1X2)
}
if(mod==2){
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X1X2<-X1*X2
fX<-cbind(Int,X1,X2,X1sq,X2sq,X1X2)
}
}
# Case for 3 variables ########################
if(kvar1==3) {
names<-list("x1","x2","x3")
names(des1)<-names
lmod<-stats::lm(y~x1+x2+x3,data=des1)
imod<-stats::lm(y~x1+x2+x3+x1:x2+x1:x3+x2:x3,data=des1)
qmod<-stats::lm(y~x1+x2+x3+I(x1^2)+I(x2^2)+I(x3^2)+x1:x2+x1:x3+x2:x3,data=des1)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X2X3<-X2*X3
fX<-cbind(Int,X1,X2,X3,X1X2,X1X3,X2X3)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X1X2<-X1*X2
X1X3<-X1*X3
X2X3<-X2*X3
fX<-cbind(Int,X1,X2,X3,X1sq,X2sq,X3sq,X1X2,X1X3,X2X3)
}
}
# Case for 4 variables ########################
if(kvar1==4) {
names<-list("x1","x2","x3","x4")
names(des1)<-names
lmod<-stats::lm(y~x1+x2+x3+x4,data=des1)
imod<-stats::lm(y~x1+x2+x3+x4+x1:x2+x1:x3+x1:x4+x2:x3+x2:x4+x3:x4,data=des1)
qmod<-stats::lm(y~x1+x2+x3+x4+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+x1:x2+x1:x3+x1:x4+x2:x3+x2:x4+x3:x4,data=des1)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X2X3<-X2*X3
X2X4<-X2*X4
X3X4<-X3*X4
fX<-cbind(Int,X1,X2,X3,X4,X1X2,X1X3,X1X4,X2X3,X2X4,X3X4)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X2X3<-X2*X3
X2X4<-X2*X4
X3X4<-X3*X4
fX<-cbind(Int,X1,X2,X3,X4,X1sq,X2sq,X3sq,X4sq,X1X2,X1X3,X1X4,X2X3,X2X4,X3X4)
}
}
# Case for 5 variables ########################
if(kvar1==5) {
names<-list("x1","x2","x3","x4","x5")
names(des1)<-names
lmod<-stats::lm(y~x1+x2+x3+x4+x5,data=des1)
imod<-stats::lm(y~x1+x2+x3+x4+x5+x1:x2+x1:x3+x1:x4+x1:x5+x2:x3+x2:x4+x2:x5+x3:x4+x3:x5+x4:x5,data=des1)
qmod<-stats::lm(y~x1+x2+x3+x4+x5+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+I(x5^2)+x1:x2+x1:x3+x1:x4+x1:x5+x2:x3+x2:x4+x2:x5+x3:x4+x3:x5+x4:x5,data=des1)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4,X5)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X3X4<-X3*X4
X3X5<-X3*X5
X4X5<-X4*X5
fX<-cbind(Int,X1,X2,X3,X4,+X5+X1X2,X1X3,X1X4,X1X5,X2X3,X2X4,X2X5,X3X4,X3X5,X4X5)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X5sq<-X5*X5
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X3X4<-X3*X4
X3X5<-X3*X5
X4X5<-X4*X5
fX<-cbind(Int,X1,X2,X3,X4,X5,X1sq,X2sq,X3sq,X4sq,X5sq,X1X2,X1X3,X1X4,X1X5,X2X3,X2X4,X2X5,X3X4,X3X5,X4X5)
}
}
# Case for 6 variables ########################
if(kvar1==6) {
names<-list("x1","x2","x3","x4","x5","x6")
names(des1)<-names
lmod<-stats::lm(y~x1+x2+x3+x4+x5+x6,data=des1)
imod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x2:x3+x2:x4+x2:x5+x2:x6+x3:x4+x3:x5+x3:x6+x4:x5+x4:x6+x5:x6,data=des1)
qmod<-stats::lm(y~x1+x2+x3+x4+x5+x6+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+I(x5^2)+I(x6^2)+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x2:x3+x2:x4+x2:x5+x2:x6+x3:x4+x3:x5+x3:x6+x4:x5+x4:x6+x5:x6,data=des1)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4,X5,X6)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X4X5<-X4*X5
X4X6<-X4*X6
X5X6<-X5*X6
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X1X2,X1X3,X1X4,X1X5,X1X6,X2X3,X2X4,X2X5,X2X6,X3X4,X3X5,X3X6,X4X5,X4X6,X5X6)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X5sq<-X5*X5
X6sq<-X6*X6
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X4X5<-X4*X5
X4X6<-X4*X6
X5X6<-X5*X6
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X1sq,X2sq,X3sq,X4sq,X5sq,X6sq,X1X2,X1X3,X1X4,X1X5,X1X6,X2X3,X2X4,X2X5,X2X6,X3X4,X3X5,X3X6,X4X5,X4X6,X5X6)
}
}
# Case for 7 variables ########################
if(kvar1==7) {
names<-list("x1","x2","x3","x4","x5","x6","x7")
names(des1)<-names
lmod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x7,data=des1)
imod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x7+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x1:x7+x2:x3+x2:x4+x2:x5+x2:x6+x2:x7+x3:x4+x3:x5+x3:x6+x3:x7+x4:x5+x4:x6+x4:x7+x5:x6+x5:x7+x6:x7,data=des1)
qmod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x7+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+I(x5^2)+I(x6^2)+I(x7^2)+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x1:x7+x2:x3+x2:x4+x2:x5+x2:x6+x2:x7+x3:x4+x3:x5+x3:x6+x3:x7+x4:x5+x4:x6+x4:x7+x5:x6+x5:x7+x6:x7,data=des1)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X7<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X7)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X7<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X1X7<-X1*X7
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X2X7<-X2*X7
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X3X7<-X3*X7
X4X5<-X4*X5
X4X6<-X4*X6
X4X7<-X4*X7
X5X6<-X5*X6
X5X7<-X5*X7
X6X7<-X6*X7
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X7,X1X2,X1X3,X1X4,X1X5,X1X6,X1X7,X2X3,X2X4,X2X5,X2X6,X2X7,X3X4,X3X5,X3X6,X3X7,X4X5,X4X6,X4X7,X5X6,X5X7,X6X7)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X7<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X5sq<-X5*X5
X6sq<-X6*X6
X7sq<-X7*X7
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X1X7<-X1*X7
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X2X7<-X2*X7
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X3X7<-X3*X7
X4X5<-X4*X5
X4X6<-X4*X6
X4X7<-X4*X7
X5X6<-X5*X6
X5X7<-X5*X7
X6X7<-X6*X7
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X7,X1sq,X2sq,X3sq,X4sq,X5sq,X6sq,X7sq,X1X2,X1X3,X1X4,X1X5,X1X6,X1X7,X2X3,X2X4,X2X5,X2X6,X2X7,X3X4,X3X5,X3X6,X3X7,X4X5,X4X6,X4X7,X5X6,X5X7,X6X7)
}
}
# get v1 and v1i
v1<-diag(fX%*%XtXI%*%t(fX))
v1i<-sort(v1)
# calculate fraction of design space
fs1<-(1:length(v1i))/length(v1i)
####################################################################
#
# Next check calculate the coordinates for des2
####################################################################
#First check to see if des2 is a matrix, and if so convert it to a data frame
if(is.matrix(des2)) des2 <- as.data.frame(des2)
#Next check to see if des2 is centered, and if not center it
chkm <-apply(des2, 2, f)
if (mean(chkm) !=0) des2<-scale(des2, scale=FALSE)
#Check to see if scaled, and if not scale it
chkm <-apply(des2, 2, mx)
if(mean(chkm) !=1) des2<-scale(des2, center=FALSE, scale=chkm)
rc <-dim(des2)
ndpts<-rc[1]
kvar1<-rc[2]
if (kvar1<2) {
stop("The minimum number of independent variables must be 2","\n")
}
if (kvar1>7) {
stop("The maximum number of independent variables allowed is 7","\n")
}
des2<-as.data.frame(des2)
# get the model
y<-stats::runif(nrow(des2),0,1)
### Case for two variables #############
if(kvar1==2) {
names<-list("x1","x2")
names(des2)<-names
lmod<-stats::lm(y~x1+x2,data=des2)
imod<-stats::lm(y~x1+x2+x1:x2,data=des2)
qmod<-stats::lm(y~x1+x2+I(x1^2)+I(x2^2)+x1:x2,data=des2)
if(mod==0){
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2)
}
if(mod==1){
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X1X2<-X1*X2
fX<-cbind(Int,X1,X2,X1X2)
}
if(mod==2){
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X1X2<-X1*X2
fX<-cbind(Int,X1,X2,X1sq,X2sq,X1X2)
}
}
# Case for 3 variables ########################
if(kvar1==3) {
names<-list("x1","x2","x3")
names(des2)<-names
lmod<-stats::lm(y~x1+x2+x3,data=des2)
imod<-stats::lm(y~x1+x2+x3+x1:x2+x1:x3+x2:x3,data=des2)
qmod<-stats::lm(y~x1+x2+x3+I(x1^2)+I(x2^2)+I(x3^2)+x1:x2+x1:x3+x2:x3,data=des2)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X2X3<-X2*X3
fX<-cbind(Int,X1,X2,X3,X1X2,X1X3,X2X3)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X1X2<-X1*X2
X1X3<-X1*X3
X2X3<-X2*X3
fX<-cbind(Int,X1,X2,X3,X1sq,X2sq,X3sq,X1X2,X1X3,X2X3)
}
}
# Case for 4 variables ########################
if(kvar1==4) {
names<-list("x1","x2","x3","x4")
names(des2)<-names
lmod<-stats::lm(y~x1+x2+x3+x4,data=des2)
imod<-stats::lm(y~x1+x2+x3+x4+x1:x2+x1:x3+x1:x4+x2:x3+x2:x4+x3:x4,data=des2)
qmod<-stats::lm(y~x1+x2+x3+x4+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+x1:x2+x1:x3+x1:x4+x2:x3+x2:x4+x3:x4,data=des2)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X2X3<-X2*X3
X2X4<-X2*X4
X3X4<-X3*X4
fX<-cbind(Int,X1,X2,X3,X4,X1X2,X1X3,X1X4,X2X3,X2X4,X3X4)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X2X3<-X2*X3
X2X4<-X2*X4
X3X4<-X3*X4
fX<-cbind(Int,X1,X2,X3,X4,X1sq,X2sq,X3sq,X4sq,X1X2,X1X3,X1X4,X2X3,X2X4,X3X4)
}
}
# Case for 5 variables ########################
if(kvar1==5) {
names<-list("x1","x2","x3","x4","x5")
names(des2)<-names
lmod<-stats::lm(y~x1+x2+x3+x4+x5,data=des2)
imod<-stats::lm(y~x1+x2+x3+x4+x5+x1:x2+x1:x3+x1:x4+x1:x5+x2:x3+x2:x4+x2:x5+x3:x4+x3:x5+x4:x5,data=des2)
qmod<-stats::lm(y~x1+x2+x3+x4+x5+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+I(x5^2)+x1:x2+x1:x3+x1:x4+x1:x5+x2:x3+x2:x4+x2:x5+x3:x4+x3:x5+x4:x5,data=des2)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4,X5)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X3X4<-X3*X4
X3X5<-X3*X5
X4X5<-X4*X5
fX<-cbind(Int,X1,X2,X3,X4,+X5+X1X2,X1X3,X1X4,X1X5,X2X3,X2X4,X2X5,X3X4,X3X5,X4X5)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X5sq<-X5*X5
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X3X4<-X3*X4
X3X5<-X3*X5
X4X5<-X4*X5
fX<-cbind(Int,X1,X2,X3,X4,X5,X1sq,X2sq,X3sq,X4sq,X5sq,X1X2,X1X3,X1X4,X1X5,X2X3,X2X4,X2X5,X3X4,X3X5,X4X5)
}
}
# Case for 6 variables ########################
if(kvar1==6) {
names<-list("x1","x2","x3","x4","x5","x6")
names(des2)<-names
lmod<-stats::lm(y~x1+x2+x3+x4+x5+x6,data=des2)
imod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x2:x3+x2:x4+x2:x5+x2:x6+x3:x4+x3:x5+x3:x6+x4:x5+x4:x6+x5:x6,data=des2)
qmod<-stats::lm(y~x1+x2+x3+x4+x5+x6+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+I(x5^2)+I(x6^2)+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x2:x3+x2:x4+x2:x5+x2:x6+x3:x4+x3:x5+x3:x6+x4:x5+x4:x6+x5:x6,data=des2)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4,X5,X6)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X4X5<-X4*X5
X4X6<-X4*X6
X5X6<-X5*X6
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X1X2,X1X3,X1X4,X1X5,X1X6,X2X3,X2X4,X2X5,X2X6,X3X4,X3X5,X3X6,X4X5,X4X6,X5X6)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X5sq<-X5*X5
X6sq<-X6*X6
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X4X5<-X4*X5
X4X6<-X4*X6
X5X6<-X5*X6
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X1sq,X2sq,X3sq,X4sq,X5sq,X6sq,X1X2,X1X3,X1X4,X1X5,X1X6,X2X3,X2X4,X2X5,X2X6,X3X4,X3X5,X3X6,X4X5,X4X6,X5X6)
}
}
# Case for 7 variables ########################
if(kvar1==7) {
names<-list("x1","x2","x3","x4","x5","x6","x7")
names(des2)<-names
lmod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x7,data=des2)
imod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x7+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x1:x7+x2:x3+x2:x4+x2:x5+x2:x6+x2:x7+x3:x4+x3:x5+x3:x6+x3:x7+x4:x5+x4:x6+x4:x7+x5:x6+x5:x7+x6:x7,data=des2)
qmod<-stats::lm(y~x1+x2+x3+x4+x5+x6+x7+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+I(x5^2)+I(x6^2)+I(x7^2)+x1:x2+x1:x3+x1:x4+x1:x5+x1:x6+x1:x7+x2:x3+x2:x4+x2:x5+x2:x6+x2:x7+x3:x4+x3:x5+x3:x6+x3:x7+x4:x5+x4:x6+x4:x7+x5:x6+x5:x7+x6:x7,data=des2)
#get the design matrix and fX
if(mod==0) {
X<-stats::model.matrix(lmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X7<-stats::runif(5000,-1,1)
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X7)
}
if(mod==1) {
X<-stats::model.matrix(imod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X7<-stats::runif(5000,-1,1)
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X1X7<-X1*X7
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X2X7<-X2*X7
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X3X7<-X3*X7
X4X5<-X4*X5
X4X6<-X4*X6
X4X7<-X4*X7
X5X6<-X5*X6
X5X7<-X5*X7
X6X7<-X6*X7
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X7,X1X2,X1X3,X1X4,X1X5,X1X6,X1X7,X2X3,X2X4,X2X5,X2X6,X2X7,X3X4,X3X5,X3X6,X3X7,X4X5,X4X6,X4X7,X5X6,X5X7,X6X7)
}
if(mod==2) {
X<-stats::model.matrix(qmod)
XtXI<-solve(t(X)%*%X)
Int<-c(rep(1,5000))
X1<-stats::runif(5000,-1,1)
X2<-stats::runif(5000,-1,1)
X3<-stats::runif(5000,-1,1)
X4<-stats::runif(5000,-1,1)
X5<-stats::runif(5000,-1,1)
X6<-stats::runif(5000,-1,1)
X7<-stats::runif(5000,-1,1)
X1sq<-X1*X1
X2sq<-X2*X2
X3sq<-X3*X3
X4sq<-X4*X4
X5sq<-X5*X5
X6sq<-X6*X6
X7sq<-X7*X7
X1X2<-X1*X2
X1X3<-X1*X3
X1X4<-X1*X4
X1X5<-X1*X5
X1X6<-X1*X6
X1X7<-X1*X7
X2X3<-X2*X3
X2X4<-X2*X4
X2X5<-X2*X5
X2X6<-X2*X6
X2X7<-X2*X7
X3X4<-X3*X4
X3X5<-X3*X5
X3X6<-X3*X6
X3X7<-X3*X7
X4X5<-X4*X5
X4X6<-X4*X6
X4X7<-X4*X7
X5X6<-X5*X6
X5X7<-X5*X7
X6X7<-X6*X7
fX<-cbind(Int,X1,X2,X3,X4,X5,X6,X7,X1sq,X2sq,X3sq,X4sq,X5sq,X6sq,X7sq,X1X2,X1X3,X1X4,X1X5,X1X6,X1X7,X2X3,X2X4,X2X5,X2X6,X2X7,X3X4,X3X5,X3X6,X3X7,X4X5,X4X6,X4X7,X5X6,X5X7,X6X7)
}
}
# get v2 and v2i
v2<-diag(fX%*%XtXI%*%t(fX))
v2i<-sort(v2)
# calculate fraction of design space
fs2<-(1:length(v2i))/length(v2i)
# make Comparison FDS Plots
maxv<-max(max(v1i),max(v2i))
plot(fs1,v1i,type='l',ylim=c(0.0,(maxv+.15)),xlab="Fraction of Space",ylab="Relative Prediction Variance")
graphics::abline(v=.5,lty=3)
graphics::abline(h=(v1i[2500]+v1i[2501])/2,lty=3)
graphics::lines(fs2,v2i,col="red",lty=2)
graphics::abline(h=(v2i[2500]+v2i[2501])/2,lty=3)
graphics::lines(c(0,.10),c(maxv,maxv),lty=1,col="black")
graphics::text(.11,maxv,name1,adj=-.1)
graphics::lines(c(0,.10),c(maxv-maxv/20,maxv-maxv/20),lty=3,col="red")
graphics::text(.11,maxv-maxv/20,name2,col="red",adj=-.1)
}
Try the Vdgraph package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.