HIV Survival data

Description

Data from a hypothetical HMO-HIV+ study shown in Table 1.1 of Hosmer, D.W. and Lemeshow, S. (1998) Applied Survival Analysis: Regression Modeling of Time to Event Data, John Wiley and Sons Inc., New York, NY

Format

A data frame with 100 observations on 7 variables.

Bayesian Cox Proportional Hazards Modelling

Description

Uses a Metropolis Hastings scheme on the proportional hazards model to draw sample from posterior. Uses a matched curvature Student's t candidate generating distribution with 4 degrees of freedom to give heavy tails.

Usage

BayesCPH(
  y,
  t,
  x,
  steps = 1000,
  priorMean = NULL,
  priorVar = NULL,
  mleMean = NULL,
  mleVar,
  startValue = NULL,
  randomSeed = NULL,
  plots = FALSE
)

Arguments

Value

A list containing the following components:

Bayesian Logistic Regression

Description

Performas Metropolis Hastings on the logistic regression model to draw sample from posterior. Uses a matched curvature Student's t candidate generating distribution with 4 degrees of freedom to give heavy tails.

Usage

BayesLogistic(
  y,
  x,
  steps = 1000,
  priorMean = NULL,
  priorVar = NULL,
  mleMean = NULL,
  mleVar,
  startValue = NULL,
  randomSeed = NULL,
  plots = FALSE
)

Arguments

Value

A list containing the following components:

Examples

data(logisticTest.df)
BayesLogistic(logisticTest.df$y, logisticTest.df$x)

Bayesian Pois Regression

Description

Performs Metropolis Hastings on the logistic regression model to draw sample from posterior. Uses a matched curvature Student's t candidate generating distribution with 4 degrees of freedom to give heavy tails.

Usage

BayesPois(
  y,
  x,
  steps = 1000,
  priorMean = NULL,
  priorVar = NULL,
  mleMean = NULL,
  mleVar,
  startValue = NULL,
  randomSeed = NULL,
  plots = FALSE
)

Arguments

Value

A list containing the following components:

Examples

data(poissonTest.df)
results = BayesPois(poissonTest.df$y, poissonTest.df$x)

Calculate the Gelman Rubin statistic

Description

Calculate the Gelman Rubin statistic

Usage

GelmanRubin(theta)

Arguments

Value

A list containing n, the between chain variance B, the within chain variance W, the estimated variance of the parameter vHat, and the Gelman Rubin statistic R = \sqrt{vHat/W}

References

Gelman, A. and Rubin, D.B. (1992) 'Inference from iterative simulations using multiple sequences with discussion.' Statistical Science 8, pp. 457-511

Examples

## take four chains sampling from a normal mixture density
theta0 = c(0,1)
theta1 = c(3,2)
p = 0.6
candidate = c(0, 3)

v1 = normMixMH(theta0, theta1, p, candidate, steps = 200)
v2 = normMixMH(theta0, theta1, p, candidate, steps = 200)
v3 = normMixMH(theta0, theta1, p, candidate, steps = 200)
v4 = normMixMH(theta0, theta1, p, candidate, steps = 200)

theta=cbind(v1,v2,v3,v4)
GelmanRubin(theta)

Metropolis Hastings sampling from a Bivariate Normal distribution

Description

This function uses the MetropolisHastings algorithm to draw a sample from a correlated bivariate normal target density using a random walk candidate and an independent candidate density respectively where we are drawing both parameters in a single draw. It can also use the blockwise Metropolis Hastings algorithm and Gibbs sampling respectively to draw a sample from the correlated bivariate normal target.

Usage

bivnormMH(rho, rho1 = 0.9, sigma = c(1.2, 1.2), steps = 1000, type = "ind")

Arguments

Value

returns a list which contains a data frame called targetSample with members x and y. These are the samples from the target density.

Examples

## independent chain
chain1.df=bivnormMH(0.9)$targetSample

## random walk chain
chain2.df=bivnormMH(0.9, type = 'r')$targetSample


## blockwise MH chain
chain3.df=bivnormMH(0.9, type = 'b')$targetSample

## Gibbs sampling chain
chain4.df=bivnormMH(0.9, type = 'g')$targetSample

oldPar = par(mfrow=c(2,2))
plot(y ~ x, type = 'l', chain1.df, main = 'Independent')
plot(y ~ x, type = 'l', chain2.df, main = 'Random Walk')
plot(y ~ x, type = 'l', chain3.df, main = 'Blockwise')
plot(y ~ x, type = 'l', chain4.df, main = 'Gibbs')
par(oldPar)

Chapter 10 Example 16 data

Description

