# This library is free software; you can redistribute it and/or # modify it under the terms of the GNU Library General Public # License as published by the Free Software Foundation; either # version 2 of the License, or (at your option) any later version. # # This library is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Library General Public License for more details. # # You should have received a copy of the GNU Library General # Public License along with this library; if not, write to the # Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, # MA 02111-1307 USA # Copyrights (C) # for this R-port: # 1999 - 2008, Diethelm Wuertz, Rmetrics Foundation, GPL # Diethelm Wuertz # www.rmetrics.org # for the code accessed (or partly included) from other R-ports: # see R's copyright and license files # for the code accessed (or partly included) from contributed R-ports # and other sources # see Rmetrics's copyright file ################################################################################ # FUNCTION: LOCATION TESTS: # locationTest Performs locations tests on two samples # .tTest Unpaired t test for differences in mean # .kw2Test Kruskal-Wallis test for differences in locations ################################################################################ locationTest <- function(x, y, method = c("t", "kw2"), title = NULL, description = NULL) { # A function implemented by Diethelm Wuertz # Description: # Correlation Tests # FUNCTION: # Test: method = match.arg(method) if (method == "t") { ans = .tTest(x, y, title = title, description = description) } if (method == "kw2") { ans = .kw2Test(x, y, title = title, description = description) } # Return Value: ans } # ------------------------------------------------------------------------------ .tTest <- function(x, y, title = NULL, description = NULL) { # A function implemented by Diethelm Wuertz # Description: # Tests if two population means are equal. # Arguments: # x, y - two numeric vector of data values or time series objects # description - a brief description of the porject of type character. # title - a character string which allows for a project title. # FUNCTION: # Call: call = match.call() # Test: test = list() # Data Set Name: DNAME = paste(deparse(substitute(x)), "and", deparse(substitute(y))) test$data.name = DNAME # Convert Type: x = as.vector(x) y = as.vector(y) # Asymptotic Test: two.sided = t.test(x = x, y = y, alternative = "two.sided", mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) less = t.test(x = x, y = y, alternative = "less", mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) greater = t.test(x = x, y = y, alternative = "greater", mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) # Assume Equal Variances: two.sided.equal = t.test(x = x, y = y, alternative = "two.sided", mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95) less.equal = t.test(x = x, y = y, alternative = "less", mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95) greater.equal = t.test(x = x, y = y, alternative = "greater", mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95) # Sample Estimates: PARAMETER = c(length(x), length(y), 0) names(PARAMETER) = c( "x Observations", "y Observations", "mu") test$parameter = PARAMETER # Sample Estimates: ESTIMATE = c(two.sided$estimate, var(x), var(y)) names(ESTIMATE) = c("Mean of x", "Mean of y", "Var of x", "Var of y") test$estimate = ESTIMATE # P Values: PVAL = c( two.sided$p.value, less$p.value, greater$p.value, two.sided.equal$p.value, less.equal$p.value, greater.equal$p.value) names(PVAL) = c( "Alternative Two-Sided", "Alternative Less", "Alternative Greater", "Alternative Two-Sided | Equal Var", "Alternative Less | Equal Var", "Alternative Greater | Equal Var") test$p.value = PVAL # Statistic: STATISTIC = c( two.sided$statistic, two.sided.equal$statistic) names(STATISTIC) = c( " T", "T | Equal Var") test$statistic = STATISTIC # Confidence Intervals: CONF.INT = cbind( a = two.sided$conf.int, b = less$conf.int, c = greater$conf.int, d = two.sided.equal$conf.int, e = less.equal$conf.int, f = greater.equal$conf.int) # For Splus compatibility use named a CONF.INT # and dimnames instead of colnames! dimnames(CONF.INT)[[2]] = c( "Two-Sided", " Less", " Greater", "Two-Sided | Equal Var", " Less | Equal Var", " Greater | Equal Var") test$conf.int = CONF.INT # Add: if (is.null(title)) title = "t Test" if (is.null(description)) description = date() # Return Value: new("fHTEST", call = call, data = list(x = x, y = y), test = test, title = as.character(title), description = as.character(description) ) } # ------------------------------------------------------------------------------ .kw2Test <- function(x, y, title = NULL, description = NULL) { # A function implemented by Diethelm Wuertz # Description: # Performs a Kruskal-Wallis rank sum test of the null that # the location parameters of the distribution of x are the # same in each group (sample). The alternative is that they # differ in at least one. # Arguments: # x, y - two numeric vector of data values or time series objects # description - a brief description of the porject of type character. # title - a character string which allows for a project title. # Note: # A function linked to "stats" # FUNCTION: # Call: call = match.call() # Test: test = list() # Data Set Name: DNAME = paste(deparse(substitute(x)), "and", deparse(substitute(y))) test$data.name = DNAME # Convert Type: x = as.vector(x) y = as.vector(y) # Sample Estimates: ESTIMATE = c(mean(x), mean(y), var(x), var(y)) names(ESTIMATE) = c("Mean of x", "Mean of y", "Var of x", "Var of y") test$estimate = ESTIMATE # Parameter: PARAMETER = c(length(x), length(y)) names(PARAMETER) = c( "x Observations", "y Observations") test$parameter = PARAMETER # Operate on Lists: x = list(x = x, y = y) if (length(x) < 2) stop("x must be a list with at least 2 elements") k = length(x) l = sapply(x, "length") g = factor(rep(1 : k, l)) x = unlist(x) # Test: n = length(x) if (n < 2) stop("not enough observations") r = rank(x) TIES = table(x) # Statistic: STATISTIC = sum(tapply(r, g, "sum")^2 / tapply(r, g, "length")) STATISTIC = ((12 * STATISTIC / (n * (n + 1)) - 3 * (n + 1)) / (1 - sum(TIES^3 - TIES) / (n^3 - n))) names(STATISTIC) = "KW chi-squared" test$statistic = STATISTIC # P Value: PVAL = 1 - pchisq(STATISTIC, 1) names(PVAL) = "" test$p.value = PVAL # Add: if(is.null(title)) title = "Kruskal-Wallis Two Sample Test" if(is.null(description)) description = date() # Return Value: new("fHTEST", call = call, data = list(x = x, y = y), test = test, title = as.character(title), description = as.character(description) ) } ################################################################################