## R
##
## R-Funktionen fuer Randomisierungs- und Bootstrap-Test
## -----------------------------------------------------

f.RandTest <- function(x, g, statistic, R = NULL,
                       alternative = c("two.sided", "less", "greater"))
{
## Purpose:    Randomisierungstest 
## -------------------------------------------------------------------------
## Arguments:            x : Daten
##                       g : Gruppencode
##               statistic : Teststatistik
##                       R : Anzahl Replikationen
##             alternative : gibt Richtung der Alternative
##                           ("greater", "less" oder "two.sided")
## -------------------------------------------------------------------------
## Author: Christian Keller, Date: 22 Aug 96
## Update: R. Frisullo: 3. Juli 02
  
        N <- length(x)
        alternative <- match.arg(alternative)
        if(!is.numeric(g) || length(g) != N || any(sort(unique(g)) != c(0, 1)))
                stop("ARG 'g' should be a vector of 0 and 1")
        if(!is.function(statistic))
                stop("ARG 'statistic' must be a function")
        foo <- function(object, code, funkt)
        {
                funkt(object[code == 0]) - funkt(object[code == 1])
        }
        
        # Wert der Teststatistik fuer Original-Stichprobe
        t.orig <- foo(object = x, code = g, funkt = statistic)
        if(is.null(R)) {
           # exakter Randomisierungstest 
                grid <- expand.grid(rep(list(c(0, 1)), N))
                anz <- apply(grid, 1, sum)
                mat <- grid[anz == sum(g),  ]
                t.star <- apply(mat, 1, foo, object = x, funkt = statistic)
        }
        else {
           # approximativer Randomisierungstest mit R Auswahlen
                if(!is.numeric(R) || length(R) > 1)
                  stop("ARG 'R' must be a number")
                t.star <- numeric(R)
                for(b in 1:R) {
                        x.star <- sample(x, N)
                        t.star[b] <- foo(object = x.star, code = g, funkt = 
                                statistic)
                }
        }
        names(t.star) <- NULL   
        p.value <- switch(alternative,
                greater = sum(t.star >= t.orig)/length(t.star),
                less = sum(t.star <= t.orig)/length(t.star),
                two.sided = sum(abs(t.star) >= abs(t.orig))/length(t.star))
        return(list(t.star = t.star, t.orig = t.orig, p.value = p.value,
                    alternative = alternative))
}


f.BootTest <- function(x, g, statistic, R,
                       alternative = c("two.sided", "less", "greater"))
{
## Purpose:    Bootstrap-Test
## -------------------------------------------------------------------------
## Arguments:            x : Daten
##                       g : Gruppencode
##               statistic : Teststatistik
##                       R : Anzahl Replikationen
##             alternative : gibt Richtung der Alternative
## -------------------------------------------------------------------------
## Author: Christian Keller, Date: 22 Aug 96
## Anpassung an R: R. Frisullo: 3. Juli 02

        N <- length(x)
        alternative <- match.arg(alternative)
        if(!is.numeric(g) || length(g) != N || any(sort(unique(g)) != c(0, 1))) 
                stop("ARG 'g' should be a vector of 0 and 1")
        if(!is.numeric(R) || length(R) > 1)
                stop("ARG 'R' must be a number")
        if(!is.function(statistic))
                stop("ARG 'statistic' must be a function")
        foo <- function(object, code, funkt)
        {
                funkt(object[code == 0]) - funkt(object[code == 1])
        }
        # Wert der Teststatistik fuer Original-Stichprobe
        t.orig <- foo(code = g, object = x, funkt = statistic)  
        # R Bootstrap-Replikationen
        f.BootStat <- function(daten, index, g, stat)
          {     foo(daten[index], code = sort(g, decreasing = TRUE),
                    funkt = stat)
          }
        t.star <- boot(x, statistic = f.BootStat, R, g=g, stat=statistic)$t
        p.value <- switch(alternative,
                greater = sum(t.star >= t.orig)/R,
                less = sum(t.star <= t.orig)/R,
                two.sided = sum(abs(t.star) >= abs(t.orig))/R)
        return(list(t.star = t.star, t.orig = t.orig, p.value = p.value,
                    alternative = alternative))
}

f.HodgesLehmann <- function(data)
{
  ## Purpose:   Hodges-Lehmann estimate and confidence interval
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 12 Aug 2002, 14:13
  .cexact <-
    c(NA,NA,NA,NA,NA,21,26,33,40,47,56,65,74,84,95,106,118,131,144,158)
  .d <- na.omit(data)
  .n <- length(.d)
  .wa <- sort(c(outer(.d,.d,"+")/2)[outer(1:.n,1:.n,"<=")])
  .c <- if (.n<=length(.cexact)) .n*(.n+1)/2+1-.cexact[.n] else
    floor(.n*(.n+1)/4-1.96*sqrt(.n*(.n+1)*(2*.n+1)/24))
  .r <- c(median(.wa), .wa[c(.c,.n*(.n+1)/2+1-.c)])
  names(.r) <- c("estimate","lower","upper")
  .r
}