NPP (version 0.6.0)

LMOMNPP_MCMC1: MCMC Sampling for Linear Regression Model of multiple historical data using Ordered Normalized Power Prior

Description

Multiple historical data are incorporated together. Conduct posterior sampling for Linear Regression Model with ordered normalized power prior. For the power parameter \(\gamma\), a Metropolis-Hastings algorithm with independence proposal is used. For the model parameters \((\beta, \sigma^2)\), Gibbs sampling is used.

Usage

LMOMNPP_MCMC1(D0, X, Y, a0, b, mu0, R, gamma_ini, prior_gamma,
              gamma_ind_prop, nsample, burnin, thin, adjust)

Value

A list of class "NPP" with four elements:

acceptrate

the acceptance rate in MCMC sampling for \(\gamma\) using Metropolis-Hastings algorithm.

beta

posterior of the model parameter \(\beta\) in vector or matrix form.

sigma

posterior of the model parameter \(\sigma^2\).

delta

posterior of the power parameter \(\delta\).

Arguments

D0

a list of \(k\) elements representing \(k\) historical data, where the \(i^{th}\) element corresponds to the \(i^{th}\) historical data named as ``D0i''.

X

a vector or matrix or data frame of covariate observed in the current data. If more than 1 covariate available, the number of rows is equal to the number of observations.

Y

a vector of individual level of the response y in the current data.

a0

a positive shape parameter for inverse-gamma prior on model parameter \(\sigma^2\).

b

a positive scale parameter for inverse-gamma prior on model parameter \(\sigma^2\).

mu0

a vector of the mean for prior \(\beta|\sigma^2\).

R

a inverse matrix of the covariance matrix for prior \(\beta|\sigma^2\).

gamma_ini

the initial value of \(\gamma\) in MCMC sampling.

prior_gamma

a vector of the hyperparameters in the prior distribution \(Dirichlet(\alpha_1, \alpha_2, ... ,\alpha_K)\) for \(\gamma\).

gamma_ind_prop

a vector of the hyperparameters in the proposal distribution \(Dirichlet(\alpha_1, \alpha_2, ... ,\alpha_K)\) for \(\gamma\).

nsample

specifies the number of posterior samples in the output.

burnin

the number of burn-ins. The output will only show MCMC samples after bunrin.

thin

the thinning parameter in MCMC sampling.

adjust

Whether or not to adjust the parameters of the proposal distribution.

Author

Qiang Zhang

Details

The outputs include posteriors of the model parameters and power parameter, acceptance rate in sampling \(\gamma\). Let \(\theta\)=\((\beta, \sigma^2)\), the normalized power prior distribution is $$\frac{\pi_0(\gamma)\pi_0(\theta)\prod_{k=1}^{K}L(\theta|D_{0k})^(\sum_{i=1}^{k}\gamma_i)}{\int \pi_0(\theta)\prod_{k=1}^{K}L(\theta|D_{0k})^(\sum_{i=1}^{k}\gamma_i)\,d\theta}.$$

Here \(\pi_0(\gamma)\) and \(\pi_0(\theta)\) are the initial prior distributions of \(\gamma\) and \(\theta\), respectively. \(L(\theta|D_{0k})\) is the likelihood function of historical data \(D_{0k}\), and \(\sum_{i=1}^{k}\gamma_i\) is the corresponding power parameter.

References

Ibrahim, J.G., Chen, M.-H., Gwon, Y. and Chen, F. (2015). The Power Prior: Theory and Applications. Statistics in Medicine 34:3724-3749.

Duan, Y., Ye, K. and Smith, E.P. (2006). Evaluating Water Quality: Using Power Priors to Incorporate Historical Information. Environmetrics 17:95-106.

See Also

LMMNPP_MCMC1; LMMNPP_MCMC2; LMOMNPP_MCMC2

Examples

if (FALSE) {
set.seed(1234)
sigsq0 = 1

n01 = 100
theta01 = c(0, 1, 1)
X01 = cbind(1, rnorm(n01, mean=0, sd=1), runif(n01, min=-1, max=1))
Y01 = X01%*%as.vector(theta01) + rnorm(n01, mean=0, sd=sqrt(sigsq0))
D01 = cbind(X01, Y01)

n02 = 70
theta02 = c(0, 2, 3)
X02 = cbind(1, rnorm(n02, mean=0, sd=1), runif(n02, min=-1, max=1))
Y02 = X02%*%as.vector(theta02) + rnorm(n02, mean=0, sd=sqrt(sigsq0))
D02 = cbind(X02, Y02)

n03 = 50
theta03 = c(0, 3, 5)
X03 = cbind(1, rnorm(n03, mean=0, sd=1), runif(n03, min=-1, max=1))
Y03 = X03%*%as.vector(theta03) + rnorm(n03, mean=0, sd=sqrt(sigsq0))
D03 = cbind(X03, Y03)

D0 = list(D01, D02, D03)
n0 = c(n01, n02, n03)

n = 100
theta = c(0, 3, 5)
X = cbind(1, rnorm(n, mean=0, sd=1), runif(n, min=-1, max=1))
Y = X%*%as.vector(theta) + rnorm(n, mean=0, sd=sqrt(sigsq0))

LMOMNPP_MCMC1(D0=D0, X=X, Y=Y, a0=2, b=2, mu0=c(0,0,0), R=diag(c(1/64,1/64,1/64)),
              gamma_ini=NULL, prior_gamma=rep(1/4,4), gamma_ind_prop=rep(1/4,4),
              nsample=5000, burnin=1000, thin=5, adjust=FALSE)
}