diff --git a/R/PPP.R b/R/PPP.R
index 6a2f11c..ec44e45 100644
--- a/R/PPP.R
+++ b/R/PPP.R
@@ -187,7 +187,8 @@ print.postPredP <- function(x, ...) {
 }
 
 matAsList <- function(matrix) {
-  lapply(apply(matrix, 1, list), function(x) x[[1]])
+  unname(split(matrix, f = seq_len(nrow(matrix))))
+  # lapply(apply(matrix, 1, list), function(x) x[[1]])
 }
 
 
diff --git a/R/plotFit.R b/R/plotFit.R
index c108d72..ad9e5ce 100644
--- a/R/plotFit.R
+++ b/R/plotFit.R
@@ -5,9 +5,13 @@
 #' predictive distribution).
 #'
 #' @inheritParams posteriorPredictive
-#' @param stat whether to plot mean frequencies (\code{"mean"}) or covariances
-#'   of individual frequencies (\code{"cov"})
-#' @param ... arguments passed to \code{\link{boxplot}}
+#' @param stat whether to plot mean frequencies (\code{"mean"}), covariances
+#'   of frequencies (\code{"cov"}), standard deviations (\code{"sd"}), or
+#'   correlations (\code{"cor"})
+#' @param main main title for plot
+#' @param ylab label for y-axis
+#' @param col color for boxplots of predicted values
+#' @param ... further arguments passed to \code{\link{boxplot}}
 #'
 #' @details If posterior predictive p-values were computed when fitting the
 #'   model (e.g., by adding the argument \code{traitMPT(...,ppp=1000)} ), the
