Comparison of medians in samples with unequal variance, size and shape

Question

What am I trying to do and how does my data look?

I am trying to perform some comparisons of medians of groups of data that, most often than not, are not-normally distributed. In some cases I am comparing some groups with 100-200 data points, but in others I'd like to compare groups with as little as 4 vs 10 data points if possible. Practically all of the compared groups have unequal sample size.

How am I doing this?

To this regard, I am doing the following:

1. Using the Kolmogorov–Smirnov test to compare the shape of the two samples (i.e. testing whether data1 exhibits a distribution that follows that of data2; In Mathematica: KolmogorovSmirnovTest[data1,data2]);
2. If shape is statistically similar, I use the Mann-Whitney U test to compare medians;
3. In most cases, groups do not have statistically similar shape. I therefore would like to perform a t-test with Welch correction (as suggested in this answer and in this wikipedia page). However, this has not been implemented in Mathematica and I currently do not have the time to do it myself.

I understand that by shifting to the Welch t-test I am no longer comparing medians but means, but I would do this if no other option is available.

Questions:

• Is there other test for comparing medians between samples with unequal variance, shape and size?
• What is the minimum sample size that I could take? I am aware that no minimal sample size exists for t-tests, because it depends on the power of the analysis. But maybe comparing shapes of distributions of 4 vs 10 points is not very good... And in such cases, the Central Limit Theorem should not have a very strong effect
• Anyone aware of such implementation in mathematica? Or maybe how to perform a Welch t-test there?

P.S. - I have searched a lot, including in this forum, but nothing shows up in this regard. I am truly sorry if this question was already answered before.

Additional Information:

In the comments it was asked some additional information. I hope this all makes sense, and I'm sorry if I cannot provide higher context into our real question.

• The information for each sample was collected from the genomic composition of many organisms. So I have two lists of independent samples in the form (0.01,0.7,0.15,...) whose distribution I do not expect to follow any particular trend (i.e. Normal, Log-normal, etc).
• For each organism we have the whole frequencies, and not only a sample, but for the groups to compare: we only have a small sample. Ideally, we would be able to collect information for some thousands of organisms, but this is not feasible or currently available. So, in our largest group we have a few hundreds of organism, in the smallest groups we have some tens.
• What we want to test is whether one population of organisms exhibits higher frequencies in one group than the other population. This because we want to get a location comparison, where the measure of comparison is representative of that population Quantile comparison, skewedness, standard deviation, etc, though informative, would possibly be not as good to get insight into how different the populations are (bearing in mind our particular biological question).So we want some general measure of comparison of both populations that permits simply saying stuff like "the two populations tend to not differ greatly in [this] overall measure" or "this population tends to exhibit higher median/mean frequency than the other".
• I am aware that means/medians per se are not informative enough about the distributions of both populations. For instance, if one is heavily tailed, they may both exhibit similar medians but one of them may exhibit higher means. However, for our purposes, we believe that medians (or means if necessary) provide enough insight.

Answer 1

I would suggest that you consider relative summary effects/relative treatment effects methods of Akritas, Arnold, Brunner, etc. The best book on the subject, generally, is Nonparametric Analysis of Longitudinal Data in Factorial Experiments by Brunner, Domhof, and Langer which is very well-written and clear.

These methods are based on modified notions of Mann-Whitney and can accommodate many of the features you mention using standard statistical software.

Answer 2

One (very conservative/low power) option would be to use Mood's median test. It does not assume normality or homogeneity of variance (or at least quite robust to unequal variance).

The test is done by making a cross tabulation of the number of observations less than or equal to the overall median vs. greater than the overall median for each group. From there, you run a chi-squared test on the table. For two group comparisons and when the groups are small, I'd use Fisher's exact test instead of chi-squared.

Edit: since I don't know if Mathematica has a function for Mood's median test, if you are willing to try using R, I have code that you can use.

Edit 2: Per Sosi's request, here's the code:

moods.median = function(v, f, exact = FALSE) {
#v is the set of Values you want to test
#f is a Factor or grouping variable
#if you want to use Fisher's exact test, specify "exact = TRUE"

if (length(v) != length(f)) {
    stop(c("v and f must have the same length"))
}



#make a new matrix data frame
m = cbind(as.character(f),v)
colnames(m) = c("group", "value")

#get the names of the factors/groups
facs = unique(f)

#count the number of factors/groups
factorN = length(unique(f))

if (factorN < 2) {
    stop("there must be at least two groups in variable 'f'")
}

#2 rows (number of values > overall median & number of values <= overall median)
#K-many columns for each level of the factor
MoodsMedianTable = matrix(NA, nrow = 2, ncol = factorN)

rownames(MoodsMedianTable) = c("> overall median", "<= overall median")
colnames(MoodsMedianTable) = c(facs[1:factorN])
colnames(MoodsMedianTable) = paste("Factor:",colnames(MoodsMedianTable))


#get the overall median
overallmedian = median(v)


#put the following into the 2 by K table:
for(j in 1:factorN){ #for each factor level

    g = facs[j] #assign a temporary "group name"

    #count the number of observations in the factor that are greater than
    #the overall median and save it to the table
    MoodsMedianTable[1,j] = sum(m[,2][ which(m[,1]==g)] > overallmedian)

    #count the number of observations in the factor that are less than
    # or equal to the overall median and save it to the table
    MoodsMedianTable[2,j] = sum(m[,2][ which(m[,1]==g)] <= overallmedian)

}

print(MoodsMedianTable)

if(exact == FALSE){return(chisq.test(MoodsMedianTable))}

if(exact == TRUE){return(list(
    chisq.test(MoodsMedianTable),
    fisher.test(MoodsMedianTable)))
}
}

Answer 3

My personal opinion is that you have the order of your tests backwards. If the Kolmogorov–Smirnov test shows that the samples come from a similar distribution, then use the t-test, otherwise use the Mann-Whitney U test.

For unequal variances, there is a procedure by Fligner and Policello (1981) who describe a methodology which can be employed for computing an adjusted U statistic if the homogeneity of variance assumption underlying the Mann-Whitney U test is violated. The reference I have for the Fligner-Policello test is Fligner, M. A. & Policello, II, G.E. (1981). Robust rank procedures for the Behrens-Fisher problem. Journal of the American Statistical Association, 76, 162-174.

the Welch's T-Test seems to be implemented in the EqualVariances->False option in the MeanDifferenceCl.

EqualVariances False whether the unknown population variances are assumed equal -- Option for MeanDifferenceCI.

• "Confidence intervals for the difference between means are also based on Student's distribution if the variances are not known. If the variances are assumed equal, MeanDifferenceCI is based on Student's distribution with Length[list1]+Length[list2]-2 degrees of freedom. If the population variances are not assumed equal, Welch's approximation for the degrees of freedom is used". See: http://reference.wolfram.com/language/HypothesisTesting/tutorial/HypothesisTesting.html