Today’s goals

• Make a simple frequency tables
• Make a more complex “Table 1”
• Run regressions in R and show the results in a table
• Introduce the framework of ggplot2 and make some simple figures

Tables

We saw that we could compute summary statistics using summarize():

nlsy |> 
  group_by(sex_cat) |> 
  summarize(prop_glasses = mean(glasses),
            mean_age_bir = mean(age_bir),
            n_vg_eye = sum(eyesight_cat == "Very Good"),
            prop_vg_eye = mean(eyesight_cat == "Very Good"),
            n_g_eye = sum(eyesight_cat == "Good"),
            prop_g_eye = mean(eyesight_cat == "Good")
  )

# A tibble: 2 × 7
  sex_cat prop_glasses mean_age_bir n_vg_eye prop_vg_eye n_g_eye prop_g_eye
  <fct>          <dbl>        <dbl>    <int>       <dbl>   <int>      <dbl>
1 Male           0.441         25.1      162       0.323      85      0.170
2 Female         0.572         22.2      223       0.317     164      0.233

Easier way to make frequency tables

The {janitor} package has nice functionality to create tables:

# install.packages("janitor")
library(janitor)
tabyl(nlsy, sex_cat)

 sex_cat   n   percent
    Male 501 0.4157676
  Female 704 0.5842324

Cross-tabulations with tabyl()

It’s easy to make a 2 x 2 (or n x m) table with tabyl()

tabyl(nlsy, sex_cat, eyesight_cat)

 sex_cat Excellent Very Good Good Fair Poor
    Male       228       162   85   21    5
  Female       246       223  164   57   14

tabyl()

We might want to see the percentages instead of counts

nlsy |> 
  tabyl(sex_cat, glasses_cat) |> 
  adorn_percentages() |> 
  adorn_pct_formatting()

 sex_cat Doesn't wear glasses Wears glasses/contacts
    Male                55.9%                  44.1%
  Female                42.8%                  57.2%

Other percentages

nlsy |>
    tabyl(sex_cat, glasses_cat) |>  
    adorn_percentages("col") |>             
    adorn_pct_formatting()

 sex_cat Doesn't wear glasses Wears glasses/contacts
    Male                48.2%                  35.4%
  Female                51.8%                  64.6%

nlsy |>
    tabyl(sex_cat, glasses_cat) |>  
    adorn_percentages("all") |>             
    adorn_pct_formatting()

 sex_cat Doesn't wear glasses Wears glasses/contacts
    Male                23.2%                  18.3%
  Female                25.0%                  33.4%

tabyl()

These are just dataframes, so we can save as csv, etc.

two_by_two <- nlsy |>
    tabyl(sex_cat, glasses_cat)

write_csv(two_by_two, 
          here::here("results", "sex-eyesight-table.csv"))

We can easily do a chi-squared test

tabyl(nlsy, sex_cat, eyesight_cat) |> 
  chisq.test()

    Pearson's Chi-squared test

data:  tabyl(nlsy, sex_cat, eyesight_cat)
X-squared = 22.738, df = 4, p-value = 0.0001428

Making more complex tables

There are lots of packages for making tables in R

One of my favorites is {gtsummary}

# install.packages("gtsummary")
library(gtsummary)

gtsummary::tbl_summary()

tbl_summary(
  nlsy,
  by = sex_cat,
  include = c(race_eth_cat,
              eyesight_cat, 
              glasses, 
              age_bir))

tbl_summary()

• Incredibly customizeable
• Really helpful with Table 1
• I often just view in the web browser and copy and paste into a Word document
• Can also be used within quarto/R Markdown
• If output is Word, I use as_flex_table() to output using the flextable package
• Make even more customizeable with the gt package with as_gt()

Regression

Regressions take a formula: y ~ x1 + x2 + x3

• Include interaction terms between x1 and x2 with y ~ x1*x2 + x3
• Main effects of x1 and x2 will be included too
• Indicator (“dummy”) variables will automatically be created for factors
• The first level will be the reference level
• If you want to include a squared term (for example), you can make the squared variable first, or wrap in I(): y ~ x1 + I(x1^2)

Regression

To fit a linear regression (by ordinary least squares), use the lm() function

To fit a generalized linear model (e.g., logistic regression, Poisson regression) use glm() and specify the family = argument

• family = gaussian() is the default: another way of fitting a linear regresion
• family = binomial() gives you logistic regression
• family = poisson() is Poisson regression

We tell R what dataset to pull the variables from with data =

Regression

linear_model <- lm(income ~ sex_cat*age_bir + race_eth_cat, 
                   data = nlsy)

