swh:1:snp:ff0951ca787d0b7f47dc2335f47fed43820a6324
Raw File
Tip revision: baa3a0f39e4ea7b973c6f69941906bb67bdcdd83 authored by Venkatraman E. Seshan on 24 March 2011, 00:00:00 UTC
version 0.9.6
Tip revision: baa3a0f
lehmann.R
power.ladesign <- function(gsize, odds.ratio, sig.level = 0.05, statistic = c("Kruskal-Wallis", "Jonckheere-Terpstra"), alternative = c("two.sided", "one.sided"), nrep=1e+6) {
  ng <- length(gsize)
  nOR <- length(odds.ratio)
  if (nOR != ng-1) {
    errmesg <- paste(" There are", ng, "groups and", nOR,"odds ratios, there should be", ng-1)
    stop(message=errmesg)
  }
  statistic <- match.arg(statistic)
  kw <- statistic=="Kruskal-Wallis"
  alternative <- match.arg(alternative)
  if (kw & alternative == "one.sided") stop("Only two-sided alternative is allowed for Kruskal-Wallis test")
  ng <- length(gsize)
  zz <- .Fortran("lehman",
                 as.integer(ng),
                 as.integer(gsize),
                 double(ng),
                 as.double(rep(1, ng)),
                 as.double(sum(gsize)),
                 double(ng),
                 as.logical(kw),
                 as.integer(nrep),
                 tstat=double(nrep))
  if (kw) {
    ldq <- quantile(zz$tstat, 1-sig.level)
  } else {
    if (alternative=="one.sided") {
      ldq <- quantile(zz$tstat, c(sig.level, 1-sig.level))
    } else {
      ldq <- quantile(zz$tstat, c(sig.level/2, 1-sig.level/2))
    }
  }
  zz <- .Fortran("lehman",
                 as.integer(ng),
                 as.integer(gsize),
                 double(ng),
                 as.double(c(1,odds.ratio)),
                 as.double(sum(gsize*c(1,odds.ratio))),
                 double(ng),
                 as.logical(kw),                 
                 as.integer(nrep),
                 tstat=double(nrep))
  if (kw) {
    ldpow <- sum(zz$tstat >= ldq)/nrep
  } else {
    if (alternative=="one.sided") {
      ldpow <- max(sum(zz$tstat <= ldq[1]), sum(zz$tstat >= ldq[2]))/nrep
    } else {
      ldpow <- (sum(zz$tstat <= ldq[1]) + sum(zz$tstat >= ldq[2]))/nrep
    }
  }
  out <- list()
  out$group.size <- gsize
  out$odds.ratio <- odds.ratio
  out$statistic <- statistic
  out$alternative <- alternative
  out$sig.level <- sig.level
  out$power <- ldpow
  class(out) <- "ladesign"
  out
}

print.ladesign <- function(x, ...) {
  cat("         Number of groups =", length(x$group.size), "\n")
  cat("               group size =", x$group.size, "\n")
  cat("odds ratios w.r.t group 1 =", round(x$odds.ratio, 3), "\n")
  cat("           test statistic =", x$statistic, "\n")
  cat("              alternative =", x$alternative, "\n")
  cat("       significance level =", x$sig.level, "\n")
  cat("                    power =", round(x$power, 3), "\n")
  invisible(x)
}
back to top