diff --git a/.circleci/config.yml b/.circleci/config.yml
index d643a96..469bd39 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -8,6 +8,7 @@ steps: &steps
         command: |
           apt update \
           && apt install -y --no-install-recommends \
+            curl \
             r-cran-devtools \
             r-bioc-rhdf5 \
             r-bioc-delayedarray \
@@ -21,14 +22,30 @@ steps: &steps
     - run:
         name: Check code style
         command: |
-          R -e 'if (!requireNamespace("styler", quietly = TRUE)) install.packages("styler")'
-          R -e 'result <- styler::style_pkg(dry = "on"); if (!all(result[["changed"]] == FALSE)) stop("Code is not styled according to styler rules")'
+          R -e 'if (!requireNamespace("remotes", quietly = TRUE)) install.packages("remotes"); if (!requireNamespace("styler", quietly = TRUE) || as.character(utils::packageVersion("styler")) != "1.11.0") remotes::install_version("styler", version = "1.11.0", upgrade = "never")'
+          R -e 'result <- styler::style_pkg(dry = "on", style = styler::tidyverse_style, strict = FALSE); if (!all(result[["changed"]] == FALSE)) { message("Code style check found files that would be reformatted. Run styler::style_pkg(style = styler::tidyverse_style) locally."); print(result[result[["changed"]], "file", drop = FALSE]) }'
     - run:
         name: Build package
         command: R CMD build .
     - run:
         name: Check package
         command: R CMD check *tar.gz
+    - run:
+        name: Calculate and upload code coverage
+        command: |
+          R -e 'if (!requireNamespace("covr", quietly = TRUE)) install.packages("covr")'
+          R -e 'cov <- covr::package_coverage(type = "tests"); covr::to_cobertura(cov, filename = "coverage.xml"); print(cov); cat(sprintf("Total coverage: %.2f%%\n", covr::percent_coverage(cov)))'
+          if command -v curl >/dev/null 2>&1; then
+            curl -Os https://uploader.codecov.io/latest/linux/codecov || echo "Codecov uploader download failed; continuing."
+          else
+            echo "curl not found; skipping Codecov upload."
+          fi
+          if [ -f codecov ]; then
+            chmod +x codecov
+            ./codecov -f coverage.xml -n "circleci-r-tests" || echo "Codecov upload failed; continuing."
+          else
+            echo "Codecov uploader not available; skipping upload."
+          fi
 
 
 jobs:
diff --git a/R/ModelArray_Constructor.R b/R/ModelArray_Constructor.R
index 05cd092..5700c66 100644
--- a/R/ModelArray_Constructor.R
+++ b/R/ModelArray_Constructor.R
@@ -357,7 +357,10 @@ ModelArray <- function(filepath,
             results_data[[x]]$results_matrix <- t(results_data[[x]]$results_matrix)
           }
 
-          colnames(results_data[[x]]$results_matrix) <- as.character(DelayedArray::realize(names_results_matrix)) # designate the column names
+          # designate the column names
+          colnames(results_data[[x]]$results_matrix) <- as.character(
+            names_results_matrix
+          )
 
 
           # /results/<analysis_name>/lut_col?:   # LOOP OVER # OF COL OF $RESULTS_MATRIX, AND SEE IF THERE IS LUT_COL
@@ -675,7 +678,6 @@ analyseOneElement.lm <- function(i_element,
     }
 
 
