# D-Optimal Design for a Causal Structure for Completely Randomized and Random Blocked Experiments

Authors: Zaher Kmail; Kent Eskridge
Journal: Journal of Probability and Statistics
 (2022)

---

## Abstract

Most experimental design literature on causal inference focuses on establishing a causal relationship between variables, but there is no literature on how to identify a design that results in the optimal parameter estimates for a structural equation model (SEM). In this research, search algorithms are used to produce aD-optimal design for a SEM for three-stage least squares and full information maximum likelihood estimators. Then, a D-optimal design for the estimate of the model parameters of a mixed-effects SEM is obtained. The efficiency of each of the D-optimal designs for SEMs is compared with univariate optimal and uniform designs. In each case, the causal relationship changed the optimal designs dramatically and the new D-optimal designs were more efficient.


---

## Body

## 1. Introduction

Most experimental design literature on causal inference focuses on establishing a causal relationship between variables, but there is no literature on how to identify a design that results in the optimal parameter estimates for a structural equation model (SEM). This work was motivated by an agricultural research study by Milander and colleagues [1] where a SEM was proposed to determine the influence of seeding rate on the yield components of waxy maize and to better determine the interrelationship between yield and yield components by estimating the direct effect of the path model based on the biological understanding of the interrelationship among the endogenous and exogenous variables. The endogenous variables in the proposed system as described in Figure 1 were the yield components rows per ear, ear length, kernels per row, kernels per ear, grain yield, ear circumference, and kernel weight. Seeding rate was the only exogenous variable.Figure 1
The proposed causal model for a maize yield component study published in Maydica by Milander et al. [1].The goodness-of-fit parameters for the model were poor, so the data were reanalyzed to obtain a model that better fit the data, with the resulting model as described in Figure2. However, it did not give any insight into the interrelationship among many of the variables, which raised the question of whether or not the design was adequate for structural equation modeling.Figure 2
The data-supported causal model for the maize yield component study published in Maydica by Milander et al. [1].The study was conducted over a period of three years and the original model was known before the data were collected [1]. Ideally, the researchers should have designed an experiment which would have allowed for the most precise estimates of the parameters, thus leading to more precise inferential statistics. Rather than uniformly applying the seeding rate levels across the experiment, the optimal design would have provided the researchers the number of replicates for each seeding rate level and which of the seeding rate levels should have been replicated more. However, they did not have the tools to obtain an optimal design for the proposed model nor were they able to take advantage of the information that was collected in the first year to adjust the design for subsequent years. There simply was no literature to support such an objective.Historically, what drove the growth of optimal design theory was the need to have a design that is optimal for all parameters that need to be estimated [2, 3]. Technological advancements made it possible to develop the search algorithms that are needed to obtain optimal designs [4–8]. The new growth in the application of experimental design came about because of the lack of theoretical results in the optimal design of blocked and split-plot experiments, which made the computerized design algorithms popular and led to three proposals for point-exchange algorithms for constructing D-optimal blocked designs [9–15].Traditional statistical analyses focus on univariate linear models, which have limitations because most random processes involve multiple dependent variables. This, in turn, led to the development of the multivariate linear model that does not address the causal relationships among the dependent (endogenous) variables, which is not the case with structural equation modeling [16, 17]. SEMs allow qualitative cause-effect information to be combined with statistical data and provide quantitative assessment of the cause-effect relationships among, and within, the endogenous and exogenous variables [17, 18].Even though there have been many algorithms developed for univariate optimal designs and univariate mixed optimal designs, there is very limited literature about multivariate optimal design with examples like the work of Soumaya and fellow authors [19] and Schwabe [20]. There is literature in experimental design in structural equation modeling but not in optimal design. The research in this area so far is to establish causal relationships (direction of the arrows) among the variables [21, 22]. This research extends to counterfactuals and experimental design where it is used to establish causal relationships [23]. However, no research exists to obtain an optimal design to estimate the parameters of a SEM.The gap in literature demonstrates the need for this current research, which comprises three objectives. For a given causal structure, the first objective is to obtain aD-optimal design that allows the most precise estimates of the endogenous and exogenous parameters of that model. Then, for a given mixed causal structure with random blocks or split-plots, the second objective is to obtain a D-optimal design for estimating the parameters of the model. Finally, the efficiency of each of the D-optimal designs for causal structures will be compared with the optimal designs for the univariate case.


## 2. Methods

For the univariate case, the linear model isy¯=Xβ¯+ε¯, where y¯ is an n×1 vector of observations, X is an n×q design matrix of rank q, β¯ is a q×1 vector of unknown coefficients that are estimable, and ε¯ is an n×1 identically independently normally distributed vector with Eε¯=0 and Varε¯=σ2In where In is the identity matrix with rank n. The least-squares estimate of β¯ is β¯∧=X′X−1X′y¯. Our objective is to choose a design “X” such that maxXVarβ∧¯=maxXX′X [7, 8]. Such a design “X” is called a D-optimal design.The same philosophy would be followed for a mixed model. The model for the data that comes from experiments with random blocks or split-plots isy¯=Xβ¯+Zu¯+ε¯ where y¯, X, β¯, and ε¯ are defined above and u¯ is a random vector with Eu¯=0 and Varu¯=σu2Ib. Z is an n×b matrix of the form Z=diag1¯k1,…,1¯kb where 1¯ki is a ki×1 vector of ones. Then, it follows that V=Vary¯=σε2In+σu2ZZ′. Now it can be shown that β˜¯=X′V−1X−1X′V−1y where Varβ˜¯=X′V−1X−1. The objective becomes to choose the design matrix “X” such that maxXX′V−1X. However, V is unknown. This presents a new challenge and shows that the more complicated the model, the more challenges can arise when choosing a D-optimal design.The first solution to the above challenge is to assume thatV=Vary¯=σε2In+σu2ZZ′ is known by making assumptions about σε2 and σu2. The second solution is to rewrite V=Vary¯=σε2In+σu2ZZ′=σε2In+ηZZ′=σε2In+ρ/1−ρZZ′, where η=σu2/σε2 and ρ=σu2/σε2+σu2 [14, 15, 24]. This allows us to rewrite X′V−1X=σε−2qX′I+ηZZ′−1X as the function to optimize a univariate linear model. It can be assumed without loss of generality σε2=1 , thus X′V−1X can be expressed in terms of variance ratio η or in terms of the correlation coefficient ρ. Therefore, the D-optimal design for a blocked or a split-plot experiment depends on the variance ratio. However, Goos and Vandebroek [14, 15] and Goos and Jones [13] argued that the dependence is minor and presented algorithms for constructing locally optimal designs for a given η or ρ. In some special cases, the optimal design is globally optimal and neither depends on the variance ratio nor the degree of correlation of which the most popular is the orthogonal design [14].The previous method, or “traditional method” as it is referred to in literature, is based on optimization ofVarβ˜¯ or X′V−1X and has two weaknesses [25]. The first weakness is that the D-optimal design depends on η or ρ and these values are not known prior to the experiment. The second weakness is that the optimal design focuses only on the fixed effect parameters and ignores the estimate of the random effect parameters. These weaknesses could have severe consequences on the estimate of the variance components, which severely affects the inferential statistics on the fixed effect parameters. Mylona et al. [25] provided scenarios for optimal designs that were constructed using the traditional D-optimality criterion and, in these designs, the random effects σu2 and σε2 are not estimable.Consequently, the two weaknesses led to a new compositeD-optimality criterion. The new composite criterion was proposed by Mylona et al. [25]. The reason for calling it composite is because this criterion takes into consideration both the fixed effect and the random effect parameters.However,D-optimality criterion and orthogonal designs are both dependent on an assumed model. If the assumed model incorrectly represents the interrelationships between the endogenous and exogenous variables, it may lead to an increased bias error. Thus, in cases where there is no knowledge of the design specification, one of the most used designs is the uniform design in which the design points are scattered over the experimental region [26]. In comparison to D-optimal and orthogonal designs, uniform designs are less sensitive for model selection [27, 28]. Consequently, uniform design does not require pre-experiment model specification [29]. Uniform design could be used to explore the relationship between variables.For the same reasons that led to the development of algorithms to account for all parameters of interest in the univariate case, our objective is to develop algorithms for SEMs that consider both the endogenous and exogenous parameters. In structural equation modeling, estimation methods are more complex and there are a variety of class estimators. This research focuses on three-stage least squares (3SLS) and full information maximum likelihood (FIML) methods because they are the only full-system estimation methods where all parameters are estimated simultaneously.


## 3.D-Optimal Design for a Causal Structure Model

### 3.1. The Information Matrix Approximation for 3SLS