@@ -28,10 +32,13 @@
 plotFit <- function(
     fittedModel,
     M = 1000,
-    stat = "mean",
+    stat = c("mean", "cov", "cor", "sd"),
+    main = NULL,
+    ylab = NULL,
+    col = "gray",
     ...
 ) {
-  stat <- match.arg(stat, c("mean", "cov"))
+  stat <- match.arg(stat)
 
   # get information about model:
   tree <- fittedModel$mptInfo$MPT$Tree
@@ -54,42 +61,98 @@ plotFit <- function(
 
   if (stat == "mean") {
     # Plot mean frequencies:
+    predicted <- t(sapply(freq.list, colMeans))[, free_cats, drop = FALSE]
+    observed  <- colMeans(dat[, free_cats, drop = FALSE])
+    if(is.null(ylab)) ylab <- "Mean frequency"
+    if(is.null(main)) main <- "Observed (red) and predicted (gray) mean frequencies"
+    x_labels <- free_cats
 
-    pred <- t(sapply(freq.list, colMeans))
-    boxplot(pred[, free_cats],
-      xaxt = "n", col = "gray",
-      main = "Observed (red) and predicted (boxplot) mean frequencies", las = 1, ...
-    )
-    axis(1, seq_along(free_cats), labels = free_cats)
 
-    xx <- by(seq_along(free_cats), tree[cats %in% free_cats], mean)
-    axis(1, xx, TreeNames, tick = F, line = NA, mgp = c(3, 2.5, 0))
-    points(1:length(free_cats), colMeans(dat)[free_cats],
-      col = "red", cex = 1.4, pch = 17
-    )
-    abline(v = cumsum(table(tree) - 1)[1:(length(TreeNames) - 1)] + .5, col = "gray")
   } else if (stat == "cov") {
-    # Plot covariance of frequencies:
 
     nams <- outer(free_cats, free_cats, paste, sep = "-")
-    sel_cov <- nams[upper.tri(nams, diag = TRUE)]
-    K <- length(sel_cov)
+    sel_cov <- nams[lower.tri(nams, diag = TRUE)]
+
 
     # observed/predicted
     c.obs <- cov(dat[, free_cats])
     c.pred <- sapply(freq.list, function(xx) {
       cc <- cov(xx[, free_cats])
-      cc[upper.tri(cc, diag = TRUE)]
+      cc[lower.tri(cc, diag = TRUE)]
     })
 
-    boxplot(t(c.pred),
-      col = "gray", ylab = "Covariance",
-      main = "Observed (red) and predicted (gray) covariances",
-      xaxt = "n", las = 1, ...
-    )
-    abline(h = 0, lty = 1, col = "gray")
-    axis(1, 1:K, labels = nams[upper.tri(nams, diag = TRUE)], las = 2)
-    points(1:K, c.obs[upper.tri(c.obs, diag = TRUE)], col = 2, pch = 17)
-    abline(v = cumsum(seq(nrow(c.obs), 2, -1)) + .5, col = "lightgray")
+    predicted <- t(c.pred)
+    observed  <- c.obs[lower.tri(c.obs, diag = TRUE)]
+    if(is.null(ylab)) ylab <- "Covariance"
+    if(is.null(main)) main <- "Observed (red) and predicted (gray) covariances"
+    x_labels <- sel_cov
+  } else if (stat == "cor") {
+    nams <- outer(free_cats, free_cats, paste, sep = "-")
+    sel_cov <- nams[lower.tri(nams, diag = FALSE)]
+
+    c.obs <- cor(dat[, free_cats])
+    c.pred <- sapply(freq.list, function(xx) {
+      cc <- cor(xx[, free_cats])
+      cc[lower.tri(cc, diag = FALSE)]
+    })
+    predicted <- t(c.pred)
+    observed  <- c.obs[lower.tri(c.obs, diag = FALSE)]
+    if(is.null(ylab)) ylab <- "Correlation"
+    if(is.null(main)) main <- "Observed (red) and predicted (gray) correlations"
+    x_labels <- sel_cov
+
+  } else if (stat == "sd") {
+    nams <- free_cats
+    s.obs <- sqrt(diag(cov(dat[, free_cats])))
+    s.pred <- sapply(freq.list, function(xx) {
+      sqrt(diag(cov(xx[, free_cats])))
+    })
+    predicted <- t(s.pred)
+    observed <- s.obs
+    if(is.null(ylab)) ylab <- "Standard deviation"
+    if(is.null(main)) main <- "Observed (red) and predicted (gray) standard deviations"
+    x_labels <- nams
+  }
+
+
+  out <- list()
+
+  out$boxplot <- boxplot(
+    x = predicted
+    , col = col
+    , ylab = ylab
+    , main = main
+    , xaxt = "n"
+    , ...
+  )
+  out$points <- points(x = seq_along(observed), y = observed, col = 2, pch = 17)
+  out$axis   <- axis(side = 1, seq_along(x_labels), labels = x_labels, las = 2)
+  abline(h = 0, lty = 1, col = "gray")
+
+
+
+
+  if(stat %in% c("mean", "sd")) {
+    # We add tree names if means or SDs are plotted
+    xx <- by(seq_along(free_cats), tree[cats %in% free_cats], mean)
+    axis(1, xx, TreeNames, tick = F, line = NA, mgp = c(3, 2.5, 0))
+    abline(v = cumsum(table(tree) - 1)[1:(length(TreeNames) - 1)] + .5, col = "gray")
+  } else {
+    # If co-variances or correlations are plotted, we add separators between blocks of columns for each tree.
+    p <- nrow(c.obs)
+
+    if (stat == "cov") {
+      block_sizes <- p:1
+    } else if (stat == "cor") {
+      block_sizes <- (p - 1):1
+    }
+
+    # positions after each completed column block
+    ends <- cumsum(block_sizes)
+
+    abline(v = ends[-length(ends)] + 0.5, col = "lightgray")
   }
+
+  # invisibly return
+  invisible(out)
 }
diff --git a/man/plotFit.Rd b/man/plotFit.Rd
index 649acbf..2727497 100644
--- a/man/plotFit.Rd
+++ b/man/plotFit.Rd
@@ -4,17 +4,32 @@
 \alias{plotFit}
 \title{Plot Posterior Predictive Mean Frequencies}
 \usage{
-plotFit(fittedModel, M = 1000, stat = "mean", ...)
+plotFit(
+  fittedModel,
+  M = 1000,
+  stat = c("mean", "cov", "cor", "sd"),
+  main = NULL,
+  ylab = NULL,
+  col = "gray",
+  ...
+)
 }
 \arguments{
 \item{fittedModel}{fitted latent-trait or beta MPT model (\code{\link{traitMPT}}, \code{\link{betaMPT}})}
 
 \item{M}{number of posterior predictive samples. As a maximum, the number of posterior samples in \code{fittedModel} is used.}
 
-\item{stat}{whether to plot mean frequencies (\code{"mean"}) or covariances
-of individual frequencies (\code{"cov"})}
+\item{stat}{whether to plot mean frequencies (\code{"mean"}), covariances
+of frequencies (\code{"cov"}), standard deviations (\code{"sd"}), or
+correlations (\code{"cor"})}
 
-\item{...}{arguments passed to \code{\link{boxplot}}}
+\item{main}{main title for plot}
+
+\item{ylab}{label for y-axis}
+
+\item{col}{color for boxplots of predicted values}
+
+\item{...}{further arguments passed to \code{\link{boxplot}}}
 }
 \description{
 Plots observed means/covariances of individual frequencies against the
diff --git a/man/plotPriorPost.Rd b/man/plotPriorPost.Rd
index cb8e823..91513d4 100644
--- a/man/plotPriorPost.Rd
+++ b/man/plotPriorPost.Rd
@@ -28,7 +28,7 @@ plotPriorPost(
 \item{nCPU}{number of CPUs used for parallel sampling. For large models and
 many participants, this may require a lot of memory.}
 
-\item{...}{arguments passed to \code{\link{boxplot}}}
+\item{...}{further arguments passed to \code{\link{boxplot}}}
 }
 \description{
 Allows to judge how much the data informed the parameter posterior