Tidy Data

Introduction

Tidy data is a standard way of organizing data values within a dataset (Wickham 2014). A dataset is tidy when:

1. Each variable forms a column.
2. Each observation forms a row.
3. Each type of observational unit forms a table.

Most real-world data arrives in an untidy form. The tidyr package provides tools to reshape data into tidy format, after which the dplyr package provides a consistent grammar for transformation and summarization. We will use both of these in this notebook to demontrate principles of tidy data analysis.

suppressPackageStartupMessages({

library(tidyr)
library(dplyr)

})

In this notebook we will work with toy datasets to practice identifying tidy vs. untidy structure, and use pivot_longer(), pivot_wider(), and core dplyr verbs to tidy and wrangle data.

Note: below we will use a function tibble() to make a type of object of the same name. A tibble is the tidyverse’s class for tabular data. It behaves like base R’s data.frame but with a few practical improvements: it prints only the first ten rows by default, always shows column types, and never silently converts strings to factors. Throughout this notebook we use tibble and data.frame interchangeably, but tibble is what the tidyverse functions return.

What does untidy data look like?

A common form of untidy data is a “wide” table where values of one variable are spread across column names. Consider a table of gene expression counts where each column is a sample. The three genes are:

• GAPDH: glycolytic enzyme expressed in nearly all cells; widely used as a housekeeping reference gene because its expression is constitutive and stable.
• IL6: pro-inflammatory cytokine that activates acute-phase responses; strongly induced by bacterial infection through NF-κB signaling.
• IRF1: transcription factor and direct target of the IFNg–JAK–STAT1 axis; a canonical marker of the interferon response.

counts_wide <- tibble(
  gene = c("GAPDH", "IL6", "IRF1"),
  control = c(520, 30, 50),
  IFNg = c(560, 90, 290),
  Salmonella = c(500, 320, 310)
)
counts_wide

# A tibble: 3 × 4
  gene  control  IFNg Salmonella
  <chr>   <dbl> <dbl>      <dbl>
1 GAPDH     520   560        500
2 IL6        30    90        320
3 IRF1       50   290        310

Here control, IFNg, and Salmonella are values of a condition variable, not variables themselves. To make this tidy we need pivot_longer().

There is nothing inherently wrong with the above presentation, and in some contexts we may prefer the data as a matrix (for example for perfoming matrix algebra). The point here is to note that different operations are easier to compute on one representation or the other.

NoteQuestion

How would you plot this data? Suppose we want counts on the y-axis and conditions (control, IFNg, Salmonella) on the x-axis.

pivot_longer()

pivot_longer() takes multiple columns and collapses them into two: one for the former column names and one for the values.

counts_long <- counts_wide |>
  pivot_longer(
    cols = -gene, # every column except 'gene'
    names_to = "condition",
    values_to = "count"
  )
counts_long

# A tibble: 9 × 3
  gene  condition  count
  <chr> <chr>      <dbl>
1 GAPDH control      520
2 GAPDH IFNg         560
3 GAPDH Salmonella   500
4 IL6   control       30
5 IL6   IFNg          90
6 IL6   Salmonella   320
7 IRF1  control       50
8 IRF1  IFNg         290
9 IRF1  Salmonella   310

NoteQuestion

How many rows does counts_long have compared to counts_wide? Why?

pivot_wider()

The inverse operation is pivot_wider(), which spreads a key-value pair back into separate columns. This is sometimes useful for display or for feeding data into a matrix-based function.

counts_long |>
  pivot_wider(
    names_from = condition,
    values_from = count
  )

# A tibble: 3 × 4
  gene  control  IFNg Salmonella
  <chr>   <dbl> <dbl>      <dbl>
1 GAPDH     520   560        500
2 IL6        30    90        320
3 IRF1       50   290        310

Exercise: pivot the penguins

The penguins dataset (from palmerpenguins) has four body measurement columns: bill_length_mm, bill_depth_mm, flipper_length_mm, and body_mass_g. This is already in a tidy format, as these are each measurements, and the rows are observtions (penguins). For more details, see the palmerpenguins package homepage, which notes:

Data were collected and made available by Dr. Kristen Gorman and the Palmer Station, Antarctica LTER, a member of the Long Term Ecological Research Network. … [the palmerpenguins package contains] data for 344 penguins. There are 3 different species of penguins in this dataset, collected from 3 islands in the Palmer Archipelago, Antarctica.

