diff --git a/NAMESPACE b/NAMESPACE
index 4b9e32c5..b87af9ea 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -17,6 +17,7 @@ export(plot_fishing_mortality)
 export(plot_indices)
 export(plot_landings)
 export(plot_natural_mortality)
+export(plot_obsvpred)
 export(plot_recruitment)
 export(plot_recruitment_deviations)
 export(plot_spawning_biomass)
diff --git a/R/plot_indices.R b/R/plot_indices.R
index e0ff78e7..102826b7 100644
--- a/R/plot_indices.R
+++ b/R/plot_indices.R
@@ -23,7 +23,7 @@ plot_indices <- function(
   unit_label = "",
   group = NULL,
   # facet always assigned to fleet since that is how indices are calc'd -- unless replaced with NULL
-  facet = "fleet",
+  facet = NULL,
   interactive = TRUE,
   module = NULL,
   focus = NULL,
@@ -63,55 +63,51 @@ plot_indices <- function(
     prepared_data <- prepared_data |>
       dplyr::filter(fleet %in% focus)
   }
+  
+  processed_data <- process_data(
+    dat = prepared_data,
+    group = group,
+    facet = facet
+  )
+  prepared_data <- processed_data[[1]]
+  group <- processed_data[[2]]
+  facet <- processed_data[[3]]
 
-  # identify if there is >1 label and create plot
-  if (length(unique(prepared_data$label)) > 1) {
-    # move label to grouping and set grouping into facet
-    if (!is.null(group)) {
-      facet <- c(facet, group)
-    }
-    # transform prep data so group_var = label
-    prepared_data <- prepared_data |>
-      dplyr::mutate(
-        group_var = label
-      )
-    # plot time series with multiple labels
-    plt <- plot_timeseries(
-      dat = prepared_data,
-      ylab = u_units,
-      group = "label",
-      facet = facet,
-      # linewidth = 1,
-      ...
-    ) +
-      # commenting out but might need this later -- not sure if this will always be true
-      ggplot2::labs(
-        linetype = "",
-        fill = ""
-      ) +
-      theme_noaa() +
-      ggplot2::scale_x_continuous(
-        breaks = ggplot2::waiver(),
-        # labels = scales::label_number(accuracy = 1),
-        guide = ggplot2::guide_axis(
-          minor.ticks = TRUE
-        )
+  # move label to grouping and set grouping into facet
+  if (!is.null(group)) {
+    facet <- c(facet, group)
+  }
+  # transform prep data so group_var = label
+  # prepared_data <- prepared_data |>
+  #   dplyr::mutate(
+  #     group_var = label
+  #   )
+  # plot time series with multiple labels
+  plt <- plot_obsvpred(
+    dat = prepared_data,
+    x = "year",
+    y = "estimate",
+    observed_label = "indices_observed",
+    predicted_label = "indices_predicted",
+    geom = "line",
+    xlab = "Year",
+    ylab = "Estimated Index",
+    group = group,
+    facet = facet
+  ) +
+    theme_noaa() +
+    ggplot2::scale_x_continuous(
+      breaks = ggplot2::waiver(),
+      # labels = scales::label_number(accuracy = 1),
+      guide = ggplot2::guide_axis(
+        minor.ticks = TRUE
       )
-    # Overwrite facets from base plot_timeseries bc scales need to be free
-    facet <- paste("~", paste(facet, collapse = " + "))
-    facet_formula <- stats::reformulate(facet)
-    plt <- plt + ggplot2::facet_wrap(facet_formula, scales = "free")
-  } else {
-    # plot time series
-    plt <- plot_error(
-      dat = prepared_data,
-      ylab = u_units,
-      group = group,
-      facet = facet,
-      ...
     )
-  }
-
+  # Overwrite facets from base plot_timeseries bc scales need to be free
+  facet <- paste("~", paste(facet, collapse = " + "))
+  facet_formula <- stats::reformulate(facet)
+  plt <- plt + ggplot2::facet_wrap(facet_formula, scales = "free")
+  
   ### Make RDA ----
   if (make_rda) {
     create_rda(
diff --git a/R/process_data.R b/R/process_data.R
index c58a351a..f30e3fb9 100644
--- a/R/process_data.R
+++ b/R/process_data.R
@@ -317,7 +317,16 @@ process_data <- function(
       group <- NULL
     }
   }
-
+  
+  # Ensure that index_variables -- group or facets are non-numeric to be plotted accurately
+  data <- data |>
+    dplyr::mutate(
+      across(
+        tidyselect::any_of(c(group, facet)),
+        as.character
+      )
+    )
+  
   # Export list of objects
   list(
     # variable,
diff --git a/R/utils_plot.R b/R/utils_plot.R
index cd630351..cc7b6a7d 100644
--- a/R/utils_plot.R
+++ b/R/utils_plot.R
@@ -973,3 +973,106 @@ check_grouping <- function(dat) {
   }
   dat_index
 }
+
+#------------------------------------------------------------------------------
+
+#' Plot observed vs. predicted data
+#'
+#' @inheritParams plot_timeseries 
+#' @param observed_label a string of the label used to filter the observed data. Default is "observed". 
+#' @param predicted_label a string of the label used to filter the predicted data. Default is "predicted".
+#'
+#' @returns Create a plot of observed vs. predicted data for a stock assessment report. 
+#' @export
+#'
+plot_obsvpred <- function(
+    dat,
+    x = "year",
+    y = "estimate",
+    observed_label = "observed",
+    predicted_label = "predicted",
+    geom = "line",
+    xlab = "Year",
+    ylab = NULL,
+    group = NULL,
+    facet = NULL,
+    ...
+) {
+  # Start plot
+  plot <- ggplot2::ggplot()
+  # make into new geom?
+  # more defaults and fxnality for ggplot
+  
+  # Add geom
+  plot <- plot +
+    ggplot2::geom_point(
+      data = dat |> dplyr::filter(grepl(observed_label, label)),
+      ggplot2::aes(
+        .data[[x]],
+        .data[[y]],
+        color = model
+      ),
+      shape = 16
+      # ...
+    ) +
+  ggplot2::geom_line(
+    data = dat |> dplyr::filter(grepl(predicted_label, label)),
+    ggplot2::aes(
+      x = .data[[x]],
+      y = .data[[y]],
+      color = model
+    ),
+    linetype = "solid"
+  )
+  
+  # Add labels to axis and legend
+  if (length(unique(dat$model)) > 1 & !is.null(group)) {
+    labs <- plot + ggplot2::labs(
+      x = xlab,
+      y = ylab
+      # color = "Model",
+      # linetype = cap_first_letter(group),
+      # fill = cap_first_letter(group),
+      # shape = cap_first_letter(group)
+    ) +
+      ggplot2::theme(legend.title = ggplot2::element_blank())
+  } else {
+    labs <- plot + ggplot2::labs(
+      x = xlab,
+      y = ylab,
+      color = "Model"
+    )
+  }
+  
+  # Remove linetype or point when there is no grouping
+  if (is.null(group) & length(unique(dat$model)) == 1) {
+    labs <- labs + ggplot2::guides(linetype = "none", shape = "none")
+  }
+  if (length(unique(dat$model)) == 1) {
+    labs <- labs + ggplot2::guides(color = "none")
+  }
+  
+  # Calc axis breaks
+  x_n_breaks <- axis_breaks(dat[[x]])
+  breaks <- ggplot2::scale_x_continuous(
+    breaks = x_n_breaks,
+    guide = ggplot2::guide_axis(
+      minor.ticks = TRUE
+    )
+  )
+  
+  # Put together final plot
+  final <- labs + breaks + ggplot2::expand_limits(y = 0) +
+    ggplot2::scale_y_continuous(
+      labels = scales::label_comma()
+    )
+  
+  # Check if facet(s) are desired
+  if (!is.null(facet) & length(facet) > 0) {
+    facet <- paste("~", paste(facet, collapse = " + "))
+    facet_formula <- stats::reformulate(facet)
+    
+    final <- final + ggplot2::facet_wrap(facet_formula)
+  }
+  final
+}
diff --git a/man/plot_indices.Rd b/man/plot_indices.Rd
index 1ef91088..8789064e 100644
--- a/man/plot_indices.Rd
+++ b/man/plot_indices.Rd
@@ -8,7 +8,7 @@ plot_indices(
   dat,
   unit_label = "",
   group = NULL,
-  facet = "fleet",
+  facet = NULL,
   interactive = TRUE,
   module = NULL,
   focus = NULL,
diff --git a/man/plot_obsvpred.Rd b/man/plot_obsvpred.Rd
new file mode 100644
index 00000000..9564ccca
--- /dev/null
+++ b/man/plot_obsvpred.Rd
@@ -0,0 +1,56 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils_plot.R
+\name{plot_obsvpred}
+\alias{plot_obsvpred}
+\title{Plot observed vs. predicted data}
+\usage{
+plot_obsvpred(
+  dat,
+  x = "year",
+  y = "estimate",
+  observed_label = "observed",
+  predicted_label = "predicted",
+  geom = "line",
+  xlab = "Year",
+  ylab = NULL,
+  group = NULL,
+  facet = NULL,
+  ...
+)
+}
+\arguments{
+\item{dat}{filtered data frame from standard output file(s) preformatted for
+the target label from \link[stockplotr]{filter_data}}
+
+\item{x}{a string of the column name of data used to plot on the x-axis (default
+is "year")}
+
+\item{y}{a string of the column name of data used to plot on the y-axis (default
+is "estimate")}
+
+\item{observed_label}{a string of the label used to filter the observed data. Default is "observed".}
+
+\item{predicted_label}{a string of the label used to filter the predicted data. Default is "predicted".}
+
+\item{geom}{type of geom to use for plotting found in ggplot2 (e.g. "point",
+"line", etc.). Default is "line". Other options are "point" and "area".}
+
+\item{xlab}{a string of the x-axis label (default is "Year")}
+
+\item{ylab}{a string of the y-axis label. If NULL, it will be set to the name
+of `y`.}
+
+\item{group}{a string of a single column that groups the data (e.g. "fleet",
+"sex", "area", etc.). Currently can only have one level of grouping.}
+
+\item{facet}{a string or vector of strings of a column that facets the data
+(e.g. "year", "area", etc.)}
+
+\item{...}{inherited arguments from internal functions from ggplot2::geom_xx}
+}
+\value{
+Create a plot of observed vs. predicted data for a stock assessment report.
+}
+\description{
+Plot observed vs. predicted data
+}