BsProb.design {FrF2}

Bayesian posterior probabilities from Box and Meyer method

Description

The function calculates Bayesian posterior probabilities according to Box and Meyer (1993) for screening experiments with 2-level factors. The function is modified from function BsProb in packge BsMD with the purpose of providing usage comfort for class design objects.

Usage

BsProb.design(design, mFac = NULL, response=NULL, select=NULL, mInt = 2, p = 0.25, g = 2,
    ng = 1, nMod = 10)

Arguments

design

an experimental design of class design with the type element of the design.info attribute containing “FrF2” or “pb” and at least one response variable

response

NULL or a character string that specifies response variable to be used, must be an element of response.names(obj); if NULL, the first response from response.names(obj) is used

mFac

integer. Maximum number of factors included in the models. The default is the number of factors in the design.

select

vector with position numbers of the factors to be included;
default: all factors.

mInt

integer <= 3. Maximum order of interactions considered in the models. This can strongly impact the result.

p

numeric. Prior probability assigned to active factors. This can strongly impact the result.

g

numeric vector. Variance inflation factor(s) gamma associated to active and interaction factors; see "Details" section

ng

integer <=20. Number of different variance inflation factors (g) used in calculations.

nMod

integer <=100. Number of models to keep with the highest posterior probability.

Details

Factor and model posterior probabilities are computed by the Box and Meyer (1993) Bayesian procedure. The design factors - or a selection of these given by column numbers in select - are considered together with the specified response or the first response of the design. The function has been adapted from function BsProb in package BsMD, and a vignette in that package (../../BsMD/doc/BsMD.pdf) explains details of the usage regarding the parameters.

If g, the variance inflation factor (VIF) gamma, is a vector of length 1, the same VIF is used for factor main effects and interactions. If the length of g is 2 and ng is 1, g[1] is used for factor main effects and g[2] for the interaction effects. If ng greater than 1, then ng values of VIFs between g[1] and g[2] are used for calculations with the same gamma value for main effects and interactions. The function calls the FORTRAN subroutine bm and captures summary results. The complete output of the FORTRAN code is save in the BsPrint.out file in the working directory. The output is a list of class BsProb for which print, plot and summary methods are available from package BsMD.

Value

cf. documentation of function BsProb

Note

This method relies on the availability of package BsMD.

Author(s)

Daniel Meyer, ported to R by Ernesto Barrios, port adapted to designs by Ulrike Groemping.

References

Barrios, E. (2013). Using the BsMD Package for Bayesian Screening and Model Discrimination. Vignette. ../../BsMD/doc/BsMD.pdf.

Box, G. E. P and R. D. Meyer (1986). An Analysis for Unreplicated Fractional Factorials. Technometrics 28, 11-18.

Box, G. E. P and R. D. Meyer (1993). Finding the Active Factors in Fractionated Screening Experiments. Journal of Quality Technology 25, 94-105.

See Also

plot.BsProb, print.BsProb, summary.BsProb, BsMD

Examples

   ### there are several success stories and recommendations for this method
   ### in the simulated example here (not fabricated, 
   ###         it was the first one that came to my mind), 
   ### the method goes wrong, at least when using mInt=2 (the default, because 
   ###         Daniel plots work quite well for pure main effects models):
   ### active factors are A to E (perhaps too many for the method to work),
   ### the method identifies F, J, and L with highest probability 
   ### (but is quite undecided)
   plan <- pb(12)
   dn <- desnum(plan)
   set.seed(8655)
   y <- dn%*%c(2,2,2,2,3,0,0,0,0,0,0) + dn[,1]*dn[,3]*2 - dn[,5]*dn[,4] + rnorm(12)/10
   plan.r <- add.response(plan, response=y)
   if (requireNamespace("BsMD", quiet=TRUE)){
   plot(bpmInt2 <- BsProb.design(plan.r), code=FALSE)
   plot(bpmInt1 <- BsProb.design(plan.r, mInt=1), code=FALSE) ## much better!
   summary(bpmInt2)
   summary(bpmInt1)
   }
   ### For comparison: A Daniel plot does not show any significant effects according 
   ### to Lenths method, but makes the right effects stick out
   DanielPlot(plan.r, half=TRUE, alpha=1)