I want to know if there is a difference between the pH values of samples measured using two different methods (group 2 vs 1). Since the groups are paired and the values in each group are not normally distributed, I used a Wilcoxon signed-rank test.
library(dplyr) library(rstatix) # test for normality dat0_normality <- dat0 |> group_by(group) |> shapiro_test(ph)
# medians and CI5-95 dat0_medians <- dat0 |> group_by(group) |> summarise(nb = n(), median = round(quantile(ph, .50, na.rm=TRUE), 4), ci5 = round(quantile(ph, .05, na.rm=TRUE), 4), ci95 = round(quantile(ph, .95, na.rm=TRUE), 4))
# wilcoxon signed-rank test dat0_wilcox <- dat0 |> pairwise_wilcox_test(ph ~ group, paired = TRUE, detailed = TRUE, exact = TRUE) |> adjust_pvalue(method = "bonferroni") |> add_significance()
Here is my interpretation, perhaps inaccurate (I am not a statistician):
The absolute median of differences (AMD) between group 2 vs 1 is +0.008458209 (the estimate), i.e. a relative median of difference (RMD) of +0.116% (i.e. 100*0.008458209/7.2800) with a CI5-95 of +0.08% to +0.15%.
Given the low adjusted p-value of 1.26e-6, these AMD and RMD may be considered as significant.
However, from a clinical and analytical perspective, these AMD and RMD are negligible (the analytical acceptability threshold is +/-0.55%, five times wider than the RMD of 0.116%). Confused by this discrepancy between the highly significant p-value and the clinically/analytically negligible variation, I have two questions:
- Is the Wilcoxon signed-rank test appropriate for the objective?
- How can we explain the significance of the test given such close medians?
Boxplot for illustration:
Note: The effect size is large, but does this explain the low p-value?
library(coin) library(libcoin) dat0_effsize <- dat0 |> wilcox_effsize(ph ~ group, paired = TRUE)
Data:
dat0 <- structure(list(group = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), levels = c("1", "2"), class = "factor"), ph = c(7.343, 7.362, 7.323, 7.336, 7.411, 7.418, 7.318, 7.304, 7.216, 7.238, 7.271, 7.277, 7.265, 7.281, 7.282, 7.297, 7.28, 7.294, 7.221, 7.23, 7.255, 7.262, 7.207, 7.211, 7.421, 7.437, 7.265, 7.28, 7.315, 7.329, 7.332, 7.352, 7.322, 7.319, 7.33, 7.318, 7.338, 7.343, 7.395, 7.401, 7.366, 7.374, 7.247, 7.255, 7.361, 7.354, 7.39, 7.391, 7.392, 7.406, 7.289, 7.293, 7.338, 7.344, 7.383, 7.386, 7.38, 7.399, 7.241, 7.265, 7.257, 7.258, 7.014, 6.999, 7.41, 7.317, 7.28, 7.277, 6.955, 6.962, 7.253, 7.256, 7.163, 7.174, 7.284, 7.29, 7.221, 7.247, 7.109, 7.128, 7.444, 7.447, 7.412, 7.425, 7.376, 7.375, 7.189, 7.194, 7.23, 7.244, 7.21, 7.213, 7.14, 7.149, 7.312, 7.322, 7.277, 7.285, 7.352, 7.359, 7.32, 7.327, 7.22, 7.246, 7.219, 7.24, 7.125, 7.134, 7.021, 7.01, 7.33, 7.337, 7.087, 7.094, 6.917, 6.927, 6.932, 6.942, 6.954, 6.987 )), row.names = c(NA, -120L), class = c("tbl_df", "tbl", "data.frame" ))
Thanks for help
