Bayesian model inference for multiple fractional polynomial models

Description

Bayesian model inference for multiple fractional polynomial models is conducted by means of either exhaustive model space evaluation or posterior model sampling.

Usage

BayesMfp(formula = formula(data), data = parent.frame(), family =
gaussian, priorSpecs = list(a = 4, modelPrior = "flat"), method =
c("ask", "exhaustive", "sampling"), subset = NULL, na.action = na.omit,
verbose = TRUE, nModels = NULL, nCache=1e9L, chainlength = 1e5L)

bfp(x, max = 2, scale = TRUE, rangeVals=NULL)

uc(x)

Arguments

Details

The formula is of the form y ~ bfp (x1, max = 4) + uc (x2 + x3), that is, the auxiliary functions bfp and uc must be used for defining the fractional polynomial and uncertain fixed form covariates terms, respectively. There must be an intercept, and no other fixed covariates are allowed. All max arguments of the bfp terms must be identical.

The prior specifications are a list:

a

hyperparameter for hyper-g prior which must be greater than 3 and is recommended to be not greater than 4 (default is 4)

modelPrior

choose if a flat model prior (default, "flat"), a model prior favoring sparse models explicitly ("sparse"), or a dependent model prior ("dependent") should be used.

If method = "ask", the user is prompted with the maximum cardinality of the model space and can then decide whether to use posterior sampling or the exhaustive model space evaluation.

Note that if you specify only one FP term, the exhaustive model search must be done, due to the structure of the model sampling algorithm. However, in reality this will not be a problem as the model space will typically be very small.

Value

Returns an object of class BayesMfp that inherits from list. It is essentially a list of models. Each model is a list and has the following components:

Subsetting the object with [.BayesMfp returns again a BayesMfp object with the same attributes, which are

Note

logNormConst may be unusable due to necessary conversion from long double to double!

Various methods for posterior summaries are available.

See Also

BayesMfp Methods, BmaSamples

Examples

## generate some data
set.seed(19)

x1 <- rnorm(n=15)
x2 <- rbinom(n=15, size=20, prob=0.5) 
x3 <- rexp(n=15)

y <- rt(n=15, df=2)

## run an exhaustive model space evaluation with a flat model prior and
## a uniform prior (a = 4) on the shrinkage factor t = g/(1 + g):
test <- BayesMfp(y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 100,
                 method="exhaustive")
test

## now the same with a *dependent* model prior:
test2 <- BayesMfp(y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 100,
		 priorSpecs = list(a = 4, modelPrior = "dependent"),
                 method="exhaustive")
test2

Other methods for BayesMfp objects

Description

Print the object (print), get fitted values (fitted) and corresponding residuals (residuals).

Usage

## S3 method for class 'BayesMfp'
print(x, ...)
## S3 method for class 'BayesMfp'
fitted(object, design = getDesignMatrix(object), post =
getPosteriorParms(object, design = design), ...) 
## S3 method for class 'BayesMfp'
residuals(object, ...)

Arguments

Author(s)

Daniel Saban\'es Bov\'e

See Also

BayesMfp, BmaSamples Methods

Examples

## generate a BayesMfp object
set.seed(19)

x1 <- rnorm(n=15)
x2 <- rbinom(n=15, size=20, prob=0.5) 
x3 <- rexp(n=15)

y <- rt(n=15, df=2)

test <- BayesMfp(y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 100,
                 method="exhaustive")

## the print method
test

## extract fitted values and corresponding residuals
fitted(test)
residuals(test)

Bayesian model averaging over multiple fractional polynomial models

Description

Draw samples from the Bayesian model average over the models in saved in a BayesMfp-object.

Usage

BmaSamples(object, sampleSize = length(object) * 10, postProbs =
posteriors(object), gridList = list(), gridSize = 203, newdata=NULL,
verbose = TRUE, includeZeroSamples=FALSE)  

Arguments

Value

Return an object of class BmaSamples, which is a list with various elements that describe the BayesMfp object over which was averaged, model frequencies in the samples, the samples themselves etc:

See Also

BmaSamples Methods, BayesMfp

Examples

## construct a BayesMfp object
set.seed(19)

x1 <- rnorm (n=15)
x2 <- rbinom (n=15, size=20, prob=0.5) 
x3 <- rexp (n=15)

y <- rt (n=15, df=2)

test <- BayesMfp (y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 200, method="exhaustive")

## now draw samples from the Bayesian model average
testBma <- BmaSamples (test)
testBma

## We can also draw predictive samples for new data points, but then
## we need to supply the new data to BmaSamples:
newdata <- data.frame(x1 = rnorm(15),
                      x2 = rbinom(n=15, size=5, prob=0.2) + 1,
                      x3 = rexp(n=15))
testBma <- BmaSamples(test, newdata=newdata)
predict(testBma)

## test that inclusion of zero samples works
testBma <- BmaSamples (test, includeZeroSamples=TRUE)
testBma

Other methods for BmaSamples objects

Description

Print the object (print), get fitted values (fitted) and corresponding residuals (residuals).

Usage

## S3 method for class 'BmaSamples'
print(x, ...)
## S3 method for class 'BmaSamples'
fitted(object, ...)
## S3 method for class 'BmaSamples'
residuals(object, ...)

Arguments

Author(s)

Daniel Saban\'es Bov\'e

See Also

predict.BmaSamples, summary.BmaSamples

Examples

## construct a BayesMfp object
set.seed(19)

x1 <- rnorm (n=15)
x2 <- rbinom (n=15, size=20, prob=0.5) 
x3 <- rexp (n=15)

y <- rt (n=15, df=2)

test <- BayesMfp (y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 200, method="exhaustive")

## now draw samples from the Bayesian model average
testBma <- BmaSamples (test)

## the print method:
testBma

## the fitted method:
fitted(testBma)

## and the corresponding residuals:
residuals(testBma)

Extract method for BayesMfp objects

Description

Extract a subset of models from a BayesMfp object.

Usage

## S3 method for class 'BayesMfp'
x[...]

Arguments

Author(s)

Daniel Saban\'es Bov\'e

See Also

BayesMfp

Examples

## generate a BayesMfp object
set.seed(19)

x1 <- rnorm(n=15)
x2 <- rbinom(n=15, size=20, prob=0.5) 
x3 <- rexp(n=15)

y <- rt(n=15, df=2)

test <- BayesMfp(y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 100,
                 method="exhaustive")

## extract the top ten models
test[1:10]

Calculate and print the summary of a BayesMfp object

Description

Calculate and print the summary of a BayesMfp object, using S3 methods for the class.

Usage

## S3 method for class 'BayesMfp'
summary(object, level=0.95, table=TRUE,
                           shrinkage=NULL, ...)
## S3 method for class 'summary.BayesMfp'
print(x, ...)

Arguments

Value

summary.BayesMfp returns a list with S3 class summary.BayesMfp, where the arguments “call”, “numVisited”, “termNames”, “shiftScaleMax”, “inclusionProbs”, “chainlength” (only for model sampling results) are copied from the attributes of the BayesMfp object, please see its help page for details.

The other elements are:

If there are multiple models in object, the list element postProbs contains the exact (for exhaustively explored model spaces) or estimated (if model sampling has been done) posterior model probabilities.

If object contains only one FP model, then this one is summarized in more detail:

Note

Note that if you extract the summary of a single model with these functions, you ignore the uncertainty about the shrinkage factor t=g/(g+1) by plugging in the number shrinkage. If you want to incorporate this uncertainty, you must run BmaSamples on this model and call the corresponding method summary.BmaSamples.

Author(s)

Daniel Saban\'es Bov\'e

See Also

summary.BmaSamples

Examples

## generate a BayesMfp object
set.seed(19)

x1 <- rnorm(n=15)
x2 <- rbinom(n=15, size=20, prob=0.5) 
x3 <- rexp(n=15)

y <- rt(n=15, df=2)

test <- BayesMfp(y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 100,
                 method="exhaustive")

## summary of multiple models:
summary(test)

## summary of just one model (no. 10):
summary(test[10])

## internal structure is usually not interesting:
str(summary(test[10]))

Calculate and print the summary of a BmaSamples object

Description

Calculate and print the summary of a BmaSamples object, using S3 methods for the class.

Usage

## S3 method for class 'BmaSamples'
summary(object, level = 0.95, hpd = TRUE, ...)
## S3 method for class 'summary.BmaSamples'
print(x, table = TRUE, ...)

Arguments

Value

The summary method returns an S3 object, where “sampleSize”, “modelData” and “modelFreqs” are copied from the BmaSamples object, please see its help page for the details. “intervalType” and “level” copy the function's parameters.

“summaryMat” contains the posterior summaries for the intercept and uncertain fixed form covariates. “sigma2Sum” and “shrinkageSum” contain the posterior summaries for the regression variance and the shrinkage factor, respectively. The summaries are always the median, mean, lower and upper credible bounds for the coefficients.

Author(s)

Daniel Saban\'es Bov\'e

See Also

summary.BayesMfp

Examples

## generate a BmaSamples object
set.seed(19)

x1 <- rnorm(n=15)
x2 <- rbinom(n=15, size=20, prob=0.5) 
x3 <- rexp(n=15)

y <- rt(n=15, df=2)

test <- BayesMfp(y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 100,
                 method="exhaustive")

testBma <- BmaSamples(test)

## look at the summary
summary(testBma)

## and its structure
str(summary(testBma))

Convert a BayesMfp object to a data frame

Description

Convert the BayesMfp object to a data frame with the saved models.

Usage

## S3 method for class 'BayesMfp'
as.data.frame(x, row.names = NULL, ..., freq = TRUE)

Arguments

Author(s)

Daniel Saban\'es Bov\'e

See Also

summary.BayesMfp

Examples

## generate a BayesMfp object
set.seed(19)

x1 <- rnorm(n=15)
x2 <- rbinom(n=15, size=20, prob=0.5) 
x3 <- rexp(n=15)

y <- rt(n=15, df=2)

test <- BayesMfp(y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 100,
                 method="exhaustive")

## get the models data frame
as.data.frame(test)

BMA prediction for new data points

Description

Make a Bayesian model averaged prediction for new data points, from those models saved in a BayesMfp object.

Usage

bmaPredict(BayesMfpObject, postProbs = posteriors(BayesMfpObject), newdata)

Arguments

Value

The predicted values as a vector.

Note

Note that this function is not an S3 predict method for BmaSamples objects, but a function working on BayesMfp objects (because we do not need BMA samples to do BMA point predictions).

Author(s)

Daniel Saban\'es Bov\'e

See Also

BmaSamples Methods

Examples

## generate a BayesMfp object
set.seed(19)

x1 <- rnorm(n=15)
x2 <- rbinom(n=15, size=20, prob=0.5) 
x3 <- rexp(n=15)

y <- rt(n=15, df=2)

test <- BayesMfp(y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 100,
                 method="exhaustive")

## predict new responses at (again random) covariates
bmaPredict(test,
           newdata = list(x1 = rnorm(n=15),
                          x2 = rbinom(n=15, size=5, prob=0.2) + 1,
                          x3 = rexp(n=15)))

Construct the model matrix for new data based on an existing BayesMfp object

Description

This is an internal function, which constructs a model matrix for new covariate data, based on the formula and scaling info in an existing BayesMfp object. The matrix can then be passed to prediction functions.

Usage

constructNewdataMatrix(BayesMfpObject, newdata)

Arguments

Value

The (uncentered!) model matrix with the FP columns shifted and scaled as for the original data.

Author(s)

Daniel Saban\'es Bov\'e

Examples

## construct a BayesMfp object
set.seed(19)

x1 <- rnorm (n=15)
x2 <- rbinom (n=15, size=20, prob=0.5) 
x3 <- rexp (n=15)

y <- rt (n=15, df=2)

test <- BayesMfp (y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 200, method="exhaustive")

## some new covariate data
newdata <- data.frame(x1 = 1:10, x2 = rbinom(n=10, size=20, prob=0.6),
x3 = 2:11)

## construct the new design matrix:
newmatrix <- bfp:::constructNewdataMatrix(BayesMfpObject=test, newdata=newdata)

## check the result:
stopifnot(identical(newmatrix,
                    structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.6,
				1.2, 1.2, 1.2,  1.3, 1.3, 1.4, 1.1,
				1.2, 1.1, 1, 2, 3, 4, 5, 6, 7, 8, 9,
				10,  2, 3, 4, 5, 6, 7, 8, 9, 10, 11),
				.Dim = c(10L, 4L),
				.Dimnames =
                                 list(c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"),
				      c("(Intercept)", "x2", "x1", "x3")),
				assign = 0:3)))

