有没有类似于R的brglm来帮助使用statsmodels Logit来处理Python中的准分离？

Question

我使用statsmodels的Logit来创建回归模型。

我得到错误：LinAlgError：奇异矩阵然后当我从我的数据集中一次删除1个变量时，我终于得到了一个不同的错误：PerfectSeparationError：检测到完美分离，结果不可用。

我怀疑原始错误（LinAlgError）与完美分离有关，因为我在R中遇到了同样的问题，并使用brglm（偏差减少的glm）绕过它。

我有一个布尔y变量和23个数字和布尔x变量。

我已经运行了一个VIF函数来删除任何具有高多重共线性分数的变量（我从26个变量开始）。

我尝试使用firth_regression.py代替完美分离，但我收到了一个内存错误：MemoryError。（https://gist.github.com/johnlees/3e06380965f367e4894ea20fbae2b90d）

我已经尝试过sklearn中的LogisticRegression，但无法得到对我不利的p值。

我甚至尝试从我的数据集中一次删除1个变量。当我得到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')