The design that we utilized in Section 4.2 to illustrate the integrative framework for geometric and hidden projection is the DSD(\(11\)) in Table 4 of the supplement. In this case study, we supposed it was of interest to augment the projection on factors \(A_1, A_2, A_3\), and \(A_4\) of the DSD(\(11\)) such that the augmented design satisfies the second-order orthogonality criterion. The code below specifies this design matrix.
# The design matrix that we consider follows below in F_design
F_design = matrix(c(0,0,1,-1,1,-1,1,-1,1,-1,0,
1,-1,0,0,-1,1,-1,1,1,-1,0,
1,-1,-1,1,0,0,1,-1,1,-1,0,
-1,1,-1,1,1,-1,0,0,1,-1,0), nrow=11, ncol=4)Our formulation of this task based on Corollary \(1\) requires the calculation of all indicator function coefficients that correspond to main effects, two-factor interactions, and linear-linear-linear three-factor interactions among \(A_1, A_2, A_3\), and \(A_4\) in the DSD(\(11\)). These are especially simple to perform with the \(\otimes\) operation for DSDs due to the fold-over structures of these designs. The code for calculating all of the indicator function coefficients involving just the factors in the projection of the DSD(\(11\)) follow below.
# First, we calculate indicator function coefficients of order 1
indicator_function_coefficients_order_1 = function(design_matrix, factor_index)
{
hadamard_matrix = matrix(c(1,1,-1,1),nrow=2,ncol=2)
indicator_function_coefficients = 2^(-1)*3^(1-ncol(design_matrix))*hadamard_matrix%*%
c(sum(design_matrix[,factor_index]==1),
sum(design_matrix[,factor_index]==-1)) -
c(0, nrow(design_matrix)/3^ncol(design_matrix))
return(indicator_function_coefficients)
}
indicator_function_coefficients_1 = indicator_function_coefficients_order_1(F_design, 1)
indicator_function_coefficients_2 = indicator_function_coefficients_order_1(F_design, 2)
indicator_function_coefficients_3 = indicator_function_coefficients_order_1(F_design, 3)
indicator_function_coefficients_4 = indicator_function_coefficients_order_1(F_design, 4)
# Second, we calculate indicator function coefficients of order 2
indicator_function_coefficients_order_2 = function(design_matrix, factor_indices)
{
# We assume that factor_indices is a vector of size 2 that is ordered lexicographically
indicator_function_coefficients_factor_1 = indicator_function_coefficients_order_1(design_matrix, factor_indices[1])
indicator_function_coefficients_factor_2 = indicator_function_coefficients_order_1(design_matrix, factor_indices[2])
hadamard_matrix = matrix(c(1,1,-1,1),nrow=2,ncol=2)
hadamard_matrix = kronecker(hadamard_matrix,hadamard_matrix)
indicator_function_coefficients = 2^(-2)*3^(2-ncol(design_matrix))*hadamard_matrix%*%
c(sum(design_matrix[,factor_indices[1]]==1 & design_matrix[,factor_indices[2]]==1),
sum(design_matrix[,factor_indices[1]]==1 & design_matrix[,factor_indices[2]]==-1),
sum(design_matrix[,factor_indices[1]]==-1 & design_matrix[,factor_indices[2]]==1),
sum(design_matrix[,factor_indices[1]]==-1 & design_matrix[,factor_indices[2]]==-1)) -
c(0,
indicator_function_coefficients_factor_1[1],
indicator_function_coefficients_factor_2[1],
indicator_function_coefficients_factor_1[2] +
indicator_function_coefficients_factor_2[2] + nrow(design_matrix)/3^ncol(design_matrix))
return(indicator_function_coefficients)
}
indicator_function_coefficients_12 = indicator_function_coefficients_order_2(F_design, c(1,2))
indicator_function_coefficients_13 = indicator_function_coefficients_order_2(F_design, c(1,3))
indicator_function_coefficients_14 = indicator_function_coefficients_order_2(F_design, c(1,4))
indicator_function_coefficients_23 = indicator_function_coefficients_order_2(F_design, c(2,3))
indicator_function_coefficients_24 = indicator_function_coefficients_order_2(F_design, c(2,4))
indicator_function_coefficients_34 = indicator_function_coefficients_order_2(F_design, c(3,4))
# Third, we calculate indicator function coefficients of order 3
indicator_function_coefficients_order_3 = function(design_matrix, factor_indices)
{
# We assume that factor_indices is a vector of size 3 that is ordered lexicographically
indicator_function_coefficients_factor_1 = indicator_function_coefficients_order_1(design_matrix, factor_indices[1])
indicator_function_coefficients_factor_2 = indicator_function_coefficients_order_1(design_matrix, factor_indices[2])
indicator_function_coefficients_factor_3 = indicator_function_coefficients_order_1(design_matrix, factor_indices[3])
indicator_function_coefficients_factor_12 = indicator_function_coefficients_order_2(design_matrix, factor_indices[c(1,2)])
indicator_function_coefficients_factor_13 = indicator_function_coefficients_order_2(design_matrix, factor_indices[c(1,3)])
indicator_function_coefficients_factor_23 = indicator_function_coefficients_order_2(design_matrix, factor_indices[c(2,3)])
hadamard_matrix = matrix(c(1,1,-1,1),nrow=2,ncol=2)
hadamard_matrix = kronecker(hadamard_matrix,hadamard_matrix)
hadamard_matrix = kronecker(matrix(c(1,1,-1,1),nrow=2,ncol=2),hadamard_matrix)
indicator_function_coefficients = 2^(-3)*3^(3-ncol(design_matrix))*hadamard_matrix%*%
c(sum(design_matrix[,factor_indices[1]]==1 & design_matrix[,factor_indices[2]]==1 & design_matrix[,factor_indices[3]]==1),
sum(design_matrix[,factor_indices[1]]==1 & design_matrix[,factor_indices[2]]==1 & design_matrix[,factor_indices[3]]==-1),
sum(design_matrix[,factor_indices[1]]==1 & design_matrix[,factor_indices[2]]==-1 & design_matrix[,factor_indices[3]]==1),
sum(design_matrix[,factor_indices[1]]==1 & design_matrix[,factor_indices[2]]==-1 & design_matrix[,factor_indices[3]]==-1),
sum(design_matrix[,factor_indices[1]]==-1 & design_matrix[,factor_indices[2]]==1 & design_matrix[,factor_indices[3]]==1),
sum(design_matrix[,factor_indices[1]]==-1 & design_matrix[,factor_indices[2]]==1 & design_matrix[,factor_indices[3]]==-1),
sum(design_matrix[,factor_indices[1]]==-1 & design_matrix[,factor_indices[2]]==-1 & design_matrix[,factor_indices[3]]==1),
sum(design_matrix[,factor_indices[1]]==-1 & design_matrix[,factor_indices[2]]==-1 & design_matrix[,factor_indices[3]]==-1)) -
c(0,
indicator_function_coefficients_factor_12[1],
indicator_function_coefficients_factor_13[1],
indicator_function_coefficients_factor_1[1] +
indicator_function_coefficients_factor_12[2] +
indicator_function_coefficients_factor_13[2],
indicator_function_coefficients_factor_23[1],
indicator_function_coefficients_factor_2[1] +
indicator_function_coefficients_factor_12[3] +
indicator_function_coefficients_factor_23[2],
indicator_function_coefficients_factor_3[1] +
indicator_function_coefficients_factor_13[3] +
indicator_function_coefficients_factor_23[3],
indicator_function_coefficients_factor_1[2] +
indicator_function_coefficients_factor_2[2] +
indicator_function_coefficients_factor_3[2] +
indicator_function_coefficients_factor_12[4] +
indicator_function_coefficients_factor_13[4] +
indicator_function_coefficients_factor_23[4] + nrow(design_matrix)/3^ncol(design_matrix))
return(indicator_function_coefficients)
}
indicator_function_coefficients_123 = indicator_function_coefficients_order_3(F_design, c(1,2,3))
indicator_function_coefficients_124 = indicator_function_coefficients_order_3(F_design, c(1,2,4))
indicator_function_coefficients_134 = indicator_function_coefficients_order_3(F_design, c(1,3,4))
indicator_function_coefficients_234 = indicator_function_coefficients_order_3(F_design, c(2,3,4))
# Fourth, we calculate indicator function coefficients of order 4
indicator_function_coefficients_order_4 = function(design_matrix, factor_indices)
{
# We assume that factor_indices is a vector of size 4 that is ordered lexicographically
indicator_function_coefficients_factor_1 = indicator_function_coefficients_order_1(design_matrix, factor_indices[1])
indicator_function_coefficients_factor_2 = indicator_function_coefficients_order_1(design_matrix, factor_indices[2])
indicator_function_coefficients_factor_3 = indicator_function_coefficients_order_1(design_matrix, factor_indices[3])
indicator_function_coefficients_factor_4 = indicator_function_coefficients_order_1(design_matrix, factor_indices[4])
indicator_function_coefficients_factor_12 = indicator_function_coefficients_order_2(design_matrix, factor_indices[c(1,2)])
indicator_function_coefficients_factor_13 = indicator_function_coefficients_order_2(design_matrix, factor_indices[c(1,3)])
indicator_function_coefficients_factor_14 = indicator_function_coefficients_order_2(design_matrix, factor_indices[c(1,4)])
indicator_function_coefficients_factor_23 = indicator_function_coefficients_order_2(design_matrix, factor_indices[c(2,3)])
indicator_function_coefficients_factor_24 = indicator_function_coefficients_order_2(design_matrix, factor_indices[c(2,4)])
indicator_function_coefficients_factor_34 = indicator_function_coefficients_order_2(design_matrix, factor_indices[c(3,4)])
indicator_function_coefficients_factor_123 = indicator_function_coefficients_order_3(design_matrix, factor_indices[c(1,2,3)])
indicator_function_coefficients_factor_124 = indicator_function_coefficients_order_3(design_matrix, factor_indices[c(1,2,4)])
indicator_function_coefficients_factor_134 = indicator_function_coefficients_order_3(design_matrix, factor_indices[c(1,3,4)])
indicator_function_coefficients_factor_234 = indicator_function_coefficients_order_3(design_matrix, factor_indices[c(2,3,4)])
hadamard_matrix = matrix(c(1,1,-1,1),nrow=2,ncol=2)
hadamard_matrix = kronecker(hadamard_matrix,hadamard_matrix)
hadamard_matrix = kronecker(matrix(c(1,1,-1,1),nrow=2,ncol=2),hadamard_matrix)
hadamard_matrix = kronecker(matrix(c(1,1,-1,1),nrow=2,ncol=2),hadamard_matrix)
indicator_function_coefficients = 2^(-4)*3^(4-ncol(design_matrix))*hadamard_matrix%*%
c(sum(design_matrix[,factor_indices[1]]==1 & design_matrix[,factor_indices[2]]==1 & design_matrix[,factor_indices[3]]==1 & design_matrix[,factor_indices[4]]==1),
sum(design_matrix[,factor_indices[1]]==1 & design_matrix[,factor_indices[2]]==1 & design_matrix[,factor_indices[3]]==1 & design_matrix[,factor_indices[4]]==-1),
sum(design_matrix[,factor_indices[1]]==1 & design_matrix[,factor_indices[2]]==1 & design_matrix[,factor_indices[3]]==-1 & design_matrix[,factor_indices[4]]==1),
sum(design_matrix[,factor_indices[1]]==1 & design_matrix[,factor_indices[2]]==1 & design_matrix[,factor_indices[3]]==-1 & design_matrix[,factor_indices[4]]==-1),
sum(design_matrix[,factor_indices[1]]==1 & design_matrix[,factor_indices[2]]==-1 & design_matrix[,factor_indices[3]]==1 & design_matrix[,factor_indices[4]]==1),
sum(design_matrix[,factor_indices[1]]==1 & design_matrix[,factor_indices[2]]==-1 & design_matrix[,factor_indices[3]]==1 & design_matrix[,factor_indices[4]]==-1),
sum(design_matrix[,factor_indices[1]]==1 & design_matrix[,factor_indices[2]]==-1 & design_matrix[,factor_indices[3]]==-1 & design_matrix[,factor_indices[4]]==1),
sum(design_matrix[,factor_indices[1]]==1 & design_matrix[,factor_indices[2]]==-1 & design_matrix[,factor_indices[3]]==-1 & design_matrix[,factor_indices[4]]==-1),
sum(design_matrix[,factor_indices[1]]==-1 & design_matrix[,factor_indices[2]]==1 & design_matrix[,factor_indices[3]]==1 & design_matrix[,factor_indices[4]]==1),
sum(design_matrix[,factor_indices[1]]==-1 & design_matrix[,factor_indices[2]]==1 & design_matrix[,factor_indices[3]]==1 & design_matrix[,factor_indices[4]]==-1),
sum(design_matrix[,factor_indices[1]]==-1 & design_matrix[,factor_indices[2]]==1 & design_matrix[,factor_indices[3]]==-1 & design_matrix[,factor_indices[4]]==1),
sum(design_matrix[,factor_indices[1]]==-1 & design_matrix[,factor_indices[2]]==1 & design_matrix[,factor_indices[3]]==-1 & design_matrix[,factor_indices[4]]==-1),
sum(design_matrix[,factor_indices[1]]==-1 & design_matrix[,factor_indices[2]]==-1 & design_matrix[,factor_indices[3]]==1 & design_matrix[,factor_indices[4]]==1),
sum(design_matrix[,factor_indices[1]]==-1 & design_matrix[,factor_indices[2]]==-1 & design_matrix[,factor_indices[3]]==1 & design_matrix[,factor_indices[4]]==-1),
sum(design_matrix[,factor_indices[1]]==-1 & design_matrix[,factor_indices[2]]==-1 & design_matrix[,factor_indices[3]]==-1 & design_matrix[,factor_indices[4]]==1),
sum(design_matrix[,factor_indices[1]]==-1 & design_matrix[,factor_indices[2]]==-1 & design_matrix[,factor_indices[3]]==-1 & design_matrix[,factor_indices[4]]==-1)) -
c(0,
indicator_function_coefficients_factor_123[1],
indicator_function_coefficients_factor_124[1],
indicator_function_coefficients_factor_12[1] +
indicator_function_coefficients_factor_123[2] +
indicator_function_coefficients_factor_124[2],
indicator_function_coefficients_factor_134[1],
indicator_function_coefficients_factor_13[1] +
indicator_function_coefficients_factor_123[3] +
indicator_function_coefficients_factor_134[2],
indicator_function_coefficients_factor_14[1] +
indicator_function_coefficients_factor_124[3] +
indicator_function_coefficients_factor_134[3],
indicator_function_coefficients_factor_1[1] +
indicator_function_coefficients_factor_12[2] +
indicator_function_coefficients_factor_13[2] +
indicator_function_coefficients_factor_14[2] +
indicator_function_coefficients_factor_123[4] +
indicator_function_coefficients_factor_124[4] +
indicator_function_coefficients_factor_134[4],
indicator_function_coefficients_factor_234[1],
indicator_function_coefficients_factor_23[1] +
indicator_function_coefficients_factor_123[5] +
indicator_function_coefficients_factor_234[2],
indicator_function_coefficients_factor_24[1] +
indicator_function_coefficients_factor_124[5] +
indicator_function_coefficients_factor_234[3],
indicator_function_coefficients_factor_2[1] +
indicator_function_coefficients_factor_12[3] +
indicator_function_coefficients_factor_23[2] +
indicator_function_coefficients_factor_24[2] +
indicator_function_coefficients_factor_123[6] +
indicator_function_coefficients_factor_124[6] +
indicator_function_coefficients_factor_234[4],
indicator_function_coefficients_factor_34[1] +
indicator_function_coefficients_factor_134[5] +
indicator_function_coefficients_factor_234[5],
indicator_function_coefficients_factor_3[1] +
indicator_function_coefficients_factor_13[3] +
indicator_function_coefficients_factor_23[3] +
indicator_function_coefficients_factor_34[2] +
indicator_function_coefficients_factor_123[7] +
indicator_function_coefficients_factor_134[6] +
indicator_function_coefficients_factor_234[6],
indicator_function_coefficients_factor_4[1] +
indicator_function_coefficients_factor_14[3] +
indicator_function_coefficients_factor_24[3] +
indicator_function_coefficients_factor_34[3] +
indicator_function_coefficients_factor_124[7] +
indicator_function_coefficients_factor_134[7] +
indicator_function_coefficients_factor_234[7],
nrow(design_matrix)/3^ncol(design_matrix) +
indicator_function_coefficients_factor_1[2] +
indicator_function_coefficients_factor_2[2] +
indicator_function_coefficients_factor_3[2] +
indicator_function_coefficients_factor_4[2] +
indicator_function_coefficients_factor_12[4] +
indicator_function_coefficients_factor_13[4] +
indicator_function_coefficients_factor_14[4] +
indicator_function_coefficients_factor_23[4] +
indicator_function_coefficients_factor_24[4] +
indicator_function_coefficients_factor_34[4] +
indicator_function_coefficients_factor_123[8] +
indicator_function_coefficients_factor_124[8] +
indicator_function_coefficients_factor_134[8] +
indicator_function_coefficients_factor_234[8])
return(indicator_function_coefficients)
}
indicator_function_coefficients_1234 = indicator_function_coefficients_order_4(F_design, c(1,2,3,4))
# Finally, we organize all of the previously calculated indicator function coefficients into one vector
b_F = c(nrow(F_design)/3^ncol(F_design),
indicator_function_coefficients_1,
indicator_function_coefficients_2,
indicator_function_coefficients_3,
indicator_function_coefficients_4,
indicator_function_coefficients_12,
indicator_function_coefficients_13,
indicator_function_coefficients_14,
indicator_function_coefficients_23,
indicator_function_coefficients_24,
indicator_function_coefficients_34,
indicator_function_coefficients_123,
indicator_function_coefficients_124,
indicator_function_coefficients_134,
indicator_function_coefficients_234,
indicator_function_coefficients_1234)The code below contains our formulation of the optimization for deriving an augmentation for the projection of the DSD(\(11\)) onto its first four factors.
# We first specify the full model matrix in Definition 6 that corresponds to the number of columns in F_design
k = ncol(F_design)
x_1_L = c(rep(-1,3^(k-1)),rep(0,3^(k-1)),rep(1,3^(k-1)))
x_1_Q = 3*x_1_L^2 - 2
x_2_L = rep(c(rep(-1,3^(k-2)),rep(0,3^(k-2)),rep(1,3^(k-2))),3^1)
x_2_Q = 3*x_2_L^2 - 2
x_3_L = rep(c(rep(-1,3^(k-3)),rep(0,3^(k-3)),rep(1,3^(k-3))),3^2)
x_3_Q = 3*x_3_L^2 - 2
x_4_L = rep(c(rep(-1,3^(k-4)),rep(0,3^(k-4)),rep(1,3^(k-4))),3^3)
x_4_Q = 3*x_4_L^2 - 2
x_12_LL = x_1_L*x_2_L
x_12_LQ = x_1_L*x_2_Q
x_12_QL = x_1_Q*x_2_L
x_12_QQ = x_1_Q*x_2_Q
x_13_LL = x_1_L*x_3_L
x_13_LQ = x_1_L*x_3_Q
x_13_QL = x_1_Q*x_3_L
x_13_QQ = x_1_Q*x_3_Q
x_14_LL = x_1_L*x_4_L
x_14_LQ = x_1_L*x_4_Q
x_14_QL = x_1_Q*x_4_L
x_14_QQ = x_1_Q*x_4_Q
x_23_LL = x_2_L*x_3_L
x_23_LQ = x_2_L*x_3_Q
x_23_QL = x_2_Q*x_3_L
x_23_QQ = x_2_Q*x_3_Q
x_24_LL = x_2_L*x_4_L
x_24_LQ = x_2_L*x_4_Q
x_24_QL = x_2_Q*x_4_L
x_24_QQ = x_2_Q*x_4_Q
x_34_LL = x_3_L*x_4_L
x_34_LQ = x_3_L*x_4_Q
x_34_QL = x_3_Q*x_4_L
x_34_QQ = x_3_Q*x_4_Q
x_123_LLL = x_1_L*x_2_L*x_3_L
x_123_LLQ = x_1_L*x_2_L*x_3_Q
x_123_LQL = x_1_L*x_2_Q*x_3_L
x_123_LQQ = x_1_L*x_2_Q*x_3_Q
x_123_QLL = x_1_Q*x_2_L*x_3_L
x_123_QLQ = x_1_Q*x_2_L*x_3_Q
x_123_QQL = x_1_Q*x_2_Q*x_3_L
x_123_QQQ = x_1_Q*x_2_Q*x_3_Q
x_124_LLL = x_1_L*x_2_L*x_4_L
x_124_LLQ = x_1_L*x_2_L*x_4_Q
x_124_LQL = x_1_L*x_2_Q*x_4_L
x_124_LQQ = x_1_L*x_2_Q*x_4_Q
x_124_QLL = x_1_Q*x_2_L*x_4_L
x_124_QLQ = x_1_Q*x_2_L*x_4_Q
x_124_QQL = x_1_Q*x_2_Q*x_4_L
x_124_QQQ = x_1_Q*x_2_Q*x_4_Q
x_134_LLL = x_1_L*x_3_L*x_4_L
x_134_LLQ = x_1_L*x_3_L*x_4_Q
x_134_LQL = x_1_L*x_3_Q*x_4_L
x_134_LQQ = x_1_L*x_3_Q*x_4_Q
x_134_QLL = x_1_Q*x_3_L*x_4_L
x_134_QLQ = x_1_Q*x_3_L*x_4_Q
x_134_QQL = x_1_Q*x_3_Q*x_4_L
x_134_QQQ = x_1_Q*x_3_Q*x_4_Q
x_234_LLL = x_2_L*x_3_L*x_4_L
x_234_LLQ = x_2_L*x_3_L*x_4_Q
x_234_LQL = x_2_L*x_3_Q*x_4_L
x_234_LQQ = x_2_L*x_3_Q*x_4_Q
x_234_QLL = x_2_Q*x_3_L*x_4_L
x_234_QLQ = x_2_Q*x_3_L*x_4_Q
x_234_QQL = x_2_Q*x_3_Q*x_4_L
x_234_QQQ = x_2_Q*x_3_Q*x_4_Q
x_1234_LLLL = x_1_L*x_2_L*x_3_L*x_4_L
x_1234_LLLQ = x_1_L*x_2_L*x_3_L*x_4_Q
x_1234_LLQL = x_1_L*x_2_L*x_3_Q*x_4_L
x_1234_LLQQ = x_1_L*x_2_L*x_3_Q*x_4_Q
x_1234_LQLL = x_1_L*x_2_Q*x_3_L*x_4_L
x_1234_LQLQ = x_1_L*x_2_Q*x_3_L*x_4_Q
x_1234_LQQL = x_1_L*x_2_Q*x_3_Q*x_4_L
x_1234_LQQQ = x_1_L*x_2_Q*x_3_Q*x_4_Q
x_1234_QLLL = x_1_Q*x_2_L*x_3_L*x_4_L
x_1234_QLLQ = x_1_Q*x_2_L*x_3_L*x_4_Q
x_1234_QLQL = x_1_Q*x_2_L*x_3_Q*x_4_L
x_1234_QLQQ = x_1_Q*x_2_L*x_3_Q*x_4_Q
x_1234_QQLL = x_1_Q*x_2_Q*x_3_L*x_4_L
x_1234_QQLQ = x_1_Q*x_2_Q*x_3_L*x_4_Q
x_1234_QQQL = x_1_Q*x_2_Q*x_3_Q*x_4_L
x_1234_QQQQ = x_1_Q*x_2_Q*x_3_Q*x_4_Q
full_model_matrix = cbind(rep(1, 3^k),
x_1_L,
x_1_Q,
x_2_L,
x_2_Q,
x_3_L,
x_3_Q,
x_4_L,
x_4_Q,
x_12_LL,
x_12_LQ,
x_12_QL,
x_12_QQ,
x_13_LL,
x_13_LQ,
x_13_QL,
x_13_QQ,
x_14_LL,
x_14_LQ,
x_14_QL,
x_14_QQ,
x_23_LL,
x_23_LQ,
x_23_QL,
x_23_QQ,
x_24_LL,
x_24_LQ,
x_24_QL,
x_24_QQ,
x_34_LL,
x_34_LQ,
x_34_QL,
x_34_QQ,
x_123_LLL,
x_123_LLQ,
x_123_LQL,
x_123_LQQ,
x_123_QLL,
x_123_QLQ,
x_123_QQL,
x_123_QQQ,
x_124_LLL,
x_124_LLQ,
x_124_LQL,
x_124_LQQ,
x_124_QLL,
x_124_QLQ,
x_124_QQL,
x_124_QQQ,
x_134_LLL,
x_134_LLQ,
x_134_LQL,
x_134_LQQ,
x_134_QLL,
x_134_QLQ,
x_134_QQL,
x_134_QQQ,
x_234_LLL,
x_234_LLQ,
x_234_LQL,
x_234_LQQ,
x_234_QLL,
x_234_QLQ,
x_234_QQL,
x_234_QQQ,
x_1234_LLLL,
x_1234_LLLQ,
x_1234_LLQL,
x_1234_LLQQ,
x_1234_LQLL,
x_1234_LQLQ,
x_1234_LQQL,
x_1234_LQQQ,
x_1234_QLLL,
x_1234_QLLQ,
x_1234_QLQL,
x_1234_QLQQ,
x_1234_QQLL,
x_1234_QQLQ,
x_1234_QQQL,
x_1234_QQQQ)
colnames(full_model_matrix) = NULL
# We now take into account which columns of the full model matrix correspond to the desired aliasing relations in the augmented design
full_model_matrix_specified = full_model_matrix[,c(2:34,42,50,58)]
full_model_matrix_unspecified = full_model_matrix[,-c(2:34,42,50,58)]
# We also take into account the specified indicator function coefficients for the augmentation that correspond to the desired aliasing relations in the augmented design
b_specified = matrix(-b_F[c(2:34,42,50,58)], ncol=1)
# We finally specify the objective function for the optimization according to the approach corresponding to equation (3)
# The fixed value of b_{\mathcal{G}, \phi, \phi} that we consider in this optimization approach is contained in b_phi_phi below
b_phi_phi = 16/3^4
objective_function = function(z)
{
x = matrix(c(b_phi_phi,z), ncol=1)
objective_value = (full_model_matrix_unspecified%*%x+
full_model_matrix_specified%*%b_specified)*
(matrix(rep(1,nrow(full_model_matrix_unspecified)), ncol=1)-
full_model_matrix_unspecified%*%x-
full_model_matrix_specified%*%b_specified)
return(sum((objective_value)^2))
}The code below contains our execution of the optimization for deriving an augmentation for the projection of the DSD(\(11\)) onto its first four factors. Note that we consider a single starting point for this particular optimization. Additional optimum augmentations can be found by executing this code for multiple starting points.
# We perform the optimization using the BFGS method in the optim function
# The starting point for this optimization is contained in starting_parameters
starting_parameters = c(-5e-02,-1e-01,0,5e-02,
0,0,2e-02,-5e-02,
5e-02,0,1e-01,0,
0,2e-02,-1e-01,5e-02,
0,5e-02,0,0,
2e-02,5e-02,1e-01,0,
5e-02,0,0,2e-02,
-2e-01,0,0,5e-02,
0,2e-02,-5e-02,0,
0,-5e-02,-2e-02,0,
-5e-02,0,0,4e-02)
optim_output = optim(starting_parameters, objective_function,
method="BFGS",
control=list(maxit=10000))
solution = optim_output$parFinally, after executing the optimization, we create the follow-up design from the resulting indicator function coefficients.
# First, we obtain the indicator function coefficients for the follow-up design
new_indicator_function_coefficient_vector = c(b_phi_phi,
b_specified[1:33],
solution[1:7],
b_specified[34],
solution[8:14],
b_specified[35],
solution[15:21],
b_specified[36],
solution[22:length(solution)])
# Second, we obtain the actual runs of the follow-up design
run_indicators = round(full_model_matrix%*%new_indicator_function_coefficient_vector, 0.00000001)
follow_up_design = full_model_matrix[run_indicators==1,c(2,4,6,8)]
# Finally, we construct the augmented design
augmented_design = rbind(F_design, follow_up_design)
augmented_design## [,1] [,2] [,3] [,4]
## [1,] 0 1 1 -1
## [2,] 0 -1 -1 1
## [3,] 1 0 -1 -1
## [4,] -1 0 1 1
## [5,] 1 -1 0 1
## [6,] -1 1 0 -1
## [7,] 1 -1 1 0
## [8,] -1 1 -1 0
## [9,] 1 1 1 1
## [10,] -1 -1 -1 -1
## [11,] 0 0 0 0
## [12,] -1 -1 0 0
## [13,] -1 -1 1 -1
## [14,] -1 0 -1 0
## [15,] -1 0 0 1
## [16,] -1 1 1 1
## [17,] 0 -1 0 1
## [18,] 0 -1 1 0
## [19,] 0 0 -1 1
## [20,] 0 0 1 -1
## [21,] 0 1 -1 0
## [22,] 0 1 0 -1
## [23,] 1 -1 -1 -1
## [24,] 1 0 0 -1
## [25,] 1 0 1 0
## [26,] 1 1 -1 1
## [27,] 1 1 0 0