diff --git a/.gitignore b/.gitignore
index 5b6a065..0fc4bfa 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,4 @@
 .Rhistory
 .RData
 .Ruserdata
+
diff --git a/DESCRIPTION b/DESCRIPTION
index 9535c52..5ffea75 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -32,7 +32,8 @@ LazyData: true
 Imports: 
     quadprog,
     Rcpp,
-    pbapply
+    pbapply,
+    ggplot2
 LinkingTo: 
     Rcpp,
     RcppArmadillo,
diff --git a/NAMESPACE b/NAMESPACE
index 8c76b79..3bf12f1 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,5 +1,11 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(plot,kliep)
+S3method(plot,lhss)
+S3method(plot,naivedensityratio)
+S3method(plot,naivesubspacedensityratio)
+S3method(plot,spectral)
+S3method(plot,ulsif)
 S3method(predict,kliep)
 S3method(predict,lhss)
 S3method(predict,naivedensityratio)
@@ -29,9 +35,32 @@ export(kmm)
 export(lhss)
 export(naive)
 export(naivesubspace)
+export(plot_bivariate)
+export(plot_univariate)
 export(spectral)
 export(ulsif)
 importFrom(Rcpp,sourceCpp)
+importFrom(ggplot2,aes)
+importFrom(ggplot2,facet_grid)
+importFrom(ggplot2,facet_wrap)
+importFrom(ggplot2,geom_histogram)
+importFrom(ggplot2,geom_hline)
+importFrom(ggplot2,geom_point)
+importFrom(ggplot2,ggplot)
+importFrom(ggplot2,ggplotGrob)
+importFrom(ggplot2,labs)
+importFrom(ggplot2,position_dodge2)
+importFrom(ggplot2,scale_color_viridis_d)
+importFrom(ggplot2,scale_colour_gradient2)
+importFrom(ggplot2,scale_fill_viridis_d)
+importFrom(ggplot2,scale_shape_manual)
+importFrom(ggplot2,scale_x_continuous)
+importFrom(ggplot2,scale_y_continuous)
+importFrom(ggplot2,theme)
+importFrom(ggplot2,theme_bw)
+importFrom(grid,grid.draw)
+importFrom(grid,grid.newpage)
+importFrom(grid,nullGrob)
 importFrom(parallel,detectCores)
 importFrom(parallel,makeCluster)
 importFrom(parallel,stopCluster)
diff --git a/R/checks.R b/R/checks.R
index 0a504d0..c9f622f 100644
--- a/R/checks.R
+++ b/R/checks.R
@@ -439,3 +439,53 @@ check.newdata <- function(object, newdata) {
   }
   newdata
 }
