Measures of central tendency and dispersion are crucial tools for characterizing data. The mean and standard deviation are commonly used in survey reports to visualize how the data is distributed. This procedure is a standard in academia and science, referring to the $\sigma$-rules of the normal distribution. We have gotten to know the arithmetic mean, the geometric equivalent to it and the related standard deviations, and variances in the previous chapters.
Following up, we will examine the proper use of these measures when reporting data.

Real-world example for the use of means: Suspension load in rivers

We will examine a data set on the suspension load concentration from real data from a gauge of the Rhine River at Diepoldsau. The data has been downloaded from the Gemstat portal. The units of our variables are: total suspension solid (tss) $mg/L$.

Loading the data

Let us read the data and save it as pandas data frame.

import pandas as pd
import numpy as np

data = pd.read_csv("https://userpage.fu-berlin.de/soga/data/workshop_data/Suspsolids.csv",sep = ";",comment = "#")
# next, we adress the variable types correctly:
data["Susp_solids_mg_L"] = data["Susp_solids_mg_L"].astype(str).str.replace(",", "")
data["Susp_solids_mg_L"] =pd.to_numeric(data["Susp_solids_mg_L"])
data["DateTime"] = pd.to_datetime(data["DateTime"],format="%d.%m.%y  %H:%M")   # here two spaces 

Next step: clean the data set, so that we can use it properly.

# There are some missing values to remove:
tss = data["Susp_solids_mg_L"][np.isfinite(data["Susp_solids_mg_L"])]
Q = data["Discharge_m3_s"][np.isfinite(data["Susp_solids_mg_L"])]

Next, we shall visualize the frequency distribution of the observed suspension loads.

tss.describe()

count     260.000000
mean      138.046154
std       195.447447
min         1.000000
25%        18.750000
50%        64.000000
75%       172.500000
max      1203.000000
Name: Susp_solids_mg_L, dtype: float64

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8,5))

ax.hist(tss, bins=30, color="lightblue",edgecolor = "grey")
ax.set_title("Frequency Distribution of suspension load of a river")
ax.set_xlabel("Suspension load mg/L")
ax.set_ylabel("Frequency")
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
plt.show()

The variable is strongly left-skewed (positive skewness). It ranges from $1$ to $1203$. The frequency of values increases the closer we get towards the lower limit of the scale. Logically, we would estimate the central tendency close to the lower part of the distribution as well.
To confirm this assumption, we want to calculate the arithmetic mean and median for this distribution.

fig, ax = plt.subplots(figsize=(8,5))

ax.hist(tss, bins=30, color="lightblue",edgecolor = "grey")
ax.set_title("Frequency Distribution of suspension load of a river")
ax.set_xlabel("Suspension load mg/L")
ax.set_ylabel("Frequency")
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.vlines(np.mean(tss),0,105, color = "red",label = "arithmetic mean" )
ax.vlines(np.median(tss),0,105, color = "blue",label = "median"  )
plt.legend(loc = 'center right')
plt.show()

Interesting! The difference between the median and the mean is quite large. Meanwhile, both values confirm that our data’s central tendency is towards the lower limit of the scale. The Median is a good reference for the mean, as it is very robust against outliers.

But the central tendency is only one half of the deal, as the data disperses around its central tendency.

While the range is already known, the dispersion can be described more detailed by the data's variance $s^2$ or its standard deviation $s$. Both values are calculated in relation to the mean.

We already know the arithmetic mean. It is most commonly used to report the central tendency in literature and textbooks. We can now calculate the standard deviation in relation to it.

The standard deviation provides further information, as it is closely linked to the Empirical rule, which we will discuss later.

By using the std() command, we can estimate the arithmetic standard deviation as a dispersion measure in relation to the arithmetic mean $\bar x$.

s = np.std(tss)

tss_mean = np.mean(tss)
tss_lower = (tss_mean - np.std(tss_mean))
tss_upper = (tss_mean + np.std(tss_mean))

