mzipBvs | R Documentation |
The function can be used to perform variable selection for multivariate zero-inflated count responses.
mzipBvs(Y, lin.pred, data, model = "generalized", offset = NULL, hyperParams, startValues, mcmcParams)
Y |
a data.frame containing q count outcomes from |
lin.pred |
a list containing three formula objects: the first formula specifies the p_z covariates for which variable selection is to be performed in the binary component of the model; the second formula specifies the p_x covariates for which variable selection is to be performed in the count part of the model; the third formula specifies the p_0 confounders to be adjusted for (but on which variable selection is not to be performed) in the regression analysis. |
data |
a data.frame containing the variables named in the formulas in |
model |
a character that specifies the type of model: A generalized multivariate Bayesian variable selection method of Lee et al.(2018) can be implemented by setting |
offset |
an optional numeric vector with an a priori known component to be included as the linear predictor in the count part of model. |
hyperParams |
a list containing lists or vectors for hyperparameter values in hierarchical models. Components include,
|
startValues |
a numeric vector containing starting values for model parameters. See Examples below. |
mcmcParams |
a list containing variables required for MCMC sampling. Components include,
|
mzipBvs
returns an object of class mzipBvs
.
Kyu Ha Lee, Brent A. Coull, Jacqueline R. Starr
Maintainer: Kyu Ha Lee <klee15239@gmail.com>
Lee, K. H., Coull, B. A., Moscicki, A.-B., Paster, B. J., Starr, J. R. (2020),
Bayesian variable selection for multivariate zero-inflated models: application to microbiome count data, Biostatistics, Volume 21, Issue 3, Pages 499-517.
## Not run: # loading a data set data(simData_mzip) Y <- simData_mzip$Y data <- simData_mzip$X n = dim(Y)[1] q = dim(Y)[2] form.bin <- as.formula(~cov.1) form.count <- as.formula(~cov.1) form.adj <- as.formula(~1) lin.pred <- list(form.bin, form.count, form.adj) Xmat0 <- model.frame(lin.pred[[1]], data=data) Xmat1 <- model.frame(lin.pred[[2]], data=data) Xmat_adj <- model.frame(lin.pred[[3]], data=data) p_adj = ncol(Xmat_adj) p0 <- ncol(Xmat0) + p_adj p1 <- ncol(Xmat1) + p_adj nonz <- rep(NA, q) for(j in 1:q) nonz[j] <- sum(Y[,j] != 0) ##################### ## Hyperparameters ## ## Generalized model ## rho0 <- q + 3 + 1 Psi0 <- diag(3, q) mu_alpha0 <- 0 mu_alpha <- rep(0, q) mu_beta0 <- 0 mu_beta <- rep(0, q) a_alpha0 <- 0.7 b_alpha0 <- 0.7 a_alpha <- rep(0.7, p0) b_alpha <- rep(0.7, p0) a_beta0 <- 0.7 b_beta0 <- 0.7 a_beta <- rep(0.7, p1) b_beta <- rep(0.7, p1) v_beta = rep(1, q) omega_beta = rep(0.1, p1-p_adj) v_alpha = rep(1, q) omega_alpha = rep(0.1, p0-p_adj) ## hyperParams.gen <- list(rho0=rho0, Psi0=Psi0, mu_alpha0=mu_alpha0, mu_alpha=mu_alpha, mu_beta0=mu_beta0, mu_beta=mu_beta, a_alpha0=a_alpha0, b_alpha0=b_alpha0, a_alpha=a_alpha, b_alpha=b_alpha, a_beta0=a_beta0, b_beta0=b_beta0, a_beta=a_beta, b_beta=b_beta, v_beta=v_beta, omega_beta=omega_beta, v_alpha=v_alpha, omega_alpha=omega_alpha) ################### ## MCMC SETTINGS ## ## Setting for the overall run ## numReps <- 100 thin <- 1 burninPerc <- 0.5 ## Settings for storage ## storeV <- TRUE storeW <- TRUE ## Tuning parameters for specific updates ## ## - Generalized model beta0.prop.var <- 0.5 alpha.prop.var <- 0.5 beta.prop.var <- 0.5 V.prop.var <- 0.05 ## mcmc.gen <- list(run=list(numReps=numReps, thin=thin, burninPerc=burninPerc), storage=list(storeV=storeV, storeW=storeW), tuning=list(beta0.prop.var=beta0.prop.var, alpha.prop.var=alpha.prop.var, beta.prop.var=beta.prop.var, V.prop.var=V.prop.var)) ##################### ## Starting Values ## ## Generalized model ## B <- matrix(0.1, p1, q, byrow = T) A <- matrix(0.1, p0, q, byrow = T) V <- matrix(rnorm(n*q, 0, 0.1), n, q) W <- matrix(rnorm(n*q, 0, 0.1), n, q) beta0 <- log(as.vector(apply(Y, 2, mean))) alpha0 <- log(nonz/n / ((n-nonz)/n)) Sigma_V <- matrix(0, q, q) diag(Sigma_V) <- 1 R <- matrix(0, q, q) diag(R) <- 1 sigSq_alpha0 <- 1 sigSq_alpha <- rep(1, p0) sigSq_beta0 <- 1 sigSq_beta <- rep(1, p1) startValues.gen <- list(B=B, A=A, V=V, W=W, beta0=beta0, alpha0=alpha0, R=R, sigSq_alpha0=sigSq_alpha0, sigSq_alpha=sigSq_alpha, sigSq_beta0=sigSq_beta0, sigSq_beta=sigSq_beta, Sigma_V=Sigma_V) ################################### ## Fitting the generalized model ## ################################### fit.gen <- mzipBvs(Y, lin.pred, data, model="generalized", offset=NULL, hyperParams.gen, startValues.gen, mcmc.gen) print(fit.gen) summ.fit.gen <- summary(fit.gen); names(summ.fit.gen) summ.fit.gen ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.