-
Notifications
You must be signed in to change notification settings - Fork 23
/
wilcox_modified.txt
142 lines (138 loc) · 4.58 KB
/
wilcox_modified.txt
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
## ===================================================
## The functions in this file were modified from the
## original code from Rand Wilcox:
## http://dornsife.usc.edu/labs/rwilcox/software/
## ===================================================
shifthd<-function(x,y,nboot=200){
#
# Compute confidence intervals for the difference between deciles
# of two independent groups. The simultaneous probability coverage is .95.
# The Harrell-Davis estimate of the qth quantile is used.
# The default number of bootstrap samples is nboot=200
#
# The results are stored and returned in a 9 by 5 matrix,
# the ith row corresponding to the i/10 quantile.
# The first column is the lower end of the confidence interval.
# The second column is the upper end.
# The third column is the estimated difference between the deciles
# (second group minus first).
#
# Modified from Rallfun-v31: GAR - University of Glasgow - 2016-06-22
# simplified inputs/outputs to make independent figures
# changed difference to x-y instead of y-x for consistency across functions
# now outputs confidence intervals, similarly to shiftdhd
# changed labels to reflect nature of data: 'groups' here, 'conditions' for shiftdhd
# now returns a data frame for use with ggplot2
x<-x[!is.na(x)]
y<-y[!is.na(y)]
crit<-80.1/(min(length(x),length(y)))^2+2.73
m<-matrix(0,9,5)
for (d in 1:9){
q<-d/10
data<-matrix(sample(x,size=length(x)*nboot,replace=TRUE),nrow=nboot)
bvec<-apply(data,1,hd,q)
sex<-var(bvec)
data<-matrix(sample(y,size=length(y)*nboot,replace=TRUE),nrow=nboot)
bvec<-apply(data,1,hd,q)
sey<-var(bvec)
m[d,1]<-hd(x,q)
m[d,2]<-hd(y,q)
m[d,3]<-m[d,1]-m[d,2]
m[d,4]<-m[d,3]-crit*sqrt(sex+sey)
m[d,5]<-m[d,3]+crit*sqrt(sex+sey)
}
out<-data.frame(m)
names(out)<-c('group1','group2','difference','ci_lower','ci_upper')
out
}
shiftdhd<-function(x,y,nboot=200){
#
# Compute confidence intervals for the difference between deciles
# of two dependent groups. The simultaneous probability coverage is .95.
# The Harrell-Davis estimate of the qth quantile is used.
# The default number of bootstrap samples is nboot=200
#
# The results are stored and returned in a 9 by 5 matrix,
# the ith row corresponding to the i/10 quantile.
# The first column is the lower end of the confidence interval.
# The second column is the upper end.
# The third column is the estimated difference between the deciles
# (second group minus first).
# The fourth column contains the estimated standard error.
#
# No missing values are allowed.
#
# If the goal is to use an alpha value different from .05,
# use the function qcomdhd or qdec2ci
#
# Modified from Rallfun-v31: GAR - University of Glasgow - 2016-06-22
# simplified inputs/outputs for use with new shift_function_figure()
# changed labels to reflect nature of data: 'conditions' here, 'groups' for shifthd
# now returns a data frame for use with ggplot2
xy=elimna(cbind(x,y))
x=xy[,1]
y=xy[,2]
crit<-37/length(x)^(1.4)+2.75
m<-matrix(0,9,5)
data<-matrix(sample(length(x),size=length(x)*nboot,replace=TRUE),nrow=nboot)
xmat<-matrix(x[data],nrow=nboot,ncol=length(x))
ymat<-matrix(y[data],nrow=nboot,ncol=length(x))
for (d in 1:9){
q<-d/10
bvec<-apply(xmat,1,hd,q)-apply(ymat,1,hd,q)
se<-sqrt(var(bvec))
m[d,1]=hd(x,q)
m[d,2]=hd(y,q)
m[d,3]<-m[d,1]-m[d,2]
m[d,4]<-m[d,3]-crit*se
m[d,5]<-m[d,3]+crit*se
}
out<-data.frame(m)
names(out)<-c('condition1','condition2','difference','ci_lower','ci_upper')
out
}
q1469<-function(x){
#
# Estimate the deciles for the data in vector x
# using the Harrell-Davis estimate of the qth quantile
# Modified from `deciles` to return only quantiles 1:4 & 6:9
# GAR, University of Glasgow, 2016-07-14
xs<-sort(x)
n<-length(x)
vecx<-seq(along=x)
xq<-0
todo<-c(1:4,6:9)
for (qi in 1:length(todo)){
q<-todo[qi]/10
m1<-(n+1)*q
m2<-(n+1)*(1-q)
wx<-pbeta(vecx/n,m1,m2)-pbeta((vecx-1)/n,m1,m2) # W sub i values
xq[qi]<-sum(wx*xs)
}
xq
}
quantiles_pbci<-function(x,q=seq(1,9)/10,nboot=2000,alpha=0.05){
# Merges `deciles` & `qcipb` from `Rallfun-v31.txt`
# Compute percentile bootstrap confidence intervals of quantiles, `quant`,
# estimated using the Harrell-Davis estimator.
# Default to deciles
# GAR, University of Glasgow, 2016-07-15
low<-round((alpha/2)*nboot)
up<-nboot-low
low<-low+1
x=elimna(x)
nq=length(q)
output=matrix(NA,ncol=4,nrow=nq)
dimnames(output)=list(NULL,c("quantile","est_q","ci.low","ci.up"))
bdata<-matrix(sample(x,size=length(x)*nboot,replace=TRUE),nrow=nboot) # bootstrap samples
for (qi in 1:nq){
output[qi,1]=q[qi]
output[qi,2]=hd(x,q=q[qi])
bvec<-apply(bdata,1,hd,q=q[qi])
bvec<-sort(bvec)
output[qi,3]=bvec[low]
output[qi,4]=bvec[up]
}
output <- data.frame(output)
print(output)
}