Handling Missing Data in Tabular Datasets

A practical guide with a clinical dataset

1 Introduction

Missing data is one of the most common — and most under-addressed — problems in applied biostatistics and epidemiology. Naively dropping incomplete rows or filling gaps with the column mean can silently bias effect estimates, shrink variance, and distort standard errors. This guide walks through a principled workflow for diagnosing, visualizing, and handling missing values, using a simulated clinical dataset (missing_data_demo.csv) that was constructed with known missingness mechanisms, so we can check whether our diagnostics correctly recover them.

The dataset contains 300 simulated patients with the following variables:

Variable	Description	Missingness mechanism injected
patient_id	Unique identifier	None
age	Age in years	None
sex	Female / Male	None
bmi	Body mass index	MNAR — depends on BMI itself
smoking_status	Never / Former / Current	MCAR
systolic_bp	Systolic blood pressure	MAR — depends on sex
hba1c	Glycated haemoglobin (%)	MAR — depends on age
income_usd	Monthly income (USD)	MNAR — depends on income itself
outcome_event	Binary outcome (e.g. disease event)	None

2 Setup

library(tidyverse)   # data wrangling & viz
library(naniar)      # missing data visualization & summaries
library(VIM)         # missingness patterns (aggr, matrixplot)
library(mice)         # multiple imputation by chained equations
library(gtsummary)    # summary tables
library(missForest)   # random-forest based imputation (optional comparison)

df <- read_csv("missing_data_demo.csv") |>
  mutate(
    sex = factor(sex),
    smoking_status = factor(smoking_status, levels = c("Never", "Former", "Current")),
    outcome_event = factor(outcome_event, labels = c("No event", "Event"))
  )

glimpse(df)

Rows: 300
Columns: 9
$ patient_id     <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
$ age            <dbl> 49.3, 30.4, 55.5, 58.2, 18.0, 26.8, 46.8, 40.6, 44.8, 3…
$ sex            <fct> Female, Female, Female, Female, Female, Male, Female, M…
$ bmi            <dbl> 35.1, 20.9, 28.2, 36.1, 18.7, 20.6, 24.9, 24.4, 31.1, 2…
$ smoking_status <fct> Current, NA, Never, Never, Never, Never, Current, Forme…
$ systolic_bp    <dbl> 147, 132, 159, 169, 125, 145, 156, NA, 158, 158, NA, 14…
$ hba1c          <dbl> 6.76, 6.68, 6.20, 7.67, 5.90, 5.70, 6.67, 6.16, 6.31, 6…
$ income_usd     <dbl> 1849, 1595, 855, 1613, NA, 436, 1374, NA, 1196, 415, 14…
$ outcome_event  <fct> No event, No event, No event, No event, No event, No ev…

3 Step 1: Quantify missingness

Always start with a simple count before doing anything else.

miss_var_summary(df)

# A tibble: 9 × 3
  variable       n_miss pct_miss
  <chr>           <int>    <num>
1 income_usd         82    27.3 
2 systolic_bp        38    12.7 
3 hba1c              36    12   
4 smoking_status     30    10   
5 bmi                17     5.67
6 patient_id          0     0   
7 age                 0     0   
8 sex                 0     0   
9 outcome_event       0     0   

pct_complete_case(df)   # % of rows with zero missing values

[1] 47.66667

n_var_miss(df)          # how many variables have any missingness

[1] 5

4 Step 2: Visualize missingness patterns

4.1 Overall missingness map

vis_miss(df, sort_miss = TRUE)

Figure 1: Missingness map across all variables

4.2 Combination patterns (which variables tend to be missing together)

gg_miss_upset(df)

Figure 2: Intersections of missingness across variables

4.3 VIM aggregation plot

aggr(df, numbers = TRUE, prop = FALSE, sortVars = TRUE,
     labels = names(df), cex.axis = 0.6, gap = 2)

 Variables sorted by number of missings: 
       Variable Count
     income_usd    82
    systolic_bp    38
          hba1c    36
 smoking_status    30
            bmi    17
     patient_id     0
            age     0
            sex     0
  outcome_event     0

Figure 3: Proportion and combinations of missing values (VIM::aggr)

5 Step 3: Investigate the missingness mechanism

This is the step people skip — but it determines which handling method is valid.

• MCAR (Missing Completely At Random): missingness unrelated to any observed or unobserved data.
• MAR (Missing At Random): missingness depends on observed variables (e.g., older patients less likely to have hba1c recorded).
• MNAR (Missing Not At Random): missingness depends on the unobserved value itself (e.g., higher earners less likely to report income_usd).

5.1 Visual check: does missingness in one variable relate to another observed variable?

ggplot(df, aes(x = age, y = is.na(hba1c))) +
  geom_jitter(height = 0.05, alpha = 0.4) +
  labs(y = "hba1c is missing", x = "Age",
       title = "Older patients more likely to be missing hba1c") +
  theme_minimal()

