diff --git a/R/classicalmetaanalysis.R b/R/classicalmetaanalysis.R
index fd926c19..f86daaa7 100644
--- a/R/classicalmetaanalysis.R
+++ b/R/classicalmetaanalysis.R
@@ -54,7 +54,31 @@ ClassicalMetaAnalysis <- function(jaspResults, dataset = NULL, options, ...) {
   "permutationTest", "permutationTestIteration", "permutationTestType", "setSeed", "seed",
   # multilevel/multivariate specific
   "randomEffects", "randomEffectsSpecification",
-  "computeCovarianceMatrix", "computeCovarianceMatrix"
+  "computeCovarianceMatrix", "computeCovarianceMatrix",
+  # multivariate effect size computation
+  "varianceCovarianceMatrixType",
+  "varianceCovarianceMatrixFile",
+  "varianceCovarianceMatrixCorrelationMatrix",
+  "varianceCovarianceMatrixSubcluster",
+  "varianceCovarianceMatrixCluster",
+  "varianceCovarianceMatrixForcePositiveDefiniteness",
+  "varianceCovarianceMatrixCheckPositiveDefiniteness",
+  "varianceCovarianceMatrixCorrelationMatrix",
+  "varianceCovarianceMatrixConstruct",
+  "varianceCovarianceMatrixConstructType",
+  "varianceCovarianceMatrixTime1",
+  "varianceCovarianceMatrixTime2",
+  "varianceCovarianceMatrixGroup1",
+  "varianceCovarianceMatrixGroup1",
+  "varianceCovarianceMatrixGroupSize1",
+  "varianceCovarianceMatrixGroupSize2",
+  "varianceCovarianceMatrixConstructCorrelationMatrix",
+  "varianceCovarianceMatrixConstructCorrelationMatrixValue",
+  "varianceCovarianceMatrixConstructCorrelationMatrixFilePath",
+  "varianceCovarianceMatrixConstructTypeCorrelationMatrix",
+  "varianceCovarianceMatrixConstructTypeCorrelationMatrixValue",
+  "varianceCovarianceMatrixConstructTypeCorrelationMatrixFilePath",
+  "varianceCovarianceMatrixTimeLag1Correlation"
 )
 .maForestPlotDependencies <- c(
   # do not forget to add variable carrying options to the .maDataPlottingDependencies
@@ -183,15 +207,25 @@ ClassicalMetaAnalysis <- function(jaspResults, dataset = NULL, options, ...) {
 
   .hasErrors(
     dataset              = dataset,
-    type                 = c("infinity", "observations", "variance"),
+    type                 = c("infinity", "observations"),
     all.target           = c(
       options[["effectSize"]],
-      options[["effectSizeStandardError"]],
-      options[["predictors"]][options[["predictors.types"]] == "scale"]
+      options[["effectSizeStandardError"]]
     ),
     observations.amount  = "< 2",
     exitAnalysisIfErrors = TRUE)
 
+  # do not check effect sizes / standard errors for 0 variance
+  otherVariable <- options[["predictors"]][options[["predictors.types"]] == "scale"]
+  if (length(otherVariable) > 0) {
+    .hasErrors(
+      dataset              = dataset,
+      type                 = c("infinity", "observations", "variance"),
+      all.target           = otherVariable,
+      observations.amount  = "< 2",
+      exitAnalysisIfErrors = TRUE)
+  }
+
   if (length(options[["effectSizeModelTerms"]]) > 0)
     .hasErrors(
       dataset              = dataset,
diff --git a/R/classicalmetaanalysiscommon.R b/R/classicalmetaanalysiscommon.R
index c33bf320..90555eef 100644
--- a/R/classicalmetaanalysiscommon.R
+++ b/R/classicalmetaanalysiscommon.R
@@ -22,8 +22,6 @@
 
 
 # TODO:
-# Forest plot
-# - allow aggregation of studies by a factor (then show simple REML aggregation within and overlaying shaded estimates)
 # AIC/BIC Model-averaging
 # Diagnostics
 # - model re-run on presence of influential cases
@@ -97,11 +95,14 @@
   if (options[["diagnosticsResidualFunnel"]])
     .maResidualFunnelPlot(jaspResults, options)
 
-
   # additional
   if (options[["showMetaforRCode"]])
     .maShowMetaforRCode(jaspResults, options)
 
+  # export the variance-covariance matrix if requested
+  if (.maIsMultilevelMultivariate(options) && options[["varianceCovarianceMatrixSaveComputedVarianceCovarianceMatrix"]] != "")
+    .mammExportVarianceCovarianceMatrix(dataset, options)
+
   return()
 }
 
@@ -154,10 +155,11 @@
       subgroupIndx  <- dataset[[options[["subgroup"]]]] == subgroupLevel
       subgroupData  <- droplevels(dataset[subgroupIndx, ])
 
-      # forward NAs information
+      # forward NAs information and additional attributes
       tempNasIds    <- attr(dataset, "NasIds")[!attr(dataset, "NasIds")]
-      attr(subgroupData, "NAs")    <- sum(tempNasIds[subgroupIndx])
-      attr(subgroupData, "NasIds") <- tempNasIds[subgroupIndx]
+      attr(subgroupData, "NAs")          <- sum(tempNasIds[subgroupIndx])
+      attr(subgroupData, "NasIds")       <- tempNasIds[subgroupIndx]
+      attr(subgroupData, "subgroupIndx") <- subgroupIndx
 
       # fit the model
       fitOutput[[paste0("subgroup", subgroupLevel)]] <- .maFitModelFun(subgroupData, options, subgroupName = as.character(subgroupLevel))
@@ -178,16 +180,17 @@
 
   # specify the effect size and outcome
   if (options[["module"]] == "metaAnalysis") {
+    # specify the univariate input
     rmaInput <- list(
       yi   = as.name(options[["effectSize"]]),
       sei  = as.name(options[["effectSizeStandardError"]]),
       data = dataset
     )
   } else if (options[["module"]] == "metaAnalysisMultilevelMultivariate") {
-    # TODO: extend to covariance matrices
+    # specify the multivariate input
     rmaInput <- list(
       yi   = as.name(options[["effectSize"]]),
-      V    = as.name("samplingVariance"), # precomputed on data load
+      V    = if (.mammVarianceCovarianceMatrixReady(options)) .mammGetVarianceCovarianceMatrix(dataset, options) else as.name("samplingVariance"),
       data = dataset
     )
   }
@@ -757,6 +760,14 @@
     )
   }
 
