5

I cannot figure out how to extract the standard error "sd(coef)" information from the logistf() regression model. These models are of class logistf and the manual states that data can be extracted this way:

The following generic methods are available for logistf`s output object: print, summary, coef, vcov, confint, anova, extractAIC, add1, drop1, profile, terms, nobs.

However, the standard error is not there.There isn't an object in str(summary(fit)) for the se(coef) and I have had a look at the source code without luck.

Any help would be appreciated!

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
MichaelChirico
  • 33,841
  • 14
  • 113
  • 198
Lel
  • 139
  • 1
  • 2
  • 13

1 Answers1

6

Indeed sd(fit) does not work and I am not sure if that usually works for other model classes.

However, the covariance matrix is available via vcov(fit) assuming the logistf model object is in fit. Then you can get the se(coef) column simply by calculating the square root of the main diagonal of the covariance matrix: 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
thie1e
  • 3,588
  • 22
  • 22