Feature enhancement: exporting to latex logitr results

[Ales-G]

Dear Helvy, thanks a lot to you for your great work and package! it is really helpful!

This may be a trivial question, but is there a way to easily extract coefficients and standard errors from logitr models and export them into latex? For instance, using logitr results to create tables with stargazer? It would be helpful to compare multiple models!

thanks a lot in advance for your help and apologies if this feature is already available!

[jhelvy]

This is a great suggestion! I do eventually plan on getting {logitr} integrated into other packages like {broom} and {stargazer}, but it's further back on the development list. For now, here's a solution that hopefully gets you close to what you need. After I estimate a model, I get the model summary table as a data frame using coef(summary(model)), then I can manually add the stars. Finally, {xtable} is a handy package for converting a data frame into LaTeX. You could also check out {kable} and {kableExtra} for getting more customizing LaTeX tables straight from the data frame.

library(logitr)
library(dplyr)
library(tibble)        
library(xtable)

# Estimate model
model <- logitr(
  data    = yogurt,
  outcome = "choice",
  obsID   = "obsID",
  pars    = c("price", "feat", "brand")
)

# Get summary table
summary_table <- coef(summary(model)) %>%
    round(3) %>%
    rownames_to_column() %>% 
    rename(
      "coefficients" = "rowname",
      "prob" = "Pr(>|z|)") %>%
    mutate(
        sig = ifelse(
            prob <= 0.001,'***', ifelse(
            prob > 0.001 & prob <= 0.01, '**', ifelse(
            prob > 0.01 & prob <= 0.05, '*', ifelse(
            prob > 0.05 & prob <= 0.1, '.', '   '))))
    )

# Convert table into LaTeX
xtable(summary_table)

[Ales-G]

Dear Helvy, thanks a lot again for your quick, and extremely helpful reply. I will use it for sure in my work!! You are really doing great work!

[etiennebacher]

Hi @jhelvy, I don't use this package but I came across this issue by chance. I just wanted to let you know that you might consider adding logitr to the list of packages supported by texreg. It is a great package to make tables in latex (texreg()), html (htmlreg()) or to print in the console (screenreg()) (there are also other *reg() functions).

This package is still active and it's super easy to add new classes. You just need to implement extract.logitr() and everything else is taken care of automatically. I did a small try with the following code:

extract.logitr <- function(model,
                           include.nobs = TRUE,
                           include.loglik = TRUE,
                           include.aic = TRUE,
                           include.bic = TRUE,
                           ...) {
  summary <- summary(model)

  coefnames <- names(summary$coefficients)
  coefs <- summary$coefTable$Estimate
  se <- summary$coefTable$`Std. Error`
  pval <- summary$coefTable$`Pr(>|z|)`

  n <- summary$n$obs
  ll <- summary$logLik

  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.nobs == TRUE) {
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Num. obs.")
    gof.decimal <- c(gof.decimal, FALSE)
  }
  if (include.loglik == TRUE) {
    gof <- c(gof, ll)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.aic == TRUE) {
    gof <- c(gof, summary$statTable["AIC", ])
    gof.names <- c(gof.names, "AIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.bic == TRUE) {
   gof <- c(gof, summary$statTable["BIC", ])
    gof.names <- c(gof.names, "BIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }

  tr <- createTexreg(
    coef.names = coefnames,
    coef = coefs,
    se = se,
    pvalues = pval,
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}

I tested with some examples on your website and it seems to work well (but I'm not sure since I don't use logitr):

library(texreg)
#> Version:  1.38.6
#> Date:     2022-04-06
#> Author:   Philip Leifeld (University of Essex)
#> 
#> Consider submitting praise using the praise or praise_interactive functions.
#> Please cite the JSS article in your publications -- see citation("texreg").
library(logitr)
#> Warning: package 'logitr' was built under R version 4.2.1

mnl_pref <- logitr(
  data    = yogurt,
  outcome = "choice",
  obsID   = "obsID",
  pars    = c("price", "feat", "brand")
)
#> Running model...
#> Done!

screenreg(mnl_pref, stars = c(0.01, 0.05, 0.1))
#> 
#> ============================
#>                 Model 1     
#> ----------------------------
#> price              -0.37 ***
#>                    (0.02)   
#> feat                0.49 ***
#>                    (0.12)   
#> brandhiland        -3.72 ***
#>                    (0.15)   
#> brandweight        -0.64 ***
#>                    (0.05)   
#> brandyoplait        0.73 ***
#>                    (0.08)   
#> ----------------------------
#> Num. obs.        2412       
#> Log Likelihood  -2656.89    
#> AIC              5323.78    
#> BIC              5352.72    
#> ============================
#> *** p < 0.01; ** p < 0.05; * p < 0.1
texreg(mnl_pref, stars = c(0.01, 0.05, 0.1))
#> 
#> \begin{table}
#> \begin{center}
#> \begin{tabular}{l c}
#> \hline
#>  & Model 1 \\
#> \hline
#> price          & $-0.37^{***}$ \\
#>                & $(0.02)$      \\
#> feat           & $0.49^{***}$  \\
#>                & $(0.12)$      \\
#> brandhiland    & $-3.72^{***}$ \\
#>                & $(0.15)$      \\
#> brandweight    & $-0.64^{***}$ \\
#>                & $(0.05)$      \\
#> brandyoplait   & $0.73^{***}$  \\
#>                & $(0.08)$      \\
#> \hline
#> Num. obs.      & $2412$        \\
#> Log Likelihood & $-2656.89$    \\
#> AIC            & $5323.78$     \\
#> BIC            & $5352.72$     \\
#> \hline
#> \multicolumn{2}{l}{\scriptsize{$^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.1$}}
#> \end{tabular}
#> \caption{Statistical models}
#> \label{table:coefficients}
#> \end{center}
#> \end{table}

mnl_wtp_unweighted <- logitr(
  data    = cars_us,
  outcome = 'choice',
  obsID   = 'obsnum',
  pars    = c(
    'hev', 'phev10', 'phev20', 'phev40', 'bev75', 'bev100', 'bev150',
    'american', 'japanese', 'chinese', 'skorean', 'phevFastcharge',
    'bevFastcharge','opCost', 'accelTime'),
  scalePar   = 'price',
  robust     = TRUE,
  # Since WTP space models are non-convex, run a multistart
  numMultiStarts =10
)
#> Setting clusterID to 'obsnum' since robust == TRUE
#> Running multistart...
#>   Iterations: 10
#>   Cores: 7
#> Done!

mnl_wtp_weighted <- logitr(
  data    = cars_us,
  outcome = 'choice',
  obsID   = 'obsnum',
  pars    = c(
    'hev', 'phev10', 'phev20', 'phev40', 'bev75', 'bev100', 'bev150',
    'american', 'japanese', 'chinese', 'skorean', 'phevFastcharge',
    'bevFastcharge','opCost', 'accelTime'),
  scalePar = 'price',
  weights  = 'weights', # This enables the weights
  robust   = TRUE,
  numMultiStarts =10
)
#> Setting clusterID to 'obsnum' since robust == TRUE
#> Running multistart...
#>   Iterations: 10
#>   Cores: 7
#> Done!

screenreg(
  list(
    mnl_wtp_unweighted, 
    mnl_wtp_weighted
  ), 
  stars = c(0.01, 0.05, 0.1),
  custom.model.names = c("Unweighted", "Weighted")
)
#> 
#> ==========================================
#>                 Unweighted    Weighted    
#> ------------------------------------------
#> scalePar            0.07 ***      0.05 ***
#>                    (0.00)        (0.00)   
#> hev                 0.81         -1.18    
#>                    (1.00)        (2.91)   
#> phev10              1.16          0.03    
#>                    (1.06)        (3.13)   
#> phev20              1.65          1.69    
#>                    (1.06)        (3.10)   
#> phev40              2.58 **       2.65    
#>                    (1.05)        (2.99)   
#> bev75             -16.05 ***    -20.14 ***
#>                    (1.25)        (3.67)   
#> bev100            -13.00 ***    -19.50 ***
#>                    (1.24)        (3.63)   
#> bev150             -9.57 ***    -13.69 ***
#>                    (1.16)        (3.49)   
#> american            2.34 ***      8.19 ***
#>                    (0.80)        (2.41)   
#> japanese           -0.38          0.93    
#>                    (0.80)        (2.36)   
#> chinese           -10.27 ***    -19.01 ***
#>                    (0.89)        (2.85)   
#> skorean            -6.03 ***     -9.51 ***
#>                    (0.85)        (2.52)   
#> phevFastcharge      2.88 ***      3.94 *  
#>                    (0.80)        (2.36)   
#> bevFastcharge       2.92 ***      3.34    
#>                    (0.92)        (2.81)   
#> opCost             -1.64 ***     -1.60 ***
#>                    (0.07)        (0.19)   
#> accelTime          -1.70 ***     -1.17 ** 
#>                    (0.16)        (0.48)   
#> ------------------------------------------
#> Num. obs.        5760          5760       
#> Log Likelihood  -4616.95      -3425.63    
#> AIC              9265.90       6883.26    
#> BIC              9372.44       6989.80    
#> ==========================================
#> *** p < 0.01; ** p < 0.05; * p < 0.1
texreg(
  list(
    mnl_wtp_unweighted, 
    mnl_wtp_weighted
  ), 
  stars = c(0.01, 0.05, 0.1),
  custom.model.names = c("Unweighted", "Weighted")
)
#> 
#> \begin{table}
#> \begin{center}
#> \begin{tabular}{l c c}
#> \hline
#>  & Unweighted & Weighted \\
#> \hline
#> scalePar       & $0.07^{***}$   & $0.05^{***}$   \\
#>                & $(0.00)$       & $(0.00)$       \\
#> hev            & $0.81$         & $-1.18$        \\
#>                & $(1.00)$       & $(2.91)$       \\
#> phev10         & $1.16$         & $0.03$         \\
#>                & $(1.06)$       & $(3.13)$       \\
#> phev20         & $1.65$         & $1.69$         \\
#>                & $(1.06)$       & $(3.10)$       \\
#> phev40         & $2.58^{**}$    & $2.65$         \\
#>                & $(1.05)$       & $(2.99)$       \\
#> bev75          & $-16.05^{***}$ & $-20.14^{***}$ \\
#>                & $(1.25)$       & $(3.67)$       \\
#> bev100         & $-13.00^{***}$ & $-19.50^{***}$ \\
#>                & $(1.24)$       & $(3.63)$       \\
#> bev150         & $-9.57^{***}$  & $-13.69^{***}$ \\
#>                & $(1.16)$       & $(3.49)$       \\
#> american       & $2.34^{***}$   & $8.19^{***}$   \\
#>                & $(0.80)$       & $(2.41)$       \\
#> japanese       & $-0.38$        & $0.93$         \\
#>                & $(0.80)$       & $(2.36)$       \\
#> chinese        & $-10.27^{***}$ & $-19.01^{***}$ \\
#>                & $(0.89)$       & $(2.85)$       \\
#> skorean        & $-6.03^{***}$  & $-9.51^{***}$  \\
#>                & $(0.85)$       & $(2.52)$       \\
#> phevFastcharge & $2.88^{***}$   & $3.94^{*}$     \\
#>                & $(0.80)$       & $(2.36)$       \\
#> bevFastcharge  & $2.92^{***}$   & $3.34$         \\
#>                & $(0.92)$       & $(2.81)$       \\
#> opCost         & $-1.64^{***}$  & $-1.60^{***}$  \\
#>                & $(0.07)$       & $(0.19)$       \\
#> accelTime      & $-1.70^{***}$  & $-1.17^{**}$   \\
#>                & $(0.16)$       & $(0.48)$       \\
#> \hline
#> Num. obs.      & $5760$         & $5760$         \\
#> Log Likelihood & $-4616.95$     & $-3425.63$     \\
#> AIC            & $9265.90$      & $6883.26$      \\
#> BIC            & $9372.44$      & $6989.80$      \\
#> \hline
#> \multicolumn{3}{l}{\scriptsize{$^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.1$}}
#> \end{tabular}
#> \caption{Statistical models}
#> \label{table:coefficients}
#> \end{center}
#> \end{table}

^{Created on 2022-07-14 by the reprex package (v2.0.1)}

I hope it's useful!

[jhelvy]

Very nice @etiennebacher! And well-timed too. I was getting ready to send another version to CRAN soon to fix a couple bugs that were just found, so I could drop in this method as well for that release. It's now on the todo list!

[etiennebacher]

Actually extract.logitr() has to be included in texreg, not in logitr so it can't be included in your next release. But you can make a PR on texreg repo to include it

[jhelvy]

Ah woops my misunderstanding. I thought it was similar to {broom} where the parent package supplies the method:

for maintainability purposes, the broom package authors now ask that requests for new methods be first directed to the parent package (i.e. the package that supplies the model object) rather than to broom.

I was planning on adding methods for broom too, so I'll try to add support for both of these approaches for extracting results from model objects.

[etiennebacher]

Well I didn't think about including extract.logitr in logitr but actually it could work. If logitr and texreg are loaded I suppose that extract.logitr will be available in the environment and texreg will use it. I didn't try this approach but curious to know if it would work

[jhelvy]

Okay I just added your extract.logitr() method in a new branch here, but it seems it won't work.

When I run this:

library(logitr)
library(texreg)

mnl_pref <- logitr(
    data    = yogurt,
    outcome = "choice",
    obsID   = "obsID",
    pars    = c("price", "feat", "brand")
)

texreg(mnl_pref, stars = c(0.01, 0.05, 0.1))

I get this error:

Error in extract(l[[i]], ...) : 
  Neither texreg nor broom supports models of class logitr.

So perhaps I need to add it to {texreg}.

[etiennebacher]

Adding

setMethod("extract", signature = className("logitr", "logitr"),
          definition = extract.logitr)

after extract.logitr makes it work (I don't know for sure that this is good practice, I copied this piece of code directly from texreg).

---

Edit: whoops maybe I spoke too quickly

> devtools::load_all(".")
ℹ Loading logitr
 Error in setMethod("extract", signature = className("logitr", "logitr"), : 
no existing definition for function ‘extract’

[etiennebacher]

Indeed, I think adding it to texreg (with some tests) would be easier, and it would ensure that everything works in case texreg code changes later

[jhelvy]

Alright folks, this should finally be resolved. Version 0.8.0 (now on CRAN) has methods added for both {gtsummary} and {texreg} support as options for printing summary tables. I also added a new vignette on how to use these tools: https://jhelvy.github.io/logitr/articles/summarizing_results.html#formatted-summary-tables