WilcoxonSRPVs <- function(X) # WilcoxonSRPVs(X) # computes for a data vector X (vector of differences) # Wilcoxon's signed-rank test statistic # T = sum_{i=1}^n sign(X[i]) R[i] , # where R(i) denotes the rank of |X[i]| among the non-zero # components of (|X[1]|, |X[2]|, ..., |X[n]|). # # In case of X[i] = 0 we set R[i] = 0, i.e. components # equal to zero are ignored. # # In case of ties among the absolute values of X's # non-zero components, we use the usual convention of # average ranks. # # In addition, the left-, right- and two-sided P-values # pvl = Pr(T_* <= T) , # pvr = Pr(T_* >= T) , # pvz = 2 * min(pvl, pvr) # are computed. Here T_* = sum_{i=1}^n S[i]*R[i] with # independent random signs S[1], S[2], ..., S[n], while # R is viewed temporarily as a fixed vector. # # Additional output arguments are # - supp : a vector containing all support points of the # distribution function F(.) := Pr(T_* <= .), # - cdf : a vector containing the numbers cdf(i) = F(supp[i]). # # Lutz Duembgen, May 20, 2011 { R <- abs(X) tmp <- (R > 0) R[tmp] <- rank(R[tmp]) T <- sum(sign(X)*R) N <- sum(tmp) R2 <- ceiling(2*sort(R[tmp])) h2 <- N*(N+1) h1 <- ceiling(h2/2) # Now we work temporarily with the test statistic # T2 = T + N(N+1)/2 # = sum_{i=1}^N 1{X[i] > 0} 2R[i] # taking values in {0,1,...,N(N+1)}: cdf <- rep(1,h2+1) m <- 1 for (i in 1:N) { m_new <- m+R2[i] cdf[(R2[i]+1):m_new] <- (cdf[1:m] + cdf[(R2[i]+1):m_new])/2 cdf[1:R2[i]] <- cdf[1:R2[i]]/2 m <- m_new } supp <- (0:h2) - h1 ind <- h1 + T + 1 pv.left <- cdf[ind] if (ind >= 2) { pv.right <- 1 - cdf[ind-1] } else { pv.right <- 1 } pv.twosided <- 2*min(pv.left,pv.right) return(list(T=T, pv.left=pv.left,pv.right=pv.right,pv.twosided=pv.twosided, supp=supp,cdf=cdf)) }