bsar: Bayesian Shape-Restricted Spectral Analysis Regression

Description

This function fits a Bayesian semiparametric regression model to estimate shape-restricted functions using a spectral analysis of Gaussian process priors.

Usage

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)

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.
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 system.time().

Arguments

formula: an object of class “formula”
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): nblow0 (1000) giving the number of initialization period for adaptive metropolis, maxmodmet (5) giving the maximum number of times to modify metropolis, nblow (10000) giving the number of MCMC in transition period, nskip (10) giving the thinning interval, smcmc (1000) giving the number of MCMC for analysis, and ndisp (1000) giving the number of saved draws to be displayed on screen (the function reports on the screen when every ndisp iterations have been carried out).
prior: a list giving the prior information. The list includes the following parameters (default values specify the non-informative prior): iflagprior choosing a smoothing prior for spectral coefficients (iflagprior=0 assigns T-Smoother prior (default), iflagprior=1 chooses Lasso-Smoother prior), theta0_m0 and theta0_s0 giving the hyperparameters for prior distribution of the spectral coefficients (theta0_m0 and theta0_s0 are used when the functions have shape-restriction), tau2_m0, tau2_s0 and w0 giving the prior mean and standard deviation of smoothing prior (When iflagprior=1, tau2_m0 is only used as the hyperparameter), beta_m0 and beta_v0 giving the hyperparameters of the multivariate normal distribution for parametric part including intercept, sigma2_m0 and sigma2_v0 giving the prior mean and variance of the inverse gamma prior for the scale parameter of response, alpha_m0 and alpha_s0 giving the prior mean and standard deviation of the truncated normal prior distribution for the constant of integration, iflagpsi determining the prior of slope for logisitic function in S or U shaped (iflagpsi=1 (default), slope $\psi$ is sampled and iflagpsi=0, $\psi$ is fixed), psifixed giving initial value (iflagpsi=1) or fixed value (iflagpsi=0) of slope, omega_m0 and omega_s0 giving the prior mean and standard deviation of the truncated normal prior distribution for the inflection point of S or U shaped function.
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 TRUE, the iteration number and the Metropolis acceptance rate are printed to the screen.

Details

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\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 normal distribution, $N(0,\sigma^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 $\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 | \sigma, \tau, \gamma \sim N(0, \sigma^2\tau^2\exp[-j\gamma]), ~ j \ge 1$$

The priors for the spectral coefficients of shape restricted functions are: $$\theta_0 | \sigma \sim N(m_{\theta_0}, \sigma v^2_{\theta_0}), \quad \theta_j | \sigma, \tau, \gamma \sim N(m_{\theta_j}, \sigma\tau^2\exp[-j\gamma]), ~ j \ge 1$$

To complete the model specification, the conjugate priors are assumed for $\beta$ and $\sigma$: $$\beta | \sigma \sim N(m_{0,\beta}, \sigma^2V_{0,\beta}), \quad \sigma^2 \sim IG\left(\frac{r_{0,\sigma}}{2}, \frac{s_{0,\sigma}}{2}\right)$$

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.

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.

Examples

Run this code

# \dontshow{
  # Toy Example
# 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 <- 20

# Fit the model with default priors and mcmc parameters
fout <- bsar(y ~ fs(x), nbasis = nbasis, shape = 'Free')

# Summary
print(fout); summary(fout)

# Trace plots
plot(fout)

# fitted values
fit <- fitted(fout)

# Plot
plot(fit, ask = FALSE)
# }

if (FALSE) {

##########################################
# 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)
}

Run the code above in your browser using DataLab