logistic_model <- glm(glasses ~ eyesight_cat + sex_cat + income, 
                      data = nlsy, family = binomial())

Regression

coef(linear_model)

                        (Intercept)                       sex_catFemale 
                          1029.3339                          -4681.2484 
                            age_bir                   race_eth_catBlack 
                           482.4650                           -299.6490 
race_eth_catNon-Black, Non-Hispanic               sex_catFemale:age_bir 
                          6418.4568                            161.5008 

confint(linear_model)

                                           2.5 %    97.5 %
(Intercept)                          -3746.58543 5805.2531
sex_catFemale                       -10553.32090 1190.8242
age_bir                                303.50836  661.4217
race_eth_catBlack                    -2448.39108 1849.0932
race_eth_catNon-Black, Non-Hispanic   4500.91244 8336.0012
sex_catFemale:age_bir                  -77.30931  400.3109

Regression

exp(coef(logistic_model))

          (Intercept) eyesight_catVery Good      eyesight_catGood 
            0.6338301             0.9172091             0.8608297 
     eyesight_catFair      eyesight_catPoor         sex_catFemale 
            0.5798278             1.1624355             1.8362460 
               income 
            1.0000174 

summary(logistic_model)

Call:
glm(formula = glasses ~ eyesight_cat + sex_cat + income, family = binomial(), 
    data = nlsy)

Coefficients:
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)           -4.560e-01  1.383e-01  -3.298 0.000975 ***
eyesight_catVery Good -8.642e-02  1.400e-01  -0.617 0.537031    
eyesight_catGood      -1.499e-01  1.604e-01  -0.934 0.350267    
eyesight_catFair      -5.450e-01  2.535e-01  -2.150 0.031585 *  
eyesight_catPoor       1.505e-01  4.786e-01   0.315 0.753121    
sex_catFemale          6.077e-01  1.210e-01   5.021 5.14e-07 ***
income                 1.745e-05  4.613e-06   3.782 0.000155 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1668.9  on 1204  degrees of freedom
Residual deviance: 1626.8  on 1198  degrees of freedom
AIC: 1640.8

Number of Fisher Scoring iterations: 4

Helpful regression packages

{broom} helps “tidy” regression results

library(broom)
tidy(linear_model)

# A tibble: 6 × 5
  term                                estimate std.error statistic  p.value
  <chr>                                  <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)                            1029.    2434.      0.423 6.72e- 1
2 sex_catFemale                         -4681.    2993.     -1.56  1.18e- 1
3 age_bir                                 482.      91.2     5.29  1.46e- 7
4 race_eth_catBlack                      -300.    1095.     -0.274 7.84e- 1
5 race_eth_catNon-Black, Non-Hispanic    6418.     977.      6.57  7.63e-11
6 sex_catFemale:age_bir                   162.     122.      1.33  1.85e- 1

broom::tidy()

tidy(logistic_model, conf.int = TRUE, exponentiate = TRUE)

# A tibble: 7 × 7
  term                  estimate  std.error statistic p.value conf.low conf.high
  <chr>                    <dbl>      <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)              0.634 0.138         -3.30  9.75e-4    0.483     0.830
2 eyesight_catVery Good    0.917 0.140         -0.617 5.37e-1    0.697     1.21 
3 eyesight_catGood         0.861 0.160         -0.934 3.50e-1    0.628     1.18 
4 eyesight_catFair         0.580 0.254         -2.15  3.16e-2    0.350     0.949
5 eyesight_catPoor         1.16  0.479          0.315 7.53e-1    0.458     3.08 
6 sex_catFemale            1.84  0.121          5.02  5.14e-7    1.45      2.33 
7 income                   1.00  0.00000461     3.78  1.55e-4    1.00      1.00 

broom::tidy() can also help other statistics

tabyl(nlsy, sex_cat, eyesight_cat) |> 
  chisq.test() |> 
  tidy()

# A tibble: 1 × 4
  statistic  p.value parameter method                    
      <dbl>    <dbl>     <int> <chr>                     
1      22.7 0.000143         4 Pearson's Chi-squared test

t.test(income ~ sex_cat, data = nlsy) |> 
  tidy()

# A tibble: 1 × 10
  estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high
     <dbl>     <dbl>     <dbl>     <dbl>   <dbl>     <dbl>    <dbl>     <dbl>
1    2397.    16690.    14292.      3.00 0.00273      963.     831.     3963.
# ℹ 2 more variables: method <chr>, alternative <chr>

You could put several together

linear_model_no_int <- lm(income ~ sex_cat + age_bir + race_eth_cat, data = nlsy)