In most systems, researchers are interested in more than one endogenous variable and these variables directly and indirectly affect each other. To model such a study, we use SEMs. To demonstrate the methodology and introduce the new approach forD-optimal designs, the following examples will be introduced. A simple example is shown in Figure 3 with two exogenous variables and two endogenous variables.Figure 3
Identified SEM where each exogenous variable directly affects only one endogenous variable.Many phenomena can be modeled based on this SEM. For example, at the University of Washington Tacoma there are two pre-calculus pathways into CalculusI, a two-course sequence (TMATH 115 and TMATH 116) and a one-course sequence (TMATH 120). In this model, x1 represents the grade in TMATH 115, y1 represents the grade in TMATH 116, x2 represents the grade in TMATH 120, and y2 represents the Calculus I grade. If we are designing a study to estimate the parameters in Figure 3, our first objective would be to determine how many students should be enrolled in each pathway to obtain the best estimates for the three parameters γ11, γ22, and b12. There are two parameters that need to be estimated through the first pathway but only one through the second pathway. In this example, assigning half of the sample size to each pathway may not be optimal because there are more parameters on the first pathway than the second. So, the fundamental question here is, how many students should be assigned to each pathway to obtain the optimal estimates for all parameters?We will focus on the model in Figure3. In this centered model, the kth observation can be expressed mathematically as yk1=γ11xk1+εk1 and yk2=b12yk1+γ22xk2+εk2 where εk1εk2∼N0,Σ for k=1…n. Some definitions are necessary to develop the model. Yj is an n×pj−1 matrix that consists of the endogenous independent variables that are in the jth equation. y¯j is an n×1 vector that consists of the responses of the endogenous dependent variable for the jth equation. Xj is an n×qj matrix that consists of the exogenous variables that are in the jth equation, b¯j are the endogenous parameters of the jth equation, and γ¯j are the exogenous parameters of the jth equation. ε¯j is an n×1 vector of random error for the jth endogenous variable. Each equation in the system can be rewritten using vector notations. So, y¯j=Yjb¯j+Xjγ¯j+ε¯j and ε¯j∼N0¯,σεj2In. Now rewrite the jth equation as y¯j=YjXjbj′γj′′+ε¯j where Wj=YjXj and δ¯j=b¯j′γ¯j′. Therefore, y¯j=Wjδ¯j+ε¯j. Multiply both sides by X′, use matrix notation, and use the Kronecker product to obtain I⊗X′y¯∗=I⊗X′W∗δ¯∗+I⊗X′ε¯∗ where y¯∗=y¯1′,…,y¯p′′, W∗=diagW1,…,Wp, δ¯∗=δ¯1′,…,δ¯p′′, and ε¯∗=ε¯1′,…,ε¯p′′ [30].Then, to estimateδ¯∗, use the Generalized Least Square (GLS) to obtain δ∧¯3SLS∗=W∗′Σ−1⊗XX′X−1X′W∗−1W∗′Σ−1⊗XX′X−1X′−1y¯∗. So, the covariance matrix of the 3SLS estimators is Varδ∧¯3SLS∗=W∧∗′Σ−1⊗XX′X−1X′W∧∗−1+O1/n where W∧∗=EW [31, 32], which is equivalent to the asymptotic information matrix for FIML [33]. Thus, the estimate of the covariance matrix would be(1)Var∧δ∧¯∗=W∧∗′Σ∧−1⊗XX′X−1X′W∧∗−1.The reason that we are interested inVar∧δ∧¯∗ is because our objective is to obtain the design X that “minimizes” Var∧δ∧¯∗. Since normality is assumed throughout the research, then Var∧δ∧¯∗−1 is called the information matrix estimate and is denoted by M∧=Var∧δ∧¯∗−1. Thus, our objective would be to obtain the D-optimal design “X” such that(2)maxXM∧=minXVarδ¯∧∗.In this paper, we assume that there are no constraints on the treatments or on the exogenous variables. However, if there are constraints on the treatments, such as if one of the treatments is limited, then the optimality criteria can be adjusted as is done in the univariate case, which is known as the optimal design of mixture experiments problem [13, 34–37].Example 1.
.Based on equation (2), the D-optimal design for the causal structure below (Figure 4) assumes that the treatments are qualitative treatments for twenty observations and assumes that the true values for γ11, γ22, and b12 are 8, 2, and 5, respectively Figure 4:
In vector notation, the model can be expressed asy¯1=γ11x¯1+ε¯1 and y¯2=b12y¯1+γ22x¯2+ε¯2 where ε¯1ε¯2∼N0¯,Σ⊗I20. In this example, y¯1 is a 20×1 vector that includes the responses of the first endogenous variable, y¯2 is a 20×1 vector that includes the responses of the second endogenous variable, x¯1 is a 20×1 vector that indicates whether or not the first treatment is applied on the kth experimental unit as denoted by xk1=1 or xk1=0, x¯2 is a 20×1 vector that indicates whether or not the second treatment is applied on the kth experimental unit as denoted by xk2=1 or xk2=0, ε¯1 is a 20×1 vector of random residuals for the first endogenous variable, ε¯2 is a 20×1 vector of random residuals for the second endogenous variable, 0¯ is a 40×1 vector of zeros, Σ=I2, and I20 is a 20×20 identity matrix. To obtain a D-optimal design, the algorithm starts with multiple random designs in order to avoid local optimality. For this specific example, 50 initial random designs were used, each design with size 20. The initial designs were constructed randomly where each xkj value had a 50–50 chance of being a 1 or 0. To improve the design, the algorithm compares each point in the design with the candidate points 0,0,1,0,0,1,1,1 and makes a simultaneous exchange with the candidate point that improves the optimality criterion from equation (2) the most. The exchanges will continue until no further improvement is achieved or the improvement is sufficiently small. This procedure will be repeated for each of the 50 initial designs. The final designs will be compared. The design X that has the smallest Var∧δ∧¯∗ from among the 50 final designs will be selected as the D-optimal design.
Using the previous algorithm, aD-optimal design for 20 observations for the model in Example 1 is given in Table 1.
In theD-optimal design for the causal structure for the model in Figure 4, the points 1,1 and 1,0 were each replicated nine times, and the point 0,1 was replicated twice. The point 0,0 was not replicated, which is expected because the model was assumed to have no intercept.
Table1 was based on the true value of W∧∗ where the true parameters are used in W∧∗. However, the true values are unknown. There are two approaches that address this issue. First, when prior data are available, estimates of the parameters can be used to estimate Var∧δ∧¯3SLS∗ instead of the true parameters. One question that arises with the use of this approach is the robustness of the design. To address the sensitivity of the true values of the parameters on the design, simulation is used to obtain the optimal design based on the estimates of the parameters.
When there is no prior data available, a second approach is a Bayesian approach which assumes a prior distribution for the endogenous variables, which is similar to the approach that is used by Mylona et al. [25] and Goos and Mylona [38] to address the variance components optimal design for blocked or split-plot experiments. This approach allows for the uncertainty of the endogenous parameters. Using the Bayesian approach to obtain a 3SLS D-optimal design for a causal structure will be the subject of future work.
Optimal design theory is currently limited to univariate or multivariate applications where one of the biggest weaknesses of these approaches is that they ignore any possible causal structure among the endogenous variables. Because of the lack of theoretical results and the algorithms to produce optimal designs for a causal structure, one approach is to use univariate optimal designs for causal models, which leads to a loss of efficiency. To obtain a univariateD-optimal design for Example 1, it is assumed that the treatment will affect both endogenous variables to avoid the situation where there is an optimal design for each equation. In this case, we obtain a D-optimal design for both treatment combinations. An optimal design in the univariate case for two qualitative treatments (i.e.y¯=γ11x¯1+γ22x¯2+ε¯) can be found in Table 2 where the four possible candidate points are 0,0,0,1,1,0, and 1,1. The three candidate points 0,1,1,0, and 1,1 are replicated seven times, seven times, and six times, respectively. For the same reasons as in the D-optimal design for a causal structure, the point 0,0 was not replicated since the intercept was not included in the model.
The ratio of the determinants of the information matrices of the univariate versus the 3SLS designs and theD-efficiency of the two designs are used to compare the optimality of the designs [9, 10]. The results from Table 3 compare the efficiency of the univariate optimal design from Table 2 to the efficiency of the new D-optimal design from Tables 2 and 3.
When comparing the univariate optimal design to the 3SLSD-optimal design for a causal structure, there was an approximately 18% increase in the determinant of the asymptotic information matrix estimates. The univariate optimal design is about 94% as D-efficient as the 3SLS D-optimal design for a causal structure in all four cases.
One criticism over the comparison of theD-optimal designs obtained above is that they are based on the asymptotic information matrices, but the D-optimal designs were obtained for small samples sizes (20 observations). Thus, it is important to ensure that the results are consistent for small samples sizes, as well. To verify the comparison, we simulated data based on the causal structure D-optimal design and data based on the univariate optimal design 100,000 times each. We then estimated the parameters for both designs using the GLS estimate δ∧¯3SLS∗ and calculated the covariance matrix of the parameter estimates for both designs, whose elements are covδ∧i,δ∧j=∑j=1105∑i=1105δi−δ∧iδj−δ∧j/105−1 where δ∧i is the ith element of δ∧¯3SLS∗. The determinants of the covariance matrices of the parameters estimates were computed. This process was repeated three times. The determinant of the covariance matrix of the parameters estimates for the 3SLS causal structure D-optimal design was consistently smaller than the determinant of the covariance matrix of the parameters estimates for the univariate optimal design, as shown in Table 4.
The results from Table4 were consistent with the results for the comparison of the asymptotic information matrices. The 3SLS causal structure D-optimal design consistently produced a smaller determinant for the covariance matrix of the parameters estimates than the 3SLS univariate optimal design. Specifically, the determinant of the covariance matrix of the parameter estimates for the 3SLS causal structure D-optimal design was about 15%–19% smaller than the determinant of the covariance matrix of the parameter estimates for the 3SLS univariate optimal design. Based on those results, the 3SLS causal structure D-optimal design was 5%–7% more D-efficient than the univariate optimal design. It is important to note that whether the true parameters or their estimates are used, the 3SLS D-optimal design for a causal structure did not change.Figure 4
The path model for Example1.Table 1
3SLS and FIMLD-optimal design for Example 1.
#1234567891011121314151617181920x111111111101111101111x200110110110101110100Table 2
Classical optimal design for univariate modely¯=γ11x¯1+γ22x¯2+ε¯.
#1234567891011121314151617181920x101101101101101101101x210110110110110110110Table 3
Parameter values, estimates, determinants of information matrices, their ratios, andD-efficiency based on a 3SLS D-optimal design and a univariate optimal design for three simulations.
Parameters and information matricesParameter valuesSimulation 1 estimatesSimulation 2 estimatesSimulation 3 estimatesγ1187.918.258.34γ2222.310.191.30b1255.015.105.03M=Varδ¯Univariate∗−1110,656108,062.4117,621.26120,228.29M3SLS=Varδ¯3SLS∗−1134,784131,624.88143,268146,443.48M/M3SLS82%82%82%82%D‡−efficiency=M/M3SLS1/394%94%94%94%‡ See [10].Table 4
Determinants of the covariance matrices for the 3SLS causal structureD-optimal design and the 3SLS univariate optimal design based on three simulations each with 100,000 data sets.
1st 100,000 data sets2nd 100,000 data sets3rd 100,000 data setsDeterminant of the covariance matrix for the causal structureD-optimal design7,613,179.57,718,001.37,290,743.0Determinant of the covariance matrix for the univariate optimal design9,008,371.88,893,254.39,035,191.3M/M3SLS85%87%81%D−efficiency=M/M3SLS1/395%95%93%