suppressPackageStartupMessages({
  library(palmerpenguins)
})

penguins

# A tibble: 344 × 8
   species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g sex     year
   <fct>   <fct>              <dbl>         <dbl>             <int>       <int> <fct>  <int>
 1 Adelie  Torgersen           39.1          18.7               181        3750 male    2007
 2 Adelie  Torgersen           39.5          17.4               186        3800 female  2007
 3 Adelie  Torgersen           40.3          18                 195        3250 female  2007
 4 Adelie  Torgersen           NA            NA                  NA          NA <NA>    2007
 5 Adelie  Torgersen           36.7          19.3               193        3450 female  2007
 6 Adelie  Torgersen           39.3          20.6               190        3650 male    2007
 7 Adelie  Torgersen           38.9          17.8               181        3625 female  2007
 8 Adelie  Torgersen           39.2          19.6               195        4675 male    2007
 9 Adelie  Torgersen           34.1          18.1               193        3475 <NA>    2007
10 Adelie  Torgersen           42            20.2               190        4250 <NA>    2007
# ℹ 334 more rows

NoteExercise

Even though the penguins dataset is already in a tidy format, we can still make the table into a longer format. Pivot those four columns to a long format with columns species, measurement, and value, keeping only species as the identifier.

We will return to the penguins dataset in a later exercise on modeling missing data.

TipHint

cols takes an unquoted column selector (like -species) because that column already exists and R can look it up by name. names_to and values_to take quoted strings because they are new column names you are inventing; they don’t exist in the data yet, so there is nothing for R to look up.

penguins |>
  select(species, bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g) |>
  pivot_longer(
    cols = ______,
    names_to = ______,
    values_to = ______
  )

Solution (Solution). -species selects every column except species, which is the identifier.

penguins |>
  select(species, bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g) |>
  pivot_longer(
    cols = -species,
    names_to = "measurement",
    values_to = "value"
  )

# A tibble: 1,376 × 3
   species measurement        value
   <fct>   <chr>              <dbl>
 1 Adelie  bill_length_mm      39.1
 2 Adelie  bill_depth_mm       18.7
 3 Adelie  flipper_length_mm  181  
 4 Adelie  body_mass_g       3750  
 5 Adelie  bill_length_mm      39.5
 6 Adelie  bill_depth_mm       17.4
 7 Adelie  flipper_length_mm  186  
 8 Adelie  body_mass_g       3800  
 9 Adelie  bill_length_mm      40.3
10 Adelie  bill_depth_mm       18  
# ℹ 1,366 more rows

gradethis::grade_this_code()

Core dplyr verbs

With tidy data in hand, dplyr provides five key verbs:

Verb	Purpose
filter()	Keep rows matching a condition
select()	Keep (or drop) columns
mutate()	Add or modify columns
summarize()	Reduce rows to a summary
arrange()	Sort rows
slice()	Select rows by position (e.g. slice_sample() for random rows)

These are usually combined with group_by() to apply operations within groups.

pandas is the most widely used Python package for tabular data manipulation and provides equivalents for all of these verbs. See the Python section below for more details and an example of sharing data between R and Python. In the table below, df is the pandas DataFrame; method calls are written without the df. prefix.

R (dplyr)	Python (pandas)
filter(count > 100)	query("count > 100")
select(gene, count)	[["gene", "count"]]
select(starts_with("bill"))	filter(regex="^bill", axis=1)
select_if(is.numeric)	select_dtypes(include="number")
mutate(log_count = log(count))	assign(log_count=lambda d: np.log(d["count"]))
group_by(condition) %>% summarize(total = sum(count))	groupby("condition").agg(total=("count", "sum"))
arrange(count)	sort_values("count")
slice_sample(prop = 0.1)	sample(frac=0.1)

Below we show an example of combining several verbs in a pipeline. The following code chunk can be read as: take counts_long, keep only rows where condition is not "control", group the remaining rows by condition, then compute the total and mean count within each group.

counts_long |>
  filter(condition != "control") |>
  group_by(condition) |>
  summarize(total = sum(count), mean_count = mean(count))

# A tibble: 2 × 3
  condition  total mean_count
  <chr>      <dbl>      <dbl>
1 IFNg         940       313.
2 Salmonella  1130       377.

Exercise: filter and summarize

NoteExercise

