Priors on model space for variable selection problems

Description

unifPrior implements a uniform prior (equal a priori probability for all models). binomPrior implements a Binomial prior. bbPrior implements a Beta-Binomial prior.

Usage

unifPrior(sel, logscale=TRUE, groups=1:length(sel),
constraints=lapply(1:length(unique(groups)), function(z) integer(0)))
binomPrior(sel, prob=.5, logscale=TRUE, probconstr=prob, groups=1:length(sel),
constraints=lapply(1:length(unique(groups)), function(z) integer(0)))
bbPrior(sel, alpha=1, beta=1, logscale=TRUE, alphaconstr=alpha,
betaconstr=beta, groups=1:length(sel),
constraints=lapply(1:length(unique(groups)), function(z) integer(0)))

Arguments

Value

Prior probability of the specified model

Author(s)

David Rossell

Examples

sel <- c(TRUE,TRUE,FALSE,FALSE)
unifPrior(sel,logscale=FALSE)
binomPrior(sel,prob=.5,logscale=FALSE)
bbPrior(sel,alpha=1,beta=1,logscale=FALSE)

Model with best AIC, BIC, EBIC or other general information criteria (getIC)

Description

Search for the regression model attaining the best value of the specified information criterion

Usage

 bestAIC(...)
 bestAIC_fast(...)
 bestBIC(...)
 bestBIC_fast(...)
 bestEBIC(...)
 bestEBIC_fast(...)
 bestIC(..., penalty)
 bestIC_fast(..., penalty)

Arguments

Details

bestAIC, bestBIC, bestEBIC, bestIC perform full model enumeration when possible and otherwise resort to MCMC to explore the models, as discussed in function modelSelection. bestAIC_fast, bestBIC_fast, bestEBIC_fast, bestIC_fast use a faster algorithm. For Gaussian and binary outcomes it uses function L0Learn.fit from package L0Learn (Hazimeh et al, 2023), which combines coordinate descent with local combinatorial search. For any other outcome type it uses only coordinte descent.

bestBIC and the other functions codumented here take similar arguments to those of modelSelection, the primary difference is that no priors on models or parameters are needed.

Let p be the total number of parameters and n the sample size. The BIC of a model k with p_k parameters is

- 2 L_k + p_k log(n)

the AIC is

- 2 L_k + p_k 2

the EBIC is

- 2 L_k + p_k log(n) + 2 log(p choose p_k)

and a general information criterion with a given model size penalty

- 2 L_k + p_k penalty

The MCMC model search is based on assigning a probability to each model, and then using MCMC to sample models from this distribution. The probability of model k is

exp(- IC_k / 2) / sum_l exp(- IC_l / 2)

where IC_k is the value of the information criterion (BIC, EBIC...)

Hence the model with best (lowest) IC_k has highest probability, which means that it is likely to be sampled by the MCMC algorithm.

Value

Object of class icfit. Use (coef, summary, confint, predict) to get inference for the top model, and help(icfit-class) for more details on the returned object.

Author(s)

David Rossell

References

H. Hazimeh, R. Mazumder, T. Nonet. L0learn: A scalable package for sparse learning using l0 regularization. Journal of Machine Learning Research 24.205 (2023): 1-8.

See Also

modelSelection to perform Bayesian model selection, and for a description of the arguments that should be passed in ...

Examples

x <- matrix(rnorm(100*3),nrow=100,ncol=3)
theta <- matrix(c(1,1,0),ncol=1)
y <- x %*% theta + rnorm(100)
ybin <- y>0
df <- data.frame(y, ybin, x)
#BIC for all models (the intercept is also selected in/out)
fit= bestBIC(y ~ X1 + X2, data=df)
fit
#Same, but setting the BIC's log(n) penalty manually
#change the penalty for other General Info Criteria
#n= nrow(x)
#fit= bestIC(y ~ X1 + X2, data=df, penalty=log(n))
summary(fit) #usual GLM summary
coef(fit) #MLE under top model
#confint(fit) #conf int under top model (requires MASS package)
#Binary outcome
fit2= bestBIC(ybin ~ X1 + X2, data=df, family='binomial')
fit2

Number of Normal mixture components under Normal-IW and Non-local priors

Description

Posterior sampling and Bayesian model selection to choose the number of components k in multivariate Normal mixtures.

bfnormmix computes posterior probabilities under non-local MOM-IW-Dir(q) priors, and also for local Normal-IW-Dir(q.niw) priors. It also computes posterior probabilities on cluster occupancy and posterior samples on the model parameters for several k.

Usage

bfnormmix(x, k=1:2, mu0=rep(0,ncol(x)), g, nu0, S0, q=3, q.niw=1,
B=10^4, burnin= round(B/10), logscale=TRUE, returndraws=TRUE, verbose=TRUE)

Arguments

Details

The likelihood is

p(x[i,] | mu,Sigma,eta)= sum_j eta_j N(x[i,]; mu_j,Sigma_j)

The Normal-IW-Dir prior is

Dir(eta; q.niw) prod_j N(mu_j; mu0, g Sigma) IW(Sigma_j; nu0, S0)

The MOM-IW-Dir prior is

d(\mu,A) Dir(\eta; q) \prod_j N(\mu_j; \mu0, g \Sigma_j) IW(\Sigma_j; \nu_0, S0)

where