Functions for the inverse gamma distribution.

Description

Functions for the inverse gamma distribution, analogues to the gamma distribution functions.

Usage

dinvGamma(x, a, b, log = FALSE, normalize = TRUE)
pinvGamma(q, a, b, lower.tail = TRUE, log.p = FALSE)
qinvGamma(p, a, b, lower.tail = TRUE, log.p = FALSE)
rinvGamma(n, a, b)

Arguments

Details

See the gamma distribution (rgamma) for the details.

Author(s)

Daniel Saban\'es Bov\'e

Construct an empirical HPD interval from samples

Description

Construct an empirical highest posterior density (HPD) interval from samples which have been drawn from the distribution of a quantity of interest.

Usage

empiricalHpd(theta, level)

Arguments

Value

A vector with the estimated lower and upper bounds of the HPD interval.

Author(s)

Daniel Saban\'es Bov\'e

Examples

## draw standard normal variates
test <- rnorm(n=1000)

## estimate the 95% HPD interval with these samples:
empiricalHpd(theta=test, level=0.95)

## compare with true HPD:
qnorm(p=c(0.025, 0.975))

Find a specific fractional polynomial model in a BayesMfp object

Description

Returns the index of the wished model if it is present in the model list, and otherwise returns NA.

Usage

findModel(model, BayesMfpObject)