fig, ax = plt.subplots(figsize=(8,5))
ax.hist(tss, bins=30, color="lightblue",edgecolor = "grey")
ax.set_title("Frequency Distribution of suspension load of a river")
ax.set_xlabel("Suspension load mg/L")
ax.set_ylabel("Frequency")
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.vlines(np.mean(tss),0,105, color = "red",label = "mean", lw= 2 )
ax.hlines(20,(tss_mean - np.std(tss)),(tss_mean + np.std(tss)), color = "green",label = r"$\pm$ 1. sd range", lw= 2)
plt.legend(loc = 'center right')
plt.show()

Hold on! According to our calculations, the data disperses around the mean in a range of $\pm$ 195.4 $mg/L$. This implies the sediment load would vary between -57.4 $mg/L$ and 333.5 $mg/L$!

The range exceeds the scale: A sediment load of 0 mg/L or less is physically not possible.

There are at least two reasons for this problem:

1. The arithmetic mean and standard deviation are meaningful exclusively for normally distributed data, which is not the case for our data.

2. The arithmetic mean and standard deviation are not suitable for data on a ratio scale, implying that the data is limited to $\mathbb{R_+}$ resp. $\mathbb{Q_+}$, which is the case for our data as well. Only multiplicative measures are unable to leave the positive branch, while additive measures can easily do so.

Therefore, we must report our data differently.

Geometric measures and real-world data

Naive data reporting using the commonly used $\bar x \pm s$ resp. $\mu \pm \sigma$ - interval may result in physically impossible values, as we have seen in our example. Assuming normality implicates interval scaled data, meaningful on the whole range of positive and negative rational/real numbers. Unfortunately, many real-world data sets are not normally distributed and are limited to the positive branch $\mathbb R_+$ or $\mathbb Q_+$. Examples include precipitation, height, weight, area, volume, and many more. The sediment concentration in rivers is just one of them.

In other words, real-world data is commonly limited to the $\mathbb{R_+}$ space (E. Limpert, 2001)

On the scale of positive real numbers $\mathbb{R_+}$ Multiplication and its inverse, division, cannot leave the set of positive real numbers. Thus, the positive real numbers together with the multiplication, denoted by $(\mathbb{R_+},*)$ forms a commutative group with 1 as neutral element. In contrast, the set of positive real numbers together with the addition $(\mathbb{R_+},+)$ has neither inverse elements nor a neutral element (because 0 is not a positive number).

The only mean exclusively designed for a multiplicative scale is the geometric mean. Together with the geometric standard deviation, they are the first moments of the log-normal distribution, which is the most common distribution for data on a ratio scale (D. McAllister, 1879).

Geometric vs. arithmetic measures

Let us recall the geometric mean: \begin{equation*} \bar{x}_{geo}=\sqrt[n]{\prod_{i=1}^N x_i}=\prod_{i=1}^N x_i^{\frac{1}{n}}= x_1^{\frac{1}{n}}\cdot x_2^{\frac{1}{n}}\cdots x_n^{\frac{1}{n}} \, . \end{equation*}

The main apparent difference between these two means is the operator: sum for the arithmetic mean and product for the geometric mean. Thus, the arithmetic mean is recommended for an additive or distance scale, and the geometric mean for a ratio scale approach. In our case, we are dealing with a variable clearly hosted on the $\mathbb{R}_+$-scale where only multiplication is permitted without limits.

But we often need an arithmetic mean for further statistical parameters: e.g.

• standard deviation
• skewness
• covariance
• correlation
• regression
• etc.

To solve this problem, we can apply the logarithm to the geometric mean: \begin{align*} log(\bar{x}_{geo}) & =log(x_1^{\frac{1}{n}}\cdot x_2^{\frac{1}{n}}\cdots x_n^{\frac{1}{n}})\\ & =\frac{1}{n}log(x_1)+\frac{1}{n}log(x_2)+\cdots+\frac{1}{n}log(x_N) \\ & =\frac{1}{n}\sum_{i=1}^{N}{log(x_i)}={\overline{log(x_i)}}_{arith} \end{align*}

!!! The logarithm of the geometric mean is the arithmetic mean of the logarithms!!!

Applying a logarithm to our sediment load data yields the following distribution:

log_tss = np.log(tss)

fig, ax = plt.subplots(figsize=(8,5))
ax.hist(log_tss, bins=30, color="lightblue",edgecolor = "grey")
ax.set_title("Rhine suspension load on log-scale")
ax.set_xlabel("log (mg/l)")
ax.set_ylabel("Frequency")
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
plt.show()

