从 R 中的logistf中提取数据答案

【问题标题】：Extracting data from logistf in R从 R 中的logistf中提取数据
【发布时间】：2014-11-23 09:07:29
【问题描述】：

我不知道如何从logistf() 回归模型中提取标准错误“sd(coef)”信息。这些模型属于 logistf 类，手册指出可以通过这种方式提取数据：

logistf 的输出对象可以使用以下通用方法：print、summary、coef、vcov、confint、anova、extractAIC、add1、drop1、profile、terms、nobs。

但是，标准错误不存在。在 str(summary(fit)) 中没有用于 se(coef) 的对象，我查看了源代码没有运气。

任何帮助将不胜感激！

logistf(formula = newdata[, i] ~ newdata[, j], data = newdata, 
    firth = TRUE)
Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood 

                      coef  se(coef) lower 0.95 upper 0.95        Chisq p
(Intercept)  -4.110874e+00 0.8276236  -6.283919  -2.824075          Inf 0
newdata[, j]  1.691332e-08 1.6689839  -4.993849   2.957865 3.552714e-15 1

Likelihood ratio test=3.552714e-15 on 1 df, p=1, n=122

【问题讨论】：

标签： r logistic-regression logistf

【解决方案1】：

确实sd(fit) 不起作用，我不确定这通常是否适用于其他模型类。

但是，假设logistf 模型对象位于fit 中，则协方差矩阵可通过vcov(fit) 获得。那么只需计算协方差矩阵的主对角线的平方根就可以得到se(coef)列：sqrt(diag(vcov(fit)))

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

> fit
logistf(formula = case ~ age + oc + vic + vicl + vis + dia, data = sex2)
Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood 

                   coef  se(coef) lower 0.95  upper 0.95       Chisq            p
(Intercept)  0.12025404 0.4855415 -0.8185574  1.07315110  0.06286298 8.020268e-01
age         -1.10598130 0.4236601 -1.9737884 -0.30742658  7.50773092 6.143472e-03
oc          -0.06881673 0.4437934 -0.9414360  0.78920202  0.02467044 8.751911e-01
vic          2.26887464 0.5484160  1.2730233  3.43543174 22.93139022 1.678877e-06
vicl        -2.11140816 0.5430824 -3.2608608 -1.11773667 19.10407252 1.237805e-05
vis         -0.78831694 0.4173676 -1.6080879  0.01518319  3.69740975 5.449701e-02
dia          3.09601263 1.6750089  0.7745730  8.03029352  7.89693139 4.951873e-03

Likelihood ratio test=49.09064 on 6 df, p=7.15089e-09, n=239

> diag(vcov(fit))
[1] 0.2357506 0.1794879 0.1969526 0.3007601 0.2949384 0.1741957 2.8056550

> sqrt(diag(vcov(fit)))
[1] 0.4855415 0.4236601 0.4437934 0.5484160 0.5430824 0.4173676 1.6750089

【讨论】：