Figure 4: hba1c missingness by age — a MAR-style pattern

df |>
  group_by(sex) |>
  summarise(pct_missing_sbp = mean(is.na(systolic_bp)) * 100)

# A tibble: 2 × 2
  sex    pct_missing_sbp
  <fct>            <dbl>
1 Female            7.14
2 Male             18.5 

5.2 Little’s MCAR test

A formal (if imperfect) test of the MCAR assumption:

mcar_test(df |> select(-patient_id, -outcome_event))

# A tibble: 1 × 4
  statistic    df p.value missing.patterns
      <dbl> <dbl>   <dbl>            <int>
1      110.    87  0.0461               18

A significant p-value suggests the data are not MCAR — consistent with how we simulated hba1c, systolic_bp, and income_usd. Note this test cannot distinguish MAR from MNAR; that requires subject-matter judgment about whether the missingness plausibly depends on the unobserved value itself.

6 Step 4: Naive approaches (and why to be cautious)

6.1 Listwise deletion (complete-case analysis)

df_cc <- df |> drop_na()
nrow(df_cc)  # how many rows survive

[1] 143

Simple, but valid only under MCAR. Here we lose a large share of rows, and since some missingness depends on age, sex, and income_usd, complete-case analysis will bias downstream estimates toward whichever subgroup has more complete records.

6.2 Simple imputation (mean / median / mode)

df_simple <- df |>
  mutate(
    bmi = if_else(is.na(bmi), median(bmi, na.rm = TRUE), bmi),
    hba1c = if_else(is.na(hba1c), mean(hba1c, na.rm = TRUE), hba1c),
    systolic_bp = if_else(is.na(systolic_bp), mean(systolic_bp, na.rm = TRUE), systolic_bp)
  )

Fast, but artificially shrinks variance and attenuates correlations between variables — avoid this for anything beyond quick sanity checks.

6.3 Missingness indicators

Sometimes the fact that a value is missing is itself informative (as with income_usd and bmi here). Adding flag columns preserves that signal for prediction models:

df_flagged <- df |>
  mutate(
    bmi_missing = is.na(bmi),
    income_missing = is.na(income_usd)
  )

7 Step 5: Multiple imputation with mice

Multiple imputation is generally the right default for inferential work (regression coefficients, confidence intervals), because it propagates the uncertainty due to missingness rather than pretending imputed values are known.

# Quick view of the default method assigned per variable
init <- mice(df |> select(-patient_id), maxit = 0)
init$method

           age            sex            bmi smoking_status    systolic_bp 
            ""             ""          "pmm"      "polyreg"          "pmm" 
         hba1c     income_usd  outcome_event 
         "pmm"          "pmm"             "" 

imp <- mice(
  df |> select(-patient_id),
  m = 5,            # number of imputed datasets
  maxit = 20,        # iterations per dataset
  method = "pmm",    # predictive mean matching — good default for skewed/continuous vars
  seed = 123
)

7.1 Diagnostics: did the chains converge?

plot(imp)

Figure 5: Convergence trace plots

Figure 6: Convergence trace plots

7.2 Compare imputed vs observed distributions

densityplot(imp, ~ hba1c + bmi + systolic_bp + income_usd)

Figure 7: Density of observed (blue) vs imputed (red) values

7.3 Fit a model across the imputed datasets and pool results

fit <- with(imp, glm(outcome_event ~ age + sex + bmi + hba1c + systolic_bp,
                      family = binomial))
pooled <- pool(fit)
summary(pooled)

         term     estimate  std.error  statistic        df      p.value
1 (Intercept) -13.02464634 3.60180819 -3.6161410 207.84956 0.0003753052
2         age   0.04267201 0.02028031  2.1041098 210.34837 0.0365566848
3     sexMale   0.44613121 0.39141813  1.1397817 291.04710 0.2553144097
4         bmi   0.07310228 0.04833025  1.5125576  75.99135 0.1345399718
5       hba1c   0.50032914 0.34145132  1.4653016 228.30658 0.1442144441
6 systolic_bp   0.02177220 0.02301983  0.9458022 112.85518 0.3462701316

The pooled standard errors correctly reflect imputation uncertainty — this is the key advantage over single (deterministic) imputation.

8 Step 6: Alternative — random forest imputation (missForest)

Useful when relationships are non-linear and you don’t need inferential (p-value) validity, e.g. for a prediction pipeline.

df_for_rf <- df |> select(-patient_id) |> as.data.frame()
rf_imp <- missForest(df_for_rf, verbose = FALSE)
rf_imp$OOBerror   # out-of-bag imputation error estimate

    NRMSE       PFC 
0.3770512 0.1814815 

9 Step 7: Comparing approaches