### 3.2. The Information Matrix Approximation for FIML

The full information maximum likelihood (FIML) is another estimation method for SEMs. In this section, we will discuss optimal designs based on the FIML estimators and will investigate the relationship between the optimal designs.Durbin [33] proposed a transformation for the maximum likelihood equations that simplified the computations and made it easier to study the properties of FIML estimators and their advantages over 3SLS estimators. However, even with the special notations, the FIML information matrix is sophisticated and takes much more work to obtain than the 3SLS information matrix. Based on Durbin’s [33] results, our objective is to obtain a design matrix “X” such that maxXMFIML=maxXQ¯′GQ¯. The derivation and the special notation are noted in Appendix A. The reason that both D-optimal designs are considered is because the 3SLS information matrix is easy to obtain, but FIML estimators have advantages over 3SLS estimators as discussed by Durbin [33].As with the 3SLSD-optimal design, the objective would be to obtain the design matrix X that will minimize the estimate of the determinant of the covariance of the FIML estimates or maximize the determinant of the FIML information matrix estimate. The performance of the D-optimal design for a causal structure will then be compared to the univariate optimal design by computing |Var∧δ∧¯FIML∗| through simulation. The design which consistently has the largest determinants of the information matrix would be the better design. In the following example, the objective will be to use the estimate of the information matrix to obtain the FIML optimal design.For the causal structure in Example1, the 3SLS D-optimal design is also the D-optimal design for the FIML estimates. Similar to 3SLS, the true parameters are unknown, so the FIML estimates are used to replace the true parameters for the three simulations. In practice, the information matrix is estimated based on the parameters estimates. This is similar to the approach of Goos and Vandebroek [14, 15] and Goos and Jones [13] in the case of optimal design for random blocks or split-plots where they suggest using the estimates of the variance parameters or a reasonable guess, arguing that the design minimally depends on those values. The results obtained for a D-optimal design for a causal structure in this research are consistent with those of Goos and Vandebroek [14, 15] and Goos and Jones [13] where the optimal design did not change when replacing the true parameters with their estimates.The results shown in Table5 show the efficiency of the D-optimal design from Table 1 as compared to the efficiency of the univariate optimal design from Table 2.Table 5
Parameter values, estimates, determinants of information matrices, their ratios, andD-efficiency based on a FIML D-optimal design and a univariate optimal design for three simulations.
Parameters and information matricesParameter valuesSimulation 1 estimatesSimulation 2 estimatesSimulation 3 estimatesγ1187.988.208.36γ2222.051.172.04b1255.034.984.98M=Varδ¯Univariate∗−1110,656110,169.65116,153.07120,709.06MFIML=Varδ¯FIML∗−1134,784133,044.67141,479.68147,029.08M/MFIML82%83%82%82%D‡−efficiency=M/MFIML1/394%94%94%94%‡ See [10].The results for the FIML causal structureD-optimal design are similar to the results for the 3SLS causal structure D-optimal design. The optimal design for the univariate case was about 6% less D-efficient than the FIML design for the causal structure. The loss of efficiency demonstrates the importance of taking the endogenous parameters into account to obtain an optimal design to study the causal structure.For both the 3SLS and FIMLD-optimal designs, the goal of the three simulations is to show that the D-optimal design did not change even though the estimate of the parameters deviated from the true value. However, it may be problematic to base a conclusion on three simulations only; therefore, 300 data sets were simulated and then the parameters were estimated. These parameters were used to obtain a D-optimal design based on the estimate of the parameters. In each of the 300 simulations, the D-optimal design did not change and was equivalent to the D-optimal design in Table 1.Another way to demonstrate the robustness of the design is to incrementally change the value of the parameter(s) from the true value, then obtain aD-optimal design based on the value of the parameters, and finally compare the D-optimal design to that which was obtained based on the true value. For the model in Figure 4, the D-optimal design depends only on γ11=8 since W∧∗ includes only y¯1=γ11x¯1+ε¯1. Table 1 is the D-optimal design for γ11=8. Now, we will change the value of γ11 incrementally based on the interval 8±6=14,2. We first start with 8 and then decrease γ11 incrementally by 0.1 (i.e.8,7.9,7.8,…,2). Then, we obtain the D-optimal design based on these values and check whether it is equivalent to the D-optimal design in Table 1. Similarly, we start with 8 and then increase γ11 incrementally by 0.1 (i.e.8,8.1,8.2,…,14). Then, again, we obtain the D-optimal design based on these values and check whether it is equivalent to the D-optimal design in Table 1. The D-optimal designs did not change even though the value of the parameters changed up to 75% from the true value. These results are consistent with Goos and Vandebroek [14, 15] and Goos and Jones [13] where they argued that the optimal design depended minimally on the specific parameter values.


## 3.1. The Information Matrix Approximation for 3SLS

