diff --git a/dev/dev.R b/dev/dev.R
index 7d62418..f42a87a 100644
--- a/dev/dev.R
+++ b/dev/dev.R
@@ -1,2 +1,8 @@
 # Development
 devtools::load_all("/Users/markolalovic/dev/latent2likert/")
+
+covr::codecov(
+  quiet = FALSE,
+  clean = FALSE,
+  token = Sys.getenv("CODECOV_TOKEN")
+)
diff --git a/tests/testthat/test_discretization.R b/tests/testthat/test_discretization.R
index 7b09a23..7511131 100644
--- a/tests/testthat/test_discretization.R
+++ b/tests/testthat/test_discretization.R
@@ -1,22 +1,21 @@
 context("Testing discretization")
 
 testthat::test_that("discretization of N(0,1) using 4 levels gives expected results", {
-    #print("TODO: test coverage")
-    res = discretize_density(density_fn = dnorm, n_levels = 4)
-    expected = list(
-        prob = c(0.16, 0.34, 0.34, 0.16), 
-        endp = c(-Inf, -0.98, 0, 0.98, Inf), 
-        repr = c(-1.51, -0.45, 0.45, 1.51), 
-        dist = 0.12)
-    testthat::expect_equal(res, expected, tolerance = 0.05)
+  res = discretize_density(density_fn = dnorm, n_levels = 4)
+  expected = list(
+      prob = c(0.16, 0.34, 0.34, 0.16), 
+      endp = c(-Inf, -0.98, 0, 0.98, Inf), 
+      repr = c(-1.51, -0.45, 0.45, 1.51), 
+      dist = 0.12)
+  testthat::expect_equal(res, expected, tolerance = 0.05)
 })
 
 testthat::test_that("discretization of N(0,1) using 5 levels gives expected results", {
-    res = discretize_density(density_fn = dnorm, n_levels = 5)
-    expected = list(
-        prob = c(0.11, 0.24, 0.30, .24, 0.11), 
-        endp = c(-Inf, -1.25, -0.38, 0.38, 1.25, Inf), 
-        repr = c(-1.73, -0.77, 0, 0.77, 1.73), 
-        dist = 0.08)
-    testthat::expect_equal(res, expected, tolerance = 0.05)
+  res = discretize_density(density_fn = dnorm, n_levels = 5)
+  expected = list(
+      prob = c(0.11, 0.24, 0.30, .24, 0.11), 
+      endp = c(-Inf, -1.25, -0.38, 0.38, 1.25, Inf), 
+      repr = c(-1.73, -0.77, 0, 0.77, 1.73), 
+      dist = 0.08)
+  testthat::expect_equal(res, expected, tolerance = 0.05)
 })
diff --git a/tests/testthat/test_estimation.R b/tests/testthat/test_estimation.R
index dfe97c0..fd92e8a 100644
--- a/tests/testthat/test_estimation.R
+++ b/tests/testthat/test_estimation.R
@@ -1,13 +1,15 @@
 context("Testing estimation of parameters")
 
