########
# bmmp #
########

naph <- read.table("data/naphthalene.dat", header=T, 
                   row.names=NULL)
attach(naph)
fit <- lm(Yn ~ (AN+Btemp+Ctime)^2 + I(AN^2) + I(Btemp^2) +
          I(Ctime^2))

nsamp <- 100
samp <- postsamp.lm(fit, nsamp = nsamp)

trellis.device() # trellis graphics device (S-PLUS only)
motif() # traditional graphics device (S-PLUS UNIX/Linux)
graphsheet() # traditional graphics device (S-PLUS Windows)
x11() # traditional graphics device (R UNIX/Linux)
windows() # traditional graphics device (R Windows)

# if you are using traditional graphics in S-PLUS, include 
# the argument usetrellis = F in the commands below; in R 
# usetrellis is F by default

# if you are using R, attach library(modreg)
library(modreg) # R only

bmmp(fit, samp = samp)

bmmp(fit, samp = samp, h = AN)
bmmp(fit, samp = samp, h = cbind(AN, Btemp))
bmmp(fit, samp = samp, h = AN, default = F)
bmmp(fit, samp = samp, default = F, random = 2)

bmmp(fit, samp = samp, sptry = c(2, 3, 4, 5))
bmmp(fit, samp = samp, spartry = c(.04, .01, .003, .001)) 
     # S-PLUS
bmmp(fit, samp = samp, spartry = c(1.3, 1.2, 1.1, 1)) # R
bmmp(fit, samp = samp, h = AN, sparams = c(3, 3))
bmmp(fit, samp = samp, h = AN, spar = c(.04, .06)) # S-PLUS
bmmp(fit, samp = samp, h = AN, spar = c(1.3, 1.1)) # R
bmmp(fit, samp = samp, spdefault = T)

bmmp(fit, samp = samp, sparams = 3, incvar = F)

bmmp(fit, samp = samp, sparams = 3, mult = .2)

bmmp(fit, samp = samp, h = AN, sparams = c(3, 3), 
     included = names(AN[AN<2])) # S-PLUS
bmmp(fit, samp = samp, h = AN, sparams = c(3, 3),
     included = as.character(seq(along=AN))[AN<2]) # R

bmmp(fit, samp = samp, usetrellis = F) # traditional graphics

bmmp(fit, samp = samp, smoother = "loess1")
bmmp(fit, samp = samp, smoother = "loess1", 
     h = cbind(AN, Btemp, Ctime), 
     sparams = c(.7, .8, .8, .8))
bmmp(fit, samp = samp, smoother = "loess2")
bmmp(fit, samp = samp, smoother = "loess2", 
     h = cbind(AN, Btemp, Ctime), 
     sparams = c(1, 1, 1, .9))
bmmp(fit, samp = samp, smoother = "kernel")
bmmp(fit, samp = samp, smoother = "kernel", 
     h = cbind(AN, Btemp, Ctime), 
     sparams = c(.2, .4, .3, .3))
bmmp(fit, samp = samp, smoother = "supsmu", incvar = F)
bmmp(fit, samp = samp, smoother = "supsmu", incvar = F,
     h = cbind(AN, Btemp, Ctime), 
     sparams = c(.5, .6, .6, .4))

x <- .397*AN + .445*Ctime + .802*Btemp
fit <- gam(Yn ~ s(x, df = 4)) # S-PLUS UNIX/Linux only

# posterior sampling requires Hastie and Tibshirani's 
# gibbs.gam() function
ggnaph <- gibbs.gam(fit, nwarm = 200, nkeep = nsamp, 
                    var.comp = F)

# alternatively download ggnaph from bmmp web-site
source("data/ggnaph") # S-PLUS UNIX/Linux only

samp <- postsamp.gam(ggnaph)
# ggnaph is the object output from gibbs.gam()

bmmp(fit, samp = samp, h = AN)
bmmp(fit, samp = samp, default = F, random = 1) 
     # adds no information when p = 1
bmmp(fit, samp = samp, h = cbind(AN, Btemp, Ctime), 
     sparams = c(4, 3, 3, 3))
bmmp(fit, samp = samp, h = AN, sparams = c(4, 3), 
     included = names(AN[AN<2]))
bmmp(fit, samp = samp, smoother = "loess1")
bmmp(fit, samp = samp, smoother = "loess1", 
     h = cbind(AN, Btemp, Ctime), 
     sparams = c(.8, .8, .8, .8))
