1 2# This library is free software; you can redistribute it and/or 3# modify it under the terms of the GNU Library General Public 4# License as published by the Free Software Foundation; either 5# version 2 of the License, or (at your option) any later version. 6# 7# This library is distributed in the hope that it will be useful, 8# but WITHOUT ANY WARRANTY; without even the implied warranty of 9# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 10# GNU Library General Public License for more details. 11# 12# You should have received a copy of the GNU Library General 13# Public License along with this library; if not, write to the 14# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 15# MA 02111-1307 USA 16 17 18################################################################################ 19# FUNCTION: DESCRIPTION: 20# scaleTest Performs scale tests on two samples 21# .ansariTest Ansari-Bradley test for differences in scale 22# .moodTest Mood test for differences in scale 23################################################################################ 24 25 26scaleTest <- 27function(x, y, method = c("ansari", "mood"), 28 title = NULL, description = NULL) 29{ 30 # A function implemented by Diethelm Wuertz 31 32 # Description: 33 # Performs Scale Tests 34 35 # FUNCTION: 36 37 # Test: 38 method = match.arg(method) 39 if (method == "ansari") { 40 ans = .ansariTest(x, y, title = title, description = description) 41 } 42 if (method == "mood") { 43 ans = .moodTest(x, y, title = title, description = description) 44 } 45 46 # Return Value: 47 ans 48} 49 50 51# ------------------------------------------------------------------------------ 52 53 54.ansariTest <- 55function(x, y, title = NULL, description = NULL) 56{ 57 # A function implemented by Diethelm Wuertz 58 59 # Description: 60 # Ansari-Bradley's test for differences in scale 61 62 # Arguments: 63 # x - a numeric vector of data values. 64 # description - a brief description of the porject of type 65 # character. 66 # title - a character string which allows for a project title. 67 68 # FUNCTION: 69 70 # Call: 71 call = match.call() 72 73 # Test: 74 test = list() 75 76 # Data Set Name: 77 DNAME = paste(deparse(substitute(x)), "and", deparse(substitute(y))) 78 test$data.name = DNAME 79 80 # Convert Type: 81 x = as.vector(x) 82 y = as.vector(y) 83 84 # Test: 85 two.sided = .ansari2Test(x = x, y = y, alternative = "two.sided", 86 exact = FALSE, conf.int = TRUE, conf.level = 0.95) 87 less = .ansari2Test(x = x, y = y, alternative = "less", 88 exact = FALSE, conf.int = TRUE, conf.level = 0.95) 89 greater = .ansari2Test(x = x, y = y, alternative = "greater", 90 exact = FALSE, conf.int = TRUE, conf.level = 0.95) 91 92 two.sided.exact = .ansari2Test(x = x, y = y, alternative = "two.sided", 93 exact = TRUE, conf.int = TRUE, conf.level = 0.95) 94 less.exact = .ansari2Test(x = x, y = y, alternative = "less", 95 exact = TRUE, conf.int = TRUE, conf.level = 0.95) 96 greater.exact = .ansari2Test(x = x, y = y, alternative = "greater", 97 exact = TRUE, conf.int = TRUE, conf.level = 0.95) 98 99 # Statistic: 100 STATISTIC = c(two.sided$statistic) 101 names(STATISTIC) = "AB" 102 test$statistic = STATISTIC 103 104 # P Values: 105 PVAL = c( 106 two.sided$p.value, 107 two.sided.exact$p.value, 108 less$p.value, 109 less.exact$p.value, 110 greater$p.value, 111 greater.exact$p.value) 112 names(PVAL) = c( 113 "Alternative Two-Sided ", 114 "Alternative Two-Sided | Exact", 115 "Alternative Less ", 116 "Alternative Less | Exact", 117 "Alternative Greater ", 118 "Alternative Greater | Exact") 119 test$p.values = PVAL 120 121 # Confidence Levels: 122 CONF.INT = cbind( 123 a = two.sided$conf.int, 124 b = two.sided.exact$conf.int, 125 c = less$conf.int, 126 d = less.exact$conf.int, 127 e = greater$conf.int, 128 f = greater.exact$conf.int) 129 # For Splus compatibility use named a CONF.INT 130 # and dimnames instead of colnames! 131 dimnames(CONF.INT)[[2]] = c( 132 "Two-Sided | Asymptotic ", 133 "Two-Sided | Exact ", 134 "Less | Asymptotic ", 135 "Less | Exact ", 136 "Greater | Asymptotic ", 137 "Greater | Exact ") 138 test$conf.int = CONF.INT 139 140 # Add: 141 if(is.null(title)) title = "Ansari-Bradley Test for Scale" 142 if(is.null(description)) description = date() 143 144 # Return Value: 145 new("fHTEST", 146 call = call, 147 data = list(x = x, y = y), 148 test = test, 149 title = as.character(title), 150 description = as.character(description) ) 151} 152 153 154# ------------------------------------------------------------------------------ 155 156 157.ansari2Test <- 158function(x, y, alternative = c("two.sided", "less", "greater"), 159 exact = TRUE, conf.int = FALSE, conf.level = 0.95, ...) 160{ 161 # A function implemented by Diethelm Wuertz 162 163 # Arguments: 164 # x - numeric vector of data values. 165 # y - numeric vector of data values. 166 # alternative - indicates the alternative hypothesis and must 167 # be one of "two.sided", "greater" or "less". You can specify 168 # just the initial letter. 169 # exact - a logical indicating whether an exact p-value should 170 # be computed. 171 # conf.int - a logical,indicating whether a confidence interval 172 # should be computed. 173 # conf.level - confidence level of the interval. 174 175 # FUNCTION: 176 177 # Return Value: 178 ansari.test(x, y, alternative, exact, conf.int, conf.level, ...) 179} 180 181 182# ------------------------------------------------------------------------------ 183 184 185.moodTest <- 186function(x, y, title = NULL, description = NULL) 187{ 188 # A function implemented by Diethelm Wuertz 189 190 # Description: 191 # Performs Mood's two-sample test for a difference in 192 # scale parameters. 193 194 # Arguments: 195 # x, y - a numeric vector of data values. 196 # description - a brief description of the porject of type 197 # character. 198 # title - a character string which allows for a project title. 199 200 # Notes: 201 # A modified copy originally from R's ctest package Version 1.8.1 202 203 # FUNCTION: 204 205 # Call: 206 call = match.call() 207 208 # Test: 209 test = list() 210 211 # Data Set Name: 212 DNAME = paste(deparse(substitute(x)), "and", deparse(substitute(y))) 213 test$data.name = DNAME 214 215 # Convert Type: 216 x = as.vector(x) 217 y = as.vector(y) 218 219 # Check Data: 220 m = length(x) 221 n = length(y) 222 s = m + n 223 if (s < 3) stop("not enough observations") 224 225 # Statistic: 226 r = rank(c(x, y)) 227 # From R: 228 # z = ((sum((r[1:length(x)] - (s + 1) / 2)^2) - m * (s^2 - 1) / 12) 229 # / sqrt(m * n * (s + 1) * (s + 2) * (s - 2) / 180) ) 230 # To run also under S-Plus use ... 231 a = sum( (r[1:length(x)]-0.5*(s+1))^2 ) - m*(s*s-1)/12 232 b = sqrt(as.numeric(m) * n * (s + 1) * (s + 2) * (s - 2) / 180) 233 STATISTIC = a/b 234 names(STATISTIC) = "Z" 235 test$statistic = STATISTIC 236 237 # P Values: 238 p = pnorm(STATISTIC) 239 PVAL = c(2 * min(p, 1 - p), p, 1 - p) 240 names(PVAL) = c( 241 "Alternative Two-Sided", 242 "Alternative Less", 243 "Alternative Greater") 244 test$p.value = PVAL 245 246 # Add: 247 if(is.null(title)) title = "Mood Two-Sample Test of Scale" 248 if(is.null(description)) description = date() 249 250 # Return Value: 251 new("fHTEST", 252 call = call, 253 data = list(x = x, y = y), 254 test = test, 255 title = as.character(title), 256 description = as.character(description) ) 257} 258 259 260################################################################################ 261 262