| Title: | Estimating Dynamic Correlation |
|---|---|
| Description: | Implementations for two different Bayesian models of differential co-expression. scdeco.cop() fits the bivariate Gaussian copula model from Zichen Ma, Shannon W. Davis, Yen-Yi Ho (2023) <doi:10.1111/biom.13701>, while scdeco.pg() fits the bivariate Poisson-Gamma model from Zhen Yang, Yen-Yi Ho (2022) <doi:10.1111/biom.13457>. |
| Authors: | Anderson Bussing [aut, cre], Yen-Yi Ho [aut, ths], Zichen Ma [aut], Zhen Yang [aut] |
| Maintainer: | Anderson Bussing <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.1.0 |
| Built: | 2024-10-03 06:39:25 UTC |
| Source: | CRAN |
Copula dynamic correlation fitting function
scdeco.cop(y, x, marginals, w = NULL, n.mcmc = 5000, burn = 1000, thin = 10)scdeco.cop(y, x, marginals, w = NULL, n.mcmc = 5000, burn = 1000, thin = 10)
y |
2-column matrix of observations |
x |
covariates |
marginals |
length-2 vector with strings of the two marginals |
w |
(optional) |
n.mcmc |
number of mcmc iterations to run |
burn |
how many of the mcmc iterations to burn |
thin |
how much to thin the mcmc iterations |
matrix with mcmc samples as rows and columns corresponding to the different parameters
n <- 1000 x.use = rnorm(n) w.use = runif(n,-1,1) eta1.use = c(-2.2, 0.7) eta2.use = c(-2, 0.8) beta1.use = c(1,0.5) beta2.use = c(1,1) alpha1.use = 7 alpha2.use = 3 tau.use = c(-0.2, .3) marginals.use <- c("ZINB", "ZIGA") y.use <- scdeco.sim.cop(marginals=marginals.use, x=x.use, eta1.true=eta1.use, eta2.true=eta2.use, beta1.true=beta1.use, beta2.true=beta2.use, alpha1.true=alpha1.use, alpha2.true=alpha2.use, tau.true=tau.use, w=w.use) mcmc.out <- scdeco.cop(y=y.use, x=x.use, marginals=marginals.use, w=w.use, n.mcmc=10, burn=0, thin=1) # n.mcmc=1000, burn=100, thin=5) lowerupper <- t(apply(mcmc.out, 2, quantile, c(0.025, 0.5, 0.975))) estmat <- cbind(lowerupper[,1], c(eta1.use, eta2.use, beta1.use, beta2.use, alpha1.use, alpha2.use, tau.use), lowerupper[,c(2,3)]) colnames(estmat) <- c("lower", "trueval", "estval", "upper") estmatn <- 1000 x.use = rnorm(n) w.use = runif(n,-1,1) eta1.use = c(-2.2, 0.7) eta2.use = c(-2, 0.8) beta1.use = c(1,0.5) beta2.use = c(1,1) alpha1.use = 7 alpha2.use = 3 tau.use = c(-0.2, .3) marginals.use <- c("ZINB", "ZIGA") y.use <- scdeco.sim.cop(marginals=marginals.use, x=x.use, eta1.true=eta1.use, eta2.true=eta2.use, beta1.true=beta1.use, beta2.true=beta2.use, alpha1.true=alpha1.use, alpha2.true=alpha2.use, tau.true=tau.use, w=w.use) mcmc.out <- scdeco.cop(y=y.use, x=x.use, marginals=marginals.use, w=w.use, n.mcmc=10, burn=0, thin=1) # n.mcmc=1000, burn=100, thin=5) lowerupper <- t(apply(mcmc.out, 2, quantile, c(0.025, 0.5, 0.975))) estmat <- cbind(lowerupper[,1], c(eta1.use, eta2.use, beta1.use, beta2.use, alpha1.use, alpha2.use, tau.use), lowerupper[,c(2,3)]) colnames(estmat) <- c("lower", "trueval", "estval", "upper") estmat
ZENCO Poisson Gamma dynamic correlation fitting function
scdeco.pg( dat, b0, b1, adapt_iter = 100, update_iter = 100, coda_iter = 1000, coda_thin = 5, coda_burnin = 100 )scdeco.pg( dat, b0, b1, adapt_iter = 100, update_iter = 100, coda_iter = 1000, coda_thin = 5, coda_burnin = 100 )
dat |
matrix containing expression values as first two columns and covariate as third column |
b0 |
intercept of zinf parameter |
b1 |
slope of zinf parameter |
adapt_iter |
number of adaptation iterations in the jags.model function |
update_iter |
update iterations in the update function |
coda_iter |
number of iterations for the coda.sample function |
coda_thin |
how much to thin the resulting MCMC output |
coda_burnin |
how many iterations to burn before beginning coda sample collection |
MCMC samples that have been adapted, burned, and thinned
phi1_use <- 4 phi2_use <- 4 phi3_use <- 1/7 mu1_use <- 15 mu2_use <- 15 mu3_use <- 7 b0_use <- -3 b1_use <- 0.1 tau0_use <- -2 tau1_use <- 0.4 simdat <- scdeco.sim.pg(N=1000, b0=b0_use, b1=b1_use, phi1=phi1_use, phi2=phi2_use, phi3=phi3_use, mu1=mu1_use, mu2=mu2_use, mu3=mu3_use, tau0=tau0_use, tau1=tau1_use) zenco_out <- scdeco.pg(dat=simdat, b0=b0_use, b1=b1_use, adapt_iter=1, # 500, update_iter=1, # 500, coda_iter=5, # 5000, coda_thin=1, # 10, coda_burnin=0) # 1000 boundsmat <- cbind(zenco_out$quantiles[,1], c(1/phi1_use, 1/phi2_use, 1/phi3_use, mu1_use, mu2_use, mu3_use, tau0_use, tau1_use), zenco_out$quantiles[,c(3,5)]) colnames(boundsmat) <- c("lower", "true", "est", "upper") boundsmatphi1_use <- 4 phi2_use <- 4 phi3_use <- 1/7 mu1_use <- 15 mu2_use <- 15 mu3_use <- 7 b0_use <- -3 b1_use <- 0.1 tau0_use <- -2 tau1_use <- 0.4 simdat <- scdeco.sim.pg(N=1000, b0=b0_use, b1=b1_use, phi1=phi1_use, phi2=phi2_use, phi3=phi3_use, mu1=mu1_use, mu2=mu2_use, mu3=mu3_use, tau0=tau0_use, tau1=tau1_use) zenco_out <- scdeco.pg(dat=simdat, b0=b0_use, b1=b1_use, adapt_iter=1, # 500, update_iter=1, # 500, coda_iter=5, # 5000, coda_thin=1, # 10, coda_burnin=0) # 1000 boundsmat <- cbind(zenco_out$quantiles[,1], c(1/phi1_use, 1/phi2_use, 1/phi3_use, mu1_use, mu2_use, mu3_use, tau0_use, tau1_use), zenco_out$quantiles[,c(3,5)]) colnames(boundsmat) <- c("lower", "true", "est", "upper") boundsmat
Simulating from copula model
scdeco.sim.cop( marginals, x, eta1.true, eta2.true, beta1.true, beta2.true, alpha1.true, alpha2.true, tau.true, w = NULL )scdeco.sim.cop( marginals, x, eta1.true, eta2.true, beta1.true, beta2.true, alpha1.true, alpha2.true, tau.true, w = NULL )
marginals |
provide vector of length 2 of which marginals to use |
x |
covariate matrix |
eta1.true |
zero-inflation parameters for marginal 1 |
eta2.true |
zero-inflation parameters for marginal 2 |
beta1.true |
mean coefficients for marginal 1 |
beta2.true |
mean coefficients for marginal 2 |
alpha1.true |
second parameter coefficients for marginal 1 |
alpha2.true |
second parameter coefficients for marginal 2 |
tau.true |
coefficients for correlation |
w |
(optional) covariate matrix for zero-inflation portion |
matrix with values simulated from copula model
n <- 2500 x.use = rnorm(n) w.use = runif(n,-1,1) eta1.use = c(-2.2, 0.7) eta2.use = c(-2, 0.8) beta1.use = c(1,0.5) beta2.use = c(1,1) alpha1.use = 7 alpha2.use = 3 tau.use = c(-0.2, .3) marginals.use <- c("ZINB", "ZIGA") y.use <- scdeco.sim.cop(marginals=marginals.use, x=x.use, eta1.true=eta1.use, eta2.true=eta2.use, beta1.true=beta1.use, beta2.true=beta2.use, alpha1.true=alpha1.use, alpha2.true=alpha2.use, tau.true=tau.use, w=w.use) y.use[1:10,]n <- 2500 x.use = rnorm(n) w.use = runif(n,-1,1) eta1.use = c(-2.2, 0.7) eta2.use = c(-2, 0.8) beta1.use = c(1,0.5) beta2.use = c(1,1) alpha1.use = 7 alpha2.use = 3 tau.use = c(-0.2, .3) marginals.use <- c("ZINB", "ZIGA") y.use <- scdeco.sim.cop(marginals=marginals.use, x=x.use, eta1.true=eta1.use, eta2.true=eta2.use, beta1.true=beta1.use, beta2.true=beta2.use, alpha1.true=alpha1.use, alpha2.true=alpha2.use, tau.true=tau.use, w=w.use) y.use[1:10,]
Simulating from ZENCO Model
scdeco.sim.pg( N, b0, b1, phi1, phi2, mu1, mu2, tau0, tau1, mu3, phi3, tau2 = NULL, tau3 = NULL, xc = NULL )scdeco.sim.pg( N, b0, b1, phi1, phi2, mu1, mu2, tau0, tau1, mu3, phi3, tau2 = NULL, tau3 = NULL, xc = NULL )
N |
size of sample to be generated |
b0 |
intercept of zinf parameter |
b1 |
slope of zinf parameter |
phi1 |
over-dispersion parameter of first marginal |
phi2 |
over-dispersion parameter of second marginal |
mu1 |
mean parameter of first marginal |
mu2 |
mean parameter of second marginal |
tau0 |
intercept of correlation |
tau1 |
slope of of correlation |
mu3 |
mean parameter of covariate vector |
phi3 |
over-dispersion parameter of covariate vector |
tau2 |
(optional) correlation coefficient on optional xc covariate vector |
tau3 |
(optional) correlation coefficient on interaction between x3 and xc |
xc |
(optional) secondary covariate to be regressed |
a matrix with expressions as first two columns and covariates as remaining columns
phi1_use <- 4 phi2_use <- 4 phi3_use <- 1/6 mu1_use <- 15 mu2_use <- 15 mu3_use <- 7 b0_use <- 0.6882 b1_use <- -0.2995 tau0_use <- 0.07 tau1_use <- 0.05 simdat <- scdeco.sim.pg(N=1000, b0=b0_use, b1=b1_use, phi1=phi1_use, phi2=phi2_use, phi3=phi3_use, mu1=mu1_use, mu2=mu2_use, mu3=mu3_use, tau0=tau0_use, tau1=tau1_use) simdat[1:10,]phi1_use <- 4 phi2_use <- 4 phi3_use <- 1/6 mu1_use <- 15 mu2_use <- 15 mu3_use <- 7 b0_use <- 0.6882 b1_use <- -0.2995 tau0_use <- 0.07 tau1_use <- 0.05 simdat <- scdeco.sim.pg(N=1000, b0=b0_use, b1=b1_use, phi1=phi1_use, phi2=phi2_use, phi3=phi3_use, mu1=mu1_use, mu2=mu2_use, mu3=mu3_use, tau0=tau0_use, tau1=tau1_use) simdat[1:10,]