#############################################################
### Analyses for JCGS paper on Bayes Marginal Model Plots ###
#############################################################

naph <- read.table("data/naphthalene.dat", header = T, row.names = NULL)
naph <- cbind(naph, greg = .397*naph[,1] + .445*naph[,3] + .802*naph[,2])
attach(naph)

fit <- lm(Yn ~ AN + Btemp + Ctime)

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

library(modreg)

x11() # open a traditional graphics device if needed (UNIX)
windows() # open a traditional graphics device if needed (MS Windows)

bmmp(fit, samp = samp, incvar = F, bdm = F)
bmmp(fit, samp = samp, sparams = 4, incvar = F, bdm = F)

fit <- lm(Yn ~ (AN + Btemp + Ctime)^2 + I(AN^2) + I(Btemp^2) + I(Ctime^2))

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

bmmp(fit, samp = samp, sparams = 4, incvar = F, bdm = F)

fit <- lm(log(Yn) ~ greg)

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

bmmp(fit, samp = samp, default = F, h = greg, sparams = 4, incvar = F, bdm = F)

#fit <- gam(Yn ~ s(greg, df = 4))

#gamsamp <- gibbs.gam(fit, nwarm = 300, nkeep = nsamp, var.comp = F)
#samp <- postsamp.gam(gamsamp)

#bmmp(fit, samp = samp, default = F, h = greg, sparams = 4, incvar = F,
#     bdm = F)

detach(naph)
rm(naph, fit, nsamp, samp)