+  # add multivariate settings notes
+  if (.maIsMultilevelMultivariate(options)) {
+    multivariateReadyNotes <- attr(.mammVarianceCovarianceMatrixReady(options), "messages")
+    for (i in seq_along(multivariateReadyNotes)) {
+      testsTable$addFootnote(multivariateReadyNotes[i])
+    }
+  }
+
   # bind and clean rows
   tests <- .maSafeRbind(tests)
   tests <- .maSafeOrderAndSimplify(tests, "test", options)
@@ -3316,10 +3327,15 @@
       data = as.name("dataset")
     )
   } else if (options[["module"]] == "metaAnalysisMultilevelMultivariate") {
-    # TODO: extend to covariance matrices
+
+    if (.mammVarianceCovarianceMatrixReady(options)) {
+      vcalcInput <-.mammGetVarianceCovarianceMatrix(NULL, options, returnCall = TRUE)
+      vcalcInput$data <- as.name("dataset")
+    }
+
     rmaInput <- list(
       yi   = as.name(options[["effectSize"]]),
-      V    = paste0(options[["effectSizeStandardError"]], "^2"), # precomputed on data load
+      V    = if (.mammVarianceCovarianceMatrixReady(options)) "effectSizeVarianceCovarianceMatrix" else paste0(options[["effectSizeStandardError"]], "^2"),
       data = as.name("dataset")
     )
   }
@@ -3401,6 +3417,19 @@
     fit <- paste0("fit <- rma(\n\t", paste(names(rmaInput), "=", rmaInput, collapse = ",\n\t"), "\n)\n")
   }
 