A random sample of size 10 from a N(\mu, \sigma^{2}) distribution where both mu and sigma are unknown parameters.

Format

A data frame with 10 observations in a single variable called y

Coronary Heart Disease Chapter 8 Example 11

Description

The age and coronory heart disease status of 100 individuals taken from Hosmer and Lemeshow (1989).

Format

A data frame with 100 observations in two columns

Calculate a credible interval from a numerically specified posterior CDF or from a sample from the posterior

Description

Calculates a lower, upper, or two-sided credible interval from the numerical posterior CDF or from a sample from the posterior.

Usage

credInt(theta, cdf = NULL, conf = 0.95, type = "twosided")

Arguments

Details

This function uses linear interpolation to calculate bounds for points that may not be specified by CDF

Value

a list containing the elements lower.bound, uppper.bound or both depending on type

Examples

## commands for calculating a numerical posterior CDF.
## In this example, the likelihood is proportional to
## \eqn{\theta^{3/2}\times \exp(-\theta/4)} and a N(6, 9) prior is used.
theta = seq(from = 0.001, to = 40, by = 0.001)
prior = dnorm(theta,6,3)
ppnLike = theta^1.5*exp(-theta/4)
ppnPost = prior*ppnLike
scaleFactor = sintegral(theta, ppnPost)$int
posterior = ppnPost/scaleFactor
cdf = sintegral(theta, posterior)$y
ci = credInt(theta, cdf)
par(mfrow=c(2,2))
plot(prior ~ theta, type = 'l',  main = 'Prior N(6, 9)')
plot(ppnLike ~ theta, type = 'l', main = 'Proportional likelihood')
plot(posterior ~ theta, type = 'l', main = 'Posterior')
abline(v=c(unlist(ci)))

## Use an inverse method to take a random sample of size 1000
## from the posterior
suppressWarnings({Finv = approxfun(cdf,theta)})
thetaSample = Finv(runif(1000))
ci = credInt(thetaSample)

Calculate a credible interval from a numerically specified posterior CDF

Description

Calculates a lower, upper, or two-sided credible interval from the numerical posterior CDF.

Usage

credIntNum(theta, cdf, conf = 0.95, type = "twosided")

Arguments

Details

This function uses linear interpolation to calculate bounds for points that may not be specified by CDF

Value

a list containing the elements lower.bound, uppper.bound or both depending on type

Examples

## commands for calculating a numerical posterior CDF.
## In this example, the likelihood is proportional to
## \eqn{\theta^{3/2}\times \exp(-\theta/4)} and a N(6, 9) prior is used.
theta = seq(from = 0.001, to = 40, by = 0.001)
prior = dnorm(theta,6,3)
ppnLike = theta^1.5*exp(-theta/4)
ppnPost = prior*ppnLike
scaleFactor = sintegral(theta, ppnPost)$int
posterior = ppnPost/scaleFactor
cdf = sintegral(theta, posterior)$y
ci=credIntNum(theta, cdf)
par(mfrow=c(2,2))
plot(prior ~ theta, type = 'l',  main = 'Prior N(6, 9)')
plot(ppnLike ~ theta, type = 'l', main = 'Proportional likelihood')
plot(posterior ~ theta, type = 'l', main = 'Posterior')
abline(v=c(unlist(ci)))

Calculate a credible interval from a numerically specified posterior CDF

Description

Calculates a lower, upper, or two-sided credible interval from the numerical posterior CDF.

Usage

credIntSamp(theta, conf = 0.95, type = "twosided")

Arguments

Details

This function uses linear interpolation to calculate bounds for points that may not be specified by CDF

Value

a list containing the elements lower.bound, uppper.bound or both depending on type

Examples

## posterior is N(0,1)
theta = rnorm(1000)
ci=credIntSamp(theta)
plot(density(theta))
abline(v=c(unlist(ci)))

Give simple descriptive statistics for a matrix or a data frame

Description

This function is designed to emulate the Minitab function DESCRIBE. It gives simple descriptive statistics for a data frame

Usage

describe(x, varNames = NULL)

Arguments

Value

A data.frame containing the following elements:

Examples

data(poissonTest.df)
describe(poissonTest.df)

Hierarchical Normal Means Regression Model

Description

fits a hierarchical normal model of the form E[y_{ij}] = \mu_{j} + \beta_{1}x_{i1}+\dots+\beta_{p}x_{ip}

Usage

hierMeanReg(
  design,
  priorTau,
  priorPsi,
  priorVar,
  priorBeta = NULL,
  steps = 1000,
  startValue = NULL,
  randomSeed = NULL
)

Arguments

Value

A data frame with variables:

Examples

