mmcm.resamp: The permuted modified maximum contrast method

Description Usage Arguments Details Value References See Also Examples

View source: R/mmcm.resamp.R

Description

This function gives P-value for the permuted modified maximum contrast method.

Usage

1
2
mmcm.resamp(x, g, contrast, alternative = c("two.sided", "less", "greater"),
  nsample = 20000, abseps = 0.001, seed = NULL)

Arguments

x

a numeric vector of data values

g

a integer vector giving the group for the corresponding elements of x

contrast

a numeric contrast coefficient matrix for permuted modified maximum contrast statistics

alternative

a character string specifying the alternative hypothesis, must be one of "two.sided" (default), "greater" or "less". You can specify just the initial letter.

nsample

specifies the number of resamples (default: 20000)

abseps

specifies the absolute error tolerance (default: 0.001)

seed

a single value, interpreted as an integer; see set.seed() function. (default: NULL)

Details

mmcm.resamp performs the permuted modified maximum contrast method that is detecting a true response pattern under the unequal sample size situation.

Y_ij (i = 1, 2, ...; j = 1, 2, ..., n_i) is an observed response for j-th individual in i-th group.

C is coefficient matrix for permuted modified maximum contrast statistics (i x k matrix, i: No. of groups, k: No. of pattern).

C = (c_1, c_2, ..., c_k)^T

c_k is coefficient vector of k-th pattern.

c_k = (c_k1, c_k2, ..., c_ki)^T (sum from i of c_ki = 0)

M_max is a permuted modified maximum contrast statistic.

Ybar_i = (sum from j of Y_ij) / n_i, Ybar = (Ybar_1, Ybar_2, ..., Ybar_i, ..., Ybar_a)^T (a x 1 vector), M_k = c_k^t Ybar / (c_k^t c_k)^(1/2),

M_max = max(M_1, M_2, ..., M_k).

Consider testing the overall null hypothesis H_0: μ_1=μ_2=…=μ_i, versus alternative hypotheses H_1 for response petterns (H_1: μ_1<μ_2<…<μ_i, μ_1=μ_2<…<μ_i, μ_1<μ_2<…=μ_i). The P-value for the probability distribution of M_max under the overall null hypothesis is

P-value = Pr(M_max > m_max | H0)

m_max is observed value of statistics. This function gives distribution of M_max by using the permutation method, follow algorithm:

1. Initialize counting variable: COUNT = 0. Input parameters: NRESAMPMIN (minimum resampling count, we set 1000), NRESAMPMAX (maximum resampling count), and epsilon (absolute error tolerance).

2. Calculate m_max that is the observed value of the test statistic.

3. Let y_ij(r) donate data, which are sampled without replacement, and independently, form observed value y_ij. Where, (r) is suffix of the resampling number (r = 1, 2, ...).

4. Calculate m(r)_max from y_ij(r). If m(r)_max > m_max, then increment the counting variable: COUNT = COUNT + 1. Calculate approximate P-value hat-p(r) = COUNT / r, and the simulation standard error [hat-p(r) (1 - hat-p(r)) / r]^(1/2).

5. Repeat 3–4, while r > 1000 and 3.5 SE(hat-p(r)) < epsilon (corresponding to 99% confidence level), or NRESAMPMAX times. Output the approximate P-value hat-p(r).

Value

statistic

the value of the test statistic with a name describing it.

p.value

the p-value for the test.

alternative

a character string describing the alternative hypothesis.

method

the type of test applied.

contrast

a character string giving the names of the data.

contrast.index

a suffix of coefficient vector of the kth pattern that gives permuted modified maximum contrast statistics (row number of the coefficient matrix).

error

estimated absolute error and,

msg

status messages.

References

Nagashima, K., Sato, Y., Hamada, C. (2011). A modified maximum contrast method for unequal sample sizes in pharmacogenomic studies Stat Appl Genet Mol Biol. 10(1): Article 41. http://dx.doi.org/10.2202/1544-6115.1560

Sato, Y., Laird, N.M., Nagashima, K., et al. (2009). A new statistical screening approach for finding pharmacokinetics-related genes in genome-wide studies. Pharmacogenomics J. 9(2): 137–146. http://www.ncbi.nlm.nih.gov/pubmed/19104505

See Also

mmcm.mvt

Examples

 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
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
## Example 1 ##
#  true response pattern: dominant model c=(1, 1, -2)
set.seed(136885)
x <- c(
  rnorm(130, mean =  1 / 6, sd = 1),
  rnorm( 90, mean =  1 / 6, sd = 1),
  rnorm( 10, mean = -2 / 6, sd = 1)
)
g <- rep(1:3, c(130, 90, 10))
boxplot(
  x ~ g,
  width = c(length(g[g==1]), length(g[g==2]), length(g[g==3])),
  main  = "Dominant model (sample data)",
  xlab  = "Genotype", ylab="PK parameter"
)

# coefficient matrix
# c_1: additive, c_2: recessive, c_3: dominant
contrast <- rbind(
  c(-1, 0, 1), c(-2, 1, 1), c(-1, -1, 2)
)
y <- mmcm.resamp(x, g, contrast, nsample = 20000, abseps = 0.01, seed = 5784324)
y

## Example 2 ##
#  for dataframe
#  true response pattern: pos = 1 dominant  model c=( 1,  1, -2)
#                               2 additive  model c=(-1,  0,  1)
#                               3 recessive model c=( 2, -1, -1)
set.seed(3872435)
x <- c(
  rnorm(130, mean =  1 / 6, sd = 1),
  rnorm( 90, mean =  1 / 6, sd = 1),
  rnorm( 10, mean = -2 / 6, sd = 1),
  rnorm(130, mean = -1 / 4, sd = 1),
  rnorm( 90, mean =  0 / 4, sd = 1),
  rnorm( 10, mean =  1 / 4, sd = 1),
  rnorm(130, mean =  2 / 6, sd = 1),
  rnorm( 90, mean = -1 / 6, sd = 1),
  rnorm( 10, mean = -1 / 6, sd = 1)
)
g   <- rep(rep(1:3, c(130, 90, 10)), 3)
pos <- rep(c("rsXXXX", "rsYYYY", "rsZZZZ"), each=230)
xx  <- data.frame(pos = pos, x = x, g = g)

# coefficient matrix
# c_1: additive, c_2: recessive, c_3: dominant
contrast <- rbind(
  c(-1, 0, 1), c(-2, 1, 1), c(-1, -1, 2)
)

mmcmtapply <- function(r) {
  mmcm.resamp(
    xx$x[xx$pos==r[1]], xx$g[xx$pos==r[1]],
    contrast, nsample = 10000, abseps = 0.01, seed = 5784324+as.numeric(r[1])
  )
}
y <- tapply(xx$pos, xx$pos, mmcmtapply)
yy <- data.frame(
  Pos      = as.vector(names(y)),
  Pval     = as.vector(sapply(y, "[[", 3)),
  Pattern  = as.vector(sapply(y, "[[", 7)),
  MC_Error = as.vector(sapply(y, "[[", 9))
)
yy

mmcm documentation built on May 29, 2017, 3:31 p.m.