examples.R
In hypergsplines: Bayesian model selection with penalised splines and hyper-g prior

### R code from vignette source 'examples.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: pima-data
###################################################
library(MASS)
pima <- rbind(Pima.tr, Pima.te)
pima$hasDiabetes <- as.numeric(pima$type == "Yes")
pima.nObs <- nrow(pima)


###################################################
### code chunk number 2: pima-setup
###################################################
library(hypergsplines)
modelData.pima <- with(pima,
                       glmModelData(y=hasDiabetes,
                                    X=
                                    cbind(npreg,
                                          glu,
                                          bp,
                                          skin,
                                          bmi,
                                          ped,
                                          age),
                                    splineType="cubic",
                                    nKnots=4L,                                    
                                    family=binomial))


###################################################
### code chunk number 3: pima-search-setup
###################################################
chainlength.pima <- 100
computation.pima <- getComputation(useOpenMP=FALSE,
                                   higherOrderCorrection=FALSE)


###################################################
### code chunk number 4: pima-search
###################################################
set.seed(93)
time.pima <-
    system.time(models.pima <-
                stochSearch(modelData=modelData.pima,
                            modelPrior="dependent",
                            chainlength=chainlength.pima,
                            nModels=chainlength.pima,
                            computation=computation.pima))


###################################################
### code chunk number 5: pima-top-models
###################################################
head(models.pima$models)
map.pima <- models.pima$models[1, 1:7]


###################################################
### code chunk number 6: pima-incprobs
###################################################
round(models.pima$inclusionProbs,2)


###################################################
### code chunk number 7: pima-sampling
###################################################
mcmc.pima <- getMcmc(samples=500L,
                     burnin=100L,
                     step=1L,
                     nIwlsIterations=2L)

set.seed(634)
map.samples.pima <- glmGetSamples(config=map.pima,
                                  modelData=modelData.pima,
                                  mcmc=mcmc.pima,
                                  computation=computation.pima)


###################################################
### code chunk number 8: pima-samples-structure
###################################################
str(map.samples.pima)


###################################################
### code chunk number 9: pima-plot-ex
###################################################
plotCurveEstimate(covName="age",
                  samples=map.samples.pima$samples,
                  modelData=modelData.pima)


###################################################
### code chunk number 10: pima-postprocess
###################################################
optim.map.pima <- postOptimize(modelData=modelData.pima,
                               modelConfig=map.pima,
                               computation=computation.pima)
optim.map.pima


###################################################
### code chunk number 11: pima-fit-samples
###################################################
fit.samples.pima <- getFitSamples(X=modelData.pima$origX[1:10,],
                                  samples=map.samples.pima$samples,
                                  modelData=modelData.pima)
obs.samples.pima <- plogis(fit.samples.pima)


###################################################
### code chunk number 12: pima-postpred-means
###################################################
rowMeans(obs.samples.pima)


###################################################
### code chunk number 13: pima-observations
###################################################
modelData.pima$Y[1:10]


###################################################
### code chunk number 14: pima-model-averaging
###################################################
average.samples.pima <- 
    with(models.pima,
         getBmaSamples(config=models[1:10,],
                       logPostProbs=models$logMargLik[1:10] + 
                       models$logPrior[1:10],
                       nSamples=500L,
                       modelData=modelData.pima,
                       mcmc=mcmc.pima,
                       computation=computation.pima))