另一种方法是brglm包。例如，使用与@chl's Answer 相同的数据/模型

data(sex2)
fm <- case ~ age + oc + vic + vicl + vis + dia
fit <- brglm(fm, data = sex2)

predict(fit, newdata = sex2[1:5, ], type = "response")

这会产生

> predict(fit, newdata = sex2[1:5, ], type = "response")
        1         2         3         4         5 
0.3389307 0.9159945 0.9159945 0.9159945 0.9159945

请注意，在这种brglm()情况下，由于函数的工作方式，您在上面看到的只是predict.glm()R 中标准函数/方法的结果。

您可能可以用很少的代数计算任何您想要的预测。让我们考虑示例数据集，

data(sex2)
fm <- case ~ age+oc+vic+vicl+vis+dia
fit <- logistf(fm, data=sex2)

一旦我们得到回归系数，设计矩阵是计算预测概率的唯一缺失部分，由下式给出

betas <- coef(fit)

因此，让我们首先尝试对观察到的数据进行预测：

X <- model.matrix(fm, data=sex2)       # add a column of 1's to sex2[,-1]
pi.obs <- 1 / (1 + exp(-X %*% betas))  # in case there's an offset, δ, it 
                                       # should be subtracted as exp(-Xβ - δ)

我们可以检查我们是否得到了正确的结果

> pi.obs[1:5]
[1] 0.3389307 0.9159945 0.9159945 0.9159945 0.9159945
> fit$predict[1:5]
[1] 0.3389307 0.9159945 0.9159945 0.9159945 0.9159945

现在，您可以在上面的设计矩阵中放入X您感兴趣的值。例如，将所有协变量设置为 1

new.x <- c(1, rep(1, 6))
1 / (1 + exp(-new.x %*% betas)) 

我们得到个体概率为 0.804，而当所有协变量都设置为 0 ( new.x <- c(1, rep(0, 6))) 时，估计概率为 0.530。