Supplement 2 for: A new u-statistic with superior
design sensitivity in matched observational studies.

This supplement contains an R session and R code.  It
is provided without guarantees of any kind.

These are the 21 matched pair differences for the 
methotrexate example in section 1.2 of the paper.

> dif
 [1]  0.03 -0.05 -0.02  0.06  0.00  0.05  0.03  0.08 -0.08
[10] -0.06  0.16 -0.03  0.45  0.08  0.25 -0.64  0.32  0.12
[19]  0.38  0.14  0.10

The function uscore computes the ranks for matched pair
differences y with the suggested adjustement for ties
in the web based supplement.

> uscore
function(y,m=5,mb=4,mB=5){
        I<-length(y)
        q<-rep(0,I)
        a<-rank(abs(y))
        for (l in mb:mB){
                q<-q+choose(a-1,l-1)*choose(I-a,m-l)
                }
        q*(abs(y)>0)
}

The following computes the rank scores with non-monotone 
ranks (m,mb,mB) = (8,6,7) as in section 5.2 of the paper.

> q<-uscore(dif,m=8,mb=6,mB=7)

These are the rank scores scaled to have maximum 1 to aid
viewing.  Notice that the largest rank 1.000 is awarded to a 
moderately large difference of 0.25.  In contrast, with monotone 
scores, mB=m, the largest |Y| receives the largest rank.
This is an illustration.  Use monotone scores except with 
long-tailed distributions.

> round(q/max(q),3)
 [1] 0.000 0.005 0.000 0.047 0.000 0.005 0.000 0.231 0.231
[10] 0.047 0.945 0.000 0.466 0.231 1.000 0.000 0.956 0.670
[19] 0.785 0.825 0.508

The following function uses expression (3) in section 2.2
to compute the approximate upper bound on the one sided 
P-value.

> approx
function(y,m=5,mb=4,mB=5,gamma=1){
        rk<-uscore(y,m=m,mb=mb,mB=mB)*(abs(y)>0)
        T<-sum((y>0)*rk)
        ET<-sum(rk)*gamma/(1+gamma)
        VT<-sum(rk*rk)*gamma/((1+gamma)^2)
        dev<-(T-ET)/sqrt(VT)
        1-pnorm(dev)    
}

> approx(dif,m=8,mb=6,mB=7,gamma=2)
[1] 0.02910648

The value above appears in Table 2 at gamma=2, 
(8,6,7).

The functions exact and lambda compute the exact upper
bound on the P-value for small sample sizes I.  The
function exact calls lambda, and the function lambda is 
defined in section 5.2.  The function lambda is a 
recursive function: it calls itself.

> exact
function(y,m=5,mb=4,mB=5,gamma=1){
        I<-length(y)
        y<-y[order(abs(y))]
        rk<-uscore(y,m=m,mb=mb,mB=mB)*(abs(y)>0)
        T<-sum((y>0)*rk)
        lambda(I,T,rk,gamma=gamma)
}

> lambda
function(I,c,rk,gamma=1){
        #calculates Pr(T>=c) where T is a signed rank statistic with ranks in the
        #vector rk, where the user calls lambda with I=length(rk), and gamma is
        #the sensitivity parameter
        kappa<-gamma/(1+gamma)
        if (I==length(rk)) rk<-sort(rk)
        if (c<=0) 1
                else if (I<=0) 0
                else lambda(I-1,c-rk[I],rk,gamma=gamma)*kappa+lambda(I-1,c,rk,gamma=gamma)*(1-kappa)
}

The following computes the exact upper bound on the P-value
for gamma=2 using (8,6,7) as reported in section 5.2

> exact(dif,m=8,mb=6,mB=7,gamma=2)
[1] 0.02611934