d(\mu,A)= [\prod_{j<l} (\mu_j-\mu_l)' A (\mu_j-\mu_l)]

and A is the average of \Sigma_1^{-1},...,\Sigma_k^{-1}. Note that one must have q>1 for the MOM-IW-Dir to define a non-local prior.

By default the prior parameter g is set such that

P( (mu[j]-mu[l])' A (mu[j]-mu[l]) < 4)= 0.05.

The reasonale when Sigma[j]=Sigma[l] and eta[j]=eta[l] then (mu[j]-mu[l])' A (mu[j]-mu[l])>4 corresponds to a bimodal density. That is, the default g focuses 0.95 prior prob on a degree of separation between components giving rise to a bimodal mixture density.

bfnormmix computes posterior model probabilities under the MOM-IW-Dir and Normal-IW-Dir priors using MCMC output. As described in Fuquene, Steel and Rossell (2018) the estimate is based on the posterior probability that one cluster is empty under each possible k.

Value

A list with elements

Author(s)

David Rossell

References

Fuquene J., Steel M.F.J., Rossell D. On choosing mixture components via non-local priors. 2018. arXiv

Examples

x <- matrix(rnorm(100*2),ncol=2)
bfnormmix(x=x,k=1:3)

Treatment effect estimation for linear models via Confounder Importance Learning using non-local priors.

Description

Treatment effect estimation for linear models in the presence of multiple treatments and a potentially high-dimensional number of controls, i.e. p \gg n can be handled.

Confounder Importance Learning (CIL) proposes an estimation framework where the importance of the relationship between treatments and controls is factored in into the establishment of prior inclusion probabilities for each of these controls on the response model. This is combined with the use of non-local priors to obtain BMA estimates and posterior model probabilities.

cil is built on modelSelection and produces objects of type cilfit. Use coef and postProb to obtain treatment effect point estimates and posterior model probabilities, respectively, on this object class.

Usage

cil(y, D, X, I = NULL, family = 'normal', familyD = 'normal',
 R = 1e4, Rinit = 500, th.search = 'EB', mod1 = 'lasso_bic',
 th.prior = 'unif', priorCoef = momprior(taustd=1),
 rho.min = NULL, rho.max = 0.95,
 th.range = NULL, max.mod = 2^20, lpen = 'lambda.1se',
 eps = 1e-10, bvs.fit0 = NULL, th.EP = NULL, center = TRUE, scale =
 TRUE, includevars, verbose = TRUE)

Arguments

Details

We estimate treatment effects for the features present in the treatment matrix D. Features in X, which may or may not be causal factors of the treatments of interest, only act as controls and, therefore, are not used as inferential subjects.

Confounder importance learning learns from data the amount of confounding in the data, based on a so-called confounding coefficient, and uses this to set covariate prior inclusion probabilities. Roughly, the coefficient measures to what extend covariates that are truly related to the outcome are also truly related (high confounding) or not related (no confounding) to the treatment(s).

In high confounding, covariates in X that appear to be related to D are assigned high inclusion probability in the regression for y. In low confounding, they're assigned low prior inclusion probability. See function plotprior to visualize these prior probabilities. Prior probabilities are regulated through a hyper-parameter theta that is set according to the method supplied to th.search. While the EB option obtains a more precise estimate, particularly when there are many covariates in X, the EP is much faster computationally and typically gives very similar results.

See references for further details on implementation and computation.

Value

Object of class cilfit, which extends a list with elements

Author(s)

Miquel Torrens

References

Torrens i Dinares M., Papaspiliopoulos O., Rossell D. Confounder importance learning for treatment effect inference. https://arxiv.org/abs/2110.00314, 2021, 1–48.

See Also

postProb to obtain posterior model probabilities.

coef for inference on the treatment parameters.

plotprior to plot the estimated prior probabilities, as a function of

Examples

# Simulate data
set.seed(1)
X <- matrix(rnorm(100 * 50), nrow = 100, ncol = 50)
beta_y <- matrix(c(rep(1, 6), rep(0, 44)), ncol = 1)
beta_d <- matrix(c(rep(1, 6), rep(0, 44)), ncol = 1)
alpha <- 1
d <- X %*% beta_d + rnorm(100)
y <- d * alpha + X %*% beta_y + rnorm(100)
# Confounder Importance Learning
fit1 <- cil(y = y, D = d, X = X, th.search = 'EP')
# BMA for treatment effects
coef(fit1)
# BMA for all covariates
head(fit1$coef)
# Estimated prior inclusion prob
# vs. treatment regression coefficients
plotprior(fit1)

Density and random draws from the asymmetric Laplace distribution

Description

dalapl evaluates the probability density function, palapl the cumulative probability function and ralapl generates random draws.

Usage

dalapl(x, th=0, scale=1, alpha=0, logscale=FALSE)
palapl(x, th=0, scale=1, alpha=0)
ralapl(n, th=0, scale=1, alpha=0)

Arguments

Details

For x<=th the asymmetric Laplace pdf is

0.5*exp(-abs(th-x)/(sqrt(scale)*(1+alpha)))/sqrt(scale)

and for x>th it is

0.5*exp(-abs(th-x)/(sqrt(scale)*(1-alpha)))/sqrt(scale)

Value

dalapl returns the density function, palapl the cumulative probability, ralapl random draws.

Author(s)

David Rossell

Examples

e <- ralapl(n=10^4, th=1, scale=2, alpha=0.5)
thseq <- seq(min(e),max(e),length=1000)
hist(e, main='', breaks=30, prob=TRUE)
lines(thseq, dalapl(thseq, th=1, scale=2, alpha=0.5), col=2)

Dirichlet density

Description

Evaluate the density of a Dirichlet distribution

Usage

ddir(x, q, logscale=TRUE)

Arguments

Value

Density of a Dirichlet(q) distribution evaluated at each row of x

Author(s)

David Rossell

Examples

x= matrix(c(1/3,2/3,.5,.5),nrow=2,byrow=TRUE)
ddir(x,q=2)

Density for Inverse Wishart distribution

Description

diwish returns the density for the inverse Wishart(nu,S) evaluated at Sigma.

Usage

diwish(Sigma, nu, S, logscale=FALSE)

Arguments

Value

Inverse Wishart(nu,S) density evaluated at Sigma

Author(s)

David Rossell

See Also

dpostNIW for the Normal-IW posterior density

Examples

Sigma= matrix(c(2,1,1,2),nrow=2)
diwish(Sigma,nu=4,S=diag(2))

Non-local prior density, cdf and quantile functions.

Description

dmom, dimom and demom return the density for the moment, inverse moment and exponential moment priors. pmom, pimom and pemom return the distribution function for the univariate moment, inverse moment and exponential moment priors (respectively). qmom and qimom return the quantiles for the univariate moment and inverse moment priors. dmomigmarg returns the marginal density implied by a MOM(x;tau*phi)*Invgamma(phi;a/2,b/2), pmomigmarg its cdf. Analogously demomigmarg and demomigmarg for eMOM(x;tau*phi)*Invgamma(phi;a/2,b/2)

Usage

dmom(x, tau, a.tau, b.tau, phi=1, r=1, V1, baseDensity='normal', nu=3,
logscale=FALSE, penalty='product')
dimom(x, tau=1, phi=1, V1, logscale=FALSE, penalty='product')
demom(x, tau, a.tau, b.tau, phi=1, logscale=FALSE)
pmom(q, V1 = 1, tau = 1)
pimom(q, V1 = 1, tau = 1, nu = 1)
pemom(q, tau, a.tau, b.tau)
qmom(p, V1 = 1, tau = 1)
qimom(p, V1 = 1, tau = 1, nu = 1)
dmomigmarg(x,tau,a,b,logscale=FALSE)
pmomigmarg(x,tau,a,b)
demomigmarg(x,tau,a,b,logscale=FALSE)
pemomigmarg(x,tau,a,b)

Arguments

Details

For type=='quadratic' the density is as follows. Define the quadratic form q(theta)= (theta-theta0)' * solve(V1) * (theta-theta0) / (tau*phi). The normal moment prior density is proportional to q(theta)*dmvnorm(theta,theta0,tau*phi*V1). The T moment prior is proportional to q(theta)*dmvt(theta,theta0,tau*phi*V1,df=nu). The inverse moment prior density is proportional to q(theta)^(-(nu+d)/2) * exp(-1/q(theta)).

pmom, pimom and qimom use closed-form expressions, while qmom uses nlminb to find quantiles numerically. Only the univariate version is implemented. In this case the product MOM is equivalent to the quadratic MOM. The same happens for the iMOM.

dmomigmarg returns the marginal density

p(x)= int MOM(x;0,tau*phi) IG(phi;a/2,b/2) dphi

Value

Prior density, cumulative distribution function or quantile.

Author(s)

David Rossell

References

Johnson V.E., Rossell D. Non-Local Prior Densities for Default Bayesian Hypothesis Tests. Journal of the Royal Statistical Society B, 2010, 72, 143-170.

Johnson V.E., Rossell D. Bayesian model selection in high-dimensional settings. Journal of the American Statistical Assocation, 2012, 107, 649-660

See http://rosselldavid.googlepages.com for technical reports.

Examples

#evaluate and plot moment and inverse moment prior densities
tau <- 1
thseq <- seq(-3, 3, length=500)
plot(thseq, dmom(thseq,tau=tau), type='l', ylab='Prior density')
lines(thseq, dimom(thseq,tau=tau), lty=2, col=2)

Posterior Normal-IWishart density

Description

dpostNIW evalutes the posterior Normal-IWishart density at (mu,Sigma). rpostNIW draws independent samples. This posterior corresponds to a Normal model for the data

x[i,] ~ N(mu,Sigma) iid i=1,...,n

under conjugate priors

mu | Sigma ~ N(mu0, g Sigma) Sigma ~ IW(nu0, S0)

Usage

dpostNIW(mu, Sigma, x, g=1, mu0=rep(0,length(mu)), nu0=nrow(Sigma)+1, S0,
 logscale=FALSE)
rpostNIW(n, x, g=1, mu0=0, nu0, S0, precision=FALSE)

Arguments

Value

dpostNIW returns the Normal-IW posterior density evaluated at (mu,Sigma).

rpostNIW returns a list with two elements. The first element are posterior draws for the mean. The second element are posterior draws for the covariance (or its inverse if precision==TRUE). Only lower-diagonal elements are returned (Sigma[lower.tri(Sigma,diag=TRUE)]).

Author(s)

David Rossell

See Also

diwish for the inverse Wishart prior density, marginalNIW for the integrated likelihood under a Normal-IW prior

Examples

#Simulate data
x= matrix(rnorm(100),ncol=2)
#Evaluate posterior at data-generating truth
mu= c(0,0)
Sigma= diag(2)
dpostNIW(mu,Sigma,x=x,g=1,nu0=4,log=FALSE)

Expectation of a product of powers of Normal or T random variables

Description

Compute the mean of prod(x)^power when x follows T_dof(mu,sigma) distribution (dof= -1 for multivariate Normal).

Usage

eprod(m, S, power = 1, dof = -1) 

Arguments

Details

The calculation is based on the computationally efficient approach by Kan (2008).

Value

Expectation of the above-mentioned product

Author(s)

John Cook

References

Kan R. From moments of sum to moments of product. Journal of Multivariate Analysis 99 (2008), 542-554.

Examples

#Check easy independence case
m <- c(0,3); S <- matrix(c(2,0,0,1),ncol=2)
eprod(m, S, power=2)
(m[1]^2+S[1][1])*(m[2]^2+S[2][2])

Class "icfit"

Description

Stores the output of the search for the model with best information criterion value, e.g. produced by bestBIC, bestBIC, bestAIC or bestIC. The class extends a list, so all usual methods for lists also work for icfit objects, e.g. accessing elements, retrieving names etc.

Methods are provided to extract coefficients, predictions, confidence intervals and summary information about the best model.

Objects from the Class

icfit objects are automatically created by a call to bestBIC or similar.

Slots

The class extends a list with elements:

topmodel

names of the variables in the top model

topmodel.fit

top model as fitted by glm

models

data frame with the information criterion for all models (when enumeration is feasible) or those visited by an MCMC model search in modelSelection (when enumeration is not feasible)

varnames

the names of all variables in the design matrix

msfit

Output of modelSelection (used to search the top model)

Methods

coef

glm fit for the top model

confint

Confidence intervals under the top model

predict

Predictions for the top model (predict.glm)

show

signature(object = "icfit"): Displays general information about the object.

Author(s)

David Rossell

See Also

See also bestBIC .

Examples

showClass("icfit")

Extract estimated inverse covariance

Description

Extract the estimated inverse covariance from an msfit_ggm object.

Bayesian model averaging is used, optionally entries with posterior probability of being non-zero below a threshold are set to 0.

Usage

icov(fit, threshold) 

Arguments

Details

The inverse covariance is obtained via Bayesian model averaging, using posterior samples of Omega. When threshold is specified, entries in the BMA estimate are set to zero, which may result in a non positive-definite matrix.

Value

Estimated inverse covariance matrix.

Author(s)

David Rossell

See Also

modelSelectionGGM , coef.msfit_ggm for Bayesian model averaging estimates and intervals.

Examples

## See modelSelectionGGM

Local variable selection

Description

Learn whether covariate effects are zero at given coordinates using Bayesian model selection or information criteria.

Use coef to extract estimates and posterior probabilities for local effects.

Usage

localnulltest(y, x, z, x.adjust, localgridsize, localgrid,
nbaseknots=20, nlocalknots=c(5,10,15), localknots,
basedegree=3, cutdegree=0,
usecutbasis=TRUE, priorCoef=normalidprior(),
priorGroup=priorCoef, priorModel=modelbbprior(),
mc.cores=min(4,length(nlocalknots)), return.mcmc=FALSE, verbose=FALSE,
...)
localnulltest_fda(y, x, z, x.adjust, function_id,
Sigma='AR/MA', localgridsize, localgrid, nbaseknots=20,
nlocalknots=c(5,10,15), localknots,
basedegree=3, cutdegree=0, usecutbasis=TRUE,
priorCoef=momprior(), priorGroup=groupmomprior(),
priorModel=modelbbprior(), mc.cores=min(4,length(nlocalknots)),
return.mcmc=FALSE, verbose=FALSE, ...)
localnulltest_givenknots(y, x, z, x.adjust, localgridsize,
localgrid, nbaseknots=20, nlocalknots=10, localknots,
basedegree=3, cutdegree=0,
usecutbasis=TRUE, priorCoef=normalidprior(),
priorGroup=priorCoef, priorModel=modelbbprior(),
verbose=FALSE, ...)
localnulltest_fda_givenknots(y, x, z, x.adjust, function_id,
Sigma='AR/MA', localgridsize, localgrid, nbaseknots=20,
nlocalknots=10, localknots,
basedegree=3, cutdegree=0, usecutbasis=TRUE,
priorCoef=momprior(), priorGroup=groupmomprior(),
priorModel=modelbbprior(), verbose=FALSE, ...)

Arguments

Details

Local variable selection considers the model

y_i= \beta_0(z_i) + sum_{j=1}^p \beta_j(z_i, x_i) + e_i

\beta_0(z_i) is the baseline mean

\beta_j(z_i,x_i) is local effect of covariate j at coordinate z_i

e_i a Gaussian error term assumed either independent or with a covariance structure given by Sigma. If assuming independence it is possible to consider alternatives to Gaussianity, e.g. set family='binomial' for logistic regression or family='poisson' for Poisson regression

Note: a sum-to-zero type constraint is set on \beta_1(z_i,x_i) so that it defines a deviation from the baseline mean \beta_0(z_i)

We model \beta_0 using B-splines of degree basedegree with nbaseknots knots. We model \beta_j using B-splines of degree cutdegree with nlocalknots. Using cutdegree=0 runs fastest is usually gives similar inference than higher degrees, and is hence recommended by default.

Value

Object of class localtest, which extends a list with elements

Author(s)

David Rossell

Examples

#Simulate outcome and 2 covariates
#Covariate 1 has local effect for z>0
#Covariate 2 has no effect for any z
truemean= function(x,z) {
 ans= double(nrow(x))
 group1= (x[,1]==1)
 ans[group1]= ifelse(z[group1] <=0, cos(z[group1]), 1)
 ans[!group1]= ifelse(z[!group1]<=0, cos(z[!group1]), 1/(z[!group1]+1)^2)
 return(ans)
}
n= 1000
x1= rep(0:1,c(n/2,n/2))
x2= x1 + rnorm(n)
x= cbind(x1,x2)
z= runif(n,-3,3)
m= truemean(x,z)
y= truemean(x,z) + rnorm(n, 0, .5)
#Run localnulltest with 10 knots
fit0= localnulltest(y, x=x, z=z, nlocalknots=10, niter=1000)
#Estimated covariate effects and posterior probabilities
b= coef(fit0)
b

Marginal (or integrated) likelihood density of the observed data for an individual model handled by modelSelection (regression, GLM, GAM, accelerated failure time, regression with Normal, two-piece Normal, Laplace or two-piece Laplace residuals

Description

The marginal density of the data, i.e. the likelihood integrated with respect to the prior distribution on the regression coefficients of the variables included in the model.

Usage

marginalLikelihood(sel, y, x, data, smoothterms, nknots=9, groups=1:ncol(x),
family="normal", priorCoef, priorGroup,
priorVar=igprior(alpha=0.01,lambda=0.01), priorSkew=momprior(tau=0.348),
neighbours, 
phi, method='auto', adj.overdisp='intercept', hess='asymp', optimMethod,
optim_maxit, initpar='none', B=10^5, logscale=TRUE, XtX, ytX)

Arguments

Details

The marginal density of the data y under a given model is

p(y | model) int p(y | theta) d P(theta | model)

where P(theta | model) is the prior distribution on the parameters included by the model.

Value

Marginal (or integrated) likelihood of the data under the specified prior.

Author(s)

David Rossell

See Also

modelSelection to perform model selection based on product non-local priors.

Examples

x <- matrix(rnorm(100*2),ncol=2)
y <- x %*% matrix(c(.5,1),ncol=1) + rnorm(nrow(x))
#Marginal likelihood for 2 models under pMOM prior
marginalLikelihood(c(TRUE,FALSE), y=y, x=x, priorCoef=momprior())
marginalLikelihood(c(TRUE, TRUE), y=y, x=x, priorCoef=momprior())
#Same, under Normal prior with diagonal covariance
marginalLikelihood(c(TRUE,FALSE), y=y, x=x, priorCoef=normalidprior())
marginalLikelihood(c(TRUE, TRUE), y=y, x=x, priorCoef=normalidprior())

Marginal likelihood under a multivariate Normal likelihood and a conjugate Normal-inverse Wishart prior.

Description

The argument z can be used to specify cluster allocations. If left missing then the usual marginal likelihood is computed, else it is computed conditional on the clusters (this is equivalent to the product of marginal likelihoods across clusters)

Usage

marginalNIW(x, xbar, samplecov, n, z, g, mu0=rep(0,ncol(x)),
nu0=ncol(x)+4, S0, logscale=TRUE)

Arguments

Details

The function computes

p(x)= int p(x | mu,Sigma) p(mu,Sigma) dmu dSigma

where p(x[i])= N(x[i]; mu,Sigma) iid i=1,...,n

p(mu | Sigma)= N(mu; mu0, g Sigma) p(Sigma)= IW(Sigma; nu0, S0)

Value

If z is missing the integrated likelihood under a Normal-IW prior. If z was specified then the product of integrated likelihoods across clusters

Author(s)

David Rossell

See Also

dpostNIW for the posterior Normal-IW density.

Examples

#Simulate data
x= matrix(rnorm(100),ncol=2)
#Integrated likelihood under correct model
marginalNIW(x,g=1,nu0=4,log=FALSE)
#Integrated likelihood under random cluster allocations
z= rep(1:2,each=25)
marginalNIW(x,z=z,g=1,nu0=4,log=FALSE)

Class "mixturebf"

Description

Stores the output of Bayesian model selection for mixture models, e.g. as produced by function bfnormmix.

Methods are provided for retrieving the posterior probability of a given number of mixture components, posterior means and posterior samples of the mixture model parameters.

Objects from the Class

Typically objects are automatically created by a call to bfnormmix.

Slots

The class has the following slots:

postprob

data.frame containing posterior probabilities for different numbers of components (k) and log-posterior probability of a component being empty (contain no individuals)

p

Number of variables in the data to which the model was fit

n

Number of observations in the data to which the model was fit

priorpars

Prior parameters used when fitting the model

postpars

Posterior parameters for a 1-component mixture, e.g. for a Normal mixture the posterior is N(mu1,Sigma/prec) IW(nu1,S1)

mcmc

For each considered value of k, posterior samples for the parameters of the k-component model are stored

Methods

coef

Computes posterior means for all parameters

show

signature(object = "mixturebf"): Displays general information about the object.

postProb

signature(object = "mixturebf"): Extracts posterior model probabilities, Bayes factors and posterior probability of a cluster being empty

postSamples

signature(object = "mixturebf"): Extracts posterior samples

Author(s)

David Rossell

References

Fuquene J., Steel M.F.J., Rossell D. On choosing mixture components via non-local priors. 2018. arXiv

See Also

See also bfnormmix

Examples

showClass("mixturebf")

Bayesian variable selection for generalized linear and generalized additive models.

Description

Bayesian model selection for linear, asymmetric linear, median and quantile regression, for GLMs and GAMs under non-local or Zellner priors. p>>n can be handled.

modelSelection enumerates all models when feasible and uses a Gibbs scheme otherwise. See coef and coefByModel for estimates and posterior intervals of regression coefficients, and rnlp for posterior samples.

modelSelection performs standard Bayesian model selection.

modelSelection_eBayes uses empirical Bayes to set prior inclusion probabilities pi. pi can be regressed on meta-covariates Z via log(pi/(1-pi))= Z w. If Z contains only an intercept, a global prior inclusion probability is learn from data.

modelsearchBlockDiag seeks the highest posterior probability model using an iterative block search.

Usage

modelSelection(y, x, data, smoothterms, nknots=9,
groups, constraints, center=TRUE, scale=TRUE,
enumerate, includevars, models,
maxvars, niter=5000, thinning=1,
burnin=round(niter/10), family='normal', priorCoef,
priorGroup, priorModel=modelbbprior(1,1),
priorConstraints,
priorVar=igprior(.01,.01),
priorSkew=momprior(tau=0.348),
neighbours, phi, deltaini,
initSearch='CDA', method='auto', adj.overdisp='intercept',
hess='asymp', optimMethod, optim_maxit, initpar='none', B=10^5,
XtXprecomp, verbose=TRUE)
modelSelection_eBayes(Z, wini, niter.mcmc= 5000, niter.mstep= 1000,
niter.eBayes= 20, priorvar.w, verbose=TRUE, ...)
modelsearchBlockDiag(y, x, priorCoef=momprior(tau=0.348),
priorModel=modelbbprior(1,1), priorVar=igprior(0.01,0.01),
blocksize=10, maxiter=10, maxvars=100, maxlogmargdrop=20,
maxenum=10, verbose=TRUE)

Arguments

method=='auto' attempts to use exact calculations when possible, otherwise ALA if available, otherwise Laplace approx.

Details

Let delta be the vector indicating inclusion/exclusion of each column of x in the model. The Gibbs algorithm sequentially samples from the posterior of each element in delta conditional on all the remaining elements in delta and the data. To do this it is necessary to evaluate the marginal likelihood for any given model. These have closed-form expression for the MOM prior, but for models with >15 variables these are expensive to compute and Laplace approximations are used instead (for the residual variance a log change of variables is used, which improves the approximation). For other priors closed forms are not available, so by default Laplace approximations are used. For the iMOM prior we also implement a Hybrid Laplace-IS which uses a Laplace approximation to evaluate the integral wrt beta and integrates wrt phi (residual variance) numerically.

It should be noted that Laplace approximations tend to under-estimate the marginal densities when the MLE for some parameter is very close to 0. That is, it tends to be conservative in the sense of excluding more variables from the model than an exact calculation would.

Finally, method=='plugin' provides a BIC-type approximation that is faster than exact or Laplace methods, at the expense of some accuracy. In non-sparse situations where models with many variables have large posterior probability method=='plugin' can be substantially faster.

For more details on the methods used to compute marginal densities see Johnson & Rossell (2012).

modelsearchBlockDiag uses the block search method described in Papaspiliopoulos & Rossell. Briefly, spectral clustering is run on X'X to cluster variables into blocks of blocksize and subsequently the Coolblock algorithm is used to define a sequence of models of increasing size. The exact integrated likelihood is evaluated for all models in this path, the best model chosen, and the scheme iteratively repeated to add and drop variables until convergence.

Value

Object of class msfit, which extends a list with elements

.

Author(s)

David Rossell

References

Johnson V.E., Rossell D. Non-Local Prior Densities for Default Bayesian Hypothesis Tests. Journal of the Royal Statistical Society B, 2010, 72, 143-170.

Johnson V.E., Rossell D. Bayesian model selection in high-dimensional settings. Journal of the American Statistical Association, 2012, 107, 649-660.

Papaspiliopoulos O., Rossell, D. Scalable Bayesian variable selection and model averaging under block orthogonal design. 2016

Rossell D., Rubio F.J. Tractable Bayesian variable selection: beyond normality. 2016

See Also

msfit-class for details on the output. postProb to obtain posterior model probabilities. coef.msfit for Bayesian model averaging estimates and intervals. predict.msfit for BMA estimates and intervals for user-supplied covariate values. plot.msfit for an MCMC diagnostic plot showing estimated marginal posterior inclusion probabilities vs. iteration number. rnlp to obtain posterior samples for the coefficients. marginalLikelihood to compute marginal densities for a given model.

Examples

#Simulate data
x <- matrix(rnorm(100*3),nrow=100,ncol=3)
theta <- matrix(c(1,1,0),ncol=1)
y <- x %*% theta + rnorm(100)
df <- data.frame(y, x)
#Specify prior parameters
priorCoef <- momprior()
priorModel <- modelunifprior()
#Alternative model space prior: 0.5 prior prob for including any covariate
priorModel <- modelbinomprior(p=0.5)
#Alternative: Beta-Binomial prior for model space
priorModel <- modelbbprior()
#Model selection
fit1 <- modelSelection(y ~ ., data=df, 
priorCoef=priorCoef, priorModel=priorModel)
coef(fit1) #BMA estimates, 95% intervals, marginal post prob
postProb(fit1) #posterior model probabilities

Bayesian variable selection for linear models via non-local priors.

Description

Bayesian model selection for linear, asymmetric linear, median and quantile regression under non-local or Zellner priors. p>>n can be handled.

modelSelection enumerates all models when feasible and uses a Gibbs scheme otherwise. See coef and coefByModel for estimates and posterior intervals of regression coefficients, and rnlp for posterior samples.

modelsearchBlockDiag seeks the highest posterior probability model using an iterative block search.

Usage

modelSelectionGGM(y, priorCoef=normalidprior(tau=1), 
priorModel=modelbinomprior(1/ncol(y)), 
priorDiag=exponentialprior(lambda=1), center=TRUE, scale=TRUE, 
global_proposal= 'regression', prob_global=0.5,
tempering= 0.5, truncratio= 100,
save_proposal= FALSE, niter=10^3, burnin= round(niter/10), 
updates_per_iter= ceiling(sqrt(ncol(y))), updates_per_column= 10,
sampler='Gibbs', pbirth=0.75, pdeath=0.5*(1-pbirth),
bounds_LIT, Omegaini='glasso-ebic', verbose=TRUE)

Arguments

Details

Let Omega be the inverse covariance matrix. A spike-and-slab prior is used. Specifically, independent priors are set on all Omega[j,k], and then a positive-definiteness truncation is added.

The prior on diagonal entries Omega[j,j] is given by priorDiag. Off-diagonal Omega[j,k] are equal to zero with probability given by modelPrior and, when non-zero, they are

Independent spike-and-slab priors are set on the off-diagonal entries of Omega, i.e. Omega[j,k]=0 with positive probability (spike) and otherwise arises from the distribution indicated in priorCoef (slab).

Inference is based on MCMC posterior sampling. All sampling algorithms proceed by updating Omega[,k] given y and Omega[,-k] (of course, Omega[k,] is also set to Omega[,k]). Omega[,k] is updated by first updating the set of non-zero entries (i.e. edges in the graphical model) using either Gibbs sampling or a proposal distribution (see below), and then the non-zero entries of Omega[,k] are updated from their exact posterior given the current set of edges.

An MCMC iteration consists of iterating over updates_per_iter columns chosen at random and, for each column, do updates_per_column proposals.

If global_proposal == "none", a local MCMC proposal is used to update what entries in Omega[,k] are zero. Local means that the proposed model is a neighbour of the current model, e.g. only one edge is added / killed. Available local kernels are Gibbs, birth-death-swap and LIT (Zhou et al 2022).

If global_proposal == "regression", a global proposal is used. A new model (set of non-zero entries in Omega[,k]) is proposed based on the posterior distribution of a linear regression of y[,k] on the other covariates. prob_global indicates the probability of using the global proposal, otherwise a local proposal is used. Proposal probabilities are tempered by raising them to the power tempering. Further, any model with probability below prob(top model) / truncratio is assigned proposal probability prob(top model) / truncratio.

Value

Posterior inference on the inverse covariance of y. Object of class msfit_ggm, which extends a list with elements

Author(s)

David Rossell

References

Zhou, Quan, et al. Dimension-free mixing for high-dimensional Bayesian variable selection. JRSS-B 2022, 84, 1751-1784

See Also

msfit_ggm-class for further details on the output. icov for the estimated precision (inverse covariance) matrix. coef.msfit_ggm for Bayesian model averaging estimates and intervals.

Examples

#Simulate data with p=3
Th= diag(3); Th[1,2]= Th[2,1]= 0.5
sigma= solve(Th)
z= matrix(rnorm(1000*3), ncol=3)
y= z 
#Obtain posterior samples
#y has few columns, initialize the MCMC at the sample precision matrix
#Else, leave the argument Omegaini in modelSelectionGGM empty
Omegaini= solve(cov(y))
fit= modelSelectionGGM(y, scale=FALSE, Omegaini=Omegaini)
#Parameter estimates, intervals, prob of non-zero
coef(fit)
#Estimated inverse covariance
icov(fit)
#Estimated inverse covariance, entries set to 0
icov(fit, threshold=0.95)
#Shows first posterior samples
head(fit$postSample)

Class "msPriorSpec"

Description

Stores the prior distributions to be used for Bayesian variable selection in normal regression models. This class can be used to specify the prior on non-zero regression coefficients, the model indicator or the nuisance parameters.

Usage

aic()
bic()
bicprior()
ic(penalty)
momprior(taustd=1, tau, tau.adj=10^6, r=1)
imomprior(tau, tau.adj=10^6)
emomprior(tau, tau.adj=10^6)
zellnerprior(taustd=1, tau, tau.adj=10^6)
normalidprior(taustd=1, tau, tau.adj=10^6)
icarplusprior(a=0.5, taustd=1, tau.adj=10^6)
exponentialprior(lambda = 1)
groupmomprior(taustd=1, tau, tau.adj=10^6)
groupimomprior(tau, tau.adj=10^6)
groupemomprior(tau, tau.adj=10^6)
groupzellnerprior(taustd=1, tau, tau.adj=10^6)
modelunifprior()
modelbinomprior(p=0.5)
modelbbprior(alpha.p=1, beta.p=1)
modelcomplexprior(c=1)
igprior(alpha=.01, lambda=.01)

Arguments

Details

DISCUSSION OF PRIOR ON PARAMETERS

Let beta=(beta_1,...,beta_p) be the regression coefficients for individual variables and delta=(delta_1,...,delta_q) those for grouped variables (e.g. factors or smooth terms in modelSelection).

momprior, emomprior, imomprior, zellnerprior, normalid and icarplus can be priors on both beta or delta. For further information see the vignette.

groupzellnerprior is the prior density on delta

p_z(\delta; \tau)= \prod_j N(\delta_j; 0, (\tau/p_j)) (X_j'X_j)^{-1}

where X_j are the design matrix columns associated to delta_j and p_j=ncol(X_j) is the number of covariates in the group (for groupmomprior, the term in the denominator is (p_j +2) instead of p_j). A default tau=n=nrow(X_j) mimics the unit information prior and implies that the ratio of variance explained by X_j / residual variance is expected to be 1 a priori. To set the dispersion in terms of unit information prior, taustd is also available.

groupmomprior adds a quadratic MOM penalty

p_m(delta; tau)= p_z(delta; tau * n) prod_j delta_j'X_j'X_jdelta_j ncol(X_j)/(tau * n * p_j / (p_j + 2))

and analogously for eMOM and iMOM. Note that unlike groupzellnerprior, the nrow(X_j) factor is already included in the code. This is done to give user introduced tau values a roughly similar meaning between momprior and groupmomprior.

DISCUSSION OF PRIOR ON MODELS

Under the uniform prior, the prior probability of any model is 1 / number of models.

Under the Binomial, Beta-Binomial and Complexity priors a model with k out of K active variables has prior probability P(Z=k) / (K choose k), where where Z ~ Binomial(K,p), Z ~ BetaBinomial(K,alpha.p,beta.p) or for the Complexity prior P(Z=k) proportional to 1/K^(c*k).

Objects from the Class

Objects can be created by calls of the form new("msPriorSpec", ...), but it is easier to use creator functions.

For priors on regression coefficients use momprior, imomprior or emomprior. For prior on model space modelunifprior, modelbinomprior modelbbprior, or modelcomplexprior. For prior on residual variance use igprior.

Slots

priorType:

Object of class "character". "coefficients" indicates that the prior is for the non-zero regression coefficients. "modelIndicator" that it is for the model indicator, and "nuisancePars" that it is for the nuisance parameteres. Several prior distributions are available for each choice of priorType, and these can be speicified in the slot priorDist.

priorDistr:

Object of class "character". If priorType=="coefficients", priorDistr can be equal to "pMOM", "piMOM", "peMOM", "zellner", "normalid", "groupMOM" or "groupzellner" (product moment, product inverse moment, product exponential moment, Zellner prior, normal prior with \Sigma=\mathbf{I}, respectively). If priorType=="modelIndicator", priorDistr can be equal to "uniform" or "binomial" to specify a uniform prior (all models equaly likely a priori) or a binomial prior, or to "complexity" for the Complexity prior of Castillo et al 2015. For a binomial prior, the prior inclusion probability for any single variable must be specified in slot priorPars['p']. For a beta-binomial prior, the Beta hyper-prior parameters must be in priorPars['alpha.p'] and priorPars['beta.p']. For the Complexity prior, the prior parameter must be in the slot priorPars['c']. If priorType=="nuisancePars", priorDistr must be equal to "invgamma". This corresponds to an inverse gamma distribution for the residual variance, with parameters specified in the slot priorPars.

priorPars:

Object of class "vector", where each element must be named. For priorDistr=='pMOM', there must be an element "r" (MOM power is 2r). For any priorDistr there must be either an element "tau" indicating the prior dispersion or elements "a.tau" and "b.tau" specifying an inverse gamma hyper-prior for "tau". Optionally, there may be an element "tau.adj" indicating the prior dispersion for the adjustment variables (i.e. not undergoing variable selection). If not defined, "tau.adj" is set to 0.001 by default. For priorDistr=='binomial', there must be either an element "p" specifying the prior inclusion probability for any single covariate, or a vector with elements "alpha.p" and "beta.p" specifying a Beta(alpha.p,beta.p) hyper-prior on p. For priorDistr=='invgamma' there must be elements "alpha" and "lambda". The prior for the residual variance is an inverse gamma with parameteres .5*alpha and .5*lambda.

Methods

No methods defined with class "msPriorSpec" in the signature.

Note

When new instances of the class are created a series of check are performed to ensure that a valid prior specification is produced.

Author(s)

David Rossell

References

Johnson VE, Rossell D. Non-Local Prior Densities for Default Bayesian Hypothesis Tests. Journal of the Royal Statistical Society B, 2010, 72, 143-170

Johnson VE, Rossell D. Bayesian model selection in high-dimensional settings. Journal of the American Statistical Association, 107, 498:649-660.

Rossell D, Abril O, Bhattacharya A. Approximate Laplace approximations for scalable model selection (2021). Journal of the Royal Statistical Society B, 83, 4, 853-879

See Also

See also modelSelection for an example of defining an instance of the class and perform Bayesian model selection.

Examples

showClass("msPriorSpec")

Class "msfit"

Description

Stores the output of Bayesian variable selection, as produced by function modelSelection. The class extends a list, so all usual methods for lists also work for msfit objects, e.g. accessing elements, retrieving names etc.

Methods are provided to compute posterior probabilities, obtaining regression coefficient estimates and posterior intervals (both via Bayesian model averaging and for individual models), and sampling from their posterior distribution, as indicated below.

Objects from the Class

Typically objects are automatically created by a call to modelSelection. Alternatively, objects can be created by calls of the form new("msfit",x) where x is a list with the adequate elements (see slots).

Slots

The class extends a list with elements:

postSample

matrix with posterior samples for the model indicator. postSample[i,j]==1 indicates that variable j was included in the model in the MCMC iteration i

postOther

postOther returns posterior samples for parameters other than the model indicator, i.e. basically hyper-parameters. If hyper-parameters were fixed in the model specification, postOther will be empty.

margpp

Marginal posterior probability for inclusion of each covariate. This is computed by averaging marginal post prob for inclusion in each Gibbs iteration, which is much more accurate than simply taking colMeans(postSample)

.

postMode

Model with highest posterior probability amongst all those visited

postModeProb

Unnormalized posterior prob of posterior mode (log scale)

postProb

Unnormalized posterior prob of each visited model (log scale)

family

Residual distribution, i.e. argument family when calling modelSelection

p

Number of variables

priors

Priors specified when calling modelSelection

ystd

For internal use. Stores the response variable, standardized if center or scale were set to TRUE

xstd

For internal use. Stores the covariates, standardized if center or scale were set to TRUE

stdconstants

For internal use. If center or scale were set to TRUE, stores the sample mean and standard deviation of the outcome and covariates

call

Stores info about the call, the formula used (if any), splines used etc

Methods

coef

Obtains posterior means and intervals via Bayesian model averaging

coefByModel

Obtains posterior means and intervals for individual models

plot

Shows estimated posterior inclusion probability for each parameter vs. number of MCMC iterations

predict

Obtains posterior means and intervals for given covariate values. These are posterior intervals for the mean, not posterior predictive intervals for the outcome

show

signature(object = "msfit"): Displays general information about the object.

postProb

signature(object = "msfit"): Extracts posterior model probabilities.

rnlp

signature(object = "msfit"): Obtain posterior samples for regression coefficients.

Author(s)

David Rossell

References

Johnson VE, Rossell D. Non-Local Prior Densities for Default Bayesian Hypothesis Tests. Journal of the Royal Statistical Society B, 2010, 72, 143-170

Johnson VE, Rossell D. Bayesian model selection in high-dimensional settings. Journal of the American Statistical Association, 107, 498:649-660.

See Also

See also modelSelection and rnlp .

Examples

showClass("msfit")

Class "msfit_ggm"

Description

Stores the output of Bayesian Gaussian graphical model selection and averaging, as produced by function modelSelectionGGM. The class extends a list, so all usual methods for lists also work for msfit_ggm objects, e.g. accessing elements, retrieving names etc.

Methods are provided to obtain parameter estimates, posterior intervals (Bayesian model averaging), and posterior probabilities of parameters being non-zero

Objects from the Class

Objects are created by a call to modelSelectionGGM.

Slots

The class extends a list with elements:

postSample

Sparse matrix (dgCMatrix) with posterior samples for the Gaussian precision (inverse covariance) parameters. Each row is a posterior sample. Within each row, only the upper-diagonal of the precision matrix is stored in a flat manner. The row and column indexes are stored in indexes

indexes

For each column in postSample, it indicates the row and column of the precision matrix

p

Number of variables

priors

Priors specified when calling modelSelection

Methods

coef

Obtain BMA posterior means, intervals and posterior probability of non-zeroes

plot

Shows estimated posterior inclusion probability for each parameter vs. number of MCMC iterations. Only up to the first 5000 parameters are shown

show

signature(object = "msfit_ggm"): Displays general information about the object.

Author(s)

David Rossell

See Also

modelSelectionGGM

Examples

showClass("msfit_ggm")

Plot estimated marginal prior inclusion probabilities

Description

Plot marginal prior inclusion probabilities as estimated by cil versus regression coefficients for the treatment(s) equation(s)

Usage

 plotprior(object, xlab, ylab, ylim=c(0,1), ...)

Arguments

Value

A plot of prior inclusion probabilities vs treatment regression coefficients (dots). The line shows the (empirical Bayes) fit

Author(s)

David Rossell

See Also

cil

Examples

#See help(cil)

Obtain posterior model probabilities

Description

Obtain posterior model probabilities after running Bayesian model selection

Usage

 postProb(object, nmax, method='norm')

Arguments

Value

A data.frame with posterior model probabilities in column pp. Column modelid indicates the indexes of the selected covariates (empty for the null model with no covariates)

For msfit_ggm objects, a list with posterior probabilities (element "pp") and the corresponding model (element "modelid") indicating what edges are non-zero

For mixturebf objects, posterior probabilities for the specified number of components

For localtest objects, posterior probabilities of a local covariate effect at various regions

Author(s)

David Rossell

See Also

modelSelection to perform model selection

Examples

#See help(modelSelection)

Extract posterior samples from an object

Description

Obtain posterior model probabilities after running Bayesian model selection

Usage

 postSamples(object)

Arguments

Value

For objects of class mixturebf, a list with one element for each considered number of mixture components.

Each element in the list contains posterior samples on the mixture weights (eta) and other component-specific parameters such as means (mu) and Cholesky decomposition of the inverse covariance matrix (cholSigmainv)

Author(s)

David Rossell

Examples

#See help(bfnormmix)

Moment and inverse moment prior elicitation

Description

priorp2g finds the g value giving priorp prior probability to the interval (-q,q).

Usage

priorp2g(priorp, q, nu=1, prior=c("iMom", "normalMom", "tMom"))

Arguments

Details

See pmom and pimom for the MOM/iMOM cumulative distribution functions.

Value

priorp2g returns g giving priorp prior probability to the interval (-q,q).

Author(s)

David Rossell rosselldavid@gmail.com

References

See http://rosselldavid.googlepages.com for technical reports.

See Also

pmom , pimom

Examples

#find g value giving 0.05 probability to interval (-.2,.2)
priorp <- .05; q <- .2
gmom <- priorp2g(priorp=priorp, q=q, prior='normalMom')
gimom <- priorp2g(priorp=priorp, q=q, prior='iMom')
gmom
gimom

Posterior sampling for regression parameters

Description

Gibbs sampler for linear, generalized linear and survival models under product non-local priors, Zellner's prior and a Normal approximation to the posterior. Both sampling conditional on a model and Bayesian model averaging are implemented (see Details).

If x and y not specified samples from non-local priors/posteriors with density proportional to d(theta) N(theta; m, V) are produced, where d(theta) is the non-local penalty term.

Usage

rnlp(y, x, m, V, msfit, outcometype, family, priorCoef, priorGroup,
priorVar, priorprec, isgroup,
niter=10^3, burnin=round(niter/10), thinning=1, pp='norm')

Arguments

Details

The algorithm is implemented for product MOM (pMOM), product iMOM (piMOM) and product eMOM (peMOM) priors. The algorithm combines an orthogonalization that provides low serial correlation with a latent truncation representation that allows fast sampling.

When y and x are specified sampling is for the linear regression posterior. When argument msfit is left missing, posterior sampling is for the full model regressing y on all covariates in x. When msfit is specified each model is drawn with probability given by postProb(msfit). In this case, a Bayesian Model Averaging estimate of the regression coefficients can be obtained by applying colMeans to the rnlp ouput matrix.

When y and x are left missing, sampling is from a density proportional to d(theta) N(theta; m,V), where d(theta) is the non-local penalty (e.g. d(theta)=prod(theta^(2r)) for the pMOM prior).

Value

Matrix with posterior samples

Author(s)

David Rossell

References

D. Rossell and D. Telesca. Non-local priors for high-dimensional estimation, 2014. http://arxiv.org/pdf/1402.5107v2.pdf

See Also

modelSelection to perform model selection and compute posterior model probabilities. For more details on prior specification see msPriorSpec-class .

Examples

#Simulate data
x <- matrix(rnorm(100*3),nrow=100,ncol=3)
theta <- matrix(c(1,1,0),ncol=1)
y <- x %*% theta + rnorm(100)
fit1 <- modelSelection(y=y, x=x, center=FALSE, scale=FALSE)
th <- rnlp(msfit=fit1, niter=100)
colMeans(th)