bsaqdpm {bsamGP}  R Documentation 
This function fits a Bayesian semiparametric quantile regression model to estimate shaperestricted functions using a spectral analysis of Gaussian process priors. The model assumes that the errors follow a Dirichlet process mixture model.
bsaqdpm(formula, xmin, xmax, p, nbasis, nint,
mcmc = list(), prior = list(), egrid, ngrid = 500,
shape = c('Free', 'Increasing', 'Decreasing', 'IncreasingConvex', 'DecreasingConcave',
'IncreasingConcave', 'DecreasingConvex', 'IncreasingS', 'DecreasingS',
'IncreasingRotatedS', 'DecreasingRotatedS', 'InvertedU', 'Ushape'),
verbose = FALSE)
formula 
an object of class “ 
xmin 
a vector or scalar giving userspecific minimum values of x. The default values are minimum values of x. 
xmax 
a vector or scalar giving userspecific maximum values of x. The default values are maximum values of x. 
p 
quantile of interest (default=0.5). 
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 noninformative 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. 
shape 
a vector giving types of shape restriction. 
verbose 
a logical variable. If 
This generic function fits a Bayesian spectral analysis quantile regression model for estimating shaperestricted functions using Gaussian process priors. For enforcing shaperestrictions, 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 shaperestricted function.
The model for estimating shaperestricted 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 shaperestricted 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 of an asymmetric Laplace distribution,
ALD_p(0,\sigma^2)
, which has the following probability density function:
\epsilon_i \sim f(\epsilon) = \int ALD_p(\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).
The prior of function without shape restriction is:
f(x) = Z(x),
where Z
is a secondorder 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 shaperestricted 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 scaleinvariant 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})
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. 
rsquarey 
correlation between 
imodmet 
the number of times to modify Metropolis. 
pmet 
proportion of 
call 
the matched call. 
mcmctime 
running time of Markov chain from 
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, 310320.
Kozumi, H. and Kobayashi, G. (2011) Gibbs sampling methods for Bayesian quantile regression. Journal of Statistical Computation and Simulation, 81(11), 15651578.
Lenk, P. and Choi, T. (2017) Bayesian Analysis of ShapeRestricted Functions using Gaussian Process Priors. Statistica Sinica, 27, 4369.
MacEachern, S. N. and M\"uller, P. (1998) Estimating mixture of Dirichlet process models. Journal of Computational and Graphical Statistics, 7, 223238.
Mukhopadhyay, S. and Gelfand, A. E. (1997) Dirichlet process mixed generalized linear models. Journal of the American Statistical Association, 92, 633639.
Neal, R. M. (2000) Markov chain sampling methods for Dirichlet process mixture models. Journal of Computational and Graphical Statistics, 9, 249265.
## Not run:
######################
# Increasingconcave #
######################
# Simulate data
set.seed(1)
n < 500
x < runif(n)
e < c(rald(n/2, scale = 0.5, p = 0.5),
rald(n/2, scale = 3, p = 0.5))
y < log(1 + 10*x) + e
# Number of cosine basis functions
nbasis < 50
# Fit the model with default priors and mcmc parameters
fout1 < bsaqdpm(y ~ fs(x), p = 0.25, nbasis = nbasis,
shape = 'IncreasingConcave')
fout2 < bsaqdpm(y ~ fs(x), p = 0.5, nbasis = nbasis,
shape = 'IncreasingConcave')
fout3 < bsaqdpm(y ~ fs(x), p = 0.75, nbasis = nbasis,
shape = 'IncreasingConcave')
# fitted values
fit1 < fitted(fout1)
fit2 < fitted(fout2)
fit3 < fitted(fout3)
# plots
plot(x, y, lwd = 2, xlab = 'x', ylab = 'y')
lines(fit1$xgrid, fit1$wbeta$mean[1] + fit1$fxgrid$mean, lwd=2, col=2)
lines(fit2$xgrid, fit2$wbeta$mean[1] + fit2$fxgrid$mean, lwd=2, col=3)
lines(fit3$xgrid, fit3$wbeta$mean[1] + fit3$fxgrid$mean, lwd=2, col=4)
legend('topleft',legend=c('1st Quartile','2nd Quartile','3rd Quartile'),
lwd=2, col=2:4, lty=1)
## End(Not run)