In most systems, researchers are interested in more than one endogenous variable and these variables directly and indirectly affect each other. To model such a study, we use SEMs. To demonstrate the methodology and introduce the new approach forD-optimal designs, the following examples will be introduced. A simple example is shown in Figure 3 with two exogenous variables and two endogenous variables.Figure 3
Identified SEM where each exogenous variable directly affects only one endogenous variable.Many phenomena can be modeled based on this SEM. For example, at the University of Washington Tacoma there are two pre-calculus pathways into CalculusI, a two-course sequence (TMATH 115 and TMATH 116) and a one-course sequence (TMATH 120). In this model, x1 represents the grade in TMATH 115, y1 represents the grade in TMATH 116, x2 represents the grade in TMATH 120, and y2 represents the Calculus I grade. If we are designing a study to estimate the parameters in Figure 3, our first objective would be to determine how many students should be enrolled in each pathway to obtain the best estimates for the three parameters γ11, γ22, and b12. There are two parameters that need to be estimated through the first pathway but only one through the second pathway. In this example, assigning half of the sample size to each pathway may not be optimal because there are more parameters on the first pathway than the second. So, the fundamental question here is, how many students should be assigned to each pathway to obtain the optimal estimates for all parameters?We will focus on the model in Figure3. In this centered model, the kth observation can be expressed mathematically as yk1=γ11xk1+εk1 and yk2=b12yk1+γ22xk2+εk2 where εk1εk2∼N0,Σ for k=1…n. Some definitions are necessary to develop the model. Yj is an n×pj−1 matrix that consists of the endogenous independent variables that are in the jth equation. y¯j is an n×1 vector that consists of the responses of the endogenous dependent variable for the jth equation. Xj is an n×qj matrix that consists of the exogenous variables that are in the jth equation, b¯j are the endogenous parameters of the jth equation, and γ¯j are the exogenous parameters of the jth equation. ε¯j is an n×1 vector of random error for the jth endogenous variable. Each equation in the system can be rewritten using vector notations. So, y¯j=Yjb¯j+Xjγ¯j+ε¯j and ε¯j∼N0¯,σεj2In. Now rewrite the jth equation as y¯j=YjXjbj′γj′′+ε¯j where Wj=YjXj and δ¯j=b¯j′γ¯j′. Therefore, y¯j=Wjδ¯j+ε¯j. Multiply both sides by X′, use matrix notation, and use the Kronecker product to obtain I⊗X′y¯∗=I⊗X′W∗δ¯∗+I⊗X′ε¯∗ where y¯∗=y¯1′,…,y¯p′′, W∗=diagW1,…,Wp, δ¯∗=δ¯1′,…,δ¯p′′, and ε¯∗=ε¯1′,…,ε¯p′′ [30].Then, to estimateδ¯∗, use the Generalized Least Square (GLS) to obtain δ∧¯3SLS∗=W∗′Σ−1⊗XX′X−1X′W∗−1W∗′Σ−1⊗XX′X−1X′−1y¯∗. So, the covariance matrix of the 3SLS estimators is Varδ∧¯3SLS∗=W∧∗′Σ−1⊗XX′X−1X′W∧∗−1+O1/n where W∧∗=EW [31, 32], which is equivalent to the asymptotic information matrix for FIML [33]. Thus, the estimate of the covariance matrix would be(1)Var∧δ∧¯∗=W∧∗′Σ∧−1⊗XX′X−1X′W∧∗−1.The reason that we are interested inVar∧δ∧¯∗ is because our objective is to obtain the design X that “minimizes” Var∧δ∧¯∗. Since normality is assumed throughout the research, then Var∧δ∧¯∗−1 is called the information matrix estimate and is denoted by M∧=Var∧δ∧¯∗−1. Thus, our objective would be to obtain the D-optimal design “X” such that(2)maxXM∧=minXVarδ¯∧∗.In this paper, we assume that there are no constraints on the treatments or on the exogenous variables. However, if there are constraints on the treatments, such as if one of the treatments is limited, then the optimality criteria can be adjusted as is done in the univariate case, which is known as the optimal design of mixture experiments problem [13, 34–37].Example 1.
.Based on equation (2), the D-optimal design for the causal structure below (Figure 4) assumes that the treatments are qualitative treatments for twenty observations and assumes that the true values for γ11, γ22, and b12 are 8, 2, and 5, respectively Figure 4:
In vector notation, the model can be expressed asy¯1=γ11x¯1+ε¯1 and y¯2=b12y¯1+γ22x¯2+ε¯2 where ε¯1ε¯2∼N0¯,Σ⊗I20. In this example, y¯1 is a 20×1 vector that includes the responses of the first endogenous variable, y¯2 is a 20×1 vector that includes the responses of the second endogenous variable, x¯1 is a 20×1 vector that indicates whether or not the first treatment is applied on the kth experimental unit as denoted by xk1=1 or xk1=0, x¯2 is a 20×1 vector that indicates whether or not the second treatment is applied on the kth experimental unit as denoted by xk2=1 or xk2=0, ε¯1 is a 20×1 vector of random residuals for the first endogenous variable, ε¯2 is a 20×1 vector of random residuals for the second endogenous variable, 0¯ is a 40×1 vector of zeros, Σ=I2, and I20 is a 20×20 identity matrix. To obtain a D-optimal design, the algorithm starts with multiple random designs in order to avoid local optimality. For this specific example, 50 initial random designs were used, each design with size 20. The initial designs were constructed randomly where each xkj value had a 50–50 chance of being a 1 or 0. To improve the design, the algorithm compares each point in the design with the candidate points 0,0,1,0,0,1,1,1 and makes a simultaneous exchange with the candidate point that improves the optimality criterion from equation (2) the most. The exchanges will continue until no further improvement is achieved or the improvement is sufficiently small. This procedure will be repeated for each of the 50 initial designs. The final designs will be compared. The design X that has the smallest Var∧δ∧¯∗ from among the 50 final designs will be selected as the D-optimal design.
Using the previous algorithm, aD-optimal design for 20 observations for the model in Example 1 is given in Table 1.
In theD-optimal design for the causal structure for the model in Figure 4, the points 1,1 and 1,0 were each replicated nine times, and the point 0,1 was replicated twice. The point 0,0 was not replicated, which is expected because the model was assumed to have no intercept.
Table1 was based on the true value of W∧∗ where the true parameters are used in W∧∗. However, the true values are unknown. There are two approaches that address this issue. First, when prior data are available, estimates of the parameters can be used to estimate Var∧δ∧¯3SLS∗ instead of the true parameters. One question that arises with the use of this approach is the robustness of the design. To address the sensitivity of the true values of the parameters on the design, simulation is used to obtain the optimal design based on the estimates of the parameters.
When there is no prior data available, a second approach is a Bayesian approach which assumes a prior distribution for the endogenous variables, which is similar to the approach that is used by Mylona et al. [25] and Goos and Mylona [38] to address the variance components optimal design for blocked or split-plot experiments. This approach allows for the uncertainty of the endogenous parameters. Using the Bayesian approach to obtain a 3SLS D-optimal design for a causal structure will be the subject of future work.
Optimal design theory is currently limited to univariate or multivariate applications where one of the biggest weaknesses of these approaches is that they ignore any possible causal structure among the endogenous variables. Because of the lack of theoretical results and the algorithms to produce optimal designs for a causal structure, one approach is to use univariate optimal designs for causal models, which leads to a loss of efficiency. To obtain a univariateD-optimal design for Example 1, it is assumed that the treatment will affect both endogenous variables to avoid the situation where there is an optimal design for each equation. In this case, we obtain a D-optimal design for both treatment combinations. An optimal design in the univariate case for two qualitative treatments (i.e.y¯=γ11x¯1+γ22x¯2+ε¯) can be found in Table 2 where the four possible candidate points are 0,0,0,1,1,0, and 1,1. The three candidate points 0,1,1,0, and 1,1 are replicated seven times, seven times, and six times, respectively. For the same reasons as in the D-optimal design for a causal structure, the point 0,0 was not replicated since the intercept was not included in the model.
The ratio of the determinants of the information matrices of the univariate versus the 3SLS designs and theD-efficiency of the two designs are used to compare the optimality of the designs [9, 10]. The results from Table 3 compare the efficiency of the univariate optimal design from Table 2 to the efficiency of the new D-optimal design from Tables 2 and 3.
When comparing the univariate optimal design to the 3SLSD-optimal design for a causal structure, there was an approximately 18% increase in the determinant of the asymptotic information matrix estimates. The univariate optimal design is about 94% as D-efficient as the 3SLS D-optimal design for a causal structure in all four cases.
One criticism over the comparison of theD-optimal designs obtained above is that they are based on the asymptotic information matrices, but the D-optimal designs were obtained for small samples sizes (20 observations). Thus, it is important to ensure that the results are consistent for small samples sizes, as well. To verify the comparison, we simulated data based on the causal structure D-optimal design and data based on the univariate optimal design 100,000 times each. We then estimated the parameters for both designs using the GLS estimate δ∧¯3SLS∗ and calculated the covariance matrix of the parameter estimates for both designs, whose elements are covδ∧i,δ∧j=∑j=1105∑i=1105δi−δ∧iδj−δ∧j/105−1 where δ∧i is the ith element of δ∧¯3SLS∗. The determinants of the covariance matrices of the parameters estimates were computed. This process was repeated three times. The determinant of the covariance matrix of the parameters estimates for the 3SLS causal structure D-optimal design was consistently smaller than the determinant of the covariance matrix of the parameters estimates for the univariate optimal design, as shown in Table 4.
The results from Table4 were consistent with the results for the comparison of the asymptotic information matrices. The 3SLS causal structure D-optimal design consistently produced a smaller determinant for the covariance matrix of the parameters estimates than the 3SLS univariate optimal design. Specifically, the determinant of the covariance matrix of the parameter estimates for the 3SLS causal structure D-optimal design was about 15%–19% smaller than the determinant of the covariance matrix of the parameter estimates for the 3SLS univariate optimal design. Based on those results, the 3SLS causal structure D-optimal design was 5%–7% more D-efficient than the univariate optimal design. It is important to note that whether the true parameters or their estimates are used, the 3SLS D-optimal design for a causal structure did not change.Figure 4
The path model for Example1.Table 1
3SLS and FIMLD-optimal design for Example 1.
#1234567891011121314151617181920x111111111101111101111x200110110110101110100Table 2
Classical optimal design for univariate modely¯=γ11x¯1+γ22x¯2+ε¯.
#1234567891011121314151617181920x101101101101101101101x210110110110110110110Table 3
Parameter values, estimates, determinants of information matrices, their ratios, andD-efficiency based on a 3SLS D-optimal design and a univariate optimal design for three simulations.
Parameters and information matricesParameter valuesSimulation 1 estimatesSimulation 2 estimatesSimulation 3 estimatesγ1187.918.258.34γ2222.310.191.30b1255.015.105.03M=Varδ¯Univariate∗−1110,656108,062.4117,621.26120,228.29M3SLS=Varδ¯3SLS∗−1134,784131,624.88143,268146,443.48M/M3SLS82%82%82%82%D‡−efficiency=M/M3SLS1/394%94%94%94%‡ See [10].Table 4
Determinants of the covariance matrices for the 3SLS causal structureD-optimal design and the 3SLS univariate optimal design based on three simulations each with 100,000 data sets.
1st 100,000 data sets2nd 100,000 data sets3rd 100,000 data setsDeterminant of the covariance matrix for the causal structureD-optimal design7,613,179.57,718,001.37,290,743.0Determinant of the covariance matrix for the univariate optimal design9,008,371.88,893,254.39,035,191.3M/M3SLS85%87%81%D−efficiency=M/M3SLS1/395%95%93%


