diff --git a/R/plssem.R b/R/plssem.R
index 0a589d19..ad036718 100644
--- a/R/plssem.R
+++ b/R/plssem.R
@@ -2022,3 +2022,45 @@ checkCSemModel <- function(model, availableVars) {
   return()
 
 }
+
+
+.plsSEMVIFhelper <- function(fit){
+  # Make VIFs into a matrix
+  # Restructure the VIFs into a table.
+  VIFspath <- cSEM::assess(.object = fit,.quality_criterion = 'vif')
+
+  idx <- which(VIFspath$VIF!=0,arr.ind = T)
+
+  if(nrow(idx)!=0){
+    VIFDf <- data.frame(Relation=paste(rownames(VIFspath$VIF)[idx[,'row']],'~',colnames(VIFspath$VIF)[idx[,'col']]),
+                        vif=VIFspath$VIF[cbind(rownames(VIFspath$VIF)[idx[,'row']],colnames(VIFspath$VIF)[idx[,'col']])])
+
+    VIFvector <-setNames(VIFDf$vif, VIFDf$Relation)
+  } else{
+    VIFvector <- NULL
+  }
+  return(VIFvector)
+}
+
+.plsSEMVIFBhelper <- function(fit){
+  VIFsweights <- cSEM::calculateVIFModeB(fit)
+
+  # If there is only one weight, cSEM::calculateVIFModeB() returns NA for that VIF
+  # therefore, replace NAs with 0
+  VIFsweights[is.na(VIFsweights)] <- 0
+
+
+  if(!is.null(VIFsweights)&sum(VIFsweights)!=0){
+    idx <- which(VIFsweights!=0,arr.ind = T)
+
+    VIFBDf <- data.frame(Relation=paste(rownames(VIFsweights)[idx[,'row']],'<~',colnames(VIFsweights)[idx[,'col']]),
+                         vif=VIFsweights[cbind(rownames(VIFsweights)[idx[,'row']],colnames(VIFsweights)[idx[,'col']])])
+
+    VIFBvector <-setNames(VIFBDf$vif, VIFBDf$Relation)
+
+  } else{
+    VIFBvector <- NULL
+  }
+  return(VIFBvector)
+
+}
diff --git a/tests/testthat/test-plssem.R b/tests/testthat/test-plssem.R
index 2e49d00b..4b8e46b4 100644
--- a/tests/testthat/test-plssem.R
+++ b/tests/testthat/test-plssem.R
@@ -1,20 +1,20 @@
 context("Partial Least Squares Structural Equation Modeling")
 
-test_that("Basic PLSSEM works", {
 
-  options <- jaspTools::analysisOptions("PLSSEM")
-  model <- "
-    ind60 =~ x1 + x2 + x3
-    dem60 =~ y1 + y2 + y3 + y4
-    dem65 =~ y5 + y6 + y7 + y8
-    dem60 ~ ind60
-    dem65 ~ ind60 + dem60
-  "
-  options$models <- list(list(name = "Model1", syntax = list(model = model, columns = c("x1", "x2", "x3", "y1", "y2", "y3", "y4", "y5", "y6", "y7", "y8"))))
-  options$group                     <- ""
-  options$innerWeightingScheme      <- "path"
-  options$convergenceCriterion      <- "absoluteDifference"
-  options$correctionFactor          <- "squaredEuclidean"
+
+options <- jaspTools::analysisOptions("PLSSEM")
+model <- "
+  ind60 =~ x1 + x2 + x3
+  dem60 =~ y1 + y2 + y3 + y4
+  dem65 =~ y5 + y6 + y7 + y8
+  dem60 ~ ind60
+  dem65 ~ ind60 + dem60
+"
+options$models <- list(list(name = "Model1", syntax = list(model = model, columns = c("x1", "x2", "x3", "y1", "y2", "y3", "y4", "y5", "y6", "y7", "y8"))))
+options$group                     <- ""
+options$innerWeightingScheme      <- "path"
+options$convergenceCriterion      <- "absoluteDifference"
+options$correctionFactor          <- "squaredEuclidean"
 
 
 results <- jaspTools::runAnalysis("PLSSEM", "poldem_grouped.csv", options)