tibble(
  method = c("Complete-case", "Mean imputation", "MICE (pooled, m=5)"),
  mean_hba1c = c(
    mean(df_cc$hba1c, na.rm = TRUE),
    mean(df_simple$hba1c, na.rm = TRUE),
    mean(complete(imp, action = "long")$hba1c, na.rm = TRUE)
  )
)

?(caption)

# A tibble: 3 × 2
  method             mean_hba1c
  <chr>                   <dbl>
1 Complete-case            6.66
2 Mean imputation          6.64
3 MICE (pooled, m=5)       6.64

Compare these to the true pre-missingness mean if you’re working with a simulation like this one — it’s a good way to sanity-check that your imputation model isn’t introducing new bias.

10 Practical recommendations

• Report the amount and pattern of missingness for every variable in any writeup — not just the analytic sample size.
• Never use single mean/mode imputation for inferential (as opposed to purely predictive) analyses.
• Default to MICE (or another proper multiple-imputation method) when you need valid standard errors and confidence intervals.
• MNAR is the hard case. No automatic imputation method fully solves it — you need a sensitivity analysis (e.g. pattern-mixture models, or mice’s delta-adjustment) and domain reasoning about why values are missing.
• Missingness indicators are cheap insurance in prediction pipelines, especially when missingness itself carries signal (as with income here).
• Always visualize before and after imputation to catch implausible values or distributional distortion.

11 Session info

sessionInfo()

R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 26100)

Matrix products: default


locale:
[1] LC_COLLATE=English_Kenya.utf8  LC_CTYPE=English_Kenya.utf8   
[3] LC_MONETARY=English_Kenya.utf8 LC_NUMERIC=C                  
[5] LC_TIME=English_Kenya.utf8    

time zone: Africa/Nairobi
tzcode source: internal

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

other attached packages:
 [1] missForest_1.6.1 gtsummary_2.5.0  mice_3.19.0      VIM_6.2.2       
 [5] colorspace_2.1-2 naniar_1.1.0     lubridate_1.9.5  forcats_1.0.1   
 [9] stringr_1.5.2    dplyr_1.2.0      purrr_1.2.2      readr_2.1.5     
[13] tidyr_1.3.1      tibble_3.3.0     ggplot2_4.0.3    tidyverse_2.0.0 

loaded via a namespace (and not attached):
 [1] Rdpack_2.6.6         gridExtra_2.3        rlang_1.1.7         
 [4] magrittr_2.0.4       e1071_1.7-17         compiler_4.3.1      
 [7] vctrs_0.7.2          pkgconfig_2.0.3      shape_1.4.6.1       
[10] crayon_1.5.3         fastmap_1.2.0        backports_1.5.1     
[13] labeling_0.4.3       utf8_1.2.6           rmarkdown_2.30      
[16] tzdb_0.5.0           nloptr_2.2.1         itertools_0.1-3     
[19] UpSetR_1.4.1         visdat_0.6.0         bit_4.6.0           
[22] xfun_0.53            glmnet_5.0           jomo_2.7-6          
[25] randomForest_4.7-1.2 jsonlite_2.0.0       pan_1.9             
[28] broom_1.0.12         parallel_4.3.1       R6_2.6.1            
[31] stringi_1.8.7        vcd_1.4-13           RColorBrewer_1.1-3  
[34] ranger_0.18.0        car_3.1-5            boot_1.3-28.1       
[37] rpart_4.1.19         lmtest_0.9-40        Rcpp_1.1.1          
[40] iterators_1.0.14     knitr_1.50           zoo_1.8-15          
[43] Matrix_1.5-4.1       splines_4.3.1        nnet_7.3-19         
[46] timechange_0.4.0     tidyselect_1.2.1     rstudioapi_0.18.0   
[49] abind_1.4-8          yaml_2.3.10          codetools_0.2-19    
[52] doRNG_1.8.6.3        plyr_1.8.9           lattice_0.21-8      
[55] withr_3.0.2          S7_0.2.0             evaluate_1.0.5      
[58] survival_3.5-5       proxy_0.4-29         norm_1.0-11.1       
[61] pillar_1.11.1        carData_3.0-6        rngtools_1.5.2      
[64] foreach_1.5.2        reformulas_0.4.4     generics_0.1.4      
[67] vroom_1.6.6          sp_2.2-1             hms_1.1.4           
[70] scales_1.4.0         laeken_0.5.3         minqa_1.2.8         
[73] class_7.3-22         glue_1.8.0           tools_4.3.1         
[76] robustbase_0.99-7    data.table_1.18.4    lme4_2.0-1          
[79] rbibutils_2.4.1      nlme_3.1-162         Formula_1.2-5       
[82] cli_3.6.5            gtable_0.3.6         DEoptimR_1.1-4      
[85] digest_0.6.37        htmlwidgets_1.6.4    farver_2.1.2        
[88] htmltools_0.5.8.1    lifecycle_1.0.5      mitml_0.4-5         
[91] bit64_4.8.2          MASS_7.3-60