bmmp(fit, samp = samp, smoother = "loess2")
bmmp(fit, samp = samp, smoother = "loess2", 
     h = cbind(AN, Btemp, Ctime), 
     sparams = c(1, 1, 1, 1))
bmmp(fit, samp = samp, smoother = "kernel")
bmmp(fit, samp = samp, smoother = "kernel", 
     h = cbind(AN, Btemp, Ctime), 
     sparams = c(.2, .4, .3, .3))
bmmp(fit, samp = samp, smoother = "supsmu", incvar = F)
bmmp(fit, samp = samp, smoother = "supsmu", incvar = F,
     h = cbind(AN, Btemp, Ctime), 
     sparams = c(.6, .6, .6, .6))

detach("naph")
rm(fit, naph, nsamp, samp, x, ggnaph) # S-PLUS UNIX/Linux
rm(fit, naph, nsamp, samp) # S-PLUS Windows, R

wbcd <- read.table("data/WBCD.dat", header = T, 
                   row.names = NULL)
attach(wbcd)
fit <- glm(Class1 ~ Adhes + BNucl + Chrom + NNucl + Thick, 
           family = binomial)

# posterior sampling requires BUGS and BOA

# alternatively download boawbcd from bmmp web-site
data.restore("data/boawbcd") # S-PLUS
load("data/boawbcd") # R

nsamp <- 100
samp <- postsamp.binlog(fit, boaobject = boawbcd, 
                        nsamp = nsamp)
# boawbcd is the object output from BOA

bmmp(fit, samp = samp, h = Mitos, sptry = c(4, 8, 12, 16), 
     incvar = F)
bmmp(fit, samp = samp, default = F, random = 2, incvar = F)
bmmp(fit, samp = samp, h = cbind(Mitos, Adhes, BNucl), 
     sparams = c(12, 4, 4, 4), incvar = F)
bmmp(fit, samp = samp, h = Mitos, sparams = c(12, 4), 
     incvar = F, included = names(Mitos[Mitos<10])) # S-PLUS
bmmp(fit, samp = samp, h = Mitos, sparams = c(12, 4), 
     included = as.character(seq(along=Mitos))[Mitos<10],
     incvar = F) # R
bmmp(fit, samp = samp, smoother = "loess1", incvar = F, 
     sptry = c(.4, .3, .2, .1))
bmmp(fit, samp = samp, smoother = "loess1", incvar = F,
     h = cbind(Mitos, Adhes, BNucl), 
     sparams = c(.1, .95, .7, .7))
bmmp(fit, samp = samp, smoother = "loess2", incvar = F, 
     sptry = c(.4, .3, .2, .1))
bmmp(fit, samp = samp, smoother = "loess2", incvar = F,
     h = cbind(Mitos, Adhes, BNucl), 
     sparams = c(.3, .9, .85, .9))
bmmp(fit, samp = samp, smoother = "kernel", incvar = F, 
     sptry = c(.4, .3, .2, .1))
bmmp(fit, samp = samp, smoother = "kernel", incvar = F,
     h = cbind(Mitos, Adhes, BNucl), 
     sparams = c(.1, .3, .3, .3))
bmmp(fit, samp = samp, smoother = "supsmu", incvar = F,
     sptry = c(.4, .3, .2, .1))
bmmp(fit, samp = samp, smoother = "supsmu", incvar = F,
     h = cbind(Mitos, Adhes, BNucl), 
     sparams = c(.1, .1, .1, .1))

detach("wbcd")
rm(fit, wbcd, boawbcd, nsamp, samp)

#######
# mmp #
#######

naph <- read.table("data/naphthalene.dat", header=T, 
                   row.names=NULL)
attach(naph)
fit <- lm(Yn ~ (AN+Btemp+Ctime)^2 + I(AN^2) + I(Btemp^2) +
          I(Ctime^2))

# if you are using traditional graphics in S-PLUS, include 
# the argument usetrellis = F in the commands below; in R 
# usetrellis is F by default

# if you are using R, attach library(modreg)
library(modreg) # R only

mmp(fit)

mmp(fit, h = AN)
mmp(fit, h = cbind(AN, Btemp))
mmp(fit, h = AN, default = F)
mmp(fit, default = F, random = 2)

mmp(fit, sptry = c(2, 3, 4, 5))
mmp(fit, spartry = c(.04, .01, .003, .001)) # S-PLUS
mmp(fit, spartry = c(1.3, 1.2, 1.1, 1)) # R
mmp(fit, h = AN, sparams = c(3, 3))
mmp(fit, h = AN, spar = c(.04, .06)) # S-PLUS
mmp(fit, h = AN, spar = c(1.3, 1.1)) # R
mmp(fit, spdefault = T)

