View source: R/BsProb.design.R
BsProb.design | R Documentation |
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.
BsProb.design(design, mFac = NULL, response=NULL, select=NULL, mInt = 2, p = 0.25, g = 2,
ng = 1, nMod = 10)
design |
an experimental design of class |
response |
NULL or a character string that specifies response variable to be used,
must be an element of |
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; |
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. |
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.
cf. documentation of function BsProb
This method relies on the availability of package BsMD.
Daniel Meyer, ported to R by Ernesto Barrios, port adapted to designs by Ulrike Groemping.
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.
plot.BsProb
, print.BsProb
,
summary.BsProb
, BsMD
### 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.