Using counts_long, keep only rows where count > 100, then compute the number of data per condition and the SD of counts per condition.

TipHint

The function n() will tell you how many rows are present, and it respects grouping by group_by().

counts_long |>
  filter(______) |>
  group_by(______) |>
  summarize(n_remain = ___, sd_count = ______)

Solution (Solution).

counts_long |>
  filter(count > 100) |>
  group_by(condition) |>
  summarize(n_remain = n(), sd_count = sd(count))

# A tibble: 3 × 3
  condition  n_remain sd_count
  <chr>         <int>    <dbl>
1 IFNg              2     191.
2 Salmonella        3     107.
3 control           1      NA 

gradethis::grade_this_code()

Plotting tidy data with plot()

One payoff of tidy data is that a single column can drive an aesthetic like color. With counts_long already in long format, we can map the condition column directly to point color.

gene_levels <- counts_wide$gene
conditions <- names(counts_wide)[-1]
colors <- palette.colors(7)[c(7,4,3)]
names(colors) <- conditions

plot(
  as.numeric(factor(counts_long$gene, levels = gene_levels)),
  counts_long$count,
  col = colors[counts_long$condition],
  pch = 16,
  cex = 1.5,
  xaxt = "n",
  xlab = "gene",
  ylab = "count",
  main = "counts by gene and condition"
)
axis(1, at = 1:3, labels = gene_levels)
legend(
  "topright", legend = conditions, col = colors, 
  pch = 16, y.intersp = 1.5, inset = 0.02
  )

The key line is col = colors[counts_long$condition]: because condition is a character vector of condition names, and colors is a named vector, this indexes colors by name to produce a color for every row.

Plotting tidy data with ggplot2

The same plot requires far less code with ggplot2, and many things that take manual effort in base R (like the legend) appear automatically as sensible defaults.

library(ggplot2)

ggplot(
  data = counts_long, 
  mapping = aes(x = gene, y = count, color = condition)
) +
  geom_point(size = 3)

The color = condition mapping is all it takes: ggplot2 assigns colors and adds a guide (legend) automatically. In base R we had to build a named color vector and call legend() by hand.

ggplot2 is built around a layered grammar of graphics. Each + adds a new layer or modification (a geometry like geom_point, a scale, a theme, a facet), so you can start with a minimal plot and layer complexity on top without rewriting anything.

NoteExercise

Extend the previous plot by adding lines to connect the points within each condition, and facet the plot by condition (make three plots, one for each condition).

TipHint

Lines can be drawn by adding geom_line() but requires telling the mapping (in aes) to know how to group the observations. Faceting can be accomplished with facet_wrap(~ variable_to_facet_on).

ggplot(
  data = counts_long,
  mapping = aes(
    x = gene, y = count,
    color = condition, group = ______
  )
) +
  geom_point(size = 3) +
  ______ +
  ______

Solution (Solution).

ggplot(
  data = counts_long,
  mapping = aes(
    x = gene, y = count,
    color = condition, group = condition
  )
) +
  geom_point(size = 3) +
  geom_line() +
  facet_wrap(~ condition)

gradethis::grade_this_code()

More sophisticated plots can be built by layering additional geoms: geom_line() connects points, geom_smooth() fits and draws a curve through the data, and stat_summary() computes statistics (such as means and standard errors) and renders them as geoms directly within the plotting call, without needing a separate dplyr summarization step. For a complete reference of available geometries, statistics, scales, and themes, see the ggplot2 documentation.

Tips for making data findable, accessible, interoperable, and re-usable

Above we considered the choice of tidy data vs other formats of storing or presenting data. Here we describe some other tips that help make data more easy to use.

A well-known set of guiding principles of working with and sharing data is called FAIR, an acronym representing the words in the section header here.

Some tips for making data both tidy and FAIR are as follows:

• Missing values are data. While typically represented with NA in R or NaN in python, one can model missingness itself, converting to a binary “shadow” matrix. See Tidy Missing Data in the naniar package vignette. Avoid -1 for missing or flagged as this can lead to downstream confusion.
• Dates, units, and categories should be explicit and consistent (no “10mg” mixed with “10”). Units can be provided in the column name (dosage_mg). It’s a good idea to provide an accompanying data dictionary defining units and categories. Use standard formats for dates such as ISO 8601: YYYY-MM-DD. For working in R, see the lubridate package.
• One thing per cell. Avoid embedding multiple values (e.g., “120/80” for blood pressure). Discourage use of highlighting, color, or style in place of data; instead use an additional column taking boolean values, e.g. flagged or highlight.
• Column names should always be machine-readable: no spaces, no special characters, unambiguous. Replace dashes with underscore. When not sure what to name a column, try to find a previously published dataset and re-use their naming scheme rather than invent a new one.
• Raw data is sacred: never overwrite it. Keep a clean separation between raw and processed data. Raw data can be excluded from version control, but should be described including provenance in an accompanying README.md. Where did this data come from, how was it collected? Ideally provide a script for downloading from stable, remote sources such as Zenodo.
• Finally, a stranger, or a future “you”, should be able to re-run everything. Comments, links, software versions (see “Session info” below), can be of great value for future data and analysis re-use. If using randomness as part of analysis (e.g. bootstrapping, sub-sampling, etc.) always set a random seed. For simulations, keep track of the seeds used to generate the multiple simulated datasets.
• Data should be findable: when depositing data in a repository (Zenodo, Figshare, GEO, etc.), always attach structured metadata describing what the dataset contains, how it was collected, and who produced it. A little effort goes a long way. Ask someone who works in your field if there are additional pieces of information they would find valuable. Where possible, use standard ontologies (e.g. the Gene Ontology, the Cell Ontology, or SNOMED CT for clinical data) to describe variables and sample attributes. Ontology terms are machine-readable identifiers that make metadata computable: other researchers and automated tools can find, filter, and integrate your data without having to interpret free-text descriptions.

Tidy data in Python

The tidyverse has counterparts in Python. pandas is the most widely used tabular data library and covers all of the dplyr verbs and pivot operations. polars is a newer alternative with a Rust backend and a lazy query engine, offering significant speed improvements on large datasets and a chaining syntax closer to dplyr pipelines. siuba ports dplyr’s verbs directly into Python using the >> pipe operator, so counts_long >> filter(_.count > 100) reads almost exactly like the R version.

One practical aspect of interoperability is file format. CSV is universal but loses column types. The parquet format stores typed columnar data efficiently and is readable by R (via the arrow package), Python (via pandas or pyarrow), and many other tools, making it a good choice when a dataset crosses language boundaries or needs long-term archiving. The following example writes counts_long from R, then reads it into Python and reproduces the grouped summary from above.

library(arrow)
write_parquet(counts_long, "counts_long.parquet")

import pandas as pd

counts = pd.read_parquet("counts_long.parquet")

# mirror: counts_long |> filter(condition != "control") |>
#   group_by(condition) |>
#   summarize(total = sum(count), mean_count = mean(count))
(
    counts
    .query('condition != "control"')
    .groupby("condition")
    .agg(total=("count", "sum"), mean_count=("count", "mean"))
    .reset_index()
)

Because parquet preserves column types (gene and condition remain strings, count remains an integer), no manual casting is needed after reading. Parquet also compresses well and can be queried column-selectively without reading the full file, which matters as datasets grow large.

Exercise: preserving raw data

A common operation in genomics is to center (and potentially scale) each gene’s expression values by subtracting the gene’s mean across conditions (perhaps after taking the logarithm). One tempting approach works directly on a matrix and reassigns the result:

counts_mat <- as.matrix(counts_wide[,-1])
rownames(counts_mat) <- counts_wide[,1,drop=TRUE]
counts_mat <- log10(counts_mat + 1)
cat("geometric mean:\n")

geometric mean:

10^(rowMeans(counts_mat))

    GAPDH       IL6      IRF1 
527.08352  96.74667 166.49755 

counts_mat <- counts_mat - rowMeans(counts_mat)
cat("\n")

cat("log-transformed and centered data:\n")

log-transformed and centered data:

counts_mat

           control        IFNg  Salmonella
GAPDH -0.005041714  0.02708342 -0.02204171
IL6   -0.494274346 -0.02659465  0.52086899
IRF1  -0.513837675  0.24248514  0.27135254

This works numerically, but the reassignment destroys information. The values in counts_mat are no longer counts; they are deviations of original counts, yet the variable name gives no indication of what has been done. More importantly, the original information about the scale of the counts per gene is gone: we can no longer recover that GAPDH had ~530 counts, while IL6 had ~100 and IRF1 had ~160, in terms of the geometric mean.

