R-Squared with logistf

[Jeffcmoore]

I'm learning R after years using SPSS. One of the reasons for the transition is access to the firth method via logistf.

I'm able to run analysis- but cannot find how to compute Pseudo R sqaured.

I also get error messages when I tried to run Loglik (saying this function ie not availablenworh logisft). In the absence of computer R sqaured I'm thinking that I could calculate McFadden R squared if I had log likelihood values.

Also - if it is possible to access log likelihood values is the somewhere I can computer the r pseudo r sqaured - I dont trust myself to do this manually.

[technocrat]

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.

[Jeffcmoore]

@technocrat - richard - thanks for the v helpful response. Can you tell me why R provides two figures for loglik?

[technocrat]

From help(logitf)

loglik
a vector of the (penalized) log-likelihood of the restricted and the full models.

[Jeffcmoore]

Thanks you - v much appreicated.