## 3.2. The Information Matrix Approximation for FIML

The full information maximum likelihood (FIML) is another estimation method for SEMs. In this section, we will discuss optimal designs based on the FIML estimators and will investigate the relationship between the optimal designs.Durbin [33] proposed a transformation for the maximum likelihood equations that simplified the computations and made it easier to study the properties of FIML estimators and their advantages over 3SLS estimators. However, even with the special notations, the FIML information matrix is sophisticated and takes much more work to obtain than the 3SLS information matrix. Based on Durbin’s [33] results, our objective is to obtain a design matrix “X” such that maxXMFIML=maxXQ¯′GQ¯. The derivation and the special notation are noted in Appendix A. The reason that both D-optimal designs are considered is because the 3SLS information matrix is easy to obtain, but FIML estimators have advantages over 3SLS estimators as discussed by Durbin [33].As with the 3SLSD-optimal design, the objective would be to obtain the design matrix X that will minimize the estimate of the determinant of the covariance of the FIML estimates or maximize the determinant of the FIML information matrix estimate. The performance of the D-optimal design for a causal structure will then be compared to the univariate optimal design by computing |Var∧δ∧¯FIML∗| through simulation. The design which consistently has the largest determinants of the information matrix would be the better design. In the following example, the objective will be to use the estimate of the information matrix to obtain the FIML optimal design.For the causal structure in Example1, the 3SLS D-optimal design is also the D-optimal design for the FIML estimates. Similar to 3SLS, the true parameters are unknown, so the FIML estimates are used to replace the true parameters for the three simulations. In practice, the information matrix is estimated based on the parameters estimates. This is similar to the approach of Goos and Vandebroek [14, 15] and Goos and Jones [13] in the case of optimal design for random blocks or split-plots where they suggest using the estimates of the variance parameters or a reasonable guess, arguing that the design minimally depends on those values. The results obtained for a D-optimal design for a causal structure in this research are consistent with those of Goos and Vandebroek [14, 15] and Goos and Jones [13] where the optimal design did not change when replacing the true parameters with their estimates.The results shown in Table5 show the efficiency of the D-optimal design from Table 1 as compared to the efficiency of the univariate optimal design from Table 2.Table 5
Parameter values, estimates, determinants of information matrices, their ratios, andD-efficiency based on a FIML D-optimal design and a univariate optimal design for three simulations.
Parameters and information matricesParameter valuesSimulation 1 estimatesSimulation 2 estimatesSimulation 3 estimatesγ1187.988.208.36γ2222.051.172.04b1255.034.984.98M=Varδ¯Univariate∗−1110,656110,169.65116,153.07120,709.06MFIML=Varδ¯FIML∗−1134,784133,044.67141,479.68147,029.08M/MFIML82%83%82%82%D‡−efficiency=M/MFIML1/394%94%94%94%‡ See [10].The results for the FIML causal structureD-optimal design are similar to the results for the 3SLS causal structure D-optimal design. The optimal design for the univariate case was about 6% less D-efficient than the FIML design for the causal structure. The loss of efficiency demonstrates the importance of taking the endogenous parameters into account to obtain an optimal design to study the causal structure.For both the 3SLS and FIMLD-optimal designs, the goal of the three simulations is to show that the D-optimal design did not change even though the estimate of the parameters deviated from the true value. However, it may be problematic to base a conclusion on three simulations only; therefore, 300 data sets were simulated and then the parameters were estimated. These parameters were used to obtain a D-optimal design based on the estimate of the parameters. In each of the 300 simulations, the D-optimal design did not change and was equivalent to the D-optimal design in Table 1.Another way to demonstrate the robustness of the design is to incrementally change the value of the parameter(s) from the true value, then obtain aD-optimal design based on the value of the parameters, and finally compare the D-optimal design to that which was obtained based on the true value. For the model in Figure 4, the D-optimal design depends only on γ11=8 since W∧∗ includes only y¯1=γ11x¯1+ε¯1. Table 1 is the D-optimal design for γ11=8. Now, we will change the value of γ11 incrementally based on the interval 8±6=14,2. We first start with 8 and then decrease γ11 incrementally by 0.1 (i.e.8,7.9,7.8,…,2). Then, we obtain the D-optimal design based on these values and check whether it is equivalent to the D-optimal design in Table 1. Similarly, we start with 8 and then increase γ11 incrementally by 0.1 (i.e.8,8.1,8.2,…,14). Then, again, we obtain the D-optimal design based on these values and check whether it is equivalent to the D-optimal design in Table 1. The D-optimal designs did not change even though the value of the parameters changed up to 75% from the true value. These results are consistent with Goos and Vandebroek [14, 15] and Goos and Jones [13] where they argued that the optimal design depended minimally on the specific parameter values.


## 4.D-Optimal Design for a Causal Structure Model with Random Blocks

In applied research, blocked designs are likely the most commonly used experimental designs since they can be used to account for variation attributable to sources other than treatments and can considerably improve the precision of the experiment. The first objective is to establish the blocked causal structure and then to obtain the 3SLS and FIML estimators for the endogenous and exogenous parameters, which will be used to estimate their information matrices. The second objective of this section is to use the algorithm described in Section1.1 in Example 1 to obtain a D-optimal design for both the endogenous and exogenous parameters. The last objective is to compare the efficiency of the D-optimal design to the classical univariate mixed model optimal design.


### 4.1. 3SLS D-Optimal Design for a Causal Structure with Random Blocks

