Input validation

Good functions check their inputs and fail with helpful messages:

my_function <- function(data, var_name) {
  
  # Check data type
  if (!is.data.frame(data)) {
    stop("'data' must be a data frame")
  }
  
  # Check if column exists
  if (!var_name %in% names(data)) {
    stop(paste("Column", var_name, "not found in data"))
  }
  
  # Check if variable is the right type
  
  # ... rest of function
}

Passing column names to functions

Two main approaches in the tidyverse:

Approach 1: Column name as string + .data[[]]

group_mean <- function(data, group_col, value_col) {
  data |> 
    group_by(.data[[group_col]]) |> 
    summarize(mean_value = mean(.data[[value_col]], na.rm = TRUE))
}

group_mean(penguins, "species", "bill_length_mm")

## # A tibble: 3 × 2
##   species   mean_value
##   <fct>          <dbl>
## 1 Adelie          38.8
## 2 Chinstrap       48.8
## 3 Gentoo          47.5

Approach 2: Bare column name + {{ }} (tidy evaluation)

group_mean2 <- function(data, group_col, value_col) {
  data |> 
    group_by({{ group_col }}) |> 
    summarize(mean_value = mean({{ value_col }}, na.rm = TRUE))
}

group_mean2(penguins, species, bill_length_mm)

## # A tibble: 3 × 2
##   species   mean_value
##   <fct>          <dbl>
## 1 Adelie          38.8
## 2 Chinstrap       48.8
## 3 Gentoo          47.5

Common plot types

library(gridExtra)

# histogram
p1 <- ggplot(penguins, aes(x = bill_length_mm)) +
  geom_histogram(binwidth = 2) +
  labs(
  title = "Histogram",
  subtitle = "Use for the distribution of a numeric variable",
  x = "Bill length (mm)",
  y = "Count"
  ) +
  theme_minimal(base_size = 10)


# boxplot
p2 <- ggplot(penguins, aes(x = species, y = body_mass_g, fill = species)) +
  geom_boxplot(alpha = 0.8) +
  labs(
  title = "Boxplot",
  subtitle = "Use to compare distributions across categories",
  x = "Species",
  y = "Body mass (g)"
  ) +
  theme_minimal(base_size = 10) +
  theme(legend.position = "none")

# scatterplot
p3 <- ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g)) +
  geom_point(alpha = 0.7) +
  labs(
  title = "Scatterplot",
  subtitle = "Use for relationships between two numeric variables",
  x = "Flipper length (mm)",
  y = "Body mass (g)"
  ) +
  theme_minimal(base_size = 10)

# another scatterplot
p4 <- ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g, color = species)) +
  geom_point(alpha = 0.7) +
  labs(
  title = "Scatterplot by Group",
  subtitle = "Use when comparing numeric relationships across \ngroups",
  x = "Flipper length (mm)",
  y = "Body mass (g)",
  color = "Species"
  ) +
  theme_minimal(base_size = 10)


grid.arrange(p1, p2, p3, p4, ncol = 2)

Faceting

Another way to show groups is to split your plot into panels. An example:

ggplot(penguins, aes(x = bill_length_mm, fill = species)) +
  geom_histogram(binwidth = 2, color = "white") +
  labs(
    title = "Faceted Histogram",
    subtitle = "Use to compare the distribution of a numeric variable across groups",
    x = "Bill length (mm)",
    y = "Count"
  ) +
  facet_wrap(~ species, ncol = 1) +
  theme_minimal() +
  theme(legend.position = "none")

Making plots standalone

A good plot should be understandable without reading surrounding text. Be mindful of labels, titles, legends, and captions.

ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g, color = species)) +
  geom_point(alpha = 0.7) +
  labs(
    title = "Penguin Body Mass vs. Flipper Length",
    subtitle = "Data from Palmer Station, Antarctica (2007-2009)",
    x = "Flipper length (mm)",
    y = "Body mass (g)",
    color = "Species",
    caption = "Source: palmerpenguins  package"
  ) +
  theme_minimal()

Running models

model <- lm(body_mass_g ~ flipper_length_mm + species, data = penguins)
summary(model)

## 
## Call:
## lm(formula = body_mass_g ~ flipper_length_mm + species, data = penguins)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -927.70 -254.82  -23.92  241.16 1191.68 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -4031.477    584.151  -6.901 2.55e-11 ***
## flipper_length_mm    40.705      3.071  13.255  < 2e-16 ***
## speciesChinstrap   -206.510     57.731  -3.577 0.000398 ***
## speciesGentoo       266.810     95.264   2.801 0.005392 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 375.5 on 338 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.7826, Adjusted R-squared:  0.7807 
## F-statistic: 405.7 on 3 and 338 DF,  p-value: < 2.2e-16

Tidying with broom

library(broom)

# coefficients as a tidy data frame
tidy(model)

## # A tibble: 4 × 5
##   term              estimate std.error statistic  p.value
##   <chr>                <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)        -4031.     584.       -6.90 2.55e-11
## 2 flipper_length_mm     40.7      3.07     13.3  1.40e-32
## 3 speciesChinstrap    -207.      57.7      -3.58 3.98e- 4
## 4 speciesGentoo        267.      95.3       2.80 5.39e- 3

# model fit statistics
glance(model)

## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic   p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>     <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1     0.783         0.781  376.      406. 1.25e-111     3 -2511. 5032. 5051.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

# add fitted values and residuals to data
augment(model) |> head()

## # A tibble: 6 × 10
##   .rownames body_mass_g flipper_length_mm species .fitted .resid    .hat .sigma
##   <chr>           <int>             <int> <fct>     <dbl>  <dbl>   <dbl>  <dbl>
## 1 1                3750               181 Adelie    3336.  414.  0.0120    375.
## 2 2                3800               186 Adelie    3540.  260.  0.00767   376.
## 3 3                3250               195 Adelie    3906. -656.  0.00833   374.
## 4 5                3450               193 Adelie    3825. -375.  0.00724   376.
## 5 6                3650               190 Adelie    3703.  -52.5 0.00662   376.
## 6 7                3625               181 Adelie    3336.  289.  0.0120    376.
## # ℹ 2 more variables: .cooksd <dbl>, .std.resid <dbl>

Presenting with modelsummary

library(modelsummary)

# fit multiple models
m1 <- lm(body_mass_g ~ flipper_length_mm, data = penguins)
m2 <- lm(body_mass_g ~ flipper_length_mm + species, data = penguins)
m3 <- lm(body_mass_g ~ flipper_length_mm * species, data = penguins)

models <- list(
  "Flipper only" = m1,
  "+ Species" = m2,
  "Interaction" = m3
)

# regression table
modelsummary(
  models,
  gof_omit = "AIC|BIC|Log.Lik|RMSE",
  output = "gt"
) |> 
  tab_options(table.align = "center")

# coefficients plot
modelplot(models, coef_omit = "Intercept") +
  theme_minimal() +
  labs(title = "Coefficients")

This script was drafted by Killian Conyngham and Carol Sobral.