The most common function used to view an overall summary of a model
is the summary() function:
library(logitr)
model <- logitr(
data = yogurt,
outcome = "choice",
obsID = "obsID",
pars = c("price", "feat", "brand")
)
summary(model)
#> =================================================
#>
#> Model estimated on: Fri Nov 01 19:42:01 2024
#>
#> Using logitr version: 1.1.2
#>
#> Call:
#> logitr(data = yogurt, outcome = "choice", obsID = "obsID", pars = c("price",
#> "feat", "brand"))
#>
#> Frequencies of alternatives:
#> 1 2 3 4
#> 0.402156 0.029436 0.229270 0.339138
#>
#> Exit Status: 3, Optimization stopped because ftol_rel or ftol_abs was reached.
#>
#> Model Type: Multinomial Logit
#> Model Space: Preference
#> Model Run: 1 of 1
#> Iterations: 21
#> Elapsed Time: 0h:0m:0.02s
#> Algorithm: NLOPT_LD_LBFGS
#> Weights Used?: FALSE
#> Robust? FALSE
#>
#> Model Coefficients:
#> Estimate Std. Error z-value Pr(>|z|)
#> price -0.366555 0.024365 -15.0441 < 2.2e-16 ***
#> feat 0.491439 0.120062 4.0932 4.254e-05 ***
#> brandhiland -3.715477 0.145417 -25.5506 < 2.2e-16 ***
#> brandweight -0.641138 0.054498 -11.7645 < 2.2e-16 ***
#> brandyoplait 0.734519 0.080642 9.1084 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Log-Likelihood: -2656.8878790
#> Null Log-Likelihood: -3343.7419990
#> AIC: 5323.7757580
#> BIC: 5352.7168000
#> McFadden R2: 0.2054148
#> Adj McFadden R2: 0.2039195
#> Number of Observations: 2412.0000000The summary prints out a table of the model coefficients as well as other information about the model, such as the log-likelihood, the number of observations, etc.
You can extract the coefficient from the summary table as a
data.frame using coef(summary(model)):
coefs <- coef(summary(model))
coefs
#> Estimate Std. Error z-value Pr(>|z|)
#> price -0.3665546 0.02436526 -15.044150 0.000000e+00
#> feat 0.4914392 0.12006175 4.093221 4.254226e-05
#> brandhiland -3.7154773 0.14541671 -25.550553 0.000000e+00
#> brandweight -0.6411384 0.05449794 -11.764450 0.000000e+00
#> brandyoplait 0.7345195 0.08064229 9.108366 0.000000e+00Another approach for extracting the model coefficients as a data
frame is to use the tidy() function from the {broom} package:
library(broom)
coefs <- tidy(model)
coefs
#> # A tibble: 5 × 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 price -0.367 0.0244 -15.0 0
#> 2 feat 0.491 0.120 4.09 0.0000425
#> 3 brandhiland -3.72 0.145 -25.6 0
#> 4 brandweight -0.641 0.0545 -11.8 0
#> 5 brandyoplait 0.735 0.0806 9.11 0The tidy() function returns a tibble and
provides a more standardized output and interfaces well with other
packages. You can also append a confidence interval to the data
frame:
coefs <- tidy(model, conf.int = TRUE, conf.level = 0.95)
coefs
#> # A tibble: 5 × 7
#> term estimate std.error statistic p.value conf.low conf.high
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 brandhiland -3.72 0.145 -25.6 0 -4.00 -3.43
#> 2 brandweight -0.641 0.0545 -11.8 0 -0.748 -0.534
#> 3 brandyoplait 0.735 0.0806 9.11 0 0.577 0.892
#> 4 feat 0.491 0.120 4.09 0.0000425 0.256 0.726
#> 5 price -0.367 0.0244 -15.0 0 -0.415 -0.319You can also extract other values of interest at the solution, such as:
The estimated coefficients
coef(model)
#> price feat brandhiland brandweight brandyoplait
#> -0.3665546 0.4914392 -3.7154773 -0.6411384 0.7345195The estimated standard errors
se(model)
#> price feat brandhiland brandweight brandyoplait
#> 0.02436526 0.12006175 0.14541671 0.05449794 0.08064229The log-likelihood
logLik(model)
#> 'log Lik.' -2656.888 (df=5)The variance-covariance matrix
vcov(model)
#> price feat brandhiland brandweight
#> price 0.0005936657 5.729619e-04 0.001851795 1.249988e-04
#> feat 0.0005729619 1.441482e-02 0.000855011 5.092011e-06
#> brandhiland 0.0018517954 8.550110e-04 0.021146019 1.490080e-03
#> brandweight 0.0001249988 5.092012e-06 0.001490080 2.970026e-03
#> brandyoplait -0.0015377721 -1.821331e-03 -0.003681036 7.779428e-04
#> brandyoplait
#> price -0.0015377721
#> feat -0.0018213311
#> brandhiland -0.0036810363
#> brandweight 0.0007779427
#> brandyoplait 0.0065031782You can also view a summary of statistics about the model using the
glance() function from the {broom} package:
glance(model)
#> # A tibble: 1 × 7
#> logLik null.logLik AIC BIC r.squared adj.r.squared nobs
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 -2657. -3344. 5324. 5353. 0.205 0.204 2412Often times you will need to create summary tables that are formatted
for publication. The {gtsummary} package
offers a convenient solution that works well with
logitrmodels. For example, a formatted summary table can be
obtained using the tbl_regression() function:
library(gtsummary)
model |>
tbl_regression()| Characteristic | Beta | 95% CI1 | p-value |
|---|---|---|---|
| brand | |||
| dannon | — | — | |
| hiland | -3.7 | -4.0, -3.4 | |
| weight | -0.64 | -0.75, -0.54 | |
| yoplait | 0.73 | 0.58, 0.89 | |
| feat | 0.49 | 0.26, 0.73 | |
| price | -0.37 | -0.41, -0.32 | |
| 1 CI = Confidence Interval | |||
The tbl_regression() function has many options
for customizing the output table. For example, you can change the
coefficient names with the label argument:
model |>
tbl_regression(
label = list(
feat = "Newspaper ad shown?",
brand = "Yogurt's brand"
)
)| Characteristic | Beta | 95% CI1 | p-value |
|---|---|---|---|
| Yogurt's brand | |||
| dannon | — | — | |
| hiland | -3.7 | -4.0, -3.4 | |
| weight | -0.64 | -0.75, -0.53 | |
| yoplait | 0.73 | 0.58, 0.89 | |
| Newspaper ad shown? | 0.49 | 0.26, 0.72 | |
| price | -0.37 | -0.41, -0.32 | |
| 1 CI = Confidence Interval | |||
The {gtsummary} package supports a wide variety of output types,
including support for LaTeX. One you create the
table with tbl_regression(), you can print it a variety of
ways. For example, once you’ve created a table x,
x <- model |>
tbl_regression()you can print it to LaTeX with any of the following ways:
as_gt(x) |> gt::as_latex()as_kable_extra(x, format = "latex")as_hux_table(x) |> to_latex()as_kable(x, format = "latex")Multiple models can also be printed in the same table:
model1 <- model
model2 <- logitr(
data = yogurt,
outcome = "choice",
obsID = "obsID",
pars = c("price*feat", "brand")
)
# Make individual tables
t1 <- tbl_regression(model1)
t2 <- tbl_regression(model2)
# Merge tables
tbl_merge(
tbls = list(t1, t2),
tab_spanner = c("**Baseline**", "**Interaction**")
)| Characteristic |
Baseline
|
Interaction
|
||||
|---|---|---|---|---|---|---|
| Beta | 95% CI1 | p-value | Beta | 95% CI1 | p-value | |
| brand | ||||||
| dannon | — | — | — | — | ||
| hiland | -3.7 | -4.0, -3.4 | -3.7 | -4.0, -3.4 | ||
| weight | -0.64 | -0.75, -0.53 | -0.64 | -0.75, -0.54 | ||
| yoplait | 0.73 | 0.57, 0.89 | 0.72 | 0.57, 0.88 | ||
| feat | 0.49 | 0.25, 0.73 | 1.2 | 0.42, 1.9 | 0.002 | |
| price | -0.37 | -0.41, -0.32 | -0.36 | -0.41, -0.31 | ||
| price * feat | -0.09 | -0.18, 0.01 | 0.068 | |||
| 1 CI = Confidence Interval | ||||||
Another option for obtaining a formatted table is to use the {texreg} package. This is particularly useful for obtaining tables formatted for use in LaTeX.
For example, you can print a summary to the screen using
screenreg():
library(texreg)
screenreg(model, 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.1Likewise, you can print the LaTeX code for a summary table using
texreg()
library(texreg)
texreg(model, 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}Similar to {gtsummary}, multiple models can be printed using
screenreg() or texreg():
model1 <- model
model2 <- logitr(
data = yogurt,
outcome = "choice",
obsID = "obsID",
pars = c("price*feat", "brand")
)
screenreg(
list(
model1,
model2
),
stars = c(0.01, 0.05, 0.1),
custom.model.names = c("Baseline", "Interaction")
)
#>
#> ==========================================
#> Baseline Interaction
#> ------------------------------------------
#> price -0.37 *** -0.36 ***
#> (0.02) (0.02)
#> feat 0.49 *** 1.16 ***
#> (0.12) (0.38)
#> brandhiland -3.72 *** -3.72 ***
#> (0.15) (0.15)
#> brandweight -0.64 *** -0.64 ***
#> (0.05) (0.05)
#> brandyoplait 0.73 *** 0.72 ***
#> (0.08) (0.08)
#> price:feat -0.09 *
#> (0.05)
#> ------------------------------------------
#> Num. obs. 2412 2412
#> Log Likelihood -2656.89 -2655.54
#> AIC 5323.78 5323.08
#> BIC 5352.72 5357.81
#> ==========================================
#> *** p < 0.01; ** p < 0.05; * p < 0.1