MATLAB Answers

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

8 views (last 30 days)
T
T on 14 Sep 2021
Commented: T on 17 Sep 2021
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

Answers (1)

Jeff Miller
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
T
T on 17 Sep 2021
The Tukey's HSD and the Steel-Dwass test are different. I could not get the correct p-value with the function.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!