How to get correct p-values for the Steel-Dwass test?

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

ans = 6×1

2.6802 2.5400 1.2826 3.7461 2.0468 3.3845

stats_steelDwass.p

ans = 6×1

0.0276 0.0320 0.1344 0.0100 0.0550 0.0138

Answer 1

[Jeff Miller]

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).

[T]

This statistical value "t" is studentized. So I would like to get the p-values from the studentized range cumulative distribution function. In this case, the number of the groups to compare is 4 (i.e. df); if I work on R, probabily the code is like the following. (But I need to work on MATLAB!)

p <- ptukey(t-value, groupNum, Inf, lower.tail = FALSE)

[Jeff Miller]

OK, then this function on file exchange is probably what you want: cdfTukey

[T]

The Tukey's HSD and the Steel-Dwass test are different. I could not get the correct p-value with the function.