Типичные задачи аналитика. Часть 2. А есть ли тренд?

Введение

В первой части статьи на Habr мы рассмотрели классические подходы к оценке изменений метрики при условии ее стационарности. В этом контексте статистические критерии, применяемые в A/B тестировании, оказались весьма эффективными.

Однако, если существует стабильный тренд, например, среднемесячная аудитория увеличивается из года в год, оценка разницы средних за два смежных периода времени может быть некорректной. В таком случае среднее значение предыдущего периода всегда будет отличаться от среднего постпериода, и это часто может быть не связано с исследуемым функционалом.

Одна из причин — тренд не всегда зависит от действий компании и часто является следствием внешних условий. Например, рост аудитории может быть связан с увеличением благосостояния населения, масштабированием бизнеса или сезонными факторами.

Таким образом, наличие или отсутствие тренда является важным аспектом анализа данных. Рассмотрим несколько успешных и неудачных подходов, которые можно применять для решения этой задачи.

Визуальный анализ

В контексте анализа временных рядов визуальный анализ играет важную роль в определении тренда, но у него есть очевидные недостатки и ограничения. Субъективность восприятия, неопределенность результатов и ограничения в обработке больших объемов данных - лишь некоторые из аспектов, которые стоит учитывать при использовании этого метода.

Визуальный анализ полезен для первичной оценки данных, выявления общих паттернов и трендов, а также для выбора подходящего метода дальнейшего анализа. Он помогает понять тип зависимости в данных, выделить основные характеристики временного ряда и сформулировать предварительные гипотезы.

Ниже рассмотрим несколько примеров:

Код Python

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns


