# Best Model Discrimination (MD) Follow-Up Experiments Best follow-up experiments based on the MD criterion are suggested to discriminate between competing models. ```r MD(X, y, nFac, nBlk = 0, mInt = 3, g = 2, nMod, p, s2, nf, facs, nFDes = 4, Xcand, mIter = 20, nStart = 5, startDes = NULL, top = 20, eps = 1e-05) ``` ## Arguments - `X`: matrix. Design matrix of the initial experiment. - `y`: vector. Response vector of the initial experiment. - `nFac`: integer. Number of factors in the initial experiment. - `nBlk`: integer >=1. The number of blocking factors in the initial experiment. They are accommodated in the first columns of matrix `X`. - `mInt`: integer. Maximum order of the interactions in the models. - `g`: vector. Variance inflation factor for main effects (`g[1]`) and interactions effects (`g[2]`). If vector length is 1 the same inflation factor is used for main and interactions effects. - `nMod`: integer. Number of competing models. - `p`: vector. Posterior probabilities of the competing models. - `s2`: vector. Competing model variances. - `nf`: vector. Factors considered in each of the models. - `facs`: matrix. Matrix [`nMod x max(nf)`] of factor numbers in the design matrix. - `nFDes`: integer. Number of runs to consider in the follow-up experiment. - `Xcand`: matrix. Candidate runs to be chosen for the follow-up design. - `mIter`: integer. If 0, then user-entered designs `startDes` are evaluated, otherwise the maximum number of iterations for each Wynn search. - `nStart`: integer. Number of starting designs. - `startDes`: matrix. Matrix `[nStart x nFDes]`. Each row has the row numbers of the user-supplied starting design. - `top`: integer. Highest MD follow-up designs recorded. - `eps`: numeric. A small number (1e-5 by default) used for computations. ## Details The MD criterion, proposed by Meyer, Steinberg and Box is used to discriminate among competing models. Random starting runs chosen from `Xcand` are used for the Wynn search of best MD follow-up designs. `nStart` starting points are tried in the search limited to `mIter` iterations. If `mIter=0` then `startDes` user-provided designs are used. Posterior probabilities and variances of the competing models are obtained from `BsProb`. The function calls the FORTRAN subroutine md and captures summary results. The complete output of the FORTRAN code is save in the MDPrint.out file in the working directory. ## Returns A list with all input and output parameters of the FORTRAN subroutine `MD`. Most of the variable names kept to match FORTRAN code. - **NSTART**: Number of starting designs. - **NRUNS**: Number of runs used in follow-up designs. - **ITMAX**: Maximum number of iterations for each Wynn search. - **INITDES**: Number of starting points. - **NO**: Numbers of runs already completed before follow-up. - **IND**: Indicator; 0 indicates the user supplied starting designs. - **X**: Matrix for initial data (`nrow(X)=N0`; `ncol(X)=COLS+BL`). - **Y**: Response values from initial experiment (`length(Y)=N0`). - **GAMMA**: Variance inflation factor. - **GAM2**: If `IND=1`, `GAM2` was used for interaction factors. - **BL**: Number of blocks (>=1) accommodated in first columns of `X` and `Xcand`. - **COLS**: Number of factors. - **N**: Number of candidate runs. - **Xcand**: Matrix of candidate runs. (`nrow(Xcand)=N`, `ncol(Xcand)=ncol(X)`). - **NM**: Number of models considered. - **P**: Models posterior probability. - **SIGMA2**: Models variances. - **NF**: Number of factors per model. - **MNF**: Maximum number of factor in models. (`MNF=max(NF)`). - **JFAC**: Matrix with the factor numbers for each of the models. - **CUT**: Maximum interaction order considered. - **MBEST**: If `INITDES=0`, the first row of the `MBEST[1,]` matrix has the first user-supplied starting design. The last row the `NSTART`-th user-supplied starting design. - **NTOP**: Number of the top best designs. - **TOPD**: The D value for the best `NTOP` designs. - **TOPDES**: Top `NTOP` design factors. - **ESP**: "Small number" provided to the md FORTRAN subroutine. 1e-5 by default. - **flag**: Indicator = 1, if the md subroutine finished properly, -1 otherwise. ## Note The function is a wrapper to call the FORTAN subroutine md , modification of Daniel Meyer's original program, MD.f , part of the `mdopt` bundle for Bayesian model discrimination of multifactor experiments. ## References Meyer, R. D., Steinberg, D. M. and Box, G. E. P. (1996). "Follow-Up Designs to Resolve Confounding in Multifactor Experiments (with discussion)". **Technometrics**, Vol. 38, No. 4, pp. 303--332. Box, G. E. P and R. D. Meyer (1993). "Finding the Active Factors in Fractionated Screening Experiments". **Journal of Quality Technology.** Vol. 25. No. 2. pp. 94--105. ## Author(s) R. Daniel Meyer. Adapted for by Ernesto Barrios. ## See Also `print.MD`, `BsProb` ## Examples ```r ### Injection Molding Experiment. Meyer et al. 1996, example 2. library(BsMD) data(BM93.e3.data,package="BsMD") X <- as.matrix(BM93.e3.data[1:16,c(1,2,4,6,9)]) y <- BM93.e3.data[1:16,10] p <- c(0.2356,0.2356,0.2356,0.2356,0.0566) s2 <- c(0.5815,0.5815,0.5815,0.5815,0.4412) nf <- c(3,3,3,3,4) facs <- matrix(c(2,1,1,1,1,3,3,2,2,2,4,4,3,4,3,0,0,0,0,4),nrow=5, dimnames=list(1:5,c("f1","f2","f3","f4"))) nFDes <- 4 Xcand <- matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, -1,-1,-1,-1,1,1,1,1,-1,-1,-1,-1,1,1,1,1, -1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1, -1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1, -1,1,1,-1,1,-1,-1,1,1,-1,-1,1,-1,1,1,-1), nrow=16,dimnames=list(1:16,c("blk","f1","f2","f3","f4")) ) injectionMolding.MD <- MD(X = X, y = y, nFac = 4, nBlk = 1, mInt = 3, g = 2, nMod = 5, p = p, s2 = s2, nf = nf, facs = facs, nFDes = 4, Xcand = Xcand, mIter = 20, nStart = 25, top = 10) summary(injectionMolding.MD) ### Reactor Experiment. Meyer et al. 1996, example 3. par(mfrow=c(1,2),pty="s") data(Reactor.data,package="BsMD") # Posterior probabilities based on first 8 runs X <- as.matrix(cbind(blk = rep(-1,8), Reactor.data[c(25,2,19,12,13,22,7,32), 1:5])) y <- Reactor.data[c(25,2,19,12,13,22,7,32), 6] reactor8.BsProb <- BsProb(X = X, y = y, blk = 1, mFac = 5, mInt = 3, p =0.25, g =0.40, ng = 1, nMod = 32) plot(reactor8.BsProb,prt=TRUE,,main="(8 runs)") # MD optimal 4-run design p <- reactor8.BsProb$ptop s2 <- reactor8.BsProb$sigtop nf <- reactor8.BsProb$nftop facs <- reactor8.BsProb$jtop nFDes <- 4 Xcand <- as.matrix(cbind(blk = rep(+1,32), Reactor.data[,1:5])) reactor.MD <- MD(X = X, y = y, nFac = 5, nBlk = 1, mInt = 3, g =0.40, nMod = 32, p = p,s2 = s2, nf = nf, facs = facs, nFDes = 4, Xcand = Xcand, mIter = 20, nStart = 25, top = 5) summary(reactor.MD) # Posterior probabilities based on all 12 runs X <- rbind(X, Xcand[c(4,10,11,26), ]) y <- c(y, Reactor.data[c(4,10,11,26),6]) reactor12.BsProb <- BsProb(X = X, y = y, blk = 1, mFac = 5, mInt = 3, p = 0.25, g =1.20,ng = 1, nMod = 5) plot(reactor12.BsProb,prt=TRUE,main="(12 runs)") ```

MD function