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

19 次查看(过去 30 天)
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

回答(1 个)

Jeff Miller
Jeff Miller 2021-9-15
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 个评论
Jeff Miller
Jeff Miller 2021-9-16
OK, then this function on file exchange is probably what you want: cdfTukey
T
T 2021-9-17
The Tukey's HSD and the Steel-Dwass test are different. I could not get the correct p-value with the function.

请先登录,再进行评论。

Community Treasure Hunt

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

Start Hunting!

Translated by