The model with random blocks is similar to the model for a completely randomized design except that the block effect needs to be added to every endogenous variable. Therefore, for theith endogenous variable as the dependent variable, the model can be written as. y¯i=Yib¯i+Xiγ¯i+Zu¯i+ε¯i where u¯iu¯j∼N0¯,⌈σui2σuijσuijσuj2⌉⊗Ib and all of the other terms have been previously defined. If the blocks have the same size, then Z=diag1k,1k,…,1k=Ib×b⊗1k, the elements of u¯i and ε¯i are assumed to be mutually independent and normally distributed with zero mean and variances σui2 and σi2, respectively.LetWi=YiXi and δ¯i=b¯i′γ¯i′′ . Then, the model with random blocks can be rewritten as y¯i=Wiδ¯i+Zu¯i+ε¯i. Use matrix notation and then the Kronecker product to rewrite the previous model as I⊗X′y¯∗=I⊗X′W∗δ¯∗+I⊗X′Zu¯∗+I⊗X′ε¯∗.Therefore,δ∧3SLS∗=I⊗X′W∗′V−1I⊗X′W∗−1I⊗X′W∗′V−1I⊗X′y¯∗ where V=VarI⊗X′Zu¯∗+I⊗X′ε¯∗. Let Varu¯∗=Σu⊗Ib and Varε¯∗=Σ⊗In where Σu is a p×p nonnegative definite matrix. Then, VarI⊗X′Zu¯∗+I⊗X′ε¯∗=I⊗X′ZVaru¯∗I⊗X′Z′+I⊗X′Varε¯∗I⊗X′′=I⊗X′ZΣu⊗IbI⊗X′Z′+I⊗X′Σ⊗InI⊗X′′=Z∗G∗Z∗′+R∗ where z∗=I⊗X′Z and G∗=Σu⊗Ib.Using the same approach as described in Section1.1, the asymptotic information matrix would be Μ3SLS=I⊗X′W∧∗′V−1I⊗X′W∧∗ where V=Z∗G∗Z∗′+R∗.Based on the information matrix for the 3SLS estimates for the causal structureD-optimal design with random blocks, our objective would be to use the algorithm described in Section 1.1 in Example 1 to obtain the design X such that maxXM3SLS=maxXI⊗X′W∧∗′V−1I⊗X′W∧∗. Since the information matrix depends on the unknown parameters through Σ∧ and W∧∗, we face the same predicament as in the case of the completely randomized D-optimal design for a causal structure where the parameters are unknown. The solutions that were proposed for the completely randomized D-optimal design for a causal structure in Section 1.1 are appropriate for the causal structure D-optimal design with random blocks. Other approaches are also possible. For example, in an experiment with random blocks that is conducted in multiple sequential stages, the data of one block (year, lab, day, etc.) could be utilized to obtain initial estimates that are then used to obtain the D-optimal design for the rest of the blocks (other years, other labs, other days, etc.). Requiring preliminary information is a standard part of univariate design. If there is no prior data, then a reasonable guess can be used similar to the approach of Goos and Vandebroek [14, 15] and Goos and Jones [13]. In many experiments, multiple factors and limitations of the experimental material require the use of split-plot designs and/or incomplete block designs. Specifically, the way the factors are applied and the nature of the experimental material give rise to different sizes of experimental units and limits on block size. For this reason, an incomplete block design with three factors is used as an example in this section. In the next example, the objective is to use the asymptotic information matrix estimates to obtain a D-optimal design for a causal structure with four random blocks where the size of each block is four.Example 2.
.Assume that there are three binary treatment factors (exogenous variables) with two endogenous variables for the causal structure given in Figure5. Here, the x¯i, i=1,2,3 each with a xkj value, had a 50–50 chance of being a 1 or −1. Also, assume that the true values for γ11, γ22, γ31, and b12 are 8, 2, 3, and 5, respectively Figure 5:
AD-optimal design for a causal structure with four random blocks with block size equal to four under the given assumptions can be found in Table 6.
The expected determinant of the information matrix is 5,222,400 and the previous design is not orthogonal. The expected determinant of the information matrix of the orthogonal design in Table7 is 4,784,128. The orthogonal design is the optimal design for the univariate mixed model with three binary treatment factors and four blocks each block with size four.
One concern for bothD-optimal and orthogonal designs is that they both require an assumed model structure. One popular alternative is the uniform design, which is “a space filling design” ([39], p. 131) that uniformly distributes experimental points across the design domain without model pre-specification [39]. Therefore, we will compare the information of the uniform design to the orthogonal design and the new D-optimal design. Using a uniform design table [39, 40], a uniform design for 16 runs with 3 factors and an experimental cube of −1,13 is shown in Tables 7 and 8.
The expected determinant of the information matrix for the uniform design is 28,612, which is significantly smaller than the determinants of theD-optimal and the orthogonal design information matrices.
The results show that the newD-optimal design was approximately 2.2% more D-efficient than the orthogonal design, whereas the D-optimal design was 72.8% more D-efficient than the uniform design. The D-optimal design performed significantly better than the uniform design but is still comparable to the orthogonal design. Although orthogonal designs are desirable because they are universally optimal for first-order linear models and are also robust for nonlinear models, they also come with concerns in that they sometimes do not exist and require a large sample size [25]. It is also important to keep in mind that in Example 2, there is only one endogenous parameter. The significance in the increase in efficiency of the determinant may not be easily identified through our simple example. In a more complex design, there will be more endogenous variables and if those variables are interconnected in a more sophisticated way, then the increase of the determinant of the information matrix would be amplified.
It is worth noting the differences between the 3SLS causal structureD-optimal design and the orthogonal design. First, the orthogonal design completely confounds the three-way interaction xi1xj2xk3 with the blocks where xi1 is the ith element of x¯1, xj2 is the jth element of x¯2, and xk3 is the kth element of x¯3. The number of replicates for each candidate point is 2. However, this is not the case in the 3SLS causal structure D-optimal design where the two-way interaction of xj2xk3 is confounded with Blocks 1 and 3, the three-way interaction of xi1xj2xk3 is confounded with Block 4, and the two-way interaction of xi1xk3 is confounded with Block 2. In the 3SLS causal structure D-optimal design, the number of candidate points is not replicated equally. For example, the candidate point −1,−1,−1 was replicated three times, but the candidate point 1,−1,−1 was replicated only once.
To understand the differences between the designs, it is important to contrast the asymptotic information matrix of the 3SLS D-optimal design for a causal structure with that of the orthogonal design. In this way, we can see where information is gained and parameter estimates are improved for the former design and where information is sacrificed and parameter estimates worsen for the latter design. Tables9 and 10 show the asymptotic information matrices for the model in Figure 5 for the orthogonal design and the 3SLS D-optimal design, respectively Table 9.
When the information of the estimates is compared, we can see that the exogenous parameters in both models have the same information except for the elements.
γ11,γ31. However, the endogenous parameters in the 3SLS D-optimal design gained more information with an increase of about 16% over the orthogonal design. These results are expected because univariate optimal designs do not take the precision of the endogenous parameters into account, whereas D-optimal designs for causal structures do. Therefore, our new approach has the advantage of giving the optimal combination of treatments for the entire model and taking into account the precision of both the exogenous and endogenous parameters.Figure 5
The identified path model for Example2.Table 6
A 3SLS causal structureD-optimal design with three binary factors and random blocks where each block size is equal to four.
Block 1Block 2Block 3Block 4x1x2x3x1x2x3x1x2x3x1x2x3111111−11−1−1−1−1−1−1−1−1−1−1−1−11−1111−1−1−11−11−111−11−1111−1111−111−1Table 7
An orthogonal design for three binary treatment factors and four blocks of size four.
Block 1Block 2Block 3Block 4x1x2x3x1x2x3x1x2x3x1x2x3−1111−1−1−11−1−1−1−11−11−1−111−1−1−11111−1111−1−111−11−1−1−1−11−111111−1Table 8
A uniform design for 16 runs with 3 factors and an experimental cube of−1,13.
#123x1x2x311127−17/15−3/152591−7/151/15−13344−11/15−9/15−9/15411225/15−13/15−13/15541415−9/1511/1513/15613769/15−3/15−5/1572612−13/15−5/157/1581413311/159/15−11/1591531413/15−11/1511/15101011133/155/159/15111610913/151/15128816−1/15−1/151139581/15−7/15−1/15147165−3/151−7/15156110−5/15−13/15161215117/1513/155/15Table 9
Asymptotic information matrix for Example2 based on the orthogonal design.
γ11γ22γ31b12γ1116000γ2201600γ3100160b120001,168Table 10
Asymptotic information matrix for Example2 based on the 3SLS D-optimal design.
γ11γ22γ31b12γ1116040γ2201600γ3140160b120001,360


### 4.2. FIML D-Optimal Design for a Causal Structure with Random Blocks

As in the completely randomizedD-optimal design for a causal structure, another important class of estimators is the maximum likelihood estimators. We will obtain the FIML estimators for a causal structure with random blocks and then obtain the covariance matrix estimates for the estimators, or the information matrix estimates. The information matrix estimates will then be used to obtain a D-optimal design for the FIML estimators.Starting with the notation from Section2.1, y¯∗−W∗δ¯∗=Z∗u¯∗+ε¯∗ where u¯∗ε¯∗∼N0¯0¯,Σu⊗Ib00Σ⊗In. Then, ℒV,δ¯∗|y¯∗,W=−np/2log2π−1/2logV−1/2y¯∗−W∗δ¯∗′V−1y¯∗−W∗δ¯∗. To maximize the likelihood, ∂ℒ/∂δ¯∗=W∗′V−1y¯∗−W∗′V−1W∗δ∧¯∗ and set the partial derivative equal to 0, which gives W∗′V−1y¯∗−W∗′V−1W∗δ∧¯∗=0. Solving for δ∧¯∗, makes δ∧¯∗=W∗′V−1W∗−1W∗′V−1y¯∗. Then, find the second partial derivative and take its expected value to obtain the information matrix, ∂2ℒ/∂δ¯∗′∂δ¯∗=−W∗′V−1W∗. Next, Var−1δ∧¯∗=EW∗′V−1W∗. However, in general, we cannot find the expected value explicitly until the model is known since W∗ includes fixed and random variables. But, for Example 2 above the information matrix is simplified to [30].Based on the information matrix above, a FIMLD-optimal design for Example 2 is equivalent to the design in Table 6. The determinant of the information matrix is 5,280,295.4. Neither the orthogonal design nor the uniform design was as efficient as the new D-optimal design. The determinants of the information matrices were 4,845,883.1 and 36,622 for the orthogonal design and uniform design, respectively. Thus, the D-optimal design was approximately 2.1% more D-efficient than the orthogonal design and 71.1% more D-efficient than the uniform design. These results are similar to the D-optimal design for Section 1.2 where both the 3SLS and the FIML estimates had the same D-optimal design.