The arithmetic measures are limited to an open space scale with normally distributed values. They are not suitable for our data, as we have seen. The sediment load data is clearly non-normal and is limited to $\mathbb{R_+}$.

Thus, we will apply the log-transformation to the data to meet the preconditions for meaningful application of linear statistics.

A log-transformation of any variable $$x\in\mathbb{R_+}: x\to x'=\log(x)$$ opens the scale to $x'\to\mathbb{R}$ and transforms the neutral element of $(\mathbb{R_+},\cdot)$, given by 1, to the neutral element of $(\mathbb{R},+)$ that is given by 0!
The field $(\mathbb{R},+,\cdot)$ provides nearly unlimited operations for statistical purposes (except e.g. $\sqrt{-1}$ or $\frac{1}{0}$).

The inverse transformation powering $$x'\to x=e^{x'}\, .$$

Let us apply the geometric mean to our data and see how it performs in comparison to the arithmetic mean and the median.

gmean_tss = np.exp(np.mean(np.log(tss)))

fig, ax = plt.subplots(figsize=(8,5))
ax.hist(tss, bins=30, color="lightblue",edgecolor = "grey")
ax.set_title("Frequency Distribution of suspension load of a river")
ax.set_xlabel("Suspension load mg/L")
ax.set_ylabel("Frequency")
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.vlines(np.mean(tss),0,105, color = "red",label = "(arithmetic) mean", lw= 2 )
ax.vlines(np.median(tss),0,105, color = "blue",label = "median", lw= 2  )
ax.vlines(gmean_tss,0,105, color = "orange",label = "geometric mean", lw= 2  )
plt.legend(loc = 'center right')
plt.show()

The geometric mean is very close to the median. It is way smaller than the arithmetic mean and lies closer to the actual maximum of the distribution. Obviously, the geometric mean seems to be a far better estimate for the central tendency.

3. Applying geometric measures with log-transformation

We can use the log-transformation to project our data to the logistic scale. This enables us to perform all calculations needed in order to get to a meaningful dispersion interval relative to the geometric mean.

• When to use: when data are on a ratio scale, and the arithmetic mean and standard deviation are not appropriate measures of central tendency and variability, respectively. This is often the case when data are positively skewed (with right-hand tail) or have a multiplicative nature (rates of change, concentrations, etc.).

• Step 1: Check the distribution of the data using histograms, boxplots, or Q-Q plots to assess skewness and identify potential outliers. If the data has a lower limit greater than zero, shift the data by adding a constant to all values to ensure that all values are positive before applying the log-transformation.

• Step 2: Apply the log-transformation to the data using the formula $x' = \log(x)$, where $x$ is the original data and $x'$ is the transformed data. This can be done using software such as R or Python. For example, in R, you can use the log() function.

• Step 3: After log-transformation, check the distribution of the transformed data to ensure that it is more symmetric and closer to a normal distribution.

• Step 4: Calculate the arithmetic mean and standard deviation of the log-transformed data, which are now appropriate measures of central tendency and variability for the transformed data.

• Step 5: Back-transform the results to the original scale by applying the exponential function to the mean and standard deviation of the log-transformed data.

After back transformation, the mean will be the geometric mean of the original data, and the standard deviation will be the geometric standard deviation of the original data. Due to the direct multiplicative nature of your original data, you can now report the 1-sigma range of the original data as the geometric mean $^*/$ geometric standard deviation. Even applying higher sigma intervals, these ranges cannot exceed $\mathbb R_+$ resp. $\mathbb Q_+$ anymore.

Application for this Recipe: Let us apply the steps above to the aforementioned

• Step 1: Checking the distribution -We already checked the distribution of the data and found a strong skewness and a large variance. The data is meaningful on a ratio scale with a lower bound of 0. There is no obvious lower boundary greater than zero; we can skip the step of shifting the data by adding a constant to all values.

• Step 2: Apply the log-transformation to the data

• Step 3: Check the log-transformed distribution

log_tss = np.log(tss)

fig, ax = plt.subplots(figsize=(8,5))
ax.hist(log_tss, bins=30, color="lightblue",edgecolor = "grey")
ax.set_title("Frequency Distribution of logarithmic values")
ax.set_xlabel("Suspension load as log(mg/l)")
ax.set_ylabel("Frequency")
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
plt.show()

