How to get correct p-values for the Steel-Dwass test?
14 views (last 30 days)
Show older comments
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 Comments
Answers (1)
Jeff Miller
on 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 Comments
Jeff Miller
on 16 Sep 2021
OK, then this function on file exchange is probably what you want: cdfTukey
See Also
Categories
Find more on Hypothesis Tests in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!