## 4.1. 3SLS D-Optimal Design for a Causal Structure with Random Blocks

The model with random blocks is similar to the model for a completely randomized design except that the block effect needs to be added to every endogenous variable. Therefore, for theith endogenous variable as the dependent variable, the model can be written as. y¯i=Yib¯i+Xiγ¯i+Zu¯i+ε¯i where u¯iu¯j∼N0¯,⌈σui2σuijσuijσuj2⌉⊗Ib and all of the other terms have been previously defined. If the blocks have the same size, then Z=diag1k,1k,…,1k=Ib×b⊗1k, the elements of u¯i and ε¯i are assumed to be mutually independent and normally distributed with zero mean and variances σui2 and σi2, respectively.LetWi=YiXi and δ¯i=b¯i′γ¯i′′ . Then, the model with random blocks can be rewritten as y¯i=Wiδ¯i+Zu¯i+ε¯i. Use matrix notation and then the Kronecker product to rewrite the previous model as I⊗X′y¯∗=I⊗X′W∗δ¯∗+I⊗X′Zu¯∗+I⊗X′ε¯∗.Therefore,δ∧3SLS∗=I⊗X′W∗′V−1I⊗X′W∗−1I⊗X′W∗′V−1I⊗X′y¯∗ where V=VarI⊗X′Zu¯∗+I⊗X′ε¯∗. Let Varu¯∗=Σu⊗Ib and Varε¯∗=Σ⊗In where Σu is a p×p nonnegative definite matrix. Then, VarI⊗X′Zu¯∗+I⊗X′ε¯∗=I⊗X′ZVaru¯∗I⊗X′Z′+I⊗X′Varε¯∗I⊗X′′=I⊗X′ZΣu⊗IbI⊗X′Z′+I⊗X′Σ⊗InI⊗X′′=Z∗G∗Z∗′+R∗ where z∗=I⊗X′Z and G∗=Σu⊗Ib.Using the same approach as described in Section1.1, the asymptotic information matrix would be Μ3SLS=I⊗X′W∧∗′V−1I⊗X′W∧∗ where V=Z∗G∗Z∗′+R∗.Based on the information matrix for the 3SLS estimates for the causal structureD-optimal design with random blocks, our objective would be to use the algorithm described in Section 1.1 in Example 1 to obtain the design X such that maxXM3SLS=maxXI⊗X′W∧∗′V−1I⊗X′W∧∗. Since the information matrix depends on the unknown parameters through Σ∧ and W∧∗, we face the same predicament as in the case of the completely randomized D-optimal design for a causal structure where the parameters are unknown. The solutions that were proposed for the completely randomized D-optimal design for a causal structure in Section 1.1 are appropriate for the causal structure D-optimal design with random blocks. Other approaches are also possible. For example, in an experiment with random blocks that is conducted in multiple sequential stages, the data of one block (year, lab, day, etc.) could be utilized to obtain initial estimates that are then used to obtain the D-optimal design for the rest of the blocks (other years, other labs, other days, etc.). Requiring preliminary information is a standard part of univariate design. If there is no prior data, then a reasonable guess can be used similar to the approach of Goos and Vandebroek [14, 15] and Goos and Jones [13]. In many experiments, multiple factors and limitations of the experimental material require the use of split-plot designs and/or incomplete block designs. Specifically, the way the factors are applied and the nature of the experimental material give rise to different sizes of experimental units and limits on block size. For this reason, an incomplete block design with three factors is used as an example in this section. In the next example, the objective is to use the asymptotic information matrix estimates to obtain a D-optimal design for a causal structure with four random blocks where the size of each block is four.Example 2.
.Assume that there are three binary treatment factors (exogenous variables) with two endogenous variables for the causal structure given in Figure5. Here, the x¯i, i=1,2,3 each with a xkj value, had a 50–50 chance of being a 1 or −1. Also, assume that the true values for γ11, γ22, γ31, and b12 are 8, 2, 3, and 5, respectively Figure 5:
AD-optimal design for a causal structure with four random blocks with block size equal to four under the given assumptions can be found in Table 6.
The expected determinant of the information matrix is 5,222,400 and the previous design is not orthogonal. The expected determinant of the information matrix of the orthogonal design in Table7 is 4,784,128. The orthogonal design is the optimal design for the univariate mixed model with three binary treatment factors and four blocks each block with size four.
One concern for bothD-optimal and orthogonal designs is that they both require an assumed model structure. One popular alternative is the uniform design, which is “a space filling design” ([39], p. 131) that uniformly distributes experimental points across the design domain without model pre-specification [39]. Therefore, we will compare the information of the uniform design to the orthogonal design and the new D-optimal design. Using a uniform design table [39, 40], a uniform design for 16 runs with 3 factors and an experimental cube of −1,13 is shown in Tables 7 and 8.
The expected determinant of the information matrix for the uniform design is 28,612, which is significantly smaller than the determinants of theD-optimal and the orthogonal design information matrices.
The results show that the newD-optimal design was approximately 2.2% more D-efficient than the orthogonal design, whereas the D-optimal design was 72.8% more D-efficient than the uniform design. The D-optimal design performed significantly better than the uniform design but is still comparable to the orthogonal design. Although orthogonal designs are desirable because they are universally optimal for first-order linear models and are also robust for nonlinear models, they also come with concerns in that they sometimes do not exist and require a large sample size [25]. It is also important to keep in mind that in Example 2, there is only one endogenous parameter. The significance in the increase in efficiency of the determinant may not be easily identified through our simple example. In a more complex design, there will be more endogenous variables and if those variables are interconnected in a more sophisticated way, then the increase of the determinant of the information matrix would be amplified.
It is worth noting the differences between the 3SLS causal structureD-optimal design and the orthogonal design. First, the orthogonal design completely confounds the three-way interaction xi1xj2xk3 with the blocks where xi1 is the ith element of x¯1, xj2 is the jth element of x¯2, and xk3 is the kth element of x¯3. The number of replicates for each candidate point is 2. However, this is not the case in the 3SLS causal structure D-optimal design where the two-way interaction of xj2xk3 is confounded with Blocks 1 and 3, the three-way interaction of xi1xj2xk3 is confounded with Block 4, and the two-way interaction of xi1xk3 is confounded with Block 2. In the 3SLS causal structure D-optimal design, the number of candidate points is not replicated equally. For example, the candidate point −1,−1,−1 was replicated three times, but the candidate point 1,−1,−1 was replicated only once.
To understand the differences between the designs, it is important to contrast the asymptotic information matrix of the 3SLS D-optimal design for a causal structure with that of the orthogonal design. In this way, we can see where information is gained and parameter estimates are improved for the former design and where information is sacrificed and parameter estimates worsen for the latter design. Tables9 and 10 show the asymptotic information matrices for the model in Figure 5 for the orthogonal design and the 3SLS D-optimal design, respectively Table 9.
When the information of the estimates is compared, we can see that the exogenous parameters in both models have the same information except for the elements.
γ11,γ31. However, the endogenous parameters in the 3SLS D-optimal design gained more information with an increase of about 16% over the orthogonal design. These results are expected because univariate optimal designs do not take the precision of the endogenous parameters into account, whereas D-optimal designs for causal structures do. Therefore, our new approach has the advantage of giving the optimal combination of treatments for the entire model and taking into account the precision of both the exogenous and endogenous parameters.Figure 5
The identified path model for Example2.Table 6
A 3SLS causal structureD-optimal design with three binary factors and random blocks where each block size is equal to four.
Block 1Block 2Block 3Block 4x1x2x3x1x2x3x1x2x3x1x2x3111111−11−1−1−1−1−1−1−1−1−1−1−1−11−1111−1−1−11−11−111−11−1111−1111−111−1Table 7
An orthogonal design for three binary treatment factors and four blocks of size four.
Block 1Block 2Block 3Block 4x1x2x3x1x2x3x1x2x3x1x2x3−1111−1−1−11−1−1−1−11−11−1−111−1−1−11111−1111−1−111−11−1−1−1−11−111111−1Table 8
A uniform design for 16 runs with 3 factors and an experimental cube of−1,13.
#123x1x2x311127−17/15−3/152591−7/151/15−13344−11/15−9/15−9/15411225/15−13/15−13/15541415−9/1511/1513/15613769/15−3/15−5/1572612−13/15−5/157/1581413311/159/15−11/1591531413/15−11/1511/15101011133/155/159/15111610913/151/15128816−1/15−1/151139581/15−7/15−1/15147165−3/151−7/15156110−5/15−13/15161215117/1513/155/15Table 9
Asymptotic information matrix for Example2 based on the orthogonal design.
γ11γ22γ31b12γ1116000γ2201600γ3100160b120001,168Table 10
Asymptotic information matrix for Example2 based on the 3SLS D-optimal design.
γ11γ22γ31b12γ1116040γ2201600γ3140160b120001,360