tbl_no_int <- tbl_regression(
  linear_model_no_int, 
  intercept = TRUE,
  label = list(
    sex_cat ~ "Sex",
    race_eth_cat ~ "Race/ethnicity",
    age_bir ~ "Age at first birth"
  ))

tbl_int <- tbl_regression(
  linear_model, 
  intercept = TRUE,
  label = list(
    sex_cat ~ "Sex",
    race_eth_cat ~ "Race/ethnicity",
    age_bir ~ "Age at first birth",
    `sex_cat:age_bir` ~ "Sex/age interaction"
  ))

You could put several together

tbl_merge(list(tbl_no_int, tbl_int), 
          tab_spanner = c("**Model 1**", "**Model 2**"))

Regression package recommendation from Xiyue: autoReg

# install.packages("autoReg")
library(autoReg)
autoReg(linear_model) |> 
  myft()

Figures in R using ggplot()

Figures in R using ggplot()

Why ggplot?

• Powerful and flexible: create complex and customized visualizations easily
• Reproducibility and efficiency: promotes reproducibility by offering consistent syntax and saves time through automation of plot creation
• Layered approach: incrementally build visualizations with multiple layers, exploring different aspects of data
• Extensive customization: essentially infinite options to tailor visualizations

ggplot builds figures by adding on pieces via a particular “grammar of graphics”

Basic structure of a ggplot

ggplot(data = {data}, 
       aes(x = {xvar}, y = {yvar}, <characteristic> = {othvar}, ...)) +
      <geom>(<characteristic> = "value", ...) + 
      ...

What are some of the layers we may need for this one?

Returning to our basic structure

ggplot(data = {data}, 
       aes(x = {xvar}, y = {yvar}, <characteristic> = {othvar}, ...)) +
      <geom>(<characteristic> = "value", ...) + 
      ...

Let’s walk through some more examples in depth

Scatterplot: geom_point()

ggplot(data = nlsy, 
       aes(x = income, y = age_bir)) +
  geom_point()

What if we want to change the color of the points?

ggplot(data = nlsy, 
       aes(x = income, y = age_bir)) +
  geom_point(color = "blue")

When we put color = outside the aes(), it means we’re giving it a specific color value that applies to all the points.

What if we want the color to correspond to values of a variable?

ggplot(data = nlsy, 
       aes(x = income, y = age_bir, 
           color = eyesight_cat)) +
  geom_point()

When we put color = inside the aes() – with no quotation marks – it means we’re telling it how it should assign colors.

Alternative specification

ggplot(data = nlsy, 
       aes(x = income, y = age_bir)) +
  geom_point(aes(color = eyesight_cat))

Note that we could also put the aes() (aesthetics) in the geom_() itself.

Let’s change the colors

ggplot(data = nlsy, 
       aes(x = income, y = age_bir, 
           color = eyesight_cat)) +
  geom_point() +
  scale_color_manual(
    values = c("blue", "purple", "red",
               "green", "yellow"))

We add on another layer to specify the colors we want.

Color palettes

ggplot(data = nlsy, 
       aes(x = income, y = age_bir, 
           color = eyesight_cat)) +
  geom_point()  +
  scale_color_brewer(palette = "Set1")

There are tons of different options in R for color palettes.

You can play around with those in the RColorBrewer package here.

You can access the scales in that package with scale_color_brewer()

Change the title on the legend

ggplot(data = nlsy, 
       aes(x = income, y = age_bir, 
           color = eyesight_cat)) +
  geom_point()  +
  scale_color_brewer(palette = "Set1", 
                  name = "Eyesight")

Change the axis scale

ggplot(data = nlsy, 
       aes(x = income, y = age_bir, 
           color = eyesight_cat)) +
  geom_point()  +
  scale_x_log10()

There are a lot of scale_x_() and scale_y_() functions for you to explore

We can label the axis better

ggplot(nlsy,
       aes(x = income, y = age_bir, 
           color = eyesight_cat)) +
  geom_point()  +
  scale_x_log10(labels = scales::dollar, 
                limits = c(100, 100000),
                breaks = c(100, 1000, 10000, 100000),
                name = "Income (USD)")

The {scales} packages contains lots of helpful number formatting functions

See the the examples in help(scale_x_log10) for some of them

Facets

• One of the most useful features of {ggplot2} is the ability to “facet” a graph by splitting it up according to the values of some variable
• You might use this to show results for a lot of outcomes or exposures at once, for example, or see how some relationship differs by something like age or geographic region

ggplot(data = nlsy, aes(x = nsibs)) +
  geom_bar() +
  labs(x = "Number of siblings")

