-
Notifications
You must be signed in to change notification settings - Fork 2
/
wmwTest.ado
250 lines (218 loc) · 5.48 KB
/
wmwTest.ado
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
*! version 1.0.0 26jan2021
cap program drop wmwTest
program wmwTest, rclass
version 15
syntax varname [if] [in], by(varname) ///
[phiNull(string) CORRect alternative(string) Level(cilevel) NOTIEadjust vmethod(string) ]
tempname glev
tempname touse
marksample `touse', novarlist
*ignore missings in by variable
cap assert !missing(`by')
if _rc {
//dis "Missings in by variable will be ignored"
cap replace ``touse''=0 if missing(`by')
}
*ignore missings in varname
cap assert !missing(`varlist')
if _rc {
//dis as text "Missings in `varlist' will be ignored"
cap replace ``touse''=0 if missing(`varlist')
}
*load moremata for roots
cap mata mata which mm_root()
if _rc {
di as error "mm_root() from -moremata- is required; type -ssc install moremata- to obtain it"
exit 499
}
*null hyopthesis for test, default: 0.5
if "`phiNull'"=="" {
local phi0 0.5
}
else {
local phi0 `phiNull'
}
*continuity correction if correct is included
if "`correct'"=="" {
local correct 0
}
else {
local correct 1
}
*variance function for the asymptotic method, LAPH or PO, very similar he first is faster
if "`vmethod'"=="" | "`vmethod'" == "LAPH" {
local vmeth 0
}
if "`vmethod'"=="PO" {
local vmeth 1
}
if "`vmethod'"!="" & !inlist("`vmethod'","LAPH","PO") {
dis as text "Invalid vmethod, default (LAPH) used"
local vmeth 0
}
assert inlist(`vmeth',0,1)
*alternative
if "`alternative'"=="" | "`alternative'"=="two-sided" {
local alt 0
}
if "`alternative'"=="greater" {
local alt 1
}
if "`alternative'"=="less" {
local alt 2
}
if "`alternative'"!="" & !inlist("`alternative'","two-sided","greater","less") {
dis as text "Invalid alternative, only two-sided, greater or less allowed. Default two-sided used."
local alt 0
}
*level
if "`level'"=="" {
local level 95
}
local alpha = (100 - `level')/100
qui levelsof `by' if ``touse'', local(`glev')
local wc: word count ``glev''
assert `wc'==2
forvalues i=1/`wc' {
local l`i': word `i' of ``glev''
qui count if `by'==`l`i'' & ``touse''
local n`i'=`r(N)'
}
qui count if ``touse''
local n=`r(N)'
assert `n'==`n1' + `n2'
tempvar r
qui egen `r'=rank(`varlist') if ``touse''
qui sum `r' if `by'==`l2' & ``touse''
local sum2 = `r(sum)'
local phi = (`sum2' - `n2' * (`n2' + 1)/2)/(`n2' * `n1')
preserve
qui keep if ``touse''
tempvar nties ntiesd
qui duplicates tag `r', gen(`nties')
qui replace `nties'=`nties'+1
qui duplicates drop `r', force
qui gen `ntiesd'=`nties'^3-`nties'
qui sum `ntiesd'
local ntiess=`r(sum)'
restore
local tiefactor=1-`ntiess'/(`n'*(`n'+1)*(`n'-1))
if "`notieadjust'" != "" {
if (`tiefactor' < 1) {
dis as error "tie factor removed, but there are ties in the data."
}
local tiefactor = 1
}
if (`tiefactor' == 0) {
local lci = 0
local uci = 1
local pv = 1
}
else {
local corrless=0
local corrgr=0
if (`correct'==1) {
local corrless = -0.5/(`n2'* `n1')
local corrgr = 0.5/(`n2' * `n1')
}
if (`vmeth' == 0) {
local v = `tiefactor' * (`phi0' * (1 - `phi0')/(`n2' * `n1')) * ///
(1 + ((`n2' + `n1' - 2)/2) * ((1 - `phi0')/(2 - `phi0') + `phi0'/(1 + `phi0')))
}
else {
//Vpo(PHI, tf, ny, nx)
//?Vpo
//more complicated not implemented at the moment
dis as error "vmethod PO not implemented at the moment"
exit, clear
}
local z0Less = (`phi' - `phi0' - `corrless')/sqrt(`v')
local z0Gr = (`phi' - `phi0' - `corrgr')/sqrt(`v')
if `alt'==0 {
local pv=2*min(normal(`z0Less'),1-normal(`z0Gr'))
}
if `alt'==1 {
local pv=1-normal(`z0Gr')
}
if `alt'==2 {
local pv=normal(`z0Less')
}
local epsilon = 10^(-8)
local phimin = `epsilon'
local phimax = 1 - `epsilon'
if `alt'==0 {
local zql=invnormal(`alpha'/2)
local zqu=invnormal(1-`alpha'/2)
}
if inlist(`alt',1,2) {
local zql=invnormal(`alpha')
local zqu=invnormal(1-`alpha')
}
tempname A CI
matrix `A'=`phi',`tiefactor',`corrgr',`corrless',`n1',`n2',`zql',`zqu',`phimin',`phimax',`epsilon'
qui mata: mci("`A'","`CI'")
if `alt'==0 {
local lci = `CI'[1,1]
local uci = `CI'[1,2]
local lcis=string(`lci',"%10.4f")
local ucis=string(`uci',"%10.4f")
}
if `alt'==1 {
local lci = `CI'[1,1]
local uci = 1
local lcis=string(`lci',"%10.4f")
local ucis=string(`uci',"%10.4f")
}
if `alt'==2 {
local lci = 0
local uci = `CI'[1,2]
local lcis=string(`lci',"%10.4f")
local ucis=string(`uci',"%10.4f")
}
}
local phis=string(`phi',"%10.4f")
//dis as text "Mann-Whitney statistic:"
//dis as text "`phis' (`lcis' to `ucis')"
//dis as text "P-value = `=string(`pv'),"%4.4f")'"
return scalar p=`pv'
return scalar uci=`uci'
return scalar lci=`lci'
return scalar estimate=`phi'
end
capture mata mata drop mci
clear mata
mata:
void mci(string scalar vname, string scalar CI)
{
A = st_matrix(vname)
PHI = A[1,1]
TF = A[1,2]
CORRGR = A[1,3]
CORRLESS = A[1,4]
nx = A[1,5]
ny = A[1,6]
ZQL = A[1,7]
ZQU = A[1,8]
PHIMIN = A[1,9]
PHIMAX = A[1,10]
EPSILON=A[1,11]
mm_root(x=., &myfunc_lci(), PHIMIN, PHIMAX, EPSILON, 1000, PHI, CORRGR, TF, nx, ny, ZQU)
lci=x
mm_root(x=., &myfunc_uci(), PHIMIN, PHIMAX, EPSILON, 1000, PHI, CORRLESS, TF, nx, ny, ZQL)
uci=x
ma=lci, uci
st_matrix(CI, ma)
}
end
mata:
function myfunc_lci(x, PHI, CORRGR, TF, nx, ny, ZQU) ///
return((PHI - x - CORRGR)/ ///
sqrt(TF * (x * (1 - x)/(ny * nx)) * (1 + ((ny + nx - 2)/2) * ((1 - x)/(2 - x) + x/(1 + x)))) ///
- ZQU)
end
mata:
function myfunc_uci(x, PHI, CORRLESS, TF, nx, ny, ZQL) ///
return((PHI - x - CORRLESS)/ ///
sqrt(TF * (x * (1 - x)/(ny * nx)) * (1 + ((ny + nx - 2)/2) * ((1 - x)/(2 - x) + x/(1 + x)))) ///
- ZQL)
end