+  if (.maIsMultilevelMultivariate(options) &&  .mammVarianceCovarianceMatrixReady(options)) {
+    if (options[["varianceCovarianceMatrixType"]] == "precomputed") {
+      fit <- paste0(
+        paste0("effectSizeVarianceCovarianceMatrix <- ", vcalcInput[["file"]], "\n"), "\n",
+        fit
+      )
+    } else {
+      fit <- paste0(
+        paste0("effectSizeVarianceCovarianceMatrix <- vcalc(\n\t", paste(names(vcalcInput), "=", vcalcInput, collapse = ",\n\t"), "\n)\n"), "\n",
+        fit
+      )
+    }
+  }
 
   # add clustering if specified
   if (options[["clustering"]] != "") {
diff --git a/R/classicalmetaanalysismultilevelmultivariate.R b/R/classicalmetaanalysismultilevelmultivariate.R
index 05af28e9..a3a65268 100644
--- a/R/classicalmetaanalysismultilevelmultivariate.R
+++ b/R/classicalmetaanalysismultilevelmultivariate.R
@@ -39,11 +39,15 @@ ClassicalMetaAnalysisMultilevelMultivariate <- function(jaspResults, dataset = N
   # random effects variables
   randomVariables <- .mammExtractRandomVariableNames(options)
 
+  # variance-covariance variables
+  varianceCovarianceVariables <- .mammExtractVarianceCovarianceMatrixNames(options)
+
   # omit NAs
   omitOnVariables <- c(
     options[["effectSize"]],
     options[["effectSizeStandardError"]],
     unlist(randomVariables),
+    varianceCovarianceVariables,
     if (options[["clustering"]] != "") options[["clustering"]],
     if (options[["subgroup"]] != "")   options[["subgroup"]],
     if (length(predictorsNominal) > 0) predictorsNominal,
@@ -65,16 +69,29 @@ ClassicalMetaAnalysisMultilevelMultivariate <- function(jaspResults, dataset = N
 
   .hasErrors(
     dataset              = dataset,
-    type                 = c("infinity", "observations", "variance"),
+    type                 = c("infinity", "observations"),
     all.target           = c(
       options[["effectSize"]],
-      options[["effectSizeStandardError"]],
-      options[["predictors"]][options[["predictors.types"]] == "scale"],
-      c(randomVariables$scale, randomVariables$ordinal)
+      options[["effectSizeStandardError"]]
     ),
     observations.amount  = "< 2",
     exitAnalysisIfErrors = TRUE)
 
+  # do not check effect sizes / standard errors for 0 variance
+  otherVariable <- c(
+    options[["predictors"]][options[["predictors.types"]] == "scale"],
+    c(randomVariables$scale, randomVariables$ordinal)
+  )
+  if (length(otherVariable) > 0) {
+    .hasErrors(
+      dataset              = dataset,
+      type                 = c("infinity", "observations", "variance"),
+      all.target           = otherVariable,
+      observations.amount  = "< 2",
+      exitAnalysisIfErrors = TRUE)
+  }
+
+
   .hasErrors(
     dataset              = dataset,
     type                 = c("modelInteractions"),
@@ -268,7 +285,352 @@ ClassicalMetaAnalysisMultilevelMultivariate <- function(jaspResults, dataset = N
 
   return(randomFormulas)
 }
-.mammExtractRandomVariableNames  <- function(options) {
+.mammGetVarianceCovarianceMatrix    <- function(dataset, options, returnCall = FALSE) {
+
+  if (options[["varianceCovarianceMatrixType"]] == "precomputed") {
+
+    # load a pre-computed correlation matrix and check the dimensions against the dataset
+    vMatrixFileName <- options[["varianceCovarianceMatrixFile"]]
+
+    if (vMatrixFileName == "")
+      return()
+
+    vMatrix <- try(as.matrix(read.csv(file = vMatrixFileName, header = FALSE)))
+
+    if (returnCall)
+      return(list(file = vMatrixFileName))
+
+    if (inherits(vMatrix, "try-error"))
+      .quitAnalysis(gettextf("Error reading the variance-covariance matrix file: %1$s", vMatrix))
+
+    # if there is only one column, try csv2 (indicates different decimals enconding)
+    if (ncol(vMatrix) == 1)
+      vMatrix <- try(as.matrix(read.csv2(file = vMatrixFileName, header = FALSE)))
+
+    if (inherits(vMatrix, "try-error"))
+      .quitAnalysis(gettextf("Error reading the variance-covariance matrix file: %1$s", vMatrix))
+
+    if (nrow(vMatrix) != ncol(vMatrix))
+      .quitAnalysis(gettextf("The variance-covariance matrix must be square. The number of rows (%1$i) does not match the number of columns (%2$i).",
+                             nrow(vMatrix), ncol(vMatrix)))
+
+    # extract the subgroup relevant section if subgroups are specified
+    if (options[["subgroup"]] != "" && !is.null(attr(dataset, "subgroupIndx")))
+      vMatrix <- vMatrix[attr(dataset, "subgroupIndx"), attr(dataset, "subgroupIndx")]
+
+    if (nrow(vMatrix) != nrow(dataset))
+      .quitAnalysis(gettextf("The variance-covariance matrix must match the dimensions of the data set. The number of the matrix rows/columns (%1$i) does not match the number of observations in the dataset (%2$i).",
+                             nrow(vMatrix), nrow(dataset)))
+
+
+    if (options[["varianceCovarianceMatrixForcePositiveDefiniteness"]]) {
+      isPD <- !any(eigen(vMatrix, symmetric = TRUE, only.values = TRUE)$values <= .Machine$double.eps)
+      if (!isPD) {
+        # only force into PD if is not PD
+        vMatrix <- Matrix::nearPD(vMatrix, corr = FALSE)
+      }
+    }
+
+    if (options[["varianceCovarianceMatrixCheckPositiveDefiniteness"]]) {
+      isPD <- !any(eigen(vMatrix, symmetric = TRUE, only.values = TRUE)$values <= .Machine$double.eps)
+      if (!isPD)
+        .quitAnalysis(gettextf("The supplied variance-covariance matrix is not positive definite."))
+    }
+
+    return(vMatrix)
+
+  } else if (options[["varianceCovarianceMatrixType"]] == "correlationMatrix") {
+
+    vcalcCal <- list(
+      vi      = as.name("samplingVariance"),
+      rvars   = as.call(c(quote(c), lapply(options[["varianceCovarianceMatrixCorrelationMatrix"]], as.name))),
+      cluster = as.name(options[["varianceCovarianceMatrixCluster"]]),
+      data    = dataset,
+      checkpd = options[["varianceCovarianceMatrixCheckPositiveDefiniteness"]],
+      nearpd  = options[["varianceCovarianceMatrixForcePositiveDefiniteness"]]
+    )
+
+    if (options[["varianceCovarianceMatrixSubcluster"]] != "")
+      vcalcCal$subgroup <- as.name(options[["varianceCovarianceMatrixSubcluster"]])
+
+    if (returnCall) {
+      return(vcalcCal)
+    }
+    vMatrix <- try(do.call(metafor::vcalc, vcalcCal))
+
+    # try cleaning the error message before returning
+    if (jaspBase::isTryError(vMatrix)) {
+      vMatrix <- gsub("'rvars'",  "'Correlation Matrix'", vMatrix)
+      .quitAnalysis(gettextf("Error computing the variance-covariance matrix: %1$s", vMatrix))
+    }
+
+    return(vMatrix)
+
+  } else if (options[["varianceCovarianceMatrixType"]] == "constructsGroupsTimes") {
+
+    vcalcCal <- list(
+      vi      = as.name("samplingVariance"),
+      data    = dataset
+    )
+
+    # add clusters and subclusters
+    if (options [["varianceCovarianceMatrixSubcluster"]] != "") {
+      vcalcCal$subgroup <- as.name(options[["varianceCovarianceMatrixSubcluster"]])
+    }
+    if (options [["varianceCovarianceMatrixCluster"]] != "") {
+      vcalcCal$cluster  <- as.name(options[["varianceCovarianceMatrixCluster"]])
+    }
+
+    # resolve construsts and subconstructs
+    if (options[["varianceCovarianceMatrixConstruct"]] != "") {
+
+      vcalcCal$obs <- as.name(options[["varianceCovarianceMatrixConstruct"]])
+
+      # get the rho
+      if (options[["varianceCovarianceMatrixConstructCorrelationMatrix"]] == "commonCorrelation") {
+        # set a single value
+        rho1 <- options[["varianceCovarianceMatrixConstructCorrelationMatrixValue"]]
+      } else if(options[["varianceCovarianceMatrixConstructCorrelationMatrix"]] == "correlationMatrix") {
+        # load a complete matrix
+        rho1File <- options[["varianceCovarianceMatrixConstructCorrelationMatrixFilePath"]]
+
+        if (rho1File != "") {
+          if (returnCall) {
+            rho1 <- rho1File
+          } else  {
+            rho1 <- try(as.matrix(read.csv(file = rho1File, row.names = 1)))
+
+            if (inherits(rho1, "try-error"))
+              .quitAnalysis(gettextf("Error reading the variance-covariance matrix file: %1$s", rho1))
+
+            # if there is only one column, try csv2 (indicates different decimals enconding)
+            if (ncol(rho1) == 1)
+              rho1 <- try(as.matrix(read.csv2(file = rho1File, row.names = 1)))
+
+            if (inherits(rho1, "try-error"))
+              .quitAnalysis(gettextf("Error reading the variance-covariance matrix file: %1$s", rho1))
+
+            if (nrow(rho1) != ncol(rho1))
+              .quitAnalysis(gettextf("The variance-covariance matrix must be square. The number of rows (%1$i) does not match the number of columns (%2$i).",
+                                     nrow(rho1), ncol(rho1)))
+          }
+        } else {
+          rho1 <- 0
+        }
+      }
+      # set rho
+      vcalcCal$rho <- rho1
+    }
+
+    if (options[["varianceCovarianceMatrixConstructType"]] != "") {
+
+      vcalcCal$type <- as.name(options[["varianceCovarianceMatrixConstructType"]])
+
+      if (options[["varianceCovarianceMatrixConstructTypeCorrelationMatrix"]] == "commonCorrelation") {
+        # set a single value
+        rho2 <- options[["varianceCovarianceMatrixConstructTypeCorrelationMatrixValue"]]
+      } else if(options[["varianceCovarianceMatrixConstructTypeCorrelationMatrix"]] == "correlationMatrix") {
+        # load a complete matrix
+        rho2File <- options[["varianceCovarianceMatrixConstructTypeCorrelationMatrixFilePath"]]
+
+        if (rho2File != "") {
+          if (returnCall) {
+            rho2 <- rho2File
+          } else  {
+            rho2 <- try(as.matrix(read.csv(file = rho2File, row.names = 1)))
+
+            if (inherits(rho2, "try-error"))
+              .quitAnalysis(gettextf("Error reading the variance-covariance matrix file: %1$s", rho2))
+
+            # if there is only one column, try csv2 (indicates different decimals enconding)
+            if (ncol(rho2) == 1)
+              rho2 <- try(as.matrix(read.csv2(file = rho2File, row.names = 1)))
+
+            if (inherits(rho2, "try-error"))
+              .quitAnalysis(gettextf("Error reading the variance-covariance matrix file: %1$s", rho2))
+
+            if (nrow(rho2) != ncol(rho2))
+              .quitAnalysis(gettextf("The variance-covariance matrix must be square. The number of rows (%1$i) does not match the number of columns (%2$i).",
+                                     nrow(rho2), ncol(rho2)))
+          }
+        } else {
+          rho2 <- 0
+        }
+      }
+      # add rho
+      if (!is.null(vcalcCal$rho)) {
+        vcalcCal$rho <- list(vcalcCal$rho, rho2)
+      } else {
+        vcalcCal$rho <- rho2
+      }
+    }
+
+    # add time lags
+    if (options[["varianceCovarianceMatrixTime1"]] != "") {
+
+      vcalcCal$time1 <- as.name(options[["varianceCovarianceMatrixTime1"]])
+
+      if (options[["varianceCovarianceMatrixTime2"]] != "") {
+        vcalcCal$time2 <- as.name(options[["varianceCovarianceMatrixTime2"]])
+      }
+
+      vcalcCal$phi <- options[["varianceCovarianceMatrixTimeLag1Correlation"]]
+    }
+
+    # add weights
+    if (options[["varianceCovarianceMatrixGroup1"]] != "") {
+      vcalcCal$grp1 <- as.name(options[["varianceCovarianceMatrixGroup1"]])
+      if (options[["varianceCovarianceMatrixGroupSize1"]] != "") {
+        vcalcCal$w1 <- as.name(options[["varianceCovarianceMatrixGroupSize1"]])
+      }
+
+      if (options[["varianceCovarianceMatrixGroup1"]] != "") {
+        vcalcCal$grp2 <- as.name(options[["varianceCovarianceMatrixGroup2"]])
+        if (options[["varianceCovarianceMatrixGroupSize2"]] != "") {
+          vcalcCal$w2 <- as.name(options[["varianceCovarianceMatrixGroupSize2"]])
+        }
+      }
+    }
+
+    # make the correlation matrix call
+    if (returnCall) {
+      return(vcalcCal)
+    }
+    vMatrix <- try(do.call(metafor::vcalc, vcalcCal))
+
+    # try cleaning the error message before returning
+    if (jaspBase::isTryError(vMatrix)) {
+      if (grepl("'rho'", vMatrix) && grepl("'obs'", vMatrix)) {
+        vMatrix <- gsub("'rho'",  "'Construct Correlation Matrix'", vMatrix)
+        vMatrix <- gsub("'obs'",  "'Construct Correlation Matrix'", vMatrix)
+      }
+      if (grepl("'rho'", vMatrix) && grepl("'type'", vMatrix)) {
+        vMatrix <- gsub("'rho'",  "'Construct Type Correlation Matrix'", vMatrix)
+        vMatrix <- gsub("'type'", "'Construct Type Correlation Matrix'", vMatrix)
+      }
+      if (grepl("'phi'", vMatrix)) {
+        vMatrix <- gsub("'time1'",  "'Time 1'", vMatrix)
+        vMatrix <- gsub("'time2'",  "'Time 2'", vMatrix)
+        vMatrix <- gsub("'phi'",  "'Time lag 1 correlation'", vMatrix)
+      }
+      if (grepl("'grp", vMatrix)) {
+        vMatrix <- gsub("'grp1'",  "'Group 1'", vMatrix)
+        vMatrix <- gsub("'grp2'",  "'Group 2'", vMatrix)
+      }
+      if (grepl("'w", vMatrix)) {
+        vMatrix <- gsub("'w1'",  "'Group Size 1'", vMatrix)
+        vMatrix <- gsub("'w2'",  "'Group Size 2'", vMatrix)
+      }
+      .quitAnalysis(gettextf("Error computing the variance-covariance matrix: %1$s", vMatrix))
+    }
+
+    return(vMatrix)
+  }
+}
+.mammVarianceCovarianceMatrixReady  <- function(options) {
+
+  ready    <- FALSE
+  messages <- c()
+
+  if (options[["varianceCovarianceMatrixType"]] == "precomputed") {
+
+    # load a pre-computed correlation matrix and check the dimensions against the dataset
+    vMatrixFileName <- options[["varianceCovarianceMatrixFile"]]
+
+    if (vMatrixFileName == "") {
+      messages <- c(messages, gettext("Please provide a file with the precomputed variance-covariance matrix."))
+    } else {
+      ready <- TRUE
+    }
+
+
+  } else if (options[["varianceCovarianceMatrixType"]] == "correlationMatrix") {
+
+    if (options[["varianceCovarianceMatrixCluster"]] == "" && length(options[["varianceCovarianceMatrixCorrelationMatrix"]]) > 0) {
+      messages <- c(messages, gettext("Please provide a clustering variable for the correlation matrix."))
+    } else if (options[["varianceCovarianceMatrixCluster"]] != "" && length(options[["varianceCovarianceMatrixCorrelationMatrix"]]) == 0) {
+      messages <- c(messages, gettext("Please provide correlation matrix variables."))
+    } else if (options[["varianceCovarianceMatrixCluster"]] != "" && length(options[["varianceCovarianceMatrixCorrelationMatrix"]]) > 0) {
+      ready <- TRUE
+    }
+
+  } else if (options[["varianceCovarianceMatrixType"]] == "constructsGroupsTimes") {
+
+    if (options[["varianceCovarianceMatrixCluster"]] == "" &&
+        (options[["varianceCovarianceMatrixSubcluster"]]    != "" ||
+         options[["varianceCovarianceMatrixConstructType"]] != "" ||
+         options[["varianceCovarianceMatrixConstruct"]]     != "" ||
+         options[["varianceCovarianceMatrixTime1"]]         != "" ||
+         options[["varianceCovarianceMatrixGroup1"]]        != "")) {
+      messages <- c(messages, gettext("Please provide a clustering variable for the correlation matrix."))
+    }
+
+    if (options[["varianceCovarianceMatrixCluster"]] != "" && options[["varianceCovarianceMatrixConstructType"]] != "") {
+      if (options[["varianceCovarianceMatrixConstructTypeCorrelationMatrix"]] == "commonCorrelation") {
+        ready <- TRUE
+        if (options[["varianceCovarianceMatrixConstructTypeCorrelationMatrixValue"]] == 0) {
+          # ready but with odd settings
+          messages <- c(messages, gettext("The value of the correlation between construct types is set to 0. This corresponds to no adjustment for dependency between construct types."))
+        }
+      } else if (options[["varianceCovarianceMatrixConstructTypeCorrelationMatrix"]] == "correlationMatrix") {
+        if (options[["varianceCovarianceMatrixConstructTypeCorrelationMatrixFilePath"]] == "") {
+          messages <- c(messages, gettext("Please provide a file with the correlation matrix for the construct types."))
+        } else {
+          ready <- TRUE
+        }
+      }
+    }
+    if (options[["varianceCovarianceMatrixCluster"]] != "" && options[["varianceCovarianceMatrixConstruct"]] != "") {
+      if (options[["varianceCovarianceMatrixConstructCorrelationMatrix"]] == "commonCorrelation") {
+        ready <- TRUE
+        if (options[["varianceCovarianceMatrixConstructCorrelationMatrixValue"]] == 0) {
+          # ready but with odd settings
+          messages <- c(messages, gettext("The value of the correlation between construct types is set to 0. This corresponds to no adjustment for dependency between constructs."))
+        }
+      } else if (options[["varianceCovarianceMatrixConstructCorrelationMatrix"]] == "correlationMatrix") {
+        if (options[["varianceCovarianceMatrixConstructCorrelationMatrixFilePath"]] == "") {
+          messages <- c(messages, gettext("Please provide a file with the correlation matrix for the construct."))
+        } else {
+          ready <- TRUE
+        }
+      }
+    }
+    if (options[["varianceCovarianceMatrixCluster"]] != "" && options[["varianceCovarianceMatrixTime1"]] != "") {
+
+      ready <- TRUE
+      if (options[["varianceCovarianceMatrixTimeLag1Correlation"]] == 0) {
+        # ready but with odd settings
+        messages <- c(messages, gettext("The value of the correlation between time lags is set to 0. This corresponds to no adjustment for dependency between time lags."))
+      }
+    }
+    if (options[["varianceCovarianceMatrixCluster"]] != "" && options[["varianceCovarianceMatrixGroup1"]] != "") {
+
+      ready <- TRUE
+      if (options[["varianceCovarianceMatrixGroupSize1"]] == "") {
+        # ready with default settings
+        messages <- c(messages, gettext("The 'Group Size' for specifying dependency due to common group comparison is unspecified. As default, the groups are assumed to be equal."))
+      }
+    }
+
+  }
+
+  attr(ready, "messages") <- messages
+  return(ready)
+}
+.mammExportVarianceCovarianceMatrix <- function(dataset, options) {
+
+  if (!.mammVarianceCovarianceMatrixReady(options))
+    return()
+
+  V <- .mammGetVarianceCovarianceMatrix(dataset, options)
+  try(write.table(V, file = options[["varianceCovarianceMatrixSaveComputedVarianceCovarianceMatrix"]], sep = ",", row.names = FALSE, col.names = FALSE))
+
+  if (jaspBase::isTryError(V))
+    .quitAnalysis(gettextf("Error writing the variance-covariance matrix file: %1$s", V))
+
+  return()
+}
+.mammExtractRandomVariableNames           <- function(options) {
 
   if (length(options[["randomEffects"]]) == 0)
     return(NULL)
@@ -351,6 +713,38 @@ ClassicalMetaAnalysisMultilevelMultivariate <- function(jaspResults, dataset = N
     ordinal = if (length(variablesOrdinal) != 0) variablesOrdinal
   ))
 }
+.mammExtractVarianceCovarianceMatrixNames <- function(options) {
+
+  if (options[["varianceCovarianceMatrixType"]] == "precomputed") {
+    return()
+  } else if (options[["varianceCovarianceMatrixType"]] == "correlationMatrix") {
+
+    # varianceCovarianceMatrixCorrelationMatrix can contain NA's as it is only lower triangular
+    variableNames <- c(
+      options[["varianceCovarianceMatrixCluster"]],
+      options[["varianceCovarianceMatrixSubluster"]]
+    )
+    variableNames <- variableNames[variableNames != ""]
+    return(variableNames)
+
+  } else if (options[["varianceCovarianceMatrixType"]] == "constructsGroupsTimes") {
+
+    variableNames <- c(
+      options[["varianceCovarianceMatrixCluster"]],
+      options[["varianceCovarianceMatrixSubcluster"]],
+      options[["varianceCovarianceMatrixConstructType"]],
+      options[["varianceCovarianceMatrixConstruct"]],
+      options[["varianceCovarianceMatrixTime1"]],
+      options[["varianceCovarianceMatrixTime2"]],
+      options[["varianceCovarianceMatrixGroup1"]],
+      options[["varianceCovarianceMatrixGroup2"]],
+      options[["varianceCovarianceMatrixGroupSize1"]],
+      options[["varianceCovarianceMatrixGroupSize2"]]
+    )
+    variableNames <- variableNames[variableNames != ""]
+    return(variableNames)
+  }
+}
 .mammRandomEstimatesTable        <- function(jaspResults, options) {
 
   # obtain the overall container
diff --git a/inst/qml/ClassicalMetaAnalysisMultilevelMultivariate.qml b/inst/qml/ClassicalMetaAnalysisMultilevelMultivariate.qml
index 01cf8efa..d3b096d4 100644
--- a/inst/qml/ClassicalMetaAnalysisMultilevelMultivariate.qml
+++ b/inst/qml/ClassicalMetaAnalysisMultilevelMultivariate.qml
@@ -115,6 +115,340 @@ Form
 		}
 	}
 
+	Section
+	{	
+		title:		qsTr("Effect Size Variance-Covariance Matrix")
+		expanded:	false
+		info: qsTr("Options for specifying the approximate variance-covariance matrix of the effect sizes. This matrix is used to account for the correlation between effect sizes when they are not independent.")
+
+		RadioButtonGroup
+		{
+			id:			varianceCovarianceMatrixType
+			name:		"varianceCovarianceMatrixType"
+			title:		qsTr("Type")
+			info: qsTr("Type of variance-covariance matrix input method.")
+			columns:	3
+
+			RadioButton
+			{
+				id:			varianceCovarianceMatrixCorrelationMatrix
+				name:		"correlationMatrix"
+				label:		qsTr("Correlation matrix")
+				checked:	true
+				info: qsTr("Use a list of variables to specify the correlation matrix of studies corresponding to the same cluster. Corresponds to the `rvars` option in the metafor's 'vcalc' function.")
+			}
+
+			RadioButton
+			{
+				id:			varianceCovarianceMatrixConstructsGroupsTimes
+				name:		"constructsGroupsTimes"
+				label:		qsTr("Constructs, groups, and times")
+				info: qsTr("Specify constructs, groups, and times of the measurement to specify the correlation matrix of studies corresponding to the same cluster.")
+			}
+
+			RadioButton
+			{
+				id:			varianceCovarianceMatrixPrecomputed
+				name:		"precomputed"
+				label:		qsTr("Precomputed")
+				info: qsTr("Load a csv file containing the precomputed variance-covariance matrix.")
+			}
+		}
+		
+		FileSelector
+		{
+			name:		"varianceCovarianceMatrixFile"
+			label:		qsTr("Effect size variance-covariance matrix file")
+			visible:	varianceCovarianceMatrixPrecomputed.checked
+			filter:		"*.csv"
+			save:		false
+			info: qsTr("CSV file containing the precomputed effect size variance-covariance matrix. The matrix needs to match the dimensions of the data set and cannot contain any other variables or names.")
+		}
+
+
+		VariablesForm
+		{
+			removeInvisibles:	true
+			preferredHeight: 	(varianceCovarianceMatrixCorrelationMatrix.checked ? 250 : 525) * preferencesModel.uiScale
+			visible:			varianceCovarianceMatrixConstructsGroupsTimes.checked || varianceCovarianceMatrixCorrelationMatrix.checked	
+
+			AvailableVariablesList
+			{
+				name:				"varianceCovarianceMatrixAllVariables"
+			}
+				
+			AssignedVariablesList
+			{
+				name:				"varianceCovarianceMatrixCorrelationMatrix"
+				title:				qsTr("Correlation Matrix")
+				allowedColumns:		["scale"]
+				allowTypeChange:	true
+				singleVariable:		false
+				visible:			varianceCovarianceMatrixCorrelationMatrix.checked
+				property bool active:	varianceCovarianceMatrixCorrelationMatrix.checked
+				onActiveChanged: 		if (!active && count > 0) itemDoubleClicked(0)
+				info: qsTr("Variable specifying the correlation between the individual estimates within cluster. The column order of the variable correspond to the row order of effect size estimates within cluster. Only the lower triangle needs to be specified. Corresponds to the `rvars` option in the metafor's 'vcalc' function.")
+			}
+
+			AssignedVariablesList
+			{
+				name:				"varianceCovarianceMatrixConstruct"
+				id:					varianceCovarianceMatrixConstruct
+				title:				qsTr("Construct")
+				allowedColumns:		["nominal"]
+				allowTypeChange:	true
+				singleVariable:		true
+				visible:			varianceCovarianceMatrixConstructsGroupsTimes.checked
+				property bool active:	varianceCovarianceMatrixConstructsGroupsTimes.checked
+				onActiveChanged: 		if (!active && count > 0) itemDoubleClicked(0)
+				info: qsTr("Variable specifying the construct measured by the effect size. Corresponds to the `obs` option in the metafor's 'vcalc' function.")
+			}
+
+			AssignedVariablesList
+			{
+				name:				"varianceCovarianceMatrixConstructType"
+				id:					varianceCovarianceMatrixConstructType
+				title:				qsTr("Construct Type")
+				allowedColumns:		["nominal"]
+				allowTypeChange:	true
+				singleVariable:		true
+				visible:			varianceCovarianceMatrixConstructsGroupsTimes.checked
+				property bool active:	varianceCovarianceMatrixConstructsGroupsTimes.checked
+				onActiveChanged: 		if (!active && count > 0) itemDoubleClicked(0)
+				info: qsTr("Variable specifying the type of construct measured by the effect size. Construct Type corresponds to a higher level grouping of the Constructs. Corresponds to the `type` option in the metafor's 'vcalc' function.")
+			}
+
+			AssignedVariablesList
+			{
+				name:				"varianceCovarianceMatrixTime1"
+				id:					varianceCovarianceMatrixTime1
+				title:				qsTr("Time 1")
+				allowedColumns:		["nominal"]
+				allowTypeChange:	true
+				singleVariable:		true
+				visible:			varianceCovarianceMatrixConstructsGroupsTimes.checked
+				property bool active:	varianceCovarianceMatrixConstructsGroupsTimes.checked
+				onActiveChanged: 		if (!active && count > 0) itemDoubleClicked(0)
+				info: qsTr("Variable specifying the time point of the effect size measurement. In case multiple time points are specified, the first time corresponds to the first codntion.")
+			}			
+
+			AssignedVariablesList
+			{
+				name:				"varianceCovarianceMatrixTime2"
+				enabled:			varianceCovarianceMatrixTime1.count != 0
+				title:				qsTr("Time 2")
+				allowedColumns:		["nominal"]
+				allowTypeChange:	true
+				singleVariable:		true
+				visible:			varianceCovarianceMatrixConstructsGroupsTimes.checked
+				property bool active:	varianceCovarianceMatrixConstructsGroupsTimes.checked
+				onActiveChanged: 		if (!active && count > 0) itemDoubleClicked(0)
+				info: qsTr("Variable specifying the time point of the effect size measurement. In case multiple time points are specified, the second time corresponds to the second condition.")
+			}			
+
+			AssignedVariablesList
+			{
+				name:				"varianceCovarianceMatrixGroup1"
+				id:					varianceCovarianceMatrixGroup1
+				title:				qsTr("Group 1")
+				allowedColumns:		["nominal"]
+				allowTypeChange:	true
+				singleVariable:		true
+				visible:			varianceCovarianceMatrixConstructsGroupsTimes.checked
+				property bool active:	varianceCovarianceMatrixConstructsGroupsTimes.checked
+				onActiveChanged: 		if (!active && count > 0) itemDoubleClicked(0)
+				info: qsTr("Variable to specify the group of the first condition when the observed effect sizes or outcomes represent contrasts between two conditions. Corresponds to the `grp1` argument in the metafor's 'vcalc' function.")
+			}
+
+			AssignedVariablesList
+			{
+				name:				"varianceCovarianceMatrixGroup2"
+				id:					varianceCovarianceMatrixGroup2
+				enabled:			varianceCovarianceMatrixGroup1.count != 0
+				title:				qsTr("Group 2")
+				allowedColumns:		["nominal"]
+				allowTypeChange:	true
+				singleVariable:		true
+				visible:			varianceCovarianceMatrixConstructsGroupsTimes.checked
+				property bool active:	varianceCovarianceMatrixConstructsGroupsTimes.checked
+				onActiveChanged: 		if (!active && count > 0) itemDoubleClicked(0)
+				info: qsTr("Variable to specify the group of the second condition when the observed effect sizes or outcomes represent contrasts between two conditions. Corresponds to the `grp1` argument in the metafor's 'vcalc' function.")
+			}
+
+			AssignedVariablesList
+			{
+				name:				"varianceCovarianceMatrixGroupSize1"
+				enabled:			varianceCovarianceMatrixGroup1.count != 0
+				title:				qsTr("Group Size 1")
+				allowedColumns:		["scale"]
+				allowTypeChange:	true
+				singleVariable:		true
+				visible:			varianceCovarianceMatrixConstructsGroupsTimes.checked
+				property bool active:	varianceCovarianceMatrixConstructsGroupsTimes.checked
+				onActiveChanged: 		if (!active && count > 0) itemDoubleClicked(0)
+				info: qsTr("Variable to specify the size of the group (or more generally, the inverse-sampling variance weight) of the first condition when the observed effect sizes or outcomes represent contrasts between two conditions. Corresponds to the `w1` argument in the metafor's 'vcalc' function.")
+			}
+
+			AssignedVariablesList
+			{
+				name:				"varianceCovarianceMatrixGroupSize2"
+				enabled:			varianceCovarianceMatrixGroup2.count != 0
+				title:				qsTr("Group Size 2")
+				allowedColumns:		["scale"]
+				allowTypeChange:	true
+				singleVariable:		true
+				visible:			varianceCovarianceMatrixConstructsGroupsTimes.checked
+				property bool active:	varianceCovarianceMatrixConstructsGroupsTimes.checked
+				onActiveChanged: 		if (!active && count > 0) itemDoubleClicked(0)
+				info: qsTr("Optional numeric vector to specify the size of the group (or more generally, the inverse-sampling variance weight) of the second condition when the observed effect sizes or outcomes represent contrasts between two conditions. Corresponds to the `w2` argument in the metafor's 'vcalc' function.")
+			}
+
+			AssignedVariablesList
+			{
+				name:				"varianceCovarianceMatrixSubcluster"
+				enabled:			varianceCovarianceMatrixCluster.count != 0
+				title:				qsTr("Subcluster")
+				singleVariable:		true
+				allowedColumns:		["nominal"]
+				info: qsTr("Variable specifying additional structure of the subgroups. Effect sizes within the same cluster with different values of the cluster variable are assumed to be independent. Note that this input corresponds to the 'subgroup' option in the metafor's 'vcalc' function and is renamed for differentiation from subgroup analysis. ")
+			}
+
+			AssignedVariablesList
+			{
+				name:				"varianceCovarianceMatrixCluster"
+				id:					varianceCovarianceMatrixCluster
+				title:				qsTr("Cluster")
+				singleVariable:		true
+				allowedColumns:		["nominal"]
+				info: qsTr("Variable specifying clustering of the effect sizes for computing the variance covariance matrix. Effect sizes with different values of the cluster variable are assumed to be independent. Note that this input differs from the 'Clustering' option in the data input which is used to specify cluster-robust standard error. In most cases however, both input should contain the same variable.")
+			}
+		}
+
+		RadioButtonGroup
+		{
+			title:		qsTr("Construct Correlation Matrix")
+			name: 		"varianceCovarianceMatrixConstructCorrelationMatrix"
+			visible:	varianceCovarianceMatrixConstructsGroupsTimes.checked
+			enabled:	varianceCovarianceMatrixConstructsGroupsTimes.checked && varianceCovarianceMatrixConstruct.count != 0
+			columns:	1
+
+			RadioButton
+			{
+				name:		"commonCorrelation"
+				label:		qsTr("Common correlation")
+				info: qsTr("Specify the correlation between the same construct levels.")
+				childrenOnSameRow: true
+				checked: true
+
+				DoubleField
+				{
+					name:				"varianceCovarianceMatrixConstructCorrelationMatrixValue"
+					defaultValue:		0
+					min: 				-1
+					max: 				1
+					inclusive: 			JASP.None
+				}
+			}
+
+			RadioButton
+			{
+				name:		"correlationMatrix"
+				label:		qsTr("Correlation matrix")
+				info: qsTr("CSV file containing the correlation matrix between the constructs levels. The first row and the first column of the file must contain names that map the matrix entries to the construct level names (the names cannot start with a number).")
+				childrenOnSameRow: true
+
+				FileSelector
+				{
+					name:		"varianceCovarianceMatrixConstructCorrelationMatrixFilePath"
+					filter:		"*.csv"
+					save:		false
+				}
+			}
+		}
+
+		RadioButtonGroup
+		{
+			title:		qsTr("Construct Type Correlation Matrix")
+			name: 		"varianceCovarianceMatrixConstructTypeCorrelationMatrix"
+			visible:	varianceCovarianceMatrixConstructsGroupsTimes.checked
+			enabled:	varianceCovarianceMatrixConstructType.count != 0
+			columns:	1
+
+			RadioButton
+			{
+				name:		"commonCorrelation"
+				label:		qsTr("Common correlation")
+				info: qsTr("Specify the correlation between the same construct levels.")
+				childrenOnSameRow: true
+				checked: true
+
+				DoubleField
+				{
+					name:				"varianceCovarianceMatrixConstructTypeCorrelationMatrixValue"
+					defaultValue:		0
+					min: 				-1
+					max: 				1
+					inclusive: 			JASP.None
+				}
+			}
+
+			RadioButton
+			{
+				name:		"correlationMatrix"
+				label:		qsTr("Correlation matrix")
+				info: qsTr("CSV file containing the correlation matrix between the constructs type level. The first row and the first column of the file must contain names that map the matrix entries to the construct level type names (the names cannot start with a number).")
+				childrenOnSameRow: true
+
+				FileSelector
+				{
+					name:		"varianceCovarianceMatrixConstructTypeCorrelationMatrixFilePath"
+					filter:		"*.csv"
+					save:		false
+				}
+			}
+		}
+
+		DoubleField
+		{
+			label:				qsTr("Time lag 1 correlation")
+			visible:			varianceCovarianceMatrixConstructsGroupsTimes.checked
+			enabled:			varianceCovarianceMatrixTime1.count != 0
+			name:				"varianceCovarianceMatrixTimeLag1Correlation"
+			defaultValue:		0
+			min: 				-1
+			max: 				1
+			inclusive: 			JASP.None
+		}
+
+		Group
+		{
+			CheckBox
+			{
+				label:		qsTr("Check positive definiteness")
+				name:		"varianceCovarianceMatrixCheckPositiveDefiniteness"
+				checked: 	true
+				info: qsTr("Check if the variance-covariance matrix is symmetric.")
+			}
+
+			CheckBox
+			{
+				label:		qsTr("Force positive definiteness")
+				name:		"varianceCovarianceMatrixForcePositiveDefiniteness"
+				checked: 	false
+				info: qsTr("Force the variance-covariance matrix to be positive definite. This option shuld be used with caution as non-positive definite matricies often indicate input misspecification.")
+			}
+		}
+
+		FileSelector
+		{
+			label:		qsTr("Save computed variance-covariance matrix")
+			visible:	varianceCovarianceMatrixConstructsGroupsTimes.checked || varianceCovarianceMatrixCorrelationMatrix.checked
+			name:		"varianceCovarianceMatrixSaveComputedVarianceCovarianceMatrix"
+			filter:		"*.csv"
+			save:		true
+		}
+	}
+
 
 	Section
 	{