Arguments

Details

See BayesMfp for the description of a model.

Value

Index of model in BayesMfpObject if it is present in the model list, otherwise NA.

Note

The searched model must have exactly the same construction as the models in BayesMfpObject. See the example below for the recommended use.

Examples

## construct a BayesMfp object
set.seed(92)

x1 <- rnorm (15)
x2 <- rbinom (n=15, size=20, prob=0.6)
x3 <- rexp (15)
y <- rt (15, df=2)

test <- BayesMfp (y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels=2000, method="exhaustive")

## copy one model 
myModel <- test[[1]]

## and modify it!
myModel$powers[["x2"]] <- c (1, 2)
stopifnot(identical(findModel (myModel, test),
                    31L))

Internal function for shifting and rescaling a fractional polynomial covariate

Description

This function is derived from mfp:::fp.scale written by Gareth Ambler and Axel Benner.

Usage

fpScale(x, scaling = TRUE)

Arguments

Extract the design matrix of a multiple FP model.

Description

Extract the (centered) design matrix of the first element (== model) of a BayesMfp object. This is an internal function not intended to be used publicly.

Usage

getDesignMatrix(x, fixedColumns=TRUE, center=TRUE)

Arguments

Value

The design matrix with an attribute shifts containing the used shifts for the optional centering of the non-intercept columns. (If center is FALSE, the shifts vector will contain only zeroes.)

Author(s)

Daniel Saban\'es Bov\'e

See Also

getFpTransforms

Examples

## construct a BayesMfp object
set.seed(19)

x1 <- rnorm (n=15)
x2 <- rbinom (n=15, size=20, prob=0.5) 
x3 <- rexp (n=15)

y <- rt (n=15, df=2)

test <- BayesMfp (y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 200, method="exhaustive")

## get the design matrix of the fifth best model
a <- bfp:::getDesignMatrix(test[5])
attr(a, "shifts")

## and once again but without centering
b <- bfp:::getDesignMatrix(test[5], center=FALSE)
stopifnot(all(attr(b, "shifts") == 0))

Internal functions to handle FP transforms.

Description

Transform a variable according to FP transformation formula and attach proper names to the resulting design matrix. The binary function bt is the Box-Tidwell transform.

Usage

x %bt% pow
getFpTransforms(vec, powers, center=TRUE)
getTransformName(name, pow)

Arguments

Author(s)

Daniel Saban\'es Bov\'e

See Also

getDesignMatrix

Examples

bfp:::getFpTransforms(cbind(x=1:10), powers=c(-2, 0, 1))

Extract log marginal likelihood from a model.

Description

Extract log marginal likelihood from a model, saved as the first element of a BayesMfp object.

Usage

getLogMargLik(x, design=getDesignMatrix(x), nObs = nrow(design), dim = ncol(design))

Arguments

Details

This function interfaces the C++ function logMargLik, and can be used to compute the marginal likelihood of a model not saved in the model list. But be careful to adjust the saved R^2 of the model, too, and not only the powers! Therefore this function is internal only... and is used e.g. in transformMfp.

Author(s)

Daniel Saban\'es Bov\'e

See Also

getLogPrior

Extract the log prior probability from a model.

Description

Extract the log prior probability from the first model saved in a BayesMfp object.

Usage

getLogPrior(x)

Arguments

Details

Note that the returned value might be only defined up to a constant, because some models might not be identifiable or have a dimension too large, and are thus implicitly assigned a zero prior probabibility.

