Group mean/median comparison

Question

I have timeseries of concentrations in tropospheric air (NO, O3, etc.) and related meteorological parameters (Temperature, Wind Speed, etc.).

Some of those time series correlate well (eg. CO2 vs NO) and none of them have a normal distribution or even symmetric distribution.

Where correlation coefficients are:

measurandkey       CO2        NO         T      WS-s
measurandkey                                        
CO2           1.000000  0.620874 -0.561119 -0.341159
NO            0.620874  1.000000 -0.237885 -0.246442
T            -0.561119 -0.237885  1.000000 -0.044186
WS-s         -0.341159 -0.246442 -0.044186  1.000000

One of my goal is to assess the impact of:

• Heating (linked with external temperature);
• Troposphere Stability (CO2 can be thought as a passive tracer linked with combustion sources and stability);
• Horizontal Shear Stress (linked with wind speed).

On concentration of a pollutants directly linked with source of combustion (mostly traffic and heating) such as Nitrogen Oxide.

With T, CO2 (and Wind Speed) timeseries, I can create few categories based on some clear cut criterion:

• Shear Stress is categorized before its distribution and mechanical criteria (definition of a calm wind, etc.), leading to four categories;
• Stability is categorized before CO2 seasonal residuals distribution (CO2 seasonal trend is removed, then split into three categories);
• Heating is categorized before Degree-Day distribution (3 days smoothing triangle window is applied on temperature);

Below a subsample of data and groups created from criteria above:

measurandkey           CO2    NO      T   WS-s Shear Stress   CO2 Trend  \
timevalue                                                                 
2012-01-09 08:00:00  423.5  12.0   7.45  2.925          low  426.265359   
2012-01-09 18:00:00  417.5   2.5   8.65  2.880          low  426.214362   
2012-01-10 04:00:00  425.0   2.5   8.50  1.645         null  426.163365   
2012-01-10 14:00:00  420.5  13.0   9.55  2.035          low  426.112368   
2012-01-11 00:00:00  416.0   NaN   8.40  1.825         null  426.061371   
2012-01-11 10:00:00  411.0   7.5  10.05  3.740          low  426.010374   
2012-01-11 20:00:00  410.5   2.0   8.90  3.480          low  425.959377   
2012-01-12 06:00:00  413.0   3.0   7.15  5.580          med  425.908380   
2012-01-12 16:00:00  410.0   2.0   9.35  4.410          low  425.857383   
2012-01-13 02:00:00  422.5   2.0   2.35  2.285          low  425.806387   

measurandkey         CO2 Residual Stability      T_deg Heating    NO Log  
timevalue                                                                 
2012-01-09 08:00:00     -2.765359       med  11.476601    high  2.484907  
2012-01-09 18:00:00     -8.714362       med  11.470853    high  0.916291  
2012-01-10 04:00:00     -1.163365       med  11.429398    high  0.916291  
2012-01-10 14:00:00     -5.612368       med  11.076427    high  2.564949  
2012-01-11 00:00:00    -10.061371       low  10.352276    high       NaN  
2012-01-11 10:00:00    -15.010374       low   9.804668     med  2.014903  
2012-01-11 20:00:00    -15.459377       low   9.374306     med  0.693147  
2012-01-12 06:00:00    -12.908380       low   9.208738     med  1.098612  
2012-01-12 16:00:00    -15.857383       low   9.279360     med  0.693147  
2012-01-13 02:00:00     -3.306387       med   9.401215     med  0.693147 

Then I am able to compute basic statistics such as mean, median, count, etc. on pair of group modalities, which leads to the following contingency matrix:

          count                   mean                      median           
Stability   low    med  high       low       med       high    low  med  high
Heating                                                                      
null        348    610    86  2.278736  2.542623   3.517442    2.0  2.0   2.0
low        2090   2793   452  2.284450  2.999821   6.263274    2.0  2.0   2.0
med        7013  11136  1223  2.667902  4.165365  23.881848    2.0  2.0  11.5
high       3393   8423  1196  2.680519  6.720052  51.109114    2.0  3.0  42.0

What I see is:

• for some combination of group modalities, there are few or no data, and therefore some statistical quantities may be meaningless;
• grouped sub samples often differ in term of distributions, mean and median (chi-by-eye on following figures).

We can render the above table as a Bar chart:

Or grouped distributions with Whisker Boxes:

Grouping can also be seen on scatter plot:

Goal

What I am looking for is: a more objective criterion in order to confirm that my appreciation (eg. mean are different among groups) is not misled by some visual effect or a priori on groups.

I don't know the distribution from which my original sample is drawn. I know that the distribution is strictly positive (because it is concentrations) and highly skewed (often the case with primary pollutants and measurement sites close to sources). I cannot assume normality.

What I want to assess is:

• Are mean (or median) of pair of group modalities significantly different?
• Are there enough observations in each pair of group modalities in order to perform reliable statistics on it?

I am looking for references and keywords pointing to procedures that could achieve those tasks, then I could implement post-processing. Thank you for your help.

As AdamO pointed out that I may look for cross correlation. I think mainly because I could work with continuous random variables instead of discretizing them. But here I decided to work with categorical variables and this why I have created this post. I am open to any constructive feedback about this decision. Anyway it could happens in real life problem that you want to compare continuous variable grouped over categorical groups.

Answer 1

I'm not sure what is meant by "a more objective criterion" as multiple metrics are possible and the decision of which one(s) to use can be motivated by a number of considerations. In other words, there are more than a few ways to analyze your information.

