Трехвариантный корреляционный анализ: классификация одной переменной по интервалам для достижения лучшей корреляции между двумя другими

Question

После этого поста и моей неспособности понять метод вычисления интервала Джонсона-Неймана [JN], я думаю о том, чтобы выйти из коробки и рассмотреть другие алгоритмы, чтобы сделать тривариант (многовариантный между тремя переменными) ) корреляционный анализ. Рассмотрим следующий Pandas фрейм данных:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import pearsonr, spearmanr, kendalltau, pointbiserialr

np.random.seed(42)

x1 = np.random.rand(10) * 10
x2 = np.random.rand(10) * 10
x3 = np.random.rand(10) * 10

x = np.concatenate((x1, x2, x3))

y1 = 1.5 * x1 + (np.random.rand(10) - 0.5) * 3.5
y2 = 0.5 * x2 + (np.random.rand(10) - 0.5) * 4.1
y3 = 3.0 * x3 + (np.random.rand(10) - 0.5) * 2.4

y = np.concatenate((y1, y2, y3))

M = [1.5, 5.5]

m1 = np.random.rand(10) * M[0]
m2 = np.random.rand(10) * (M[1] - M[0]) + M[0]
m3 = np.random.rand(10) * 2 + M[1]

m = np.concatenate((m1, m2, m3))

df = pd.DataFrame({'X': x, 'Y': y, 'M': m})
df = df.sample(frac=1).reset_index(drop=True)

Теперь мне нужно иметь функцию:

classify_moderator(dataframe=df, predictor="x", response="y", moderator="m")

, которая возвращает список из одного или нескольких значений для m, который разбивает фрейм данных на две или более значимые корреляционные группы. Например, для приведенного выше фрейма данных мы должны получить M = [1.5, 5.5], в котором в группах df["m"] < M[0], M[0] < df["m"] < M[1] и M[1] < df["m"] ассоциация между подгруппами df["x"] и df["y"] сильнее суммарных данных фрейм с точки зрения двухфакторных алгоритмов корреляционного анализа, таких как методы Пирсона, Кендалла, Спирмена и Пойнт-Бизерана. Важно учитывать, что мы не знаем точно, сколько существует интересующих разделов. Это может быть два или более. Но если это делает метод излишне сложным, мы можем принять и число.

Я предпочитаю обычные методы (например, модуль scipy.optimize), а не heuristi c, поскольку они помогают мне лучше понять основную математику, но в остальном также приветствуются. Спасибо за вашу поддержку заранее.

PS1. Здесь моя первая неудачная попытка сгруппироваться в две части:

def groupedCor(modValue, dataframe, predictor, response, moderator):

    df1 = dataframe.loc[dataframe[moderator] < modValue]
    df2 = dataframe.loc[modValue <= dataframe[moderator]]

    return pearsonr(df1[predictor], df1[response])[1] + pearsonr(df2[predictor], df2[response])[1]

M1 = (df["M"].max() - df["M"].min()) / 2 + df["M"].min()

M1 = fmin(groupedCor, x0=M1, args=(df, "X", "Y", "M"))

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

ValueError: длины должны совпадать, чтобы сравнивать

, тогда я, вероятно, мог бы решить проблему самостоятельно.

PS2. Я понял, что ни один из методов двумерного корреляционного анализа не является хорошим представителем. Поэтому я использовал r-квадрат линейной регрессии:

from sklearn.linear_model import LinearRegression

def groupedScore(modValue, dataframe, predictor, response, moderator):
    df1 = dataframe[dataframe[moderator] < modValue]
    df2 = dataframe[modValue <= dataframe[moderator]]

    x1 = df1[predictor].values
    x1_ = x1.reshape((-1, 1))
    y1 = df1[response].values
    model1 = LinearRegression().fit(x1_, y1)

    x2 = df2[predictor].values
    x2_ = x2.reshape((-1, 1))
    y2 = df2[response].values
    model2 = LinearRegression().fit(x2_, y2)

    return (model1.score(x1_, y1) + model2.score(x2_, y2)) / 2

def groupedScore_v(x, dataframe, predictor, response, moderator):
    return [[groupedScore(xi, dataframe, predictor, response, moderator) for xi in x[2:-1]], x[2:-1]]

scores = groupedScore_v(np.sort(df["M"].values).tolist(), df, "X", "Y", "M")

plt.plot(np.asarray(scores[1]), np.asarray(scores[0]))

, что приводит к:

              <img src="https://i.stack.imgur.com/U7qFq.png" width="400">               

, что является очень странным результатом!