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