Experiment Data

The experiment data arguments include what is commonly referred to in regression as the “design matrix”, the response, and the number of observations that have been seen in the design matrix so far.

The design matrix is a \(n\) by \(p\) object (where \(p\) is the number of model parameters), which can have the structure of either a dataframe or a matrix. The response contains all of the observations seen to date, which may be in the form of a vector, a dataframe, or a matrix. The response may either have the number of observations that have been seen, or it may have \(n\) observations where the observations that have been seen have “NA” entries.

For this example, we have data loaded in an file named exData. In this data set, the first column is the response variable, and the remaining 13 columns are the entries of design matrix. The purpose of this simulated experiment is to determine if the mean number of miles traveled on one charge is greater than or equal to 400 miles (for a simulated electric semi-truck). We will refer to the mean number of miles traveled as the parameter of interest (\(\phi\)), and the threshold value the parameter of interest most obtain as \(\phi_0\) (coded as phi.0).

The experiment selected was a \(2^4\) full factorial with five replicates that included main effects and two-way interactions (excluding the uncontrollable factor, \(\omega_k\)) for the experimental design, resulting in 80 test events. This design has a power of 80%, with an 80% confidence level, to detect a difference of 50 miles with a standard deviation of 100 miles.

This example assumes that the first 75 observations of the 80 observations have been seen, and seeks to use predictive probability to determine if the last five observations need to be seen to make a determination about the mean number of miles traveled, or if the experiment has provided enough information to stop testing early.

Of note, there are limitations to what this R program can support.

1. Experimental designs that can be used include fully randomized designs (e.g., factorial designs and optimal designs) and designs that block on replicates.

2. The priors on all model parameters are assumed to be normal, and the prior on the precision is assumed to be gamma. The user may select the the appropriate values in these priors; however, these priors were selected to ensure computational efficiency.

3. The data is assumed to follow an reference cell ANOVA data model, where the model parameters are mutually independent. The constraint selected is that the first level of each factor (or an interaction that contains a factor at the first level) are all equal to 0. This gives \(\eta\), the baseline parameter, as the first level of each factor. For example, the data model for this data set is:

\[y_{ijklmp|\mu_{ijklm}\sim N(\mu_{ijklm}, \frac{1}{\tau})}\]

where

\[\mu_{ijklm} = \eta + \alpha_i + \beta_j + \omega_k + \gamma_l + \delta+m + \alpha\beta_{(ij)} + \alpha\gamma_{(il)} + \alpha\delta_{(im)} + \beta\gamma_{(jl)} + \beta\delta_{(jm)} + \delta\gamma_{(lm)}\]

and

\[ i=1, 2\\ j=1, 2 \\ k=1, 2, 3 \text{ (uncontrolled factor)}) \\ l=1,2 \\ m=1,2 \]

# Design Matrix
X <- exData[,c(2:14)]
head(X)

##      V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14
## [1,]  1  0  0  1  0  1  1  0   0   0   0   0   1
## [2,]  1  1  1  0  0  0  0  1   0   0   0   0   0
## [3,]  1  0  0  0  0  0  1  0   0   0   0   0   0
## [4,]  1  0  1  1  0  0  0  0   0   0   0   0   0
## [5,]  1  1  1  0  0  1  1  1   1   1   1   1   1
## [6,]  1  1  1  0  0  0  0  1   0   0   0   0   0

# Response Column
Y <- as.matrix(exData[,1], ncol=1) # could also be a vector e.g., Y <- exData[,1]
# Observations Seen
n.seen <- 75

Priors

This program uses a reference cell ANOVA model that constrains the first level of each model parameter (or an interaction that contains the first level of a factor) to 0. Therefore, degenerate priors at 0 were selected for these constrained parameters and are not included in the analysis. The (weakly informative) priors for the unconstrained parameters are below:

\[ p(\eta) \sim N(400, 100^2) \\ p(\alpha) \sim N(50, 100^2) \\ p(\beta) \sim N(50, 50^2) \\ p(\omega_2) \sim N(-25, 50^2) \\ p(\omega_3) \sim N(-50, 50^2) \\ p(\gamma) \sim N(100, 100^2) \\ p(\delta) \sim N(100, 100^2) \\ p(\tau) ~ \sim \text{Gam}(0.0001, 0.0001) \]

The prior on all two way interactions is \(N(0, 100^2)\).

The set of unconstrained model parameters is captured by \(\beta\) (not to be confused with \(\beta_j\)). The mean of all the priors (excluding \(\tau\)) are assigned to beta.mean, which can either be a matrix (as below) or a vector. The precision of the prior is assigned to beta.precision, and may also be either a matrix (as below) or a vector. This results in a prior on \(\beta\) that is a multivariate normal distribution with a mean vector equal to the mean of the priors and a precision matrix with a diagonal of the prior precisions and off-diagonals of 0 (to represent independent factors).

The sampler will initiate the MCMC chain at the mean of the priors selected. For more specific details on the priors, refer to Sieck and Christensen (2021).

# Means for betas
beta.mean <- matrix(
  c(400, 50, 50, -25, -50, 100, 100, rep(0, 6)),
  ncol = 1
) 
# could also be vector: beta.mean <- c(400, 50, 50, -25, -50, 100, 100, rep(0, 6))

