-
Notifications
You must be signed in to change notification settings - Fork 0
/
pvalRW.m
80 lines (73 loc) · 2.26 KB
/
pvalRW.m
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
% Calculating the efficient p-values based on Bootstap based on the paper
% of Romano and Wolf (2016) paper link:
% http://www.sciencedirect.com/science/article/pii/S0167715216000389
% Script developed by Arman Hassanniakalager for the 3rd chapter of PhD,
% Created on 04 Jul 2017 17:48 BST,
% Last modified 07 Jul 2017 14:59 BST.
% ****************************************
%% Input arguments:
% -tststat: test statistic of the hyposethes
% -tststatB: Bootstrap test statistic of the hyposethes
% -mondj: adjustment for monotonicity
function pvaladj=pvalRW(tststat,tststatBUN,monadj)
if nargin<3
monadj=1;
end
tststatBCNT=tststatBUN-tststat;
[toptststat,rperm]=sort(tststat,'descend');
topboots=tststatBCNT(:,rperm);
M=size(tststatBCNT,1);
S=size(tststatBCNT,2);
for m=1:M
for j=1:S
max_t(m,j)=max(topboots(m,j:S));
end
end
pvalrperm=nan(1,numel(tststat));
excesstimes1=sum(max_t(:,1)>=toptststat(1));
excesstimes2=sum(max_t(:,1)<=toptststat(1));
excesstimes=min(excesstimes1,excesstimes2);
pvalrperm(1)=(2*excesstimes+1)/(M+1);
for s=2:S
excesstimes1=sum(max_t(:,s)>=toptststat(s));
excesstimes2=sum(max_t(:,s)<=toptststat(s));
excesstimes=min(excesstimes1,excesstimes2);
pvalrperm(s)=(2*excesstimes+1)/(M+1);
pvalrperm(s)=monadj*max(pvalrperm(s),pvalrperm(s-1))...
+(1-monadj)*pvalrperm(s);
end
pvaladj(rperm)=pvalrperm;
pvaladj(isnan(tststat)) = nan;
end
% function pvaladj=pvalRW(tststat,tststatBUN,monadj)
% if nargin<3
% monadj=1;
% end
% tststatB=tststatBUN-tststat;
% [toptststat,rperm]=sort(tststat,'descend');
% topboots=tststatB(:,rperm);
% M=size(tststatB,1);
% S=size(tststatB,2);
% for m=1:M
% for j=1:S
% max_t(m,j)=max(topboots(m,j:S));
% end
% end
% pvalrperm=nan(1,numel(tststat));
% excesstimes=sum(max_t(:,1)>=toptststat(1));
% if excesstimes==0
% excesstimes=nan;
% end
% pvalrperm(1)=(excesstimes+1)/(M+1);
% for s=2:S
% excesstimes=sum(max_t(:,s)>=toptststat(s));
% if excesstimes==0
% excesstimes=nan;
% end
% pvalrperm(s)=(excesstimes+1)/(M+1);
% pvalrperm(s)=monadj*max(pvalrperm(s),pvalrperm(s-1))...
% +(1-monadj)*pvalrperm(s);
% end
% pvaladj(rperm)=pvalrperm;
%
% end