log_tss.describe()

count    260.000000
mean       4.083419
std        1.391661
min        0.000000
25%        2.930922
50%        4.158883
75%        5.150385
max        7.092574
Name: Susp_solids_mg_L, dtype: float64

The distribution of the log-transformed data is more symmetric, but now appears to be bi-modal. The latter characteristic was hidden in the original data histogram. However, we should expect a bimodal distribution in such an alpine river with a second-order run-off regime. We will keep in mind that our data is not perfectly log-normal distributed, so it does not form a bell-shaped normal distribution after log-transformation (instead, it is platykurtic and bimodal, but it is symmetric in log-space). We will continue with the log-transformed data as it is the best solution for our problem.

• Step 4: Calculate the arithmetic mean and standard deviation on the logistic scale Note: This time we are using the log-transformed values to estimate $\bar{x}$ and $s$.

m_log = np.mean(log_tss)
med_log = np.median(log_tss)
sd_log = np.std(log_tss)
onesigma_l_tss = m_log - sd_log
onesigma_u_tss = m_log + sd_log

• Step 5: Transform these parameters back by powering

geo_mean_tss = np.exp(m_log)
cor_median_tss = np.exp(med_log)
geo_sd_tss = np.exp(sd_log)
geo_onesigma_l_tss = np.exp(onesigma_l_tss)
geo_onesigma_u_tss = np.exp(onesigma_u_tss)
print(f"geomean = {geo_mean_tss:.1f}, geo_sd = {geo_sd_tss:.1f}, median = {cor_median_tss:.1f}, 1-sigma-interval = [{geo_onesigma_l_tss:.1f},{geo_onesigma_u_tss:.1f}]")

geomean = 59.3, geo_sd = 4.0, median = 64.0, 1-sigma-interval = [14.8,238.0]

Now, we can present our results, as we tried before:

Results:

1. Geometric mean $\bar x$ = 59 is close to the median (which has not changed by transformation)

2. Geometric standard deviation $\sigma_g$ = 4 which has to be interpreted as a multiplicative factor for the geometric mean, thus the 1-$\sigma$-interval is given by $\bar x_{geo} / s_{geo}$ and $\bar x_{geo} \cdot \sigma_{geo}$

3. The result of 2. can also be reported as 1- $\sigma$-interval: $\mu \;{^*/} \,\sigma$= [15, 239]

Thus, we not only get a reasonable mean, but the 1- $\sigma$-bounds are asymmetric as they should be for skewed distributions.

Finally, let us plot and report the moments correctly:

fig, ax = plt.subplots(figsize=(8,5))
ax.hist(tss, bins=30, color="lightgrey",edgecolor = "grey")
ax.set_title("Frequency Distribution of suspension load of a river")
ax.set_xlabel("Suspension load mg/L")
ax.set_ylabel("Frequency")
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.vlines(geo_mean_tss,0,105, color = "orange",label = "geometric mean" , lw= 2)
ax.hlines(20,geo_onesigma_l_tss,geo_onesigma_u_tss, color = "green",label = "*/ 1-sd range", lw= 2)
plt.text(150,25,f"{geo_mean_tss:.1f} */ {geo_sd_tss:.1f} mg/L \n1-$\\sigma$ interval = [{geo_onesigma_l_tss:.1f}, {geo_onesigma_u_tss:.1f}]")
plt.legend(loc = 'center right')
plt.show()

Perfect, we now have a meaningful mean and a meaningful 1-$\sigma$-interval for our data, which are both in the original units of measurement and can be interpreted in a meaningful way. The mean represents the central tendency of the data, while the 1-$\sigma$-interval represents the range of values that are within one standard deviation of the mean, accounting for the skewness of the original data.

4. The Empirical Rule

The need for geometric measures is clear by now. Nevertheless, we should discuss why the arithmetic $\bar x \pm s$ interval is still a useful measure for normally distributed variables and why it is used so commonly in literature. One of the main reasons is the empirical rule.

The empirical rule applies only to a specific type of distribution called a bell-shaped distribution or normal distribution

There are 3 rules:

For a bell-shaped distribution, approximately