Author(s)

Daniel Saban\'es Bov\'e

See Also

getLogMargLik

Extract posterior expected g and shrinkage factor (g/(1+g)) from a model

Description

For a valid BayesMfp object, extract posterior expected g and shrinkage factor g/(1+g)

Usage

getPostExpectedg(x, design=getDesignMatrix(x), nObs = nrow(design),
dim = ncol(design))
getPostExpectedShrinkage(x, design=getDesignMatrix(x),
nObs =nrow(design), dim = ncol(design))

Arguments

Value

The posterior expected g or shrinkage factor g/(1+g) from the model.

Author(s)

Daniel Saban\'es Bov\'e

Extract updated posterior parameters for the normal inverse gamma distribution from a model, given a shrinkage factor.

Description

Conditional on a fixed shrinkage factor t=g/(g+1), the posterior joint distribution of the effects and the regression variance is normal inverse gamma. With this function, you can compute the parameters of this distribution.

Usage

getPosteriorParms(x, shrinkage=x[[1]]$postExpectedShrinkage,  
                  design = getDesignMatrix(x))

Arguments

Value

A list with four parameters:

Author(s)

Daniel Saban\'es Bov\'e

Examples

## construct a BayesMfp object
set.seed(19)

x1 <- rnorm (n=15)
x2 <- rbinom (n=15, size=20, prob=0.5) 
x3 <- rexp (n=15)

y <- rt (n=15, df=2)

test <- BayesMfp (y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 200, method="exhaustive")

## now get the posterior parameters of the third best model
getPosteriorParms(test[3])

Compute (model averaged) posterior variable inclusion probabilites

Description

Compute (model averaged) posterior inclusion probabilites for the uncertain variables (including FP variables) based on a BayesMfp object.

Usage

inclusionProbs(BayesMfpObject, postProbs = posteriors(BayesMfpObject, ind = 1))

Arguments

Value

Named numeric vector with the estimated variable inclusion probabilities. Note that these can differ noticeably from the “global” inclusion probabilities computed from all discovered models, from which only the best were saved in the BayesMfp object.

Author(s)

Daniel Saban\'es Bov\'e

Examples

## construct a BayesMfp object
set.seed(19)

x1 <- rnorm (n=15)
x2 <- rbinom (n=15, size=20, prob=0.5) 
x3 <- rexp (n=15)

y <- rt (n=15, df=2)

test <- BayesMfp (y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 200, method="exhaustive")

## now get the local inclusion probabilities
local <- inclusionProbs(test)

## they can be compared with the global inclusion probabilities
local - attr(test, "inclusionProbs")

Ozone data from Breiman and Friedman, 1985

Description

This is the Ozone data discussed in Breiman and Friedman (JASA, 1985, p. 580). These data are for 330 days in 1976. All measurements are in the area of Upland, CA, east of Los Angeles.

Usage

data(ozone)

Format

A data frame with 366 observations on the following 13 variables.

month

month of the year

day

day of the month

weekday

day of the week: a factor with levels Monday, Tueday, Wednesday, Thursday, Friday, Saturday, Sunday

hourAverageMax

maximum 1-hour average ozone level [ppm]

pressure500Height

500 millibar pressure height [meters]

windSpeed

wind speed [mph]

humidity

relative humidity [%]

tempSandburg

temperature at Sandberg, CA [degrees F]

tempElMonte

temperature at El Monte, CA [degrees F]

inversionBaseHeight

inversion base height [feet]

pressureGradientDaggett

pressure gradient from LAX to Daggett, CA [mm Hg]

inversionBaseTemp

inversion base temperature [degrees F]

visibility

visibility [miles]

Source

Breiman, L and Friedman, J. (1985), “Estimating Optimal Transformations for Multiple Regression and Correlation”, Journal of the American Statistical Association, 80, 580-598.

Generic function for plotting a fractional polynomial curve estimate

Description

Plot a fractional polynomial curve estimate for either a single model or a Bayesian model average over BayesMfp objects. Optionally, credible intervals and / or bands can be added to the plot.

Usage

plotCurveEstimate(model, termName, plevel = 0.95, slevel = plevel,
plot = TRUE, legendPos = "topleft", rug = FALSE, partialResids=TRUE, 
hpd=TRUE,..., main = NULL) 

Arguments

Details

Further arguments for application on a BayesMfp object:

grid

