Есть ли что-то похожее на brglm R, чтобы помочь справиться с квазоразделением в Python, используя statsmodels Logit?

Question

Я использую Logit из statsmodels для создания регрессионной модели.

Я получаю ошибку: LinAlgError: Сингулярная матрица, а затем, когда я удаляю 1 переменную за раз из своего набора данных, я, наконец, получаю другую ошибку: PerfectSeparationError: Обнаружено идеальное разделение, результаты недоступны.

Я подозреваю, что исходная ошибка (LinAlgError) связана с идеальным разделением, потому что у меня была та же проблема в R, и я обошел ее, используя brglm (смещение, уменьшенное glm).

У меня есть логическая переменная y и 23 числовые и логические переменные x.

Я уже запустил функцию VIF, чтобы удалить все переменные, которые имеют высокие оценки мультиколлинеарности (я начал с 26 переменных).

Я попытался использовать firth_regression.py вместо этого, чтобы учесть идеальное разделение, но я получил ошибку памяти: MemoryError. (https://gist.github.com/johnlees/3e06380965f367e4894ea20fbae2b90d)

Я пробовал LogisticRegression от sklearn, но не могу получить значения p, что мне не подходит.

Я даже пытался удалить по одной переменной за раз из моего набора данных. Когда у меня осталось 4 переменные (у меня было 23), я получил PerfectSeparationError: Обнаружено идеальное разделение, результаты недоступны.

Кто-нибудь сталкивался с этим и как с этим справиться?

Ценю любой совет!

    X = df.loc[:, df.columns != 'VehicleMake']
    y = df.iloc[:,0]
    # Split data
    X_train, X_test, y_train, y_test = skl.model_selection.train_test_split(X, y, test_size=0.3)

Код вопроса:

    # Perform logistic regression and get p values
    logit_model = sm.Logit(y_train, X_train.astype(float))
    result = logit_model.fit()

Я попытался использовать это firth_regression, которое привело к ошибке памяти:

# For the firth_regression
import sys
import warnings
import math
import statsmodels
from scipy import stats
import statsmodels.formula.api as smf


def firth_likelihood(beta, logit):
    return -(logit.loglike(beta) + 0.5*np.log(np.linalg.det(-logit.hessian(beta))))

step_limit=1000
convergence_limit=0.0001

logit_model = smf.Logit(y_train, X_train.astype(float))

start_vec = np.zeros(X.shape[1])

beta_iterations = []
beta_iterations.append(start_vec)
for i in range(0, step_limit):
    pi = logit_model.predict(beta_iterations[i])
    W = np.diagflat(np.multiply(pi, 1-pi))
    var_covar_mat = np.linalg.pinv(-logit_model.hessian(beta_iterations[i]))

    # build hat matrix
    rootW = np.sqrt(W)
    H = np.dot(np.transpose(X_train), np.transpose(rootW))
    H = np.matmul(var_covar_mat, H)
    H = np.matmul(np.dot(rootW, X), H)

    # penalised score
    U = np.matmul(np.transpose(X_train), y - pi + np.multiply(np.diagonal(H), 0.5 - pi))
    new_beta = beta_iterations[i] + np.matmul(var_covar_mat, U)

    # step halving
    j = 0
    while firth_likelihood(new_beta, logit_model) > firth_likelihood(beta_iterations[i], logit_model):
        new_beta = beta_iterations[i] + 0.5*(new_beta - beta_iterations[i])
        j = j + 1
        if (j > step_limit):
            sys.stderr.write('Firth regression failed\n')
            None

    beta_iterations.append(new_beta)
    if i > 0 and (np.linalg.norm(beta_iterations[i] - beta_iterations[i-1]) < convergence_limit):
        break

return_fit = None
if np.linalg.norm(beta_iterations[i] - beta_iterations[i-1]) >= convergence_limit:
    sys.stderr.write('Firth regression failed\n')
else:
# Calculate stats
    fitll = -firth_likelihood(beta_iterations[-1], logit_model)
    intercept = beta_iterations[-1][0]
    beta = beta_iterations[-1][1:].tolist()
    bse = np.sqrt(np.diagonal(-logit_model.hessian(beta_iterations[-1])))

    return_fit = intercept, beta, bse, fitll
#print(return_fit)

Answer 1

Я исправил свою проблему, изменив метод по умолчанию в logit-регрессии на method = 'bfgs'. ​​

result = logit_model.fit(method = 'bfgs')