How to get correct p-values for the Steel-Dwass test?
10 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
I would like to create a code file for the Steel-Dwass test. I use a following smaple data with R (from: https://www.trifields.jp/introducing-steel-dwass-in-r-1632).
----------------------------------------------------------------------------------------------------------------------
> data <- c(
6.9, 7.5, 8.5, 8.4, 8.1, 8.7, 8.9, 8.2, 7.8, 7.3, 6.8,
9.6, 9.4, 9.5, 8.5, 9.4, 9.9, 8.7, 8.1, 7.8, 8.8,
5.7, 6.4, 6.8, 7.8, 7.6, 7.0, 7.7, 7.5, 6.8, 5.9,
7.6, 8.7, 8.5, 8.5, 9.0, 9.2, 9.3, 8.0, 7.2, 7.9, 7.8
)
> group <- rep(1:4, c(11, 10, 10, 11))
> Steel.Dwass(data, group)
t p
1:2 2.680234 0.036960431
1:3 2.539997 0.053980573
1:4 1.282642 0.574011771
2:3 3.746076 0.001031145
2:4 2.046776 0.170965537
3:4 3.384456 0.003976894
----------------------------------------------------------------------------------------------------------------------
I got the same t-values as bottom, but I could not get the same p-value. Could you help me?
Here is my code:
----------------------------------------------------------------------------------------------------------------------
data = [6.9 9.6 5.7 7.6; 7.5, 9.4, 6.4, 8.7; 8.5, 9.5, 6.8, 8.5; 8.4, 8.5, 7.8, 8.5; 8.1, 9.4, 7.6, 9; 8.7, 9.9, 7, 9.2; 8.9, 8.7, 7.7, 9.3;...
8.2, 8.1, 7.5, 8; 7.8, 7.8, 6.8, 7.2; 7.3, 8.8, 5.9, 7.9; 6.8, NaN, NaN, 7.8];
groupNum = size(data, 2);
N = zeros(groupNum, 1);
for i = 1:groupNum
N(i, 1) = length(data(:, i)) - sum(isnan(data(:, i)));
end
clear i
c = nchoosek(1:groupNum, 2);
t_steelDwass = zeros(length(c), 1);
p_steelDwass = zeros(length(c), 1);
stats_steelDwass = [];
for i = 1:length(c)
d1 = data(:, c(i, 1)); % group1 data
d2 = data(:, c(i, 2)); % group2 data
n1 = N(c(i, 1), 1); % group1 num
n2 = N(c(i, 2), 1); % group2 num
R = tiedrank([d1; d2]);
R2 = R.^2;
R2_sum = [sum(R2(1:length(d1), 1), 'omitnan'); sum(R2(length(d1)+1:length(R2), 1), 'omitnan')];
E = n1 * (n1 + n2 + 1) / 2;
x = (n1*n2) / ((n1 + n2)^2 - (n1 + n2));
y = sum(R2_sum) - (((n1 + n2) * (n1 + n2 + 1)^2)/4);
V = x * y;
t_steelDwass(i, 1) = abs((sum(R(1:length(d1), 1), 'omitnan') - E) / sqrt(V)); %%%%%%%%% t-value %%%%%%%%%
p_steelDwass(i, 1) = 1 - tcdf(t_steelDwass(i, 1), Inf); %%%%%%%%% p-value %%%%%%%%%
clear d1 d2 n1 n2 R R2 R2_sum E x y V
end
clear i
stats_steelDwass.df = groupNum;
stats_steelDwass.p = p_steelDwass;
stats_steelDwass.t = t_steelDwass;
stats_steelDwass.c = c;
clear p_steelDwass t_steelDwass c
stats_steelDwass.t
stats_steelDwass.p
0 Kommentare
Antworten (1)
Jeff Miller
am 15 Sep. 2021
The problem is with the 'inf' value given as the degrees of freedom in this line
p_steelDwass(i, 1) = 1 - tcdf(t_steelDwass(i, 1), Inf); %%%%%%%%% p-value %%%%%%%%%
I'm not familiar with the steel-dwass test, but I am pretty sure the df value should be n1+n2-2 (or something close to that) instead of Inf. If n1+n2-2 doesn't give you the right p values, you should be able to work out the correct df formula by changing that quantity (e.g., maybe n1+n2-1 or -3).
3 Kommentare
Jeff Miller
am 16 Sep. 2021
OK, then this function on file exchange is probably what you want: cdfTukey
Siehe auch
Kategorien
Mehr zu Hypothesis Tests finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!