bsaq | R Documentation |
This function fits a Bayesian semiparametric quantile regression model to estimate shape-restricted functions using a spectral analysis of Gaussian process priors.
bsaq(formula, xmin, xmax, p, nbasis, nint, mcmc = list(), prior = list(), shape = c('Free', 'Increasing', 'Decreasing', 'IncreasingConvex', 'DecreasingConcave', 'IncreasingConcave', 'DecreasingConvex', 'IncreasingS', 'DecreasingS', 'IncreasingRotatedS','DecreasingRotatedS','InvertedU','Ushape', 'IncMultExtreme','DecMultExtreme'), nExtreme = NULL, marginal.likelihood = TRUE, spm.adequacy = FALSE, verbose = FALSE)
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. |
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 non-informative prior):
|
shape |
a vector giving types of shape restriction. |
nExtreme |
a vector of extreme points for 'IncMultExtreme', 'DecMultExtreme' shape restrictions. |
marginal.likelihood |
a logical variable indicating whether the log marginal likelihood is calculated. The methods of Gelfand and Dey (1994) and Newton and Raftery (1994) are used. |
spm.adequacy |
a logical variable indicating whether the log marginal likelihood of linear model is calculated. The marginal likelihood gives the values of the linear regression model excluding the nonlinear parts. |
verbose |
a logical variable. If |
This generic function fits a Bayesian spectral analysis quantile regression model for estimating shape-restricted functions using Gaussian process priors. For enforcing shape-restrictions, the model assumed that the derivatives of the functions are squares of Gaussian processes.
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β + ∑_{k=1}^K f_k(x_{i,k}) + ε_i, ~ i=1,…,n,
where f_k is an unknown shape-restricted function of the scalar x_{i,k} \in [0,1] and the error terms \{ε_i\} are a random sample from an asymmetric Laplace distribution, ALD_p(0,σ^2), which has the following probability density function:
ALD_p(ε; μ, σ^2) = \frac{p(1-p)}{σ^2}\exp\Big(-\frac{(x-μ)[p - I(x ≤ μ)]}{σ^2}\Big),
where 0 < p < 1 is the skew parameter, σ^2 > 0 is the scale parameter, -∞ < μ < ∞ is the location parameter, and I(\cdot) is the indication function.
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 ν(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) = ∑_{j=0}^∞ θ_j\varphi_j(x)
\varphi_0(x) = 1 ~~ \code{and} ~~ \varphi_j(x) = √{2}\cos(π j x), ~ j ≥ 1, ~ 0 ≤ x ≤ 1
The shape-restricted functions are modeled by assuming the qth derivatives of f are squares of Gaussian processes:
f^{(q)}(x) = δ Z^2(x)h(x), ~~ δ \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[ψ(x - ω)]}{1 + \exp[ψ(x - ω)]}, ~~ ψ > 0, ~~ 0 < ω < 1
For the spectral coefficients of functions without shape constraints, the scale-invariant prior is used (The intercept is included in β):
θ_j | σ, τ, γ \sim N(0, σ^2τ^2\exp[-jγ]), ~ j ≥ 1
The priors for the spectral coefficients of shape restricted functions are:
θ_0 | σ \sim N(m_{θ_0}, σ v^2_{θ_0}), \quad θ_j | σ, τ, γ \sim N(m_{θ_j}, στ^2\exp[-jγ]), ~ j ≥ 1
To complete the model specification, the conjugate priors are assumed for β and σ:
β | σ \sim N(m_{0,β}, σ^2V_{0,β}), \quad σ^2 \sim IG≤ft(\frac{r_{0,σ}}{2}, \frac{s_{0,σ}}{2}\right)
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. |
lmarg.lm |
log marginal likelihood for linear quantile regression model. |
lmarg.gd |
log marginal likelihood using Gelfand-Dey method. |
lmarg.nr |
log marginal likelihood using Netwon-Raftery method, which is biased. |
rsquarey |
correlation between y and \hat{y}. |
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, 310-320.
Lenk, P. and Choi, T. (2017) Bayesian Analysis of Shape-Restricted Functions using Gaussian Process Priors. Statistica Sinica, 27: 43-69.
Gelfand, A. E. and Dey, K. K. (1994) Bayesian model choice: asymptotics and exact calculations. Journal of the Royal Statistical Society. Series B - Statistical Methodology, 56, 501-514.
Kozumi, H. and Kobayashi, G. (2011) Gibbs sampling methods for Bayesian quantile regression. Journal of Statistical Computation and Simulation, 81(11), 1565-1578.
Newton, M. A. and Raftery, A. E. (1994) Approximate Bayesian inference with the weighted likelihood bootstrap (with discussion). Journal of the Royal Statistical Society. Series B - Statistical Methodology, 56, 3-48.
bsar
, gbsar
## Not run: ###################### # Increasing-concave # ###################### # Simulate data set.seed(1) n <- 200 x <- runif(n) y <- log(1 + 10*x) + rald(n, scale = 0.5, p = 0.5) # Number of cosine basis functions nbasis <- 50 # Fit the model with default priors and mcmc parameters fout1 <- bsaq(y ~ fs(x), p = 0.25, nbasis = nbasis, shape = 'IncreasingConcave') fout2 <- bsaq(y ~ fs(x), p = 0.5, nbasis = nbasis, shape = 'IncreasingConcave') fout3 <- bsaq(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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.