The tidy alternative is to leave counts_long untouched and add new columns with mutate(), keeping the original count alongside intermediate values. A reader (or a future you) can immediately see both the absolute level and the deviation, and the operation is self-documenting in the column names.

NoteExercise

Starting from counts_long, add two new columns without overwriting anything:

• log_count: log10 of the count plus one
• ctr_log10_count: log_count centered by subtracting the per-gene mean

counts_long |>
  ______(log_count = ______) |>
  ______(______) |>
  ______(ctr_log10_count = ______)

Solution (Solution).

counts_long |>
  mutate(log_count = log10(count + 1)) |>
  group_by(gene) |>
  mutate(ctr_log10_count = log_count - mean(log_count))

# A tibble: 9 × 5
# Groups:   gene [3]
  gene  condition  count log_count ctr_log10_count
  <chr> <chr>      <dbl>     <dbl>           <dbl>
1 GAPDH control      520      2.72        -0.00504
2 GAPDH IFNg         560      2.75         0.0271 
3 GAPDH Salmonella   500      2.70        -0.0220 
4 IL6   control       30      1.49        -0.494  
5 IL6   IFNg          90      1.96        -0.0266 
6 IL6   Salmonella   320      2.51         0.521  
7 IRF1  control       50      1.71        -0.514  
8 IRF1  IFNg         290      2.46         0.242  
9 IRF1  Salmonella   310      2.49         0.271  

gradethis::grade_this_code()

Exercise: modeling missingness

The FAIR tips above note that missing values are data and that missingness itself can be modeled. Here we construct a dataset where missingness depends on a fully observed covariate, then explore whether we can recover that structure from the data.

We start from the complete cases of the penguins dataset and simulate a realistic mechanism: suppose penguins with flippers longer than 215 mm often cannot fit through the door to the room where the scale is located, so their body mass is unrecorded 75% of the time. Because missingness depends on flipper_length_mm, which is always observed, this is an example of missing at random (MAR) given the covariates. A companion indicator column records whether each body mass value is missing.

library(palmerpenguins)

set.seed(42)
penguins_complete <- penguins |> drop_na()

penguins_obs <- penguins_complete |>
  mutate(
    body_mass_missing = flipper_length_mm >= 215 & runif(n()) < 0.75,
    body_mass_obs = if_else(body_mass_missing, NA, body_mass_g)
  )

names(penguins_obs)

 [1] "species"           "island"            "bill_length_mm"    "bill_depth_mm"    
 [5] "flipper_length_mm" "body_mass_g"       "sex"               "year"             
 [9] "body_mass_missing" "body_mass_obs"    

Note the difference in flipper length between penguins with and without a recorded body mass. To show this we introduce slice_sample(), which allows for random sampling within groups.

set.seed(5) # to fix the random sampling
penguins_obs |>
  select(species, flipper_length_mm, body_mass_g, body_mass_obs, body_mass_missing) |>
  group_by(body_mass_missing) |>
  slice_sample(n = 5)

# A tibble: 10 × 5
# Groups:   body_mass_missing [2]
   species   flipper_length_mm body_mass_g body_mass_obs body_mass_missing
   <fct>                 <int>       <int>         <int> <lgl>            
 1 Gentoo                  230        5800          5800 FALSE            
 2 Chinstrap               196        3900          3900 FALSE            
 3 Gentoo                  214        4925          4925 FALSE            
 4 Chinstrap               197        3750          3750 FALSE            
 5 Chinstrap               197        3300          3300 FALSE            
 6 Gentoo                  222        5350            NA TRUE             
 7 Gentoo                  230        5500            NA TRUE             
 8 Gentoo                  217        4900            NA TRUE             
 9 Gentoo                  219        5200            NA TRUE             
10 Gentoo                  224        5350            NA TRUE             

Because missingness depends on flipper length, which we always observe, we can model it directly. Flipper length is also highly correlated with body mass, such that we will tend to not observe heavier penguins. Let’s confirm the correlation using our full dataset:

penguins_obs |>
  ggplot(aes(flipper_length_mm, body_mass_g)) + 
  geom_point()

Again, this correlation implies we can model the binary missingness indicator as a function of flipper length using logistic regression. Logistic regression is implemented in R via glm() with family = binomial, or equivalently family = "binomial" (quotes not needed as binomial is an object in base R, a special function called a family). Technically the model is a Bernoulli regression, a special case of the binomial where the number of trials is 1, and the response must be in {0, 1} (or equivalently TRUE/FALSE).

