Here I provide a exact one sample or paired sample test for mean(s), which is the exact version of permutation test to compare a single sample against a mean or to compare paired samples’ differences against a mean.
The algorithm in this function is based on Bryan F. J. Manly. 1997.
Randomization, bootstrap and Monte Carlo methods in biology. 2nd edition. pp. 91-97.
That is, all possible permutations (i.e. all possible cases of exchange between $x_i$ and $\mu$ in one-sample test or all possible cases of exchange between $x_{1i}$ and $x_{2i}$) in paired-sample test are processed to calculate a exact p-value.
R code
Arguments
The usage of this function exactOneOrPairedSampleTest()
is very similar to the R built-in function t.test()
. There are four arguments:
x1
: a numeric vector specifying $x_{1i}$
x2 = NULL
: a numeric vector specifying $x_{2i}$ (in paired-sample case)
mu = 0
: a number specifying true mean or true difference $\mu$
alternative = c("t","g","l")
a single character specifying $H_0$: true mean of $x_1$ or mean of $x_{1i} - x_{2i}$ is equal, greater or less than $\mu$, respectively.
Example 1: one-sample test
Let $x_i = \{43,67,64,64,51,53,53,26,36,48,34,48,6 \}$, $i=1 \ldots 13$:
> x <- c(43,67,64,64,51,53,53,26,36,48,34,48,6)
Now we can call the above function exactOneOrPairedSampleTest()
to test $H_0$: true mean of $x_i = 56$ against $H_A$: true mean of $x_i \neq 56$:
> test1 <- exactOneOrPairedSampleTest(x, alternative="t", mu=56)
> test1
Exact one sample test
Alternative hypothesis: true mean is not equal to 56
mean(x) - mu = -10.3846153846154
Number of total permutation = 8192
Number of rejected permutation = 364
P-value = 0.04443359
The results show that:
- $\sum_{i=1}^{13} x_i / 13 - \mu = 45.61538 - 56 = -10.38462$;
- Totally $8192 = 2^{13}$ permutations are processed;
- 364 permutations do not support $H_0$;
- $P = 364/2^{13} = 0.0444$.
We can also call the function hist.exactOneOrPairedSampleTest()
:
> hist(test1)
As you can see, the vertical lines $x=\pm 10.38462$ shows the critical boundary of rejecting $H_0$. Note that it is a two-tail test.
Finally, you may call the returned object:
> str(test1)
List of 14
$ x1 : num [1:13] 43 67 64 64 51 53 53 26 36 48 ...
$ x2 : NULL
$ x1.name : chr "x"
$ x2.name : chr "NULL"
$ n : int 13
$ mu : num 56
$ test.0 : num -10.4
$ is.onesample: logi TRUE
$ alternative : chr "t"
$ test.perm : num [1:8192] -10.38 -8.38 -12.08 -10.08 -11.62 ...
$ DF :'data.frame': 13 obs. of 3 variables:
..$ x1 : num [1:13] 43 67 64 64 51 53 53 26 36 48 ...
..$ mu : num [1:13] 56 56 56 56 56 56 56 56 56 56 ...
..$ diff: num [1:13] -13 11 8 8 -5 -3 -3 -30 -20 -8 ...
$ N : num 8192
$ p.value : num 0.0444
$ rejected.N : int 364
- attr(*, "class")= chr "exactOneOrPairedSampleTest"
to get more details of the results if you need them.
Example 2: paired-sample test
Let $x_{1i} = \{92, 0,72,80,57,76,81,67,50,77,90\}$ and $x_{2i} = \{43,67,64,64,51,53,53,26,36,48,34\}$, where each $i$ is paired and $i=1 \ldots 11$:
> x1 <- c(92, 0,72,80,57,76,81,67,50,77,90)
> x2 <- c(43,67,64,64,51,53,53,26,36,48,34)
Now we can call the function exactOneOrPairedSampleTest()
to test $H_0$: true mean of $x_{1i} - x_{2i} \leq 10$ against $H_A$: true mean of $x_{1i} - x_{2i} > 10$:
> test2 <- exactOneOrPairedSampleTest(x1, x2, alternative="g", mu=10)
> # equivalent to test2 <- exactOneOrPairedSampleTest(x1 - x2, alternative="g", mu=10)
> test2
Exact paired sample test
Alternative hypothesis: means of x1 - x2 is greater than 10
mean((x1)-(x2)) - mu = 8.45454545454546
Number of total permutation = 2048
Number of rejected permutation = 445
P-value = 0.2172852
The results show that:
- $\sum_{i=1}^{11} (x_{1i}-x_{2i}) / 11 - 10 = 8.45455$;
- Totally $2048 = 2^{11}$ permutations are processed;
- 445 permutations do not support $H_0$;
- $P = 445/2^{11} = 0.2172852$.
We can also call the function hist.exactOneOrPairedSampleTest()
:
> hist(test2)
As you can see, the vertical lines $x=8.45455$ shows the critical boundary of rejecting $H_0$. Note that this is a right-tail test.
Again, you may call the returned object:
> str(test2)
List of 14
$ x1 : num [1:11] 92 0 72 80 57 76 81 67 50 77 ...
$ x2 : num [1:11] 43 67 64 64 51 53 53 26 36 48 ...
$ x1.name : chr "x1"
$ x2.name : chr "x2"
$ n : int 11
$ mu : num 10
$ test.0 : num 8.45
$ is.onesample: logi FALSE
$ alternative : chr "g"
$ test.perm : num [1:2048] 8.45 1.36 22.45 15.36 8.82 ...
$ DF :'data.frame': 11 obs. of 3 variables:
..$ x1 : num [1:11] 92 0 72 80 57 76 81 67 50 77 ...
..$ x2 : num [1:11] 43 67 64 64 51 53 53 26 36 48 ...
..$ diff: num [1:11] 39 -77 -2 6 -4 13 18 31 4 19 ...
$ N : num 2048
$ p.value : num 0.217
$ rejected.N : int 445
- attr(*, "class")= chr "exactOneOrPairedSampleTest"
to get more details of the results if you need them.