bsardpm {bsamGP} | R Documentation |
Bayesian Shape-Restricted Spectral Analysis Regression with Dirichlet Process Mixture Errors
Description
This function fits a Bayesian semiparametric regression model to estimate shape-restricted functions using a spectral analysis of Gaussian process priors. The model assumes that the errors follow a Dirichlet process mixture model.
Usage
bsardpm(formula, xmin, xmax, nbasis, nint,
mcmc = list(), prior = list(), egrid, ngrid, location = TRUE,
shape = c('Free', 'Increasing', 'Decreasing', 'IncreasingConvex', 'DecreasingConcave',
'IncreasingConcave', 'DecreasingConvex', 'IncreasingS', 'DecreasingS',
'IncreasingRotatedS','DecreasingRotatedS','InvertedU','Ushape'),
verbose = FALSE)
Arguments
formula |
an object of class “ |
xmin |
a vector or scalar giving user-specific minimum values of x. The default values are minimum values of x. |
xmax |
a vector or scalar giving user-specific maximum values of x. The default values are maximum values of x. |
nbasis |
number of cosine basis functions. |
nint |
number of grid points where the unknown function is evaluated for plotting. The default is 200. |
mcmc |
a list giving the MCMC parameters.
The list includes the following integers (with default values in parentheses):
|
prior |
a list giving the prior information. The list includes the following parameters
(default values specify the non-informative prior):
|
egrid |
a vector giving grid points where the residual density estimate is evaluated. The default range is from -10 to 10. |
ngrid |
a vector giving number of grid points where the residual density estimate is evaluated. The default value is 500. |
location |
a logical value. If it is true, error density is modelled using location-scale mixture. |
shape |
a vector giving types of shape restriction. |
verbose |
a logical variable. If |
Details
This generic function fits a Bayesian spectral analysis regression model for estimating shape-restricted functions using Gaussian process priors. For enforcing shape-restrictions, the model assumes that the derivatives of the functions are squares of Gaussian processes. The model also assumes that the errors follow a Dirichlet process mixture model.
Let y_i
and w_i
be the response and the vector of parametric predictors, respectively.
Further, let x_{i,k}
be the covariate related to the response through an unknown shape-restricted function.
The model for estimating shape-restricted functions is as follows.
y_i = w_i^T\beta + \sum_{k=1}^K f_k(x_{i,k}) + \epsilon_i, ~ i=1,\ldots,n,
where f_k
is an unknown shape-restricted function of the scalar x_{i,k} \in [0,1]
and
the error terms \{\epsilon_i\}
are a random sample from a Dirichlet process mixture model,
1. scale mixture :
\epsilon_i \sim f(\epsilon) = \int N(\epsilon; 0,\sigma^2)dG(\sigma^2),
G \sim DP(M,G0), ~~ G0 = Ga\left(\sigma^{-2}; \frac{r_{0,\sigma}}{2},\frac{s_{0,\sigma}}{2}\right).
2. location-scale mixture :
\epsilon_i \sim f(\epsilon) = \int N(\epsilon; \mu,\sigma^2)dG(\mu,\sigma^2),
G \sim DP(M,G0), ~~ G0 = N\left(\mu;\mu_0,\kappa\sigma^2\right)Ga\left(\sigma^{-2}; \frac{r_{0,\sigma}}{2},\frac{s_{0,\sigma}}{2}\right).
The prior of function without shape restriction is:
f(x) = Z(x),
where Z
is a second-order Gaussian process with mean function equal to zero and covariance function
\nu(s,t) = E[Z(s)Z(t)]
for s, t \in [0, 1]
. The Gaussian process is expressed with
the spectral representation based on cosine basis functions:
Z(x) = \sum_{j=0}^\infty \theta_j\varphi_j(x)
\varphi_0(x) = 1 ~~ \code{and} ~~ \varphi_j(x) = \sqrt{2}\cos(\pi j x), ~ j \ge 1, ~ 0 \le x \le 1
The shape-restricted functions are modeled by assuming the q
th derivatives of f
are squares of Gaussian processes:
f^{(q)}(x) = \delta Z^2(x)h(x), ~~ \delta \in \{1, -1\}, ~~ q \in \{1, 2\},
where h
is the squish function. For monotonic, monotonic convex, and concave functions, h(x)=1
, while
for S
and U
shaped functions, h
is defined by
h(x) = \frac{1 - \exp[\psi(x - \omega)]}{1 + \exp[\psi(x - \omega)]}, ~~ \psi > 0, ~~ 0 < \omega < 1
For the spectral coefficients of functions without shape constraints, the scale-invariant prior is used
(The intercept is included in \beta
):
\theta_j | \tau, \gamma \sim N(0, \tau^2\exp[-j\gamma]), ~ j \ge 1
The priors for the spectral coefficients of shape restricted functions are:
\theta_0 \sim N(m_{\theta_0}, v^2_{\theta_0}), \quad
\theta_j | \tau, \gamma \sim N(m_{\theta_j}, \tau^2\exp[-j\gamma]), ~ j \ge 1
To complete the model specification, the popular normal prior is assumed for \beta
:
\beta | \sim N(m_{0,\beta}, V_{0,\beta})
Value
An object of class bsam
representing the Bayesian spectral analysis model fit.
Generic functions such as print
, fitted
and plot
have methods to show the results of the fit.
The MCMC samples of the parameters in the model are stored in the list mcmc.draws
,
the posterior samples of the fitted values are stored in the list fit.draws
, and
the MCMC samples for the log marginal likelihood are saved in the list loglik.draws
.
The output list also includes the following objects:
post.est |
posterior estimates for all parameters in the model. |
lpml |
log pseudo marginal likelihood using Mukhopadhyay and Gelfand method. |
imodmet |
the number of times to modify Metropolis. |
pmet |
proportion of |
call |
the matched call. |
mcmctime |
running time of Markov chain from |
References
Jo, S., Choi, T., Park, B. and Lenk, P. (2019). bsamGP: An R Package for Bayesian Spectral Analysis Models Using Gaussian Process Priors. Journal of Statistical Software, 90, 310-320.
Lenk, P. and Choi, T. (2017) Bayesian Analysis of Shape-Restricted Functions using Gaussian Process Priors. Statistica Sinica, 27, 43-69.
MacEachern, S. N. and Müller, P. (1998) Estimating mixture of Dirichlet process models. Journal of Computational and Graphical Statistics, 7, 223-238.
Mukhopadhyay, S. and Gelfand, A. E. (1997) Dirichlet process mixed generalized linear models. Journal of the American Statistical Association, 92, 633-639.
Neal, R. M. (2000) Markov chain sampling methods for Dirichlet process mixture models. Journal of Computational and Graphical Statistics, 9, 249-265.
See Also
Examples
## Not run:
#####################
# Increasing-convex #
#####################
# Simulate data
set.seed(1)
n <- 200
x <- runif(n)
e <- c(rnorm(n/2, sd = 0.5), rnorm(n/2, sd = 3))
y <- exp(6*x - 3) + e
# Number of cosine basis functions
nbasis <- 50
# Fit the model with default priors and mcmc parameters
fout <- bsardpm(y ~ fs(x), nbasis = nbasis, shape = 'IncreasingConvex')
# Summary
print(fout); summary(fout)
# fitted values
fit <- fitted(fout)
# Plot
plot(fit, ask = TRUE)
## End(Not run)