+
+check.var.names <- function(vars, data){
+  if(!all(vars %in% names(data))) {
+  stop("Indicated variable (s) are not present in object. Check whether variable names are correct")
+  }
+}
+
+check.object.type <- function(object) {
+  models <- c(
+    "ulsif",
+    "kliep",
+    "lhss",
+    "spectral",
+    "naivedensityratio",
+    "naivesubspacedensityratio"
+  )
+
+  if (!attr(object, "class") %in% models) {
+    stop("Only densityratio objects can be plotted.")
+  }
+}
+
+check.logscale <- function(ext, logscale, tol){
+  if (logscale) {
+    # Convert values lower than tolerance to tol
+    negdr <- ext$dr < tol
+    ext$dr[negdr] <- tol
+    if (any(negdr)) {
+      warning(
+        paste(
+          "Negative estimated density ratios for", sum(negdr),
+          "observation(s) converted to",tol,
+          "before applying logarithmic transformation"
+        ),
+        call. = FALSE
+      )
+    }
+
+    # Apply log transformation
+    ext$dr <- log(ext$dr)
+
+    # Set y axis intercept to 0
+    ext$yintercept <- 0
+  } else {
+    # Set y axis intercept to 1
+    ext$yintercept <- 1
+  }
+  return(ext)
+}
+
diff --git a/R/plot.R b/R/plot.R
new file mode 100644
index 0000000..731c434
--- /dev/null
+++ b/R/plot.R
@@ -0,0 +1,511 @@
+
+#' A histogram of density ratio estimates
+#'
+#' Creates a histogram of the density ratio estimates. Useful to understand the
+#' distribution of estimated density ratios in each sample, or compare it among
+#' samples. It is the default plotting method for density ratio objects.
+#'
+#' @param x Density ratio object created with e.g., [kliep()], [ulsif()],
+#'   or [naive()]
+#' @param samples Character string indicating whether to plot the 'numerator',
+#'   'denominator', or 'both' samples. Default is 'both'.
+#' @param tol Numeric indicating the tolerance: values below this value will be set to the tolerance value, for legibility of the plots
+#' @param logscale Logical indicating whether to plot the density ratio
+#' estimates on a log scale. Default is TRUE.
+#' @param binwidth Numeric indicating the width of the bins, passed on to
+#'   `ggplot2`.
+#' @param bins Numeric indicating the number of bins. Overriden by binwidth, and
+#'   passed on to `ggplot2`.
+#' @param ... Additional arguments passed on to `predict()`.
+#'
+#' @return A histogram of density ratio estimates.
+#' @importFrom ggplot2 ggplot
+#' @importFrom ggplot2 aes
+#' @importFrom ggplot2 position_dodge2
+#' @importFrom ggplot2 geom_histogram
+#' @importFrom ggplot2 scale_fill_viridis_d
+#' @importFrom ggplot2 theme_bw
+#' @importFrom ggplot2 labs
+#'
+#'
+dr.histogram <- function(x,
+                         samples = "both",
+                         logscale = TRUE,
+                         binwidth = NULL,
+                         bins = NULL,
+                         tol = 10e-3,
+                         ...) {
+
+  # Check object type
+  check.object.type(x)
+
+
+  # Create data object and estimate density ratio
+  data <- rbind(
+    check.datatype(x$df_numerator),
+    check.datatype(x$df_denominator)
+  )
+  ext <- data.frame(dr = predict(x, newdata = data, ...),
+                    sample = c(rep("numerator", nrow(x$df_numerator)),
+                               rep("denominator", nrow(x$df_denominator))))
+
+  # If logscale = TRUE, transform density ratio estimates
+  ext <- check.logscale(ext, logscale, tol)
+
+  # Assign x-axis label
+  if (logscale) {
+    x_lab <- "Log (Density Ratio)"
+  } else {
+    x_lab <- "Density Ratio"
+  }
+
+  # Filter correct subset
+  samples <- match.arg(samples, c("both", "numerator", "denominator"))
+  if(samples != "both"){
+  data <- subset(data, ext$sample == samples)
+  ext  <- subset(ext, ext$sample == samples)
+  }
+
+  # Plot
+  plot <-
+    ggplot2::ggplot(data, ggplot2::aes(x = ext$dr)) +
+    ggplot2::geom_histogram(ggplot2::aes(fill = ext$sample),
+                   alpha = .75,
+                   color = "black",
+                   binwidth = binwidth,
+                   bins = bins,
+                   position = ggplot2::position_dodge2(preserve = "single",
+                                              padding = 0.2,
+                                              reverse = TRUE)) +
+
+    ggplot2::scale_fill_viridis_d(option = "cividis",
+                        breaks = c("numerator", "denominator"),
+                        labels = c("Numerator", "Denominator")) +
+    ggplot2::theme_bw()  +
+    ggplot2::labs(
+      x = x_lab,
+      y = "Count",
+      title = "Distribution of density ratio estimates",
+      fill = "Sample")
+
+  return(plot)
+}
+
+
+
+#' @rdname dr.histogram
+#' @return A histogram of density ratio estimates.
+#' @method plot ulsif
+#' @export
+#'
+plot.ulsif <- function(x, samples = "both", logscale = TRUE, binwidth = NULL,
+                       bins = NULL, tol = 10e-3, ...) {
+  dr.histogram(x, samples = samples, logscale = logscale, binwidth = binwidth,
+               bins = bins, tol = tol)
+}
+
+
+#' @rdname dr.histogram
+#' @return A histogram of density ratio estimates.
+#' @method plot kliep
+#' @export
+#'
+
+plot.kliep <- function(x, samples = "both", logscale = TRUE, binwidth = NULL,
+                       bins = NULL, tol = 10e-3, ...) {
+  dr.histogram(x, samples = samples, logscale = logscale, binwidth = binwidth,
+               bins = bins, tol = tol)
+}
+
+#' @rdname dr.histogram
+#' @return A histogram of density ratio estimates.
+#' @method plot spectral
+#' @export
+#'
+
+plot.spectral <- function(x, samples = "both", logscale = TRUE, binwidth = NULL,
+                       bins = NULL, tol = 10e-3, ...) {
+  dr.histogram(x, samples = samples, logscale = logscale, binwidth = binwidth,
+               bins = bins, tol = tol)
+}
+
+#' @rdname dr.histogram
+#' @return A histogram of density ratio estimates.
+#' @method plot lhss
+#' @export
+#'
+plot.lhss <- function(x, samples = "both", logscale = TRUE, binwidth = NULL,
+                       bins = NULL, tol = 10e-3, ...) {
+  dr.histogram(x, samples = samples, logscale = logscale, binwidth = binwidth,
+               bins = bins, tol = tol)
+}
+
+#' @rdname dr.histogram
+#' @return A histogram of density ratio estimates.
+#' @method plot naivedensityratio
+#' @export
+#'
+plot.naivedensityratio <- function(x, samples = "both", logscale = TRUE,
+                                   binwidth = NULL, bins = NULL, tol = 10e-3, ...) {
+  dr.histogram(x, samples = samples, logscale = logscale, binwidth = binwidth,
+               bins = bins, tol = tol)
+}
+
+#' @rdname dr.histogram
+#' @return A histogram of density ratio estimates.
+#' @method plot naivesubspacedensityratio
+#' @export
+#'
+plot.naivesubspacedensityratio <- function(x, samples = "both", logscale = TRUE,
+                                           binwidth = NULL, bins = NULL,
+                                           tol = 10e-3, ...) {
+  dr.histogram(x, samples = samples, logscale = logscale, binwidth = binwidth,
+               bins = bins, tol = tol)
+}
+
+#' Indivual univariate plot
+#'
+#' Scatterplot of individual values and density ratio estimates. Used internally in [create_univariate_plot()]
+#' @param data Data frame with the individual values and density ratio estimates
+#' @param var Name of the variable to be plotted on the x-axis
+#' @param y_lab Name of the y-axis label, typically ("Density Ratio" or "Log Density Ratio")
+#' @param ext Data frame with the density ratio estimates and sample indicator
+#' @param sample.facet Logical indicating whether to facet the plot by sample. Default is TRUE.
+#'
+#' @return A scatterplot of variable values and density ratio estimates.
+#'
+#' @importFrom ggplot2 ggplot
+#' @importFrom ggplot2 aes
+#' @importFrom ggplot2 geom_point
+#' @importFrom ggplot2 theme_bw
+#' @importFrom ggplot2 labs
+#' @importFrom ggplot2 scale_color_viridis_d
+#' @importFrom ggplot2 scale_y_continuous
+#' @importFrom ggplot2 facet_wrap
+#' @importFrom ggplot2 geom_hline
+#'
+#'
+create_univariate_plot <- function(data, ext, var, y_lab, sample.facet = TRUE){
+
+  y_max <- max(2, ext$dr)
+  y_min <- min(-2, ext$dr)
+
+  plot <-
+    ggplot2::ggplot(data, ggplot2::aes(x = .data[[var]], y = ext$dr)) +
+    ggplot2::geom_point(ggplot2::aes(col = ext$sample),
+               alpha = .6) +
+    ggplot2::theme_bw() +
+    ggplot2::labs(title = "Scatter plot of individual values and density ratio",
+         color = "Sample",
+         y = y_lab) +
+    ggplot2::scale_color_viridis_d(option = "cividis",
+                         breaks = c("numerator", "denominator"),
+                         labels = c("Numerator", "Denominator")) +
+    ggplot2::scale_y_continuous(limits = c(y_min, y_max))
+
+  if(sample.facet){
+    plot <- plot + ggplot2::facet_wrap(~ext$sample)
+  }
+
+  plot <- plot +
+    ggplot2::geom_hline(yintercept = ext$yintercept, linetype = "dashed")
+
+  return(plot)
+}
+
+#' Scatter plot of density ratios and individual variables
+#'
+#' A scatter plot showing the relationship between estimated density ratios and individual variables.
+#'
+#' @inheritParams dr.histogram
+#' @param vars Character vector of variable names to be plotted.
+#' @param grid Logical indicating whether output should be a list of individual plots ("individual"), or one facetted plot with all variables ("assembled"). Defaults to "individual".
+#' @param sample.facet Logical indicating whether to facet the plot by sample, i.e, showing plots separate for each sample, and side to side. Defaults to FALSE.
+#' @param nrow Integer indicating the number of rows in the assembled plot. If NULL, the number of rows is automatically calculated.
+#' @param ... Additional arguments passed to the predict() function.
+#'
+#' @return Scatter plot of density ratios and individual variables.
+#' @export
+#'
+#' @importFrom ggplot2 ggplot
+#' @importFrom ggplot2 aes
+#' @importFrom ggplot2 geom_point
+#' @importFrom ggplot2 geom_hline
+#' @importFrom ggplot2 theme_bw
+#' @importFrom ggplot2 labs
+#' @importFrom ggplot2 scale_color_viridis_d
+#' @importFrom ggplot2 scale_y_continuous
+#' @importFrom ggplot2 facet_wrap
+#' @importFrom ggplot2 facet_grid
+#'
+
+plot_univariate <- function(x, vars = NULL, samples = "both", logscale = TRUE,
+                            grid = FALSE, sample.facet = FALSE,
+                            nrow = NULL, tol = 10e-3, ...) {
+
+  nu <- check.datatype(x$df_numerator)
+  de <- check.datatype(x$df_denominator)
+
+  if(is.null(vars)){
+    vars <- names(nu)
+  }
+  # Check object type
+  check.object.type(x)
+
+  # Create data object, and external object with density ratio and sample indicators
+  data <- rbind(nu, de)
+  ext <- data.frame(dr = predict(x, newdata = data, ...),
+                    sample = c(rep("numerator", nrow(x$df_numerator)),
+                               rep("denominator", nrow(x$df_denominator))))
+  # Check variable names
+  check.var.names(vars, data)
+
+  # Filter appropriate subset
+  samples <- match.arg(samples, c("both", "numerator", "denominator"))
+  if(samples != "both"){
+    data <- subset(data, ext$sample == samples)
+    ext  <- subset(ext, ext$sample == samples)
+  }
+
+
+  # If logscale is TRUE, transform density ratio estimates
+  ext <- check.logscale(ext, logscale, tol)
+
+  # Set y axis label
+  if(logscale) {
+    y_lab <- "Log (Density Ratio)"
+  }  else {
+    y_lab <- "Density Ratio"
+  }
+
+  # Plot either individual plots in a list, or a grid of individual plots
+  if(!grid){
+
+  plot <- lapply(vars, function(var) create_univariate_plot(data, ext, var, y_lab, sample.facet))
+
+  } else {
+
+    values <- data[, vars] |> unlist(use.names = FALSE)
+    variable <- rep(vars, each = length(values)/length(vars))
+    dr <- rep(ext$dr, length(vars))
+    sample <- rep(ext$sample, length(vars))
+    ext <- data.frame(values = values, variable = variable,
+                      dr = dr, sample = sample, yintercept = rep(ext$yintercept, length(vars)))
+
+    # Maximum scale for y
+    y_max <- max(1,ext$dr)
+    y_min <- min(-1, ext$dr)
+
+    plot <- ggplot2::ggplot(ext) +
+      ggplot2::geom_point(ggplot2::aes(x = values, y = dr, col = sample),
+                 alpha = .6) +
+      ggplot2::theme_bw() +
+      ggplot2::labs(title = "Scatter plot of individual values and density ratio",
+           color = "Sample",
+           y = "Density ratio") +
+      ggplot2::scale_color_viridis_d(option = "cividis",
+                            breaks = c("numerator", "denominator"),
+                            labels = c("Numerator", "Denominator")) +
+      ggplot2::scale_y_continuous(limits = c(y_min, y_max))
+
+    if(sample.facet){
+
+      plot <- plot +
+        ggplot2::facet_grid(rows = ggplot2::vars(sample),
+                            cols = ggplot2::vars(variable),
+                   scales = "free_x")
+
+      } else {
+        plot <- plot +
+          ggplot2::facet_wrap(ggplot2::vars(variable), scales = "free_x", nrow = nrow)
+      }
+
+  plot <- plot +
+    ggplot2::geom_hline(yintercept = ext$yintercept, linetype = "dashed")
+  }
+
+  return(plot)
+}
+
+
+#' Bivariate plot
+#'
+#' @inheritParams create_univariate_plot
+#' @param show.sample Logical indicating whether to give different shapes to observations, depending on the sample they come from (numerator or denominator). Defaults to FALSE.
+#' @param vars Character vector of variable names to be plotted.
+#' @param logscale Logical indicating whether the density ratio should be plotted in log scale. Defaults to TRUE.
+#'
+#' @return Bivariate plot
+#'
+#'
+#' @importFrom ggplot2 ggplot
+#' @importFrom ggplot2 aes
+#' @importFrom ggplot2 geom_point
+#' @importFrom ggplot2 scale_colour_gradient2
+#' @importFrom ggplot2 scale_shape_manual
+#' @importFrom ggplot2 theme_bw
+#' @importFrom ggplot2 labs
+#'
+create_bivariate_plot  <- function(data, ext, vars, logscale, show.sample){
+
+
+  dr_max <- ifelse(logscale, max(2, ext$dr), max(exp(2), ext$dr))
+  dr_min <- ifelse(logscale, min(-2, ext$dr), min(exp(-2), ext$dr))
+
+  plot <-
+    ggplot2::ggplot(data, mapping = ggplot2::aes(x = .data[[vars[1]]], y = .data[[vars[2]]])) +
+    ggplot2::geom_point(ggplot2::aes(colour = ext$dr, shape = show.sample),
+               alpha = 1,
+               size = 2.0) +
+    ggplot2::scale_colour_gradient2(low = "#00204DFF",
+                           high = "#7D0000",
+                           mid = "lightyellow",
+                           midpoint = 0,
+                           limits = c(dr_min, dr_max)) +
+    ggplot2::theme_bw() +
+    ggplot2::labs(title = "Scatter plot, with density ratio mapped to colour",
+         colour = "Log (Density ratio)") +
+    ggplot2::scale_shape_manual(values = c(21, 24))
+
+  return(plot)
+}
+
+#' Densityratio in bidimensional plot
+#'
+#' Plots a scatterplot of two variables, with densityratio mapped to the colour
+#' scale.
+#'
+#' @inheritParams plot_univariate
+#' @param logscale Logical indicating whether to plot the density ratio
+#'   estimates on a log scale. Default is TRUE.
+#' @param show.sample Logical indicating whether to give different shapes to
+#' observations, depending on the sample they come from (numerator or
+#' denominator). Defaults to FALSE.
+#' @param vars Character vector of variable names for which all pairwise
+#' bivariate plots are created
+#'
+#' @return Bivariate scatter plots of all combinations of variables in vars.
+#' @export
+#'
+#' @importFrom ggplot2 ggplot
+#' @importFrom ggplot2 aes
+#' @importFrom ggplot2 geom_point
+#' @importFrom ggplot2 facet_grid
+#' @importFrom ggplot2 scale_colour_gradient2
+#' @importFrom ggplot2 scale_y_continuous
+#' @importFrom ggplot2 scale_x_continuous
+#' @importFrom ggplot2 theme_bw
+#' @importFrom ggplot2 theme
+#' @importFrom ggplot2 labs
+#' @importFrom ggplot2 scale_shape_manual
+#' @importFrom ggplot2 ggplotGrob
+#' @importFrom grid grid.draw
+#' @importFrom grid grid.newpage
+#' @importFrom grid nullGrob
+#'
+plot_bivariate <- function(x, vars, samples = "both", grid = FALSE,
+                           logscale = TRUE, show.sample = NULL, tol = 10e-3, ...) {
+
+  if (length(vars) <= 1) {
+    stop("At least two variables are required for a bivariate plot.")
+  }
+
+  # Check object type
+  check.object.type(x)
+
+  # Create data object and estimate density ratio
+
+  data <- rbind(
+    check.datatype(x$df_numerator),
+    check.datatype(x$df_denominator)
+  )
+
+  # Check variable names
+  check.var.names(vars, data)
+  var_combinations <- expand.grid(vars, vars)
+
+  ext <- data.frame(dr = predict(object, newdata = data, ...),
+                    sample = c(rep("numerator", nrow(x$df_numerator)),
+                               rep("denominator", nrow(x$df_denominator))))
+
+  # Check if logscale is TRUE, then change ext
+  ext <- check.logscale(ext, logscale, tol)
+
+  # Assign correct labels depending on logscale
+  if(logscale){
+    colour_label <- "Log (Density Ratio)"
+  } else {
+    colour_label <- "Density Ratio"
+  }
+
+
+  # Select data
+  samples <- match.arg(samples, c("both", "numerator", "denominator"))
+  if(samples != "both"){
+    data <- subset(data, ext$sample == samples)
+    ext  <- subset(ext, ext$sample == samples)
+  }
+
+  if(!grid){
+
+  plot <- lapply(var_combinations, function(vars) create_bivariate_plot(data, ext, vars, logscale, show.sample))
+  return(plot)
+
+  } else {
+
+    ext <- data.frame(data, dr = ext$dr, sample = ext$sample)
+    datlist <- lapply(var_combinations, \(x) data.frame(values.x = ext[,x[1]],
+                                              values.y = ext[,x[2]],
+                                              xvar = rep(x[1], nrow(ext)),
+                                              yvar = rep(x[2], nrow(ext)),
+                                              sample = ext[, "sample"],
+                                              dr = ext[, "dr"]))
+
+    plot_data <- do.call(datlist, what = rbind)
+
+    dr_max <- max(1, ext$dr)
+    dr_min <- min(-1, ext$dr)
+
+    plot <-
+      ggplot2::ggplot(plot_data, mapping = ggplot2::aes(x = values.x, y = values.y,
+                                      shape = show.sample)) +
+      ggplot2::geom_point(ggplot2::aes(colour = dr),
+                 size = 2.0) +
+      ggplot2::facet_grid(rows = ggplot2::vars(yvar), cols = ggplot2::vars(xvar), scales = "free",
+                 switch = "both") +
+      ggplot2::scale_colour_gradient2(low = "#00204DFF",
+                             high = "#7D0000",
+                             mid = "lightyellow",
+                             midpoint = 0,
+                             limits = c(dr_min, dr_max)) +
+      ggplot2::scale_y_continuous(position = "left") +
+      ggplot2::scale_x_continuous(position = "bottom") +
+      ggplot2::theme_bw() +
+      ggplot2::theme(strip.placement = "outside") +
+      ggplot2::labs(title = "Scatter plots, with density ratio mapped to colour",
+          x = NULL,
+          y = NULL,
+          colour = colour_label,
+          shape = show.sample) +
+      ggplot2::scale_shape_manual(values = c(21, 24))
+
+  # Erase upper diagonal
+  ## Create plot into a grob
+  grob <- ggplot2::ggplotGrob(plot)
+
+  ## Create name of empty panels in the upper diagonal
+  empty_panels <- expand.grid(seq(1:length(vars)), seq(1:length(vars))) |>
+    subset(Var2 > Var1)
+  empty_panels$panel <- paste0("panel-", empty_panels$Var1, "-", empty_panels$Var2)
+
+  # Delete panels in upper diagonal, based in their index
+  idx <- which(grob$layout$name %in% empty_panels$panel)
+  for (i in idx) grob$grobs[[i]] <- grid::nullGrob()
+
+  out <- grob
+
+  grid::grid.newpage()
+  grid::grid.draw(out)
+
+  }
+}
diff --git a/man/create_bivariate_plot.Rd b/man/create_bivariate_plot.Rd
new file mode 100644
index 0000000..9e1bf6b
--- /dev/null
+++ b/man/create_bivariate_plot.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot.R
+\name{create_bivariate_plot}
+\alias{create_bivariate_plot}
+\title{Bivariate plot}
+\usage{
+create_bivariate_plot(data, ext, vars, logscale, show.sample)
+}
+\arguments{
+\item{data}{Data frame with the individual values and density ratio estimates}
+
+\item{ext}{Data frame with the density ratio estimates and sample indicator}
+
+\item{vars}{Character vector of variable names to be plotted.}
+
+\item{logscale}{Logical indicating whether the density ratio should be plotted in log scale. Defaults to TRUE.}
+
+\item{show.sample}{Logical indicating whether to give different shapes to observations, depending on the sample they come from (numerator or denominator). Defaults to FALSE.}
+}
+\value{
+Bivariate plot
+}
+\description{
+Bivariate plot
+}
diff --git a/man/create_univariate_plot.Rd b/man/create_univariate_plot.Rd
new file mode 100644
index 0000000..0248bc9
--- /dev/null
+++ b/man/create_univariate_plot.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot.R
+\name{create_univariate_plot}
+\alias{create_univariate_plot}
+\title{Indivual univariate plot}
+\usage{
+create_univariate_plot(data, ext, var, y_lab, sample.facet = TRUE)
+}
+\arguments{
+\item{data}{Data frame with the individual values and density ratio estimates}
+
+\item{ext}{Data frame with the density ratio estimates and sample indicator}
+
+\item{var}{Name of the variable to be plotted on the x-axis}
+
+\item{y_lab}{Name of the y-axis label, typically ("Density Ratio" or "Log Density Ratio")}
+
+\item{sample.facet}{Logical indicating whether to facet the plot by sample. Default is TRUE.}
+}
+\value{
+A scatterplot of variable values and density ratio estimates.
+}
+\description{
+Scatterplot of individual values and density ratio estimates. Used internally in \code{\link[=create_univariate_plot]{create_univariate_plot()}}
+}
diff --git a/man/dr.histogram.Rd b/man/dr.histogram.Rd
new file mode 100644
index 0000000..a7f5645
--- /dev/null
+++ b/man/dr.histogram.Rd
@@ -0,0 +1,122 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot.R
+\name{dr.histogram}
+\alias{dr.histogram}
+\alias{plot.ulsif}
+\alias{plot.kliep}
+\alias{plot.spectral}
+\alias{plot.lhss}
+\alias{plot.naivedensityratio}
+\alias{plot.naivesubspacedensityratio}
+\title{A histogram of density ratio estimates}
+\usage{
+dr.histogram(
+  x,
+  samples = "both",
+  logscale = TRUE,
+  binwidth = NULL,
+  bins = NULL,
+  tol = 0.01,
+  ...
+)
+
+\method{plot}{ulsif}(
+  x,
+  samples = "both",
+  logscale = TRUE,
+  binwidth = NULL,
+  bins = NULL,
+  tol = 0.01,
+  ...
+)
+
+\method{plot}{kliep}(
+  x,
+  samples = "both",
+  logscale = TRUE,
+  binwidth = NULL,
+  bins = NULL,
+  tol = 0.01,
+  ...
+)
+
+\method{plot}{spectral}(
+  x,
+  samples = "both",
+  logscale = TRUE,
+  binwidth = NULL,
+  bins = NULL,
+  tol = 0.01,
+  ...
+)
+
+\method{plot}{lhss}(
+  x,
+  samples = "both",
+  logscale = TRUE,
+  binwidth = NULL,
+  bins = NULL,
+  tol = 0.01,
+  ...
+)
+
+\method{plot}{naivedensityratio}(
+  x,
+  samples = "both",
+  logscale = TRUE,
+  binwidth = NULL,
+  bins = NULL,
+  tol = 0.01,
+  ...
+)
+
+\method{plot}{naivesubspacedensityratio}(
+  x,
+  samples = "both",
+  logscale = TRUE,
+  binwidth = NULL,
+  bins = NULL,
+  tol = 0.01,
+  ...
+)
+}
+\arguments{
+\item{x}{Density ratio object created with e.g., \code{\link[=kliep]{kliep()}}, \code{\link[=ulsif]{ulsif()}},
+or \code{\link[=naive]{naive()}}}
+
+\item{samples}{Character string indicating whether to plot the 'numerator',
+'denominator', or 'both' samples. Default is 'both'.}
+
+\item{logscale}{Logical indicating whether to plot the density ratio
+estimates on a log scale. Default is TRUE.}
+
+\item{binwidth}{Numeric indicating the width of the bins, passed on to
+\code{ggplot2}.}
+
+\item{bins}{Numeric indicating the number of bins. Overriden by binwidth, and
+passed on to \code{ggplot2}.}
+
+\item{tol}{Numeric indicating the tolerance: values below this value will be set to the tolerance value, for legibility of the plots}
+
+\item{...}{Additional arguments passed on to \code{predict()}.}
+}
+\value{
+A histogram of density ratio estimates.
+
+A histogram of density ratio estimates.
+
+A histogram of density ratio estimates.
+
+A histogram of density ratio estimates.
+
+A histogram of density ratio estimates.
+
+A histogram of density ratio estimates.
+
+A histogram of density ratio estimates.
+}
+\description{
+Creates a histogram of the density ratio estimates. Useful to understand the
+distribution of estimated density ratios in each sample, or compare it among
+samples. It is the default plotting method for density ratio objects.
+}
diff --git a/man/plot_bivariate.Rd b/man/plot_bivariate.Rd
new file mode 100644
index 0000000..69d500d
--- /dev/null
+++ b/man/plot_bivariate.Rd
@@ -0,0 +1,47 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot.R
+\name{plot_bivariate}
+\alias{plot_bivariate}
+\title{Densityratio in bidimensional plot}
+\usage{
+plot_bivariate(
+  x,
+  vars,
+  samples = "both",
+  grid = FALSE,
+  logscale = TRUE,
+  show.sample = NULL,
+  tol = 0.01,
+  ...
+)
+}
+\arguments{
+\item{x}{Density ratio object created with e.g., \code{\link[=kliep]{kliep()}}, \code{\link[=ulsif]{ulsif()}},
+or \code{\link[=naive]{naive()}}}
+
+\item{vars}{Character vector of variable names for which all pairwise
+bivariate plots are created}
+
+\item{samples}{Character string indicating whether to plot the 'numerator',
+'denominator', or 'both' samples. Default is 'both'.}
+
+\item{grid}{Logical indicating whether output should be a list of individual plots ("individual"), or one facetted plot with all variables ("assembled"). Defaults to "individual".}
+
+\item{logscale}{Logical indicating whether to plot the density ratio
+estimates on a log scale. Default is TRUE.}
+
+\item{show.sample}{Logical indicating whether to give different shapes to
+observations, depending on the sample they come from (numerator or
+denominator). Defaults to FALSE.}
+
+\item{tol}{Numeric indicating the tolerance: values below this value will be set to the tolerance value, for legibility of the plots}
+
+\item{...}{Additional arguments passed to the predict() function.}
+}
+\value{
+Bivariate scatter plots of all combinations of variables in vars.
+}
+\description{
+Plots a scatterplot of two variables, with densityratio mapped to the colour
+scale.
+}
diff --git a/man/plot_univariate.Rd b/man/plot_univariate.Rd
new file mode 100644
index 0000000..3519476
--- /dev/null
+++ b/man/plot_univariate.Rd
@@ -0,0 +1,46 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot.R
+\name{plot_univariate}
+\alias{plot_univariate}
+\title{Scatter plot of density ratios and individual variables}
+\usage{
+plot_univariate(
+  x,
+  vars = NULL,
+  samples = "both",
+  logscale = TRUE,
+  grid = FALSE,
+  sample.facet = FALSE,
+  nrow = NULL,
+  tol = 0.01,
+  ...
+)
+}
+\arguments{
+\item{x}{Density ratio object created with e.g., \code{\link[=kliep]{kliep()}}, \code{\link[=ulsif]{ulsif()}},
+or \code{\link[=naive]{naive()}}}
+
+\item{vars}{Character vector of variable names to be plotted.}
+
+\item{samples}{Character string indicating whether to plot the 'numerator',
+'denominator', or 'both' samples. Default is 'both'.}
+
+\item{logscale}{Logical indicating whether to plot the density ratio
+estimates on a log scale. Default is TRUE.}
+
+\item{grid}{Logical indicating whether output should be a list of individual plots ("individual"), or one facetted plot with all variables ("assembled"). Defaults to "individual".}
+
+\item{sample.facet}{Logical indicating whether to facet the plot by sample, i.e, showing plots separate for each sample, and side to side. Defaults to FALSE.}
+
+\item{nrow}{Integer indicating the number of rows in the assembled plot. If NULL, the number of rows is automatically calculated.}
+
+\item{tol}{Numeric indicating the tolerance: values below this value will be set to the tolerance value, for legibility of the plots}
+
+\item{...}{Additional arguments passed to the predict() function.}
+}
+\value{
+Scatter plot of density ratios and individual variables.
+}
+\description{
+A scatter plot showing the relationship between estimated density ratios and individual variables.
+}