Like lm(), a logistic regression fits a linear combination of predictor variables (here, to the log-odds of the response being 1), but instead of minimizing squared residuals, it estimates coefficients by maximizing the likelihood of the observed 0/1 outcomes under the model. The fitted log-odds can be converted to a probability with type = "response" in predict(). The question we are asking: does a longer flipper predict a higher probability of the body mass being missing?

NoteExercise

Fit a logistic regression predicting body_mass_missing from flipper_length_mm in penguins_obs. You will then see the estimated coefficients with the tidy() function in the broom package.

TipHint

Modeling one variable on another uses the formula syntax y ~ x, where y is the response variable and x is the predictor variable.

fit <- glm(
  ______ ~ ______,
  data = penguins_obs,
  family = ______
)
broom::tidy(fit)

Solution (Solution).

fit <- glm(
  body_mass_missing ~ flipper_length_mm,
  data = penguins_obs,
  family = binomial
)
broom::tidy(fit)

# A tibble: 2 × 5
  term              estimate std.error statistic  p.value
  <chr>                <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)        -53.4      7.60       -7.04 1.97e-12
2 flipper_length_mm    0.246    0.0352      6.99 2.80e-12

gradethis::grade_this_code()

The coefficient on flipper_length_mm is positive and significant: penguins with longer flippers are indeed more likely to have their body mass missing, which recovers the MAR mechanism we built in.

In practice, missingness models are rarely this simple. A more realistic model would include multiple predictors, and the relationship between a predictor and the log-odds of missingness may not be exactly linear; semi-parametric approaches such as splines allow the shape of the relationship to be estimated from the data rather than assumed. Beyond predicting missingness, if the goal is causal estimation (understanding how one variable affects another while accounting for confounding), it becomes important to diagram the relationships among all variables in a directed acyclic graph (DAG). DAGs make explicit which variables are confounders, mediators, or colliders, and guide which variables should or should not be conditioned on. A comprehensive and freely available introduction to this framework is What If by Hernán and Robins (miguelhernan.org/whatifbook).

Now let’s attach the model’s predicted missingness probability back to each row. First we will run the fitting code from the exercise again so that the fitted model is available.

fit <- glm(
  body_mass_missing ~ flipper_length_mm,
  data = penguins_obs,
  family = binomial
)
broom::glance(fit) # statistical model summaries

# A tibble: 1 × 8
  null.deviance df.null logLik   AIC   BIC deviance df.residual  nobs
          <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int> <int>
1          314.     332  -68.0  140.  148.     136.         331   333

NoteExercise

Add a column prob_missing to penguins_obs containing the predicted probability of missingness from fit, then assign the result to the dataset penguins_modeled.

TipHint

The predict() function generates predicted values from a fitted model. For a logistic regression, use type = "response" to return predicted probabilities rather than log-odds.

penguins_modeled <- penguins_obs |>
  mutate(______ = ______(fit, type = ______))

Solution (Solution).

penguins_modeled <- penguins_obs |>
  mutate(prob_missing = predict(fit, type = "response"))

gradethis::grade_this_code()

One application of a fitted missingness model is inverse probability weighting (IPW). The intuition is as follows: penguins with long flippers are under-represented in our observed body mass values, so we compensate by giving those observed penguins extra weight in any downstream analysis (here “weight” refers to the influence of a datum on the estimate, not to the weight of the penguin!).

Specifically, each observed value receives a weight of 1 / (1 - prob_missing): a penguin with a 75% predicted probability of being missing, but whose value we do happen to observe, is up-weighted in the equations used to estimate parameters by a factor of 4, because it is standing in for roughly four penguins like itself.

Because our mechanism is MAR (missingness depends only on flipper length, which is fully observed), this re-weighting can recover unbiased estimates of quantities computed on the full data. If the mechanism were not at random (MNAR), for example if very heavy penguins were missing regardless of flipper length, but due to unobserved values including the missing variable itself, our estimates of the missingness probability would not be sufficient to recover unbiased estimates.

Now let’s plot our new variable prob_missing against true body mass, with the point shape indicating whether the value was missing (empty circle) or not (filled circle). This will let us see how well the model separates the two groups. Generally this plot wouldn’t be possible because the missing values would be missing!

