Description Usage Arguments Details Value Author(s) References Examples
Density, distribution function and quantile function for the distribution of one and two sample permutation tests using the ShiftAlgorithm by Streitberg & R\"ohmel.
1 2 3 4 5 6 7 8  dperm(x, scores, m, paired=NULL, tol = 0.01, fact=NULL, density=FALSE,
simulate=FALSE, B=10000)
pperm(q, scores, m, paired=NULL, tol = 0.01, fact=NULL,
alternative=c("less", "greater", "two.sided"), pprob=FALSE,
simulate=FALSE, B=10000)
qperm(p, scores, m, paired=NULL, tol = 0.01, fact=NULL,
simulate=FALSE, B=10000)
rperm(n, scores, m)

x, q 
vector of quantiles. 
p 
vector of probabilities. 
scores 
arbitrary scores of the observations
of the 
m 
sample size of the 
paired 
logical. Indicates if paired observations are used. Needed
to discriminate between a paired problem and the distribution
of the total sum of the scores (which has mass 1 at the
point 
.
tol 
real. Real valued scores are mapped into integers by rounding
after multiplication with an appropriate factor.
Make sure that the absolute difference between the
each possible test statistic for the original scores and the
rounded scores is less than 
fact 
real. If 
n 
number of random observations to generate. 
alternative 
character indicating whether the probability
P(T ≤ q) ( 
pprob 
logical. Indicates if the probability P(T = q) should be computed additionally. 
density 
logical. When 
simulate 
logical. Use conditional MonteCarlo to compute the distribution. 
B 
number of MonteCarlo replications to be used. 
The exact distribution of the sum of the first m
scores is
evaluated using the ShiftAlgorithm by Streitberg & R\"ohmel under the
hypothesis of exchangeability (or, equivalent, the hypothesis that all
permutations of the scores are equally likely). The algorithm is able
to deal with tied scores, so the conditional distribution can be
evaluated.
The algorithm is defined for positive integer valued scores only.
There are two ways dealing with real valued scores.
First, one can try to find integer valued scores that lead to statistics
which differ not more than tol
from the statistics computed for the original scores. This can be done as
follows.
Without loss of generality let a_i > 0 denote real valued scores in
reverse ordering and
f a positive factor (this is the fact
argument).
Let R_i = f \cdot a_i  round(f \cdot a_i). Then
∑_{i=1}^m f \cdot a_i = ∑_{i=1}^m round(f \cdot a_i)  R_i.
Clearly, the maximum difference between 1/f ∑_{i=1}^m f \cdot a_i and 1/f ∑_{i=1}^n round(f \cdot a_i) is given by ∑_{i=1}^m R_i. Therefore one searches for f with
∑_{i=1}^m R_i ≤ ∑_{i=1}^m R_i ≤ tol.
If f induces more that 100.000 columns in the ShiftAlgorithm by Streitberg & R\"ohmel, f is restricted to the largest integer that does not.
The second idea is to map the scores into integers by taking the
integer part of a_i N / \max(a_i) (Hothorn & Lausen, 2002).
This induces additional ties, but the shape of the
scores is very similar. That means we do not try to approximate something
but use a different test (with integer valued scores), serving for the
same purpose
(due to a similar shape of the scores). However, this has to be done prior
to calling pperm
(see the examples).
Exact twosided pvalues are computed as suggested in the StatXact5 manual, page 225, equation (9.31) and equation (8.18), p. 179 (paired case). In detail: For the paired case the twosided pvalue is just twice the onesided one. For the independent sample case the two sided pvalue is defined as
p_2 = P( T  E(T) ≥  q  E(T) )
where q is the quantile passed to pperm
.
dperm
gives the density, pperm
gives the distribution
function and qperm
gives the quantile function. If pprob
is
true, pperm
returns a list with elements
PVALUE 
the probability specified by 
PPROB 
the probability P(T = q). 
rperm
is a wrapper to sample
.
Torsten Hothorn <Torsten.Hothorn@rzmail.unierlangen.de>
Bernd Streitberg & Joachim R\"ohmel (1986), Exact distributions for permutations and rank tests: An introduction to some recently published algorithms. Statistical Software Newsletter 12(1), 10–17.
Bernd Streitberg & Joachim R\"ohmel (1987), Exakte Verteilungen f\"ur Rang und Randomisierungstests im allgemeinen $c$Stichprobenfall. EDV in Medizin und Biologie 18(1), 12–19.
Torsten Hothorn (2001), On exact rank tests in R. R News 1(1), 11–12.
Cyrus R. Mehta & Nitin R. Patel (2001), StatXact5 for Windows. Manual, Cytel Software Cooperation, Cambridge, USA
Torsten Hothorn & Berthold Lausen (2003), On the exact distribution of maximally selected rank statistics. Computational Statistics \& Data Analysis, 43(2), 121137.
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 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81  # exact onesided pvalue of the Wilcoxon test for a tied sample
x < c(0.5, 0.5, 0.6, 0.6, 0.7, 0.8, 0.9)
y < c(0.5, 1.0, 1.2, 1.2, 1.4, 1.5, 1.9, 2.0)
r < cscores(c(x,y), type="Wilcoxon")
pperm(sum(r[seq(along=x)]), r, 7)
# Compare the exact algorithm as implemented in ctest and the
# ShiftAlgorithm by Streitberg & Roehmel for untied samples
# Wilcoxon:
n < 10
x < rnorm(n, 2)
y < rnorm(n, 3)
r < cscores(c(x,y), type="Wilcoxon")
# exact distribution using the ShiftAlgorithm
dwexac < dperm((n*(n+1)/2):(n^2 + n*(n+1)/2), r, n)
sum(dwexac) # should be something near 1 :)
# exact distribution using dwilcox
dw < dwilcox(0:(n^2), n, n)
# compare the two distributions:
plot(dw, dwexac, main="Wilcoxon", xlab="dwilcox", ylab="dperm")
# should give a "perfect" line
# Wilcoxon signed rank test
n < 10
x < rnorm(n, 5)
y < rnorm(n, 5)
r < cscores(abs(x  y), type="Wilcoxon")
pperm(sum(r[x  y > 0]), r, length(r))
wilcox.test(x,y, paired=TRUE, alternative="less")
psignrank(sum(r[x  y > 0]), length(r))
# AnsariBradley
n < 10
x < rnorm(n, 2, 1)
y < rnorm(n, 2, 2)
# exact distribution using the ShiftAlgorithm
sc < cscores(c(x,y), type="Ansari")
dabexac < dperm(0:(n*(2*n+1)/2), sc, n)
sum(dabexac)
# real scores are allowed (but only result in an approximation)
# e.g. v.d. Waerden test
n < 10
x < rnorm(n)
y < rnorm(n)
scores < cscores(c(x,y), type="NormalQuantile")
X < sum(scores[seq(along=x)]) # < v.d. Waerden normal quantile statistic
# critical value, twosided test
abs(qperm(0.025, scores, length(x)))
# pvalues
p1 < pperm(X, scores, length(x), alternative="two.sided")
# generate integer valued scores with the same shape as normal quantile
# scores, this no longer v.d.Waerden, but something very similar
scores < cscores(c(x,y), type="NormalQuantile", int=TRUE)
X < sum(scores[seq(along=x)])
p2 < pperm(X, scores, length(x), alternative="two.sided")
# compare p1 and p2
p1  p2

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.