I want to compare a parameter between two independent groups of volunteers (see MWE below). All volunteers had baseline reading (V1) and a subsequent reading by a later visit (V2) for intervention (
Interv; Yes/No). The parameters are gene expression data and expressed in fold change (i.e. V2/V1). For brevity, I mentioned only one
Par1. I came from here, but I didn't find what I want though how to report was updated in the question but it was not addressed yet I guess, so it might deserve a separate question like this one here.
To report the central tendency, I am intending to use the
median as it is less affected by outliers, may be with
MAD (Median Absolute Deviation), though not very often used, but I think is the equivalent to
SD for the
Subjects without intervention were (
n=15), whereas those with intervention were (
n=19). I am not a statistician, but I know this from my field that when it comes to gene expression data it is always good to log transform the data before doing analysis (to make the data mimic the normal distribution). So to report the
median I used the fold change reading because it is more intuitive to the reader and makes more sense than reporting log values.
To measure the difference between the two groups as inferential statistics, I am more inclined to use the
permutation techniques available in
R, like the
coin package, than to use
t-test, thinking that is more robust and less demanding of assumptions (confession: biologists do not respect these aliens, they will meet anybody except assumptions).
After log transformation, I did permutation test, wilcox test (Mann-Whitney-U-test), and t-test just for sake of comparison and learning.
The dataset (use
dget(<code below>) to get the dataset:
> dput(dfFold) structure(list(Interv = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Yes", "No"), class = "factor"), variable = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Par1", "Par2"), class = "factor"), Fold = c(181.552149430047, 365.287611324048, 1.37615543475921, 29.4749766533613, 7.30460900092265, 0.55729573732714, 0.856344617183982, 12.2151438678328, 8.22023432258123, 7.52091826273063, 11.3960901635291, 141.40637778502, 1.19646833351009, 4234.12741125762, 9410.18608208175, 0.11188317488791, 3.39142462873992, 1.54986235136038, 0.138010480189387, 0.150720706883034, 8.58182709551405, 1.65995914129576, 1.02664060360863, 40.9872498215584, 51.7412085216083, 3.64546611583102, 5.64546161492277, 8.92993632885674, 4.71129389070826, 1497.44850394305, 0.0209326180515642, 2.90945271783875, 5.09753127386217, 9.85421224939179 ), Log = c(5.20154293711853, 5.90068501952248, 0.31929369443907, 3.38354165434038, 1.98850551891357, -0.584659233229834, -0.15508239368987, 2.50267648260455, 2.10659871505853, 2.01768823991267, 2.43327032846195, 4.95163785700696, 0.179374162082163, 8.35093254365055, 9.14954800731289, -2.19030003343294, 1.22125007758044, 0.438166121413757, -1.98042565330713, -1.89232677812372, 2.14964683902473, 0.506792988388275, 0.0262919218996518, 3.71326103837802, 3.9462545340732, 1.29348423542785, 1.73085196838103, 2.18940926483575, 1.54996258170508, 7.31151794138033, -3.86644666427115, 1.0679649940249, 1.62875635858491, 2.28789900330443)), .Names = c("Interv", "variable", "Fold", "Log"), class = "data.frame", row.names = c(NA, -34L))
How it looks like:
> head(dfFold) Interv variable Fold Log 1 Yes Par1 181.5521494 5.2015429 2 Yes Par1 365.2876113 5.9006850 3 Yes Par1 1.3761554 0.3192937 4 Yes Par1 29.4749767 3.3835417 5 Yes Par1 7.3046090 1.9885055 6 Yes Par1 0.5572957 -0.5846592
To report the central tendency (I included
> median(dfFold[dfFold$Interv=="Yes",]$Fold,)  11.39609 > median(dfFold[dfFold$Interv=="No",]$Fold,)  3.645466 > mean(dfFold[dfFold$Interv=="Yes",]$Fold,)  960.8452 > mean(dfFold[dfFold$Interv=="No",]$Fold,)  86.71587
Now the the tests,
wilcox test from
> pvalue(oneway_test(Log~Interv, data=dfFold ,distribution=approximate(B=9999),conf.int=TRUE))  0.03590359 99 percent confidence interval: 0.03128298 0.04097196 Warning message: In independence_test.IndependenceProblem(object, check = check, : additional arguments conf.int will be ignored > set.seed(1234) > pvalue(wilcox_test(Log~Interv, data=dfFold,distribution=approximate(B=9999),conf.int=TRUE))  0.05890589 99 percent confidence interval: 0.05300214 0.06523398 > t.test (Log~Interv, data=dfFold) Welch Two Sample t-test data: Log by Interv t = 2.1519, df = 27.43, p-value = 0.04036 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 0.09772892 4.04392270 sample estimates: mean in group Yes mean in group No 3.183037 1.112211
I got a warning message from the
oneway_test(), anybody can help explain that?
So permutation test with approximation (B=9999) the difference was significant at 0.05 level, in agreement with the
t-test, but not with the
Now the question:
If I were to report this difference based on the permutation test what should I mention? Do I have to include the confidence intervals of 99% of
oneway_test()? (see the output), how to change it to 95%? and what is the difference in concept of statistic between
permutation a computer-intensive technique that can be used with many statistical tests? but when I do it with
oneway_test() what kind of test this function eventually uses? Any help to raise my statistical knowledge would be appreciated. Also few lines of comment on the result as one would be expected to write in a report would be much appreciated alike.