3
$\begingroup$

I am running a multiple comparison using the non-parametric Kruskal Wallis test (K-W), using the ggpubr library and I am a bit confused about the results.

When i just run the KW-test using "base R" the result is different and I am not sure if there is an issue with the data or I am doing something wrong.

In case of ggpubr im getting significant result, i mean the p value is much lower than the critical threshold, but in case of base R this does not occer. I am posting the code and results figure.

Any suggestion would be highly appreciated why results are varying .. Data K-W test output 

library(ggpubr)
df1 = read.csv("NEW_RBP/CORPLOT/RBP_DISEASE_CLUSTER/RBP_DISEASE_C1_C4_AVG.txt",header = TRUE,sep = "\t")
head(df1)
my_comparisons <- list( c("HSC", "LSC"), c("LSC", "Blast"), c("HSC", "Blast") )

ex <- melt(df1, id.vars=c("gene"))

head(ex)

ggboxplot(ex, x = "variable", y = "value",ggtheme = theme_bw(base_size = 30),
          color = "variable", palette = "jco", add = "jitter")+ 
  stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
  stat_compare_means(label.y = 50)     # Add global p-value




kruskal.test(HSC ~ LSC, data = df1)

    Kruskal-Wallis rank sum test

data:  HSC by LSC
Kruskal-Wallis chi-squared = 416, df = 416, p-value = 0.4908
Improve this question
$\endgroup$
2
  • $\begingroup$ Could you provide the output of the base R solution? $\endgroup$
    – llrs
    Jan 2, 2019 at 8:19
  • 1
    $\begingroup$ i have added that to my question please have a look $\endgroup$
    – kcm
    Jan 2, 2019 at 9:02

3 Answers 3

Reset to default
0
$\begingroup$

The OP says KW is giving different results between "base R" and ggpubr. Background information: KW test is a non-parametric (assumption free) version of 1-way ANOVA.

The analysis you have presented provides a clear result: HSC and Blast are of an approximately equivalent distribution, however LSC is an outlier to both HSC and Blast.

The global value is very highly significant (because of LSC). In other words the global value (HSC, Blast and LSC used in a single analysis) is a much lower probability than the pairwise tests because it will amplify the discrepancy between LSC-HSC and LSC-Blast.

I personally don't see any discrepancy in behaviour within this analysis.

I suspect that HSC and Blast might be significantly different in different analyses because the result is quite low p=0.14, although it isn't borderline (e.g. if p=0.5 I would be happy). The basic issue in my opinion is that non-parametric tests are not particularly powerful and p=0.14 isn't hugely above p=0.05.

If this concerns you, I would examine the HSC and Blast groups using 1-way ANOVA and see if this sheds light on the problem, i.e. the global test will also be rejected by 1-way ANOVA.

1-way ANOVA requires the assumption of normality, but would be supportive of the KW-result. If 1-way ANOVA resulted in HSC and Blast being not significantly different from each other, that is your answer because 1-way ANOVA is a more powerful test (i.e, statistical power), thus it is more likely to reject a relationship between HSC and Blast.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ okay ,but im still not getting why i do see different result between the base R and ggpubr with the KW test? $\endgroup$
    – kcm
    Jan 2, 2019 at 5:11
  • 1
    $\begingroup$ Like I said, "I suspect that HSC and Blast might be significantly different in different analyses ....". I didn't know which analysis was in "base R". @llrs identified you have run a single data compartment in a multi-compartment test. If you are able to meaningfully compartmentalise your data (I don't know the data set), the HW-test will be more powerful than the Wilcox test. IMO its worth investigating because a lot of variables significantly different from each other isn't cool, but a hierarchy of differences is much more interesting. $\endgroup$
    – M__ ♦
    Jan 2, 2019 at 12:41
  • $\begingroup$ sounds more clear as my concepts about these are not so good ..i will read more about it $\endgroup$
    – kcm
    Jan 2, 2019 at 12:48
2
$\begingroup$

Your problem comes from using an incorrect test, the formula parameter says:

a formula of the form response ~ group where response gives the data values and group a vector or factor of the corresponding groups.

When you apply the kruskal test, you don't have groups in the LSC variables! You aren't comparing the mean of several groups based on LSC! When you use stat_compare_means it is doing a wilcox.test (it hints to it in the help page "a list of additional arguments used for the test method. For example one might use method.args = list(alternative = "greater") for wilcoxon test.").