priorTau = list(tau0 = 0, v0 = 1000)
priorPsi = list(psi0 = 500, eta0 = 1)
priorVar = list(s0 = 500, kappa0 = 1)
priorBeta = list(b0 = c(0,0), bMat = matrix(c(1000,100,100,1000), ncol = 2))

data(hiermeanRegTest.df)
data.df = hiermeanRegTest.df
design = list(y = data.df$y, group = data.df$group,
               x = as.matrix(data.df[,3:4]))
r=hierMeanReg(design, priorTau, priorPsi, priorVar, priorBeta)

oldPar = par(mfrow = c(3,3))
plot(density(r$tau))
plot(density(r$psi))
plot(density(r$mu.1))
plot(density(r$mu.2))
plot(density(r$mu.3))
plot(density(r$beta.1))
plot(density(r$beta.2))
plot(density(r$sigmaSq))
par(oldPar)

## example with no covariates
priorTau = list(tau0 = 0, v0 = 1000)
priorPsi = list(psi0 = 500, eta0 = 1)
priorVar = list(s0 = 500, kappa0 = 1)

data(hiermeanRegTest.df)
data.df = hiermeanRegTest.df
design = list(y = data.df$y, group = data.df$group, x = NULL)
r=hierMeanReg(design, priorTau, priorPsi, priorVar)

oldPar = par(mfrow = c(3,2))
plot(density(r$tau))
plot(density(r$psi))
plot(density(r$mu.1))
plot(density(r$mu.2))
plot(density(r$mu.3))
plot(density(r$sigmaSq))
par(oldPar)

Test data for hiermeanReg

Description

Data for testing hiermeanReg which uses Gibbs sampling on a hierarchical normal mean model with regression on covariates

Format

A data frame with 30 observations on 4 variables.

See Also

hiermeanReg

Test data for bayesLogistic

Description

A test data set for bayesLogisticReg

Format

A data frame with 100 observations on 6 variables.

See Also

bayesLogistic

Draw a sample from a posterior distribution of data with an unknown mean and variance using Gibbs sampling

Description

normGibbs draws a Gibbs sample from the posterior distribution of the parameters given the data fron normal distribution with unknown mean and variance. The prior for \mu given var is prior mean m0 and prior variance var/n0 . That means n0 is the 'equivalent sample size.' The prior distribution of the variance is s0 times an inverse chi-squared with kappa0 degrees of freedom. The joint prior is the product g(var)g(mu|var).

Usage

normGibbs(y, steps = 1000, type = "ind", ...)

Arguments

Value

A data frame containing three variables

Author(s)

James M. Curran

Examples

## firstly generate some random data
mu = rnorm(1)
sigma = rgamma(1,5,1)
y = rnorm(100, mu, sigma)

## A \eqn{N(10,3^2)} prior for \eqn{\mu} and a 25 times inverse chi-squared
## with one degree of freedom prior for \eqn{\sigma^2}
MCMCSampleInd = normGibbs(y, steps = 5000, priorMu = c(10,3),
                           priorVar = c(25,1))


## We can also use a joint conjugate prior for \eqn{\mu} and \eqn{\sigma^2}.
## This will be a \emph{normal}\eqn{(m,\sigma^2/n_0)} prior for \eqn{\mu} given
## the variance \eqn{\sigma^2}, and an \eqn{s0} times an \emph{inverse
## chi-squared} prior for \eqn{\sigma^2}.
MCMCSampleJoint = normGibbs(y, steps = 5000, type = 'joint',
                             priorMu = c(10,3), priorVar = c(25,1))

## Now plot the results
oldPar = par(mfrow=c(2,2))

plot(density(MCMCSampleInd$mu),xlab=expression(mu), main =
'Independent')
abline(v=mu)
plot(density(MCMCSampleInd$sig),xlab=expression(sig), main =
'Independent')
abline(v=sigma)

plot(density(MCMCSampleJoint$mu),xlab=expression(mu), main =
'Joint')
abline(v=mu)
plot(density(MCMCSampleJoint$sig),xlab=expression(sig), main =
'Joint')
abline(v=sigma)

Sample from a normal mixture model using Metropolis-Hastings

Description

normMixMH uses the Metropolis-Hastings algorithm to draw a sample from a univariate target distribution that is a mixture of two normal distributions using an independent normal candidate density or a random walk normal candidate density.

Usage

normMixMH(
  theta0,
  theta1,
  p,
  candidate,
  steps = 1000,
  type = "ind",
  randomSeed = NULL,
  startValue = NULL
)

Arguments

Value

A vector containing a sample from the normal mixture distribution.

Examples

## Set up the normal mixture
theta0 = c(0,1)
theta1 = c(3,2)
p = 0.8

## Sample from an independent N(0,3^2) candidate density
candidate = c(0, 3)
MCMCsampleInd = normMixMH(theta0, theta1, p, candidate)


## If we wish to use the alternative random walk N(0, 0.5^2)
## candidate density
candidate = c(0, 0.5)
MCMCsampleRW = normMixMH(theta0, theta1, p, candidate, type = 'rw')

Test a one sided hypothesis from a numerically specified posterior CDF or from a sample from the posterior

Description

Calculates the probability of a one sided null hypothesis from a numerically calculated posterior CDF or from a sample from the posterior.

Usage

pNull(theta0, theta, cdf = NULL, type = "upper")

Arguments

Details

This function uses linear interpolation to calculate bounds for points that may not be specified by CDF

Value

a list containing the element prob which will be the upper or lower tail probability depending on type

Examples

## commands for calculating a numerical posterior CDF.
## In this example, the likelihood is proportional to
## \eqn{\theta^{3/2}\times \exp(-\theta/4)} and a N(6, 9) prior is used.
theta = seq(from = 0.001, to = 40, by = 0.001)
prior = dnorm(theta,6,3)
ppnLike = theta^1.5*exp(-theta/4)
ppnPost = prior*ppnLike
scaleFactor = sintegral(theta, ppnPost)$int
posterior = ppnPost/scaleFactor
cdf = sintegral(theta, posterior)$y
pNull(15, theta, cdf)

## Use an inverse method to take a random sample of size 1000
## from the posterior
suppressWarnings({Finv = approxfun(cdf, theta)})
thetaSample = Finv(runif(1000))
pNull(15, thetaSample)

Test a one sided hypothesis from a numerically specified posterior CDF

Description

Calculates the probability of a one sided null hypothesis from a numerically calculated posterior CDF.

Usage

pnullNum(theta0, theta, cdf, type = "upper")

Arguments

Details

This function uses linear interpolation to calculate bounds for points that may not be specified by CDF

Value

a list containing the element prob which will be the upper or lower tail probability depending on type

Examples

## commands for calculating a numerical posterior CDF.
## In this example, the likelihood is proportional to
## \eqn{\theta^{3/2}\times \exp(-\theta/4)} and a N(6, 9) prior is used.
theta = seq(from = 0.001, to = 40, by = 0.001)
prior = dnorm(theta,6,3)
ppnLike = theta^1.5*exp(-theta/4)
ppnPost = prior*ppnLike
scaleFactor = sintegral(theta, ppnPost)$int
posterior = ppnPost/scaleFactor
cdf = sintegral(theta, posterior)$y
pnullNum(1, theta, cdf)

Test a one sided hypothesis using a sample from a posterior density

Description

Calculates the probability of a one sided null hypothesis from a sample from a posterior density.

Usage

pnullSamp(theta, theta0 = 0, type = "upper")

Arguments

Details

This function uses linear interpolation to calculate bounds for points that may not be specified by CDF

Value

a list containing the element prob which will be the upper or lower tail probability depending on type

Examples

## The posterior density is N(3,1)
theta = rnorm(1000,3)

## test whether the true mean is greater than 0 (it is obviously!)
pnullSamp(theta)

A test data set for bayesPois

Description

A test data set for bayesPois. The data come from the equation \log(\lambda_{i}) = 1 + 5x_{i} + \epsilon_{i} where \epsilon_{i} comes from N(0,0.01).

Format

A data frame with 100 observations on 5 variables.

See Also

bayesPois

Numerical integration using Simpson's Rule

Description

Takes a vector of x values and a corresponding set of postive f(x)=y values and evaluates the area under the curve:

\int{f(x)dx}

.

Usage

sintegral(x, fx, n.pts = 256)

Arguments

Value

returns a list with the following elements

Examples

## integrate the normal density from -3 to 3
x=seq(-3,3,length=100)
fx=dnorm(x)
estimate=sintegral(x,fx)$int
true.val=diff(pnorm(c(-3,3)))
cat(paste('Absolute error :',round(abs(estimate-true.val),7),'\n'))
cat(paste('Relative percentage error :', 100*round((abs(estimate-true.val)/true.val),6),'%\n'))

Thin an MCMC sample

Description

Thins the output from an MCMC process

Usage

thin(x, k)

Arguments

Details

Note this function does not check to see if k is sensible.

Value

A thinned vector, matrix or data frame containing every kth element of x.

Examples

## A blockwise Metropolis-Hastings chain of 1000 elements, thinned to
## 5th element
##

MCMCSampleBW = bivnormMH(0.9, type = 'block')
MCMCSampleBW = thin(MCMCSampleBW, 5)