| Title: | Multi-Objective Optimal Design of Experiments |
|---|---|
| Description: | Provides functionality to generate compound optimal designs for targeting the multiple experimental objectives directly, ensuring that the full set of research questions is answered as economically as possible. Designs can be found using point or coordinate exchange algorithms combining estimation, inference and lack-of-fit criteria that account for model inadequacy. Details and examples are given by Koutra et al. (2024) <doi:10.48550/arXiv.2412.17158>. |
| Authors: | Vasiliki Koutra [aut, cre, cph] (ORCID: <https://orcid.org/0000-0002-1117-4554>), Olga Egorova [aut, cph], Steven Gilmour [aut, cph], Luzia Trinca [aut, cph] |
| Maintainer: | Vasiliki Koutra <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 1.1.0 |
| Built: | 2026-05-28 08:15:07 UTC |
| Source: | https://github.com/vkstats/moode |
This function forms the full extended candidate set with all polynomial terms, with labels, not orthonormalised.
candidate_set_full(cand.trt, K)candidate_set_full(cand.trt, K)
cand.trt |
Candidate set of treatments, the first column contains treatment labels. Usually obtained as output from the candidate_trt_set function. |
K |
Number of factors. |
The full extended candidate set: the column of treatment labels, then named columns with polynomial terms up to the 4th order.
For example, "x12" stands for , and "x1x2 stands for , and "x23x4" for .
# Full extended candidate set for two 3-level factors K<-2; Levels <- rep(list(1:3),K); candidate_set_full(candidate_trt_set(Levels, K), K)# Full extended candidate set for two 3-level factors K<-2; Levels <- rep(list(1:3),K); candidate_set_full(candidate_trt_set(Levels, K), K)
This function forms the full extended orthonormalised candidate set of primary and potential terms, with an intercept column and labels.
candidate_set_orth(cand.full, primary.terms, potential.terms)candidate_set_orth(cand.full, primary.terms, potential.terms)
cand.full |
Candidate set containing terms up to 4th order, with labels in the first column. |
primary.terms |
Character vector identifying primary model terms. |
potential.terms |
Character vector identifying potential model terms. |
The orthonormalised full candidate set containing primary and potential terms, with labels.
# Full extended orthonormalised candidate set for two 4-level factors, # full quadratic polynomial model as primary model and all three-order terms as potential. K<-2; Levels <- rep(list(1:4),K) cand.trt <- candidate_trt_set(Levels, K) cand.full <- candidate_set_full(cand.trt, K) prime.terms <- colnames(cand.full)[2:7] poten.terms <- colnames(cand.full)[8:11] Parameters <- c(1, rep(1,K), rep(1,K), K*(K-1)/2) candidate_set_orth(cand.full, prime.terms, poten.terms)# Full extended orthonormalised candidate set for two 4-level factors, # full quadratic polynomial model as primary model and all three-order terms as potential. K<-2; Levels <- rep(list(1:4),K) cand.trt <- candidate_trt_set(Levels, K) cand.full <- candidate_set_full(cand.trt, K) prime.terms <- colnames(cand.full)[2:7] poten.terms <- colnames(cand.full)[8:11] Parameters <- c(1, rep(1,K), rep(1,K), K*(K-1)/2) candidate_set_orth(cand.full, prime.terms, poten.terms)
This function forms the candidate set of treatments from the factors' levels, adds labels, with optional spherical transformation of the coordinates.
candidate_trt_set(Levels, K, Hypercube = TRUE)candidate_trt_set(Levels, K, Hypercube = TRUE)
Levels |
Levels of each factor. |
K |
Number of factors. |
Hypercube |
Indicates if the experimental region is a hypercube ('TRUE') or spherical ('FALSE'). |
Matrix of candidate set of treatments, with treatment labels contained in the first column.
# Candidate treatment set for five 3-level factors K<-5; Levels <- rep(list(1:3),K); candidate_trt_set(Levels, K)# Candidate treatment set for five 3-level factors K<-5; Levels <- rep(list(1:3),K); candidate_trt_set(Levels, K)
This function evaluates the Generalised Ds-criterion (Goos et al. 2005) for given primary and potential model matrices.
criteria.GD(X1, X2, search.object, eps = 1e-23)criteria.GD(X1, X2, search.object, eps = 1e-23)
X1 |
The primary model matrix, with the first column containing the labels of treatments, and the second – the intercept term. |
X2 |
The matrix of potential terms, with the first column containing the labels of treatments. |
search.object |
Object of class |
eps |
Computational tolerance, the default value is 10^-23 |
A list of values: indicator of whether the evaluation was successful ("eval"), Ds-criterion value – intercept excluded ("Ds"), Lack-of-fit criterion value ("LoF"), the bias component value ("bias"), the number of pure error degrees of freedom ("df") and the value of the compound criterion ("compound").
Goos P, Kobilinsky A, O'Brien TE, Vandebroek M (2005). “Model-Robust and Model-Sensitive Designs.” Computational Statistics and Data Analysis, 49, 201-216.
#Experiment: one 5-level factor, primary model -- full quadratic, one potential (cubic) term # setting up the example ex.mood <- mood(K = 1, Levels = 5, Nruns = 7, criterion.choice = "GDP", kappa = list(kappa.Ds = 1./3, kappa.LoF = 1./3, kappa.bias = 1./3), model_terms = list(primary.model = "second_order", potential.model = "cubic_terms")) # Generating candidate set: orthonormalised K <- ex.mood$K Levels <- ex.mood$Levels cand.not.orth <- candidate_set_full(candidate_trt_set(Levels, K), K) cand.full.orth <- candidate_set_orth(cand.not.orth, ex.mood$primary.terms, ex.mood$potential.terms) # Choosing a design index <- c(rep(1, 2), 3, 4, rep(5, 3)) X.primary <- cand.full.orth[index, c(1, match(ex.mood$primary.terms, colnames(cand.full.orth)))] X.potential <- cand.full.orth[index, (c(1, match(ex.mood$potential.terms, colnames(cand.full.orth))))] # Evaluating a compound GD-criterion criteria.GD(X1 = X.primary, X2 = X.potential, ex.mood) # Output: eval = 1, Ds = 0.7334291, LoF = 0.7212544, bias = 1.473138, df = 3, compound = 0.9202307#Experiment: one 5-level factor, primary model -- full quadratic, one potential (cubic) term # setting up the example ex.mood <- mood(K = 1, Levels = 5, Nruns = 7, criterion.choice = "GDP", kappa = list(kappa.Ds = 1./3, kappa.LoF = 1./3, kappa.bias = 1./3), model_terms = list(primary.model = "second_order", potential.model = "cubic_terms")) # Generating candidate set: orthonormalised K <- ex.mood$K Levels <- ex.mood$Levels cand.not.orth <- candidate_set_full(candidate_trt_set(Levels, K), K) cand.full.orth <- candidate_set_orth(cand.not.orth, ex.mood$primary.terms, ex.mood$potential.terms) # Choosing a design index <- c(rep(1, 2), 3, 4, rep(5, 3)) X.primary <- cand.full.orth[index, c(1, match(ex.mood$primary.terms, colnames(cand.full.orth)))] X.potential <- cand.full.orth[index, (c(1, match(ex.mood$potential.terms, colnames(cand.full.orth))))] # Evaluating a compound GD-criterion criteria.GD(X1 = X.primary, X2 = X.potential, ex.mood) # Output: eval = 1, Ds = 0.7334291, LoF = 0.7212544, bias = 1.473138, df = 3, compound = 0.9202307
This function evaluates the Generalised DPs-criterion for given primary and potential model matrices. Components: Ds-, DPs-, LoF(DP)- and Bias(D)-optimality.
criteria.GDP(X1, X2, search.object, eps = 10^-23)criteria.GDP(X1, X2, search.object, eps = 10^-23)
X1 |
The primary model matrix, with the first column containing the labels of treatments, and the second – the intercept term. |
X2 |
The matrix of potential terms, with the first column containing the labels of treatments. |
search.object |
Object of class |
eps |
Computational tolerance, the default value is 10^-23 |
A list of values: indicator of whether the evaluation was successful ("eval"), Ds-criterion value – intercept excluded ("Ds"), DPs-criterion value – intercept excluded ("DPs"), Lack-of-fit(DP) criterion value ("LoF"), the bias component value ("bias"), the number of pure error degrees of freedom ("df") and the value of the compound criterion ("compound").
# Experiment: one 5-level factor, primary model -- full quadratic, X^3 and X^4 potential terms. ex.mood <- mood(K = 1, Levels = 5, Nruns = 8, criterion.choice = "GDP", kappa = list(kappa.Ds = .25, kappa.LoF = .25, kappa.bias = .25, kappa.DP = .25), model_terms = list(primary.model = "second_order", potential.terms = "x14")) # Generating candidate sets: primary and full orthonormalised ones K <- ex.mood$K Levels <- ex.mood$Levels cand.not.orth <- candidate_set_full(candidate_trt_set(Levels, K), K) cand.full.orth <- candidate_set_orth(cand.not.orth, ex.mood$primary.terms, ex.mood$potential.terms) # Choosing a design index <- c(rep(1,2),3,rep(4,2),rep(5,3)) X.primary <- cand.full.orth[index, c(1, match(ex.mood$primary.terms, colnames(cand.full.orth)))] X.potential <- cand.full.orth[index, (c(1, match(ex.mood$potential.terms, colnames(cand.full.orth))))] # Evaluating a compound GDP-criterion criteria.GDP(X1 = X.primary, X2 = X.potential, ex.mood) # Output: eval = 1, Ds = 0.6884783, DP = 4.4538023, LoF = 3.895182, # bias = 1.03807, df = 4, compound = 2.465318# Experiment: one 5-level factor, primary model -- full quadratic, X^3 and X^4 potential terms. ex.mood <- mood(K = 1, Levels = 5, Nruns = 8, criterion.choice = "GDP", kappa = list(kappa.Ds = .25, kappa.LoF = .25, kappa.bias = .25, kappa.DP = .25), model_terms = list(primary.model = "second_order", potential.terms = "x14")) # Generating candidate sets: primary and full orthonormalised ones K <- ex.mood$K Levels <- ex.mood$Levels cand.not.orth <- candidate_set_full(candidate_trt_set(Levels, K), K) cand.full.orth <- candidate_set_orth(cand.not.orth, ex.mood$primary.terms, ex.mood$potential.terms) # Choosing a design index <- c(rep(1,2),3,rep(4,2),rep(5,3)) X.primary <- cand.full.orth[index, c(1, match(ex.mood$primary.terms, colnames(cand.full.orth)))] X.potential <- cand.full.orth[index, (c(1, match(ex.mood$potential.terms, colnames(cand.full.orth))))] # Evaluating a compound GDP-criterion criteria.GDP(X1 = X.primary, X2 = X.potential, ex.mood) # Output: eval = 1, Ds = 0.6884783, DP = 4.4538023, LoF = 3.895182, # bias = 1.03807, df = 4, compound = 2.465318
This function evaluates the Generalised L-criterion (Goos et al. 2005) for given primary and potential model matrices.
criteria.GL(X1, X2, search.object, eps = 10^-23)criteria.GL(X1, X2, search.object, eps = 10^-23)
X1 |
The primary model matrix, with the first column containing the labels of treatments, and the second – the intercept term. |
X2 |
The matrix of potential terms, with the first column containing the labels of treatments. |
search.object |
Object of class |
eps |
Computational tolerance, the default value is 10^-23 |
A list of values: indicator of whether the evaluation was successful ("eval"), Ls-criterion value – intercept excluded ("Ls"), Lack-of-fit criterion value ("LoF"), the bias component value ("bias"), the number of pure error degrees of freedom ("df") and the value of the compound criterion ("compound").
Goos P, Kobilinsky A, O'Brien TE, Vandebroek M (2005). “Model-Robust and Model-Sensitive Designs.” Computational Statistics and Data Analysis, 49, 201-216.
#Experiment: one 5-level factor, primary model -- full quadratic, one potential (cubic) term # setting up the example ex.mood <- mood(K = 1, Levels = 5, Nruns = 7, criterion.choice = "GL", kappa = list(kappa.L = 1./3, kappa.LoF = 1./3, kappa.bias = 1./3), model_terms = list(primary.model = "second_order", potential.model = "cubic_terms")) # Generating candidate set: orthonormalised K <- ex.mood$K Levels <- ex.mood$Levels cand.not.orth <- candidate_set_full(candidate_trt_set(Levels, K), K) cand.full.orth <- candidate_set_orth(cand.not.orth, ex.mood$primary.terms, ex.mood$potential.terms) # Choosing a design index <- c(rep(1, 2), 3, 4, rep(5, 3)) X.primary <- cand.full.orth[index, c(1, match(ex.mood$primary.terms, colnames(cand.full.orth)))] X.potential <- cand.full.orth[index, (c(1, match(ex.mood$potential.terms, colnames(cand.full.orth))))] # Evaluating a compound GD-criterion criteria.GL(X1 = X.primary, X2 = X.potential, ex.mood) # Output: eval = 1, L = 0.3118626, LoF = 0.7212544, bias = 1.473138, df = 3, compound = 0.6919878#Experiment: one 5-level factor, primary model -- full quadratic, one potential (cubic) term # setting up the example ex.mood <- mood(K = 1, Levels = 5, Nruns = 7, criterion.choice = "GL", kappa = list(kappa.L = 1./3, kappa.LoF = 1./3, kappa.bias = 1./3), model_terms = list(primary.model = "second_order", potential.model = "cubic_terms")) # Generating candidate set: orthonormalised K <- ex.mood$K Levels <- ex.mood$Levels cand.not.orth <- candidate_set_full(candidate_trt_set(Levels, K), K) cand.full.orth <- candidate_set_orth(cand.not.orth, ex.mood$primary.terms, ex.mood$potential.terms) # Choosing a design index <- c(rep(1, 2), 3, 4, rep(5, 3)) X.primary <- cand.full.orth[index, c(1, match(ex.mood$primary.terms, colnames(cand.full.orth)))] X.potential <- cand.full.orth[index, (c(1, match(ex.mood$potential.terms, colnames(cand.full.orth))))] # Evaluating a compound GD-criterion criteria.GL(X1 = X.primary, X2 = X.potential, ex.mood) # Output: eval = 1, L = 0.3118626, LoF = 0.7212544, bias = 1.473138, df = 3, compound = 0.6919878
This function evaluates the Generalised LP-criterion for given primary and potential model matrices. Components: L-, LP-, LoF(LP)- and Bias(L)-optimality.
criteria.GLP(X1, X2, search.object, eps = 1e-23)criteria.GLP(X1, X2, search.object, eps = 1e-23)
X1 |
The primary model matrix, with the first column containing the labels of treatments, and the second – the intercept term. |
X2 |
The matrix of potential terms, with the first column containing the labels of treatments. |
search.object |
Object of class |
eps |
Computational tolerance, the default value is 10^-23 |
A list of values: indicator of whether the evaluation was successful ("eval"), Ls-criterion value – intercept excluded ("L"), LP-criterion value – intercept excluded ("LP"), Lack-of-fit(LP) criterion value ("LoF"), the bias component value ("bias"), the number of pure error degrees of freedom ("df") and the value of the compound criterion ("compound").
#' # Experiment: one 5-level factor, primary model -- full quadratic, X^3 and X^4 potential terms. ex.mood <- mood(K = 1, Levels = 5, Nruns = 8, criterion.choice = "GLP", kappa = list(kappa.L = .25, kappa.LoF = .25, kappa.bias = .25, kappa.LP = .25), model_terms = list(primary.model = "second_order", potential.terms = "x14")) # Generating candidate sets: primary and full orthonormalised ones K <- ex.mood$K Levels <- ex.mood$Levels cand.not.orth <- candidate_set_full(candidate_trt_set(Levels, K), K) cand.full.orth <- candidate_set_orth(cand.not.orth, ex.mood$primary.terms, ex.mood$potential.terms) # Choosing a design index <- c(rep(1,2),3,rep(4,2),rep(5,3)) X.primary <- cand.full.orth[index, c(1, match(ex.mood$primary.terms, colnames(cand.full.orth)))] X.potential <- cand.full.orth[index, (c(1, match(ex.mood$potential.terms, colnames(cand.full.orth))))] # Evaluating a compound GDP-criterion criteria.GLP(X1 = X.primary, X2 = X.potential, ex.mood) # Output: eval = 1, L = 0.2952603, LP = 4.584705, LoF = 3.895182, # bias = 1.03807, df = 4, compound = 1.529564#' # Experiment: one 5-level factor, primary model -- full quadratic, X^3 and X^4 potential terms. ex.mood <- mood(K = 1, Levels = 5, Nruns = 8, criterion.choice = "GLP", kappa = list(kappa.L = .25, kappa.LoF = .25, kappa.bias = .25, kappa.LP = .25), model_terms = list(primary.model = "second_order", potential.terms = "x14")) # Generating candidate sets: primary and full orthonormalised ones K <- ex.mood$K Levels <- ex.mood$Levels cand.not.orth <- candidate_set_full(candidate_trt_set(Levels, K), K) cand.full.orth <- candidate_set_orth(cand.not.orth, ex.mood$primary.terms, ex.mood$potential.terms) # Choosing a design index <- c(rep(1,2),3,rep(4,2),rep(5,3)) X.primary <- cand.full.orth[index, c(1, match(ex.mood$primary.terms, colnames(cand.full.orth)))] X.potential <- cand.full.orth[index, (c(1, match(ex.mood$potential.terms, colnames(cand.full.orth))))] # Evaluating a compound GDP-criterion criteria.GLP(X1 = X.primary, X2 = X.potential, ex.mood) # Output: eval = 1, L = 0.2952603, LP = 4.584705, LoF = 3.895182, # bias = 1.03807, df = 4, compound = 1.529564
This function evaluates the MSE DPs-criterion for given primary and potential model matrices. Candidate full model matrices do not have to be orthonormalised. Components: DPs-, LoF(DP)- and MSE(D)-optimality.
criteria.mseD(X1, X2, search.object, eps = 1e-23)criteria.mseD(X1, X2, search.object, eps = 1e-23)
X1 |
The primary model matrix, with the first column containing the labels of treatments, and the second – the intercept term. |
X2 |
The matrix of potential terms, with the first column containing the labels of treatments. |
search.object |
Object of class |
eps |
Computational tolerance, the default value is 10^-23 |
A list of values: indicator of whether the evaluation was successful ("eval"), DPs-criterion value – intercept excluded ("DP"), Lack-of-fit(DP) criterion value ("LoF"), the MSE(D) component value ("mse"), the number of pure error degrees of freedom ("df") and the value of the compound criterion ("compound").
# Experiment: one 5-level factor, primary model -- full quadratic, X^3 and X^4 potential terms. set.seed(20210930) ex.mood <- mood(K = 1, Levels = 5, Nruns = 8, criterion.choice = "MSE.D", kappa = list(kappa.DP = 1./3, kappa.LoF = 1./3, kappa.mse = 1./3), control = list(Biter = 1000), model_terms = list(primary.model = "second_order", potential.terms = "x14")) # Generating candidate sets: primary and full orthonormalised ones K <- ex.mood$K Levels <- ex.mood$Levels cand.not.orth <- candidate_set_full(candidate_trt_set(Levels, K), K) cand.full.orth <- candidate_set_orth(cand.not.orth, ex.mood$primary.terms, ex.mood$potential.terms) # Choosing a design index <- c(rep(1,2),3,rep(4,2),rep(5,3)) X.primary <- cand.full.orth[index, c(1, match(ex.mood$primary.terms, colnames(cand.full.orth)))] X.potential <- cand.full.orth[index, (c(1, match(ex.mood$potential.terms, colnames(cand.full.orth))))] # Evaluating a compound GDP-criterion criteria.mseD(X.primary, X.potential, ex.mood) # Output: eval = 1, DP = 4.538023, LoF = 3.895182, mse = 0.6986903, df = 4, compound = 2.310728# Experiment: one 5-level factor, primary model -- full quadratic, X^3 and X^4 potential terms. set.seed(20210930) ex.mood <- mood(K = 1, Levels = 5, Nruns = 8, criterion.choice = "MSE.D", kappa = list(kappa.DP = 1./3, kappa.LoF = 1./3, kappa.mse = 1./3), control = list(Biter = 1000), model_terms = list(primary.model = "second_order", potential.terms = "x14")) # Generating candidate sets: primary and full orthonormalised ones K <- ex.mood$K Levels <- ex.mood$Levels cand.not.orth <- candidate_set_full(candidate_trt_set(Levels, K), K) cand.full.orth <- candidate_set_orth(cand.not.orth, ex.mood$primary.terms, ex.mood$potential.terms) # Choosing a design index <- c(rep(1,2),3,rep(4,2),rep(5,3)) X.primary <- cand.full.orth[index, c(1, match(ex.mood$primary.terms, colnames(cand.full.orth)))] X.potential <- cand.full.orth[index, (c(1, match(ex.mood$potential.terms, colnames(cand.full.orth))))] # Evaluating a compound GDP-criterion criteria.mseD(X.primary, X.potential, ex.mood) # Output: eval = 1, DP = 4.538023, LoF = 3.895182, mse = 0.6986903, df = 4, compound = 2.310728
This function evaluates the MSE LP-criterion for given primary and potential model matrices. Candidate full model matrices do not have to be orthonormalised. Components: LP-, LoF(LP)- and MSE(L)-optimality.
criteria.mseL(X1, X2, search.object, eps = 10^-23)criteria.mseL(X1, X2, search.object, eps = 10^-23)
X1 |
The primary model matrix, with the first column containing the labels of treatments, and the second – the intercept term. |
X2 |
The matrix of potential terms, with the first column containing the labels of treatments. |
search.object |
Object of class |
eps |
Computational tolerance, the default value is 10^-23 |
A list of values: indicator of whether the evaluation was successful ("eval"), LP-criterion value – intercept excluded ("LP"), Lack-of-fit(LP) criterion value ("LoF"), the MSE(L) component value ("mse"), the number of pure error degrees of freedom ("df") and the value of the compound criterion ("compound").
#'# Experiment: one 5-level factor, primary model -- full quadratic, X^3 and X^4 potential terms. ex.mood <- mood(K = 1, Levels = 5, Nruns = 8, criterion.choice = "MSE.L", kappa = list(kappa.LP = 1./3, kappa.LoF = 1./3, kappa.mse = 1./3), model_terms = list(primary.model = "second_order", potential.terms = "x14")) # Generating candidate sets: primary and full orthonormalised ones K <- ex.mood$K Levels <- ex.mood$Levels cand.not.orth <- candidate_set_full(candidate_trt_set(Levels, K), K) cand.full.orth <- candidate_set_orth(cand.not.orth, ex.mood$primary.terms, ex.mood$potential.terms) # Choosing a design index <- c(rep(1,2),3,rep(4,2),rep(5,3)) X.primary <- cand.full.orth[index, c(1, match(ex.mood$primary.terms, colnames(cand.full.orth)))] X.potential <- cand.full.orth[index, (c(1, match(ex.mood$potential.terms, colnames(cand.full.orth))))] # Evaluating a compound GDP-criterion criteria.mseL(X.primary, X.potential, ex.mood) # Output: eval = 1, LP = 4.584705, LoF = 3.895182, mse = 0.3926842, df = 4, compound = 1.914084#'# Experiment: one 5-level factor, primary model -- full quadratic, X^3 and X^4 potential terms. ex.mood <- mood(K = 1, Levels = 5, Nruns = 8, criterion.choice = "MSE.L", kappa = list(kappa.LP = 1./3, kappa.LoF = 1./3, kappa.mse = 1./3), model_terms = list(primary.model = "second_order", potential.terms = "x14")) # Generating candidate sets: primary and full orthonormalised ones K <- ex.mood$K Levels <- ex.mood$Levels cand.not.orth <- candidate_set_full(candidate_trt_set(Levels, K), K) cand.full.orth <- candidate_set_orth(cand.not.orth, ex.mood$primary.terms, ex.mood$potential.terms) # Choosing a design index <- c(rep(1,2),3,rep(4,2),rep(5,3)) X.primary <- cand.full.orth[index, c(1, match(ex.mood$primary.terms, colnames(cand.full.orth)))] X.potential <- cand.full.orth[index, (c(1, match(ex.mood$potential.terms, colnames(cand.full.orth))))] # Evaluating a compound GDP-criterion criteria.mseL(X.primary, X.potential, ex.mood) # Output: eval = 1, LP = 4.584705, LoF = 3.895182, mse = 0.3926842, df = 4, compound = 1.914084
This function evaluates the MSE DPs-criterion for given primary and potential model matrices, using point MSE(D)-component estimation. Candidate full model matrices do not have to be orthonormalised. Components: DPs-, LoF(DP)- and MSE(D)-optimality.
criteria.mseP(X1, X2, search.object, eps = 10^-23)criteria.mseP(X1, X2, search.object, eps = 10^-23)
X1 |
The primary model matrix, with the first column containing the labels of treatments, and the second – the intercept term. |
X2 |
The matrix of potential terms, with the first column containing the labels of treatments. |
search.object |
Object of class |
eps |
Computational tolerance, the default value is 10^-23 |
A list of values: indicator of whether the evaluation was successful ("eval"), DPs-criterion value – intercept excluded ("DP"), Lack-of-fit(DP) criterion value ("LoF"), the MSE(D) component value ("mse"), the number of pure error degrees of freedom ("df") and the value of the compound criterion ("compound").
# Experiment: one 5-level factor, primary model -- full quadratic, X^3 and X^4 potential terms. ex.mood <- mood(K = 1, Levels = 5, Nruns = 8, criterion.choice = "MSE.P", kappa = list(kappa.DP = 1./3, kappa.LoF = 1./3, kappa.mse = 1./3), model_terms = list(primary.model = "second_order", potential.terms = "x14")) # Generating candidate sets: primary and full orthonormalised K <- ex.mood$K Levels <- ex.mood$Levels cand.not.orth <- candidate_set_full(candidate_trt_set(Levels, K), K) cand.full.orth <- candidate_set_orth(cand.not.orth, ex.mood$primary.terms, ex.mood$potential.terms) # Choosing a design index <- c(rep(1,2),3,rep(4,2),rep(5,3)) X.primary <- cand.full.orth[index, c(1, match(ex.mood$primary.terms, colnames(cand.full.orth)))] X.potential <- cand.full.orth[index, (c(1, match(ex.mood$potential.terms, colnames(cand.full.orth))))] # Evaluating a compound GDP-criterion criteria.mseP(X.primary, X.potential, ex.mood) # Output: eval = 1, DP = 4.538023, LoF = 3.895182, mse = 0.6992699, df = 4, compound = 2.312135# Experiment: one 5-level factor, primary model -- full quadratic, X^3 and X^4 potential terms. ex.mood <- mood(K = 1, Levels = 5, Nruns = 8, criterion.choice = "MSE.P", kappa = list(kappa.DP = 1./3, kappa.LoF = 1./3, kappa.mse = 1./3), model_terms = list(primary.model = "second_order", potential.terms = "x14")) # Generating candidate sets: primary and full orthonormalised K <- ex.mood$K Levels <- ex.mood$Levels cand.not.orth <- candidate_set_full(candidate_trt_set(Levels, K), K) cand.full.orth <- candidate_set_orth(cand.not.orth, ex.mood$primary.terms, ex.mood$potential.terms) # Choosing a design index <- c(rep(1,2),3,rep(4,2),rep(5,3)) X.primary <- cand.full.orth[index, c(1, match(ex.mood$primary.terms, colnames(cand.full.orth)))] X.potential <- cand.full.orth[index, (c(1, match(ex.mood$potential.terms, colnames(cand.full.orth))))] # Evaluating a compound GDP-criterion criteria.mseP(X.primary, X.potential, ex.mood) # Output: eval = 1, DP = 4.538023, LoF = 3.895182, mse = 0.6992699, df = 4, compound = 2.312135
Calculating values of determinant- and trace-based components of Generalized D-, DP-, L- and LP- criteria for an output of a search object, with model and control parameters set in a mood object.
criteria.values.G(search.obj, mood.obj, eps = 10^-23)criteria.values.G(search.obj, mood.obj, eps = 10^-23)
search.obj |
Output of the ‘Search’ function |
mood.obj |
Output of the ‘mood’ function |
eps |
Computational tolerance, default 10^-20 |
List of the calculated values:
df pure error degrees of freedom
Ds Ds-criterion value, intercept excluded
DP DPs-criterion value, intercept excluded
LoFD LoF(D)-criterion value from the GD-criterion
LoFDP LoF(DP)-criterion value from the GDP-criterion
biasD bias(D)-criterion value from the GD-criterion
Ls Ls-criterion value, intercept excluded
LP LPs-criterion value, intercept excluded
LoFL LoF(L)-criterion value from the GL-criterion
LoFLP LoF(LP)-criterion value from the GLP-criterion
biasL bias(L)-criterion value from the GL-criterion
Egorova,~O. (2017).
Optimal Design of Experiments for Multiple Objectives.
Ph.D. thesis, University of Southampton.
Goos P, Kobilinsky A, O'Brien TE, Vandebroek M (2005).
“Model-Robust and Model-Sensitive Designs.”
Computational Statistics and Data Analysis, 49, 201-216.
Calculating values of determinant- and trace-based components of MSE(D)- and MSE(L)- criteria for an output of a search object, with model and control parameters set in a mood object.
criteria.values.mse(search.obj, mood.obj, eps = 10^-20, Biter = 1000)criteria.values.mse(search.obj, mood.obj, eps = 10^-20, Biter = 1000)
search.obj |
Output of the ‘Search’ function |
mood.obj |
Output of the ‘mood’ function |
eps |
Computational tolerance, default 10^-20 |
Biter |
MC sample size for evaluating the mse(D)-component |
List of the calculated values:
df pure error degrees of freedom
Ds Ds-criterion value, intercept excluded
DP DPs-criterion value, intercept excluded
LoFDP LoF(DP)-criterion value
mseD mse(D)-criterion value, obtained via MC sampling
mseP mse(D)-criterion value, obtained using point prior
L L-criterion value, intercept excluded
LP LP-criterion value, intercept excluded
LoFLP LoF(LP)-criterion value
mseL mse(L)-criterion value
Creates an object containing the parameters of the experiment, compound optimality criterion with the weights and parameters of the search.
mood( K, Levels, Nruns, criterion.choice = c("GD", "GL", "GDP", "GLP", "MSE.D", "MSE.L", "MSE.P"), kappa = list(), control = list(), prob = list(), model_terms = list(primary.model = "first_order") )mood( K, Levels, Nruns, criterion.choice = c("GD", "GL", "GDP", "GLP", "MSE.D", "MSE.L", "MSE.P"), kappa = list(), control = list(), prob = list(), model_terms = list(primary.model = "first_order") )
K |
Number of factors. |
Levels |
Either (a) a common number of levels for each factor or (b) a list of length K of the vectors containing levels of each factor. |
Nruns |
Number of runs of the experiment. |
criterion.choice |
Compound criterion to be used for the optimal design search or evaluation. Possible values are:
|
kappa |
List specifying the weights used in the compound criterion. Each named entry must be between 0 and 1.
|
control |
Named list specifying control parameters for the design search.
|
prob |
Named list specifying confidence levels for DP- ( |
model_terms |
A list specifying the primary (fitted) and potential (biased) models with the following named elements (see Details).
|
The function provides different ways of specifying the levels of the factors and the models. Although some default options are provided
for, e.g., criterion.choice and kappa.* values, specification of these input parameters should be carefully chosen to reflect the aims of the experiment and available prior information.
Specifying the factors and levels
If all K factors have the same number of levels, Levels parameter is used to input that number.
Otherwise, Levels is set to be a list of vectors containing the values of the factors, e.g.
list(1:3, 1:2, 1:4) for 3 factors with equally spaced and levels respectively.
Specifying the fitted model and the potential terms
There are two ways to describe the primary and potential sets of model terms via the model_terms list.
Named elements primary.model and potential.model can be used to specify the fitted model and the potential terms via a string form.
They are used to generate the sets of primary.terms and potential.terms which alternatively can be input directly.
Possible values of primary.model are:
main_effects – main effects for all the factors (default for all criteria)
first_order – main effects and linear interactions
second_order – full second order polynomial
third_order – full second order polynomial model and all interactions of degree 3
cubic – third order polynomial model with cubic terms
The intercept is always included as a primary term.
Possible elements of the potential.model vector argument:
linear_interactions – linear interactions among the factors (default for MSE criteria)
quadratic_terms – quadratic terms for all the factors
third_order_terms – all interactions of degree 3: linear-by-linear-by-linear and quadratic-by-linear terms
cubic_terms – cubic terms for all the factors
fourth_order_terms – all interactions of degree 4, similar to third_order_terms
primary.terms and potential.terms arguments are used to specify the fitted model and the potential terms explicitly – up to the total power of 4.
Single factor powers, are coded as a string starting with with "x" and followed by
the index of the factor and the power: "x32".
For example, is coded as "x32"; "x22" stands for , and "x4" stands for the linear term .
Interaction of factors' powers are coded by merging the individual factors'
powers, ordered by the factors' indexes. For example, is coded as "x1x22",
– as "x12x3x4".
Object of class mood specifying parameters of the experiment, compound criterion of choice, and primary and potential model terms.
K Number of factors.
Klev Number of levels of each factor, if all factors have the same number of levels.
Levels List of length K of the vectors containing values of the factors.
Nruns Number of runs of the experiment.
criterion.choice Compound criterion to be used for the optimal design search or evaluation.
Nstarts The number of randomly generated start designs of the search algorithm.
Biter Number of samples for evaluating the MSE determinant-based component criterion.
tau2 The variance scaling parameter for the prior distribution of the potential terms.
tau The square root of tau2
Cubic Whether the experimental region is cubic (TRUE) or spherical (FALSE).
MC Indicator of the multiple comparison (Bonferroni) correction for trace-based criteria.
prob.DP, prob.LP, prob.LoF, prob.LoFL Confidence levels for the DP-, LP-, lack of fit determinant- and trace-based criteria.
alpha.DP, alpha.LP, alpha.LoF, alpha.LoFL Significance levels for the DP-, LP-, lack of fit determinant- and trace-based criteria.
orth Whether the candidate sets are orthonormalised (TRUE) or not (FALSE).
Z0 Z0 matrix.
W Weight matrix for Ls criterion.
primary.terms Fitted (primary) model terms.
potential.terms Potential terms.
P The number of terms in the fitted model (including intercept).
Q The number of potential terms.
kappa.Ds, kappa.DP, kappa.L, kappa.LP, kappa.LoF, kappa.bias, kappa.mse Compound criterion weights.
warning.msg Warning messages.
Gilmour SG, Trinca LA (2012).
“Optimum Design of Experiments for Statistical Inference (with discussion).”
Journal of the Royal Statistical Society C, 61, 345-401.
Goos P, Kobilinsky A, O'Brien TE, Vandebroek M (2005).
“Model-Robust and Model-Sensitive Designs.”
Computational Statistics and Data Analysis, 49, 201-216.
example1 <- mood(K = 5, Levels = 3, Nruns = 40, criterion.choice = "GDP", kappa = list(kappa.Ds = 1./3, kappa.DP = 1./3, kappa.LoF = 1./3), control = list(Nstarts = 50, tau2 = 0.1), model_terms = list(primary.model = "second_order", potential.terms = c("x12x2", "x22x3", "x32x4", "x42x5"))) example1 example2 <- mood(K = 3, Nruns = 12, Levels = list(1:3, 1:2, 1:2), criterion.choice = "MSE.L", kappa = list(kappa.LP = 1./2, kappa.LoF = 1./4, kappa.mse = 1./4), control = list(Nstarts = 50, tau2 = 1), model_terms = list(primary.terms = "first_order", potential.terms = c("x12", "x12x2", "x12x3"))) example2example1 <- mood(K = 5, Levels = 3, Nruns = 40, criterion.choice = "GDP", kappa = list(kappa.Ds = 1./3, kappa.DP = 1./3, kappa.LoF = 1./3), control = list(Nstarts = 50, tau2 = 0.1), model_terms = list(primary.model = "second_order", potential.terms = c("x12x2", "x22x3", "x32x4", "x42x5"))) example1 example2 <- mood(K = 3, Nruns = 12, Levels = list(1:3, 1:2, 1:2), criterion.choice = "MSE.L", kappa = list(kappa.LP = 1./2, kappa.LoF = 1./4, kappa.mse = 1./4), control = list(Nstarts = 50, tau2 = 1), model_terms = list(primary.terms = "first_order", potential.terms = c("x12", "x12x2", "x12x3"))) example2
Performing search for a (nearly) optimum factorial design, optimising with respect to a specified compound criterion.
Search( mood.object, algorithm = c("ptex", "coordex"), parallel = FALSE, verbose = TRUE )Search( mood.object, algorithm = c("ptex", "coordex"), parallel = FALSE, verbose = TRUE )
mood.object |
Object of class |
algorithm |
Parameter specifying the search algorithm. If |
parallel |
If |
verbose |
If |
Search takes the mood object as an input with all the parameters of the experiment. Runs a point-exchange or a coordinate-exchange algorithm, returns design and model matrices, computation time and criteria values.
See (Koutra et al. 2024) for examples of using parallel = TRUE.
Updated object of class mood containing the outputs generated by the search:
X.design Design matrix.
df The number of pure error degrees of freedom.
X1 Primary model matrix for the found (nearly-) optimum design.
X2 Model matrix of potential terms for the found (nearly-) optimum design.
compound.value The compound criterion value of the (nearly-) optimum design.
criteria.values Component criteria values of the (nearly-) optimum design.
path The "path" of compound criterion values of the optimum designs obtained after for each random start.
time Computation time.
algorithm Point exchange or coordinate exchange used to find the design?
parallel Were different runs of the algorithm performed across different CPU cores (TRUE/FALSE)
Koutra V, Egorova O, Gilmour SG, Trinca LA (2024). “MOODE: An R Package for Multi-Objective Optimal Design of Experiments.” arXiv:2412.17158, https://arxiv.org/abs/2412.17158.
example1 <- mood(K = 2, Levels = 3, Nruns = 10, criterion.choice = "GDP", kappa = list(kappa.Ds = 1./3, kappa.DP = 1./3, kappa.LoF = 1./3), control = list(tau2 = 0.1), model_terms = list(primary.model = "first_order", potential.terms = c("x12", "x22", "x1x2"))) # Using point exchange Search_point <- Search(example1, algorithm = 'ptex') Search_point # Using coordinate exchange (the default for K>4) Search_coord <- Search(example1) Search_coordexample1 <- mood(K = 2, Levels = 3, Nruns = 10, criterion.choice = "GDP", kappa = list(kappa.Ds = 1./3, kappa.DP = 1./3, kappa.LoF = 1./3), control = list(tau2 = 0.1), model_terms = list(primary.model = "first_order", potential.terms = c("x12", "x22", "x1x2"))) # Using point exchange Search_point <- Search(example1, algorithm = 'ptex') Search_point # Using coordinate exchange (the default for K>4) Search_coord <- Search(example1) Search_coord