-testthat::test_that("`mean` and `sd` are well estimated given `pk` for normal case", {
+testthat::test_that("`mean` and `sd` are well estimated given `prob`
+  for normal case", {
   actual <- c(-1, 0.5) # actual mean and sd
   prob <- c("1" = 0.313, "2" = 0.579, "3" = 0.105, "4" = 0.003)
   estimates <- estimate_mean_and_sd(prob, 5)
   testthat::expect_equal(estimates, actual, tolerance = 0.1)
 })
 
-testthat::test_that("`mean` and `sd` are well estimated given `pk` for skew case", {
+testthat::test_that("`mean` and `sd` are well estimated given `prob`
+  for skew case", {
   actual <- c(0.146858, 1.084341) # actual mean and sd
   skew <- -0.4565873 # skewness
   prob <- c(
@@ -17,3 +19,85 @@ testthat::test_that("`mean` and `sd` are well estimated given `pk` for skew case
   estimates <- estimate_mean_and_sd(prob, 7, skew)
   testthat::expect_equal(estimates, actual, tolerance = 0.05)
 })
+
+testthat::test_that("estimate_params returns accurate estimates for 
+  latent parameters when data is a vector", {
+  set.seed(12345) # for reproducibility
+  
+  # Generate test data for a single item (vector data)
+  data_vector <- rlikert(size = 1000,
+                         n_items = 1,
+                         n_levels = 5,
+                         mean = 0.5, 
+                         sd   = 1.2,
+                         skew = 0)
+  
+  # Estimate parameters
+  estimates <- estimate_params(data_vector, n_levels = 5, skew = 0)
+  
+  # Extract estimated means and standard deviations
+  estimated_mean <- as.numeric(estimates["mean"])
+  estimated_sd <- as.numeric(estimates["sd"])
+  
+  # Actual means and standard deviations
+  actual_mean <- 0.5
+  actual_sd <- 1.2
+  
+  # Check if estimated means are close to actual means
+  testthat::expect_equal(estimated_mean, actual_mean, tolerance = 0.1)
+  
+  # Check if estimated standard deviations are close to actual standard deviations
+  testthat::expect_equal(estimated_sd, actual_sd, tolerance = 0.1)
+})
+
+testthat::test_that("estimate_params returns accurate estimates for 
+latent parameters", {
+  set.seed(12345) # for reproducibility
+  
+  # Generate test data
+  corr <- matrix(c(1.00, -0.63, -0.39, 
+                   -0.63, 1.00, 0.41, 
+                   -0.39, 0.41, 1.00), nrow = 3)
+  data <- rlikert(size = 1000,
+                  n_items = 3,
+                  n_levels = c(4, 5, 6),
+                  mean = c(0, -1, 0), 
+                  sd   = c(0.8, 1, 1),
+                  corr = corr)
+  
+  # Estimate parameters
+  estimates <- estimate_params(data, n_levels = c(4, 5, 6), skew = 0)
+  
+  # Extract estimated means and standard deviations
+  estimated_means <- as.numeric(estimates["mean", ])
+  estimated_sds <- as.numeric(estimates["sd", ])
+  
+  # Actual means and standard deviations
+  actual_means <- c(0, -1, 0)
+  actual_sds <- c(0.8, 1, 1)
+  
+  # Check if estimated means are close to actual means
+  testthat::expect_equal(estimated_means, actual_means, tolerance = 0.1)
+  
+  # Check if estimated sds are close to actual sds
+  testthat::expect_equal(estimated_sds, actual_sds, tolerance = 0.1)
+})
+
+testthat::test_that("plot_contour executes without errors", {
+  # Define a simple objective function for testing
+  test_fn <- function(x, endp, prob, cdf_X) {
+    u <- x[1]
+    v <- x[2]
+    y <- cdf_X(v * endp - u * v)
+    return(matrix(utils::tail(y, -1) - utils::head(y, -1) - prob))
+  }
+  
+  # Example parameters for the test
+  endp <- c(-Inf, -1, 0, 1, Inf)
+  prob <- c(0.1, 0.2, 0.4, 0.3)
+  cdf_X <- stats::pnorm
+  trace <- matrix(c(rep(0, 100), rep(1, 100)), nrow = 2)
+  
+  # Run the plotting function
+  testthat::expect_silent(plot_contour(test_fn, endp, prob, cdf_X, trace))
+})
diff --git a/tests/testthat/test_simulation.R b/tests/testthat/test_simulation.R
index 03c16ce..95ee455 100644
--- a/tests/testthat/test_simulation.R
+++ b/tests/testthat/test_simulation.R
@@ -1,103 +1,34 @@
 context("Testing simulation")
 
-testthat::test_that("proportions of generated responses match actual probabilities", {
-    size <- 10^6
-    n_levels <- 5
-    n_items <- 3
-
-    mean <- c(-1, 0, 1)
-    sd <- c(1, 1, 0.5)
-    skew <- c(0.5, 0.5, 0.5)
-    corr <- 0.5
-    
-    data <- rlikert(size, n_items, n_levels, mean, sd, skew, corr)
-    for (i in seq_len(n_items)) {
-        cp <- c("mu"=mean[i], "sd"=sd[i], "skew"=skew[i])
-        prob <- simulate_likert(n_levels, cp)
-        data_prop <- response_prop(data[,i], n_levels)
-        testthat::expect_equal(prob, data_prop, tolerance = 0.05)
-    }
+testthat::test_that("simulate_likert using 4 levels and N(0,1) gives 
+  expected result", {
+  n_levels <- 4
+  cp <- c("mu"=0, "sd"=1, "skew"=0)
+  prob <- simulate_likert(n_levels, cp)
+  prob <- as.numeric(prob)
+  expected_prob <- c(0.163, 0.337, 0.337, 0.163)
+  testthat::expect_equal(prob, expected_prob, tolerance = 0.05)
 })
 
-testthat::test_that("proportions of generated responses match actual probabilities, 
-                     when using a random corr matrix", {
-    size <- 10^6
-    n_levels <- 5
-    n_items <- 3
-
-    mean <- c(0, -1, -1)
-    sd <- c(1, 1, 0.5)
-    skew <- c(0, 0, 0)
-    corr <- "random"
-    
-    data <- rlikert(size, n_items, n_levels, mean, sd, skew, corr)
-    for (i in seq_len(n_items)) {
-        cp <- c("mu"=mean[i], "sd"=sd[i], "skew"=skew[i])
-        prob <- simulate_likert(n_levels, cp)
-        data_prop <- response_prop(data[,i], n_levels)
-        testthat::expect_equal(prob, data_prop, tolerance = 0.05)
-    }
-})
-
-testthat::test_that("correlations of random responses match actual correlations", {
-    size <- 10^6
-    n_levels <- 5
-    n_items <- 3
-    
-    mean <- c(0, -1, -1)
-    sd <- c(1, 1, 0.5)
-    skew <- c(0.5, 0.5, 0.5)
-    corr <- 0.5
-    
-    data <- rlikert(size, n_items, n_levels, mean, sd, skew, corr)
-    data_corr_matrix <- cor(data)
-    
-    actual_corr_matrix <- generate_corr_matrix(corr, 3, n_items)
-    dimnames(actual_corr_matrix) <- dimnames(data_corr_matrix)
-
-    testthat::expect_equal(actual_corr_matrix, 
-                           data_corr_matrix, tolerance = 0.05)
-})
-
-testthat::test_that("correlations of random responses match actual correlations,
-                     harder case", {
-    size <- 10^6
-    n_levels <- 6
-    n_items <- 3
-    
-    mean <- c(-0.5, 0, 0.5)
-    sd <- c(0.5, 0.5, 0.5)
-    skew <- c(-0.3, -0.4, -0.5)
-    corr <- 0.7
-    
-    data <- rlikert(size, n_items, n_levels, mean, sd, skew, corr)
-    data_corr_matrix <- cor(data)
-    
-    actual_corr_matrix <- generate_corr_matrix(corr, 3, n_items)
-    dimnames(actual_corr_matrix) <- dimnames(data_corr_matrix)
-
-    testthat::expect_equal(actual_corr_matrix, 
-                           data_corr_matrix, tolerance = 0.11)
-})
-
-testthat::test_that("simulate_likert using 5 levels and N(0,1) gives expected result", {
-    n_levels <- 5
-    cp <- c("mu"=0, "sd"=1, "skew"=0)
-    prob <- simulate_likert(n_levels, cp)
-    prob <- as.numeric(prob)
-    expected_prob <- c(0.106, 0.244, 0.298, 0.244, 0.106)
-    testthat::expect_equal(prob, expected_prob, tolerance = 0.05)
+testthat::test_that("simulate_likert using 5 levels and N(0,1) gives 
+  expected result", {
+  n_levels <- 5
+  cp <- c("mu"=0, "sd"=1, "skew"=0)
+  prob <- simulate_likert(n_levels, cp)
+  prob <- as.numeric(prob)
+  expected_prob <- c(0.106, 0.244, 0.298, 0.244, 0.106)
+  testthat::expect_equal(prob, expected_prob, tolerance = 0.05)
 })
 
 testthat::test_that("Correlation input is handled correctly", {
-    corr_matrix <- generate_rand_corr_matrix(3)
-    corr_inputs <- list(0, "random", 0.5, corr_matrix)
-    
-    for (i in seq_along(corr_inputs)) {
-        corr <- corr_inputs[[i]]
-        corr_case <- handle_corr_case(corr)
-        testthat::expect_equal(i, corr_case)
-    }
+  corr_matrix <- generate_rand_corr_matrix(3)
+  corr_inputs <- list(0, "random", 0.5, corr_matrix)
+  
+  for (i in seq_along(corr_inputs)) {
+      corr <- corr_inputs[[i]]
+      corr_case <- handle_corr_case(corr)
+      testthat::expect_equal(i, corr_case)
+  }
 })
 
 testthat::test_that("Invalid corr input raises error", {
@@ -105,7 +36,8 @@ testthat::test_that("Invalid corr input raises error", {
   testthat::expect_equal(class(res), "try-error")
 })
 
-testthat::test_that("generate_corr_matrix returns a matrix that resembles a correlation matrix", {
+testthat::test_that("generate_corr_matrix returns a matrix that 
+  resembles a correlation matrix", {
   n_items <- 3
   corr_matrix <- generate_rand_corr_matrix(n_items)
   corr_inputs <- list("random", 0.5, corr_matrix)
@@ -115,7 +47,134 @@ testthat::test_that("generate_corr_matrix returns a matrix that resembles a corr
     corr_case <- handle_corr_case(corr)
     res <- generate_corr_matrix(corr, corr_case, n_items)
 
-    testthat::expect_true(all(diag(res) == 1), info = paste("Failed on input:", i))
-    testthat::expect_true(isSymmetric(res), info = paste("Failed on input:", i))
+    testthat::expect_true(all(diag(res) == 1), 
+                          info = paste("Failed on input:", i))
+    testthat::expect_true(isSymmetric(res), 
+                          info = paste("Failed on input:", i))
+  }
+})
+
+testthat::test_that("proportions of generated responses match actual 
+  probabilities for single item", {
+  set.seed(12345) # for reproducibility
+  size <- 10^6
+  n_items <- 1
+  n_levels <- 5
+  mean <- 0
+  sd <- 1
+  skew <- 0
+  
+  data <- rlikert(size, n_items, n_levels, mean, sd, skew)
+
+  cp <- c("mu"=mean, "sd"=sd, "skew"=skew)
+  prob <- simulate_likert(n_levels, cp)
+  data_prop <- response_prop(data, n_levels)
+  testthat::expect_equal(prob, data_prop, tolerance = 0.05)
+})
+
+testthat::test_that("proportions of generated responses match actual 
+  probabilities for multiple items without corr", {
+  set.seed(12345) # for reproducibility
+  size <- 10^6
+  n_levels <- 5
+  n_items <- 3
+
+  mean <- c(-1, 0, 1)
+  sd <- c(0.5, 1, 0.5)
+  skew <- c(0.5, 0.5, 0.5)
+  corr <- 0
+  
+  data <- rlikert(size, n_items, n_levels, mean, sd, skew, corr)
+  for (i in seq_len(n_items)) {
+      cp <- c("mu"=mean[i], "sd"=sd[i], "skew"=skew[i])
+      prob <- simulate_likert(n_levels, cp)
+      data_prop <- response_prop(data[,i], n_levels)
+      testthat::expect_equal(prob, data_prop, tolerance = 0.05)
+  }
+})
+
+testthat::test_that("proportions of generated responses match actual 
+  probabilities for multiple items with corr", {
+  set.seed(12345) # for reproducibility
+  size <- 10^6
+  n_levels <- 5
+  n_items <- 3
+
+  mean <- c(-1, 0, 1)
+  sd <- c(1, 1, 0.5)
+  skew <- c(0.5, 0.5, 0.5)
+  corr <- 0.5
+  
+  data <- rlikert(size, n_items, n_levels, mean, sd, skew, corr)
+  for (i in seq_len(n_items)) {
+      cp <- c("mu"=mean[i], "sd"=sd[i], "skew"=skew[i])
+      prob <- simulate_likert(n_levels, cp)
+      data_prop <- response_prop(data[,i], n_levels)
+      testthat::expect_equal(prob, data_prop, tolerance = 0.05)
+  }
+})
+
+testthat::test_that("proportions of generated responses match actual 
+  probabilities, when using a random corr matrix", {
+  set.seed(12345) # for reproducibility    
+  size <- 10^6
+  n_levels <- 5
+  n_items <- 3
+
+  mean <- c(0, -1, -1)
+  sd <- c(1, 1, 0.5)
+  skew <- c(0, 0, 0)
+  corr <- "random"
+  
+  data <- rlikert(size, n_items, n_levels, mean, sd, skew, corr)
+  for (i in seq_len(n_items)) {
+      cp <- c("mu"=mean[i], "sd"=sd[i], "skew"=skew[i])
+      prob <- simulate_likert(n_levels, cp)
+      data_prop <- response_prop(data[,i], n_levels)
+      testthat::expect_equal(prob, data_prop, tolerance = 0.05)
   }
 })
+
+testthat::test_that("correlations of random responses match actual 
+  correlations", {
+  set.seed(12345) # for reproducibility    
+  size <- 10^6
+  n_levels <- 5
+  n_items <- 3
+  
+  mean <- c(0, -1, -1)
+  sd <- c(1, 1, 0.5)
+  skew <- c(0.5, 0.5, 0.5)
+  corr <- 0.5
+  
+  data <- rlikert(size, n_items, n_levels, mean, sd, skew, corr)
+  data_corr_matrix <- cor(data)
+  
+  actual_corr_matrix <- generate_corr_matrix(corr, 3, n_items)
+  dimnames(actual_corr_matrix) <- dimnames(data_corr_matrix)
+
+  testthat::expect_equal(actual_corr_matrix, 
+                        data_corr_matrix, tolerance = 0.05)
+})
+
+testthat::test_that("correlations of random responses match actual 
+  correlations, harder case", {
+  set.seed(12345) # for reproducibility
+  size <- 10^6
+  n_levels <- 6
+  n_items <- 3
+  
+  mean <- c(-0.5, 0, 0.5)
+  sd <- c(0.5, 0.5, 0.5)
+  skew <- c(-0.3, -0.4, -0.5)
+  corr <- 0.7
+  
+  data <- rlikert(size, n_items, n_levels, mean, sd, skew, corr)
+  data_corr_matrix <- cor(data)
+  
+  actual_corr_matrix <- generate_corr_matrix(corr, 3, n_items)
+  dimnames(actual_corr_matrix) <- dimnames(data_corr_matrix)
+
+  testthat::expect_equal(actual_corr_matrix, 
+                          data_corr_matrix, tolerance = 0.11)
+})
diff --git a/tests/testthat/test_utilities.R b/tests/testthat/test_utilities.R
index 4f84b9b..b6c209a 100644
--- a/tests/testthat/test_utilities.R
+++ b/tests/testthat/test_utilities.R
@@ -8,7 +8,8 @@ testthat::test_that("pad_levels gives the correct result", {
   testthat::expect_equal(padded_pr, actual_pr)
 })
 
