Inspection of a fitted model object shows all the content
library(logistf)
data(sex2)
fit<-logistf(case ~ age+oc+vic+vicl+vis+dia, data=sex2)
# summary(fit)
# nobs(fit)
# drop1(fit)
# plot(profile(fit,variable="dia"))
# extractAIC(fit)
str(fit)
#> List of 28
#> $ coefficients : Named num [1:7] 0.1203 -1.106 -0.0688 2.2689 -2.1114 ...
#> ..- attr(*, "names")= chr [1:7] "(Intercept)" "age" "oc" "vic" ...
#> $ alpha : num 0.05
#> $ terms : chr [1:7] "(Intercept)" "age" "oc" "vic" ...
#> $ var : num [1:7, 1:7] 0.2358 -0.0274 -0.1952 -0.1211 -0.0289 ...
#> $ df : num 6
#> $ loglik : num [1:2] -157 -133
#> $ iter : int 9
#> $ n : num 239
#> $ y : Named int [1:239] 1 1 1 1 1 1 1 1 1 1 ...
#> ..- attr(*, "names")= chr [1:239] "1" "2" "3" "4" ...
#> $ formula :Class 'formula' language case ~ age + oc + vic + vicl + vis + dia
#> .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
#> $ call : language logistf(formula = case ~ age + oc + vic + vicl + vis + dia, data = sex2)
#> $ conv : Named num [1:3] 3.44e-12 4.54e-07 4.86e-06
#> ..- attr(*, "names")= chr [1:3] "LL change" "max abs score" "beta change"
#> $ firth : logi TRUE
#> $ linear.predictors: num [1:239] -0.668 2.389 2.389 2.389 2.389 ...
#> $ predict : num [1:239] 0.339 0.916 0.916 0.916 0.916 ...
#> $ hat.diag : num [1:239] 0.071 0.0226 0.0226 0.0226 0.0226 ...
#> $ method : chr "Penalized ML"
#> $ conflev : num 0.95
#> $ ci.lower : Named num [1:7] -0.819 -1.974 -0.941 1.273 -3.261 ...
#> ..- attr(*, "names")= chr [1:7] "(Intercept)" "age" "oc" "vic" ...
#> $ ci.upper : Named num [1:7] 1.073 -0.307 0.789 3.435 -1.118 ...
#> ..- attr(*, "names")= chr [1:7] "(Intercept)" "age" "oc" "vic" ...
#> $ prob : Named num [1:7] 8.02e-01 6.14e-03 8.75e-01 1.68e-06 1.24e-05 ...
#> ..- attr(*, "names")= chr [1:7] "(Intercept)" "age" "oc" "vic" ...
#> $ method.ci : chr [1:7] "Profile Likelihood" "Profile Likelihood" "Profile Likelihood" "Profile Likelihood" ...
#> $ pl.iter : num [1:7, 1:2] 19 9 13 19 10 8 18 29 19 19 ...
#> $ betahist :List of 2
#> ..$ lower:List of 7
#> .. ..$ : num [1:19, 1:7] -0.356 -0.588 -0.703 -0.761 -0.79 ...
#> .. ..$ : num [1:9, 1:7] 0.247 0.24 0.236 0.234 0.233 ...
#> .. ..$ : num [1:13, 1:7] 0.551 0.985 0.986 0.986 0.986 ...
#> .. ..$ : num [1:19, 1:7] 0.337 0.451 0.51 0.54 0.554 ...
#> .. ..$ : num [1:10, 1:7] 0.225 0.222 0.222 0.222 0.222 ...
#> .. ..$ : num [1:8, 1:7] 0.339 0.352 0.352 0.352 0.352 ...
#> .. ..$ : num [1:18, 1:7] 0.136 0.155 0.166 0.172 0.176 ...
#> ..$ upper:List of 7
#> .. ..$ : num [1:29, 1:7] 0.596 1.073 1.073 1.073 1.073 ...
#> .. ..$ : num [1:19, 1:7] 0.056812 0.020358 0.000966 -0.009094 -0.014263 ...
#> .. ..$ : num [1:19, 1:7] -0.311 -0.52 -0.625 -0.676 -0.702 ...
#> .. ..$ : num [1:10, 1:7] -0.313 -0.269 -0.268 -0.267 -0.267 ...
#> .. ..$ : num [1:19, 1:7] 0.068 0.0439 0.0322 0.0264 0.0235 ...
#> .. ..$ : num [1:19, 1:7] 0.0111 -0.0398 -0.065 -0.0776 -0.084 ...
#> .. ..$ : num [1:14, 1:7] 0.0887 0.11 0.1099 0.1098 0.1098 ...
#> $ pl.conv : num [1:7, 1:4] 7.15e-06 3.99e-10 1.63e-06 6.99e-06 1.04e-06 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : chr [1:7] "1" "2" "3" "4" ...
#> .. ..$ : chr [1:4] "lower, loglik" "lower, beta" "upper, loglik" "upper, beta"
#> $ flic : logi FALSE
#> $ control :List of 7
#> ..$ maxit : num 25
#> ..$ maxhs : num 5
#> ..$ maxstep : num 5
#> ..$ lconv : num 1e-05
#> ..$ gconv : num 1e-05
#> ..$ xconv : num 1e-05
#> ..$ collapse: logi TRUE
#> $ model :'data.frame': 239 obs. of 7 variables:
#> ..$ case: int [1:239] 1 1 1 1 1 1 1 1 1 1 ...
#> ..$ age : int [1:239] 0 0 0 0 0 0 0 0 0 0 ...
#> ..$ oc : int [1:239] 0 0 0 0 0 0 0 0 0 0 ...
#> ..$ vic : int [1:239] 0 1 1 1 1 1 1 1 1 1 ...
#> ..$ vicl: int [1:239] 0 0 0 0 0 0 0 0 0 0 ...
#> ..$ vis : int [1:239] 1 0 0 0 0 0 0 0 0 0 ...
#> ..$ dia : int [1:239] 0 0 0 0 0 0 0 0 0 0 ...
#> ..- attr(*, "terms")=Classes 'terms', 'formula' language case ~ age + oc + vic + vicl + vis + dia
#> .. .. ..- attr(*, "variables")= language list(case, age, oc, vic, vicl, vis, dia)
#> .. .. ..- attr(*, "factors")= int [1:7, 1:6] 0 1 0 0 0 0 0 0 0 1 ...
#> .. .. .. ..- attr(*, "dimnames")=List of 2
#> .. .. .. .. ..$ : chr [1:7] "case" "age" "oc" "vic" ...
#> .. .. .. .. ..$ : chr [1:6] "age" "oc" "vic" "vicl" ...
#> .. .. ..- attr(*, "term.labels")= chr [1:6] "age" "oc" "vic" "vicl" ...
#> .. .. ..- attr(*, "order")= int [1:6] 1 1 1 1 1 1
#> .. .. ..- attr(*, "intercept")= int 1
#> .. .. ..- attr(*, "response")= int 1
#> .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
#> .. .. ..- attr(*, "predvars")= language list(case, age, oc, vic, vicl, vis, dia)
#> .. .. ..- attr(*, "dataClasses")= Named chr [1:7] "numeric" "numeric" "numeric" "numeric" ...
#> .. .. .. ..- attr(*, "names")= chr [1:7] "case" "age" "oc" "vic" ...
#> - attr(*, "class")= chr "logistf"
Created on 2020-09-18 by the reprex package (v0.3.0)
So that to obtain log likelihood values
> fit$loglik
[1] -157.0847 -132.5394
A S/O thread has useful answers for calculating flavors of R^2.
The best way to proceed is to identify the pieces and and implement the algorithm as a function if doing it oneself.