## 4.2. FIML D-Optimal Design for a Causal Structure with Random Blocks

As in the completely randomizedD-optimal design for a causal structure, another important class of estimators is the maximum likelihood estimators. We will obtain the FIML estimators for a causal structure with random blocks and then obtain the covariance matrix estimates for the estimators, or the information matrix estimates. The information matrix estimates will then be used to obtain a D-optimal design for the FIML estimators.Starting with the notation from Section2.1, y¯∗−W∗δ¯∗=Z∗u¯∗+ε¯∗ where u¯∗ε¯∗∼N0¯0¯,Σu⊗Ib00Σ⊗In. Then, ℒV,δ¯∗|y¯∗,W=−np/2log2π−1/2logV−1/2y¯∗−W∗δ¯∗′V−1y¯∗−W∗δ¯∗. To maximize the likelihood, ∂ℒ/∂δ¯∗=W∗′V−1y¯∗−W∗′V−1W∗δ∧¯∗ and set the partial derivative equal to 0, which gives W∗′V−1y¯∗−W∗′V−1W∗δ∧¯∗=0. Solving for δ∧¯∗, makes δ∧¯∗=W∗′V−1W∗−1W∗′V−1y¯∗. Then, find the second partial derivative and take its expected value to obtain the information matrix, ∂2ℒ/∂δ¯∗′∂δ¯∗=−W∗′V−1W∗. Next, Var−1δ∧¯∗=EW∗′V−1W∗. However, in general, we cannot find the expected value explicitly until the model is known since W∗ includes fixed and random variables. But, for Example 2 above the information matrix is simplified to [30].Based on the information matrix above, a FIMLD-optimal design for Example 2 is equivalent to the design in Table 6. The determinant of the information matrix is 5,280,295.4. Neither the orthogonal design nor the uniform design was as efficient as the new D-optimal design. The determinants of the information matrices were 4,845,883.1 and 36,622 for the orthogonal design and uniform design, respectively. Thus, the D-optimal design was approximately 2.1% more D-efficient than the orthogonal design and 71.1% more D-efficient than the uniform design. These results are similar to the D-optimal design for Section 1.2 where both the 3SLS and the FIML estimates had the same D-optimal design.


## 5. Discussion

Optimal design theory and its applications are generally limited to univariate designs. However, most natural phenomena include multiple endogenous variables which are related in complicated ways. There are two weaknesses concerning univariate designs. First, univariate designs are incapable of addressing more than one endogenous variable that could be affected by different exogenous variables. The second weakness is that univariate designs ignore the endogenous parameters. Our new approach overcomes both weaknesses. TheD-optimality criteria were developed for both the 3SLS and FIML estimates for a causal structure and for a causal structure with random blocks. In addition, our approach is broadly applicable since it generalizes to any linear causal model with both fixed and random effects, and with any type of nesting and factorial structure.Our results show that theD-optimal designs were the same for both FIML and 3SLS estimates. The reader may wonder why they should care about 3SLS D-optimal design since FIML estimators have traditionally been used in structural equation modeling since 3SLS estimators are not based on optimization. However, 3SLS information matrices are mathematically easier to obtain than FIML information matrices. The interest here is in optimal design, not estimation. Since there was no difference between the FIML and 3SLS D-optimal designs for both examples and, in general, since the FIML and 3SLS information matrices are asymptotically equivalent [33], we use the 3SLS information matrix both for its simplicity in obtaining a D-optimal design and to ease hesitancy for scientists to utilize the more complicated FIML information matrix.In the example of the completely randomized causal structure, the newD-optimal design increased the information matrix determinant by at least 17% over the univariate optimal design. Because of this increase, the univariate optimal design was only 94% as D-efficient as the new D-optimal design. The information matrix of a causal structure depends on parameter values that are unknown. This raises the concern of whether or not the D-optimal design is robust. Based on the simulation results from Example 1, the D-optimal design is robust because it did not change when the parameter values changed.The fact that the information matrix of a causal structure depends on unknown parameter values is not unique to causal structureD-optimal designs. The same issue arises for D-optimal designs for random effects, or the covariance parameters, in a split-plot or a blocked design as well where the values of the parameters are needed to obtain the D-optimal design. The results in our research are consistent with those of Goos and Vandebroek [14, 15] and Goos and Jones [13] where the optimal design did not change when replacing the true parameters for an estimate.Similarly, for a causal structure model with random blocks, the endogenous parameters had a significant impact on theD-optimal design. In the univariate mixed model, orthogonal designs are universally optimal [14, 41]. However, it is not optimal for a causal structure model with random blocks. As shown in Example 2, the new D-optimal design increased the determinant of the information matrix by at least 9%, meaning that the orthogonal design was less efficient than the new D-optimal design. The orthogonal design was consistently about 97.8% and 97.9% as D-efficient as the new 3SLS and FIML D-optimal designs, respectively.While it may be a disadvantage thatD-optimal designs require an assumed model structure, the results show that the uniform design, which focuses on distributing the design points uniformly over the experimental domain, is not optimal. The uniform design was about 27.2% as D-efficient as the new 3SLS D-optimal design and 28.9% as D-efficient as the new FIML D-optimal design. In both cases, the uniform designs consistently lost more than 70% of the information in contrast to the new D-optimal design. Such a large loss of information will result in inflating the variance and has serious consequences on inferential statistics and the power of the statistical tests.The designs that were produced using the new criteria show that endogenous parameters have an important effect onD-optimal designs. The results further show the importance of taking advantage of prior knowledge of the interrelationships of endogenous and exogenous parameters to develop a D-optimal design using the algorithms and methodology developed here. An application for the use of our algorithms could be in research similar to the heart study by Mi and colleagues [42] who examined the complex interrelationships among physiological, genetic, and environmental factors and their interaction effects on coronary heart disease (CHD), as shown in Figure 6. The data in that study were collected and analysis was performed based on the data alone; there was no pre-planned design. To obtain the most precise estimates to describe the system, causal structure D-optimal design can be used to determine replication of the levels of each factor such as how many men and women to include in the study, how many of the subjects should be smokers, how many should be obese, or how to determine an optimal sample of volunteers based on budgetary constraints.Figure 6
Path estimates of SEM of gene-by-environment interaction in the development of CHD from annals of human genetics by Mi et al. [42].Another practical use of our algorithms and methodology comes from the field of agriculture through research conducted by Campbell and colleagues [43], and Dhungana and colleagues [44] where the objective was to understand how gene-environment interactions (GEI) influence interrelated and complex traits such as grain yield. This subject of future work would help to answer questions like which genotypes should be chosen for the study and what is the combination of those genotypes to obtain optimal estimates. As a multi-stage experiment which took place over a period of three years, our methods could have been implemented to use the preliminary results from the first year of the study to answer those questions if the model was pre-determined.


### 5.1. Limitations

The fact that the algorithm requires the model to be predetermined is a challenge because in most natural phenomena the interrelationships in a system could be complex. In univariate designs, the challenge is reduced because there is only one endogenous variable that is directly affected by an exogenous variable. However, in structural equation modeling, not all exogenous variables affect endogenous variables directly and the interrelationship among the endogenous variables needs to be predetermined, as well. An additional challenge to obtain aD-optimal design through our algorithm is that some of the estimates of the endogenous parameters are required before an experiment is conducted, which is unknown prior to the experiment. Addressing these limitations is also the subject of future work by using a prior distribution and a Bayesian approach to allow for the uncertainty of the unknown endogenous and exogenous parameters. Finally, how our technique works with other designs such as A-optimal, C-optimal, or E-optimal designs is yet to be explored.


## 5.1. Limitations

The fact that the algorithm requires the model to be predetermined is a challenge because in most natural phenomena the interrelationships in a system could be complex. In univariate designs, the challenge is reduced because there is only one endogenous variable that is directly affected by an exogenous variable. However, in structural equation modeling, not all exogenous variables affect endogenous variables directly and the interrelationship among the endogenous variables needs to be predetermined, as well. An additional challenge to obtain aD-optimal design through our algorithm is that some of the estimates of the endogenous parameters are required before an experiment is conducted, which is unknown prior to the experiment. Addressing these limitations is also the subject of future work by using a prior distribution and a Bayesian approach to allow for the uncertainty of the unknown endogenous and exogenous parameters. Finally, how our technique works with other designs such as A-optimal, C-optimal, or E-optimal designs is yet to be explored.



---

*Source: 7299086-2022-08-30.xml*