MCMC for normally distributed responses with additive model for the mean and variance functions achieved via spike-slab prior for variable selection. See `Details' section for a full description of the model.
mvrm(formula,formula.v=~1,data,sweeps,burn,thin=1,seed,StorageDir,
c.betaPrior="IG(0.5,0.5*n)",c.alphaPrior="IG(1.1,1.1)",
pi.muPrior="Beta(1,1)",pi.sigmaPrior="Beta(1,1)",sigmaPrior="HN(2)",...)a formula defining the response and the covariates in the mean model e.g. y ~ x or
for smooth effects y ~ sm(x).
a formula defining the variance model e.g. ~ x or
for smooth effects ~ sm(x).
a required data frame.
total number of posterior samples, including those discarded in burn-in period
(see argument burn) and those discarded by the thinning process
(see argument thin).
length of burn-in period.
thinning parameter.
optional seed for the random generator.
a required directory to store files with the posterior samples of models parameters.
The parameters of the inverse Gamma prior of \(c_{\beta}\). The default is "IG(0.5,0.5*n)", that is, an inverse Gamma with parameters \(1/2\) and \(n/2\), where \(n\) is the sample size.
The parameters of the inverse Gamma prior of \(c_{\alpha}\). The default is "IG(1.1,1.1)".
The parameters of the Beta prior of \(\pi_{\mu}\). The default is "Beta(1,1)".
The parameters of the Beta prior of \(\pi_{\sigma}\). The default is "Beta(1,1)".
The prior of \(\sigma^2\). The default is "HN(2)", a half-normal with variance equal to two, that is \(|\sigma| \sim N(0,2)\). Inverse Gamma prior is also available.
Other options that will be ignored.
Function mvrm returns the following:
the matched call.
mean model formula.
variance model formula.
the seed that was used (in case replication of the results is needed).
the dataset
the mean model design matrix.
the variance model design matrix.
the length of the vector of indicators \(\gamma\).
the length of the vector of indicators \(\delta\).
the MCMC parameters: length of burn in period, total number of samples, thinning period.
total number of posterior samples
the storage directory
contains samples from the posterior of vector \(\alpha\). Rows represent posterior samples and columns represent the regression coefficient, and they are in the same order as the columns of design matrix Z.
contains samples from the posterior of vector \(\beta\). Rows represent posterior samples and columns represent the regression coefficients, and they are in the same order as the columns of design matrix X.
contains samples from the posterior of the vector of the indicators \(\gamma\). Rows represent posterior samples and columns represent the indicator variables, and they are in the same order as the columns of design matrix X.
contains samples from the posterior of the vector of the indicators \(\delta\). Rows represent posterior samples and columns represent the indicator variables, and they are in the same order as the columns of design matrix Z.
contains samples from the posterior of the error variance \(\sigma^2\).
contains samples from the posterior of \(c_{\beta}\).
contains samples from the posterior of \(c_{\alpha}\).
Function mvrm returns samples from the posterior distributions of the parameters of
a regression model with normally distributed responses and mean and variance functions modeled
in terms of covariates. For instance, in the presence of two covariates in the mean model
(\(u_1, u_2\)) and two in the variance model (\(w_1, w_2\)), we may choose to fit
$$
\mu_u = \beta_0 + \beta_1 u_1 + f_{\mu}(u_2),$$
$$
\log(\sigma^2_W) = \alpha_0 + \alpha_1 w_1 + f_{\sigma}(w_2),$$
parametrically modelling the effects of \(u_1\) and \(w_1\) and non-parametrically
the effects of \(u_2\) and \(w_2\). Smooth functions, such as \(f_{\mu}\)
and \(f_{\sigma}\), are represented by basis function expansion. For instance
$$
f_{\mu}(u_2) = \sum_{j} \beta_{j} \phi_{j}(u_2),
$$
where \(\phi\) are the basis functions and \(\beta\) are the associated coefficients.
The variance model can equivalently be expressed as $$ \sigma^2_W = \exp(\alpha_0) \exp(\alpha_1 w_1 + f_{\sigma}(w_2)) = \sigma^2 \exp(\alpha_1 w_1 + f_{\sigma}(w_2)),$$ where \(\sigma^2 = \exp(\alpha_0)\). This is the parameterization that we adopt in this implementation.
Positive prior probability that the regression coefficients in the mean model are exactly zero is achieved by defining binary variables \(\gamma\) that take value \(\gamma=1\) if the associated coefficient \(\beta \neq 0\) and \(\gamma = 0\) if \(\beta = 0\). Indicators \(\delta\) that take value \(\delta=1\) if the associated coefficient \(\alpha \neq 0\) and \(\delta = 0\) if \(\alpha = 0\) for the variance function are defined analogously. We note that all coefficients in the mean and variance functions are subject to selection except the intercepts, \(\beta_0\) and \(\alpha_0\).
Prior specification:
For the vector of non-zero regression coefficients \(\beta_{\gamma}\) we specify a g-prior $$ \beta_{\gamma} | c_{\beta}, \sigma^2, \gamma, \alpha, \delta \sim N(0,c_{\beta} \sigma^2 (\tilde{X}_{\gamma}^{\top} \tilde{X}_{\gamma} )^{-1}). $$ where \(\tilde{X}\) is a scaled version of design matrix \(X\) of the mean model.
For the vector of non-zero regression coefficients \(\alpha_{\delta}\) we specify a normal prior $$ \alpha_{\delta} | c_{\alpha}, \delta \sim N(0,c_{\alpha} I). $$
Independent priors are specified for the indicators variables \(\gamma\) and \(\delta\) as \(P(\gamma = 1 | \pi_{\mu}) = \pi_{\mu}\) and \(P(\delta = 1 | \pi_{\sigma}) = \pi_{\sigma}\). Further, Beta priors are specified for \(\pi_{\mu}\) and \(\pi_{\sigma}\) $$ \pi_{\mu} \sim Beta(c_{\mu},d_{\mu}), \pi_{\sigma} \sim Beta(c_{\sigma},d_{\sigma}). $$ We note that blocks of regression coefficients associated with distinct covariate effects have their own probability of selection (\(\pi_{\mu}\) or \(\pi_{\sigma}\)) and this probability has its own prior distribution.
Further, we specify inverse Gamma priors for \(c_{\beta}\) and \(c_{\alpha}\) $$ c_{\beta} \sim IG(a_{\beta},b_{\beta}), c_{\alpha} \sim IG(a_{\alpha},b_{\alpha})$$
Lastly, for \(\sigma^2\) we consider inverse Gamma and half-normal priors $$ \sigma^2 \sim IG(a_{\sigma},b_{\sigma}), |\sigma| \sim N(0,\phi^2_{\sigma}). $$
Chan, D., Kohn, R., Nott, D., & Kirby, C. (2006). Locally adaptive semiparametric estimation of the mean and variance functions in regression models. Journal of Computational and Graphical Statistics, 15(4), 915-936.
# NOT RUN {
# Fit a mean/variance regression model on the cps71 dataset from package np
require(np)
require(ggplot2)
data(cps71)
meanModel<-logwage~sm(age,nknots=30,bs="rd")
varianceModel<-~sm(age,nknots=30,bs="rd")
DIR<-getwd()
m1<-mvrm(formula=meanModel,formula.v=varianceModel,data=cps71,sweeps=10,burn=5,thin=1,
seed=1,StorageDir=DIR)
# }
# NOT RUN {
m2 <- mvrm(formula=meanModel,formula.v=varianceModel,data=cps71,sweeps=10000,burn=5000,thin=2,
seed=1,StorageDir=DIR)
# }
# NOT RUN {
alpha<-mvrm2mcmc(mvrmObj=m1,"alpha")
summary(alpha)
plot(alpha)
wagePlotOptions<-list(geom_point(data=cps71,aes(x=age,y=logwage)))
plot(x=m1,model="mean",term="sm(age)",plotOptions=wagePlotOptions)
# }
Run the code above in your browser using DataLab