We’ll introduce bar graphs at the same time!

Notice how we only need an x = argument - the y-axis is automatically the count with this geom.

ggplot(data = nlsy, aes(x = nsibs)) +
  geom_bar() +
  labs(x = "Number of siblings") +
  facet_grid(cols = vars(region_cat))

The facet_grid() function splits up the data according to a variable(s).

Here we’ve split it by region into columns.

ggplot(data = nlsy, aes(x = nsibs)) +
  geom_bar() +
  labs(x = "Number of siblings") +
  facet_grid(rows = vars(region_cat))

Since that was hard to read, we’ll probably want to split by rows instead.

ggplot(data = nlsy, aes(x = nsibs)) +
  geom_bar() +
  labs(x = "Number of siblings") +
  facet_grid(rows = vars(region_cat),
             margins = TRUE)

We can also add a row for all of the data combined.

ggplot(data = nlsy, aes(x = nsibs)) +
  geom_bar() +
  labs(x = "Number of siblings") +
  facet_grid(rows = vars(region_cat),
             margins = TRUE,
             scales = "free_y")

That squishes the other rows though! We can allow them all to have their own axis limits with the scales = argument.

Other options are “free_x” if we want to allow the x-axis scale to vary, or just “free” to allow for both.

ggplot(data = nlsy, aes(x = nsibs)) +
  geom_bar() +
  labs(x = "Number of siblings") +
  facet_wrap(vars(region_cat))

We can use facet_wrap() instead, if we want to use both multiple rows and columns for all the values of a variable.

ggplot(data = nlsy, aes(x = nsibs)) +
  geom_bar() +
  labs(x = "Number of siblings") +
  facet_wrap(vars(region_cat), 
             ncol = 3)

It tries to make a good decision, but you can override how many columns you want!

Wait, these look like histograms!

When we have a variable with a lot of possible values, we may want to bin them with a histogram

ggplot(nlsy, aes(x = income)) +
  geom_histogram()

stat_bin() using bins = 30. Pick better value with binwidth.

We used discrete values with geom_bar(), but with geom_histogram() we’re combining values: the default is into 30 bins.

This is one of the most common warning messages I get in R!

ggplot(nlsy, aes(x = income)) +
  geom_histogram(bins = 10)

We can use bins = instead, if we want!

ggplot(nlsy, aes(x = income)) +
  geom_histogram(binwidth = 1000)

Sometimes the bin width actually has some meaning so we want to specify that

Here, each bin is $1000 – you can see the $5000 and $10000 increments

Themes to make our plots prettier

You probably recognize the ggplot theme. But did you know you can trick people into thinking you made your figures in Stata?

p +
  theme_minimal()

We can change the overall theme

Since we stored the plot as p, it’s easy to add on / try different things

p +
  ggthemes::theme_fivethirtyeight()

Other packages contain themes, too.

p +
  ggthemes::theme_excel_new()

In case you miss Excel….

p +
  louisahstuff::my_theme()

You can even make your own!

Finally, save it!

If your data changes, you can easily run the whole script again:

library(tidyverse)
dataset <- read_csv("dataset.csv")
ggplot(dataset) + 
  geom_point(aes(x = xvar, y = yvar))
ggsave(filename = "scatterplot.pdf")

The ggsave() function will automatically save the most recent plot in your output.

To be safe, you can store your plot, e.g., p <- ggplot(...) + ... and then

ggsave(filename = "scatterplot.pdf", plot = p)

More resources

• Cheat sheet: https://www.rstudio.com/resources/cheatsheets/#ggplot2
• Catalog: http://shiny.stat.ubc.ca/r-graph-catalog/
• Cookbook: http://www.cookbook-r.com/Graphs/
• Official package reference: https://ggplot2.tidyverse.org/index.html
• List of themes and instructions to make your own: https://www.datanovia.com/en/blog/ggplot-themes-gallery/
• Book (includes theory behind ggplot): https://ggplot2-book.org/

Today’s summary

• We learned how to use janitor and gtsummary packages to make tables
• We learned how to fit linear regressions and generlized linear models
• We learned how to use broom to tidy regression results
• We learned the basics of ggplot2

Today’s functions

• tabyl(): create frequency tables
• tbl_summary(): create summary tables with a lot of covariates
• lm(): fit linear regression models
• glm(): fit generalized linear models
• tbl_regression(): create regression tables
• tidy(): tidy regression results
• ggplot(): create figures
• geom_...(): add a geometric feature to a figure
• scale_...(): change the scale of an axis or aesthetic
• facet_...(): split a figure into panels
• theme_...(): change the look of a figure