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

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 回归中的默认方法更改为方法 ='bfgs'.

解决了我的问题

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

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

Is there something similar to R's brglm to help deal with quasi-separation in Python using statsmodels Logit?

python

statsmodels