Description Usage Arguments Details Value Examples
MAP test
1 |
est_result |
estimation_fun result: |
dd |
True index to check the true FDR |
Type |
Type of hypothesis
|
nn |
Number of QMC nodes used to estimate the test statistics |
The vector of read counts for gene g, treatment group i, replicate j, at time point t,Y_{gij}(t), follows a Negative Binomial distribution parameterized mean λ_{gi} and φ_g, where E[Y_{gij}(t)] = λ_{gi}(t). λ_{gi}(t) is further modeled as λ_{gi}(t) = S_{ij} \exp[η_{g1}I_{i = 2} + B'(t)η_{g2}I_{i = 2} + B'(t)τ_{g}]. We have B'(t) are design matrix, which is constructed by 2 orthorgonal polynomial bases.
t = 1,..., n_treat (or n_control if control group);
j = 1,..., n_rep;
g = 1,...,G; and
[η_{g1}, η_{g2}, τ_{g}] ~ 4-component gausssian mixture model. After the parameter estimation we run the MAP test to detect the DE genes.
A list of results for DE analysis.
Type Type of hypothesis of interest
FDR_final FDR estimation of each hypothesis
FDR_real The true FDR if dd is known
power Power if dd is known
T_x Likelihood result for each component
test_stat Test statistics for each hypothesis
ct_final_all Test statistics cut-off value at each estimated FDR level
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 | library(matlib)
n_basis = 2
n_control = 10
n_treat = 10
n_rep = 3
tt_treat = c(1:n_treat)/n_treat
nt = length(tt_treat)
ind_t = sort(sample(c(1:nt), n_control))
tt = tt_treat[ind_t]
tttt = c(rep(tt, n_rep), rep(tt_treat, n_rep))
z = x = matrix(0, length(tttt), n_basis)
z[,1] = 1.224745*tttt
z[,2] = -0.7905694 + 2.371708*tttt^2
x[,1] = z[,1] - Proj(z[,1], rep(1, length(tttt)))
x[,2] = z[,2] - Proj(z[,2], rep(1, length(tttt))) - Proj(z[,2], x[,1])
Y1 = data_generation(G = 100,
n_control = n_control,
n_treat = n_treat,
n_rep = n_rep,
k_real = 4,
sigma2_r = rep(1, 2),
sigma1_2_r = 1,
sigma2_2_r = c(3,2),
mu1_r = 4,
phi_g_r = rep(1,100),
p_k_real = c(0.7, 0.1, 0.1, 0.1),
x = x)
aaa <- proc.time()
est_result <- estimation_fun(n_control = n_control,
n_treat = n_treat,
n_rep = n_rep,
x = x,
Y1 = Y1,
nn = 300,
k = 4,
phi = NULL,
type = 2,
tttt = tttt)
aaa1 <- proc.time()
aaa1 - aaa
G <- 100
k_real <- 4
p_k_real <- c(0.7, 0.1, 0.1, 0.1)
dd = rep(c(0:(k_real-1)), p_k_real * G)
result = MAP_test(est_result = est_result, Type = c(1:6), dd = dd, nn = 300)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.