bsar | R Documentation |
This function fits a Bayesian semiparametric regression model to estimate shape-restricted functions using a spectral analysis of Gaussian process priors.
bsar(formula, xmin, xmax, 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. |
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 regression model (Lenk and Choi, 2015) for estimating shape-restricted functions using Gaussian process priors. For enforcing shape-restrictions, they 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 a normal distribution, N(0,σ^2).
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 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.
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.
bsardpm
## Not run: ########################################## # Increasing Convex to Concave (S-shape) # ########################################## # simulate data f <- function(x) 5*exp(-10*(x - 1)^4) + 5*x^2 set.seed(1) n <- 100 x <- runif(n) y <- f(x) + rnorm(n, sd = 1) # Number of cosine basis functions nbasis <- 50 # Fit the model with default priors and mcmc parameters fout <- bsar(y ~ fs(x), nbasis = nbasis, shape = 'IncreasingConvex', spm.adequacy = TRUE) # Summary print(fout); summary(fout) # Trace plots plot(fout) # fitted values fit <- fitted(fout) # Plot plot(fit, ask = TRUE) ############################################# # Additive Model # # Monotone-Increasing and Increasing-Convex # ############################################# # Simulate data f1 <- function(x) 2*pi*x + sin(2*pi*x) f2 <- function(x) exp(6*x - 3) n <- 200 x1 <- runif(n) x2 <- runif(n) x <- cbind(x1, x2) y <- 5 + f1(x1) + f2(x2) + rnorm(n, sd = 0.5) # Number of cosine basis functions nbasis <- 50 # MCMC parameters mcmc <- list(nblow0 = 1000, nblow = 10000, nskip = 10, smcmc = 5000, ndisp = 1000, maxmodmet = 10) # Prior information xmin <- apply(x, 2, min) xmax <- apply(x, 2, max) xrange <- xmax - xmin prior <- list(iflagprior = 0, theta0_m0 = 0, theta0_s0 = 100, tau2_m0 = 1, tau2_v0 = 100, w0 = 2, beta_m0 = numeric(1), beta_v0 = diag(100,1), sigma2_m0 = 1, sigma2_v0 = 1000, alpha_m0 = 3, alpha_s0 = 50, iflagpsi = 1, psifixed = 1000, omega_m0 = (xmin + xmax)/2, omega_s0 = (xrange)/8) # Fit the model with user specific priors and mcmc parameters fout <- bsar(y ~ fs(x1) + fs(x2), nbasis = nbasis, mcmc = mcmc, prior = prior, shape = c('Increasing', 'IncreasingS')) # Summary print(fout); summary(fout) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.