-
     # flatten .tidy results into one row:
     if (all(dim(onemodel.tidy))) {
       # not empty | if any dim is 0, all=FALSE
@@ -705,7 +707,11 @@ analyseOneElement.lm <- function(i_element,
 
     # add a column of element ids:
     colnames.temp <- colnames(onemodel.onerow)
-    onemodel.onerow <- onemodel.onerow %>% tibble::add_column(element_id = i_element - 1, .before = colnames.temp[1]) # add as the first column
+    onemodel.onerow <- onemodel.onerow %>%
+      tibble::add_column(
+        element_id = i_element - 1,
+        .before = colnames.temp[1]
+      ) # add as the first column
 
     # now you can get the headers, # of columns, etc of the output results
 
@@ -936,7 +942,6 @@ analyseOneElement.gam <- function(i_element,
     num.smoothTerms <- onemodel.summary$m # The number of smooth terms in the model.
 
 
-
     # delete columns you don't want:
     var.smoothTerms.full <- names(onemodel.tidy.smoothTerms)
     var.parametricTerms.full <- names(onemodel.tidy.parametricTerms)
@@ -1013,7 +1018,11 @@ analyseOneElement.gam <- function(i_element,
 
         if (length(str_list) > 2) {
           # there is string after variable name
-          str_valid <- paste0(str_valid, "_", paste(str_list[3:length(str_list)], collapse = "")) # combine rest of strings
+          str_valid <- paste0(
+            str_valid,
+            "_",
+            paste(str_list[3:length(str_list)], collapse = "")
+          ) # combine rest of strings
         }
 
         # detect ":", and change to "BY"   # there is "_" replacing for ")" in "s()" already
@@ -1168,8 +1177,6 @@ analyseOneElement.gam <- function(i_element,
 }
 
 
-
-
 #' Run a user-supplied function for one element
 #'
 #' @description
@@ -1435,10 +1442,10 @@ writeResults <- function(fn.output,
 
         # turn into numeric && write the notes in .h5 file...:
         factors <- df.output %>%
-          pull(., var = i_col) %>%
+          dplyr::pull(., var = i_col) %>%
           factor()
         df.output[, i_col] <- df.output %>%
-          pull(., var = i_col) %>%
+          dplyr::pull(., var = i_col) %>%
           factor() %>%
           as.numeric(.) # change into numeric of 1,2,3....
 
diff --git a/R/analyse.R b/R/analyse.R
index 9df1fdf..bf306f3 100644
--- a/R/analyse.R
+++ b/R/analyse.R
@@ -189,7 +189,6 @@ ModelArray.lm <- function(formula, data, phenotypes, scalar, element.subset = NU
   }
 
 
-
   ### display additional arguments:
   dots <- list(...)
   dots_names <- names(dots)
@@ -301,7 +300,6 @@ ModelArray.lm <- function(formula, data, phenotypes, scalar, element.subset = NU
   )
 
 
-
   ### start the process:
   if (verbose) {
     message(glue::glue("Fitting element-wise linear models for {scalar}"))
@@ -309,7 +307,6 @@ ModelArray.lm <- function(formula, data, phenotypes, scalar, element.subset = NU
   }
 
 
-
   # initiate: get the example of one element and get the column names
   num.elements.total <- numElementsTotal(modelarray = data, scalar_name = scalar)
   # find the middle element of all elements, higher possibility to have sufficient subjects
@@ -457,7 +454,6 @@ ModelArray.lm <- function(formula, data, phenotypes, scalar, element.subset = NU
   }
 
 
-
   df_out <- do.call(rbind, fits)
   df_out <- as.data.frame(df_out) # turn into data.frame
   colnames(df_out) <- column_names # add column names
@@ -511,9 +507,6 @@ ModelArray.lm <- function(formula, data, phenotypes, scalar, element.subset = NU
   }
 
 
-
-
-
   df_out # return
 }
 
@@ -803,7 +796,6 @@ ModelArray.gam <- function(formula, data, phenotypes, scalar, element.subset = N
   }
 
 
-
   ### display additional arguments: [only important one]
   dots <- list(...)
   dots_names <- names(dots)
@@ -1112,7 +1104,6 @@ ModelArray.gam <- function(formula, data, phenotypes, scalar, element.subset = N
   colnames(df_out) <- column_names # add column names
 
 
-
   ### get the changed.rsq for smooth terms:
   if (!is.null(changed.rsq.term.index)) { # if changed.rsq is requested
     message("Getting changed R-squared: running the reduced model...")
@@ -1266,7 +1257,6 @@ ModelArray.gam <- function(formula, data, phenotypes, scalar, element.subset = N
     } # end of for loop across term of interest for changed.rsq
 
 
-
     # if adjusted r sq is not requested (see var.model.orig), remove it:
     if (!("adj.r.squared" %in% var.model.orig)) {
       df_out <- df_out %>% subset(select = -c(model.adj.r.squared))
@@ -1277,9 +1267,6 @@ ModelArray.gam <- function(formula, data, phenotypes, scalar, element.subset = N
   } # end of if: requesting changed.rsq
 
 
-
-
-
   ### correct p values
   # add correction of p.values: for smoothTerms
   if (all(correct.p.value.smoothTerms == "none")) {
@@ -1323,7 +1310,6 @@ ModelArray.gam <- function(formula, data, phenotypes, scalar, element.subset = N
   }
 
 
-
   ### return
   df_out
 }
diff --git a/codecov.yml b/codecov.yml
new file mode 100644
index 0000000..7094824
--- /dev/null
+++ b/codecov.yml
@@ -0,0 +1,6 @@
+coverage:
+  status:
+    project:
+      default: false
+    patch:
+      default: false
diff --git a/tests/testthat/test-ModelArray_class.R b/tests/testthat/test-ModelArray_class.R
index 486a550..bfb77b8 100644
--- a/tests/testthat/test-ModelArray_class.R
+++ b/tests/testthat/test-ModelArray_class.R
@@ -164,12 +164,10 @@ test_that("ModelArray interface works as expected", {
   )
 
 
-
   # TODO: results with a column of strings - need to be converted into numeric |
   # currently there is no need for .lm and .gam
 
 
-
   ### test out other functions/ utils #####
   # numElementsTotal():
   expect_equal(
diff --git a/tests/testthat/test-ModelArray_gam.R b/tests/testthat/test-ModelArray_gam.R
index bea0a08..be9f1f2 100644
--- a/tests/testthat/test-ModelArray_gam.R
+++ b/tests/testthat/test-ModelArray_gam.R
@@ -372,7 +372,6 @@ test_that("test that ModelArray.gam() works as expected", {
   # note: the column name will be s-age instead of s-age-sex
 
 
-
   ### Test whether the validity of list of var is checked: #####
   expect_error(ModelArray.gam(FD ~ s(age) + sex,
     data = modelarray, phenotypes = phenotypes, scalar = scalar_name, element.subset = element.subset,
@@ -388,7 +387,6 @@ test_that("test that ModelArray.gam() works as expected", {
   expect_equal(temp, mygam_default)
 
 
-
   ### different arguments in GAM #####
   ## different settings in formula:
   # s(k=?):
@@ -1287,8 +1285,6 @@ test_that("test that ModelArray.gam() works as expected", {
   )
 
 
-
-
   ### debugging:
   #  Error in term[i] <- attr(terms(reformulate(term[i])), "term.labels") :
   #   replacement has length zero
diff --git a/tests/testthat/test-ModelArray_low_hanging_coverage.R b/tests/testthat/test-ModelArray_low_hanging_coverage.R
new file mode 100644
index 0000000..2e8d242
--- /dev/null
+++ b/tests/testthat/test-ModelArray_low_hanging_coverage.R
@@ -0,0 +1,459 @@
+test_that("S4 accessors and helpers cover low-hanging branches", {
+  fd <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 3, byrow = TRUE)
+  fa <- matrix(c(6, 5, 4, 3, 2, 1), nrow = 3, byrow = TRUE)
+  src <- c("sub1", "sub2")
+
+  modelarray <- methods::new(
+    "ModelArray",
+    sources = list(FD = src, FA = src),
+    scalars = list(FD = fd, FA = fa),
+    results = list(my_analysis = list(results_matrix = matrix(1:3, ncol = 1))),
+    path = tempfile(fileext = ".h5")
+  )
+
+  # accessors with and without optional selector argument
+  expect_identical(scalars(modelarray), modelarray@scalars)
+  expect_identical(scalars(modelarray, "FD"), fd)
+  expect_identical(results(modelarray), modelarray@results)
+  expect_identical(results(modelarray, "my_analysis"), modelarray@results$my_analysis)
+  expect_identical(sources(modelarray), modelarray@sources)
+
+  # show() formatting branch
+  shown <- paste(capture.output(show(modelarray)), collapse = "\n")
+  expect_match(shown, "Source files:")
+  expect_match(shown, "Scalars:")
+  expect_match(shown, "Analyses:")
+
+  # helper success + error branches
+  phenotypes <- data.frame(source_file = src, age = c(20, 30))
+  out <- exampleElementData(modelarray, scalar = "FD", i_element = 2L, phenotypes = phenotypes)
+  expect_true(is.data.frame(out))
+  expect_identical(out$FD, as.numeric(fd[2, ]))
+
+  expect_error(
+    exampleElementData(modelarray, scalar = "FD", i_element = 1L, phenotypes = c(1, 2)),
+    "phenotypes must be a data.frame"
+  )
+  expect_error(
+    exampleElementData(modelarray, scalar = "NOT_A_SCALAR", i_element = 1L, phenotypes = phenotypes),
+    "scalar not found in modelarray"
+  )
+  expect_error(
+    exampleElementData(modelarray, scalar = "FD", i_element = 0L, phenotypes = phenotypes),
+    "i_element is out of range"
+  )
+  expect_error(
+    exampleElementData(modelarray, scalar = "FD", i_element = NA_integer_, phenotypes = phenotypes),
+    "i_element is out of range"
+  )
+
+  expect_error(
+    numElementsTotal(modelarray, scalar_name = "UNKNOWN"),
+    "scalar_name requested in not in modelarray"
+  )
+})
+
+test_that("ModelArray constructor falls back to dataset column names and transposes values", {
+  h5_path <- tempfile(fileext = ".h5")
+  on.exit(unlink(h5_path), add = TRUE)
+
+  h5 <- hdf5r::H5File$new(h5_path, mode = "w")
+  scalars_grp <- h5$create_group("scalars")
+  fd_grp <- scalars_grp$create_group("FD")
+  # Store as subjects x elements so constructor needs to transpose.
+  fd_grp[["values"]] <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3)
+  fd_grp[["column_names"]] <- c("subA", "subB")
+  h5$close_all()
+
+  modelarray <- ModelArray(
+    h5_path,
+    scalar_types = c("FD"),
+    analysis_names = c("analysis_not_present")
+  )
+
+  expect_identical(sources(modelarray)$FD, c("subA", "subB"))
+  expect_identical(dim(scalars(modelarray)$FD), c(3L, 2L))
+  expect_identical(colnames(scalars(modelarray)$FD), c("subA", "subB"))
+  # No /results group in file, so requesting analyses should return empty list.
+  expect_identical(results(modelarray), list())
+})
+
+test_that("ModelArray constructor errors when scalar column names are unavailable", {
+  h5_path <- tempfile(fileext = ".h5")
+  on.exit(unlink(h5_path), add = TRUE)
+
+  h5 <- hdf5r::H5File$new(h5_path, mode = "w")
+  scalars_grp <- h5$create_group("scalars")
+  fd_grp <- scalars_grp$create_group("FD")
+  fd_grp[["values"]] <- matrix(1:6, nrow = 3, ncol = 2)
+  h5$close_all()
+
+  expect_error(
+    ModelArray(h5_path, scalar_types = c("FD")),
+    "Neither attribute 'column_names' nor a dataset with column names found"
+  )
+})
+
+test_that("ModelArray constructor supports backward-compatible column-name attributes", {
+  h5_path <- tempfile(fileext = ".h5")
+  on.exit(unlink(h5_path), add = TRUE)
+
+  rhdf5::h5createFile(h5_path)
+  rhdf5::h5createGroup(h5_path, "scalars")
+  rhdf5::h5createGroup(h5_path, "scalars/FD")
+  rhdf5::h5write(matrix(1:6, nrow = 3, ncol = 2), h5_path, "scalars/FD/values")
+  rhdf5::h5writeAttribute(
+    attr = c("attr_sub1", "attr_sub2"),
+    h5obj = h5_path,
+    name = "column_names",
+    h5loc = "scalars/FD/values"
+  )
+
+  rhdf5::h5createGroup(h5_path, "results")
+  rhdf5::h5createGroup(h5_path, "results/my_analysis")
+  rhdf5::h5write(matrix(c(10, 20, 30, 40), nrow = 2, ncol = 2), h5_path, "results/my_analysis/results_matrix")
+  rhdf5::h5writeAttribute(
+    attr = c("beta", "p.value"),
+    h5obj = h5_path,
+    name = "colnames",
+    h5loc = "results/my_analysis/results_matrix"
+  )
+
+  modelarray <- ModelArray(
+    h5_path,
+    scalar_types = c("FD"),
+    analysis_names = c("my_analysis")
+  )
+
+  expect_identical(sources(modelarray)$FD, c("attr_sub1", "attr_sub2"))
+  expect_identical(colnames(scalars(modelarray)$FD), c("attr_sub1", "attr_sub2"))
+  expect_identical(
+    colnames(results(modelarray)$my_analysis$results_matrix),
+    c("beta", "p.value")
+  )
+})
+
+test_that("ModelArray constructor prefers attributes over dataset column-name fallbacks", {
+  h5_path <- tempfile(fileext = ".h5")
+  on.exit(unlink(h5_path), add = TRUE)
+
+  rhdf5::h5createFile(h5_path)
+  rhdf5::h5createGroup(h5_path, "scalars")
+  rhdf5::h5createGroup(h5_path, "scalars/FD")
+  rhdf5::h5write(matrix(1:6, nrow = 3, ncol = 2), h5_path, "scalars/FD/values")
+  rhdf5::h5writeAttribute(
+    attr = c("attr_sub1", "attr_sub2"),
+    h5obj = h5_path,
+    name = "column_names",
+    h5loc = "scalars/FD/values"
+  )
+  # Conflicting fallback path should be ignored when attribute exists.
+  rhdf5::h5write(c("fallback_sub1", "fallback_sub2"), h5_path, "scalars/FD/column_names")
+
+  rhdf5::h5createGroup(h5_path, "results")
+  rhdf5::h5createGroup(h5_path, "results/my_analysis")
+  rhdf5::h5write(matrix(c(10, 20, 30, 40), nrow = 2, ncol = 2), h5_path, "results/my_analysis/results_matrix")
+  rhdf5::h5writeAttribute(
+    attr = c("attr_beta", "attr_p.value"),
+    h5obj = h5_path,
+    name = "colnames",
+    h5loc = "results/my_analysis/results_matrix"
+  )
+  # Conflicting fallback path should be ignored when attribute exists.
+  rhdf5::h5write(c("fallback_beta", "fallback_p"), h5_path, "results/my_analysis/column_names")
+
+  modelarray <- ModelArray(
+    h5_path,
+    scalar_types = c("FD"),
+    analysis_names = c("my_analysis")
+  )
+
+  expect_identical(sources(modelarray)$FD, c("attr_sub1", "attr_sub2"))
+  expect_identical(
+    colnames(results(modelarray)$my_analysis$results_matrix),
+    c("attr_beta", "attr_p.value")
+  )
+})
+
+test_that("analyseOneElement.lm covers predictor/collision/error branches", {
+  src <- c("s1", "s2", "s3", "s4")
+  modelarray <- methods::new(
+    "ModelArray",
+    sources = list(FD = src, FA = src),
+    scalars = list(
+      FD = matrix(c(1, 2, 3, 4, 2, 3, 4, 5), nrow = 2, byrow = TRUE),
+      FA = matrix(c(1, 2, NaN, 4, 5, 6, 7, 8), nrow = 2, byrow = TRUE)
+    ),
+    results = list(),
+    path = tempfile(fileext = ".h5")
+  )
+  phen <- data.frame(source_file = src, age = c(20, 30, 40, 50))
+
+  expect_error(
+    analyseOneElement.lm(
+      i_element = 1, formula = FD ~ age + FA, modelarray = modelarray,
+      phenotypes = transform(phen, FA = 1:4), scalar = "FD",
+      var.terms = c("estimate"), var.model = c("adj.r.squared"),
+      num.subj.lthr = 0, flag_initiate = TRUE
+    ),
+    "Column name collision"
+  )
+
+  modelarray_bad_src <- methods::new(
+    "ModelArray",
+    sources = list(FD = src, FA = c("x1", "x2", "x3", "x4")),
+    scalars = modelarray@scalars,
+    results = list(),
+    path = tempfile(fileext = ".h5")
+  )
+  expect_error(
+    analyseOneElement.lm(
+      i_element = 1, formula = FD ~ age + FA, modelarray = modelarray_bad_src,
+      phenotypes = phen, scalar = "FD",
+      var.terms = c("estimate"), var.model = c("adj.r.squared"),
+      num.subj.lthr = 0, flag_initiate = TRUE
+    ),
+    "do not match phenotypes\\$source_file"
+  )
+
+  # intersection threshold branch (finite FD but not enough finite FA)
+  res_insuf <- analyseOneElement.lm(
+    i_element = 1, formula = FD ~ age + FA, modelarray = modelarray,
+    phenotypes = phen, scalar = "FD",
+    var.terms = c("estimate"), var.model = c("adj.r.squared"),
+    num.subj.lthr = 3, num.stat.output = 3, flag_initiate = FALSE
+  )
+  expect_identical(res_insuf[1], 0)
+  expect_true(all(is.nan(res_insuf[-1])))
+
+  expect_warning(
+    res_skip <- analyseOneElement.lm(
+      i_element = 1, formula = FD ~ missingVar, modelarray = modelarray,
+      phenotypes = phen, scalar = "FD",
+      var.terms = c("estimate"), var.model = c("adj.r.squared"),
+      num.subj.lthr = 0, num.stat.output = 3, flag_initiate = FALSE,
+      on_error = "skip"
+    )
+  )
+  expect_identical(res_skip[1], 0)
+  expect_true(all(is.nan(res_skip[-1])))
+
+  expect_warning(
+    init_skip <- analyseOneElement.lm(
+      i_element = 1, formula = FD ~ missingVar, modelarray = modelarray,
+      phenotypes = phen, scalar = "FD",
+      var.terms = c("estimate"), var.model = c("adj.r.squared"),
+      num.subj.lthr = 0, flag_initiate = TRUE, on_error = "skip"
+    )
+  )
+  expect_true(is.nan(init_skip$column_names)[1])
+})
+
+test_that("analyseOneElement.gam covers predictor/collision/error branches", {
+  src <- c("s1", "s2", "s3", "s4", "s5")
+  modelarray <- methods::new(
+    "ModelArray",
+    sources = list(FD = src, FA = src),
+    scalars = list(
+      FD = matrix(c(1, 2, 3, 4, 5, 2, 3, 4, 5, 6), nrow = 2, byrow = TRUE),
+      FA = matrix(c(1, 2, NaN, 4, 5, 5, 6, 7, 8, 9), nrow = 2, byrow = TRUE)
+    ),
+    results = list(),
+    path = tempfile(fileext = ".h5")
+  )
+  phen <- data.frame(source_file = src, age = c(20, 30, 40, 50, 60))
+
+  expect_error(
+    analyseOneElement.gam(
+      i_element = 1, formula = FD ~ age + FA, modelarray = modelarray,
+      phenotypes = transform(phen, FA = 1:5), scalar = "FD",
+      var.smoothTerms = c("p.value"), var.parametricTerms = c("estimate"),
+      var.model = c("dev.expl"),
+      num.subj.lthr = 0, flag_initiate = TRUE
+    ),
+    "Column name collision"
+  )
+
+  modelarray_bad_src <- methods::new(
+    "ModelArray",
+    sources = list(FD = src, FA = c("x1", "x2", "x3", "x4", "x5")),
+    scalars = modelarray@scalars,
+    results = list(),
+    path = tempfile(fileext = ".h5")
+  )
+  expect_error(
+    analyseOneElement.gam(
+      i_element = 1, formula = FD ~ age + FA, modelarray = modelarray_bad_src,
+      phenotypes = phen, scalar = "FD",
+      var.smoothTerms = c("p.value"), var.parametricTerms = c("estimate"),
+      var.model = c("dev.expl"),
+      num.subj.lthr = 0, flag_initiate = TRUE
+    ),
+    "do not match phenotypes\\$source_file"
+  )
+
+  res_insuf <- analyseOneElement.gam(
+    i_element = 1, formula = FD ~ age + FA, modelarray = modelarray,
+    phenotypes = phen, scalar = "FD",
+    var.smoothTerms = c("p.value"), var.parametricTerms = c("estimate"),
+    var.model = c("dev.expl"),
+    num.subj.lthr = 4, num.stat.output = 3, flag_initiate = FALSE
+  )
+  expect_identical(res_insuf[1], 0)
+  expect_true(all(is.nan(res_insuf[-1])))
+
+  expect_warning(
+    res_skip <- analyseOneElement.gam(
+      i_element = 1, formula = FD ~ missingVar, modelarray = modelarray,
+      phenotypes = phen, scalar = "FD",
+      var.smoothTerms = c("p.value"), var.parametricTerms = c("estimate"),
+      var.model = c("dev.expl"),
+      num.subj.lthr = 0, num.stat.output = 3, flag_initiate = FALSE,
+      on_error = "skip"
+    )
+  )
+  expect_identical(res_skip[1], 0)
+  expect_true(all(is.nan(res_skip[-1])))
+})
+
+test_that("analyseOneElement.wrap covers coercion and error branches", {
+  src <- c("s1", "s2", "s3", "s4")
+  modelarray <- methods::new(
+    "ModelArray",
+    sources = list(FD = src, FA = src),
+    scalars = list(
+      FD = matrix(c(1, 2, 3, 4, 2, 3, 4, 5), nrow = 2, byrow = TRUE),
+      FA = matrix(c(11, 12, 13, 14, 21, 22, 23, 24), nrow = 2, byrow = TRUE)
+    ),
+    results = list(),
+    path = tempfile(fileext = ".h5")
+  )
+  phen <- data.frame(source_file = src, age = c(20, 30, 40, 50))
+
+  list_fun <- function(data) list(m = mean(data$FD), a = mean(data$FA))
+  init <- analyseOneElement.wrap(
+    i_element = 1, user_fun = list_fun, modelarray = modelarray,
+    phenotypes = phen, scalar = "FD",
+    num.subj.lthr = 0, flag_initiate = TRUE
+  )
+  expect_true(all(c("element_id", "m", "a") %in% init$column_names))
+
+  atomic_fun <- function(data) c(mean(data$FD), mean(data$FA))
+  out_atomic <- analyseOneElement.wrap(
+    i_element = 1, user_fun = atomic_fun, modelarray = modelarray,
+    phenotypes = phen, scalar = "FD",
+    num.subj.lthr = 0, num.stat.output = 3, flag_initiate = FALSE
+  )
+  expect_equal(length(out_atomic), 3)
+
+  expect_error(
+    analyseOneElement.wrap(
+      i_element = 1, user_fun = function(data) data.frame(a = 1:2), modelarray = modelarray,
+      phenotypes = phen, scalar = "FD",
+      num.subj.lthr = 0, flag_initiate = FALSE, num.stat.output = 2
+    ),
+    "must return a one-row data.frame/tibble"
+  )
+
+  expect_error(
+    analyseOneElement.wrap(
+      i_element = 1, user_fun = function(data) new.env(), modelarray = modelarray,
+      phenotypes = phen, scalar = "FD",
+      num.subj.lthr = 0, flag_initiate = FALSE, num.stat.output = 2
+    ),
+    "Unsupported return type"
+  )
+
+  expect_warning(
+    out_skip <- analyseOneElement.wrap(
+      i_element = 1, user_fun = function(data) stop("boom"), modelarray = modelarray,
+      phenotypes = phen, scalar = "FD",
+      num.subj.lthr = 0, flag_initiate = FALSE, num.stat.output = 3, on_error = "skip"
+    )
+  )
+  expect_identical(out_skip[1], 0)
+  expect_true(all(is.nan(out_skip[-1])))
+
+  expect_error(
+    analyseOneElement.wrap(
+      i_element = 1, user_fun = list_fun, modelarray = modelarray,
+      phenotypes = transform(phen, FD = 1:4), scalar = "FD",
+      num.subj.lthr = 0, flag_initiate = TRUE
+    ),
+    "Column name collision"
+  )
+})
+
+test_that("writeResults handles validation, non-numeric LUT, and overwrite=false", {
+  h5_path <- tempfile(fileext = ".h5")
+  on.exit(unlink(h5_path), add = TRUE)
+
+  expect_error(
+    writeResults(h5_path, df.output = c(1, 2, 3), analysis_name = "bad"),
+    "must be data of type `data.frame`"
+  )
+
+  df1 <- data.frame(
+    element_id = 0:1,
+    score = c(1.1, 2.2),
+    group = c("A", "B")
+  )
+  writeResults(h5_path, df.output = df1, analysis_name = "my_res", overwrite = TRUE)
+
+  h5 <- hdf5r::H5File$new(h5_path, mode = "r")
+  expect_true(h5[["results/my_res"]]$exists("results_matrix"))
+  expect_true(h5[["results/my_res"]]$exists("column_names"))
+  expect_true(h5[["results/my_res"]]$exists("lut_forcol3"))
+  mat1 <- h5[["results/my_res/results_matrix"]]$read()
+  h5$close_all()
+
+  df2 <- data.frame(element_id = 0:2, score = c(10, 20, 30))
+  expect_warning(
+    writeResults(h5_path, df.output = df2, analysis_name = "my_res", overwrite = FALSE),
+    "exists but not to overwrite"
+  )
+
+  h5b <- hdf5r::H5File$new(h5_path, mode = "r")
+  mat2 <- h5b[["results/my_res/results_matrix"]]$read()
+  h5b$close_all()
+  expect_identical(dim(mat2), dim(mat1))
+})
+
+test_that("ModelArray.wrap validation branches are exercised", {
+  src <- c("s1", "s2", "s3", "s4")
+  modelarray <- methods::new(
+    "ModelArray",
+    sources = list(FD = src),
+    scalars = list(FD = matrix(c(1, 2, 3, 4, 2, 3, 4, 5), nrow = 2, byrow = TRUE)),
+    results = list(),
+    path = tempfile(fileext = ".h5")
+  )
+  phen <- data.frame(source_file = src, age = c(20, 30, 40, 50))
+  simple_fun <- function(data) data.frame(m = mean(data$FD))
+
+  expect_error(
+    ModelArray.wrap(simple_fun, data = list(), phenotypes = phen, scalar = "FD", element.subset = as.integer(1)),
+    "data's class is not ModelArray"
+  )
+  expect_error(
+    ModelArray.wrap(simple_fun, data = modelarray, phenotypes = phen, scalar = "FD", element.subset = c(1, 2)),
+    "Please enter integers for element.subset"
+  )
+  expect_error(
+    ModelArray.wrap(simple_fun, data = modelarray, phenotypes = subset(phen, select = -source_file), scalar = "FD", element.subset = as.integer(1)),
+    "Did not find column 'source_file'"
+  )
+  expect_error(
+    ModelArray.wrap(simple_fun, data = modelarray, phenotypes = rbind(phen, phen[1, ]), scalar = "FD", element.subset = as.integer(1)),
+    "not the same as that in ModelArray 'data'"
+  )
+
+  # non-identical but matchable source_file order should be handled
+  phen_swapped <- phen[c(2, 1, 3, 4), ]
+  out <- ModelArray.wrap(
+    simple_fun, data = modelarray, phenotypes = phen_swapped, scalar = "FD",
+    element.subset = as.integer(c(1, 2)), n_cores = 1, pbar = FALSE, verbose = FALSE,
+    num.subj.lthr.abs = 0
+  )
+  expect_equal(out$element_id, c(0, 1))
+})
diff --git a/tests/testthat/test-ModelArray_subj_specific_masks.R b/tests/testthat/test-ModelArray_subj_specific_masks.R
index f7276f2..3c9ecd5 100644
--- a/tests/testthat/test-ModelArray_subj_specific_masks.R
+++ b/tests/testthat/test-ModelArray_subj_specific_masks.R
@@ -40,7 +40,6 @@ test_that("ModelArray handles subject-specific masks as expected", {
   }
 
 
-
   df.config <- data.frame(
     nsubj = c(60, 40, 30),
     specials = c("", "", "allInvalidValues")
diff --git a/vignettes/doc_for_developer.Rmd b/vignettes/doc_for_developer.Rmd
index b681e83..df2a9d3 100644
--- a/vignettes/doc_for_developer.Rmd
+++ b/vignettes/doc_for_developer.Rmd
@@ -141,6 +141,33 @@ devtools::load_all()
 
 Then click "Run Test" button in RStudio to run the test file. Check if there is anything failed.
 
+### Running code styling with `styler`
+Use `styler` to format all package code before committing changes.
+
+Install `styler` if needed:
+
+```{r, eval=FALSE}
+if (!requireNamespace("styler", quietly = TRUE)) {
+  install.packages("styler")
+}
+```
+
+Apply code styling to the package using the tidyverse style guide:
+
+```{r, eval=FALSE}
+styler::style_pkg(style = styler::tidyverse_style)
+```
+
+Optionally, run a dry check to confirm no files would be reformatted:
+
+```{r, eval=FALSE}
+result <- styler::style_pkg(
+  dry = "on",
+  style = styler::tidyverse_style
+)
+all(result[["changed"]] == FALSE)
+```
+
 For more details on how to write unit tests + test out, please check out `Testing` chapters in the book ["R Packages" written by Hadley Wickham, Jennifer Bryan](https://r-pkgs.org/index.html)
 
 ## Building the Docker image for ModelArray + ModelArrayIO
diff --git a/vignettes/exploring-h5.Rmd b/vignettes/exploring-h5.Rmd
index fd057e8..6c21974 100644
--- a/vignettes/exploring-h5.Rmd
+++ b/vignettes/exploring-h5.Rmd
@@ -44,8 +44,8 @@ the file path, how many sources, which scalars are loaded, and which analyses ha
 
 ```{r sources, eval=FALSE}
 src <- sources(modelarray)[["FDC"]]
-length(src)     # number of sources
-head(src)       # first few filenames
+length(src) # number of sources
+head(src) # first few filenames
 ```
 
 ```{.console}
@@ -189,6 +189,6 @@ age_pvals <- as.numeric(lm_results[, col_idx])
 sum(age_pvals < 0.05, na.rm = TRUE)
 
 # Top 10 most significant element IDs
-top_elements <- order(age_pvals)[1:10] - 1  # convert to 0-based element IDs
+top_elements <- order(age_pvals)[1:10] - 1 # convert to 0-based element IDs
 top_elements
 ```
diff --git a/vignettes/modelling.Rmd b/vignettes/modelling.Rmd
index e499fdb..65dd68b 100644
--- a/vignettes/modelling.Rmd
+++ b/vignettes/modelling.Rmd
@@ -170,7 +170,7 @@ result <- gen_gamFormula_fxSmooth(
   fx = TRUE, k = 4
 )
 formula <- result$formula
-phenotypes <- result$phenotypes  # may contain new ordered factor column
+phenotypes <- result$phenotypes # may contain new ordered factor column
 ```
 
 **Continuous interaction** (`gen_gamFormula_contIx`):