First of all, with strictly positive information truncated at zero, gaussian assumptions are clearly not relevant. Given that, the arithmetic mean is not appropriate. Next, that the pdf(s) of your data are skewed or fat-tailed also has important consequences. Putting these two observations together leads one to an initial assumption that the information may be lognormally distributed. Given that, estimating a lognormal mean would be more appropriate. This wiki article discusses approaches for that (https://en.wikipedia.org/wiki/Log-normal_distribution). It may also be the case that with time series data such as this, the lognormal geometric mean would be a good estimator (also discussed in the wiki link).

However, lognormality should only be an initial assumption and subject to some confirmation. One useful visual heuristic supporting (disconfirming) lognormality is the extent to which the empirical data fits a 45 degree, straight line in a Q-Q plot. The key thing to look for is the fit (or lack of it) in the extreme tails. If there is evidence for "lack of fit," deeper analyses are possible to help ferret out just how skewed or fat-tailed your data is. For instance, estimating a tail exponent would provide empirical evidence for fitting distributions other than lognormal, e.g., fat-tailed, extreme valued, power law distributions. A simple approach to estimating a tail index is discussed in this paper by Xavier Gabaix, Rank 1/2: A Simply Way to Improve the OLS Estimation of Tail Exponents (http://www.nber.org/papers/t0342.pdf). Basically, by taking the log of the metric as well as the log of its rank order (an ordinal, integer value), and using log-log OLS regression of the one on the other, the absolute value of the resulting beta coefficient is a heuristic proxy for a tail exponent. The closer that exponent is to 1.0, the greater the likelihood that the data is lognormally distributed. For power law distributions more extreme than lognormal, the "Examples" table from this wiki article discussing Tweedie distributions (https://en.wikipedia.org/wiki/Tweedie_distribution#Examples) can be used to identify and fit a variety alternative distributions.

However, if the tail exponent from this simple heuristic is much greater than 1.0, then consequences need to be considered. First, as the magnitude of the exponent diverges from 1.0, the less reliable simple log-log OLS estimation becomes. In other words, as your data becomes increasingly extreme valued, lognormal (log-log) assumptions become increasingly unreliable estimators of the real magnitude or nature of these extreme values. At some point, more rigorous, accurate maximum likelihood methods for tail exponent estimation such as Hill's or Pickands' methods, need to be employed. Next, as the information conforms more and more to extreme valued, power law assumptions, it may be the case that the moments (mean, std dev, etc.) are infinite and/or undefined. For instance, the Cauchy distribution is one example of a symmetric, extreme value, power law distribution with a single parameter -- the mean -- its std dev is undefined and/or infinite. Two things about the Cauchy, it's not well represented by an arithmetic mean and, next, since it is symmetric around that mean, as Gelman has suggested for use with data like yours, employing a "half-Cauchy" to capture its greater than zero pdf would be recommended.

All of this is likely confusing to you, at least wrt estimating the mean of your strictly positive, skewed data. Hopefully, your data provides an adequate fit to lognormal assumptions, making your life much easier. If not and given that mean estimation is so fraught with complications, it would be much simpler and more straightforward to estimate the median(s) and employ nonparametric methods for between group comparisons. The Mann-Whitney U test is for use when comparing two groups while the Kruskal-Wallis test, a generalization of the Mann-Whitney, can be used with multiple groups. Unfortunately, neither test was designed for use with time series data and would, therefore, require some simplifying assumptions wrt the presence of the classic issues wrt time series modeling: nonstationarity, autocorrelation, cointegration, and so on, to make them appropriate tests.

Other CV participants may have better recommendations and suggestions wrt dealing with all of the statistically focused, time series issues your data presents, particularly in the case where lognormal assumptions no longer apply. I would be interested in those suggestions. My recommendation would be to leverage information-theoretic approaches that are not based or rooted in the moments of any statistical distribution and are appropriate for use with time series data. One good example of such an approach is permutation distribution clustering (PDC). There are a number of articles and R modules out there which discuss this method, most authored by Andreas Brandmaier at the Max Planck Institute, e.g., https://cran.r-project.org/web/packages/pdc/pdc.pdf. Here is the abstract to this paper:

Permutation Distribution Clustering (pdc) represents a complexity-based approach to clustering time series. Clustering comprises methods that recover similarities in a dataset and represent the findings in group structures. Important applications of clustering include the creation of taxonomies, the discovery of anomalies, or the the discovery of reliably different subgroups for differential analysis or treatment. A crucial parameter in clustering is the choice of the similarity measure between objects. Permutation Distribution Clustering finds similarity in time series based on differences in their permutation distribution as a proxy for differences in their complexity. The permutation distribution is obtained by counting the frequency of distinct order patterns in an m-embedding of the original time series. An embedding of dimension m allows for m! different order patterns. The choice of the embedding dimension crucially influences the clustering result. A small embedding dimension might lead to a permutation distribution with a low representational power, while a large embedding dimension quickly leads to a large permutation distribution that cannot reliably be estimated. With the Minimum Entropy Heuristic (MinE), the embedding dimension can automatically be chosen, thus making the algorithm a parameter-free clustering approach. For clustering timeseries, the similarity between two time-series is defined as the divergence between two permutation distributions. PDC is particularly apt for the analysis of psychophysiological time-series because it is efficient (the time complexity is linear in the time-series length), it is robust to drift, time-series of differing length can be compared, and it is invariant to differences in mean and variance of the time-series (choosing a normalization is not essential).

PDC creates a matrix of comparisons based on the Kullback-Leibler divergence metric as a measure of similarity. If one so chooses, at that point, standard statistical tests of multiple group differences can be employed, e.g., Student-Newman-Keuls or some other multiple group, ANOVA-type test for contrasts between matrix elements.