diff --git a/DESCRIPTION b/DESCRIPTION
index f9745be7..9b7ead7a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: faux
 Title: Simulation for Factorial Designs
-Version: 1.2.3.9000
-Date: 2025-09-25
+Version: 1.2.3.9001
+Date: 2025-10-20
 Authors@R: c(
     person(
         given = "Lisa", 
diff --git a/NEWS.md b/NEWS.md
index 94c6d449..0e91522c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,7 @@
 # faux (development version)
 
+* Data frames created with `sim_design()` now have column labels that can be used by ggplot2 4.0 for figure labels.
+
 # faux 1.2.3 (2025-09-25)
 
 * Fixed failing tests due to ggplot2 4.0 ggplot object structure changes
diff --git a/R/long2wide.R b/R/long2wide.R
index b691a579..c688c65e 100644
--- a/R/long2wide.R
+++ b/R/long2wide.R
@@ -58,6 +58,7 @@ long2wide <- function(data, within = c(), between = c(), dv = "y", id = "id", se
   }
   class(d1) <- c("faux", "data.frame")
   rownames(d1) <- NULL
+  d1 <- add_labels(d1)
   
   d1
 }
diff --git a/R/sim_design.R b/R/sim_design.R
index 445f5856..e6214b1c 100644
--- a/R/sim_design.R
+++ b/R/sim_design.R
@@ -213,6 +213,7 @@ sim_data <- function(design, empirical = FALSE,
   }
   
   class(df_return) <- c("faux", "data.frame")
+  rownames(df_return) <- c() # get rid of row names
   
   if (rep == 1) {
     df_return$.rep. <- NULL
@@ -220,7 +221,7 @@ sim_data <- function(design, empirical = FALSE,
     # nest reps
     df_rep <- by(df_return, df_return$.rep., function(x) {
       x$.rep. <- NULL
-      x
+      add_labels(x, design)
     })
     df_return <- data.frame(rep = 1:rep)
     df_return$data <- df_rep # can't assign list in data.frame
@@ -230,6 +231,68 @@ sim_data <- function(design, empirical = FALSE,
     df_return <- df_return[order(df_return$rep),] 
   }
   
-  rownames(df_return) <- c() # get rid of row names
+  df_return <- add_labels(df_return, design) # add column labels
+  
   return(df_return)
 }
+
+
+#' Add Labels to Data Tables
+#'
+#' @param data a data table
+#' @param design a design specification (if NULL, is read from data)
+#'
+#' @returns a data frame with labelled columns
+#' @keywords internal
+add_labels <- function(data, design = NULL) {
+  if ("design" %in% names(attributes(data))) {
+    # get parameters from design
+    design <- get_design(data)
+  } 
+  
+  if (is.null(design)) return(data)
+  
+  # nested
+  if ("rep" %in% names(data)) 
+    attr(data$rep, "label") <- "replicate index"
+  
+  if ("data" %in% names(data)) {
+    data$data <- lapply(data$data, add_labels, design = design)
+    attr(data$data, "label") <- "data"
+    
+    return(data)
+  }
+  
+  # id column
+  id <- names(design$id)
+  attr(data[[id]], "label") <- design$id[[1]]
+  
+  # dv / within
+  dv <- names(design$dv)
+  if (dv %in% names(data)) {
+    attr(data[[dv]], "label") <- design$dv[[1]]
+    
+    # long within
+    for (w in names(design$within)) {
+      attr(data[[w]], "label") <- design$vardesc[[w]]
+    }
+  } else if (length(design$within)) {  # wide within
+    exp <- expand.grid(rev(design$within))
+    w_labels <- apply(exp, 1, function(x) { 
+      paste(rev(x), collapse = ":") 
+    })
+    within_cols <- cell_combos(design$within, sep = design$sep) 
+    names(w_labels) <- within_cols
+    
+    for (nm in within_cols) {
+      attr(data[[nm]], "label") <- w_labels[[nm]]
+    }
+  }
+  
+  # between
+  for (b in names(design$between)) {
+    attr(data[[b]], "label") <- design$vardesc[[b]]
+  }
+  
+  return(data)
+}
diff --git a/R/wide2long.R b/R/wide2long.R
index b1c8d192..d7090063 100644
--- a/R/wide2long.R
+++ b/R/wide2long.R
@@ -65,6 +65,7 @@ wide2long <- function(data, within_factors = c(), within_cols = c(),
   attributes(longdat)$design <- design
   class(longdat) <- c("faux", "data.frame")
   rownames(longdat) <- NULL
+  longdat <- add_labels(longdat)
   
   longdat
 }
diff --git a/man/add_labels.Rd b/man/add_labels.Rd
new file mode 100644
index 00000000..ff0de44a
--- /dev/null
+++ b/man/add_labels.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sim_design.R
+\name{add_labels}
+\alias{add_labels}
+\title{Add Labels to Data Tables}
+\usage{
+add_labels(data, design = NULL)
+}
+\arguments{
+\item{data}{a data table}
+
+\item{design}{a design specification (if NULL, is read from data)}
+}
+\value{
+a data frame with labelled columns
+}
+\description{
+Add Labels to Data Tables
+}
+\keyword{internal}
diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf
index f78d8c39..4644dc08 100644
Binary files a/tests/testthat/Rplots.pdf and b/tests/testthat/Rplots.pdf differ
diff --git a/tests/testthat/test-add_labels.R b/tests/testthat/test-add_labels.R
new file mode 100644
index 00000000..664763ab
--- /dev/null
+++ b/tests/testthat/test-add_labels.R
@@ -0,0 +1,144 @@
+test_that("sim_design", {  
+  between <- list(
+    B = c(B1 = "Level 1B", B2 = "Level 2B")
+  )
+  within <- list(
+    W = c(W1 = "Level 1W", W2 = "Level 2W")
+  )
+  
+  vardesc <- list(B = "Between-Subject Factor",
+                  W = "Within-Subject Factor")
+  
+  # wide ----
+  dat <- sim_design(within, between, vardesc = vardesc, plot = FALSE)
+  design <- get_design(dat)
+  expect_mapequal(design$vardesc, vardesc)
+  
+  obsID <- attr(dat$id, "label")
+  expect_equal(obsID, "id")
+  
+  obsB <- attr(dat$B, "label")
+  expect_equal(obsB, vardesc$B)
+  
+  obsW1 <- attr(dat$W1, "label")
+  expect_equal(obsW1, within$W[['W1']])
+  
+  obsW2 <- attr(dat$W2, "label")
+  expect_equal(obsW2, within$W[['W2']])
+
+  # long ----
+  dat <- sim_design(within, between, 
+                    vardesc = vardesc, 
+                    long = TRUE, plot = FALSE)
+  design <- get_design(dat)
+  expect_mapequal(design$vardesc, vardesc)
+  
+  obsID <- attr(dat$id, "label")
+  expect_equal(obsID, "id")
+  
+  obsY <- attr(dat$y, "label")
+  expect_equal(obsY, "value")
+  
+  obsB <- attr(dat$B, "label")
+  expect_equal(obsB, vardesc$B)
+  
+  obsW <- attr(dat$W, "label")
+  expect_equal(obsW, vardesc$W)
+})
+
+test_that("sim_design wide 2-within", {  
+  within <- list(
+    A = c(A1 = "Level 1A", A2 = "Level 2A"),
+    B = c(B1 = "Level 1B", B2 = "Level 2B")
+  )
+  
+  vardesc <- list(A = "First Factor",
+                  B = "Second Factor")
+  
+  dat <- sim_design(within, vardesc = vardesc, plot = FALSE)
+  
+  obsID <- attr(dat$id, "label")
+  expect_equal(obsID, "id")
+  
+  obs <- attr(dat$A1_B1, "label")
+  expect_equal(obs, "Level 1A:Level 1B")
+  
+  obs <- attr(dat$A1_B2, "label")
+  expect_equal(obs, "Level 1A:Level 2B")
+  
+  obs <- attr(dat$A2_B1, "label")
+  expect_equal(obs, "Level 2A:Level 1B")
+  
+  obs <- attr(dat$A2_B2, "label")
+  expect_equal(obs, "Level 2A:Level 2B")
+})
+
+test_that("sim_design no within", {  
+  dat <- sim_design(id = c(stim_id = "Stimulus ID"),
+                    dv = c(age = "Age"), 
+                    plot = FALSE)
+  
+  obsID <- attr(dat$stim_id, "label")
+  expect_equal(obsID, "Stimulus ID")
+  
+  obsY <- attr(dat$age, "label")
+  expect_equal(obsY, "Age")
+})
+
+test_that("sim_design rep", {  
+  # nested ----
+  dat <- sim_design(n = 20, rep = 2, plot = FALSE)
+  
+  obsID <- attr(dat$data[[1]]$id, "label")
+  expect_equal(obsID, "id")
+  
+  obsY <- attr(dat$data[[1]]$y, "label")
+  expect_equal(obsY, "value")
+  
+  expect_equal(attr(dat$data, "label"), "data")
+  expect_equal(attr(dat$rep, "label"), "replicate index")
+  
+  # unnested ----
+  dat <- sim_design(n = 20, rep = 2, nested = FALSE, plot = FALSE)
+  
+  expect_equal(attr(dat$id,  "label"), "id")
+  expect_equal(attr(dat$y,   "label"), "value")
+  expect_equal(attr(dat$rep, "label"), "replicate index")
+})
+
+test_that("long2wide/wide2long", {
+  between <- list(
+    B = c(B1 = "Level 1B", B2 = "Level 2B")
+  )
+  within <- list(
+    W = c(W1 = "Level 1W", W2 = "Level 2W")
+  )
+  vardesc <- list(B = "Between-Subject Factor",
+                  W = "Within-Subject Factor")
+  
+  # wide2long
+  dat_w <- sim_design(within, between, vardesc = vardesc, 
+                      long = FALSE, plot = FALSE)
+  l <- wide2long(dat_w)
+  
+  obsB <- attr(l$B, "label")
+  expect_equal(obsB, vardesc$B)
+  
+  obsW <- attr(l$W, "label")
+  expect_equal(obsW, vardesc$W)
+  
+  # long2wide
+  dat_l <- sim_design(within, between, vardesc = vardesc, 
+                      long = TRUE, plot = FALSE)
+  w <- long2wide(dat_l)
+  
+  obsB <- attr(w$B, "label")
+  expect_equal(obsB, vardesc$B)
+  
+  obsW1 <- attr(w$W1, "label")
+  expect_equal(obsW1, within$W[['W1']])
+  
+  obsW2 <- attr(w$W2, "label")
+  expect_equal(obsW2, within$W[['W2']])
+  
+})
diff --git a/tests/testthat/test-long2wide.R b/tests/testthat/test-long2wide.R
index 19b9a493..68850d20 100644
--- a/tests/testthat/test-long2wide.R
+++ b/tests/testthat/test-long2wide.R
@@ -132,4 +132,5 @@ test_that("from design", {
   expect_equal(names(dwide2), c("rater_id", "rater_sex", "black", "east_asian", "west_asian", "white"))
 })
 
+
 faux_options(plot = TRUE)
\ No newline at end of file
diff --git a/tests/testthat/test-norta.R b/tests/testthat/test-norta.R
index d3b39af9..e5db814b 100644
--- a/tests/testthat/test-norta.R
+++ b/tests/testthat/test-norta.R
@@ -52,7 +52,7 @@ test_that("fh_bounds", {
 })
 
 test_that("fh_bounds variation check", {
-  skip("long checking function")
+  skip("long simulation")
   
   x <- lapply(1:50, function(x) {
     fh_bounds("pois", "binom",  
diff --git a/tests/testthat/test-sim_design.R b/tests/testthat/test-sim_design.R
index bdeabd1b..1a79adfe 100644
--- a/tests/testthat/test-sim_design.R
+++ b/tests/testthat/test-sim_design.R
@@ -596,7 +596,7 @@ test_that("unnested reps", {
   n <- 10
   df <- sim_design(2, n = n, rep = rep, nested = FALSE, plot = FALSE)
   expect_equal(nrow(df), rep*n)
-  expect_equal(df$rep, rep(1:rep, each = n))
+  expect_equivalent(df$rep, rep(1:rep, each = n), )
 })
 
 # empirical ----
@@ -706,7 +706,7 @@ test_that("vardesc", {
   )
   
   vardesc <- list(B = "Between-Subject Factor",
-                       W = "Within-Subject Factor")
+                  W = "Within-Subject Factor")
   
   expect_silent(dat <- sim_design(within, between, vardesc = vardesc))
   design <- get_design(dat)