home *** CD-ROM | disk | FTP | other *** search
/ Big Green CD 8 / BGCD_8_Dev.iso / NEXTSTEP / UNIX / Educational / R-0.49-MI / R-0.49-I / library / ctest < prev    next >
Encoding:
Text File  |  1997-09-14  |  24.4 KB  |  843 lines
  1. bartlett.test <- function(x, g)
  2. {
  3.   if (is.list(x)) {
  4.     if (length(x) < 2)
  5.       stop("x must be a list with at least 2 elements")
  6.     DNAME <- deparse(substitute(x))
  7.     x <- lapply(x, function(x) x <- x[is.finite(x)])
  8.     k <- length(x)
  9.   }
  10.   else {
  11.     if (length(x) != length(g))
  12.       stop("x and g must have the same length")
  13.     DNAME <- paste(deparse(substitute(x)), "and", deparse(substitute(g)))
  14.     OK <- complete.cases(x, g)
  15.     x <- x[OK]
  16.     g <- as.factor(g[OK])
  17.     k <- nlevels(g)
  18.     if (k < 2)
  19.       stop("all observations are in the same group")
  20.     x <- split(x, g)
  21.   }
  22.  
  23.   n <- sapply(x, "length") - 1
  24.   if (any(n <= 0))
  25.     stop("there must be at least 2 observations in each group")
  26.   v <- sapply(x, "var")
  27.   n.total <- sum(n)
  28.   v.total <- sum(n * v) / n.total
  29.   STATISTIC <- ((n.total * log(v.total) - sum(n * log(v))) /
  30.         (1 + (sum(1 / n) - 1 / n.total) / (3 * (k - 1))))
  31.   names(STATISTIC) <- "Bartlett's K-square"
  32.   PARAMETER <- k - 1
  33.   names(PARAMETER) <- "df"
  34.   
  35.   RVAL <- list(statistic = STATISTIC,
  36.            parameter = PARAMETER,
  37.            p.value = 1 - pchisq(STATISTIC, PARAMETER),
  38.            data.name = DNAME,
  39.            method = "Bartlett test for homogeneity of variances")
  40.   class(RVAL) <- "htest"
  41.   return(RVAL)
  42. }
  43. binom.test <- function(x, n, p = 0.5, alternative = "two.sided")
  44. {
  45.   if ((length(n) > 1) || is.na(n) || (n < 1) || (n != round(n)))
  46.     stop("n must be a positive integer")
  47.   if ((length(x) > 1) || is.na(x) || (x < 0) || (x > n) || (x != round(x)))
  48.     stop("x must be an integer between 0 and n")
  49.   if (!missing(p) && (length(p) > 1 || is.na(p) || p < 0 || p > 1))
  50.     stop ("p must be a single number between 0 and 1")
  51.  
  52.   CHOICES <- c("two.sided", "less", "greater")
  53.   alternative <- CHOICES[pmatch(alternative, CHOICES)]
  54.   if (length(alternative) > 1 || is.na(alternative))
  55.     stop ("alternative must be \"two.sided\", \"less\" or \"greater\"")
  56.  
  57.   DNAME <- paste(deparse(substitute(x)), "and", deparse(substitute(n)))
  58.   if (alternative == "less")
  59.     PVAL <- pbinom(x, n, p)
  60.   else if (alternative == "greater")
  61.     PVAL <- 1 - pbinom(x - 1, n, p)
  62.   else {
  63.     eps <- 10^(-6)
  64.     if (x/n < p) {
  65.       PVAL <- pbinom(x, n, p)
  66.       v <- 1 - pbinom(0:n, n, p)
  67.     }
  68.     else {
  69.       PVAL <- 1 - pbinom(x - 1, n, p)
  70.       v <- pbinom(0:n, n, p)
  71.     }
  72.     PVAL <- min(1, PVAL + max(v[v <= (1 + eps) * PVAL]))
  73.   }
  74.   
  75.   names(x) <- "number of successes"    # or simply "x" ??
  76.   names(n) <- "number of trials"    # or simply "n" ??
  77.   names(p) <- "probability of success"    # or simply "p" ??
  78.  
  79.   RVAL <- list(statistic = x,
  80.            parameter = n,
  81.            p.value = PVAL,
  82.            null.value = p,
  83.            alternative = alternative,
  84.            method = "Exact binomial test",
  85.            data.name = DNAME)
  86.   class(RVAL) <- "htest"
  87.   return(RVAL)
  88. }
  89. chisq.test <- function(x, y = NULL, correct = TRUE,
  90.                p = rep(1 / length(x), length(x)))
  91. {
  92.   DNAME <- deparse(substitute(x))
  93.   if (is.matrix(x)) {
  94.     if (min(dim(x)) == 1)
  95.       x <- as.vector(x)
  96.   }
  97.   if (!is.matrix(x) && !is.null(y)) {
  98.     if (length(x) != length(y))
  99.       stop("x and y must have the same length")
  100.     DNAME <- paste(DNAME, "and", deparse(substitute(y)))
  101.     OK <- complete.cases(x, y)
  102.     x <- as.factor(x[OK])
  103.     y <- as.factor(y[OK])
  104.     if ((nlevels(x) < 2) || (nlevels(y) < 2))
  105.       stop("x and y must have at least 2 levels")
  106.     x <- table(x, y)
  107.   }
  108.  
  109.   if (any(x < 0) || any(is.na(x)))
  110.     stop("all entries of x must be nonnegative and finite")
  111.  
  112.   if (is.matrix(x)) {
  113.     METHOD <- "Pearson's Chi-square test"    
  114.     E <- outer(apply(x, 1, sum), apply(x, 2, sum), "*") / sum(x)  
  115.     if (correct && nrow(x) == 2 && ncol(x) == 2) {
  116.       YATES <- .5
  117.       METHOD <- paste(METHOD, "with Yates' continuity correction")
  118.     }
  119.     else
  120.       YATES <- 0
  121.     dimnames(E) <- dimnames(x)
  122.     STATISTIC <- sum((abs(x - E) - YATES)^2 / E)
  123.     PARAMETER <- (nrow(x) - 1) * (ncol(x) - 1)
  124.   }
  125.   else {
  126.     if (length(x) == 1)
  127.       stop("x must at least have 2 elements")
  128.     if (length(x) != length(p))
  129.       stop("x and p must have the same number of elements")
  130.     METHOD <- "Chi-square test for given probabilities"
  131.     E <- sum(x) * p
  132.     names(E) <- names(x)
  133.     STATISTIC <- sum((x - E) ^ 2 / E)
  134.     PARAMETER <- length(x) - 1
  135.   }
  136.  
  137.   names(STATISTIC) <- "X-squared"
  138.   names(PARAMETER) <- "df"
  139.   if (any(E < 5))
  140.     warning("Chi-square approximation may be incorrect") 
  141.   PVAL <- 1 - pchisq(STATISTIC, PARAMETER)
  142.   RVAL <- list(statistic = STATISTIC,
  143.            parameter = PARAMETER,
  144.            p.value = PVAL,
  145.            method = METHOD,
  146.            data.name = DNAME,
  147.            observed = x,
  148.            expected = E)
  149.   class(RVAL) <- "htest"
  150.   return(RVAL)
  151. }
  152. cor.test <- function(x, y, alternative = "two.sided", method = "pearson") 
  153. {
  154.   CHOICES <- c("two.sided", "less", "greater")
  155.   alternative <- CHOICES[pmatch(alternative, CHOICES)]
  156.   if (length(alternative) > 1 || is.na(alternative)) 
  157.     stop("alternative must be \"two.sided\", \"less\" or \"greater\"")
  158.  
  159.   CHOICES <- c("pearson", "kendall", "spearman")
  160.   method <- CHOICES[pmatch(method, CHOICES)]
  161.   if (length(method) > 1 || is.na(method)) 
  162.     stop("method must be \"pearson\", \"kendall\" or \"spearman\"")
  163.  
  164.   DNAME <- paste(deparse(substitute(x)), "and", deparse(substitute(y)))
  165.  
  166.   if (length (x) != length (y))
  167.     stop ("x and y must have the same length")
  168.   OK <- complete.cases(x, y)
  169.   x <- x[OK]
  170.   n <- length(x)
  171.   if (n < 3)
  172.     stop("not enough finite observations")
  173.   else
  174.     y <- y[OK]
  175.  
  176.   NVAL <- 0
  177.  
  178.   if (method == "pearson") {
  179.     method <- "Pearson's product-moment correlation"
  180.     names(NVAL) <- "correlation"
  181.     r <- cor(x, y)
  182.     ESTIMATE <- r
  183.     names(ESTIMATE) <- "cor"
  184.     PARAMETER <- n - 2
  185.     names(PARAMETER) <- "df"
  186.     STATISTIC <- sqrt(PARAMETER) * r / sqrt(1 - r^2)
  187.     names(STATISTIC) <- "t"
  188.     PVAL <- pt(STATISTIC, PARAMETER)
  189.   }
  190.   else {
  191.     if (method == "kendall") {
  192.       method <- "Kendall's rank correlation tau"
  193.       names(NVAL) <- "tau"
  194.       x <- rank(x)
  195.       y <- rank(y)
  196.       ESTIMATE <- cor (c(sign(outer(x, x, "-"))),
  197.                c(sign(outer(y, y, "-"))))
  198.       names(ESTIMATE) <- "tau"
  199.       STATISTIC <- ESTIMATE / sqrt ((4*n + 10) / (9 * n * (n-1)))
  200.     }
  201.     else {
  202.       method <- "Spearman's rank correlation rho"
  203.       names(NVAL) <- "rho"
  204.       ESTIMATE <- cor (rank(x), rank(y))
  205.       names(ESTIMATE) <- "rho"
  206.       STATISTIC <- sqrt (n-1) * (ESTIMATE - 6 / (n^3 - n))
  207.     }
  208.     PARAMETER <- NULL
  209.     names(STATISTIC) <- "z"
  210.     PVAL <- pnorm(STATISTIC)
  211.   }
  212.  
  213.   if (alternative == "two.sided")
  214.     PVAL <- 2 * min (PVAL, 1 - PVAL)
  215.   else if (alternative == "greater")
  216.     PVAL <- 1 - PVAL
  217.  
  218.   RVAL <- list(statistic = STATISTIC,
  219.            parameter = PARAMETER,
  220.            p.value = PVAL,
  221.            estimate = ESTIMATE,
  222.            null.value = NVAL,
  223.            alternative = alternative,
  224.            method = method,
  225.            data.name = DNAME)
  226.   class(RVAL) <- "htest"
  227.   return(RVAL)
  228. }
  229. fisher.test <- function(x, y = NULL, alternative = "two.sided")
  230. {
  231.   DNAME <- deparse(substitute(x))  
  232.   if (is.matrix(x)) {
  233.     if (any(dim(x) < 2))
  234.       stop("x must have at least 2 rows and columns")
  235.     if (any(x < 0) || any(is.na(x))) 
  236.       stop("all entries of x must be nonnegative and finite")
  237.   }
  238.   else {
  239.     if (is.null(y)) 
  240.       stop("if x is not a matrix, y must be given")
  241.     if (length(x) != length(y)) 
  242.       stop("x and y must have the same length")
  243.     DNAME <- paste(DNAME, "and", deparse(substitute(y)))
  244.     OK <- complete.cases(x, y)
  245.     x <- as.factor(x[OK])
  246.     y <- as.factor(y[OK])
  247.     if ((nlevels(x) < 2) || (nlevels(y) < 2)) 
  248.       stop("x and y must have at least 2 levels")
  249.     x <- table(x, y)
  250.   }
  251.  
  252.   if (any(dim(x) != c(2, 2)))
  253.     stop("Sorry, only 2 by 2 tables are currently implemented")
  254.  
  255.   CHOICES <- c("two.sided", "less", "greater")
  256.   alternative <- CHOICES[pmatch(alternative, CHOICES)]
  257.   if (length(alternative) > 1 || is.na(alternative)) 
  258.     stop("alternative must be \"two.sided\", \"less\" or \"greater\"")  
  259.  
  260.   m <- sum(x[, 1])
  261.   n <- sum(x[, 2])
  262.   k <- sum(x[1, ])
  263.   x <- x[1, 1]
  264.   PVAL <- switch(alternative,
  265.          less = phyper(x, m, n, k),
  266.          greater = 1 - phyper(x - 1, m, n, k),
  267.          two.sided = {
  268.            eps <- 10^(-6)
  269.            if ((PVAL <- phyper(x, m, n, k)) < .5)
  270.              v <- phyper(0:k, m, n, k)
  271.            else {
  272.              v <- 1 - phyper(0:k, m, n, k)
  273.              PVAL <- 1 - phyper(x - 1, m, n, k)
  274.            }
  275.            min(1, PVAL + max(v[v <= (1 + eps) * PVAL]))
  276.          })
  277.   RVAL <- list(p.value = PVAL,
  278.            alternative = alternative,
  279.            method = "Fisher's Exact Test for Count Data",
  280.            data.name = DNAME)
  281.   class(RVAL) <- "htest"
  282.   return(RVAL)
  283. }
  284. friedman.test <- function(y, groups, blocks)
  285. {
  286.   DNAME <- deparse(substitute(y))
  287.   if (is.matrix(y)) {
  288.     groups <- as.factor(c(col(y)))
  289.     blocks <- as.factor(c(row(y)))
  290.   }
  291.   else {
  292.     if (any(is.na(c(groups, blocks))))
  293.       stop("NA's are not allowed in groups or blocks")
  294.     if (any(diff(c(length(y), length(groups), length(blocks)))))
  295.       stop("y, groups and blocks must have the same length")
  296.     DNAME <- paste(DNAME, ", ", deparse(substitute(groups)), " and ",
  297.            deparse(substitute(blocks)), sep = "")
  298.     if (any(table(groups, blocks) != 1))
  299.       stop("Not an unreplicated complete block design")
  300.     groups <- as.factor(groups)
  301.     blocks <- as.factor(blocks)
  302.   }
  303.  
  304.   k <- nlevels(groups)
  305.   y <- matrix(unlist(split(y, blocks)), ncol = k, byrow = T)
  306.   y <- y[complete.cases(y), ]
  307.   n <- nrow(y)
  308.   r <- t(apply(y, 1, rank))
  309.   TIES <- apply(r, 1, table)
  310.   STATISTIC <- ((12 * sum((apply(r, 2, sum) - n * (k + 1) / 2)^2)) /
  311.         (n * k * (k + 1)
  312.          - (sum(unlist(lapply(TIES, function (u) {u^3 - u}))) /
  313.             (k - 1))))
  314.   PARAMETER <- k - 1
  315.   names(STATISTIC) <- "Friedman chi-square"
  316.   names(PARAMETER) <- "df"
  317.   RVAL <- list(statistic = STATISTIC,
  318.            parameter = PARAMETER,
  319.            p.value = 1 - pchisq(STATISTIC, PARAMETER),
  320.            method = "Friedman rank sum test",
  321.            data.name = DNAME)
  322.   class(RVAL) <- "htest"
  323.   return(RVAL)
  324. }
  325. kruskal.test <- function(x, g)
  326. {
  327.   if (is.list(x)) {
  328.     if (length(x) < 2)
  329.       stop("x must be a list with at least 2 elements")
  330.     DNAME <- deparse(substitute(x))
  331.     x <- lapply(x, function(x) x <- x[is.finite(x)])
  332.     k <- length(x)
  333.     g <- as.factor(rep(1 : k, sapply(x, "length")))
  334.     x <- unlist(x)
  335.   }
  336.   else {
  337.     if (length(x) != length(g))
  338.       stop("x and g must have the same length")
  339.     DNAME <- paste(deparse(substitute(x)), "and", deparse(substitute(g)))
  340.     if (!all(is.finite(g)))
  341.       stop("all group levels must be finite")
  342.     OK <- !is.na(x)
  343.     x <- x[OK]
  344.     g <- as.factor(g[OK])
  345.     k <- nlevels(g)
  346.     if (k < 2)
  347.       stop("all observations are in the same group")
  348.   }
  349.  
  350.   n <- length(x)
  351.   if (n < 2)
  352.     stop("not enough observations")
  353.   r <- rank(x)
  354.   TIES <- table(x)
  355.   STATISTIC <- sum(tapply(r, g, "sum")^2 / tapply(r, g, "length"))
  356.   STATISTIC <- ((12 * STATISTIC / (n * (n + 1)) - 3 * (n + 1)) /
  357.         (1 - sum(TIES^3 - TIES) / (n^3 - n)))
  358.   names(STATISTIC) <- "Kruskal-Wallis chi-square"
  359.   PARAMETER <- k - 1
  360.   names(PARAMETER) <- "df"
  361.  
  362.   RVAL <- list(statistic = STATISTIC,
  363.            parameter = PARAMETER,
  364.            p.value = 1 - pchisq(STATISTIC, PARAMETER),
  365.            method = "Kruskal-Wallis rank sum test",
  366.            data.name = DNAME)
  367.   class(RVAL) <- "htest"
  368.   return(RVAL)
  369. }
  370. ks.test <- function(x, y, ..., alternative = "two.sided")
  371. {
  372.   CHOICES <- c("two.sided", "less", "greater")
  373.   alternative <- CHOICES[pmatch(alternative, CHOICES)]
  374.   if (length(alternative) > 1 || is.na(alternative)) 
  375.     stop("alternative must be \"two.sided\", \"less\" or \"greater\"")
  376.  
  377.   DNAME <- deparse(substitute(x))      
  378.   x <- x[!is.na(x)]
  379.   n <- length(x)
  380.   if (n < 1)
  381.     stop("Not enough x data")
  382.   
  383.   if (is.numeric(y)) {
  384.     DNAME <- paste(DNAME, "and", deparse(substitute(y)))
  385.     y <- y[!is.na(y)]
  386.     n.x <- n
  387.     n.y <- length(y)
  388.     if (n.y < 1)
  389.       stop("Not enough y data")
  390.     METHOD <- "Two-sample Kolmogorov-Smirnov test"
  391.     n <- n.x * n.y / (n.x + n.y)
  392.     z <- ifelse (order(c(x, y)) <= n.x, 1/n.x, -1/n.y)
  393.     STATISTIC <- switch(alternative,
  394.             "two.sided" = max(abs(cumsum(z))),
  395.             "greater" = max(cumsum(z)),
  396.             "less" = - min(cumsum(z)))
  397.   }
  398.   else {
  399.     if (is.character(y))
  400.       y <- get(y, mode="function")
  401.     if (mode(y) != "function")
  402.       stop("y must be numeric or a string naming a valid function")
  403.     METHOD <- "One-sample Kolmogorov-Smirnov test"
  404.     n <- length(x)
  405.     x <- y(sort(x), ...) - (0:(n-1))/n
  406.     STATISTIC <- switch(alternative,
  407.             "two.sided" = max(abs(c(x, x-1/n))),
  408.             "greater" = max(c(x, x-1/n)),
  409.             "less" = - min(c(x, x-1/n)))
  410.   }
  411.  
  412.   names(STATISTIC) <- switch(alternative,
  413.                  "two.sided" = "D",
  414.                  "greater" = "D^+",
  415.                  "less" = "D^-")
  416.   PVAL <- ifelse(alternative == "two.sided",
  417.          1 - pks(sqrt(n) * STATISTIC),
  418.          exp(- 2 * n * STATISTIC^2))
  419.   
  420.   RVAL <- list(statistic = STATISTIC,
  421.            p.value = PVAL,
  422.            alternative = alternative,
  423.            method = METHOD,
  424.            data.name = DNAME)
  425.   class(RVAL) <- "htest"
  426.   return(RVAL)
  427. }
  428.  
  429.  
  430. mantelhaen.test <- function(x, y = NULL, z = NULL, correct = TRUE)
  431. {
  432.   DNAME <- deparse(substitute(x))
  433.   if (is.array(x)) {
  434.     if (length(dim(x)) == 3) {
  435.       if (any(is.na(x)))
  436.     stop("NAs are not allowed")
  437.       if (dim(x)[1:2] != c(2, 2))
  438.     stop("table for each stratum must be 2 by 2")
  439.     }
  440.     else
  441.       stop("x must be a 3-dimensional array")
  442.   }
  443.   else {
  444.     if (is.null(y))
  445.       stop("If x is not an array, y must be given")
  446.     if (is.null(z))
  447.       stop("If x is not an array, z must be given")
  448.     if (any(diff(c(length(x), length(y), length(z)))))
  449.       stop("x, y, and z must have the same length")
  450.     DNAME <- paste(DNAME, "and", deparse(substitute(y)), "and",
  451.            deparse(substitute(z)))
  452.     OK <- complete.cases(x, y, z)
  453.     x <- as.factor(x[OK])
  454.     y <- as.factor(y[OK])
  455.     if ((nlevels(x) != 2) || (nlevels(y) != 2))
  456.       stop("x and y must be dichotomous")
  457.     else
  458.       x <- table(x, y, z[OK])
  459.   }
  460.  
  461.   s.x <- apply(x, c(1, 3), sum)
  462.   s.y <- apply(x, c(2, 3), sum)
  463.   n <- apply(x, 3, sum)
  464.   if (any(n < 2))
  465.     stop("sample size in each stratum must be > 1")
  466.   DELTA <- abs(sum(x[1, 1, ] - s.x[1, ] * s.y[1, ] / n))
  467.   YATES <- ifelse(correct && (DELTA >= .5), .5, 0)
  468.   STATISTIC <- ((DELTA - YATES)^2 /
  469.         sum(apply(rbind(s.x, s.y), 2, prod) / (n^2 * (n - 1))))
  470.   PARAMETER <- 1
  471.   names(STATISTIC) <- "Mantel-Haenszel X-square"
  472.   names(PARAMETER) <- "df"
  473.  
  474.   RVAL <- list(statistic = STATISTIC,
  475.            parameter = PARAMETER,
  476.            p.value = 1 - pchisq(STATISTIC, PARAMETER),
  477.            method = paste("Mantel-Haenszel chi-square test",
  478.            ifelse(YATES, "with", "without"),
  479.            "continuity correction"),
  480.            data.name = DNAME)
  481.   class(RVAL) <- "htest"
  482.   return(RVAL)
  483. }
  484. mcnemar.test <- function(x, y = NULL, correct = TRUE)
  485. {
  486.   if (is.matrix(x)) {
  487.     r <- nrow(x)
  488.     if ((r < 2) || (ncol (x) != r))
  489.       stop("x must be square with at least two rows and columns")
  490.     if (any(x < 0) || any(is.na(x)))
  491.       stop("all entries of x must be nonnegative and finite")
  492.     DNAME <- deparse(substitute(x))
  493.   }
  494.   else {
  495.     if (is.null(y))
  496.       stop("if x is not a matrix, y must be given")
  497.     if (length(x) != length(y))
  498.       stop("x and y must have the same length")
  499.     DNAME <- paste(deparse(substitute(x)), "and", deparse(substitute(y)))
  500.     OK <- complete.cases(x, y)
  501.     x <- as.factor(x[OK])
  502.     y <- as.factor(y[OK])
  503.     r <- nlevels(x)
  504.     if ((r < 2) || (nlevels(y) != r))
  505.       stop("x and y must have the same number of levels (minimum 2)")
  506.     x <- table(x, y)
  507.   }
  508.  
  509.   PARAMETER <- r * (r-1) / 2
  510.   names(PARAMETER) <- "df"
  511.   METHOD <- "McNemar's Chi-square test"
  512.  
  513.   if (correct && (r == 2) && any(x - t(x))) {
  514.     y <- (abs(x - t(x)) - 1)
  515.     METHOD <- paste(METHOD, "with continuity correction")
  516.   }
  517.   else
  518.     y <- x - t(x)
  519.   x <- x + t(x)
  520.  
  521.   STATISTIC <- sum(y[upper.tri(x)]^2 / x[upper.tri(x)])
  522.   names(STATISTIC) <- "McNemar's chi-square"
  523.   PVAL <- 1 - pchisq(STATISTIC, PARAMETER)  
  524.     
  525.   RVAL <- list(statistic = STATISTIC,
  526.            parameter = PARAMETER,
  527.            p.value = PVAL,
  528.            method = METHOD,
  529.            data.name = DNAME)
  530.   class(RVAL) <- "htest"
  531.   return(RVAL)
  532. }
  533. pks <- function(x, tol=10^(-6))
  534. {
  535.   if (is.numeric(x))
  536.     x <- as.vector(x)
  537.   else
  538.     stop("Argument x must be numeric")
  539.  
  540.   PKS <- rep(0, length(x))
  541.   PKS[is.na(x)] <- NA
  542.   IND <- which(!is.na(x) & (x > 0))
  543.   if (length(IND) > 0) {
  544.     k <- 1 : ceiling(sqrt(-log(tol)/2) / min(x[IND]))
  545.     y <- outer(x[IND]^2, k,
  546.            function (t, k) { (-1)^k * exp(-2 * t * k^2) })
  547.     PKS[IND] <- 1 + 2 * apply(y, 1, "sum")
  548.   }
  549.   return(PKS)
  550. }
  551. prop.test <- function(x, n, p = NULL, alternative = "two.sided",
  552.               conf.level = 0.95, correct = TRUE)
  553. {
  554.   DNAME <- paste(deparse(substitute(x)), "out of", deparse(substitute(n)))
  555.  
  556.   if ((l <- length(x)) != length(n))
  557.     stop("x and n must have the same length")
  558.   OK <- complete.cases(x, n)
  559.   x <- x[OK]
  560.   n <- n[OK]
  561.   if ((k <- length(x)) < 1)
  562.     stop("Not enough data")
  563.   if (any(n <= 0))
  564.     stop("Elements of n must be positive")
  565.   if (any(x < 0))
  566.     stop("Elements of x must be nonnegative")
  567.   if (any(x > n))
  568.     stop("Elements of x must not be greater than those of n")
  569.  
  570.   if (is.null(p) && (k == 1))
  571.     p <- .5
  572.   if (!is.null(p)) {
  573.     DNAME <- paste(DNAME, ", null ",
  574.            ifelse(k == 1, "probability ", "probabilities "),
  575.            deparse(substitute(p)), sep = "")
  576.     if (length(p) != l)
  577.       stop("p must have the same length as x and n")
  578.     p <- p[OK]    
  579.     if (any((p <= 0) | (p >= 1)))
  580.       stop("Elements of p must be in (0,1)")
  581.   }
  582.  
  583.   CHOICES <- c("two.sided", "less", "greater")
  584.   alternative <- CHOICES[pmatch(alternative, CHOICES)]
  585.   if (length(alternative) > 1 || is.na(alternative)) 
  586.     stop("alternative must be \"two.sided\", \"less\" or \"greater\"")
  587.   if ((k > 2) || (k == 2) && !is.null(p))
  588.     alternative <- "two.sided"
  589.  
  590.   if ((length(conf.level) != 1) || is.na(conf.level) ||
  591.       (conf.level <= 0) || (conf.level >= 1))
  592.     stop("conf.level must be a single number between 0 and 1")
  593.  
  594.   correct <- as.logical(correct)
  595.   
  596.   ESTIMATE <- x/n
  597.   names(ESTIMATE) <- if (k == 1) "p" else paste("prop", 1:l)[OK]
  598.   NVAL <- p
  599.   CINT <- NULL
  600.   YATES <- ifelse(correct && (k <= 2), .5, 0)
  601.  
  602.   if (k == 1) {
  603.     z <- ifelse(alternative == "two.sided",
  604.         qnorm((1 + conf.level) / 2),
  605.         qnorm(conf.level))
  606.     YATES <- min(YATES, abs(x - n * p))
  607.     p.c <- ESTIMATE + YATES / n
  608.     p.u <- ((p.c + z^2 / (2 * n) 
  609.          + z * sqrt(p.c * (1 - p.c) / n + z^2 / (4 * n^2)))
  610.         / (1 + z^2 / n))
  611.     p.c <- ESTIMATE - YATES / n
  612.     p.l <- ((p.c + z^2 / (2 * n) 
  613.          - z * sqrt(p.c * (1 - p.c) / n + z^2 / (4 * n^2)))
  614.         / (1 + z^2 / n))
  615.     CINT <- switch(alternative,
  616.            "two.sided" = c(max(p.l, 0), min(p.u, 1)),
  617.            "greater" = c(max(p.l, 0), 1),
  618.            "less" = c(0, min(p.u, 1)))
  619.   }
  620.   else if ((k == 2) & is.null(p)) {
  621.     DELTA <- ESTIMATE[1] - ESTIMATE[2]
  622.     YATES <- min(YATES, abs(DELTA) / sum(1/n))
  623.     WIDTH <- (switch(alternative,
  624.              "two.sided" = qnorm((1 + conf.level) / 2),
  625.              qnorm(conf.level))
  626.           * sqrt(sum(ESTIMATE * (1 - ESTIMATE) / n))
  627.           + YATES * sum(1/n))
  628.     CINT <- switch(alternative,
  629.            "two.sided" = c(max(DELTA - WIDTH, -1),
  630.                            min(DELTA + WIDTH, 1)),
  631.            "greater" = c(max(DELTA - WIDTH, -1), 1),
  632.            "less" = c(-1, min(DELTA + WIDTH, 1)))
  633.   }
  634.   if (!is.null(CINT))
  635.     attr(CINT, "conf.level") <- conf.level
  636.  
  637.   METHOD <- paste(ifelse(k == 1,
  638.              "1-sample proportions test",
  639.              paste(k, "-sample test for ",
  640.                    ifelse(is.null(p), "equality of", "given"),
  641.                    " proportions", sep = "")),
  642.           ifelse(YATES, "with", "without"),
  643.           "continuity correction")
  644.  
  645.   if (is.null(p)) {
  646.     p <- sum(x)/sum(n)
  647.     PARAMETER <- k - 1
  648.   }
  649.   else {
  650.     PARAMETER <- k
  651.     names(NVAL) <- names(ESTIMATE)
  652.   }
  653.   names(PARAMETER) <- "df"
  654.  
  655.   x <- cbind(x, n - x)
  656.   E <- cbind(n * p, n * (1 - p))
  657.   if (any(E < 5))
  658.     warning("Chi-square approximation may be incorrect")
  659.   STATISTIC <- sum((abs(x - E) - YATES)^2 / E)
  660.   names(STATISTIC) <- "X-squared"
  661.  
  662.   if (alternative == "two.sided")
  663.     PVAL <- 1 - pchisq(STATISTIC, PARAMETER)
  664.   else {
  665.     if (k == 1)
  666.       z <- sign(ESTIMATE - p) * sqrt(STATISTIC)
  667.     else
  668.       z <- sign(DELTA) * sqrt(STATISTIC)
  669.     if (alternative == "greater")
  670.       PVAL <- 1 - pnorm(z)
  671.     else
  672.       PVAL <- pnorm(z)
  673.   }
  674.  
  675.   RVAL <- list(statistic = STATISTIC,
  676.            parameter = PARAMETER,
  677.            p.value = PVAL,
  678.            estimate = ESTIMATE,
  679.            null.value = NVAL,
  680.            conf.int = CINT,
  681.            alternative = alternative,
  682.            method = METHOD,
  683.            data.name = DNAME)
  684.   class(RVAL) <- "htest"
  685.   return(RVAL)
  686. }
  687. var.test <- function(x, y, ratio = 1, alternative = "two.sided",
  688.              conf.level = 0.95)  
  689. {
  690.   if (!((length(ratio) == 1) && is.finite(ratio) && (ratio > 0)))
  691.     stop("ratio must be a single positive number")
  692.   
  693.   alternative <- char.expand(alternative,
  694.                  c("two.sided", "less", "greater"))
  695.   if ((length(alternative) > 1) || is.na(alternative)) 
  696.     stop("alternative must be \"two.sided\", \"less\" or \"greater\"")
  697.  
  698.   if (!((length(conf.level) == 1) && is.finite(conf.level) &&
  699.     (conf.level > 0) && (conf.level < 1)))
  700.     stop("conf.level must be a single number between 0 and 1")
  701.  
  702.   DNAME <- paste(deparse(substitute(x)), "and", deparse(substitute(y)))
  703.  
  704.   x <- x[is.finite(x)]
  705.   DF.x <- length(x) - 1
  706.   if (DF.x < 1) 
  707.     stop("not enough x observations")
  708.   y <- y[is.finite(y)]
  709.   DF.y <- length(y) - 1
  710.   if (DF.y < 1) 
  711.     stop("not enough y observations")
  712.  
  713.   V.x <- var(x)
  714.   V.y <- var(y)
  715.   ESTIMATE <- V.x / V.y
  716.   STATISTIC <- ESTIMATE / ratio
  717.   PARAMETER <- c(DF.x, DF.y)
  718.  
  719.   PVAL <- pf(STATISTIC, DF.x, DF.y)
  720.   if (alternative == "two.sided") {
  721.     PVAL <- 2 * min(PVAL, 1 - PVAL)
  722.     BETA <- (1 - conf.level) / 2
  723.     CINT <- c(ESTIMATE / qf(1 - BETA, DF.x, DF.y),
  724.           ESTIMATE / qf(BETA, DF.x, DF.y))
  725.   }
  726.   else if (alternative == "greater") {
  727.     PVAL <- 1 - PVAL
  728.     CINT <- c(ESTIMATE / qf(conf.level, DF.x, DF.y), NA)
  729.   }
  730.   else
  731.     CINT <- c(0, ESTIMATE / qf(1 - conf.level, DF.x, DF.y))
  732.   names(STATISTIC) <- "F"
  733.   names(PARAMETER) <- c("num df", "denom df")
  734.   names(ESTIMATE) <- names(ratio) <- "ratio of variances"
  735.   attr(CINT, "conf.level") <- conf.level
  736.   RVAL <- list(statistic = STATISTIC,
  737.            parameter = PARAMETER,
  738.            p.value = PVAL,
  739.            conf.int = CINT,
  740.            estimate = ESTIMATE,
  741.            null.value = ratio,
  742.            alternative = alternative,
  743.            method = "F test to compare two variances",
  744.            data.name = DNAME)
  745.   attr(RVAL, "class") <- "htest"
  746.   return(RVAL)
  747. }
  748. wilcox.test <- function(x, y = NULL, alternative = "two.sided", mu = 0,
  749.             paired = FALSE, exact = FALSE, correct = TRUE) 
  750. {
  751.   CHOICES <- c("two.sided", "less", "greater")
  752.   alternative <- CHOICES[pmatch(alternative, CHOICES)]
  753.   if (length(alternative) > 1 || is.na(alternative)) 
  754.     stop("alternative must be \"two.sided\", \"less\" or \"greater\"")
  755.  
  756.   if (!missing(mu) && ((length(mu) > 1) || !is.finite(mu)))
  757.     stop("mu must be a single number")
  758.  
  759.   if (!is.null(y)) {
  760.     DNAME <- paste(deparse(substitute(x)), "and", deparse(substitute(y)))
  761.     if (paired) {
  762.       if (length(x) != length(y))
  763.     stop("x and y must have the same length")
  764.       OK <- complete.cases(x, y)
  765.       x <- x[OK] - y[OK]
  766.       y <- NULL
  767.     }
  768.     else {
  769.       x <- x[is.finite(x)]
  770.       y <- y[is.finite(y)]
  771.     }
  772.   }
  773.   else {
  774.     DNAME <- deparse(substitute(x))
  775.     if (paired)
  776.       stop("y missing for paired test")
  777.     x <- x[is.finite(x)]
  778.   }
  779.  
  780.   if (length(x) < 1)
  781.     stop("not enough x observations")
  782.  
  783.   if (exact)
  784.     warning("exact Wilcoxon tests are currently not implemented\n")
  785.  
  786.   PARAMETER <- NULL
  787.   
  788.   CORRECTION <- 0  
  789.  
  790.   if (is.null(y)) {
  791.     METHOD <- "Wilcoxon signed rank test"
  792.     x <- x - mu
  793.     r <- rank(abs(x))
  794.     STATISTIC <- sum(r[x > 0])
  795.     names(STATISTIC) <- "V"
  796.     n <- length(x)
  797.     TIES <- table(r)
  798.     z <- STATISTIC - n*(n+1)/4
  799.     SIGMA <- sqrt(n*(n+1)*(2*n+1)/24 - sum(TIES^3-TIES)/48)
  800.   }
  801.   else {
  802.     if (length(y) < 1)
  803.       stop("not enough y observations")
  804.     METHOD <- "Wilcoxon rank sum test"
  805.     r <- rank(c(x - mu, y))
  806.     STATISTIC <- sum(r[seq(along = x)])
  807.     names(STATISTIC) <- "W"
  808.     n.x <- length(x)
  809.     n.y <- length(y)
  810.     TIES <- table(r)
  811.     z <- STATISTIC - n.x*(n.x+n.y+1)/2
  812.     SIGMA <- sqrt((n.x*n.y/12) *
  813.           ((n.x+n.y+1)
  814.            - sum(TIES^3-TIES) / ((n.x+n.y)*(n.x+n.y+1))))
  815.   }
  816.  
  817.   if (correct) {
  818.     CORRECTION <- switch ("two.sided", sign(z) * 0.5,
  819.               "greater", 0.5,
  820.               "less", -0.5)
  821.     METHOD <- paste(METHOD, "with continuity correction")
  822.   }
  823.   
  824.   PVAL <- pnorm((z - CORRECTION) / SIGMA)
  825.   if (alternative == "two.sided")
  826.     PVAL <- 2 * min (PVAL, 1 - PVAL)
  827.   else if (alternative == "greater")
  828.     PVAL <- 1 - PVAL
  829.  
  830.   NVAL <- mu
  831.   names(NVAL) <- "mu"
  832.  
  833.   RVAL <- list(statistic = STATISTIC,
  834.            parameter = PARAMETER,
  835.            p.value = PVAL,
  836.            null.value = NVAL,
  837.            alternative = alternative,
  838.            method = METHOD,
  839.            data.name = DNAME)
  840.   class(RVAL) <- "htest"
  841.   return(RVAL)
  842. }
  843.