diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml
index 2c5bb50..f190d51 100755
--- a/.github/workflows/test-coverage.yaml
+++ b/.github/workflows/test-coverage.yaml
@@ -44,7 +44,7 @@ jobs:
 
       - name: Upload test results
         if: failure()
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: coverage-test-failures
           path: ${{ runner.temp }}/package
diff --git a/tests/testthat/test-simulations.R b/tests/testthat/test-simulations.R
index 926ae56..026f191 100755
--- a/tests/testthat/test-simulations.R
+++ b/tests/testthat/test-simulations.R
@@ -1,9 +1,11 @@
 expected_bre_fit <- list(
-  est.table = structure(
-    c(-0.60955243747617, -0.443501946207444, 0.736943004541776,
-      0.0104576873063189, 0.000414042620671129, 0.0123799088684966),
-    .Dim = 3:2,
-    .Dimnames = list(
+  est.table = matrix(
+    c(
+      -0.60955243747617, -0.443501946207444, 0.736943004541776,
+      0.0104576873063189, 0.000414042620671129, 0.0123799088684966
+    ),
+    nrow = 3,
+    dimnames = list(
       c("intersept", "log.Time", "Treatment"),
       c("Estimate", "Pvalue")
     )
@@ -14,32 +16,26 @@ expected_bre_fit <- list(
 
 test_that("simulate_beta_data works", {
   beta_data <- simulate_beta_random_effect_data()
-
-  # betareg can not fit random effect model so set the s2 to a small value (small random effect)
-  tdata <- data.frame(Y = beta_data$Y, beta_data$Z, SID = beta_data$subject_ind)
-  gy <- betareg::betareg(beta_data$Y ~ log.Time + as.factor(Treatment), data = tdata, type = "ML")
-  #print(summary(gy))
-  gy <- betareg::betareg(beta_data$Y ~ log.Time + as.factor(Treatment), data = tdata, type = "BC")
-  #print(summary(gy))
-  gy <- betareg::betareg(beta_data$Y ~ log.Time + as.factor(Treatment), data = tdata, type = "BR")
-  #print(summary(gy))
-
-  res <- fit_beta_random_effect(Z = beta_data$Z,
-                                Y = beta_data$Y,
-                                subject.ind = beta_data$subject_ind,
-                                time.ind = beta_data$time_ind,
-                                quad.n = 30,
-                                verbose = FALSE)
+  res <- fit_beta_random_effect(
+    Z = beta_data$Z,
+    Y = beta_data$Y,
+    subject.ind = beta_data$subject_ind,
+    time.ind = beta_data$time_ind,
+    quad.n = 30,
+    verbose = FALSE
+  )
 
   expect_equal(res, expected_bre_fit, tolerance = 1e-3)
 })
 
 expected_lre_fit <- list(
-  est.table = structure(
-    c(0.306795205483351, 0.280454254105502,
-      -0.667159092633291, 0.431194361789462, 0.266160176079596, 0.123214729652536),
-    .Dim = 3:2,
-    .Dimnames = list(c("intersept", "log.Time", "Treatment"), c("Estimate", "Pvalue"))
+  est.table = matrix(
+    c(
+      0.306795205483351, 0.280454254105502,
+      -0.667159092633291, 0.431194361789462, 0.266160176079596, 0.123214729652536
+    ),
+    nrow = 3,
+    dimnames = list(c("intersept", "log.Time", "Treatment"), c("Estimate", "Pvalue"))
   ),
   s1.est = 1.12114408844058
 )
@@ -47,19 +43,23 @@ expected_lre_fit <- list(
 test_that("simulate_logistic_data works", {
   logistic_data <- simulate_logistic_data()
 
-  tdata <- data.frame(Y = logistic_data$Y, logistic_data$X, SID = logistic_data$subject_ind)
-  lme.fit <- lme4::glmer(as.factor(logistic_data$Y) ~ log.Time + as.factor(Treatment) + (1 | SID),
-                         data = tdata,
-                         family = binomial,
-                         control = lme4::glmerControl(optimizer = "bobyqa"),
-                         nAGQ = 10)
+  tdata <- data.frame(
+    Y = logistic_data$Y, logistic_data$X, SID = logistic_data$subject_ind
+  )
+  lme.fit <- lme4::glmer(
+    as.factor(logistic_data$Y) ~ log.Time + as.factor(Treatment) + (1 | SID),
+    data = tdata,
+    family = binomial,
+    control = lme4::glmerControl(optimizer = "bobyqa"),
+    nAGQ = 10
+  )
   summary(lme.fit)
 
-  res <- fit_logistic_random_effect(X = logistic_data$X,
-                                    Y = logistic_data$Y,
-                                    subject.ind = logistic_data$subject_ind,
-                                    time.ind = logistic_data$time_ind)
-  dput(res)
-
+  res <- fit_logistic_random_effect(
+    X = logistic_data$X,
+    Y = logistic_data$Y,
+    subject.ind = logistic_data$subject_ind,
+    time.ind = logistic_data$time_ind
+  )
   expect_equal(res, expected_lre_fit, tolerance = 1e-3)
 })
diff --git a/tests/testthat/test-zibr.R b/tests/testthat/test-zibr.R
index f2efa6d..dead8a3 100755
--- a/tests/testthat/test-zibr.R
+++ b/tests/testthat/test-zibr.R
@@ -1,17 +1,17 @@
 real_expected <- list(
-  logistic_est_table = structure(
+  logistic_est_table = matrix(
     c(2.69732340911363, 0.200018071642455, 6.17105397782147e-06, 0.887710710500198),
-    .Dim = c(2L, 2L),
-    .Dimnames = list(
+    nrow = 2,
+    dimnames = list(
       c("intersept", "Treatment"),
       c("Estimate", "Pvalue")
     )
   ),
   logistic_s1_est = 3.27876357265477,
-  beta_est_table = structure(
+  beta_est_table = matrix(
     c(-2.78605402671509, -0.32334269957996, 0, 0.206300170582503),
-    .Dim = c(2L, 2L),
-    .Dimnames = list(
+    nrow = 2,
+    dimnames = list(
       c("intersept", "Treatment"),
       c("Estimate", "Pvalue")
     )
@@ -30,13 +30,15 @@ test_that("zibr main function works on package data", {
     time_ind = ibd$Time
   )
 
-  expected_names <- c("logistic_est_table",
-                      "logistic_s1_est",
-                      "beta_est_table",
-                      "beta_s2_est",
-                      "beta_v_est",
-                      "loglikelihood",
-                      "joint_p")
+  expected_names <- c(
+    "logistic_est_table",
+    "logistic_s1_est",
+    "beta_est_table",
+    "beta_s2_est",
+    "beta_v_est",
+    "loglikelihood",
+    "joint_p"
+  )
   expect_equal(names(zibr_fit), expected_names)
   expect_equal(zibr_fit, real_expected, tolerance = 1e-3)
 })
@@ -55,19 +57,19 @@ sim <- simulate_zero_inflated_beta_random_effect_data(
 )
 
 sim_expected <- list(
-  logistic_est_table = structure(
+  logistic_est_table = matrix(
     c(-0.571900346260645, 0.828072571379781, 0.0069577928092327, 0.00547006539016526),
-    .Dim = c(2L, 2L),
-    .Dimnames = list(
+    nrow = 2,
+    dimnames = list(
       c("intersept", "var1"),
       c("Estimate", "Pvalue")
     )
   ),
   logistic_s1_est = 1.06801391077711,
-  beta_est_table = structure(
+  beta_est_table = matrix(
     c(-0.593090577296489, 0.591745582740125, 4.18162741963046e-05, 0.00169945318503251),
-    .Dim = c(2L, 2L),
-    .Dimnames = list(
+    nrow = 2,
+    dimnames = list(
       c("intersept", "var1"),
       c("Estimate", "Pvalue")
     )
@@ -79,19 +81,19 @@ sim_expected <- list(
 )
 
 sim_expected_cov <- list(
-  logistic_est_table = structure(
+  logistic_est_table = matrix(
     c(-0.571900346260645, 0.828072571379781, 0.0069577928092327, 0.00547006539016526),
-    .Dim = c(2L, 2L),
-    .Dimnames = list(
+    nrow = 2,
+    dimnames = list(
       c("intersept", "var1"),
       c("Estimate", "Pvalue")
     )
   ),
   logistic_s1_est = 1.06801391077711,
-  beta_est_table = structure(
+  beta_est_table = matrix(
     c(-0.593090577296489, 0.591745582740125, 4.18162741963046e-05, 0.00169945318503251),
-    .Dim = c(2L, 2L),
-    .Dimnames = list(
+    nrow = 2,
+    dimnames = list(
       c("intersept", "var1"),
       c("Estimate", "Pvalue")
     )
@@ -103,21 +105,25 @@ sim_expected_cov <- list(
 )
 
 test_that("zibr main function works on simulated data", {
-  zibr_fit <- zibr(logistic_cov = sim$X,
-                   beta_cov = sim$Z,
-                   Y = sim$Y,
-                   subject_ind = sim$subject_ind,
-                   time_ind = sim$time_ind)
+  zibr_fit <- zibr(
+    logistic_cov = sim$X,
+    beta_cov = sim$Z,
+    Y = sim$Y,
+    subject_ind = sim$subject_ind,
+    time_ind = sim$time_ind
+  )
 
   expect_equal(zibr_fit, sim_expected, tolerance = 1e-3)
 })
 
 test_that("zibr main function works on simulated data with same covariates", {
-  zibr_fit <- zibr(logistic_cov = sim$X,
-                   beta_cov = sim$X,
-                   Y = sim$Y,
-                   subject_ind = sim$subject_ind,
-                   time_ind = sim$time_ind)
+  zibr_fit <- zibr(
+    logistic_cov = sim$X,
+    beta_cov = sim$X,
+    Y = sim$Y,
+    subject_ind = sim$subject_ind,
+    time_ind = sim$time_ind
+  )
 
   expect_equal(zibr_fit, sim_expected_cov, tolerance = 1e-3)
 })