mmp(fit, sparams = 3, adjvar = T)
mmp(fit, sparams = 3, incvar = F)

mmp(fit, sparams = 3, mult = .1)

mmp(fit, h = AN, sparams = c(3, 3), 
    included = names(AN[AN<2])) # S-PLUS
mmp(fit, h = AN, sparams = c(3, 3),
    included = as.character(seq(along=AN))[AN<2]) # R

mmp(fit, usetrellis = F) # traditional graphics

mmp(fit, smoother = "loess1")
mmp(fit, smoother = "loess1", h = cbind(AN, Btemp, Ctime),
    sparams = c(.7, .8, .8, .8))
mmp(fit, smoother = "loess2")
mmp(fit, smoother = "loess2", h = cbind(AN, Btemp, Ctime),
    sparams = c(1, 1, 1, .9))
mmp(fit, smoother = "kernel")
mmp(fit, smoother = "kernel", h = cbind(AN, Btemp, Ctime),
    sparams = c(.2, .4, .3, .3))
mmp(fit, smoother = "supsmu", incvar = F)
mmp(fit, smoother = "supsmu", h = cbind(AN, Btemp, Ctime),
    sparams = c(.5, .6, .6, .4), incvar = F)

x <- .397*AN + .445*Ctime + .802*Btemp
fit <- gam(Yn ~ s(x, df = 4)) # S-PLUS only

mmp(fit, h = AN)
mmp(fit, h = cbind(AN, Btemp, Ctime), 
    sparams = c(4, 3, 3, 3))
mmp(fit, h = AN, sparams = c(4, 3), 
    included = names(AN[AN<2]))
mmp(fit, smoother = "loess1")
mmp(fit, smoother = "loess1", h = cbind(AN, Btemp, Ctime),
    sparams = c(.8, .8, .8, .8))
mmp(fit, smoother = "loess2")
mmp(fit, smoother = "loess2", h = cbind(AN, Btemp, Ctime),
    sparams = c(1, 1, 1, 1))
mmp(fit, smoother = "kernel")
mmp(fit, smoother = "kernel", h = cbind(AN, Btemp, Ctime),
    sparams = c(.2, .4, .3, .3))
mmp(fit, smoother = "supsmu", incvar = F)
mmp(fit, smoother = "supsmu", h = cbind(AN, Btemp, Ctime),
    sparams = c(.6, .6, .6, .6), incvar = F)

detach("naph")
rm(fit, naph, x) # S-PLUS
rm(fit, naph) # R

wbcd <- read.table("data/WBCD.dat", header = T, 
                   row.names = NULL)
attach(wbcd)
fit <- glm(Class1 ~ Adhes + BNucl + Chrom + NNucl + Thick, 
           family = binomial)

mmp(fit, h = Mitos, sptry = c(4, 8, 12, 16), incvar = F)
mmp(fit, default = F, random = 2, incvar = F)
mmp(fit, h = cbind(Mitos, Adhes, BNucl), 
    sparams = c(12, 4, 4, 4), incvar = F)
mmp(fit, h = Mitos, sparams = c(12, 4), incvar = F,
    included = names(Mitos[Mitos<10])) # S-PLUS
mmp(fit, h = Mitos, sparams = c(12, 4), incvar = F,
    included = as.character(seq(along=Mitos))[Mitos<10]) # R
mmp(fit, smoother = "loess1", sptry = c(.4, .3, .2, .1), 
    incvar = F)
mmp(fit, smoother = "loess1", h = cbind(Mitos, Adhes, BNucl),
    sparams = c(.1, .95, .7, .7), incvar = F)
mmp(fit, smoother = "loess2", sptry = c(.4, .3, .2, .1), 
    incvar = F)
mmp(fit, smoother = "loess2", h = cbind(Mitos, Adhes, BNucl),
    sparams = c(.3, .9, .85, .9), incvar = F)
mmp(fit, smoother = "kernel", sptry = c(.4, .3, .2, .1), 
    incvar = F)
mmp(fit, smoother = "kernel", h = cbind(Mitos, Adhes, BNucl),
    sparams = c(.1, .3, .3, .3), incvar = F)
mmp(fit, smoother = "supsmu", sptry = c(.4, .3, .2, .1), 
    incvar = F) 
mmp(fit, smoother = "supsmu", h = cbind(Mitos, Adhes, BNucl),
    sparams = c(.1, .1, .1, .1), incvar = F)

detach("wbcd")
rm(fit, wbcd)