When you use the wilcox.test we get the same p-value:

wilcox.test(df1$HSC, df1$LSC)

    Wilcoxon rank sum test with continuity correction

data:  df1$HSC and df1$LSC
W = 110980, p-value = 4.849e-12
alternative hypothesis: true location shift is not equal to 0

In addition, probably you should read about the paired argument, if this two values where measured on the same samples you should use it (but I can't tell it for sure because I don't know what does it mean each variable name).

Improve this answer
$\endgroup$
4
  • $\begingroup$ thanks for pointing out the blunder i as performing , i read in the manual that by default for multiple comparison it use kruskal test which i think i misunderstood due to lack of my proper statistical knowledge, you asked about variable name so HSC is my stem cell, LSC leukemic stem cell and Blast is the blast cell type and im trying to see if there is a significant difference or not ,may be this information would help,but why does it prints the kruskal test in the graph output im not sure though $\endgroup$
    – kcm
    Jan 2, 2019 at 12:07
  • 1
    $\begingroup$ Because the topleft value is comparing the three means where each group is each variable, which indeed requires the kruskal test, but not when comparing between two groups IMHO. $\endgroup$
    – llrs
    Jan 2, 2019 at 12:16
  • $\begingroup$ okay so the value reported in the figure is kruskal or Wilcoxon im confused or its hardcoded in the library ? $\endgroup$
    – kcm
    Jan 2, 2019 at 12:39
  • 1
    $\begingroup$ The one on the topright is kruskal, but to compare between two is wilcoxon. Read about them in the documentation of the test and the help page of the command. $\endgroup$
    – llrs
    Jan 2, 2019 at 13:12
2
$\begingroup$

Just to mention the KW test is appropriate here and not an error as suggested.

Each variable is a compartment and you have 3 compartments. There is no satisfactory equivalent to a non-parametric 2-way ANOVA (possible Friedman test), which is what is being implied.

Wilcoxon will only compare 2 samples, there is no restriction in the number under KW.

If you are comparing alot of samples, KW is much preferred because of type 2 error which will result from Wilcoxon.

I would present both results,KW and Wilcoxon. The only error was that we didn't know what test was run in base R.

Furthermore if you can split each sample into say two compartments, then you would run KW on each compartment (ie run it twice), at this point Wilcoxon becomes erroneous because you would need to run it 6 times at 0.05 critical probability and the chance of a false positive (two data sets are different when they are not) becomes statistically significant.

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ so what would you suggest ? mathematically and statistically to perform when i have comparison like this KW is the appropriate test ..meanwhile i wonder did you run code and data . $\endgroup$
    – kcm
    Jan 4, 2019 at 6:31
  • 1
    $\begingroup$ To keep it simple, both your code and results are fine. You present the results of both tests (pairwise Wilcoxon and KW). The key issue is that you must explain each test and a brief description of its rationale. You might add a statement that there is marginal support for LSC being an outlier. If a reviewer asks you to verify that last statement then post-back here. $\endgroup$
    – M__ ♦
    Jan 5, 2019 at 14:09
  • $\begingroup$ actually its kind of ongoing process my guide who is not a statistician by training so i have to make them get it first as they are not really concern about most of these mostly validation and they have not much idea so i have to get it right first then downstream $\endgroup$
    – kcm
    Jan 5, 2019 at 16:06
  • $\begingroup$ "The key issue is that you must explain each test and a brief description of its rationale. You might add a statement that there is marginal support for LSC being an outlier" stuff like this wont come to my head as i know how to use it without much idea about the ratioale ..immense help from guys like you makes my understanding better translate number into biological insight $\endgroup$
    – kcm
    Jan 5, 2019 at 16:09
  • 1
    $\begingroup$ Rationale, just text book theory (or wiki's), but you must mention "non-parametric". You know about statistics because you picked the right test, and that can't have been coincidence. Instead of using the term "outlier" just use the term "statistical difference". For example, "The distributions of LSC, HBS and Blast were very highly significantly different under the non-parametric KW test (state p value [describe KW in Methods]). However, using pairwise Wilcoxon tests a statistical difference was only significant for the "LSC" data against both the 'HBS' and 'Blast' groups (state p values)." $\endgroup$
    – M__ ♦
    Jan 6, 2019 at 2:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.