## 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 (Werte 0 und 1)
##               statistic : Teststatistik
##                       R : Anzahl Replikationen
##             alternative : gibt Richtung der Alternative
##                           ("greater", "less" oder "two.sided")
##                           greater: x[g==0] > x[g==1]
## -------------------------------------------------------------------------
## Author: Christian Keller, Date: 22 Aug 96
## Update: R. Frisullo: 3. Juli 02
## Update: B. Jaggi, C. Buser: 15. September 04
## Update: C. Birrer, 27. Oktober 06
    
        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

        ## alte Berechnung der P-Werte
        ##-         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)

        ## neue Berechnung der P-Werte, 22.9.04
        t.wosided <- sum(abs(t.star) >= abs(t.orig))/length(t.star)
        p.value <- switch(alternative,
                two.sided = t.wosided,
                greater = ifelse(t.orig >= 0, t.wosided/2, 1-t.wosided/2),
                less = ifelse(t.orig <= 0, t.wosided/2, 1-t.wosided/2))
        z <- list(t.star = t.star, t.orig = t.orig, p.value = p.value,
                      alternative = alternative)
        attr(z, "test") <- "Randomisation"
        class(z) <- "f.RandBootTest"
        z
}

print.f.RandBootTest <- function(x, samp.dist = FALSE, ...)
## Purpose:    Print-Methode für f.RandTest und f.BootTest
## -------------------------------------------------------------------------
## Arguments:  x : Resultat von f.RandTest oder f.BootTest
## -------------------------------------------------------------------------
## Author:     C. Birrer, 27. Oktober 06
{
  cat("\n       ", attr(x, "test"), " Test\n", sep = "")
  if (samp.dist) {
    cat("\nSampling Distribution:\n")
    print(structure(quantile(x$t.star),
                    names = c("Min", "1Q", "Median", "3Q", "Max")))
  }
  cat("\nTest statistic of original Data: ",
      x$t.orig, ", p-value: ", x$p.value, "\nalternative hypothesis: ",
      x$alternative,"\n", sep = "")
  invisible(x)
} 

summary.f.RandBootTest <- function(x, ...)
## Purpose: summary-Methode für f.RandTest und f.BootTest (ruft nur print auf)
## -------------------------------------------------------------------------
## Arguments:  x : Resultat von f.RandTest oder f.BootTest
## -------------------------------------------------------------------------
## Author:     C. Birrer, 27. Oktober 06
{
  print(x, samp.dist = TRUE)
}

plot.f.RandBootTest <- function(x, ...)
## Purpose:    Plot-Methode für f.RandTest und f.BootTest
## -------------------------------------------------------------------------
## Arguments:  x : Resultat von f.RandTest oder f.BootTest
## -------------------------------------------------------------------------
## Author:     C. Birrer, 27. Oktober 06
{
  hist(x$t.star, cex.lab = 1.4, freq = FALSE, main = "Sampling Distribution",
       xlab = "x")
  abline(v = x$t.orig, col = "red")
  axis(1, at = x$t.orig, padj = -1.1, col = 2)
  invisible()
}

f.BootTest <- function(x, g, statistic, R,
                       alternative = c("two.sided", "less", "greater"))
{
## Purpose:    Bootstrap-Test
## -------------------------------------------------------------------------
## Arguments:            x : Daten
##                       g : Gruppencode (Werte 0 und 1)
##               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
## Update: C. Birrer, 27. Oktober 06

        library(boot)
        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)
        z <- list(t.star = t.star, t.orig = t.orig, p.value = p.value,
                  alternative = alternative)
        attr(z, "test") <- "Bootstrap"
        class(z) <- "f.RandBootTest"
        z
}

f.HodgesLehmann <- function(data)
{
  ## Purpose:   Hodges-Lehmann estimate and confidence interval
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## Remark: function changed so that CI covers >= 95%, before it was too
  ##         small (9/22/04)
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 12 Aug 2002, 14:13
  ## Update: Beat Jaggi, Date: 22 Sept 2004  
  .cexact <-
    # c(NA,NA,NA,NA,NA,21,26,33,40,47,56,65,74,84,95,107,119,131,144,158)
      c(NA,NA,NA,NA,NA,22,27,34,41,48,57,66,75,85,96,108,120,132,145,159)  
  .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
}