-testthat::test_that("response_prop gives the correct result, univariate case", {
+testthat::test_that("response_prop gives the correct result, 
+  univariate case", {
   data <- rep(c(1, 2, 3, 4), each = 2)
   tab <- response_prop(data, n_levels = 4)
   correct_tab <- rep(0.25, 4)
@@ -16,7 +17,8 @@ testthat::test_that("response_prop gives the correct result, univariate case", {
   testthat::expect_true(identical(tab, correct_tab))
 })
 
-testthat::test_that("response_prop gives the correct result, multivariate case", {
+testthat::test_that("response_prop gives the correct result, 
+  multivariate case", {
   y <- rep(c(1, 2, 3, 4), each = 2)
   data <- cbind(y, y)
   tab <- response_prop(data, n_levels = 4)
@@ -63,3 +65,55 @@ testthat::test_that("density_sn gives the same results as sn::dsn", {
 
   testthat::expect_equal(y1, y2)
 })
+
+
+testthat::test_that("plot_likert_transform runs without errors", {
+  testthat::expect_error(
+    plot_likert_transform(n_items = 3, n_levels = c(3, 4, 5)), 
+    NA
+  )
+  testthat::expect_error(
+    plot_likert_transform(n_items = 3, n_levels = 5, mean = c(0, 1, 2)),
+    NA
+  )
+  testthat::expect_error(
+    plot_likert_transform(n_items = 3, n_levels = 5, sd = c(0.8, 1, 1.2)),
+    NA
+  )
+  testthat::expect_error(
+    plot_likert_transform(n_items = 3, n_levels = 5, skew = c(-0.5, 0, 0.5)),
+    NA
+  )
+})
+
+testthat::test_that("delta_skew_normal returns correct value", {
+  alpha <- 1
+  result <- delta_skew_normal(alpha)
+  expected_result <- 0.71
+  testthat::expect_equal(result, expected_result, tolerance = 0.05)
+})
+
+testthat::test_that("mean_skew_normal returns correct value", {
+  alpha <- 1
+  result <- mean_skew_normal(alpha)
+  expected_result <- 0.56
+  testthat::expect_equal(result, expected_result, tolerance = 0.05)
+})
+
+testthat::test_that("var_skew_normal returns correct value", {
+  alpha <- 1
+  result <- var_skew_normal(alpha)
+  expected_result <- 0.68
+  testthat::expect_equal(result, expected_result, tolerance = 0.05)
+})
+
+testthat::test_that("scale_and_shift returns correct values", {
+  dp <- c(xi = 1, omega = 2, alpha = 1)
+  x <- c(-1, 0, 1)
+  result <- scale_and_shift(x, dp)
+  expected_result <- c(-0.7179052, -0.2179052, 0.2820948)
+  testthat::expect_equal(result, expected_result, tolerance = 0.05)
+})
+
+
+