def generate_data(
    n_samples: int,
    loc: int = 100,
    scale: int = 3,
    trend_max: int = 20,
    hetero: float = 1,
    n_outliers: int = 0,
) -> pd.DataFrame:
    time = pd.date_range("2024-01-01", periods=n_samples, freq="1H")
    metric_values = (
        loc
        + np.random.normal(scale=scale, size=n_samples)
        * np.linspace(1, hetero, n_samples)
        + np.linspace(0, trend_max, n_samples)
    )

    for i in range(n_outliers):
        metric_values[i * n_samples // n_outliers] = np.max(metric_values) * 1.02

    return pd.DataFrame({"time": time, "metric": metric_values})


def rolling_data(data, window_size=10):
    return data.rolling(window=window_size).mean()


ONE_MONTH = 24 * 30
test_data_1 = generate_data(ONE_MONTH, hetero=1, trend_max=20, scale=3)
test_data_2 = generate_data(ONE_MONTH, hetero=1, trend_max=20, scale=10)
smoothed_data = rolling_data(test_data_2["metric"])

fig, axes = plt.subplots(3, 1, figsize=(10, 12))

sns.lineplot(x="time", y="metric", data=test_data_1, ax=axes[0], alpha=0.5)
axes[0].set_ylim(0, 200)
axes[0].set_title("Low variance")

sns.lineplot(x="time", y="metric", data=test_data_2, ax=axes[1], alpha=0.5)
axes[1].set_ylim(0, 200)
axes[1].set_title("High variance")

axes[2].plot(test_data_2["time"], smoothed_data, alpha=0.5, color="red")
axes[2].set_ylim(0, 200)
axes[2].set_title("High variance and rolling window")

plt.subplots_adjust(hspace=0.5)
plt.show()

В случае, когда дисперсия в данных низкая, обнаружение тренда становится относительно простым, так как данные имеют менее выраженные колебания, что упрощает его визуализацию. Однако в более типичных сценариях, когда данные шумные, как на втором графике, выделение определенного тренда становится проблематичным. В таких ситуациях можно попытаться уменьшить влияние шума с помощью методов сглаживания, таких как скользящее среднее или экспоненциальное сглаживание. Эти методы также могут иметь свои недостатки и привести к ошибкам в анализе. Например, скользящее среднее, основываясь на предыдущих значениях, может пропустить недавно возникший тренд или изменение, что снизит его чувствительность к новым данным.

Стоит помнить, что визуальный анализ должен использоваться в сочетании с другими методами анализа, так как только на его основе делать окончательные выводы может быть недостаточно надежно.

Линейная регрессия

Линейная регрессия это один из наиболее популярных и эффективных методов для определения и оценки тренда. Она является неотъемлемой частью курсов по математической статистике и анализу данных.

Одно из главных преимуществ линейной регрессии является ее простота, интерпретируемость, малые требования к вычислительным ресурсам, а также ее эффективность при линейных взаимосвязях.

Однако последний плюс в тоже время является и минусом - линейная регрессия может быть недостаточно гибкой, чтобы хорошо аппроксимировать сложные взаимосвязи в данных. Если тренд имеет нелинейную форму, линейная регрессия может построить неверную модель.

Линейная регрессия предполагает независимость ошибок и их одинаковое распределение (гомоскедастичность). Нарушение этих предположений может привести к некорректным выводам.

Построение регрессии чаще всего происходит по следующим этапам

• Подготовка датасета.

• Применение линейной регрессии.

• Проверка гомоскедастичности.

  • Гомоскедастичность означает, что дисперсия ошибок модели постоянна для всех значений предикторов. Иными словами, разброс ошибок не зависит от уровня предсказываемой переменной. Почему это важно делать можно подробнее почитать тут. В данном примере линейная зависимость, поэтому будем использовать Breusch-Pagan test.

• Вычислить R² и доверительный интервал угла наклона.

  • R² в данном случае показывает качество модели. Чем он ближе к 1, тем большее кол-во данных объясняется моделью. Но не всегда низкие значения R² свидетельствуют о несостоятельности модели. В случае когда данные очень шумные и имеются выбросы, R² будет как правило довольно низким, но это не отрицает наличие линейных зависимостей или явного тренда. Детальней рассмотрим ниже.

• Проверка на независимость ошибок

  • Независимость ошибок в модели линейной регрессии важна для получения корректных оценок параметров. Для этого может применяться статистика, такая как тест Дарбина-Уотсона, который проверяет наличие автокорреляции в остатках модели. Подробней с ней можно познакомиться здесь. Обычно, когда статистика Дарбина-Уотсона находится в пределах 1.5 до 2.5, это считается приемлемым показателем отсутствия автокорреляции. Однако, для точной интерпретации, можно использовать таблицы критических значений для теста Дарбина-Уотсона или статистические пакеты, которые автоматически вычисляют p-значение для данной статистики.

• Вычисление доверительного интервала наклона с помощью бутстрапа.

  • О бутстрапе подробно было описано в предыдущей статье. Здесь лишь уточним что такой метод нужен для того чтобы избежать случайных результатов в данных с широким интервалом. Построение 1000 моделей с помощью бутсрапа позволит оценить интервал в рамках которого может изменять угол наклона нашего тренда.

Рассмотрим далее несколько примеров применения линейной регрессии для определения тренда на Python и Scipy , а также основные трудности с которым можно столкнуться при построении линейной регрессии.

Построим для примера модель в которой выполняются все условие.

Код Python

import matplotlib.pyplot as plt

import numpy as np

import pandas as pd
import scipy.stats as stats

import seaborn as sns
import statsmodels.api as sm


def generate_data(
    n_samples: int,
    loc: int = 100,
    scale: int = 10,
    trend_max: int = 10,
    hetero: float = 2,
) -> pd.DataFrame:
    time = pd.date_range("2024-01-01", periods=n_samples, freq="1H")

    metric_values = (
        loc
        + np.random.normal(scale=scale, size=n_samples)
        * np.linspace(1, hetero, n_samples)
        + np.linspace(0, trend_max, n_samples)
    )

    return pd.DataFrame({"time": time, "metric": metric_values})


ONE_MONTH = 24 * 30
test_data = generate_data(ONE_MONTH)

# plot data

sns.lineplot(x="time", y="metric", data=test_data, alpha=0.5)
# fit linear regression
slope, intercept, rvalue, pvalue, stderr = stats.linregress(
    test_data.index, test_data.metric
)

# test for Homoscedasticity of the residuals


def test_homoscedasticity(test_data: pd.DataFrame):
    # prepare y and X for statsmodels
    y = test_data.metric
    X_sm = test_data.index
    X_sm = sm.add_constant(X_sm)
    # fit linear regression using statsmodels
    mod_sm = sm.OLS(y, X_sm)
    res_sm = mod_sm.fit()
    # test for Homoscedasticity
    bp_lm, bp_lm_pvalue, bp_fvalue, bp_f_pvalue = sm.stats.diagnostic.het_breuschpagan(
        res_sm.resid, res_sm.model.exog
    )

    return bp_f_pvalue


bp_f_pvalue = test_homoscedasticity(test_data)
if bp_f_pvalue < 0.05:
    print("Residuals are not Homoscedastic. Linear regression may be invalid")
else:
    print("Residuals are Homoscedastic")


# test autocorrelation
def test_autocorrelation(residuals):
    dw_statistic = sm.stats.durbin_watson(residuals)
    return dw_statistic


# Assuming you have the 'residuals' variable defined elsewhere in your code
# Conduct the test for autocorrelation


# Check the p-value and print the result
residuals = test_data.metric - (intercept + slope * test_data.index)
dw_statistic = test_autocorrelation(residuals)

if dw_statistic > 1.5 and 2.5:
    print("Errors are independent (no signs of autocorrelation)")
else:
    print(
        "Аutocorrelation in the residuals. Independence of errors is not satisfactory."
    )

# bootstrap linear regression slope
n_iterations = 1000
slopes = np.zeros(n_iterations)
intercepts = np.zeros(n_iterations)
r_square = np.zeros(n_iterations)
for i in range(n_iterations):
    sample = test_data.sample(frac=1, replace=True)
    slope, intercept, rvalue, pvalue, stderr = stats.linregress(
        sample.index, sample.metric
    )
    slopes[i] = slope
    intercepts[i] = intercept
    r_square[i] = rvalue**2
# calculate 95% confidence interval
slopes_quantiles = np.percentile(slopes, [2.5, 97.5])
intercept_quantiles = np.percentile(intercepts, [2.5, 97.5])
r_square_quantiles = np.percentile(r_square, [2.5, 97.5])
print(f"r_square is between {r_square_quantiles}")
print(f"slope is between {slopes_quantiles}")


# plot two lines for the confidence interval slopes
plt.plot(
    test_data.time,
    intercept_quantiles[0] + slopes_quantiles[0] * test_data.index,
    label="lower bound",
    color="green",
    linestyle="--",
    linewidth=2,
)

plt.plot(
    test_data.time,
    intercept_quantiles[1] + slopes_quantiles[1] * test_data.index,
    label="upper bound",
    color="green",
    linestyle="--",
    linewidth=2,
)

# plot linear regression

plt.plot(
    test_data.time,
    intercept + slope * test_data.index,
    label="fitted line",
    color="red",
    linewidth=2,
)

plt.legend()
plt.show()

Полученные результаты указывают на присутствие тренда в данных, с углом наклона (slope) в интервале от ~0.012 до 0.018. Другими словами, каждый день метрика увеличивается от 0.012 до 0.018 единиц. Проверка гомоскедастичности подтверждает выполнение условий для построения регрессии, а анализ независимости ошибок показывает, что различия не являются статистически значимыми.

Таким образом, полученные результаты дают представление о характере тренда и его величине, что затруднительно сделать только с помощью визуального анализа.

Проверка на гомоскедастичность.

Наличие высокой гетероскедастичности создает препятствия для построения модели. Давайте рассмотрим, какие последствия могут возникнуть при построении линейной регрессии при игнорировании этой проверки.

Код Python

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.stats as stats
import seaborn as sns
import statsmodels.api as sm


def generate_data(
    n_samples: int,
    loc: int = 100,
    scale: int = 10,
    trend_max: int = 10,
    hetero: float = 30,
) -> pd.DataFrame:
    time = pd.date_range("2024-01-01", periods=n_samples, freq="1H")

    metric_values = (
        loc
        + np.random.normal(scale=scale, size=n_samples)
        * np.linspace(1, hetero, n_samples)
        + np.linspace(0, trend_max, n_samples)
    )

    return pd.DataFrame({"time": time, "metric": metric_values})


ONE_MONTH = 24 * 30
test_data = generate_data(ONE_MONTH)

# plot data

sns.lineplot(x="time", y="metric", data=test_data, alpha=0.5)
# fit linear regression
slope, intercept, rvalue, pvalue, stderr = stats.linregress(
    test_data.index, test_data.metric
)

# test for Homoscedasticity of the residuals


def test_homoscedasticity(test_data: pd.DataFrame):
    # prepare y and X for statsmodels
    y = test_data.metric
    X_sm = test_data.index
    X_sm = sm.add_constant(X_sm)
    # fit linear regression using statsmodels
    mod_sm = sm.OLS(y, X_sm)
    res_sm = mod_sm.fit()
    # test for Homoscedasticity
    bp_lm, bp_lm_pvalue, bp_fvalue, bp_f_pvalue = sm.stats.diagnostic.het_breuschpagan(
        res_sm.resid, res_sm.model.exog
    )

    return bp_f_pvalue


bp_f_pvalue = test_homoscedasticity(test_data)
if bp_f_pvalue < 0.05:
    print("Residuals are not Homoscedastic. Linear regression may be invalid")
else:
    print("Residuals are Homoscedastic")


# test autocorrelation
def test_autocorrelation(residuals):
    dw_statistic = sm.stats.durbin_watson(residuals)
    return dw_statistic


# Assuming you have the 'residuals' variable defined elsewhere in your code
# Conduct the test for autocorrelation


# Check the p-value and print the result
residuals = test_data.metric - (intercept + slope * test_data.index)
dw_statistic = test_autocorrelation(residuals)

if dw_statistic > 1.5 and 2.5:
    print("Errors are independent (no signs of autocorrelation)")
else:
    print(
        "Аutocorrelation in the residuals. Independence of errors is not satisfactory."
    )

# bootstrap linear regression slope
n_iterations = 1000
slopes = np.zeros(n_iterations)
intercepts = np.zeros(n_iterations)
r_square = np.zeros(n_iterations)

for i in range(n_iterations):
    sample = test_data.sample(frac=1, replace=True)
    slope, intercept, rvalue, pvalue, stderr = stats.linregress(
        sample.index, sample.metric
    )
    slopes[i] = slope
    intercepts[i] = intercept
    r_square[i] = rvalue**2
# calculate 95% confidence interval
slopes_quantiles = np.percentile(slopes, [2.5, 97.5])
intercept_quantiles = np.percentile(intercepts, [2.5, 97.5])
r_square_quantiles = np.percentile(r_square, [2.5, 97.5])
print(f"r_square is between {r_square_quantiles}")
print(f"slope is between {slopes_quantiles}")


# plot two lines for the confidence interval slopes
plt.plot(
    test_data.time,
    intercept_quantiles[0] + slopes_quantiles[0] * test_data.index,
    label="lower bound",
    color="green",
    linestyle="--",
    linewidth=2,
)

plt.plot(
    test_data.time,
    intercept_quantiles[1] + slopes_quantiles[1] * test_data.index,
    label="upper bound",
    color="green",
    linestyle="--",
    linewidth=2,
)


# plot linear regression

plt.plot(
    test_data.time,
    intercept + slope * test_data.index,
    label="fitted line",
    color="red",
    linewidth=2,
)

plt.legend()
plt.show()

Тест Breusch-Pagan сигналиризует о том что разброс ошибок зависим от уровня предсказываемой переменной

Результаты могут быть примерно такими:

Получаемые результаты при анализе данных могут быть непостоянными, и один и тот же код может демонстрировать как положительный, так и отрицательный тренд. В таких случаях выводы о наличии тренда остаются на усмотрение случайных факторов. Более подробную информацию о том, как бороться с гетероскедастичностью, можно найти здесь, однако решение этой проблемы часто бывает сложным, и определить явное наличие тренда в данных может быть затруднительно.

Оценка R²

Внимательный читатель мог заметить, что в первом примере была пропущена оценка коэффициента детерминации R². Возвращаясь к этому, результаты анализа показывают, что R² находится в диапазоне между 0.06 и 0.14. Это свидетельствует о том, что модель довольно слабо описывает данные, что, в целом, подтверждает наше предположение о том, что такая модель характерна для данных с высоким уровнем шума.

Ниже рассмотрим два дополнительных примера, иллюстрирующих, почему низкий R² не всегда свидетельствует о плохом качестве модели, а может быть следствием особенностей данных.

Код Python

import matplotlib.pyplot as plt

import numpy as np

import pandas as pd
import scipy.stats as stats

import seaborn as sns
import statsmodels.api as sm


def generate_data(
        n_samples: int,
        loc: int = 100,
        scale: int = 30,
        trend_max: int = 100,
        hetero: float = 1
) -> pd.DataFrame:
    time = pd.date_range("2024-01-01", periods=n_samples, freq="1H")

    metric_values = (
            loc
            + np.random.normal(scale=scale, size=n_samples)
            * np.linspace(1, hetero, n_samples)
            + np.linspace(0, trend_max, n_samples)
    )

    return pd.DataFrame({"time": time, "metric": metric_values})


ONE_MONTH = 24 * 30
test_data = generate_data(ONE_MONTH)

# plot data

sns.lineplot(x="time", y="metric", data=test_data, alpha=0.5)
# fit linear regression
slope, intercept, rvalue, pvalue, stderr = stats.linregress(
    test_data.index, test_data.metric
)


# test for Homoscedasticity of the residuals

def test_homoscedasticity(test_data: pd.DataFrame):
    # prepare y and X for statsmodels
    y = test_data.metric
    X_sm = test_data.index
    X_sm = sm.add_constant(X_sm)
    # fit linear regression using statsmodels
    mod_sm = sm.OLS(y, X_sm)
    res_sm = mod_sm.fit()
    # test for Homoscedasticity
    bp_lm, bp_lm_pvalue, bp_fvalue, bp_f_pvalue = sm.stats.diagnostic.het_breuschpagan(
        res_sm.resid, res_sm.model.exog
    )

    return bp_f_pvalue


bp_f_pvalue = test_homoscedasticity(test_data)
if bp_f_pvalue < 0.05:
    print("Residuals are not Homoscedastic. Linear regression may be invalid")
else:
    print("Residuals are Homoscedastic")


# test autocorrelation
def test_autocorrelation(residuals):
    dw_statistic = sm.stats.durbin_watson(residuals)
    return dw_statistic


# Assuming you have the 'residuals' variable defined elsewhere in your code
# Conduct the test for autocorrelation


# Check the p-value and print the result
residuals = test_data.metric - (intercept + slope * test_data.index)
dw_statistic = test_autocorrelation(residuals)

if dw_statistic > 1.5 and 2.5:
    print("Errors are independent (no signs of autocorrelation)")
else:
    print("Аutocorrelation in the residuals. Independence of errors is not satisfactory.")

# bootstrap linear regression slope
n_iterations = 1000
slopes = np.zeros(n_iterations)
intercepts = np.zeros(n_iterations)
r_square = np.zeros(n_iterations)
for i in range(n_iterations):
    sample = test_data.sample(frac=1, replace=True)
    slope, intercept, rvalue, pvalue, stderr = stats.linregress(
        sample.index, sample.metric
    )
    slopes[i] = slope
    intercepts[i] = intercept
    r_square[i] = rvalue ** 2
# calculate 95% confidence interval
slopes_quantiles = np.percentile(slopes, [2.5, 97.5])
intercept_quantiles = np.percentile(intercepts, [2.5, 97.5])
r_square_quantiles = np.percentile(r_square, [2.5, 97.5])
print(f'r_square is between {r_square_quantiles}')
print(f'slope is between {slopes_quantiles}')

# plot two lines for the confidence interval slopes
plt.plot(
    test_data.time,
    intercept_quantiles[0] + slopes_quantiles[0] * test_data.index,
    label="lower bound",
    color="green",
    linestyle='--',
    linewidth=2
)

plt.plot(
    test_data.time,
    intercept_quantiles[1] + slopes_quantiles[1] * test_data.index,
    label="upper bound",
    color="green",
    linestyle='--',
    linewidth=2
)

# plot linear regression

plt.plot(
    test_data.time,
    intercept + slope * test_data.index,
    label="fitted line",
    color="red",
    linewidth=2
)

plt.legend()
plt.show()

За период среднее значение выросло в 2 раза trend_max = 100, но из за большого кол-ва выбросов и шума в данных R² довольно низкий. В таких случаях важно понимать что R², вопреки расхожим заблуждениям не говорит о качестве модели, а лишь о том насколько она точно аппроксимирует нашу фукнцию, и в случае с большим количеством выбросов показатель всегда будет достаточно низкий. Важно в таком случае обращать внимание на другие тесты и на slope.

Также мы не упомянули важную слабость линейной регрессии - ее неустойчивость к выбросам. Это выражается в широких доверительных интервалах для коэффициентов наклона. Для улучшения качества модели часто требуется удаление выбросов или применение более устойчивых методов, способных лучше справляться с наличием аномалий в данных. Рассмотрим некоторые из них далее.

Mann-Kendall test

Метод Манн-Кендалла представляет собой непараметрический тест, который определяет наличие или отсутствие тренда на рангах данных. Он не зависит от конкретных значений метрики, а лишь учитывает их порядок, то есть знак разности между последовательными наблюдениями.

Среди преимуществ этого теста можно отметить его независимость от выбросов и способность выявлять монотонные зависимости даже в случае нелинейной формы зависимости. Однако следует проявлять осторожность при его использовании в практических задачах, поскольку уже небольшой тренд может быть признан статистически значимым. Это может привести к неправильным выводам с точки зрения бизнес-логики, где незначительные изменения могут быть недостаточно значимыми.

Давайте рассмотрим следующий пример для наглядного представления.

Код Python

import numpy as np
import pandas as pd
import scipy.stats as stats


def generate_data(
        n_samples: int,
        loc: int = 100,
        scale: int = 0.1,
        trend_max: int = 0.5,
        hetero: float = 1,
        outlier: bool = False,
) -> pd.DataFrame:
    time = pd.date_range("2024-01-01", periods=n_samples, freq="1H")
    metric_values = (
            loc
            + np.random.normal(scale=scale, size=n_samples)
            * np.linspace(1, hetero, n_samples)
            + np.linspace(0, trend_max, n_samples)
    )
    if outlier:
        metric_values[n_samples // 2] = -np.max(metric_values) * 1000
    return pd.DataFrame({"time": time, "metric": metric_values})


ONE_MONTH = 24 * 30
test_data = generate_data(ONE_MONTH)
# test mann-kendall trend
n_iterations = 100
p_values = np.zeros(n_iterations)
taus = np.zeros(n_iterations)
for i in range(n_iterations):
    test_data = generate_data(ONE_MONTH)
    tau, pvalue = stats.mstats.kendalltau(test_data.index, test_data.metric)
    taus[i] = np.round(tau, 2)
    p_values[i] = pvalue

tau_quantiles = np.percentile(taus, [2.5, 97.5])
# number of times we reject the null hypothesis
print(f'trend is significant in {np.sum(p_values < 0.05)} cases from {n_iterations}')
print(f'tau_quantiles interval {tau_quantiles}')

Таким образом, тест Манна-Кендалла, примененный к данному примеру, покажет наличие тренда в 100 из 100 случаев, даже если этот тренд является незначительным.

Для данного примера даже бустрапированная метрика Манна-Кендалла всегда будет указывать на наличие тренда. Помимо значения p-value, тест Манна-Кендалла возвращает статистику tau, которая меняется от -1 до 1. Значение 1 указывает на положительную монотонность, -1 на отрицательную, а 0 на отсутствие монотонности. Однако эта метрика сложна для интерпретации, поскольку не предоставляет информацию о величине тренда. В текущем примере значение tau около 0.6, что указывает на положительную монотонность, но не даёт информации о конкретной величине тренда.

Таким образом, после положительного заключения о наличии тренда в тесте Ман-Кендалла необходимо дополнительно проводить оценку угла наклона другими тестами. Линейная регрессия на тех же данные показывает что slope минимальный.

Если есть возможность оценить угол наклона на втором шаге, то теряется смысл в самом тесте Ман-Кендалла. Более того, частые срабатывания теста на незначительных трендах могут привести к неправильным выводам аналитика, что поднимает вопросы о его применимости в реальных сценариях.

Несмотря на это, к достоинствам теста Манна-Кендалла следует отнести его способность быстро обнаруживать нелинейные тренды, что может быть полезно в некоторых ситуациях.

Spearman’s Rank Correlation Test

Так же, как и тест Манн-Кендалла, тест Спирмена основан на анализе рангов. В тесте Спирмена рассчитывается коэффициент корреляции Спирмена, который измеряет степень монотонной связи между двумя переменными.

Код Python

import numpy as np
import pandas as pd
import scipy.stats as stats


def generate_data(
        n_samples: int,
        loc: int = 100,
        scale: int = 0.1,
        trend_max: int = 0.1,
        hetero: float = 1,
        num_outliers: int = 50,
) -> pd.DataFrame:
    time = pd.date_range("2024-01-01", periods=n_samples, freq="1H")
    metric_values = (
            loc
            + np.random.normal(scale=scale, size=n_samples)
            * np.linspace(1, hetero, n_samples)
            + np.linspace(0, trend_max, n_samples)
    )
    if num_outliers > 0:
        indices_outliers = np.random.choice(n_samples, num_outliers, replace=False)
        metric_values[indices_outliers] = np.max(metric_values) * 100

    return pd.DataFrame({"time": time, "metric": metric_values})


ONE_MONTH = 24 * 30
test_data = generate_data(ONE_MONTH)

# test mann-kendall trend
n_iterations = 100

p_values_mk = np.zeros(n_iterations)
taus = np.zeros(n_iterations)

p_values_sp = np.zeros(n_iterations)
spearman_coefficients = np.zeros(n_iterations)

for i in range(n_iterations):
    test_data = generate_data(ONE_MONTH)
    tau, p_value_mk = stats.mstats.kendalltau(test_data.index, test_data.metric)
    spearman_coefficient, p_value_sp = stats.spearmanr(test_data.index, test_data.metric)
    taus[i] = np.round(tau, 2)
    p_values_mk[i] = p_value_mk

    spearman_coefficients[i] = round(spearman_coefficient, 2)
    p_values_sp[i] = p_value_sp

tau_quantiles = np.percentile(taus, [2.5, 97.5])
spearman_quantiles = np.percentile(spearman_coefficients, [2.5, 97.5])
# number of times we reject the null hypothesis

print(f'MK trend is significant in {np.sum(p_values_mk < 0.05)} cases from {n_iterations}')
print(f'SP trend is significant in {np.sum(p_values_sp < 0.05)} cases from {n_iterations}')

print(f'tau_quantiles interval {tau_quantiles}')
print(f'spearman_quantiles interval {spearman_quantiles}')

Недостатки и преимущества теста Ман-Кендалла также переносятся и на тест Спирмана, в том числе чувствительность к низким значениям углов тренда, что также ставит под сомнение его применимость во многих задачах где цель определить наличие или отсутствие тренда и понять его характер.

Repeated median regression (RMR)

Проблема неустойчивости линейной регрессии к выбросам привела к развитию целого семейства методов, робастных к редким большим и маленьким значениям. Идея проста: для каждой точки в выборке строятся линии, проходящая через неё и все остальные точки, и затем вычисляется угол наклона этой линии. Затем из всех полученных углов наклона выбирается медианное значение. Этот процесс повторяется для каждой точки, и затем снова вычисляется медиана всех медиан.

Когда он хорошо работает:

• С выбросами в данных: Метод робастен к наличию выбросов и аномалий, что делает его полезным в ситуациях, где классическая линейная регрессия может дать искаженные результаты из-за выбросов.

• Нелинейные зависимости: Этот метод может быть эффективен при аппроксимации нелинейных зависимостей, так как он не ограничен предположением о линейности отношения между переменными.

• Работа с разнородными данными: Метод может быть эффективен при работе с данными, содержащими шум или неоднородные структуры, так как он оценивает углы наклона, учитывая различия в данных.

Важно отметить, что, этот метод может потребовать дополнительных вычислительных ресурсов, особенно при обработке больших объемов данных.

Код Python

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.stats as stats
import seaborn as sns
import statsmodels.api as sm


def generate_data(
    n_samples: int,
    loc: int = 100,
    scale: int = 10,
    trend_max: int = 20,
    hetero: float = 10,
    n_outliers: int = 30,
) -> pd.DataFrame:
    time = pd.date_range("2024-01-01", periods=n_samples, freq="1H")
    metric_values = (
        loc
        + np.random.normal(scale=scale, size=n_samples)
        * np.linspace(1, hetero, n_samples)
        + np.linspace(0, trend_max, n_samples)
    )

    for i in range(n_outliers):
        metric_values[i * n_samples // n_outliers] = -np.max(metric_values) * 10
    return pd.DataFrame({"time": time, "metric": metric_values})


ONE_MONTH = 24 * 30
test_data = generate_data(ONE_MONTH)


# bootstrap linear regression slope
n_iterations = 100
slopes_RMR = np.zeros(n_iterations)
intercepts_RMR = np.zeros(n_iterations)
slopes_LR = np.zeros(n_iterations)
intercepts_LR = np.zeros(n_iterations)
for i in range(n_iterations):
    sample = test_data.sample(frac=1, replace=True)
    res = stats.siegelslopes(sample.metric, sample.index)
    slopes_RMR[i] = res.slope
    intercepts_RMR[i] = res.intercept
    res_reg = stats.linregress(sample.metric, sample.index)
    slopes_LR[i] = res_reg.slope
    intercepts_LR[i] = res_reg.intercept
# calculate 95% confidence interval
slopes_quantiles_RMR = np.percentile(slopes_RMR, [2.5, 97.5])
intercept_quantiles_RMR = np.percentile(intercepts_RMR, [2.5, 97.5])

slopes_quantiles_LR = np.percentile(slopes_LR, [2.5, 97.5])
intercept_quantiles_LR = np.percentile(intercepts_LR, [2.5, 97.5])

# results RMR and LR
print(f"RMR  95% confidence interval for slope: {slopes_quantiles_RMR}")
print(f"RMR  95% confidence interval for intercept: {intercept_quantiles_RMR}")

print(f"LR 95% confidence interval for slope: {slopes_quantiles_LR}")
print(f"LR 95% confidence interval for intercept: {intercept_quantiles_LR}")

Исходя из результатов, видно, что этот метод не только устойчив к выбросам, но и эффективно работает в связке с бутстрапом. Другим преимуществом является то, что этот метод менее чувствителен к нарушениям гомоскедастичности и может давать хорошие результаты, даже если дисперсия ошибок не является постоянной на всех уровнях предиктора. Однако важно понимать, что он менее чувствителен, и при высоких значениях дисперсии качество модели также ухудшается. Ниже приведен результат выполнения при hetero=30.

Как видно RMR выдает значение приближенные к реальным, в то время как у классической линейной регрессии интервал становится крайне широким.

Заключение

Мы рассмотрели несколько популярных методов для определения тренда. Каждый из них имеет свои сильные и слабые стороны и выбор конкретного метода зависит от целей анализа, характера данных и требуемой точности оценки тренда.
В данной статье мы практически не затронули методы которые используется при анализе тренда в нелинейных взаимосвязях таких как полиномиальная регрессия, экспоненциальная регрессия, логистическая регрессия и другие. Или более сложные алгоритмы такие как ARIMA, SARIMA и другие.

Дополнительные проверки на автокоррелляцию и гетероскедастичность можно найти https://www.statsmodels.org/stable/diagnostic.html.

О том чем это чревата неосторожная очистка выбросов можно почитать в статье https://habr.com/ru/articles/781060/

Авторы

Roman Rudnitskiy

Marat Yuldashev

Mikhail Tretyakov