# Precisions for betas
beta.precision <- matrix(
  c(1/10000, 1/10000, 1/10000, 1/2500, 1/2500, 1/10000, 1/10000, rep(1/10000, 6)), 
  ncol = 1
) 
# could also be vector: beta.precision <- c(1/10000, 1/10000, 1/10000, 1/2500, 1/2500, 1/10000, 1/10000, rep(1/10000, 6))

# Precision for tau
shape <- 0.0001
rate <- 0.0001

Simulation Parameters

These parameters specify the amount of burn-in and number of posterior samples for the conditional draws (i.e., the nested sampler that draws observations of the model parameters, conditioned on the seen and a set of unseen data as defined by the outer sampler) and non-conditional draws (i.e., the outer sampler that draws a set of unseen data for the remaining observations). See Sieck and Christensen (2021) for more details, to include flow charts of the process.

# Conditional Draws
b.sim <- 20000
b.burnin <- 2000

# Non-Conditional Draws
n.sim <- 1000
y.burnin <- 100

Analysis Parameters

To determine if the parameter of interest exceeds the threshold value, the probability of this occurring can be calculated. If this probability is high enough, then the question can be evaluated as met (or not met if low enough). This threshold probability is theta.t, a probability which is used to express how much certainty is required before stating that the parameter of interest exceeds the threshold value.

As defined above, phi.0 is the threshold value that must be obtained by the parameter of interest for the question to be evaluated as met.

# Threshold Value
theta.t <- 0.8

# Metric Threshold Value
phi.0 <- 400

The prob object is used to support a Bayesian Mission Mean analysis. This analysis uses a Bayesian Mission Mean (BMM; Sieck and Christensen 2021). Each combination of factor levels creates a “mission set” for the system being evaluated that represent operational environments proposed for use. Using a BMM provides a method to obtain a summary statistic based on performance across operational environments without artificially reducing the operational variance. This is accomplished by considering the joint distribution of mission sets and mission means, and then marginalizing over mission sets to obtain a mixture distribution of mission means.

Mission sets are drawn based on the selected anticipated likelihood of encountering levels. Once the mission sets and the conditional posterior distribution on the model parameters are obtained from the nested sampler, they can be used to induce a distribution on our parameter of interest. This is accomplished by using the mission sets to mix the mission means from the posterior distribution. This induced distribution can then be used to calculate the direct probability that the parameter of interest is greater than the threshold value across all mission sets.

You can use the prob.creator function to construct the necessary prob matrix for the gibbs.sampler function. For prob.creator: num.factors is the number of factors in the model (e.g., \(\alpha_i\) with \(i=1,2\) is 1 factor; from above, this example has 5 factors). num.factor.levels contains the number of levels for each factor (e.g., \(\alpha_i\) with \(i=1,2\) has 2 levels), and may be a vector, a matrix, or a dataframe.

The likelihood of encountering are captured in the object likelihood.encountering, and can be a vector, dataframe, or matrix. This object, unlike the priors, must include the likelihood of encountering a level for all model parameters (constrained or unconstrained). As an example, below indicates that the likelihood of encountering \(\alpha_1\) is 1/2 and the likelihood of encountering \(\alpha_2\) is 1/2. Of note, the likelihood of encountering the levels must be in the same order as defined by the factor in the design matrix.

num.factors <- 5
num.factor.levels <- c(2, 2, 3, 2, 2)
likelihood.encountering <- c(1/2, 1/2, 4/9, 5/9, 1/3, 1/3, 1/3, 1/2, 1/2, 1/2, 1/2)

prob <- prob.creator(num.factors, num.factor.levels, likelihood.encountering, print.result = TRUE)

##       Factor Likelihood of Encountering
##  [1,]      1                  0.5000000
##  [2,]      1                  0.5000000
##  [3,]      2                  0.4444444
##  [4,]      2                  0.5555556
##  [5,]      3                  0.3333333
##  [6,]      3                  0.3333333
##  [7,]      3                  0.3333333
##  [8,]      4                  0.5000000
##  [9,]      4                  0.5000000
## [10,]      5                  0.5000000
## [11,]      5                  0.5000000

You may also write your own prob matrix, as shown below.

# Mission Probabilities (likelihood of encountering a factor level) 
prob <- matrix(
  c(1, 1, 2, 2, 3, 3, 3, 4, 4, 5, 5, 1/2, 1/2, 4/9, 5/9, 1/3, 1/3, 1/3, 1/2, 1/2, 1/2, 1/2), 
  ncol = 2, 
  dimnames = list(NULL, c("factor", "probability"))
)

Optional Arguments

By default, the model will utilize two-way interactions from all the main effects defined in the prob object. One option the user has is to remove a factor from consideration within the two-way interactions. For example, we will remove the third model parameter from consideration because it is an uncontrolled factor in the experiment. Note that factor.no.2way should be provided factor numbers to drop two-way interactions involving that factor.

factor.no.2way <- c(3)

By default, the function will assign the column names from the design matrix to the result posterior dataframe. Alternatively, users can specify their own column names for the result dataframe.

colnames.pick <- c(
  "eta", "alpha", "beta", "omega2", "omega3", 
  "theta", "gamma", "alphabeta", "alphatheta", 
  "alphagamma", "betatheta", "betagamma", "thetagamma", "tau"
)