Getting started
To begin, load the modelsummary package and download data from the Rdatasets archive:
library(modelsummary)
library(tinytable)
url <- 'https://vincentarelbundock.github.io/Rdatasets/csv/HistData/Guerry.csv'
dat <- read.csv(url)
dat$Small <- dat$Pop1831 > median(dat$Pop1831)
dat <- dat[,
c("Donations", "Literacy", "Commerce", "Crime_pers", "Crime_prop", "Clergy", "Small")
]Data Summaries
Quick overview of the data:
datasummary_skim(dat)| Unique | Missing Pct. | Mean | SD | Min | Median | Max | Histogram | |
|---|---|---|---|---|---|---|---|---|
| Donations | 85 | 0 | 7075.5 | 5834.6 | 1246.0 | 5020.0 | 37015.0 |  |
| Literacy | 50 | 0 | 39.3 | 17.4 | 12.0 | 38.0 | 74.0 |  |
| Commerce | 84 | 0 | 42.8 | 25.0 | 1.0 | 42.5 | 86.0 |  |
| Crime_pers | 85 | 0 | 19754.4 | 7504.7 | 2199.0 | 18748.5 | 37014.0 |  |
| Crime_prop | 86 | 0 | 7843.1 | 3051.4 | 1368.0 | 7595.0 | 20235.0 |  |
| Clergy | 85 | 0 | 43.4 | 25.0 | 1.0 | 43.5 | 86.0 |  |
| Small | N | % | ||||||
| FALSE | 43 | 50.0 | ||||||
| TRUE | 43 | 50.0 |
Balance table (aka “Table 1”) with differences in means by subgroups:
datasummary_balance(~Small, dat)| FALSE (N=43) | TRUE (N=43) | |||||
|---|---|---|---|---|---|---|
| Mean | Std. Dev. | Mean | Std. Dev. | Diff. in Means | Std. Error | |
| Donations | 7258.5 | 6194.1 | 6892.6 | 5519.0 | -365.9 | 1265.2 |
| Literacy | 37.9 | 19.1 | 40.6 | 15.6 | 2.7 | 3.8 |
| Commerce | 42.7 | 24.6 | 43.0 | 25.7 | 0.3 | 5.4 |
| Crime_pers | 18040.6 | 7638.4 | 21468.2 | 7044.3 | 3427.7 | 1584.6 |
| Crime_prop | 8422.5 | 3406.7 | 7263.7 | 2559.3 | -1158.8 | 649.8 |
| Clergy | 39.1 | 26.7 | 47.7 | 22.7 | 8.6 | 5.3 |
Correlation table:
datasummary_correlation(dat)| Donations | Literacy | Commerce | Crime_pers | Crime_prop | Clergy | |
|---|---|---|---|---|---|---|
| Donations | 1 | . | . | . | . | . |
| Literacy | -.13 | 1 | . | . | . | . |
| Commerce | .30 | -.58 | 1 | . | . | . |
| Crime_pers | -.04 | -.04 | .05 | 1 | . | . |
| Crime_prop | -.13 | -.37 | .41 | .27 | 1 | . |
| Clergy | .09 | -.17 | -.12 | .26 | -.07 | 1 |
Two variables and two statistics, nested in subgroups:
datasummary(Literacy + Commerce ~ Small * (mean + sd), dat)| FALSE | TRUE | |||
|---|---|---|---|---|
| mean | sd | mean | sd | |
| Literacy | 37.88 | 19.08 | 40.63 | 15.57 |
| Commerce | 42.65 | 24.59 | 42.95 | 25.75 |
Model Summaries
Estimate a linear model and display the results:
mod <- lm(Donations ~ Crime_prop, data = dat)
modelsummary(mod)| (1) | |
|---|---|
| (Intercept) | 9065.287 |
| (1738.926) | |
| Crime_prop | -0.254 |
| (0.207) | |
| Num.Obs. | 86 |
| R2 | 0.018 |
| R2 Adj. | 0.006 |
| AIC | 1739.0 |
| BIC | 1746.4 |
| Log.Lik. | -866.516 |
| F | 1.505 |
| RMSE | 5749.29 |
Now estimate five regression models, display the results side-by-side, and use the group_tt() function from the tinytable package to add column labels:
library(tinytable)
models <- list(
"I" = lm(Donations ~ Literacy + Clergy, data = dat),
"II" = lm(Crime_pers ~ Literacy + Clergy, data = dat),
"III" = lm(Crime_prop ~ Literacy + Clergy, data = dat),
"IV" = glm(Crime_pers ~ Literacy + Commerce, family = poisson, data = dat),
"V" = glm(Donations ~ Literacy + Commerce, family = poisson, data = dat)
)
modelsummary(models, stars = TRUE, gof_omit = "IC|Adj|F|RMSE|Log") |>
group_tt(j = list("Linear" = 2:4, "Poisson" = 5:6))| Linear | Poisson | ||||
|---|---|---|---|---|---|
| I | II | III | IV | V | |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |||||
| (Intercept) | 7948.667*** | 16259.384*** | 11243.544*** | 9.876*** | 8.241*** |
| (2078.276) | (2611.140) | (1011.240) | (0.003) | (0.006) | |
| Literacy | -39.121 | 3.680 | -68.507*** | 0.000*** | 0.003*** |
| (37.052) | (46.552) | (18.029) | (0.000) | (0.000) | |
| Clergy | 15.257 | 77.148* | -16.376 | ||
| (25.735) | (32.334) | (12.522) | |||
| Commerce | 0.001*** | 0.011*** | |||
| (0.000) | (0.000) | ||||
| Num.Obs. | 86 | 86 | 86 | 86 | 86 |
| R2 | 0.020 | 0.065 | 0.152 | ||
Now, save it to a Microsoft Word document:
modelsummary(models, output = "table.docx")And draw a coefficient plot:
ols <- models[1:3]
modelplot(ols, coef_omit = "Intercept")