1. 68% of the observations lie within one standard deviation of the mean.
2. 95% of the observations lie within two standard deviations of the mean.
3. 99.7% of the observations lie within three standard deviations of the mean.

Since we have sufficient coding abilities by now, we will try to test if the three rules are valid. We will explore the numpy.random.normal function to generate normally distributed data and display these rules using the generated data. The normal distribution belongs to the family of continuous distributions.

We can sample n values from a normal distribution with a given mean (default is 0) and standard deviation (default is 1) using the random.normal function: numpy.random.normal(loc=0.0, scale=1.0,size = None). Let us give it a try:

y_norm = np.random.normal(loc=0.0, scale=1.0,size = 100000)

If we plot a histogram of these numbers, we see the eponymous bell-shaped distribution.

fig, ax = plt.subplots(figsize=(8,5))

ax.hist(y_norm, bins=100, color="lightgrey",edgecolor = "grey")
ax.set_title("Normal distribution")
ax.set_ylabel("Frequency")
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
plt.show()

We already know the mean and the standard deviation of the y_norm vector, as we explicitly called the function numpy.random.normal() with loc = 0 and scale = 1. So, we just have to count the number of elements of the y_norm vector bigger than 1, and respectively smaller than -1, bigger than 2, respectively -2, and 3, respectively -3, and relate them to the length of the vector, in our case len(y_norm), to validate the three rules claimed above.

sd1 = sum((y_norm > -1) & (y_norm < 1)) / len(y_norm) * 100
sd2 = sum((y_norm > -2) & (y_norm < 2)) / len(y_norm) * 100
sd3 = sum((y_norm > -3) & (y_norm < 3)) / len(y_norm) * 100

print(f"1-sd: {sd1}\n2-sd: {sd2}\n3-sd: {sd3}")

1-sd: 68.193
2-sd: 95.50800000000001
3-sd: 99.725

Perfect match! The three empirical rules are obviously valid. To visualize our findings, we replot the histogram and add some annotations. Please note that in the hist() function, we set the argument density = True. As a consequence, the resulting histogram does not show counts on the y-axis anymore, but the density values (normalized count divided by bin width), which means that the bar areas sum to 1.

fig, ax = plt.subplots(figsize=(8, 5))

# create histogram
n, bins, patches = ax.hist(
    y_norm,
    bins=100,
    density=True,
    edgecolor="black")

# colour of bins
for patch, left, right in zip(patches, bins[:-1], bins[1:]):
    if right <= -3 or left >= 3:
        patch.set_facecolor("white")
    elif right <= -2 or left >= 2:
        patch.set_facecolor("purple")
    elif right <= -1 or left >= 1:
        patch.set_facecolor("green")
    else:
        patch.set_facecolor("orange")

# horizontal Lines
ax.plot([-2, 2], [0.48, 0.48], color="green", linewidth=3)
ax.plot([-3, 3], [0.55, 0.55], color="purple", linewidth=3)
ax.plot([-1, 1], [0.41, 0.41], color="orange", linewidth=3)

# vertikcal lines
for x, y, color in [
    (1, 0.41, "orange"), (-1, 0.41, "orange"),
    (2, 0.48, "green"), (-2, 0.48, "green"),
    (3, 0.55, "purple"), (-3, 0.55, "purple")
]:
    ax.plot([x, x], [0, y], color=color, linewidth=3)

# text
ax.text(0, 0.42, "68 %", fontsize=14, color="orange", ha="center")
ax.text(0, 0.49, "95 %", fontsize=14, color="green", ha="center")
ax.text(0, 0.56, "99.7 %", fontsize=14, color="purple", ha="center")

ax.set_title("Normal distribution")
ax.set_ylabel("Density")
ax.set_ylim(0, 0.6)
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
plt.show()

---

Citation

The E-Learning project SOGA-Py was developed at the Department of Earth Sciences by Annette Rudolph, Joachim Krois and Kai Hartmann. You can reach us via mail by soga[at]zedat.fu-berlin.de.

You may use this project freely under the Creative Commons Attribution-ShareAlike 4.0 International License.

Please cite as follow: Rudolph, A., Krois, J., Hartmann, K. (2023): Statistics and Geodata Analysis using Python (SOGA-Py). Department of Earth Sciences, Freie Universitaet Berlin.