gbsar {bsamGP} | R Documentation |
Bayesian Shape-Restricted Spectral Analysis for Generalized Partial Linear Models
Description
This function fits a Bayesian generalized partial linear regression model to estimate shape-restricted functions using a spectral analysis of Gaussian process priors.
Usage
gbsar(formula, xmin, xmax, family, link, nbasis, nint, mcmc = list(), prior = list(),
shape = c('Free','Increasing','Decreasing','IncreasingConvex','DecreasingConcave',
'IncreasingConcave','DecreasingConvex','IncreasingS','DecreasingS',
'IncreasingRotatedS','DecreasingRotatedS','InvertedU','Ushape'),
marginal.likelihood = TRUE, algorithm = c('AM', 'KS'), 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. |
family |
a description of the error distribution to be used in the model: The family contains bernoulli (“bernoulli”), poisson (“poisson”), negative-binomial (“negative.binomial”), poisson-gamma mixture (“poisson.gamma”). |
link |
a description of the link function to be used in the model. |
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. |
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. |
algorithm |
a description of the algorithm to be used in the fitting of the logistic model:
The algorithm contains the Gibbs sampler based on the Kolmogorov-Smirnov distribution ( |
verbose |
a logical variable. If |
Details
This generic function fits a Bayesian generalized partial linear regression models 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 and
be the response and the vector of parametric predictors, respectively.
Further, let
be the covariate related to the response through an unknown shape-restricted function.
The model for estimating shape-restricted functions is as follows.
where is a link function and
is an unknown nonlinear function of the scalar
.
The prior of function without shape restriction is:
where is a second-order Gaussian process with mean function equal to zero and covariance function
for
. The Gaussian process is expressed with
the spectral representation based on cosine basis functions:
The shape-restricted functions are modeled by assuming the th derivatives of
are squares of Gaussian processes:
where is the squish function. For monotonic, monotonic convex, and concave functions,
, while
for
S
and U
shaped functions, is defined by
For the spectral coefficients of functions without shape constraints, the following prior is used
(The intercept is included in ):
The priors for the spectral coefficients of shape restricted functions are:
To complete the model specification, the following prior is assumed for :
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.gd |
log marginal likelihood using Gelfand-Dey method. |
lmarg.nr |
log marginal likelihood using Netwon-Raftery method, which is biased. |
family |
the family object used. |
link |
the link object used. |
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.
Roberts, G. O. and Rosenthal, J. S. (2009) Examples of Adaptive MCMC. Journal of Computational and Graphical Statistics, 18, 349-367.
Holmes, C. C. and Held, L. (2006) Bayesian Auxiliary Variables Models for Binary and Multinomial Regression. Bayesian Analysis, 1, 145-168.
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.
Albert, J. H. and Chib, S. (1993) Bayesian Analysis of Binary and Polychotomous Response Data. Journal of the American Statistical Association, 88, 669-679.
See Also
Examples
## Not run:
###########################
# Probit Regression Model #
###########################
# Simulate data
set.seed(1)
f <- function(x) 1.5 * sin(pi * x)
n <- 1000
b <- c(1,-1)
rho <- 0.7
u <- runif(n, min = -1, max = 1)
x <- runif(n, min = -1, max = 1)
w1 <- runif(n, min = -1, max = 1)
w2 <- round(f(rho * x + (1 - rho) * u))
w <- cbind(w1, w2)
y <- w %*% b + f(x) + rnorm(n)
y <- (y > 0)
# Number of cosine basis functions
nbasis <- 50
# Fit the model with default priors and mcmc parameters
fout <- gbsar(y ~ w1 + w2 + fs(x), family = "bernoulli", link = "probit",
nbasis = nbasis, shape = 'Free')
# Summary
print(fout); summary(fout)
# fitted values
fit <- fitted(fout)
# Plot
plot(fit, ask = TRUE)
######################################
# Logistic Additive Regression Model #
######################################
# Wage-Union data
data(wage.union); attach(wage.union)
race[race==1 | race==2]=0
race[race==3]=1
y <- union
w <- cbind(race,sex,south)
x <- cbind(wage,education,age)
# mcmc parameters
mcmc <- list(nblow0 = 10000,
nblow = 10000,
nskip = 10,
smcmc = 1000,
ndisp = 1000,
maxmodmet = 10)
foutGBSAR <- gbsar(y ~ race + sex + south + fs(wage) + fs(education) + fs(age),
family = 'bernoulli', link = 'logit', nbasis = 50, mcmc = mcmc,
shape = c('Free','Decreasing','Increasing'))
# fitted values
fitGBSAR <- fitted(foutGBSAR)
# Plot
plot(fitGBSAR, ask = TRUE)
## End(Not run)