vector of unscaled abscissae, default is a length gridSize grid over the observed range specified by providing the argument NULL.

post

list with posterior parameters of the model, which may be provided manually to accelerate plotting in a loop

gridSize

default number of grid points used when no grid is supplied (default is 201)

numSim

number of simulations for estimation of the SCB (default is 500)

Value

a list of various plotting information:

See Also

BayesMfp, BmaSamples

Examples

## construct a BayesMfp object
set.seed(19)
x1 <- rnorm (n=15)
x2 <- rbinom (n=15, size=20, prob=0.5) 
x3 <- rexp (n=15)
y <- rt (n=15, df=2)

test <- BayesMfp (y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 100,
method="exhaustive")

## plot the x2 curve estimate for the 20-th best model
p1 <- plotCurveEstimate (test[20], "x2")

## look at the returned list
str(p1)

## plot the BMA curve estimate for the same covariate
testBma <- BmaSamples (test)
p2 <- plotCurveEstimate (testBma, "x2")

## look at the returned list
str(p2)

## try the new options:
plotCurveEstimate (testBma, "x2", partialResids=FALSE, hpd=FALSE)

Extract posterior model probability estimates from BayesMfp objects

Description

Extract posterior model probability estimates (either normalized estimates or sampling frequencies) from BayesMfp objects.

Usage

posteriors(BayesMfpObject, ind = 1)

Arguments

Value

The vector of probability estimates.

Author(s)

Daniel Saban\'es Bov\'e

Examples

## construct a BayesMfp object
set.seed(19)
x1 <- rnorm (n=15)
x2 <- rbinom (n=15, size=20, prob=0.5) 
x3 <- rexp (n=15)
y <- rt (n=15, df=2)

test <- BayesMfp (y ~ bfp (x1, max = 2) + bfp (x2, max = 2) + uc (x3), nModels = 100,
		  method="exhaustive")

## this works:
posteriors(test)


## only if we do model sampling there are model frequencies:
test2 <- BayesMfp (y ~ bfp (x1, max = 2) + bfp (x2, max = 2) + uc (x3), nModels = 100,
         method="sampling")
posteriors(test2, ind=2)

Predict method for BayesMfp objects

Description

Predict new responses from a single multiple FP model.

Usage

## S3 method for class 'BayesMfp'
predict(object, newdata, ...)

Arguments

Author(s)

Daniel Saban\'es Bov\'e

See Also

bmaPredict

Examples

## generate a BayesMfp object
set.seed(19)

x1 <- rnorm(n=15)
x2 <- rbinom(n=15, size=20, prob=0.5) 
x3 <- rexp(n=15)

y <- rt(n=15, df=2)

test <- BayesMfp(y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 100,
                 method="exhaustive")

## predict new responses at (again random) covariates
predict(test,
        newdata = list(x1 = rnorm (15),
                       x2 = rbinom (n=15, size=5, prob=0.2) + 1,
                       x3 = rexp (15)))

Predict method to extract point and interval predictions from BmaSamples objects

Description

Predict new responses from a Bayesian model average over FP models, from which predictive samples have already been produced.

Usage

## S3 method for class 'BmaSamples'
predict(object, level=0.95, hpd=TRUE, ...)
## S3 method for class 'predict.BmaSamples'
print(x, ...)

Arguments

Details

This function summarizes the predictive samples saved in the BmaSamples object. Using these functions, one can obtain predictive credible intervals, as opposed to just using the function bmaPredict, which only gives the means of the predictive distributions.

Value

A list of class predict.BmaSamples, which has then a separate print method. The elements of the list are:

Author(s)

Daniel Saban\'es Bov\'e

See Also

bmaPredict

Examples

## generate a BmaSamples object
set.seed(19)

x1 <- rnorm(n=15)
x2 <- rbinom(n=15, size=20, prob=0.5) 
x3 <- rexp(n=15)

y <- rt(n=15, df=2)

test <- BayesMfp(y ~ bfp (x2, max = 4) + uc (x1 + x3), nModels = 100,
                 method="exhaustive")

## predict new responses at (again random) covariates with BMA:
testBma <- BmaSamples(test,
                      newdata=data.frame(x1 = rnorm (15),
                                         x2 = rbinom (n=15, size=5, prob=0.2) + 1,
                                         x3 = rexp (15)))
predict(testBma)

Multivariate Student Random Deviates

Description

This function provides a random vector generator for the multivariate Student distribution with mean vector mean, scale matrix sigma and degrees of freedom df.

