Enhance com2links to subset matrices and improve error handling

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: BGmisc
 Title: An R Package for Extended Behavior Genetics Analysis
-Version: 1.7.0.1
+Version: 1.7.0.1.1
 Authors@R: c(
     person("S. Mason", "Garrison", , "garrissm@wfu.edu", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-4804-6003")),

NEWS.md

@@ -4,6 +4,8 @@
 
 * Optimized sliceFamilies to be more abstract
 * Created `.require_openmx()` to make it easier to use OpenMx functions without making OpenMx a dependency
+* Smarter string ID handling for ped2id
+* Fixed how different-sized matrices are handled by `com2links()`
 
 # BGmisc 1.7.0.0
 * Fixed bug in parList

R/makeLinks.R

@@ -79,22 +79,50 @@ com2links <- function(
 
   # Extract individual IDs from the first available matrix.
   ids <- NULL
+  # Find the smallest matrix by ncol (avoids extracting IDs from large matrices).
+  mat_refs <- list()
+  if (!is.null(ad_ped_matrix)) mat_refs[["ad"]] <- ncol(ad_ped_matrix)
+  if (!is.null(mit_ped_matrix)) mat_refs[["mt"]] <- ncol(mit_ped_matrix)
+  if (!is.null(cn_ped_matrix)) mat_refs[["cn"]] <- ncol(cn_ped_matrix)
+
+  if (length(mat_refs) == 0L) {
+    stop("At least one relationship matrix must be provided.")
+  }
 
+  smallest <- names(which.min(unlist(mat_refs)))
+  guide_mat <- switch(smallest,
+    ad = ad_ped_matrix,
+    mt = mit_ped_matrix,
+    cn = cn_ped_matrix
+  )
 
-  if (!is.null(cn_ped_matrix)) {
-    ids <- as.numeric(dimnames(cn_ped_matrix)[[1]])
-    nc <- ncol(cn_ped_matrix)
-  } else if (!is.null(ad_ped_matrix)) {
-    ids <- as.numeric(dimnames(ad_ped_matrix)[[1]])
-    nc <- ncol(ad_ped_matrix)
-  } else if (!is.null(mit_ped_matrix)) {
-    ids <- as.numeric(dimnames(mit_ped_matrix)[[1]])
-    nc <- ncol(mit_ped_matrix)
+  # Extract IDs only from the smallest matrix.
+  guide_ids <- dimnames(guide_mat)[[1]]
+  if (is.null(guide_ids) || length(guide_ids) == 0L) {
+    stop("Could not extract IDs from the smallest matrix.")
+  }
+  ids <- suppressWarnings(as.numeric(guide_ids))
+  if (anyNA(ids)) {
+    warning(
+      "Matrix dimnames used as IDs should be strictly numeric for 'com2links()'. Found non-numeric or non-coercible IDs in the smallest matrix."
+    )
   }
+  nc <- length(ids)
 
-  if (is.null(ids)) {
-    stop("Could not extract IDs from the provided matrices.")
+  # Subset only the larger matrices to match the smallest matrix's IDs and ordering.
+  if (!is.null(ad_ped_matrix) && ncol(ad_ped_matrix) > nc) {
+    if (verbose) message("Subsetting ad_ped_matrix from ", ncol(ad_ped_matrix), " to ", nc, " IDs.")
+    ad_ped_matrix <- ad_ped_matrix[guide_ids, guide_ids, drop = FALSE]
+  }
+  if (!is.null(mit_ped_matrix) && ncol(mit_ped_matrix) > nc) {
+    if (verbose) message("Subsetting mit_ped_matrix from ", ncol(mit_ped_matrix), " to ", nc, " IDs.")
+    mit_ped_matrix <- mit_ped_matrix[guide_ids, guide_ids, drop = FALSE]
   }
+  if (!is.null(cn_ped_matrix) && ncol(cn_ped_matrix) > nc) {
+    if (verbose) message("Subsetting cn_ped_matrix from ", ncol(cn_ped_matrix), " to ", nc, " IDs.")
+    cn_ped_matrix <- cn_ped_matrix[guide_ids, guide_ids, drop = FALSE]
+  }
+
 
   # --- matrix_case construction and switch dispatch ---
   matrix_case <- paste(sort(c(

R/segmentPedigree.R

@@ -11,6 +11,7 @@
 #' @param dadID character.  Name of the column in ped for the father ID variable
 #' @param famID character.  Name of the column to be created in ped for the family ID variable
 #' @param twinID character.  Name of the column in ped for the twin ID variable, if applicable
+#' @param overwrite logical.  If TRUE, will overwrite existing famID variable if it exists.  Default is TRUE.
 #' @param ... additional arguments to be passed to \code{\link{ped2com}}
 #' @details
 #' The general idea of this function is to use person ID, mother ID, and father ID to
@@ -31,18 +32,20 @@
 ped2fam <- function(ped, personID = "ID",
                     momID = "momID", dadID = "dadID", famID = "famID",
                     twinID = "twinID",
+                    overwrite = TRUE,
                     ...) {
   # Call to wrapper function
   .ped2id(
     ped = ped, personID = personID, momID = momID, dadID = dadID, famID = famID, twinID = twinID,
-    type = "parents"
+    type = "parents",
+    overwrite = overwrite
   )
 }
 
 .ped2id <- function(ped,
                     personID = "ID", momID = "momID", dadID = "dadID",
                     famID = "famID", twinID = "twinID",
-                    type,
+                    type, overwrite = TRUE,
                     ...) {
   # Turn pedigree into family
   pg <- ped2graph(
@@ -55,23 +58,44 @@
 
   # Create famID data.frame
   # Convert IDs to numeric, with warning if coercion collapses IDs
+
   uniques <- suppressWarnings(unique(as.numeric(names(wcc$membership))))
+  keep_string <- FALSE
 
   if (length(uniques) == 1L && is.na(uniques)) {
     warning("After converting IDs to numeric, all IDs became NA. This indicates ID coercion collapsed IDs. Please ensure IDs aren't character or factor variables.")
-
+ keep_string <- TRUE
+     } else if (length(uniques) < length(wcc$membership)) {
+     warning("After converting IDs to numeric, some IDs became NA. This indicates ID coercion collapsed some IDs. Please ensure IDs aren't character or factor variables.")
+  keep_string <- TRUE
+}
+if(keep_string==TRUE) {
     fam <- data.frame(
       V1 = names(wcc$membership),
       V2 = wcc$membership
     )
-  } else {
+    } else {
     fam <- data.frame(
       V1 = as.numeric(names(wcc$membership)),
       V2 = wcc$membership
     )
   }
 
   names(fam) <- c(personID, famID)
+
+  if(famID %in% names(ped)) {
+    if(overwrite) {
+    overwrite_message <-  "be overwritten."
+    ped[[famID]] <- NULL
+    } else {
+    overwrite_message <-  "not be overwritten."
+    fam[[famID]] <- NULL
+   }
+
+    warning(sprintf("The famID variable '%s' already exists in the pedigree. The existing variable will %s", famID, overwrite_message))
+
+   }
+
   ped2 <- merge(fam, ped,
     by = personID, all.x = FALSE, all.y = TRUE
   )

tests/testthat/test-makeLinks.R

@@ -324,3 +324,82 @@ test_that("com2links garbage collection does not affect output, using two compon
 
   expect_equal(result_gc, result_no_gc)
 })
+
+
+test_that("com2links handles mismatched matrix dimensions by subsetting to smallest", {
+  data(hazard)
+  subset_ids <- hazard$ID[seq_len(ceiling(nrow(hazard) / 2))]
+
+  ad_small <- ped2add(hazard, sparse = TRUE, keep_ids = subset_ids)
+  mit_ped_matrix <- ped2mit(hazard, sparse = TRUE)
+  cn_ped_matrix <- ped2cn(hazard, sparse = TRUE)
+
+  # All three matrices, ad is smaller
+  result_mismatch <- com2links(
+    ad_ped_matrix = ad_small,
+    mit_ped_matrix = mit_ped_matrix,
+    cn_ped_matrix = cn_ped_matrix,
+    writetodisk = FALSE
+  )
+
+  # Reference: all three matrices built from the same subset
+  result_ref <- com2links(
+    ad_ped_matrix = ad_small,
+    mit_ped_matrix = ped2mit(hazard, sparse = TRUE, keep_ids = subset_ids),
+    cn_ped_matrix = ped2cn(hazard, sparse = TRUE, keep_ids = subset_ids),
+    writetodisk = FALSE
+  )
+
+  expect_equal(result_mismatch, result_ref)
+
+  # Only IDs from the smaller matrix should appear
+  all_output_ids <- unique(c(result_mismatch$ID1, result_mismatch$ID2))
+  expect_true(all(all_output_ids %in% as.numeric(dimnames(ad_small)[[1]])))
+})
+
+test_that("com2links mismatched dimensions with two matrices", {
+  data(hazard)
+  subset_ids <- hazard$ID[seq_len(ceiling(nrow(hazard) / 2))]
+
+  ad_ped_matrix <- ped2add(hazard, sparse = TRUE)
+  cn_small <- ped2cn(hazard, sparse = TRUE, keep_ids = subset_ids)
+
+  # cn is smaller than ad
+  result_mismatch <- com2links(
+    ad_ped_matrix = ad_ped_matrix,
+    cn_ped_matrix = cn_small,
+    writetodisk = FALSE
+  )
+
+  result_ref <- com2links(
+    ad_ped_matrix = ped2add(hazard, sparse = TRUE, keep_ids = subset_ids),
+    cn_ped_matrix = cn_small,
+    writetodisk = FALSE
+  )
+
+  expect_equal(result_mismatch, result_ref)
+})
+
+test_that("com2links mismatched dimensions with mit smaller", {
+  data(hazard)
+  subset_ids <- hazard$ID[seq_len(ceiling(nrow(hazard) / 2))]
+
+  ad_ped_matrix <- ped2add(hazard, sparse = TRUE)
+  mit_small <- ped2mit(hazard, sparse = TRUE, keep_ids = subset_ids)
+
+  result_mismatch <- com2links(
+    ad_ped_matrix = ad_ped_matrix,
+    mit_ped_matrix = mit_small,
+    writetodisk = FALSE
+  )
+
+  result_ref <- com2links(
+    ad_ped_matrix = ped2add(hazard, sparse = TRUE, keep_ids = subset_ids),
+    mit_ped_matrix = mit_small,
+    writetodisk = FALSE
+  )
+
+  expect_equal(result_mismatch, result_ref)
+  expect_true(all(unique(c(result_mismatch$ID1, result_mismatch$ID2)) %in% as.numeric(dimnames(mit_small)[[1]])))
+  expect_true(all(unique(c(result_mismatch$ID1, result_mismatch$ID2)) %in% as.numeric(dimnames(ad_ped_matrix)[[1]])))
+})

tests/testthat/test-segmentPedigree.R

@@ -42,6 +42,40 @@ test_that("ped2fam gets the right families for inbreeding data", {
   expect_equal(ds$famID, ds$newFamID)
 })
 
+test_that("ped2fam overwrite = TRUE warns and replaces existing famID", {
+  data(hazard)
+  # Add a fake existing famID column with wrong values
+  hazard$famID <- 999L
+
+  # Should warn that famID will be overwritten
+  expect_warning(
+    ds <- ped2fam(hazard, famID = "famID", overwrite = TRUE),
+    regexp = "already exists.*be overwritten"
+  )
+
+  # The result should have the correctly computed famID, not the placeholder
+  expect_false(all(ds$famID == 999L))
+  # No duplicate columns
+  expect_equal(sum(names(ds) == "famID"), 1L)
+})
+
+test_that("ped2fam overwrite = FALSE warns and preserves existing famID", {
+  data(hazard)
+  # Add a known existing famID column
+  hazard$famID <- 999L
+
+  # Should warn that famID will NOT be overwritten
+  expect_warning(
+    ds <- ped2fam(hazard, famID = "famID", overwrite = FALSE),
+    regexp = "already exists.*not be overwritten"
+  )
+
+  # The result should preserve the original famID values
+  expect_true(all(ds$famID == 999L))
+  # No duplicate columns
+  expect_equal(sum(names(ds) == "famID"), 1L)
+})
+
 test_that("ped2graph produces a graph for hazard data with mothers", {
   expect_silent(data(hazard))
   g <- ped2graph(hazard, adjacent = "mothers")

vignettes/articles/tutorialmanuscript.Xmd

@@ -7,7 +7,7 @@ author:
     corresponding: true
     email: "garrissm@wfu.edu"
 abstract: |
-  Twin studies remain the dominant design in behavior genetics, yet most twin half-siblings, cousins, and multi-generational relatives whose distinct kinship coefficients jointly identify a richer set of variance components than any MZ/DZ comparison alone. We demonstrate how to fit extended pedigree models using the BGmisc package and OpenMx. 
+  Twin studies remain the dominant design in behavior genetics, yet most twin half-siblings, cousins, and multi-generational relatives whose distinct kinship coefficients jointly identify a richer set of variance components than any MZ/DZ comparison alone. We demonstrate how to fit extended pedigree models using the BGmisc package and OpenMx.
   We apply the extended pedigree model to mutiple datasets
   of Youth (a large human panel study with researcher-linked kinship), the Kluane Red Squirrel Project
   (a multi-generational animal field study), and a children-of-twins dataset.
@@ -31,8 +31,8 @@ vignette: >
   %\VignetteIndexEntry{Extended Family Modeling with BGmisc}
   %\VignetteEncoding{UTF-8}
   %\VignetteEngine{knitr::rmarkdown}
-editor_options: 
-  markdown: 
+editor_options:
+  markdown:
     wrap: 100
 ---
 
@@ -65,7 +65,8 @@ studies, either intentionally (e.g., twin registries that also include siblings
 byproduct of large panel studies (e.g., the National Longitudinal Survey of Youth, which includes
 researcher-linked kinship). In most cases, the additional relatives are excluded from analysis, and
 the twin design is applied to a subset of the data, even though these relatives carry independent
-information about the genetic and environmental architecture of the phenotype.
+information about the genetic and environmental architecture of the phenotype. For example, many twin registries also include triplets, siblings, children, and parents, all of whom provide additional information about the genetic and environmental architecture of the phenotype.
+<!-- cite https://helda.helsinki.fi/server/api/core/bitstreams/f0b6dc08-69df-449e-a8fe-e2c78abf7f60/content -->
 
 The extended pedigree model, which we have introduced elsewhere (see ETC), leverages the full range
 of kinship coefficients in a pedigree to identify a richer set of variance components than the
@@ -87,15 +88,15 @@ value provides independent leverage for disentangling genetic from environmental
 the number of distinct kinship types increases, so does the number of identifiable variance
 components.
 
-Extended pedigree designs have been used in behavior genetics since at least the 1970s [@eaves1978; @fulker_multiple_1988], but they have remained a minority practice. Partially over concerns about model identification and power (Wilson, 1982, 1989), the complexity of fitting these models, and the relative costs of collecting twin data compared to extended family data. 
+Extended pedigree designs have been used in behavior genetics since at least the 1970s [@eaves1978; @fulker_multiple_1988], but they have remained a minority practice. Partially over concerns about model identification and power (Wilson, 1982, 1989), the complexity of fitting these models, and the relative costs of collecting twin data compared to extended family data.
 
-<! -- https://onlinelibrary.wiley.com/doi/10.1002/bimj.4710310511 
+<!-- https://onlinelibrary.wiley.com/doi/10.1002/bimj.4710310511 -->
 
 but also because the twin design has been so successful and widely adopted. The twin design is often seen as the "gold standard" in behavior genetics, and many researchers may be hesitant to deviate from this established approach. Additionally, many human datasets simply do not include the necessary family structure to fit extended pedigree models, which may limit their applicability in certain contexts.
 
 <the reasosn are numerous for why this is the case, but a key factor is that many human datasets simply do not include the necessary family structure to fit these models. And the twin design is often the default analytic approach, even when more complex family data are available.
 
-Deriving  
+
 
 In contrast, similar
 models are common in plant and animal breeding, where pedigree data is more routinely collected and