Again, we need to add prob_missing from the exercise, for use below.

penguins_modeled <- penguins_obs |>
  mutate(prob_missing = predict(fit, type = "response"))

ggplot(
  data = penguins_modeled,
  mapping = aes(x = body_mass_g, y = prob_missing, shape = body_mass_missing)
) +
  geom_point(size = 2) +
  scale_shape_manual(values = c("FALSE" = 16, "TRUE" = 1)) +
  labs(x = "body mass (g)", y = "predicted P(missing)", shape = "missing")

NoteQuestion

The open circles (missing values) are concentrated at high body mass, but the predicted probability on the y-axis is not a simple increasing function of body mass. Why is this so?

Concept review

• Tidy and untidy data are both useful, depending on the task. Wide (matrix) format is natural for linear algebra and many statistical model inputs; long (tidy) format is natural for arbitrary filtering, grouping, and plotting. pivot_longer() and pivot_wider() let you move between them.

• Tidy data makes data manipulation and visualization easier because each variable occupies exactly one column, so a single column can directly drive a grouping, a color aesthetic, or a facet without any reshaping inside the plotting call.

• Good data practices make data useful to others (and to your future self): keep raw data untouched, use machine-readable column names, record units and provenance, set random seeds when necessary, log software versions, and use structured metadata with standard ontologies when depositing data publicly.

• dplyr and ggplot2 are powerful because of their composability. Both are built around a small set of simple, well-named operations that chain together cleanly. A filter() followed by a group_by() and a summarize(), or a geom_point() followed by a facet_wrap(), reads almost like a description of what you want. Complex operations built from simple pieces tend to stay readable even as they grow.

• The same ideas carry across languages. pandas in Python provides close equivalents to nearly every dplyr verb, and both benefit from the same underlying conceptual framework: apply composable, single-purpose operations to tidy data. Learning this way of data analytic thinking, and this way of structuring datasets in one language makes it much easier to pick up the other.

• Missingness is data. Modeling the probability that a value is missing, and understanding whether that probability depends on observed or unobserved quantities (MAR vs. MNAR), is an important part of any honest analysis.

Session info

sessionInfo() prints the R version, operating system, locale, and the name and version of every package loaded in the current session. Including it at the end of a script or notebook is a lightweight form of provenance: if a result changes or a bug appears months later, the session info is often the first place to look; it tells you exactly what software produced the output. It is also essential for reproducibility reports and as supplementary information in publications.

sessionInfo()

R version 4.6.0 (2026-04-24)
Platform: aarch64-apple-darwin23
Running under: macOS Tahoe 26.3.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Uzhgorod
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
[1] ggplot2_4.0.3        palmerpenguins_0.1.1 dplyr_1.2.1          tidyr_1.3.2         

loaded via a namespace (and not attached):
 [1] gtable_0.3.6       jsonlite_2.0.0     compiler_4.6.0     tidyselect_1.2.1   stringr_1.6.0     
 [6] scales_1.4.0       yaml_2.3.12        fastmap_1.2.0      R6_2.6.1           labeling_0.4.3    
[11] generics_0.1.4     knitr_1.51         htmlwidgets_1.6.4  backports_1.5.1    tibble_3.3.1      
[16] pillar_1.11.1      RColorBrewer_1.1-3 rlang_1.2.0        utf8_1.2.6         stringi_1.8.7     
[21] broom_1.0.13       xfun_0.58          S7_0.2.2           otel_0.2.0         cli_3.6.6         
[26] withr_3.0.2        magrittr_2.0.5     digest_0.6.39      grid_4.6.0         lifecycle_1.0.5   
[31] vctrs_0.7.3        evaluate_1.0.5     glue_1.8.1         farver_2.1.2       stats4_4.6.0      
[36] rmarkdown_2.31     purrr_1.2.2        tools_4.6.0        pkgconfig_2.0.3    htmltools_0.5.9   

A popular alternative is sessionInfo() from the sessioninfo package, which formats packages in a cleaner table and records the source of each package: CRAN, Bioconductor, or a GitHub repository with its commit SHA, which matters when packages installed from GitHub are not on a numbered release. devtools depends on sessioninfo and re-exports the same function as devtools::session_info(), so if devtools is already loaded you can call it either way.

Use of generative AI

Portions of this tutorial were developed with the assistance of Claude Code (Anthropic).