Usage

rmvt(n, sigma = diag(2), mu = rep(0, 2), df = 1)

Arguments

Value

The vector samples, stacked in the rows of a matrix.

Note

Note that the covariance matrix is only defined for df > 2, and is given by sigma * df / (df - 2).

Author(s)

Daniel Saban\'es Bov\'e

Examples

## samples from the multivariate Cauchy distribution:
bfp:::rmvt(20)

## here the covariance exists:
sigma <- matrix(c(1, 0.5, 0.5, 1), nrow=2)
df <- 10

## theoretical covariance:
sigma * df / (df - 2)

## this should be close:
cov(bfp:::rmvt(n=100000, sigma=sigma, df=df))

Sample from the model-specific posterior of the shrinkage factor t = g / (1 + g)

Description

Sample from the model-specific posterior of the shrinkage factor t = g / (1 + g), using inverse sampling. This function is internal only because the arguments are not very user-friendly.

Usage

rshrinkage(n=1, R2, nObs, p, alpha)

Arguments

Value

n posterior shrinkage factor samples

Author(s)

Daniel Saban\'es Bov\'e

Examples

## construct a BayesMfp object
set.seed(29)
nObs <- 15

x1 <- rnorm (n=nObs)
x2 <- rbinom (n=nObs, size=20, prob=0.5) 
x3 <- rexp (n=nObs)

y <- rt (n=nObs, df=2)

test <- BayesMfp (y ~ bfp (x2, max = 3) + bfp(x3, max=3), nModels = 200, method="exhaustive")

## get the best found model
map <- test[[1]]

## produce many samples from the shrinkage factor
samples <- bfp:::rshrinkage(n=1e+5,
			    R2=map$R2,
		            nObs=nObs,
			    p=length(unlist(map$powers)),
			    alpha=4)
hist(samples,
     nclass=50, prob=TRUE)

## compare the mean with the analytical mean
abline(v=c(sampled <- mean(samples),
	   true <- map$postExpectedShrinkage),
       col=c(1, 2))
stopifnot(abs(sampled - true) < 0.0001)

Simultaneous credible band computation (Besag, Green et al algorithm)

Description

Simultaneous credible band computation

Usage

scrBesag(samples, level=0.95)

Arguments

Details

Calculates a series of simultaneous credible bounds for one parameter type, following section 6.3 in the seminal paper "Bayesian computation and stochastic systems". The corresponding algorithm was invented by Peter Green.

Value

matrix with ‘lower’ and ‘upper’ rows

Author(s)

Thomas Kneib

References

J. Besag, P. Green, D. Higdon, K. Mengersen (1995): Bayesian computation and stochastic systems, Statistical Science 10/1, 3–41

Calculate an SCB from a samples matrix

Description

Calculate an SCB from a samples matrix, which minimizes the absolute distances of the contained samples to a mode vector, at each gridpoint. Therefore the SCB might be considered an “HPD SCB”.

Usage

scrHpd(samples, mode = apply(samples, 2, median), level = 0.95)

Arguments

Details

This function first computes the matrix of absolute distances of the samples to the mode vector. Then based on this distance matrix, a one-sided SCB as described in Besag et al. (1995) is computed, which is then mapped back to the samples.

Value

A matrix with rows “lower” and “upper”, with the lower and upper SCB bounds.

Author(s)

Daniel Saban\'es Bov\'e

References

Besag, J.; Green, P.; Higdon, D. and Mengersen, K. (1995): “Bayesian computation and stochastic systems (with discussion)”, Statistical Science, 10, 3-66.

See Also

empiricalHpd

Examples

## create some samples
time <- 1:10
nSamples <- 50
samples <- t(replicate(nSamples,
                       time * rnorm(1) + rexp(1))) +
           rnorm(length(time) * nSamples)
matplot(time, t(samples), type="l", lty=1, col=1,
        xlab="time", ylab="response")

## now test the function: 50% credible band
scb <- scrHpd(samples, level=0.5)
matlines(time, t(scb), col=2, lwd=2, lty=1)

Transform a fitted mfp model into a BayesMfp model

Description

Transform a fitted mfp model into a BayesMfp model with the correct powers etc. to compare with other (true) BayesMfp models fitted to the same (!) data.

Usage

transformMfp(mfpObject, BayesMfpObject)

Arguments

Value

A BayesMfp object with the converted model.