# NOT RUN {
# We consider a simple simulated data set to illustrate the use of the DQ function.
# We have simulated the following situation.
# We have two groups, for example, 5 patients with tumor 1 and 5 patients
# with tumor 2. For each patient 100 variables are measured, for example,
# gene expression levels. Besides, the distributions of 30 of the variables
# are different in the two groups, and the differences are in location.
# We consider a collection of densities {f1=N(0,1), f2=N(0,1/4), f3=N(2,1), f4=N(2,1/4)}. In
# the first group (tumor 1) the sample of each variable (gene) comes from one of
# the four densities with the same probability 1/4. On the other hand, in the second
# group the sample of each variable comes from the same density as in the first
# group except for 30 randomly selected variables for which the density changes
# producing location differences. Specifically, if the variable follows
# f1 in the first group, its density, in the second group, is f3 producing a
# change on its location parameter.The situation for the other cases is as follows:
# the density f2 (f3 or f4) in group 1 leads to density f4 (f1 or f2, respectively)
# in the second one.
set.seed(123)
p <- 100
n = m = 5
inds <- sample(1:4, p, replace = TRUE)
X <- matrix(rep(0, n * p), ncol = n)
for (j in 1:p){
if (inds[j] == 1){
X[j, ] <- rnorm(n)}
if (inds[j] == 2){
X[j, ] <- rnorm(n, sd = sqrt(1/4))
}
if (inds[j] == 3){
X[j, ]<-rnorm(n, mean = 2)
}
if (inds[j] == 4){
X[j, ]<-rnorm(n, mean = 2, sd = sqrt(1/4))
}
}
rho <- 0.3
ind <- sample(1:p, rho * p)
li <- length(ind)
indsy <- inds
for (l in 1:li){
if (indsy[ind[l]] == 1){indsy[ind[l]] = 3} else{
if (indsy[ind[l]] == 2){indsy[ind[l]] = 4} else {
if (indsy[ind[l]] == 3){indsy[ind[l]] = 1}
else{indsy[ind[l]] = 2}}}}
Y <- matrix(rep(0, m * p), ncol = m)
for (j in 1:p){
if (indsy[j] == 1){
Y[j, ] <- rnorm(m)}
if (indsy[j] == 2){
Y[j, ] <- rnorm(m, sd = sqrt(1/4))
}
if (indsy[j] == 3){
Y[j, ] <- rnorm(m, mean = 2)
}
if (indsy[j] == 4){
Y[j, ]<- rnorm(m, mean = 2, sd = sqrt(1/4))
}
}
# We can see which are the variables with different distributions in the two data sets.
dif <- which(inds != indsy)
# Cross table for (X,Y) density indexes:
table(inds,indsy)
# Our interest is to identify which variables have a different distribution in the two groups.
# Hence, since the differences between the distributions are in location, we applied
# Wilcoxon-Mann-Whitney test to verify for each variable the equality of distribution
# in the two groups.
# }
# NOT RUN {
library(exactRankTests)
library(coin)
# We compute the P-values
p <- nrow(X)
pv <- 1:p
for (i in 1:p){
pv[i] <-wilcox.exact(X[i, ], Y[i, ])$p.value
}
# When the sample size is small, in this case n=m=5, the distribution of
# the Wilcoxon's statistic is calibrated using an exact permutation test. Hence,
# the P-values are homogeneous discrete uniform distributed with support points
# of such distribution:
ss <- c(1, 2, 4, 7, 12, 19, 28, 39, 53,69, 87, 106, 126) / 126
# When the number of P-values is large enough "ss" is equal to:
sort(unique(pv))
# For details about the Wilcoxon-Mann-Whitney test and its exact distribution
# see Section 9.2 of Gibbons and Chakraborti (1992).
# We apply Chen method:
R <- DQ(pv, ss = ss, method = "Chen")
# The estimate of the proportion of null P-values:
R$pi0
# Summary of the vector of the estimated q-values:
summary(R$q.values)
# How many null hypotheses are rejected?
alpha <- 0.05
sum(R$q.values < alpha)
# Which variables correspond to such null hypotheses?
which(R$q.values < alpha)
# Classification table (Decision at nominal level alpha vs. reality):
table(R$q.values < alpha,inds != indsy)
# The conclusion from the previous table is that Chen method reports
# 21 true positives and 9 false negatives.
# We can also apply Liang and SS methods as follows.
RLiang <- DQ(pv, ss = ss, method = "Liang")
RSS <- DQ(pv, ss = ss, method = "SS")
# The next graphic help us to see that Liang method (for discrete P-values)
# is more powerful than SS method (only suitable for continuous P-values).
plot(RSS$q.values,RLiang$q.values)
abline(a = 0, b = 1, col = 2, lty = 2)
# We consider a poisson setting to show a case where the support of the discrete
# distribution of the P-values is infinte.
# We generate 100 values of a poisson distribution with event rate 10.
# Then we compute the probability that each of the values come from a
# a poisson distribution with event rate 10. This set of probabilities
# are considered as our set of P-values.
p <- 100
N <- rpois(p, lambda = 10)
pv <- 1:p
for(i in 1:p){
pv[i] <- ppois(N[i], lambda = 10)
}
# It is well know that the support of a poisson distribution is infinite
# and equal to the natural numbers. Hence to know the support of the P-values
# defined above, we compute for 1,..., 50 their corresponding P-value.
# We only considered 50 values because for large values than 50 the P-value is 1 again.
nn_0 <- 50
ss <- 1:(nn_0 + 1)
for (i in 1:(nn_0 + 1)){
ss[i] <- ppois(i - 1, lambda = 10)
}
# We eliminate repeated values.
ss <- unique(ss)
# For Chen method the relevant support points are only the values below tau[100] = 0.5.
# We define the support ss as such values. Then, we can apply Chen method. Of
# course s_inf = TRUE.
indicator <- which(ss <= 0.5)
ss <- ss[indicator]
R <- DQ(pv, ss = ss, s_inf = "TRUE", method = "Chen")
# For Liang method the relevant support values are also the values below 0.5, hence
# ss defined above is suitable.
R <- DQ(pv, ss = ss, s_inf = "TRUE", method = "Liang")
# For Rand method the relevant support values are those ones with match with
# the P-values, and also the largest support points smaller than each one of such P-values.
# Hence, ss only includes the relevant points, as we can see below.
pv_unique <- unique(pv)
p_u <- length(pv_unique)
ind <- 1:p_u
for(i in 1:p_u){
ind[i] <- which(ss == pv_unique[i])
}
p_u == length(ind)
ind_minus <- ind - 1
ind_final <- unique(sort(c(ind, ind_minus)))
ss <- ss[ind_final]
# Now, we can apply Rand method.
R <- DQ(pv, ss = ss, s_inf = "TRUE", method = "Rand")
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab