diff --git a/NEWS.md b/NEWS.md index 92292f4a..68aeebd8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,10 @@ # BGmisc NEWS # Development version: 1.6.0.9000 +Add option for MZ twins in the additive genetic matrix +Add option to select sex for MZ twin generation. +Add option to tweak pedigree with one id provided + # BGmisc 1.6.0.1 * Add helper functions for checkParents etc * fixed incorrect direction so that parents are pointing to children in the graphs diff --git a/R/buildComponent.R b/R/buildComponent.R index 38a96522..9c86a8fa 100644 --- a/R/buildComponent.R +++ b/R/buildComponent.R @@ -19,6 +19,9 @@ #' @param isChild_method character. The method to use for computing the isChild matrix. Options are "classic" or "partialparent" #' @param adjBeta_method numeric The method to use for computing the building the adjacency_method matrix when using the "beta" build #' @param compress logical. If TRUE, use compression when saving the checkpoint files. Defaults to TRUE. +#' @param mz_twins logical. If TRUE, merge MZ co-twin columns in the r2 matrix before tcrossprod so that MZ twins are coded with relatedness 1 instead of 0.5. Twin pairs are identified from the \code{twinID} column. When a \code{zygosity} column is also present, only pairs where both members have \code{zygosity == "MZ"} are used; otherwise all \code{twinID} pairs are assumed to be MZ. Defaults to FALSE. +#' @param mz_method character. The method to handle MZ twins. Options are "merging" (default) or "addtwins". "addtwins" adds the twin2 column to the twin1 column before tcrossprod so that all relatedness flows through a single source, then leaves the twin2 column as zero and relies on the fact that the row/col names are the same to copy the values back to twin2 after tcrossprod. "merging" merges the twin2 column into the twin1 column before tcrossprod and then copies the values back to twin2 after tcrossprod so that both twins appear in the final matrix. +#' @param beta logical. Used for benchmarking #' @param ... additional arguments to be passed to \code{\link{ped2com}} #' @details The algorithms and methodologies used in this function are further discussed and exemplified in the vignette titled "examplePedigreeFunctions". For more advanced scenarios and detailed explanations, consult this vignette. #' @export @@ -32,7 +35,7 @@ ped2com <- function(ped, component, standardize_colnames = TRUE, transpose_method = "tcrossprod", adjacency_method = "direct", - isChild_method = "classic", + isChild_method = "partialparent", saveable = FALSE, resume = FALSE, save_rate = 5, @@ -42,6 +45,9 @@ ped2com <- function(ped, component, save_path = "checkpoint/", adjBeta_method = NULL, compress = TRUE, + mz_twins = TRUE, + mz_method = "addtwins", + beta = FALSE, ...) { #------ # Check inputs @@ -121,12 +127,37 @@ ped2com <- function(ped, component, ped <- standardizeColnames(ped, verbose = config$verbose) } + mz_row_pairs <- NULL + mz_id_pairs <- NULL + + if (mz_twins == TRUE && "twinID" %in% colnames(ped)) { + df_mz <- findMZtwins(ped, + verbose = config$verbose, + returnIDs = TRUE, + returnRows = TRUE, + returnAsList = TRUE, + beta = beta + ) + mz_row_pairs <- df_mz$pair_rows + mz_id_pairs <- df_mz$pair_ids + } + + # Load final result if computation was completed if (config$resume == TRUE && file.exists(checkpoint_files$final_matrix)) { if (config$verbose == TRUE) cat("Loading final computed matrix...\n") return(readRDS(checkpoint_files$final_matrix)) } + if (mz_method %in% c("merging") && mz_twins == TRUE && !is.null(mz_row_pairs) && length(mz_row_pairs) > 0 && + config$component %in% c("additive")) { + # replace all MZ twin IDs with the first twin's ID in each pair so they are merged for the path tracing and all subsequent steps. We will copy the values back to the second twin at the end. + ped <- fuseTwins(ped = ped, mz_row_pairs = mz_row_pairs, mz_id_pairs = mz_id_pairs, config = config, beta = beta) + if (config$verbose == TRUE) { + message("Merged ", length(mz_row_pairs), " MZ twin pair(s) in pedigree dataset for path tracing") + } + } + #------ # Algorithm @@ -141,6 +172,7 @@ ped2com <- function(ped, component, cat(paste0("Family Size = ", config$nr, "\n")) } + # # Step 1: Construct parent-child adjacency matrix ## A. Resume from Checkpoint if Needed @@ -195,6 +227,8 @@ ped2com <- function(ped, component, config = config, compress = config$compress ) + + # TODO merge twin columns # --- Step 2: Compute Relatedness Matrix --- @@ -226,7 +260,7 @@ ped2com <- function(ped, component, # r is I + A + A^2 + ... = (I-A)^-1 from RAM # could trim, here ## it keeps going until it explains all of the relatedness with themselves (i.e., mtSum == 0) - # some of this precision is articifuial because we literally get to the point that the condon is eaither there or not. probabiliticy + # some of this precision is artificial because we literally get to the point that the condon is eaither there or not. probabiliticy # how much percision do we need to get unbiased estimates @@ -281,6 +315,29 @@ ped2com <- function(ped, component, compress = config$compress ) + if (mz_method == "addtwins" && mz_twins == TRUE && !is.null(mz_row_pairs) && length(mz_row_pairs) > 0) { + if (config$verbose == TRUE) { + message("MZ twin merging enabled: Will merge MZ twin columns in r2 before tcrossprod") + } + + # --- Step 3b: Add --- + # MZ twins share the same genetic source. We absorb twin2's column into + # twin1's before tcrossprod so all path-traced relatedness flows through a + # single source. After tcrossprod we copy twin1's row/col back to twin2. + if (!is.null(mz_row_pairs) && length(mz_row_pairs) > 0 && config$component %in% c("additive")) { + # Extract all indices at once for batch operations + pairs_mat <- do.call(rbind, mz_row_pairs) + idx1_all <- pairs_mat[, 1] + idx2_all <- pairs_mat[, 2] + # Batch: absorb all twin2 columns into twin1 columns, then zero twin2 + r2[, idx1_all] <- r2[, idx1_all, drop = FALSE] + r2[, idx2_all, drop = FALSE] + r2[, idx2_all] <- 0 + + if (config$verbose == TRUE) { + message("Added ", length(mz_row_pairs), " MZ twin pair column(s) in r2") + } + } + } # --- Step 4: T crossproduct --- if (config$resume == TRUE && file.exists(checkpoint_files$tcrossprod_checkpoint) && @@ -300,12 +357,68 @@ ped2com <- function(ped, component, } } + if (mz_method %in% c("merging", "addtwins") && mz_twins == TRUE && config$component %in% c("additive") && !is.null(mz_row_pairs) && length(mz_row_pairs) > 0) { + # --- Step 4b: Restore MZ twins --- + # Copy twin1's row/col to twin2 so both twins appear in the final matrix. + if (config$sparse == FALSE) { + r <- as.matrix(r) + rnames <- rownames(r) + ids_mat <- do.call(rbind, mz_id_pairs) + idx1_all <- match(ids_mat[, 1], rnames) + idx2_all <- match(ids_mat[, 2], rnames) + # Batch copy: twin1 rows/cols -> twin2 rows/cols + r[idx2_all, ] <- r[idx1_all, ] + + r[, idx2_all] <- r[, idx1_all] + } else { + # TODO this is really slow. Can we do it without coercing to dense? Maybe by doing row/col replacement on the sparse matrix directly? Or by constructing a sparse matrix with the twin2 values and adding it to r? + # r <- df_add + + rnames <- r@Dimnames[[1]] + + ids_mat <- do.call(rbind, mz_id_pairs) + # needs to use sparse indexing to avoid coercion to dense + idx1_all <- match(ids_mat[, 1], rnames) + idx2_all <- match(ids_mat[, 2], rnames) + + twin1_rows <- r[idx1_all, , drop = FALSE] + twin1_cols <- r[, idx1_all, drop = FALSE] + twin1_rows@Dimnames[[1]] <- rnames[idx2_all] + twin1_cols@Dimnames[[2]] <- rnames[idx2_all] + twin1_self <- r[idx1_all, idx1_all, drop = FALSE] + twin1_self@Dimnames[[1]] <- rnames[idx2_all] + + r[idx2_all, ] <- twin1_rows + r[, idx2_all] <- twin1_cols + r[idx2_all, idx2_all] <- twin1_self + + # Batch copy: twin1 rows/cols -> twin2 rows/cols + + # Row/column replacement on a dsCMatrix (symmetric) causes Matrix to + # coerce to dgCMatrix (general), doubling stored entries. Convert back + + r <- Matrix::drop0(r) + + # so both mz_method paths return the same sparse class. + if (methods::is(r, "CsparseMatrix") && !methods::is(r, "symmetricMatrix")) { + r <- Matrix::forceSymmetric(r) + } + } + if (config$verbose == TRUE) { + message("Restored ", length(mz_row_pairs), " MZ twin pair(s) in relatedness matrix") + } + } + + if (config$component %in% c("mitochondrial", "mtdna", "mitochondria")) { r@x <- rep(1, length(r@x)) # Assign 1 to all nonzero elements for mitochondrial component } - if (config$sparse == FALSE) { + # Remove explicit zeros so that both mz_method paths produce + # structurally identical sparse matrices + + if (config$sparse == FALSE && !methods::is(r, "matrix")) { r <- as.matrix(r) } # flattens diagonal if you don't want to deal with inbreeding @@ -335,6 +448,8 @@ ped2add <- function(ped, max_gen = 25, sparse = TRUE, verbose = FALSE, save_rate_parlist = 100000 * save_rate, save_path = "checkpoint/", compress = TRUE, + mz_twins = FALSE, + mz_method = "addtwins", ...) { ped2com( ped = ped, @@ -353,6 +468,8 @@ ped2add <- function(ped, max_gen = 25, sparse = TRUE, verbose = FALSE, save_rate_parlist = save_rate_parlist, save_path = save_path, compress = compress, + mz_twins = mz_twins, + mz_method = mz_method, ... ) } diff --git a/R/constructAdjacency.R b/R/constructAdjacency.R index 3053fb5b..fe3f7bbc 100644 --- a/R/constructAdjacency.R +++ b/R/constructAdjacency.R @@ -70,7 +70,8 @@ .adjIndexed <- function(ped, component, saveable, resume, save_path, verbose, lastComputed, checkpoint_files, update_rate, - parList, lens, save_rate_parlist, config, compress = config$compress) { + parList, lens, save_rate_parlist, + config, compress = config$compress) { # Loop through each individual in the pedigree # Build the adjacency matrix for parent-child relationships # Is person in column j the parent of the person in row i? .5 for yes, 0 for no. @@ -207,7 +208,7 @@ return(list_of_adjacency) } -#' Construct Adjacency Matrix for Parent-Child Relationships Using Beta Method +#' Construct Adjacency Matrix for Parent-Child Relationships Using Beta Methods #' This function constructs an adjacency matrix for parent-child relationships #' using a method in beta testing. It identifies parent-child pairs based on the #' specified component of relatedness. diff --git a/R/helpTwins.R b/R/helpTwins.R new file mode 100644 index 00000000..b8a2c33b --- /dev/null +++ b/R/helpTwins.R @@ -0,0 +1,278 @@ +#' Determine isTwin Status +#' @param ped pedigree data frame +#' @return isTwin 'S' matrix +#' @keywords internal + + +isTwin <- function(ped) { + isTwin <- apply(ped[, c("twinID")], 1, function(x) { + !is.na(x) + }) +} + + +#' Find MZ twin pair_rows in a pedigree +#' +#' Identifies MZ twin pair_rows from the \code{twinID} column and returns their +#' row indices. These indices are used later to merge the twins' columns in +#' the \code{r2} matrix before \code{tcrossprod}, which correctly produces +#' relatedness 1 between MZ co-twins with no diagonal or downstream artifacts. +#' +#' @param ped A pedigree data.frame with columns \code{ID} and \code{twinID}. +#' Optionally a \code{zygosity} column; when present only pair_rows where both +#' members have \code{zygosity == "MZ"} are used. +#' @param verbose logical. If TRUE, print progress messages. +#' @param returnIDs logical. If TRUE, return the IDs of the twin pair_rows instead of row indices. +#' @param returnRows logical. If TRUE, return the row indices of the twin pair_rows instead of IDs. +#' @param returnAsList logical. If TRUE, return results as a list of vectors +#' (default). If FALSE, return results as a data.frame with separate columns for each twin's ID and row index. +#' @param beta logical. If TRUE, use an optimized approach with O(1) lookups for large pedigrees. If FALSE (default), use a simpler approach that may be less efficient for large pedigrees. +#' @return A list of length-2 integer vectors \code{c(idx1, idx2)} giving the +#' row indices of each MZ pair in the pedigree, or \code{NULL} if none found. +#' @keywords internal +findMZtwins <- function(ped, verbose = FALSE, returnRows = TRUE, + returnIDs = FALSE, returnAsList = TRUE, + beta = FALSE) { + if (!"twinID" %in% colnames(ped)) { + return(NULL) + } + + twin_rows <- which(!is.na(ped$twinID)) + + # If zygosity column exists, restrict to MZ pair_rows + if ("zygosity" %in% colnames(ped)) { + twin_rows <- twin_rows[!is.na(ped$zygosity[twin_rows]) & + ped$zygosity[twin_rows] %in% c("mz", "MZ")] + } + + if (length(twin_rows) == 0) { + return(NULL) + } + + # Build ID-to-row lookup for O(1) resolution instead of which() per pair + id_to_row <- seq_len(nrow(ped)) + names(id_to_row) <- as.character(ped$ID) + + # Use environment as hash set for O(1) membership checks + processed <- new.env(hash = TRUE, parent = emptyenv()) + + pair_rows <- vector("list", length(twin_rows)) + if (returnIDs) { + pair_ids <- vector("list", length(twin_rows)) + } + n_pairs <- 0L + + for (idx in twin_rows) { + twin_id <- ped$ID[idx] + co_twin_id <- ped$twinID[idx] + + twin_id_chr <- as.character(twin_id) + co_twin_id_chr <- as.character(co_twin_id) + + # Skip if already processed this pair (O(1) lookup) + if (exists(twin_id_chr, envir = processed, inherits = FALSE) || + exists(co_twin_id_chr, envir = processed, inherits = FALSE)) { + next + } + + # O(1) row lookup via named vector + idx1 <- id_to_row[twin_id_chr] + idx2 <- id_to_row[co_twin_id_chr] + + if (is.na(idx1) || is.na(idx2)) next + + # Always put the lower index first for consistency + if (idx1 > idx2) { + tmp <- idx1 + idx1 <- idx2 + idx2 <- tmp + } + + # O(1) insert into hash set + assign(twin_id_chr, TRUE, envir = processed) + assign(co_twin_id_chr, TRUE, envir = processed) + + n_pairs <- n_pairs + 1L + pair_rows[[n_pairs]] <- c(idx1, idx2) + if (returnIDs) { + pair_ids[[n_pairs]] <- c(twin_id, co_twin_id) + } + if (verbose) { + message( + "MZ twin pair found: ", twin_id, " (row ", idx1, + ") and ", co_twin_id, " (row ", idx2, ")" + ) + } + } + + # Trim pre-allocated lists to actual size + if (n_pairs == 0L) { + return(NULL) + } + pair_rows <- pair_rows[seq_len(n_pairs)] + if (returnIDs) { + pair_ids <- pair_ids[seq_len(n_pairs)] + } + + if (returnIDs == TRUE && returnRows == FALSE) { + if (returnAsList == TRUE) { + return(pair_ids) + } else { + data.frame( + twin1_id = vapply(pair_ids, `[`, numeric(1), 1L), + twin2_id = vapply(pair_ids, `[`, numeric(1), 2L) + ) + } + } else if (returnRows == TRUE && returnIDs == FALSE) { + if (returnAsList == TRUE) { + return(pair_rows) + } else { + data.frame( + twin1_row = vapply(pair_rows, `[`, integer(1), 1L), + twin2_row = vapply(pair_rows, `[`, integer(1), 2L) + ) + } + } else if (returnIDs == TRUE && returnRows == TRUE) { + if (returnAsList == TRUE) { + return(list(pair_rows = pair_rows, pair_ids = pair_ids)) + } else { + return(data.frame( + twin1_id = vapply(pair_ids, `[`, numeric(1), 1L), + twin2_id = vapply(pair_ids, `[`, numeric(1), 2L), + twin1_row = vapply(pair_rows, `[`, integer(1), 1L), + twin2_row = vapply(pair_rows, `[`, integer(1), 2L) + )) + } + } else { + stop("Invalid combination of returnRows and returnIDs parameters") + } +} +# replace all MZ twin IDs with the first twin's ID in each pair so they are merged for the path tracing and all subsequent steps. We will copy the values back to the second twin at the end. + +#' Fuse MZ twin pairs in a pedigree dataset for path tracing +#' This function identifies MZ twin pairs in the pedigree dataset and merges their IDs for path tracing purposes. The second twin in each pair is made a founder (with NA parents), and all children of the second twin are redirected to the first twin. This allows for correct relatedness calculations without diagonal or downstream artifacts. +#' @param ped A pedigree data.frame with columns \code{ID}, \code{momID}, \code{dadID}, and optionally \code{twinID} and \code{zygosity}. The function will look for MZ twin pairs based on the \code{twinID} column and optionally restrict to MZ pairs if a \code{zygosity} column is present. +#' @param mz_id_pairs Optional list of length-2 character vectors specifying the IDs of MZ twin pairs to fuse. If provided, this will be used instead of automatically identifying MZ twins from the \code{twinID} column. Each element should be a character vector of length 2, e.g. \code{list(c("ID1", "ID2"), c("ID3", "ID4"))}. +#' @param mz_row_pairs Optional list of length-2 integer vectors specifying the row indices of MZ twin pairs to fuse. If provided, this will be used instead of automatically identifying MZ twins from the \code{twinID} column. Each element should be an integer vector of length 2, e.g. \code{list(c(1, 2), c(3, 4))}. +#' @param test_df_twins logical. If TRUE, return the data frame of twin pairs instead of the modified pedigree. Default is FALSE. +#' @param df_twins Optional data frame with columns \code{twin1_id}, \code{twin2_id}, \code{twin1_row}, and \code{twin2_row} specifying the IDs and row indices of MZ twin pairs to fuse. If provided, this will be used instead of automatically identifying MZ twins from the \code{twinID} column. If this parameter is provided, it takes precedence over \code{mz_id_pairs} and \code{mz_row_pairs}. If \code{test_df_twins} is TRUE, this data frame will be returned for testing purposes instead of performing the fusion. +#' @param beta logical. If TRUE, use an optimized approach with O(1) lookups for large pedigrees when identifying MZ twins. Default is FALSE. +#' @param config A list of configuration options. +#' @return A modified version of the input pedigree data.frame with MZ twin pairs fused for path tracing. If \code{test_df_twins} is TRUE, returns the data frame of identified twin pairs instead. +# 70% of the time is here +fuseTwins <- function(ped, + df_twins = NULL, + mz_id_pairs = NULL, + mz_row_pairs = NULL, + config = list(verbose = FALSE), + test_df_twins = FALSE, + beta = FALSE) { + # make df_twins if not already made, and test it if requested, before proceeding with the fusion. This allows users to provide their own mz_id_pairs or mz_row_pairs and have them converted to df_twins for testing before the fusion is attempted. + + if (!is.null(df_twins)) { + if (!all(c("twin1_id", "twin2_id", "twin1_row", "twin2_row") %in% colnames(df_twins))) { + stop("df_twins must have columns twin1_id, twin2_id, twin1_row, and twin2_row") + } + if (test_df_twins == TRUE) { + return(df_twins) + } + } else if (!is.null(mz_id_pairs) && !is.null(mz_row_pairs) && length(mz_id_pairs) == length(mz_row_pairs)) { + df_twins <- lapply(1:length(mz_id_pairs), function(i) { + twin1_id <- mz_id_pairs[[i]][1] + twin2_id <- mz_id_pairs[[i]][2] + twin1_row <- mz_row_pairs[[i]][1] + twin2_row <- mz_row_pairs[[i]][2] + data.frame( + twin1_id = twin1_id, twin2_id = twin2_id, + twin1_row = twin1_row, twin2_row = twin2_row + ) + }) + df_twins <- do.call(rbind, df_twins) + rownames(df_twins) <- NULL + if (test_df_twins == TRUE) { + return(df_twins) + } + } else if (is.null(mz_id_pairs) && is.null(mz_row_pairs)) { + df_twins <- findMZtwins(ped, + verbose = config$verbose, + returnRows = TRUE, returnIDs = TRUE, returnAsList = FALSE + ) + if (test_df_twins == TRUE) { + return(df_twins) + } + } else if (!is.null(mz_id_pairs) && is.null(mz_row_pairs)) { + df_twins <- lapply(mz_id_pairs, function(pair) { + twin1_row <- which(ped$ID == pair[1]) + twin2_row <- which(ped$ID == pair[2]) + data.frame( + twin1_id = pair[1], twin2_id = pair[2], + twin1_row = twin1_row, twin2_row = twin2_row + ) + }) + df_twins <- do.call(rbind, df_twins) + rownames(df_twins) <- NULL + if (test_df_twins == TRUE) { + return(df_twins) + } + } else if (is.null(mz_id_pairs) && !is.null(mz_row_pairs)) { + df_twins <- lapply(mz_row_pairs, function(row) { + twin1_id <- ped$ID[row[1]] + twin2_id <- ped$ID[row[2]] + data.frame( + twin1_id = twin1_id, + twin2_id = twin2_id, + twin1_row = row[1], + twin2_row = row[2] + ) + }) + df_twins <- do.call(rbind, df_twins) + # remove row names + rownames(df_twins) <- NULL + if (test_df_twins == TRUE) { + return(df_twins) + } + } else { + stop("Invalid input: must provide either mz_id_pairs, mz_row_pairs, or df_twins") + } + + + fuseattemptable <- !is.null(df_twins) || (!is.null(mz_id_pairs) && length(mz_id_pairs) > 0) || (!is.null(mz_row_pairs) && length(mz_row_pairs) > 0) + if (fuseattemptable == TRUE) { + if (config$verbose == TRUE) { + message("MZ twin pairs identified for fusion") + } + + twin1s_id <- df_twins$twin1_id + twin2s_id <- df_twins$twin2_id + twin2s_row <- df_twins$twin2_row + + # Make twin2s founders + # ped$momID[twin2s_row] <- NA + # ped$dadID[twin2s_row] <- NA + + twin2s_as_mom <- which(ped$momID %in% twin2s_id) + twin2s_as_dad <- which(ped$dadID %in% twin2s_id) + # Now redirect all children of twin2 to twin1 + ped$momID[twin2s_as_mom] <- twin1s_id[match(ped$momID[twin2s_as_mom], twin2s_id)] + ped$dadID[twin2s_as_dad] <- twin1s_id[match(ped$dadID[twin2s_as_dad], twin2s_id)] + + if ("spouseID" %in% colnames(ped)) { + twin2s_as_spouse <- which(ped$spouseID %in% twin2s_id) + ped$spouseID[twin2s_as_spouse] <- twin1s_id[match(ped$spouseID[twin2s_as_spouse], twin2s_id)] + } + if ("spID" %in% colnames(ped)) { + twin2s_as_spID <- which(ped$spID %in% twin2s_id) + ped$spID[twin2s_as_spID] <- twin1s_id[match(ped$spID[twin2s_as_spID], twin2s_id)] + } + + if (config$verbose == TRUE) { + message("Merged ", length(twin1s_id), " MZ twin pair(s) in pedigree dataset for path tracing") + } + } else { + if (config$verbose == TRUE) { + message("No MZ twin pair_rows found in pedigree dataset") + } + } + + return(ped) +} diff --git a/R/segmentPedigree.R b/R/segmentPedigree.R index 3179dbfd..d74abfb3 100644 --- a/R/segmentPedigree.R +++ b/R/segmentPedigree.R @@ -10,6 +10,7 @@ #' @param momID character. Name of the column in ped for the mother ID variable #' @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 ... 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 @@ -29,17 +30,25 @@ #' ped2fam <- function(ped, personID = "ID", momID = "momID", dadID = "dadID", famID = "famID", + twinID = "twinID", ...) { # Call to wrapper function - .ped2id(ped = ped, personID = personID, momID = momID, dadID = dadID, famID = famID, type = "parents") + .ped2id( + ped = ped, personID = personID, momID = momID, dadID = dadID, famID = famID, twinID = twinID, + type = "parents" + ) } .ped2id <- function(ped, personID = "ID", momID = "momID", dadID = "dadID", - famID = "famID", type, + famID = "famID", twinID = "twinID", + type, ...) { # Turn pedigree into family - pg <- ped2graph(ped = ped, personID = personID, momID = momID, dadID = dadID, adjacent = type) + pg <- ped2graph( + ped = ped, personID = personID, momID = momID, dadID = dadID, twinID = twinID, + adjacent = type + ) # Find weakly connected components of graph wcc <- igraph::components(pg) @@ -97,6 +106,7 @@ ped2graph <- function(ped, personID = "ID", momID = "momID", dadID = "dadID", + twinID = "twinID", directed = TRUE, adjacent = c("parents", "mothers", "fathers"), ...) { @@ -196,11 +206,14 @@ ped2graph <- function(ped, #' ped2maternal <- function(ped, personID = "ID", momID = "momID", dadID = "dadID", - matID = "matID", ...) { + matID = "matID", + twinID = "twinID", + ...) { # Call to wrapper function .ped2id( ped = ped, personID = personID, momID = momID, - dadID = dadID, famID = matID, type = "mothers" + dadID = dadID, famID = matID, twinID = twinID, + type = "mothers" ) } @@ -219,10 +232,14 @@ ped2maternal <- function(ped, personID = "ID", #' ped2paternal <- function(ped, personID = "ID", momID = "momID", dadID = "dadID", - patID = "patID", ...) { + patID = "patID", + twinID = "twinID", + ...) { # Call to wrapper function .ped2id( ped = ped, personID = personID, momID = momID, - dadID = dadID, famID = patID, type = "fathers" + dadID = dadID, famID = patID, + twinID = twinID, + type = "fathers" ) } diff --git a/R/tweakPedigree.R b/R/tweakPedigree.R index 88f96128..1c8d4e6a 100644 --- a/R/tweakPedigree.R +++ b/R/tweakPedigree.R @@ -9,6 +9,7 @@ #' @param verbose logical. If TRUE, print progress through stages of algorithm #' @param gen_twin A vector of \code{generation} of the twin to be imputed. #' @param zygosity A character string indicating the zygosity of the twins. Default is "MZ" for monozygotic twins. +#' @param twin_sex A character string indicating the sex of the twins. Default is randomly assigned ("R"). If specified, it should be either "M" or "F" #' @return Returns a \code{data.frame} with MZ twins information added as a new column. #' @export @@ -18,7 +19,8 @@ makeTwins <- function(ped, ID_twin1 = NA_integer_, ID_twin2 = NA_integer_, gen_twin = 2, verbose = FALSE, - zygosity = "MZ") { + zygosity = "MZ", + twin_sex = "R") { # Check if the ped is the same format as the output of simulatePedigree if (paste0(colnames(ped), collapse = "") != paste0(c( "famID", "ID", "gen", @@ -30,37 +32,57 @@ makeTwins <- function(ped, ID_twin1 = NA_integer_, } # stop("The input pedigree is not in the same format as the output of simulatePedigree") } + ped$MZtwin <- NA_integer_ - ped$zygosity <- NA_character_ + ped$MZzygosity <- NA_character_ # Check if the two IDs are provided - if (is.na(ID_twin1) || is.na(ID_twin2)) { + + # Build ID-to-row index for O(1) attribute lookups + id_row_map <- seq_len(nrow(ped)) + names(id_row_map) <- as.character(ped$ID) + + if (is.na(ID_twin1) && is.na(ID_twin2)) { # Check if the generation is provided if (is.na(gen_twin)) { stop("You should provide either the IDs of the twins or the generation of the twins") } else { # Check if the generation is valid if (gen_twin < 2 || gen_twin > max(ped$gen)) { - stop("The generation of the twins should be an integer between 2 and the maximum generation in the pedigree") + warning("The generation of the twins should be an integer between 2 and the maximum generation in the pedigree") + # remove the MZtwin and zygosity columns + ped$MZzygosity <- NULL + ped$MZtwin <- NULL + return(ped) } else { idx <- nrow(ped[ped$gen == gen_twin & !is.na(ped$dadID), ]) usedID <- c() # randomly loop through all the individuals in the generation until find an individual who is the same sex and shares the same dadID and momID with another individual + # Pre-compute generation mask once to avoid repeated subsetting + gen_mask <- ped$gen == gen_twin & !is.na(ped$dadID) + gen_ids <- ped$ID[gen_mask] + # Build ID-to-row index for O(1) attribute lookups + id_row_map <- seq_len(nrow(ped)) + names(id_row_map) <- as.character(ped$ID) for (i in 1:idx) { - # cat("loop", i, "\n") # check if i is equal to the number of individuals in the generation usedID <- c(usedID, ID_twin1) # message(usedID) if (i < idx) { # randomly select one individual from the generation - ID_twin1 <- resample(ped$ID[ped$gen == gen_twin & !(ped$ID %in% usedID) & !is.na(ped$dadID)], 1) - # cat("twin1", ID_twin1, "\n") + ID_twin1 <- resample(ped$ID[gen_mask & !(ped$ID %in% usedID)], 1) + # Cache twin1 attributes via O(1) row lookup + twin1_row <- id_row_map[as.character(ID_twin1)] + twin1_sex <- ped$sex[twin1_row] + twin1_dad <- ped$dadID[twin1_row] + twin1_mom <- ped$momID[twin1_row] # find one same sex sibling who has the same dadID and momID as the selected individual - if (zygosity %in% c("MZ", "SS")) { - twin2_Pool <- ped$ID[ped$ID != ID_twin1 & ped$gen == gen_twin & ped$sex == ped$sex[ped$ID == ID_twin1] & ped$dadID == ped$dadID[ped$ID == ID_twin1] & ped$momID == ped$momID[ped$ID == ID_twin1]] - } else if (zygosity == "DZ") { - twin2_Pool <- ped$ID[ped$ID != ID_twin1 & ped$gen == gen_twin & ped$dadID == ped$dadID[ped$ID == ID_twin1] & ped$momID == ped$momID[ped$ID == ID_twin1]] - } else if (zygosity == "OS") { - twin2_Pool <- ped$ID[ped$ID != ID_twin1 & ped$gen == gen_twin & ped$sex != ped$sex[ped$ID == ID_twin1] & ped$dadID == ped$dadID[ped$ID == ID_twin1] & ped$momID == ped$momID[ped$ID == ID_twin1]] + sib_mask <- ped$ID != ID_twin1 & gen_mask & ped$dadID == twin1_dad & ped$momID == twin1_mom + if (zygosity %in% c("mz", "MZ", "SS", "ss")) { + twin2_Pool <- ped$ID[sib_mask & ped$sex == twin1_sex] + } else if (zygosity %in% c("DZ", "dz")) { + twin2_Pool <- ped$ID[sib_mask] + } else if (zygosity %in% c("OS", "os")) { + twin2_Pool <- ped$ID[sib_mask & ped$sex != twin1_sex] } else { stop("The zygosity should be either 'MZ', 'DZ', or 'OS'") } @@ -81,7 +103,14 @@ makeTwins <- function(ped, ID_twin1 = NA_integer_, } else { # randomly select all males or females in the generation and put them in a vector if (zygosity %in% c("MZ", "SS")) { - selectGender <- ped$ID[ped$gen == gen_twin & ped$sex == resample(c("M", "F"), 1) & !is.na(ped$dadID) & !is.na(ped$momID)] + if (twin_sex == "R") { + twin_sex_select <- resample(c("M", "F"), 1) + } else { + twin_sex_select <- twin_sex + } + + selectGender <- ped$ID[ped$gen == gen_twin & ped$sex == twin_sex_select & !is.na(ped$dadID) & !is.na(ped$momID)] + notselectGender <- ped$ID[ped$gen == gen_twin & ped$sex != twin_sex_select & !is.na(ped$dadID) & !is.na(ped$momID)] } else if (zygosity %in% c("DZ")) { selectGender <- ped$ID[ped$gen == gen_twin & !is.na(ped$dadID) & !is.na(ped$momID)] } else if (zygosity %in% c("OS")) { @@ -91,9 +120,13 @@ makeTwins <- function(ped, ID_twin1 = NA_integer_, } # message(selectGender) - if (length(selectGender) < 2) { + if (length(selectGender) < 2 && length(notselectGender) < 2) { stop("There are no available same-sex people in the generation to make twins") + } else if (twin_sex == "R" && length(selectGender) < 2 && length(notselectGender) >= 2) { + selectGender <- notselectGender } + + # randomly select two individuals from the vector ID_DoubleTwin <- resample(selectGender, 2) # message(ID_DoubleTwin) @@ -112,18 +145,84 @@ makeTwins <- function(ped, ID_twin1 = NA_integer_, # Set the zygosity of the twins } } - } else { + } else if (!is.na(ID_twin1) && !is.na(ID_twin2)) { # Impute the IDs of the twin in the MZtwin column ped$MZtwin[ped$ID == ID_twin1] <- ID_twin2 ped$MZtwin[ped$ID == ID_twin2] <- ID_twin1 + } else if (!is.na(ID_twin1) && is.na(ID_twin2)) { + + twin1_row <- id_row_map[as.character(ID_twin1)] + twin1_sex <- ped$sex[twin1_row] + twin1_dad <- ped$dadID[twin1_row] + twin1_mom <- ped$momID[twin1_row] + sib_mask <- ped$ID != ID_twin1 & !is.na(ped$dadID) & ped$dadID == twin1_dad & !is.na(ped$momID) & ped$momID == twin1_mom + if (zygosity %in% c("mz", "MZ", "SS", "ss")) { + twin2_Pool <- ped$ID[sib_mask & ped$sex == twin1_sex] + } else if (zygosity %in% c("DZ", "dz")) { + twin2_Pool <- ped$ID[sib_mask] + } else if (zygosity %in% c("OS", "os")) { + twin2_Pool <- ped$ID[sib_mask & ped$sex != twin1_sex] + } else { + stop("The zygosity should be either 'MZ', 'DZ', or 'OS'") + } + twin2_Pool <- twin2_Pool[!is.na(twin2_Pool)] + if (length(twin2_Pool) == 0) { + stop("No suitable sibling found for ID_twin1 = ", ID_twin1, " with zygosity '", zygosity, "'") + } + ID_twin2 <- resample(twin2_Pool, 1) + if (verbose) cat("Auto-selected ID_twin2 =", ID_twin2, "\n") + } else if (is.na(ID_twin1) && !is.na(ID_twin2)) { + # Mirror: find a match for twin2 + twin2_row <- id_row_map[as.character(ID_twin2)] + twin2_sex <- ped$sex[twin2_row] + twin2_dad <- ped$dadID[twin2_row] + twin2_mom <- ped$momID[twin2_row] + sib_mask <- ped$ID != ID_twin2 & !is.na(ped$dadID) & ped$dadID == twin2_dad & !is.na(ped$momID) & ped$momID == twin2_mom + if (zygosity %in% c("mz", "MZ", "SS", "ss")) { + twin1_Pool <- ped$ID[sib_mask & ped$sex == twin2_sex] + } else if (zygosity %in% c("DZ", "dz")) { + twin1_Pool <- ped$ID[sib_mask] + } else if (zygosity %in% c("OS", "os")) { + twin1_Pool <- ped$ID[sib_mask & ped$sex != twin2_sex] + } else { + stop("The zygosity should be either 'MZ', 'DZ', or 'OS'") + } + + + twin1_Pool <- twin1_Pool[!is.na(twin1_Pool)] + if (length(twin1_Pool) == 0) { + stop("No suitable sibling found for ID_twin2 = ", ID_twin2, " with zygosity '", zygosity, "'") + } + ID_twin1 <- resample(twin1_Pool, 1) + if (verbose) cat("Auto-selected ID_twin1 =", ID_twin1, "\n") } + + if (verbose == TRUE) { cat("twin1", ID_twin1, "\n") cat("twin2", ID_twin2, "\n") } - names(ped)[names(ped) == "MZtwin"] <- "twinID" - ped$zygosity[ped$ID == ID_twin1] <- zygosity - ped$zygosity[ped$ID == ID_twin2] <- zygosity + ped$MZtwin[ped$ID == ID_twin1] <- ID_twin2 + ped$MZtwin[ped$ID == ID_twin2] <- ID_twin1 + + if ("twinID" %in% colnames(ped)) { + + ped$twinID[!is.na(ped$MZtwin)] <- ped$MZtwin[!is.na(ped$MZtwin)] + ped$MZtwin <- NULL + } else { + names(ped)[names(ped) == "MZtwin"] <- "twinID" + + } + if ("zygosity" %in% colnames(ped)) { + ped$zygosity[ped$ID == ID_twin1] <- zygosity + ped$zygosity[ped$ID == ID_twin2] <- zygosity + ped$MZzygosity <- NULL + } else { + names(ped)[names(ped) == "MZzygosity"] <- "zygosity" + ped$zygosity[ped$ID == ID_twin1] <- zygosity + ped$zygosity[ped$ID == ID_twin2] <- zygosity + } + return(ped) } @@ -139,6 +238,7 @@ makeTwins <- function(ped, ID_twin1 = NA_integer_, #' @param verbose logical. If TRUE, print progress through stages of algorithm #' @param gen_inbred A vector of \code{generation} of the twin to be imputed. #' @param type_inbred A character vector indicating the type of inbreeding. "sib" for sibling inbreeding and "cousin" for cousin inbreeding. +#' @param prefer_unmated A logical indicating whether to prefer unmated siblings when automatically selecting inbred mates. Default is FALSE, which means the function will consider all siblings regardless of their mating status. #' @return Returns a \code{data.frame} with some inbred mates. #' @details #' This function creates inbred mates in the simulated pedigree \code{data.frame}. This function's purpose is to evaluate the effect of inbreeding on model fitting and parameter estimation. In case it needs to be said, we do not condone inbreeding in real life. But we recognize that it is a common practice in some fields to create inbred strains for research purposes. @@ -150,7 +250,9 @@ makeInbreeding <- function(ped, ID_mate2 = NA_integer_, verbose = FALSE, gen_inbred = 2, - type_inbred = "sib") { + type_inbred = "sib", + prefer_unmated = FALSE +) { # check if the ped is the same format as the output of simulatePedigree if (paste0(colnames(ped), @@ -175,8 +277,55 @@ makeInbreeding <- function(ped, stop("The type of inbreeding should be either 'sib' or 'cousin'") return() } - # check if the two IDs are provided - if (is.na(ID_mate1) || is.na(ID_mate2)) { + + if (!is.na(ID_mate1) && is.na(ID_mate2)) { + mate1_sex <- ped$sex[ped$ID == ID_mate1] + mate1_dad <- ped$dadID[ped$ID == ID_mate1] + mate1_mom <- ped$momID[ped$ID == ID_mate1] + if (!is.na(mate1_dad) && !is.na(mate1_mom)) { + # prefer unmated opposite-sex sibling + if (prefer_unmated == TRUE) { + pool <- makePool(ped = ped, mate_id = ID_mate1, mate_sex = mate1_sex, + mate_dad = mate1_dad, mate_mom = mate1_mom, prefer_unmated = TRUE) + + if (length(pool) == 0) { + # fall back to mated opposite-sex sibling + pool <- makePool(ped = ped, mate_id = ID_mate1, mate_sex = mate1_sex, + mate_dad = mate1_dad, mate_mom = mate1_mom, prefer_unmated = FALSE) + } + } else { + pool <- makePool(ped = ped, mate_id = ID_mate1, mate_sex = mate1_sex, + mate_dad = mate1_dad, mate_mom = mate1_mom, prefer_unmated = FALSE) + } + + if (length(pool) > 0) { + ID_mate2 <- resample(pool, 1) + if (verbose) message("Auto-selected ID_mate2 = ", ID_mate2) + } + } + } else if (is.na(ID_mate1) && !is.na(ID_mate2)) { + mate2_sex <- ped$sex[ped$ID == ID_mate2] + mate2_dad <- ped$dadID[ped$ID == ID_mate2] + mate2_mom <- ped$momID[ped$ID == ID_mate2] + if (!is.na(mate2_dad) && !is.na(mate2_mom)) { + if (prefer_unmated) { + pool <- makePool(ped = ped, mate_id = ID_mate2, mate_sex = mate2_sex, + mate_dad = mate2_dad, mate_mom = mate2_mom, prefer_unmated = TRUE) + # fall back + if (length(pool) == 0) { + pool <- makePool(ped = ped, mate_id = ID_mate2, mate_sex = mate2_sex, + mate_dad = mate2_dad, mate_mom = mate2_mom, prefer_unmated = FALSE) + } + } else { + pool <- makePool(ped = ped, mate_id = ID_mate2, mate_sex = mate2_sex, + mate_dad = mate2_dad, mate_mom = mate2_mom, prefer_unmated = FALSE) + } + if (length(pool) > 0) { + ID_mate1 <- resample(pool, 1) + if (verbose) message("Auto-selected ID_mate1 = ", ID_mate1) + } + } # check if the two IDs are provided + } else if (is.na(ID_mate1) && is.na(ID_mate2)) { # Check if the generation is provided if (is.na(gen_inbred)) { stop("You should provide either the IDs of the inbred mates or the generation of the inbred mates") @@ -195,15 +344,35 @@ makeInbreeding <- function(ped, if (!is.na(ID_mate2)) { break } - ID_pool_mate1 <- ped$ID[ped$gen == gen_inbred & !is.na(ped$dadID) & !is.na(ped$momID) & is.na(ped$spID) & !(ped$ID %in% usedID)] - # if the pool is empty, find all individuals who have the same dadID and momID as the selected individual but mated - if (length(ID_pool_mate1) == 0) { + + if (prefer_unmated) { + # try to find one opposite + ID_pool_mate1 <- ped$ID[ped$gen == gen_inbred & !is.na(ped$dadID) & !is.na(ped$momID) & is.na(ped$spID) & !(ped$ID %in% usedID)] + + + + # if the pool is empty, find all individuals who have the same dadID and momID as the selected individual but mated + if (length(ID_pool_mate1) == 0) { + ID_pool_mate1 <- ped$ID[ped$gen == gen_inbred & !is.na(ped$dadID) & !is.na(ped$momID) & !(ped$ID %in% usedID)] + } + } else { ID_pool_mate1 <- ped$ID[ped$gen == gen_inbred & !is.na(ped$dadID) & !is.na(ped$momID) & !(ped$ID %in% usedID)] } + ID_mate1 <- resample(ID_pool_mate1, 1) usedID <- c(usedID, ID_mate1) # try to find one opposite-sex individual who has the same dadID and momID as the selected individual, preferalbly not mated - ID_pool_mate2 <- ped$ID[ped$gen == gen_inbred & ped$sex != ped$sex[ped$ID == ID_mate1] & ped$dadID == ped$dadID[ped$ID == ID_mate1] & ped$momID == ped$momID[ped$ID == ID_mate1] & is.na(ped$spID)] + if (prefer_unmated) { + ID_pool_mate2 <- ped$ID[ped$gen == gen_inbred & ped$sex != ped$sex[ped$ID == ID_mate1] & ped$dadID == ped$dadID[ped$ID == ID_mate1] & ped$momID == ped$momID[ped$ID == ID_mate1] & is.na(ped$spID)] + # back up + if (length(ID_pool_mate1) == 0) { + ID_pool_mate2 <- ped$ID[ped$gen == gen_inbred & ped$sex != ped$sex[ped$ID == ID_mate1] & ped$dadID == ped$dadID[ped$ID == ID_mate1] & ped$momID == ped$momID[ped$ID == ID_mate1]] + } + } else { + # ID_pool_mate2 <- ped$ID[ped$gen == gen_inbred & ped$sex != ped$sex[ped$ID == ID_mate1] & ped$dadID == ped$dadID[ped$ID == ID_mate1] & ped$momID == ped$momID[ped$ID == ID_mate1]] + ID_pool_mate2 <- makePool(ped = ped, mate_id = ID_mate1, mate_sex = ped$sex[ped$ID == ID_mate1], mate_dad = ped$dadID[ped$ID == ID_mate1], mate_mom = ped$momID[ped$ID == ID_mate1], prefer_unmated = prefer_unmated, gen_inbred = gen_inbred) + } + # if the pool is not empty, randomly select one individual from the pool if (length(ID_pool_mate2) > 0) { ID_mate2 <- resample(ID_pool_mate2, 1) @@ -237,20 +406,18 @@ makeInbreeding <- function(ped, # change the spouseID of ID_mate1 and ID_mate2 to each other ped$spID[ped$ID == ID_mate1] <- ID_mate2 ped$spID[ped$ID == ID_mate2] <- ID_mate1 - # change the individuals in next generation whose dadID and momID are ID_mate1 and ID_mate2's former mates to ID_mate1 and ID_mate2 - for (j in seq_len(nrow(ped))) { - if (!is.na(ped$dadID[j]) & !is.na(ID_mate1_former_mate) & ped$dadID[j] == ID_mate1_former_mate) { - ped$dadID[j] <- ID_mate2 - } - if (!is.na(ped$momID[j]) & !is.na(ID_mate1_former_mate) & ped$momID[j] == ID_mate1_former_mate) { - ped$momID[j] <- ID_mate2 - } - if (!is.na(ped$dadID[j]) & !is.na(ID_mate2_former_mate) & ped$dadID[j] == ID_mate2_former_mate) { - ped$dadID[j] <- ID_mate1 - } - if (!is.na(ped$momID[j]) & !is.na(ID_mate2_former_mate) & ped$momID[j] == ID_mate2_former_mate) { - ped$momID[j] <- ID_mate1 - } + # Vectorized replacement of former mates' IDs in children's parent columns + if (!is.na(ID_mate1_former_mate)) { + dad_match1 <- which(!is.na(ped$dadID) & ped$dadID == ID_mate1_former_mate) + if (length(dad_match1) > 0) ped$dadID[dad_match1] <- ID_mate2 + mom_match1 <- which(!is.na(ped$momID) & ped$momID == ID_mate1_former_mate) + if (length(mom_match1) > 0) ped$momID[mom_match1] <- ID_mate2 + } + if (!is.na(ID_mate2_former_mate)) { + dad_match2 <- which(!is.na(ped$dadID) & ped$dadID == ID_mate2_former_mate) + if (length(dad_match2) > 0) ped$dadID[dad_match2] <- ID_mate1 + mom_match2 <- which(!is.na(ped$momID) & ped$momID == ID_mate2_former_mate) + if (length(mom_match2) > 0) ped$momID[mom_match2] <- ID_mate1 } return(ped) } @@ -265,13 +432,28 @@ makeInbreeding <- function(ped, #' @param gen_drop the generation in which the randomly dropped person is. Will work if `ID_drop` is not specified. #' @param sex_drop the biological sex of the randomly dropped person. #' @param n_drop the number of times the mutation happens. +#' @param verbose logical. If TRUE, print progress through stages of algorithm #' @return a pedigree with the dropped person's `dadID` and `momID` set to NA. #' @export dropLink <- function(ped, ID_drop = NA_integer_, gen_drop = 2, sex_drop = NA_character_, - n_drop = 1) { + n_drop = 1, + verbose = FALSE) { + # Standardize column names for consistency with other functions + + if (paste0(colnames(ped), + collapse = "" + ) != paste0( + c("famID", "ID", "gen", "dadID", "momID", "spID", "sex"), + collapse = "" + )) { + ped <- standardizeColnames(ped, verbose = verbose) + if (verbose == TRUE) { + message("The input pedigree is not in the same format as the output of simulatePedigree\n") + } + } # check if the ID_drop is specified if (is.na(ID_drop)) { # check if the sex_drop is specified @@ -406,3 +588,29 @@ addPersonToPed <- function(ped, name = NULL, return(rbind(ped, new_row)) } } + + +makePool <- function(ped, mate_id, mate_sex, mate_dad, mate_mom, prefer_unmated = FALSE, + gen_inbred = NULL, usedID = NULL +) { + if (!is.null(gen_inbred)) { + ped <- ped[ped$gen == gen_inbred, ] + } + if (!is.null(usedID)) { + ped <- ped[!(ped$ID %in% usedID), ] + } + + + # should we prefer unmated siblings when automatically selecting inbred mates? If yes, we will only consider unmated siblings. If no, we will consider all siblings regardless of their mating status. + if (prefer_unmated == TRUE) { + mated_mask <- !is.na(ped$spID) + + } else { + mated_mask <- rep(TRUE, nrow(ped)) + } + + pool <- ped$ID[ped$ID != mate_id & ped$sex != mate_sex & + !is.na(ped$dadID) & ped$dadID == mate_dad & + !is.na(ped$momID) & ped$momID == mate_mom & mated_mask] + return(pool) +} diff --git a/data-raw/optimizing.R b/data-raw/optimizing.R index c279bf79..1a11914f 100644 --- a/data-raw/optimizing.R +++ b/data-raw/optimizing.R @@ -3,19 +3,11 @@ library(microbenchmark) library(tidyverse) set.seed(1164127) Ngen <- 3 -kpc <- 6 +kpc <- 4 sexR <- .50 # sometimes fails above .5 -marR <- .9 +marR <- .8 reps <- 15 -if (FALSE) { - profvis({ - simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = beta_F) - }) - profvis({ - simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = beta_T) - }) -} # mz_method_opts <- c("addtwins", "merging") beta_method_opts <- c(TRUE, FALSE) beta_F <- T @@ -29,114 +21,189 @@ df_gen1 <- simulatePedigree( ) %>% makeTwins(gen_twin = 1) -df_lowgen <- simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = F) %>% +df_lowgen <- simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = TRUE) %>% + makeTwins(gen_twin = gen_twin) + +df_midgen <- simulatePedigree(kpc = kpc, Ngen = Ngen * 2, sexR = sexR, marR = marR, beta = TRUE) %>% + makeTwins(gen_twin = gen_twin) + +df_highgen <- simulatePedigree(kpc = kpc, Ngen = Ngen * 3, sexR = sexR, marR = marR, beta = TRUE) %>% + makeTwins(gen_twin = gen_twin) + +# r_mz1 <- df_midgen %>% +# ped2com(mz_method = "merging", mz_twins = TRUE) +# r_mz2 <- df_midgen %>% +# ped2com(mz_method = "addtwins", mz_twins = TRUE) +# expect_equal(length(r_mz1@i), length(r_mz2@i)) +# expect_equal(length(r_mz1@x), length(r_mz2@x)) +# expect_equal(length(r_mz1@p), length(r_mz2@p)) + +df_gen1 <- simulatePedigree( + kpc = kpc, Ngen = 1, sexR = sexR, marR = marR, + beta = TRUE +) %>% + makeTwins(gen_twin = 1) + +df_lowgen <- simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = T) %>% makeTwins(gen_twin = gen_twin) -df_midgen <- simulatePedigree(kpc = kpc, Ngen = Ngen * 2, sexR = sexR, marR = marR, beta = F) %>% +df_midgen <- simulatePedigree(kpc = kpc, Ngen = Ngen * 2, sexR = sexR, marR = marR, beta = T) %>% makeTwins(gen_twin = gen_twin) -df_highgen <- simulatePedigree(kpc = kpc, Ngen = Ngen * 2 + 1, sexR = sexR, marR = marR, beta = F) %>% +df_highgen <- simulatePedigree(kpc = kpc, Ngen = Ngen * 3, sexR = sexR, marR = marR, beta = T) %>% makeTwins(gen_twin = gen_twin) -r_mz1 <- df_midgen %>% - ped2add(mz_method = "merging", mz_twins = TRUE) -r_mz2 <- df_midgen %>% - ped2add(mz_method = "addtwins", mz_twins = TRUE) +#r_mz1 <- df_midgen %>% +# ped2com(mz_method = "merging", mz_twins = TRUE) +#r_mz2 <- df_midgen %>% +# ped2com(mz_method = "addtwins", mz_twins = TRUE) # expect_equal(length(r_mz1@i), length(r_mz2@i)) # expect_equal(length(r_mz1@x), length(r_mz2@x)) # expect_equal(length(r_mz1@p), length(r_mz2@p)) benchmark_results <- microbenchmark( - beta_null_1gen = { + null_1gen = { df_gen1 %>% - ped2add(mz_twins = F) + ped2com(mz_twins = F,component = "additive") }, - beta_false_1gen = { + + addtwins_1gen = { df_gen1 %>% - ped2add(mz_method = "addtwins", mz_twins = TRUE) + ped2com(mz_method = "addtwins", mz_twins = TRUE,component = "additive") }, - beta_true_1gen = { + merging_1gen = { df_gen1 %>% - ped2add(mz_method = "merging", mz_twins = TRUE) + ped2com(mz_method = "merging", mz_twins = TRUE,component = "additive") }, - beta_null_lowgen = { + null_beta_1gen = { + df_gen1 %>% + ped2com(mz_twins = F,component = "additive", beta = TRUE) + }, + addtwins_beta_1gen = { + df_gen1 %>% + ped2com(mz_method = "addtwins", mz_twins = TRUE,component = "additive",beta=TRUE) + }, + merging_beta_1gen = { + df_gen1 %>% + ped2com(mz_method = "merging", mz_twins = TRUE,component = "additive",beta=TRUE) + }, + null_lowgen = { + df_lowgen %>% + ped2com(mz_twins = F,component = "additive") + }, + addtwins_lowgen = { + df_lowgen %>% + ped2com(mz_method = "addtwins", mz_twins = TRUE,component = "additive") + }, + merging_lowgen = { + df_lowgen %>% + ped2com(mz_method = "merging", mz_twins = TRUE,component = "additive") + }, + null_beta_lowgen = { df_lowgen %>% - ped2add(mz_twins = F) + ped2com(mz_twins = F,component = "additive", beta=TRUE) }, - beta_false_lowgen = { + addtwins_beta_lowgen = { df_lowgen %>% - ped2add(mz_method = "addtwins", mz_twins = TRUE) + ped2com(mz_method = "addtwins", mz_twins = TRUE,component = "additive", beta=TRUE) }, - beta_true_lowgen = { + merging_beta_lowgen = { df_lowgen %>% - ped2add(mz_method = "merging", mz_twins = TRUE) + ped2com(mz_method = "merging", mz_twins = TRUE,component = "additive", beta=TRUE) }, - beta_null_midgen = { + null_midgen = { df_midgen %>% - ped2add(mz_twins = F) + ped2com(mz_twins = F,component = "additive") }, - beta_false_midgen = { + addtwins_midgen = { df_midgen %>% - ped2add(mz_method = "addtwins", mz_twins = TRUE) + ped2com(mz_method = "addtwins", mz_twins = TRUE,component = "additive") }, - beta_true_midgen = { + merging_midgen = { df_midgen %>% - ped2add(mz_method = "merging", mz_twins = TRUE) + ped2com(mz_method = "merging", mz_twins = TRUE,component = "additive") }, - beta_null_highgen = { - df_highgen %>% - ped2add(mz_twins = F) + null_beta_midgen = { + df_midgen %>% + ped2com(mz_twins = F,component = "additive", beta=TRUE) }, - beta_false_highgen = { - df_highgen %>% - ped2add(mz_method = "addtwins", mz_twins = TRUE) + addtwins_beta_midgen = { + df_midgen %>% + ped2com(mz_method = "addtwins", mz_twins = TRUE,component = "additive", beta=TRUE) }, - beta_true_highgen = { - df_highgen %>% - ped2add(mz_method = "merging", mz_twins = TRUE) + merging_beta_midgen = { + df_midgen %>% + ped2com(mz_method = "merging", mz_twins = TRUE,component = "additive", beta=TRUE) }, + # null_highgen = { + # df_highgen %>% + # ped2com(mz_twins = F,component = "additive") + # }, + # addtwins_highgen = { + # df_highgen %>% + # ped2com(mz_method = "addtwins", mz_twins = TRUE,component = "additive") + # }, + # merging_highgen = { + # df_highgen %>% + # ped2com(mz_method = "merging", mz_twins = TRUE,component = "additive") + # }, + # null_beta_highgen = { + # df_highgen %>% + # ped2com(mz_twins = F,component = "additive", beta=TRUE) + # }, + # addtwins_beta_highgen = { + # df_highgen %>% + # ped2com(mz_method = "addtwins", mz_twins = TRUE,component = "additive", beta=TRUE) + # }, + # merging_beta_highgen = { + # df_highgen %>% + # ped2com(mz_method = "merging", mz_twins = TRUE,component = "additive", beta=TRUE) + # }, times = reps # Run each method 10 times ) - - benchmark_results <- benchmark_results %>% mutate( - beta_factor = factor(case_when( - grepl("beta_true", expr) ~ "TRUE", - grepl("beta_false", expr) ~ "FALSE", - grepl("beta_null", expr) ~ "NULL", - grepl("beta_indexed", expr) ~ "indexed" + twin_factor = factor(case_when( + grepl("merging", expr) ~ "merging", + grepl("addtwins", expr) ~ "addtwins", + grepl("null", expr) ~ "NULL" )), - beta = ifelse(grepl("beta_false", expr), FALSE, TRUE), + twin_factor = factor(twin_factor, levels = c("NULL", "addtwins", "merging")), + beta = ifelse(grepl("beta_", expr), TRUE,FALSE), gen_num = case_when( grepl("1gen", expr) ~ 1, grepl("lowgen", expr) ~ Ngen, grepl("midgen", expr) ~ Ngen * 2, - grepl("highgen", expr) ~ Ngen * 2 + 1 + grepl("highgen", expr) ~ Ngen * 3 ), gen_factor = factor(gen_num, levels = c(1, Ngen, Ngen * 2, Ngen * 2 + 1)) ) summary(benchmark_results) -lm(benchmark_results$time ~ benchmark_results$beta_factor * benchmark_results$gen_num) %>% +lm(benchmark_results$time ~ benchmark_results$beta * benchmark_results$twin_factor * benchmark_results$gen_num) %>% + summary() +lm(benchmark_results$time ~ benchmark_results$beta * benchmark_results$gen_num) %>% summary() -lm(benchmark_results$time ~ benchmark_results$beta_factor) %>% +lm(benchmark_results$time ~ benchmark_results$twin_factor * benchmark_results$beta) %>% summary() # log transform time for better visualization -ggplot(benchmark_results, aes(x = gen_factor, y = time / 1e6, color = beta_factor)) + - geom_boxplot() + +ggplot(benchmark_results, aes(x = gen_factor, y = time / 1e6, color = twin_factor, fill = beta)) + + geom_boxplot(notch = TRUE, position = position_dodge(width = 0.8)) + labs( title = "Benchmarking simulatePedigree() with and without beta parameter", x = "Generation Size", y = "Execution Time (ms)", - color = "Beta Parameter" + color = "Twin Method" ) + theme_minimal() + - scale_y_log10() + scale_y_log10() + + scale_fill_manual(values = c("TRUE" = "black", "FALSE" = "white"), name = "Beta Parameter") + + scale_color_manual(values = c("NULL" = "gray", "addtwins" = "skyblue3", + "merging" = "tomato2"), name = "Twin Method") @@ -149,6 +216,7 @@ kpc <- 4 sexR <- .50 # sometimes fails above .5 marR <- .7 reps <- 10 + if (FALSE) { profvis({ simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = FALSE) @@ -158,60 +226,93 @@ if (FALSE) { simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = TRUE) }) } -if (FALSE) { - benchmark_results <- microbenchmark( - beta_false_1gen = { - simulatePedigree(kpc = kpc, Ngen = 1, sexR = sexR, marR = marR, beta = FALSE) - }, - beta_true_1gen = { - simulatePedigree(kpc = kpc, Ngen = 1, sexR = sexR, marR = marR, beta = TRUE) - }, - beta_false_lowgen = { - simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = FALSE) - }, - beta_true_lowgen = { - simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = TRUE) - }, - beta_false_midgen = { - simulatePedigree(kpc = kpc, Ngen = Ngen * 2, sexR = sexR, marR = marR, beta = FALSE) - }, - beta_true_midgen = { - simulatePedigree(kpc = kpc, Ngen = Ngen * 2, sexR = sexR, marR = marR, beta = TRUE) - }, - beta_false_highgen = { - simulatePedigree(kpc = kpc, Ngen = Ngen * 3, sexR = sexR, marR = marR, beta = FALSE) - }, - beta_true_highgen = { - simulatePedigree(kpc = kpc, Ngen = Ngen * 3, sexR = sexR, marR = marR, beta = TRUE) - }, - times = reps # Run each method 10 times - ) - - benchmark_results <- benchmark_results %>% - mutate( - beta_factor = factor(case_when( - grepl("beta_true", expr) ~ "TRUE", - grepl("beta_false", expr) ~ "FALSE", - grepl("beta_indexed", expr) ~ "indexed" - )), - beta = ifelse(grepl("beta_false", expr), FALSE, TRUE), - gen_num = case_when( - grepl("1gen", expr) ~ 1, - grepl("lowgen", expr) ~ Ngen, - grepl("midgen", expr) ~ Ngen * 2, - grepl("highgen", expr) ~ Ngen * 3 - ), - gen_factor = factor(gen_num, levels = c(1, Ngen, Ngen * 2, Ngen * 3)) + if (FALSE) { + benchmark_results <- microbenchmark( + beta_false_1gen = { + simulatePedigree(kpc = kpc, Ngen = 1, sexR = sexR, marR = marR, beta = FALSE) + }, + beta_true_1gen = { + simulatePedigree(kpc = kpc, Ngen = 1, sexR = sexR, marR = marR, beta = TRUE) + }, + beta_false_lowgen = { + simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = FALSE) + }, + beta_true_lowgen = { + simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = TRUE) + }, + beta_false_midgen = { + simulatePedigree(kpc = kpc, Ngen = Ngen * 2, sexR = sexR, marR = marR, beta = FALSE) + }, + beta_true_midgen = { + simulatePedigree(kpc = kpc, Ngen = Ngen * 2, sexR = sexR, marR = marR, beta = TRUE) + }, + beta_false_highgen = { + simulatePedigree(kpc = kpc, Ngen = Ngen * 3, sexR = sexR, marR = marR, beta = FALSE) + }, + beta_true_highgen = { + simulatePedigree(kpc = kpc, Ngen = Ngen * 3, sexR = sexR, marR = marR, beta = TRUE) + }, + times = reps # Run each method 10 times ) - summary(benchmark_results) - lm(benchmark_results$time ~ benchmark_results$beta * benchmark_results$gen_num) %>% - summary() + benchmark_results <- benchmark_results %>% + mutate( + beta_factor = factor(case_when( + grepl("beta_true", expr) ~ "TRUE", + grepl("beta_false", expr) ~ "FALSE", + grepl("beta_indexed", expr) ~ "indexed" + )), + beta = ifelse(grepl("beta_false", expr), FALSE, TRUE), + gen_num = case_when( + grepl("1gen", expr) ~ 1, + grepl("lowgen", expr) ~ Ngen, + grepl("midgen", expr) ~ Ngen * 2, + grepl("highgen", expr) ~ Ngen * 3 + ), + gen_factor = factor(gen_num, levels = c(1, Ngen, Ngen * 2, Ngen * 3)) + ) + benchmark_results <- benchmark_results %>% + mutate( + beta_factor = factor(case_when( + grepl("beta_true", expr) ~ "TRUE", + grepl("beta_false", expr) ~ "FALSE", + grepl("beta_indexed", expr) ~ "indexed" + )), + beta = ifelse(grepl("beta_false", expr), FALSE, TRUE), + gen_num = case_when( + grepl("1gen", expr) ~ 1, + grepl("lowgen", expr) ~ Ngen, + grepl("midgen", expr) ~ Ngen * 2, + grepl("highgen", expr) ~ Ngen * 3 + ), + gen_factor = factor(gen_num, levels = c(1, Ngen, Ngen * 2, Ngen * 3)) + ) + + summary(benchmark_results) + lm(benchmark_results$time ~ benchmark_results$beta * benchmark_results$gen_num) %>% + summary() + summary(benchmark_results) + lm(benchmark_results$time ~ benchmark_results$beta * benchmark_results$gen_num) %>% + summary() - lm(benchmark_results$time ~ benchmark_results$beta) %>% - summary() - # log transform time for better visualization + lm(benchmark_results$time ~ benchmark_results$beta) %>% + summary() + # log transform time for better visualization + lm(benchmark_results$time ~ benchmark_results$beta) %>% + summary() + # log transform time for better visualization + ggplot(benchmark_results, aes(x = gen_factor, y = time / 1e6, color = beta_factor)) + + geom_boxplot() + + labs( + title = "Benchmarking simulatePedigree() with and without beta parameter", + x = "Generation Size", + y = "Execution Time (ms)", + color = "Beta Parameter" + ) + + theme_minimal() + + scale_y_log10() + } ggplot(benchmark_results, aes(x = gen_factor, y = time / 1e6, color = beta_factor)) + geom_boxplot() + labs( diff --git a/data-raw/optimizing_simulations.R b/data-raw/optimizing_simulations.R new file mode 100644 index 00000000..e824760d --- /dev/null +++ b/data-raw/optimizing_simulations.R @@ -0,0 +1,117 @@ +library(profvis) +library(microbenchmark) +library(tidyverse) +set.seed(1667) +Ngen <- 3 +kpc <- 4 +sexR <- .50 # sometimes fails above .5 +marR <- .7 +reps <- 10 + +if (FALSE) { + profvis({ + simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = FALSE) + }) + + profvis({ + simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = TRUE) + }) +} + if (FALSE) { + benchmark_results <- microbenchmark( + beta_false_1gen = { + simulatePedigree(kpc = kpc, Ngen = 1, sexR = sexR, marR = marR, beta = FALSE) + }, + beta_true_1gen = { + simulatePedigree(kpc = kpc, Ngen = 1, sexR = sexR, marR = marR, beta = TRUE) + }, + beta_false_lowgen = { + simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = FALSE) + }, + beta_true_lowgen = { + simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = TRUE) + }, + beta_false_midgen = { + simulatePedigree(kpc = kpc, Ngen = Ngen * 2, sexR = sexR, marR = marR, beta = FALSE) + }, + beta_true_midgen = { + simulatePedigree(kpc = kpc, Ngen = Ngen * 2, sexR = sexR, marR = marR, beta = TRUE) + }, + beta_false_highgen = { + simulatePedigree(kpc = kpc, Ngen = Ngen * 3, sexR = sexR, marR = marR, beta = FALSE) + }, + beta_true_highgen = { + simulatePedigree(kpc = kpc, Ngen = Ngen * 3, sexR = sexR, marR = marR, beta = TRUE) + }, + times = reps # Run each method 10 times + ) + + benchmark_results <- benchmark_results %>% + mutate( + beta_factor = factor(case_when( + grepl("beta_true", expr) ~ "TRUE", + grepl("beta_false", expr) ~ "FALSE", + grepl("beta_indexed", expr) ~ "indexed" + )), + beta = ifelse(grepl("beta_false", expr), FALSE, TRUE), + gen_num = case_when( + grepl("1gen", expr) ~ 1, + grepl("lowgen", expr) ~ Ngen, + grepl("midgen", expr) ~ Ngen * 2, + grepl("highgen", expr) ~ Ngen * 3 + ), + gen_factor = factor(gen_num, levels = c(1, Ngen, Ngen * 2, Ngen * 3)) + ) + benchmark_results <- benchmark_results %>% + mutate( + beta_factor = factor(case_when( + grepl("beta_true", expr) ~ "TRUE", + grepl("beta_false", expr) ~ "FALSE", + grepl("beta_indexed", expr) ~ "indexed" + )), + beta = ifelse(grepl("beta_false", expr), FALSE, TRUE), + gen_num = case_when( + grepl("1gen", expr) ~ 1, + grepl("lowgen", expr) ~ Ngen, + grepl("midgen", expr) ~ Ngen * 2, + grepl("highgen", expr) ~ Ngen * 3 + ), + gen_factor = factor(gen_num, levels = c(1, Ngen, Ngen * 2, Ngen * 3)) + ) + + summary(benchmark_results) + lm(benchmark_results$time ~ benchmark_results$beta * benchmark_results$gen_num) %>% + summary() + summary(benchmark_results) + lm(benchmark_results$time ~ benchmark_results$beta * benchmark_results$gen_num) %>% + summary() + + lm(benchmark_results$time ~ benchmark_results$beta) %>% + summary() + # log transform time for better visualization + lm(benchmark_results$time ~ benchmark_results$beta) %>% + summary() + # log transform time for better visualization + + ggplot(benchmark_results, aes(x = gen_factor, y = time / 1e6, color = beta_factor)) + + geom_boxplot() + + labs( + title = "Benchmarking simulatePedigree() with and without beta parameter", + x = "Generation Size", + y = "Execution Time (ms)", + color = "Beta Parameter" + ) + + theme_minimal() + + scale_y_log10() + } + ggplot(benchmark_results, aes(x = gen_factor, y = time / 1e6, color = beta_factor)) + + geom_boxplot() + + labs( + title = "Benchmarking simulatePedigree() with and without beta parameter", + x = "Generation Size", + y = "Execution Time (ms)", + color = "Beta Parameter" + ) + + theme_minimal() + + scale_y_log10() +} diff --git a/data-raw/optimizing_twins.R b/data-raw/optimizing_twins.R new file mode 100644 index 00000000..511f3207 --- /dev/null +++ b/data-raw/optimizing_twins.R @@ -0,0 +1,379 @@ +library(profvis) +library(microbenchmark) +library(tidyverse) +devtools::load_all(".") + + +# --------------------------- +# 0) Config +# --------------------------- +cfg <- list( + seed = 1164127, + Ngen_base = 3, + reps = 10, + all_scenarios = FALSE, # set to TRUE to run all scenarios defined below + include_highgen = TRUE, + include_1gen = FALSE, + include_lowgen = FALSE +) +cfg$gen_twin <- ceiling(cfg$Ngen_base - 1) + +set.seed(cfg$seed) + +# --------------------------- +# 1) Levels (edit here to extend) +# --------------------------- +levels <- list( + ped = tibble( + ped_label = c( + if (cfg$include_1gen) "1gen", + if (cfg$include_lowgen) "lowgen", + "midgen", + if (cfg$include_highgen) "highgen" + ), + Ngen_total = c( + if (cfg$include_1gen) 1, + if (cfg$include_lowgen) cfg$Ngen_base, + cfg$Ngen_base * 2, + if (cfg$include_highgen) cfg$Ngen_base * 3 + ), + gen_twin = c( + if (cfg$include_1gen) 1, + if (cfg$include_lowgen) cfg$gen_twin, + cfg$gen_twin, + if (cfg$include_highgen) cfg$gen_twin + ) + # Add highgen row whenever you want + ), + + # Simulation-side factors (simulatePedigree) + kpc = 3, # set to c(2, 3, 4) to vary + sexR = 0.50, # sometimes fails above .5 + marR = c(0.8), # set to c(0.6, 0.8, 0.9) to vary + sim_beta = TRUE, # set to c(TRUE, FALSE) if you ever want to vary + + # Conversion-side factors (ped2com) + component = c("additive"), + twin_method = c("NULL", "addtwins", "merging"), + beta = c(FALSE, TRUE), + sparse_matrix = c(FALSE, TRUE) # user-facing name, translated to ped2com's `sparse` +) + +# Which columns define a unique simulation vs conversion condition +# If you add a new factor, put its name in the right scope here. + +scopes <- list( + sim = c("ped_label", "Ngen_total", "gen_twin", "kpc", "sexR", "marR", "sim_beta"), + conv = c("component", "twin_method", "beta", "sparse_matrix") +) + +# --------------------------- +# 2) Scenarios (edit here to control crossing) +# --------------------------- +# Each scenario says what to vary and what to fix. No other code changes. +scenarios <- list( + full = list( + vary = c("ped", "marR", "twin_method", "beta", "sparse_matrix"), + fixed = list() + ), + quick = list( + vary = c("ped", "twin_method"), + fixed = list(marR = 0.8, beta = FALSE, sparse_matrix = TRUE, component = "additive") + ) + + # Add more scenarios whenever you want: + # e.g., "marR_only" = list(vary=c("marR"), fixed=list(ped = levels$ped[levels$ped$ped_label=="midgen",], ...)) +) + +if (!cfg$all_scenarios) { + scenarios <- scenarios[c("full")] # order control; also, comment out any you don't want to run +} +# --------------------------- +# 3) Generic design builder (do not edit to add factors) +# --------------------------- +`%||%` <- function(x, y) if (!is.null(x)) x else y + +level_tbl <- function(name, x) { + if (inherits(x, "data.frame")) { + return(x) + } + tibble(!!name := x) +} + +default_value <- function(x) { + if (inherits(x, "data.frame")) { + return(x[1, , drop = FALSE]) + } + x[[1]] +} + +expand_scenario <- function(levels, vary, fixed = list(), scenario = "scenario") { + df <- tibble(.dummy = 1) + + # expand varied factors + for (nm in vary) { + df <- tidyr::crossing(df, level_tbl(nm, levels[[nm]])) + } + + # add fixed/default factors not varied + not_varied <- setdiff(names(levels), vary) + for (nm in not_varied) { + lv <- levels[[nm]] + + if (inherits(lv, "data.frame")) { + const <- fixed[[nm]] %||% default_value(lv) + if (!inherits(const, "data.frame") || nrow(const) != 1) const <- default_value(lv) + df <- bind_cols(df, const[rep(1, nrow(df)), , drop = FALSE]) + } else { + df[[nm]] <- fixed[[nm]] %||% default_value(lv) + } + } + + df %>% + select(-.dummy) %>% + mutate(scenario = scenario, .before = 1) +} + +build_design <- function(levels, scenarios) { + purrr::imap_dfr( + scenarios, + ~ expand_scenario(levels, vary = .x$vary, fixed = .x$fixed, scenario = .y) + ) +} + +design <- build_design(levels, scenarios) + +# Create a stable, joinable label (also becomes the microbenchmark expr name) +design <- design %>% + mutate( + label = paste0( + scenario, + "|ped=", ped_label, + "|marR=", marR, + "|twin=", twin_method, + "|beta=", beta, + "|sparse=", sparse_matrix, + "|comp=", component + ) + ) + +# --------------------------- +# 4) Simulation cache (simulate once per unique sim condition) +# --------------------------- +simulate_one <- function(Ngen_total, gen_twin, kpc, sexR, marR, sim_beta) { + simulatePedigree( + kpc = kpc, Ngen = Ngen_total, sexR = sexR, marR = marR, + beta = sim_beta + ) %>% + makeTwins(gen_twin = gen_twin) +} + +sim_tbl <- design %>% + distinct(across(all_of(scopes$sim))) %>% + mutate( + sim_id = row_number(), + ped = pmap( + list(Ngen_total, gen_twin, kpc, sexR, marR, sim_beta, sim_id), + function(Ngen_total, gen_twin, kpc, sexR, marR, sim_beta, sim_id) { + set.seed(cfg$seed + sim_id) + simulate_one(Ngen_total, gen_twin, kpc, sexR, marR, sim_beta) + } + ) + ) %>% + select(-sim_id) + +design <- design %>% left_join(sim_tbl, by = scopes$sim) + +# Put peds and args in keyed lists so benchmark expressions stay tiny +peds_by_label <- setNames(design$ped, design$label) + +# --------------------------- +# 5) Conversion arg translation (edit only if you add non-1:1 args) +# --------------------------- +# Columns that map directly to ped2com arg names go through automatically. +# Anything else gets a translation rule here. +special_to_args <- list( + twin_method = function(v) { + if (is.null(v) || length(v) == 0 || is.na(v) || v == "NULL") { + list(mz_twins = FALSE) + } else { + list(mz_twins = TRUE, mz_method = v) + } + }, + beta = function(v) list(beta = TRUE), + sparse_matrix = function(v) list(sparse = v) +) + +make_conv_args <- function(row, conv_cols) { + # row is already a named list of scalar values when called correctly + direct_cols <- setdiff(conv_cols, names(special_to_args)) + args <- row[direct_cols] + + for (nm in intersect(names(special_to_args), conv_cols)) { + args <- c(args, special_to_args[[nm]](row[[nm]])) + } + + args +} + +make_conv_args_row <- function(...) { + row <- list(...) + make_conv_args(row, scopes$conv) +} + +# Correct per-row args creation +args_by_label <- design %>% + mutate( + conv_args = pmap(select(., all_of(scopes$conv)), make_conv_args_row) + ) %>% + select(label, conv_args) %>% + deframe() + +# --------------------------- +# 6) One microbenchmark call with all expressions (correct behavior) +# --------------------------- +bench_exprs <- lapply(names(peds_by_label), function(lbl) { + bquote( + do.call(ped2com, c(list(ped = peds_by_label[[.(lbl)]]), args_by_label[[.(lbl)]])) + ) +}) +names(bench_exprs) <- names(peds_by_label) + +write_csv(design, "data-raw/ped2com_benchmark_design.csv") +# write start time +write.table(Sys.time(), "data-raw/ped2com_benchmark_start_time.txt", row.names = FALSE, col.names = FALSE) + + +benchmark_results <- do.call( + microbenchmark::microbenchmark, + c(bench_exprs, list(times = cfg$reps)) +) + +write.table(Sys.time(), "data-raw/ped2com_benchmark_end_time.txt", row.names = FALSE, col.names = FALSE) + +results <- as_tibble(benchmark_results) %>% + mutate(label = as.character(expr)) %>% + left_join( + design %>% select(-ped), + by = "label" + ) +# notes: sparse with addtwins is slow (878 vs 250), but sparse with merging is way times slower (14257 vs 241). This is a huge difference, and suggests that the merging method is not compatible with sparse matrices in its current form. The addtwins method is slower than NULL in both cases, but the difference is much more pronounced when using sparse matrices. The means of adding or substituting twins in the pedigree may interact with the way sparse matrices are constructed or handled in ped2com, leading to increased computational overhead. The NULL method is the fastest in both cases, which is expected since it does not involve any additional processing for twins. Overall, these results suggest that while the addtwins method can be used with sparse matrices, it may not be the most efficient choice, and the merging method may not be suitable for use with sparse matrices at all. +# --------------------------- +# 7) Analysis/plot +# --------------------------- +results <- results %>% + mutate( + twin_method = factor(twin_method, levels = c("NULL", "addtwins", "merging")), + ped_label = factor(ped_label, levels = levels$ped$ped_label), + gen_factor = factor(ped_label, levels = levels$ped$ped_label, labels = paste0(levels$ped$Ngen_total, " gen")) + ) +write_csv(results, "data-raw/ped2com_benchmark_results.csv") + +summary(results) +results %>% + group_by(ped_label, twin_method, sparse_matrix) %>% + summarise(median_time_ms = median(time / 1e6), + mean_time_ms = mean(time / 1e6), + var_time_ms = var(time / 1e6), + sd_time_ms = sd(time / 1e6), + se_time_ms = sd(time / 1e6) / sqrt(n()), + n_time = n(), + .groups = "drop_last") %>% + arrange(ped_label, twin_method, sparse_matrix) %>% + write_csv("data-raw/ped2com_benchmark_summary.csv") + +# ped_label twin_method sparse_matrix median_time_ms +# +# 7 highgen NULL FALSE 222. +# 8 highgen NULL TRUE 198. +# 9 highgen addtwins FALSE 250. +# 10 highgen addtwins TRUE 878. +# 11 highgen merging FALSE 241. +# 12 highgen merging TRUE 14257. + +if (cfg$reps > 8) { + notch <- FALSE +} else { + notch <- FALSE +} + +results %>% + # dplyr::filter(!ped_label %in% c("1gen", "lowgen", "midgen")) %>% + mutate( + beta_sparse = paste0("beta=", beta, ", sparse=", sparse_matrix), + beta_sparse = factor(beta_sparse, levels = c( + "beta=FALSE, sparse=FALSE", + "beta=FALSE, sparse=TRUE", + "beta=TRUE, sparse=FALSE", + "beta=TRUE, sparse=TRUE" + )) + ) %>% + ggplot( + aes( + x = ped_label, + y = time / 1e6, + fill = twin_method, + color = sparse_matrix # , + # shape = beta + ) + ) + + geom_boxplot( + notch = notch, position = position_dodge(width = 0.8), outlier.size = 0.5, + linewidth = 0.25 + ) + + scale_y_log10() + + facet_grid(~scenario) + + labs( + title = "Benchmarking ped2com() by twin handling and beta option", + x = "Pedigree", y = "Execution time (ms)", color = "sparse_matrix", fill = "Twin method" + ) + + theme_minimal() #+ +# scale_fill_manual(values = c( +# "beta=FALSE, sparse=FALSE" = "lightgray", +# "beta=FALSE, sparse=TRUE" = "gray8", +# "beta=TRUE, sparse=FALSE" = "lightcoral", +# "beta=TRUE, sparse=TRUE" = "red2" +# )) + +# scale_color_manual(values = c("NULL" = "gray", "addtwins" = "skyblue3", "merging" = "tomato2")) +# df <- sim_tbl$ped[[2]] + +# df_add <- df %>% ped2add() + +if (FALSE) { + profvis( + { + sim_tbl$ped[[1]] %>% ped2com(component = "additive", mz_twins = TRUE, mz_method = "addtwins", beta = TRUE, sparse = TRUE) + }, + interval = 0.01) + + + profvis( + { + sim_tbl$ped[[1]] %>% ped2com(component = "additive", mz_twins = TRUE, mz_method = "merging", beta = TRUE, sparse = FALSE) + }, + interval = 0.01) + # could it be the interaction between NA and gc?? + test_df_twins <- sim_tbl$ped[[1]] %>% fuseTwins(test_df_twins = TRUE) + foundreturnedall <- findMZtwins(sim_tbl$ped[[1]], + verbose = F, + returnRows = TRUE, returnIDs = TRUE, returnAsList = T + ) + foundreturnedrows <- findMZtwins(sim_tbl$ped[[1]], + verbose = F, + returnRows = TRUE, returnIDs = FALSE, returnAsList = T + ) + foundreturnedids <- findMZtwins(sim_tbl$ped[[1]], + verbose = F, + returnRows = FALSE, returnIDs = TRUE, returnAsList = T + ) + +} +if (FALSE) { + # Example of how to run a single benchmark for one condition + profvis( + { + sim_tbl$ped[[1]] %>% fuseTwins(mz_id_pairs = NULL, mz_row_pairs = NULL) + }, + interval = 0.01 + ) +} diff --git a/data-raw/ped2com_benchmark_design.csv b/data-raw/ped2com_benchmark_design.csv new file mode 100644 index 00000000..593138ce --- /dev/null +++ b/data-raw/ped2com_benchmark_design.csv @@ -0,0 +1,25 @@ +scenario,ped_label,Ngen_total,gen_twin,marR,twin_method,beta,sparse_matrix,kpc,sexR,sim_beta,component,label,ped +full,highgen,9,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive, +full,highgen,9,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive, +full,highgen,9,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive, +full,highgen,9,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive, +full,highgen,9,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive, +full,highgen,9,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive, +full,highgen,9,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive, +full,highgen,9,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive, +full,highgen,9,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive, +full,highgen,9,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive, +full,highgen,9,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive, +full,highgen,9,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive, +full,midgen,6,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive, +full,midgen,6,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive, +full,midgen,6,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive, +full,midgen,6,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive, +full,midgen,6,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive, +full,midgen,6,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive, +full,midgen,6,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive, +full,midgen,6,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive, +full,midgen,6,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive, +full,midgen,6,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive, +full,midgen,6,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive, +full,midgen,6,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive, diff --git a/data-raw/ped2com_benchmark_end_time.txt b/data-raw/ped2com_benchmark_end_time.txt new file mode 100644 index 00000000..ae7f71ae --- /dev/null +++ b/data-raw/ped2com_benchmark_end_time.txt @@ -0,0 +1 @@ +2026-02-16 10:47:50.923385 diff --git a/data-raw/ped2com_benchmark_results.csv b/data-raw/ped2com_benchmark_results.csv new file mode 100644 index 00000000..c26e7302 --- /dev/null +++ b/data-raw/ped2com_benchmark_results.csv @@ -0,0 +1,241 @@ +expr,time,label,scenario,ped_label,Ngen_total,gen_twin,marR,twin_method,beta,sparse_matrix,kpc,sexR,sim_beta,component,gen_factor +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,1097588100,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,1427187500,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,208481200,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,6189200,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,200689000,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,27103200,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,40300200,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,6026000,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,227563000,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,7928900,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,205996300,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,7277200,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,7109900,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,238252400,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,8026300,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,286128500,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,1445191500,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,1254251600,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,6127200,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,10656600,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,9625200,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,8027400,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,11509500,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,7426000,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,443970200,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,241512000,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,452599400,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,11105000,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,7619100,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,244079300,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,902782300,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,454902600,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,234666900,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,7181500,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,8527000,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,212780500,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,1100157100,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,231629800,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,245438700,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,9059000,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,214617000,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,6712400,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,445690600,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,235071900,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,253957100,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,1028123600,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,10129800,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,232238800,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,399899300,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,266978400,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,12141300,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,7296500,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,262951300,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,8368100,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,8787900,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,13507800,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,7670000,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,261880300,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,420694500,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,6956400,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,10659100,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,6674800,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,415774000,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,7061000,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,8168300,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,12491500,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,859762100,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,10746300,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,7435100,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,10309000,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,1154719100,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,686766200,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,8156100,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,1250916100,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,720033600,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,227784000,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,210911700,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,8514700,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,209204900,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,456020500,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,14974900,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,13465800,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,194586400,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,226957000,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,182487000,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,8091400,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,666622800,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,7738000,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,7558800,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,506063000,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,428734800,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,8621300,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,11728100,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,12375600,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,8772300,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,12003600,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,824795700,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,7323800,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,245685100,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,417243300,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,9535600,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,7172700,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,7072900,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,8495100,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,238015500,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,7356900,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,515113200,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,703907500,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,7065900,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,7029700,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,223874400,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,234624100,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,10971500,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,9505400,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,8480900,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,701840500,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,8603400,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,1366937400,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,12457400,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,7429200,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,221505800,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,8338300,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,7587900,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,260293200,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,6970500,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,387071800,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,7560700,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,7024200,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,7572400,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,679198200,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,7570400,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,8952100,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,220476700,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,1324076900,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,870241700,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,236028200,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,241456900,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,6418900,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,193984200,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,174968700,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,9338700,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,237393900,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,8107300,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,7594300,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,214633800,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,223070000,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,1162039400,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,6900500,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,9390800,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,6822900,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,6129800,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,419321700,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,239979200,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,246306700,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,413186000,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,945828200,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,212045300,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,516998700,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,7050200,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,185658200,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,195559600,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,10013400,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,244554400,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,10914100,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,260740500,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,696082600,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,991407800,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,183159200,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,250327800,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,11682700,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,972089100,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,7223900,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,7554000,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,248708600,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,414717900,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,234125500,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,185366000,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,1164331700,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,192493700,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,1004205400,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,409725900,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,9967800,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,7054100,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,7931100,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,9482300,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,240506300,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,178059100,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,9036700,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,8792900,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,515438000,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,7878200,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,6172500,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,678646600,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,12139500,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,705609800,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,10250500,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,8626400,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,669870100,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,8373900,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,7465900,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,7766100,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,8029300,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,9360000,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,8011500,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,5451100,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,181073500,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,708658700,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,458115900,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,212310700,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,6463400,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,1134949700,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,10318000,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,6561300,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,415959500,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,10837000,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,253422200,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,11325500,full|ped=midgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,498453800,full|ped=highgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,6816400,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,9532900,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,7317400,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,815704200,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,8478700,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,7668400,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,189231300,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,970275100,full|ped=highgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,FALSE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,224921500,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,823923000,full|ped=highgen|marR=0.8|twin=addtwins|beta=TRUE|sparse=TRUE|comp=additive,full,highgen,9,2,0.8,addtwins,TRUE,TRUE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,8553200,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,203815400,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,8441500,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,239606300,full|ped=highgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,434145700,full|ped=highgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,7041200,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,7063600,full|ped=midgen|marR=0.8|twin=merging|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,206283600,full|ped=highgen|marR=0.8|twin=NULL|beta=TRUE|sparse=FALSE|comp=additive,full,highgen,9,2,0.8,NULL,TRUE,FALSE,3,0.5,TRUE,additive,9 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,6966300,full|ped=midgen|marR=0.8|twin=NULL|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,NULL,FALSE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,6370100,full|ped=midgen|marR=0.8|twin=NULL|beta=TRUE|sparse=TRUE|comp=additive,full,midgen,6,2,0.8,NULL,TRUE,TRUE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,6395900,full|ped=midgen|marR=0.8|twin=merging|beta=TRUE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,merging,TRUE,FALSE,3,0.5,TRUE,additive,6 gen +full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,7256300,full|ped=midgen|marR=0.8|twin=addtwins|beta=FALSE|sparse=FALSE|comp=additive,full,midgen,6,2,0.8,addtwins,FALSE,FALSE,3,0.5,TRUE,additive,6 gen diff --git a/data-raw/ped2com_benchmark_start_time.txt b/data-raw/ped2com_benchmark_start_time.txt new file mode 100644 index 00000000..91c3edb0 --- /dev/null +++ b/data-raw/ped2com_benchmark_start_time.txt @@ -0,0 +1 @@ +2026-02-16 10:46:51.927423 diff --git a/data-raw/ped2com_benchmark_summary.csv b/data-raw/ped2com_benchmark_summary.csv new file mode 100644 index 00000000..5fd17eb4 --- /dev/null +++ b/data-raw/ped2com_benchmark_summary.csv @@ -0,0 +1,13 @@ +ped_label,twin_method,sparse_matrix,median_time_ms,mean_time_ms,var_time_ms,sd_time_ms,se_time_ms,n_time +midgen,NULL,FALSE,7.08545,7.065175,0.38096159250000006,0.6172208620097024,0.13801478045847118,20 +midgen,NULL,TRUE,7.010350000000001,6.97651,0.3269079230526316,0.5717586230680143,0.1278491147901759,20 +midgen,addtwins,FALSE,8.371,8.478245,1.852214922605263,1.360961029054566,0.30432013756940757,20 +midgen,addtwins,TRUE,11.5961,12.25093,13.742667178000001,3.7071103541707524,0.8289350752019123,20 +midgen,merging,FALSE,7.703200000000001,7.761935,0.6603589255526314,0.81262471384559265,0.18170841003550597,20 +midgen,merging,TRUE,9.5205,11.184985,48.470337688710515,6.962064183035841,1.556764877698468,20 +highgen,NULL,FALSE,225.3165,243.54288,3892.258395726947,62.38796675423031,13.950373464045589,20 +highgen,NULL,TRUE,194.2853,214.36891,3924.8415531188416,62.648555874168736,14.008642962683505,20 +highgen,addtwins,FALSE,261.08675,315.74309,9358.837543218842,96.74108508394374,21.631964246478915,20 +highgen,addtwins,TRUE,1081.53665,1098.186765,40060.529971101336,200.15126772294332,44.75518404112608,20 +highgen,merging,FALSE,253.68965,307.289225,7715.235275696711,87.83641201515867,19.640818816557406,20 +highgen,merging,TRUE,691.4244,704.852615,29254.359534471863,171.03905850557018,38.24549616260185,20 diff --git a/data-raw/twin_methods_sparse.png b/data-raw/twin_methods_sparse.png new file mode 100644 index 00000000..d69e5a0a Binary files /dev/null and b/data-raw/twin_methods_sparse.png differ diff --git a/man/adjustKidsPerCouple.Rd b/man/adjustKidsPerCouple.Rd index bf5b8192..dbe7e500 100644 --- a/man/adjustKidsPerCouple.Rd +++ b/man/adjustKidsPerCouple.Rd @@ -21,7 +21,18 @@ value is 3. Returns an error when kpc equals 1.} generated from a poisson distribution with mean kpc. If FALSE, the number of kids per mate will be fixed at kpc.} -\item{beta}{logical. If TRUE, use the optimized version of the algorithm.} +\item{beta}{logical or character. Controls which algorithm version to use: +\itemize{ + \item{\code{FALSE}, \code{"base"}, or \code{"original"} (default): Use the original algorithm. + Slower but ensures exact reproducibility with set.seed().} + \item{\code{TRUE} or \code{"optimized"}: Use the optimized algorithm with 4-5x speedup. + Produces statistically equivalent results but not identical to base version + due to different random number consumption. Recommended for large simulations + where speed matters more than exact reproducibility.} +} +Note: Both versions are mathematically correct and produce valid pedigrees with the +same statistical properties (sex ratios, mating rates, etc.). The optimized version +uses vectorized operations instead of loops, making it much faster for large pedigrees.} } \value{ A numeric vector with the generated or adjusted number of kids per couple. diff --git a/man/buildBetweenGenerations.Rd b/man/buildBetweenGenerations.Rd index fa87604e..f28622b2 100644 --- a/man/buildBetweenGenerations.Rd +++ b/man/buildBetweenGenerations.Rd @@ -62,7 +62,18 @@ kids per mate will be fixed at kpc.} \item{code_female}{The value to use for females. Default is "F"} -\item{beta}{logical. If TRUE, use the optimized version of the algorithm.} +\item{beta}{logical or character. Controls which algorithm version to use: +\itemize{ + \item{\code{FALSE}, \code{"base"}, or \code{"original"} (default): Use the original algorithm. + Slower but ensures exact reproducibility with set.seed().} + \item{\code{TRUE} or \code{"optimized"}: Use the optimized algorithm with 4-5x speedup. + Produces statistically equivalent results but not identical to base version + due to different random number consumption. Recommended for large simulations + where speed matters more than exact reproducibility.} +} +Note: Both versions are mathematically correct and produce valid pedigrees with the +same statistical properties (sex ratios, mating rates, etc.). The optimized version +uses vectorized operations instead of loops, making it much faster for large pedigrees.} } \value{ The function updates the `df_Fam` data frame in place, adding or modifying columns related to parental and offspring status, diff --git a/man/buildWithinGenerations.Rd b/man/buildWithinGenerations.Rd index 905e42ad..cbddf27a 100644 --- a/man/buildWithinGenerations.Rd +++ b/man/buildWithinGenerations.Rd @@ -20,7 +20,18 @@ buildWithinGenerations( ) } \arguments{ -\item{beta}{logical. If TRUE, use the optimized version of the algorithm.} +\item{beta}{logical or character. Controls which algorithm version to use: +\itemize{ + \item{\code{FALSE}, \code{"base"}, or \code{"original"} (default): Use the original algorithm. + Slower but ensures exact reproducibility with set.seed().} + \item{\code{TRUE} or \code{"optimized"}: Use the optimized algorithm with 4-5x speedup. + Produces statistically equivalent results but not identical to base version + due to different random number consumption. Recommended for large simulations + where speed matters more than exact reproducibility.} +} +Note: Both versions are mathematically correct and produce valid pedigrees with the +same statistical properties (sex ratios, mating rates, etc.). The optimized version +uses vectorized operations instead of loops, making it much faster for large pedigrees.} \item{sizeGens}{A numeric vector containing the sizes of each generation within the pedigree.} diff --git a/man/dot-adjBeta.Rd b/man/dot-adjBeta.Rd index 5b962868..3c572806 100644 --- a/man/dot-adjBeta.Rd +++ b/man/dot-adjBeta.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/constructAdjacency.R \name{.adjBeta} \alias{.adjBeta} -\title{Construct Adjacency Matrix for Parent-Child Relationships Using Beta Method +\title{Construct Adjacency Matrix for Parent-Child Relationships Using Beta Methods This function constructs an adjacency matrix for parent-child relationships using a method in beta testing. It identifies parent-child pairs based on the specified component of relatedness.} @@ -60,7 +60,7 @@ specified component of relatedness.} \item{...}{additional arguments to be passed to \code{\link{ped2com}}} } \description{ -Construct Adjacency Matrix for Parent-Child Relationships Using Beta Method +Construct Adjacency Matrix for Parent-Child Relationships Using Beta Methods This function constructs an adjacency matrix for parent-child relationships using a method in beta testing. It identifies parent-child pairs based on the specified component of relatedness. diff --git a/man/dropLink.Rd b/man/dropLink.Rd index 668d9e82..2753587d 100644 --- a/man/dropLink.Rd +++ b/man/dropLink.Rd @@ -9,7 +9,8 @@ dropLink( ID_drop = NA_integer_, gen_drop = 2, sex_drop = NA_character_, - n_drop = 1 + n_drop = 1, + verbose = FALSE ) } \arguments{ @@ -22,6 +23,8 @@ dropLink( \item{sex_drop}{the biological sex of the randomly dropped person.} \item{n_drop}{the number of times the mutation happens.} + +\item{verbose}{logical. If TRUE, print progress through stages of algorithm} } \value{ a pedigree with the dropped person's `dadID` and `momID` set to NA. diff --git a/man/findMZtwins.Rd b/man/findMZtwins.Rd new file mode 100644 index 00000000..39a361fd --- /dev/null +++ b/man/findMZtwins.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpTwins.R +\name{findMZtwins} +\alias{findMZtwins} +\title{Find MZ twin pair_rows in a pedigree} +\usage{ +findMZtwins( + ped, + verbose = FALSE, + returnRows = TRUE, + returnIDs = FALSE, + returnAsList = TRUE, + beta = FALSE +) +} +\arguments{ +\item{ped}{A pedigree data.frame with columns \code{ID} and \code{twinID}. +Optionally a \code{zygosity} column; when present only pair_rows where both +members have \code{zygosity == "MZ"} are used.} + +\item{verbose}{logical. If TRUE, print progress messages.} + +\item{returnRows}{logical. If TRUE, return the row indices of the twin pair_rows instead of IDs.} + +\item{returnIDs}{logical. If TRUE, return the IDs of the twin pair_rows instead of row indices.} + +\item{returnAsList}{logical. If TRUE, return results as a list of vectors +(default). If FALSE, return results as a data.frame with separate columns for each twin's ID and row index. +@param beta logical. If TRUE, use an optimized approach with O(1) lookups for large pedigrees. If FALSE (default), use a simpler approach that may be less efficient for large pedigrees.} +} +\value{ +A list of length-2 integer vectors \code{c(idx1, idx2)} giving the + row indices of each MZ pair in the pedigree, or \code{NULL} if none found. +} +\description{ +Identifies MZ twin pair_rows from the \code{twinID} column and returns their +row indices. These indices are used later to merge the twins' columns in +the \code{r2} matrix before \code{tcrossprod}, which correctly produces +relatedness 1 between MZ co-twins with no diagonal or downstream artifacts. +} +\keyword{internal} diff --git a/man/fuseTwins.Rd b/man/fuseTwins.Rd new file mode 100644 index 00000000..42eb49fa --- /dev/null +++ b/man/fuseTwins.Rd @@ -0,0 +1,39 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpTwins.R +\name{fuseTwins} +\alias{fuseTwins} +\title{Fuse MZ twin pairs in a pedigree dataset for path tracing +This function identifies MZ twin pairs in the pedigree dataset and merges their IDs for path tracing purposes. The second twin in each pair is made a founder (with NA parents), and all children of the second twin are redirected to the first twin. This allows for correct relatedness calculations without diagonal or downstream artifacts.} +\usage{ +fuseTwins( + ped, + df_twins = NULL, + mz_id_pairs = NULL, + mz_row_pairs = NULL, + config = list(verbose = FALSE), + test_df_twins = FALSE, + beta = FALSE +) +} +\arguments{ +\item{ped}{A pedigree data.frame with columns \code{ID}, \code{momID}, \code{dadID}, and optionally \code{twinID} and \code{zygosity}. The function will look for MZ twin pairs based on the \code{twinID} column and optionally restrict to MZ pairs if a \code{zygosity} column is present.} + +\item{df_twins}{Optional data frame with columns \code{twin1_id}, \code{twin2_id}, \code{twin1_row}, and \code{twin2_row} specifying the IDs and row indices of MZ twin pairs to fuse. If provided, this will be used instead of automatically identifying MZ twins from the \code{twinID} column. If this parameter is provided, it takes precedence over \code{mz_id_pairs} and \code{mz_row_pairs}. If \code{test_df_twins} is TRUE, this data frame will be returned for testing purposes instead of performing the fusion.} + +\item{mz_id_pairs}{Optional list of length-2 character vectors specifying the IDs of MZ twin pairs to fuse. If provided, this will be used instead of automatically identifying MZ twins from the \code{twinID} column. Each element should be a character vector of length 2, e.g. \code{list(c("ID1", "ID2"), c("ID3", "ID4"))}.} + +\item{mz_row_pairs}{Optional list of length-2 integer vectors specifying the row indices of MZ twin pairs to fuse. If provided, this will be used instead of automatically identifying MZ twins from the \code{twinID} column. Each element should be an integer vector of length 2, e.g. \code{list(c(1, 2), c(3, 4))}.} + +\item{config}{A list of configuration options.} + +\item{test_df_twins}{logical. If TRUE, return the data frame of twin pairs instead of the modified pedigree. Default is FALSE.} + +\item{beta}{logical. If TRUE, use an optimized approach with O(1) lookups for large pedigrees when identifying MZ twins. Default is FALSE.} +} +\value{ +A modified version of the input pedigree data.frame with MZ twin pairs fused for path tracing. If \code{test_df_twins} is TRUE, returns the data frame of identified twin pairs instead. +} +\description{ +Fuse MZ twin pairs in a pedigree dataset for path tracing +This function identifies MZ twin pairs in the pedigree dataset and merges their IDs for path tracing purposes. The second twin in each pair is made a founder (with NA parents), and all children of the second twin are redirected to the first twin. This allows for correct relatedness calculations without diagonal or downstream artifacts. +} diff --git a/man/isTwin.Rd b/man/isTwin.Rd new file mode 100644 index 00000000..e2f3e887 --- /dev/null +++ b/man/isTwin.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpTwins.R +\name{isTwin} +\alias{isTwin} +\title{Determine isTwin Status} +\usage{ +isTwin(ped) +} +\arguments{ +\item{ped}{pedigree data frame} +} +\value{ +isTwin 'S' matrix +} +\description{ +Determine isTwin Status +} +\keyword{internal} diff --git a/man/makeInbreeding.Rd b/man/makeInbreeding.Rd index 1540b170..d43200c5 100644 --- a/man/makeInbreeding.Rd +++ b/man/makeInbreeding.Rd @@ -10,7 +10,8 @@ makeInbreeding( ID_mate2 = NA_integer_, verbose = FALSE, gen_inbred = 2, - type_inbred = "sib" + type_inbred = "sib", + prefer_unmated = FALSE ) } \arguments{ @@ -25,6 +26,8 @@ makeInbreeding( \item{gen_inbred}{A vector of \code{generation} of the twin to be imputed.} \item{type_inbred}{A character vector indicating the type of inbreeding. "sib" for sibling inbreeding and "cousin" for cousin inbreeding.} + +\item{prefer_unmated}{A logical indicating whether to prefer unmated siblings when automatically selecting inbred mates. Default is FALSE, which means the function will consider all siblings regardless of their mating status.} } \value{ Returns a \code{data.frame} with some inbred mates. diff --git a/man/makeTwins.Rd b/man/makeTwins.Rd index cb9733e5..dc4f3ee5 100644 --- a/man/makeTwins.Rd +++ b/man/makeTwins.Rd @@ -10,7 +10,8 @@ makeTwins( ID_twin2 = NA_integer_, gen_twin = 2, verbose = FALSE, - zygosity = "MZ" + zygosity = "MZ", + twin_sex = "R" ) } \arguments{ @@ -25,6 +26,8 @@ makeTwins( \item{verbose}{logical. If TRUE, print progress through stages of algorithm} \item{zygosity}{A character string indicating the zygosity of the twins. Default is "MZ" for monozygotic twins.} + +\item{twin_sex}{A character string indicating the sex of the twins. Default is randomly assigned ("R"). If specified, it should be either "M" or "F"} } \value{ Returns a \code{data.frame} with MZ twins information added as a new column. diff --git a/man/ped2add.Rd b/man/ped2add.Rd index dc2fee88..0a08cd09 100644 --- a/man/ped2add.Rd +++ b/man/ped2add.Rd @@ -21,6 +21,8 @@ ped2add( save_rate_parlist = 1e+05 * save_rate, save_path = "checkpoint/", compress = TRUE, + mz_twins = FALSE, + mz_method = "addtwins", ... ) } @@ -57,6 +59,10 @@ ped2add( \item{compress}{logical. If TRUE, use compression when saving the checkpoint files. Defaults to TRUE.} +\item{mz_twins}{logical. If TRUE, merge MZ co-twin columns in the r2 matrix before tcrossprod so that MZ twins are coded with relatedness 1 instead of 0.5. Twin pairs are identified from the \code{twinID} column. When a \code{zygosity} column is also present, only pairs where both members have \code{zygosity == "MZ"} are used; otherwise all \code{twinID} pairs are assumed to be MZ. Defaults to FALSE.} + +\item{mz_method}{character. The method to handle MZ twins. Options are "merging" (default) or "addtwins". "addtwins" adds the twin2 column to the twin1 column before tcrossprod so that all relatedness flows through a single source, then leaves the twin2 column as zero and relies on the fact that the row/col names are the same to copy the values back to twin2 after tcrossprod. "merging" merges the twin2 column into the twin1 column before tcrossprod and then copies the values back to twin2 after tcrossprod so that both twins appear in the final matrix.} + \item{...}{additional arguments to be passed to \code{\link{ped2com}}} } \description{ diff --git a/man/ped2com.Rd b/man/ped2com.Rd index f34d6022..dc57cdb7 100644 --- a/man/ped2com.Rd +++ b/man/ped2com.Rd @@ -15,7 +15,7 @@ ped2com( standardize_colnames = TRUE, transpose_method = "tcrossprod", adjacency_method = "direct", - isChild_method = "classic", + isChild_method = "partialparent", saveable = FALSE, resume = FALSE, save_rate = 5, @@ -25,6 +25,9 @@ ped2com( save_path = "checkpoint/", adjBeta_method = NULL, compress = TRUE, + mz_twins = TRUE, + mz_method = "addtwins", + beta = FALSE, ... ) } @@ -69,6 +72,12 @@ ped2com( \item{compress}{logical. If TRUE, use compression when saving the checkpoint files. Defaults to TRUE.} +\item{mz_twins}{logical. If TRUE, merge MZ co-twin columns in the r2 matrix before tcrossprod so that MZ twins are coded with relatedness 1 instead of 0.5. Twin pairs are identified from the \code{twinID} column. When a \code{zygosity} column is also present, only pairs where both members have \code{zygosity == "MZ"} are used; otherwise all \code{twinID} pairs are assumed to be MZ. Defaults to FALSE.} + +\item{mz_method}{character. The method to handle MZ twins. Options are "merging" (default) or "addtwins". "addtwins" adds the twin2 column to the twin1 column before tcrossprod so that all relatedness flows through a single source, then leaves the twin2 column as zero and relies on the fact that the row/col names are the same to copy the values back to twin2 after tcrossprod. "merging" merges the twin2 column into the twin1 column before tcrossprod and then copies the values back to twin2 after tcrossprod so that both twins appear in the final matrix.} + +\item{beta}{logical. Used for benchmarking} + \item{...}{additional arguments to be passed to \code{\link{ped2com}}} } \description{ diff --git a/man/ped2fam.Rd b/man/ped2fam.Rd index 3052e568..b3cf1e07 100644 --- a/man/ped2fam.Rd +++ b/man/ped2fam.Rd @@ -10,6 +10,7 @@ ped2fam( momID = "momID", dadID = "dadID", famID = "famID", + twinID = "twinID", ... ) } @@ -24,6 +25,8 @@ ped2fam( \item{famID}{character. Name of the column to be created in ped for the family ID variable} +\item{twinID}{character. Name of the column in ped for the twin ID variable, if applicable} + \item{...}{additional arguments to be passed to \code{\link{ped2com}}} } \value{ diff --git a/man/ped2graph.Rd b/man/ped2graph.Rd index 5e7ac7b1..022af3c7 100755 --- a/man/ped2graph.Rd +++ b/man/ped2graph.Rd @@ -9,6 +9,7 @@ ped2graph( personID = "ID", momID = "momID", dadID = "dadID", + twinID = "twinID", directed = TRUE, adjacent = c("parents", "mothers", "fathers"), ... @@ -23,6 +24,8 @@ ped2graph( \item{dadID}{character. Name of the column in ped for the father ID variable} +\item{twinID}{character. Name of the column in ped for the twin ID variable, if applicable} + \item{directed}{Logical scalar. Default is TRUE. Indicates whether or not to create a directed graph.} \item{adjacent}{Character. Relationship that defines adjacency in the graph: parents, mothers, or fathers} diff --git a/man/ped2maternal.Rd b/man/ped2maternal.Rd index 03e02311..1f3bcb2e 100755 --- a/man/ped2maternal.Rd +++ b/man/ped2maternal.Rd @@ -10,6 +10,7 @@ ped2maternal( momID = "momID", dadID = "dadID", matID = "matID", + twinID = "twinID", ... ) } @@ -24,6 +25,8 @@ ped2maternal( \item{matID}{Character. Maternal line ID variable to be created and added to the pedigree} +\item{twinID}{character. Name of the column in ped for the twin ID variable, if applicable} + \item{...}{additional arguments to be passed to \code{\link{ped2com}}} } \description{ diff --git a/man/ped2paternal.Rd b/man/ped2paternal.Rd index e893ec03..7fac27b9 100755 --- a/man/ped2paternal.Rd +++ b/man/ped2paternal.Rd @@ -10,6 +10,7 @@ ped2paternal( momID = "momID", dadID = "dadID", patID = "patID", + twinID = "twinID", ... ) } @@ -24,6 +25,8 @@ ped2paternal( \item{patID}{Character. Paternal line ID variable to be created and added to the pedigree} +\item{twinID}{character. Name of the column in ped for the twin ID variable, if applicable} + \item{...}{additional arguments to be passed to \code{\link{ped2com}}} } \description{ diff --git a/man/simulatePedigree.Rd b/man/simulatePedigree.Rd index fa87632e..d56ab72c 100644 --- a/man/simulatePedigree.Rd +++ b/man/simulatePedigree.Rd @@ -76,7 +76,18 @@ current version.} \item{fam_shift}{An integer to shift the person ID. Default is 1L. This is useful when simulating multiple pedigrees to avoid ID conflicts.} -\item{beta}{logical. If TRUE, use the optimized version of the algorithm.} +\item{beta}{logical or character. Controls which algorithm version to use: +\itemize{ + \item{\code{FALSE}, \code{"base"}, or \code{"original"} (default): Use the original algorithm. + Slower but ensures exact reproducibility with set.seed().} + \item{\code{TRUE} or \code{"optimized"}: Use the optimized algorithm with 4-5x speedup. + Produces statistically equivalent results but not identical to base version + due to different random number consumption. Recommended for large simulations + where speed matters more than exact reproducibility.} +} +Note: Both versions are mathematically correct and produce valid pedigrees with the +same statistical properties (sex ratios, mating rates, etc.). The optimized version +uses vectorized operations instead of loops, making it much faster for large pedigrees.} \item{...}{Additional arguments to be passed to other functions.} } diff --git a/tests/testthat/test-buildComponent.R b/tests/testthat/test-buildComponent.R index b00fa039..d863b0f7 100644 --- a/tests/testthat/test-buildComponent.R +++ b/tests/testthat/test-buildComponent.R @@ -1,3 +1,245 @@ +test_that("MZ twins coded at relatedness 1 via twinID column", { + # Simple pedigree: two parents and two MZ twin children + ped <- potter + + mz_method_opts <- c("addtwins", "merging") + + for (mz_method in mz_method_opts) { + # mz_method <- "merging" # "addtwins" + # Without mz_twins: siblings get 0.5 + r_no_mz <- ped2add(ped, mz_twins = FALSE, sparse = FALSE, mz_method = mz_method) + expect_equal(r_no_mz["12", "13"], 0.5) + expect_equal(r_no_mz["13", "12"], 0.5) + + # With mz_twins: MZ twins get 1.0 + r_mz <- ped2add(ped, mz_twins = TRUE, sparse = FALSE, mz_method = mz_method) + expect_equal(r_mz["12", "13"], 1.0) + expect_equal(r_mz["13", "12"], 1.0) + + # Self-relatedness should still be 1 + expect_equal(r_mz["12", "12"], 1.0) + expect_equal(r_mz["13", "13"], 1.0) + + # Parent-child relatedness unchanged + expect_equal(r_mz["12", "9"], 0.5) + expect_equal(r_mz["13", "9"], 0.5) + expect_equal(r_mz["12", "10"], 0.5) + expect_equal(r_mz["13", "10"], 0.5) + } + ped_kids <- potter + + # Add a child to one of the MZ twins + ped_kids <- addPersonToPed(ped_kids, sex = 0, momID = NA, dadID = NA, personID = 31) + ped_kids <- addPersonToPed(ped_kids, sex = 0, momID = NA, dadID = NA, personID = 32) + ped_kids <- addPersonToPed(ped_kids, sex = 0, momID = 31, dadID = 12, personID = 33) + ped_kids <- addPersonToPed(ped_kids, sex = 0, momID = 32, dadID = 13, personID = 34) + ped_kids <- addPersonToPed(ped_kids, sex = 0, momID = 31, dadID = 13, personID = 35) + + for (mz_method in mz_method_opts) { + r_kids <- ped2add(ped_kids, mz_twins = TRUE, sparse = FALSE, mz_method = mz_method) + # Child of twin1 (ID=31) should be 0.5 to twin1 (parent) + expect_equal(r_kids["33", "12"], 0.5) + # Child of twin1 should ALSO be 0.5 to twin2 (genetically identical to parent) + expect_equal(r_kids["33", "13"], 0.5) + # Child of twin2 (ID=32) should be 0.5 to twin + expect_equal(r_kids["34", "13"], 0.5) + # Child of twin2 should ALSO be 0.5 to twin1 (genetically identical to parent) + expect_equal(r_kids["34", "12"], 0.5) + + # different moms should be 0.25 with different mz twin dads + expect_equal(r_kids["34", "33"], 0.25) + expect_equal(r_kids["34", "35"], 0.25) + + # same mom, different mz twin dads should be 0.5 + expect_equal(r_kids["33", "35"], 0.5) + } + + r_mz1 <- ped2add(ped_kids, mz_twins = TRUE, sparse = FALSE, mz_method = mz_method_opts[1]) + + r_mz2 <- ped2add(ped_kids, mz_twins = TRUE, sparse = FALSE, mz_method = mz_method_opts[2]) + + expect_equal(r_mz1, r_mz2) +}) + + +test_that("MZ twins coded at relatedness 1 via twinID column (complex pedigree)", { + set.seed(1667) + Ngen <- 5 + kpc <- 4 + sexR <- .50 # sometimes fails above .5 + marR <- 1 + + gen_twin <- 3 # max(c(floor(Ngen / 2) - 1, 2)) + + # create base pedigree with twins at gen_twin + + df_midgen_base <- simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR, beta = TRUE) |> + makeTwins(gen_twin = gen_twin) + + twinIDS <- findMZtwins(df_midgen_base, returnRows = FALSE, + returnIDs = TRUE, returnAsList = FALSE) + + moms_of_twins <- df_midgen_base$momID[df_midgen_base$ID %in% c(twinIDS$twin1_id, twinIDS$twin2_id)] %>% unique() + dads_of_twins <- df_midgen_base$dadID[df_midgen_base$ID %in% c(twinIDS$twin1_id, twinIDS$twin2_id)] %>% unique() + parents_of_twins <- unique(c(moms_of_twins, dads_of_twins)) + + grandmothers_of_twins <- df_midgen_base$momID[df_midgen_base$ID %in% parents_of_twins] %>% unique() + grandfathers_of_twins <- df_midgen_base$dadID[df_midgen_base$ID %in% parents_of_twins] %>% unique() + grandfathers_of_twins <- grandfathers_of_twins[!is.na(grandfathers_of_twins)] + grandmothers_of_twins <- grandmothers_of_twins[!is.na(grandmothers_of_twins)] + grandparents_of_twins <- unique(c(grandmothers_of_twins, grandfathers_of_twins)) + + female_children_of_twins <- df_midgen_base$ID[(df_midgen_base$momID %in% c(twinIDS$twin1_id, twinIDS$twin2_id) | df_midgen_base$dadID %in% c(twinIDS$twin1_id, twinIDS$twin2_id)) & df_midgen_base$sex == "F"] %>% unique() + + male_children_of_twins <- df_midgen_base$ID[(df_midgen_base$momID %in% c(twinIDS$twin1_id, twinIDS$twin2_id) | df_midgen_base$dadID %in% c(twinIDS$twin1_id, twinIDS$twin2_id)) & df_midgen_base$sex == "M"] %>% unique() + + children_of_twins <- c(male_children_of_twins, female_children_of_twins) + + df_midgen_below <- df_midgen_base |> makeInbreeding(ID_mate1 = male_children_of_twins[1], ID_mate2 = female_children_of_twins[length(female_children_of_twins)]) + + df_midgen_above <- df_midgen_base + + df_midgen_above$momID[df_midgen_above$ID %in% parents_of_twins] <- grandmothers_of_twins[1] + df_midgen_above$dadID[df_midgen_above$ID %in% parents_of_twins] <- grandfathers_of_twins[1] + + if (FALSE) { + df_midgen_below %>% rename(personID = ID) %>% + ggpedigree::ggPedigreeInteractive(config = list(code_male = "M", focal_fill_personID = twinIDS$twin1_id, + focal_fill_include = TRUE, + sex_color_include = FALSE)) + + } + for (df_midgen in list(df_midgen_base, df_midgen_below, df_midgen_above)) { + r_mz1 <- df_midgen |> + ped2add(mz_method = "merging", mz_twins = TRUE) + r_mz2 <- df_midgen |> + ped2add(mz_method = "addtwins", mz_twins = TRUE) + + r_mz3 <- df_midgen |> + ped2add(mz_twins = FALSE) + + # which rows are the twins + twin_rows <- which(!is.na(df_midgen$twinID)) + child_rows <- which(df_midgen$momID %in% df_midgen$ID[twin_rows] | df_midgen$dadID %in% df_midgen$ID[twin_rows]) + + family_rows <- unique(c(twin_rows, child_rows)) + + expect_equal(sum(as.matrix(r_mz1[family_rows, family_rows]) - as.matrix(r_mz2[family_rows, family_rows])), 0) + + expect_gt(sum(as.matrix(r_mz1[family_rows, family_rows])), sum(as.matrix(r_mz3[family_rows, family_rows]))) + + expect_gt(sum(as.matrix(r_mz2[family_rows, family_rows])), sum(as.matrix(r_mz3[family_rows, family_rows]))) + + + r_mz1_ordered <- r_mz1[order(rownames(r_mz1)), order(colnames(r_mz1))] + r_mz2_ordered <- r_mz2[order(rownames(r_mz2)), order(colnames(r_mz2))] + + expect_equal(sum(r_mz1_ordered - r_mz2_ordered), 0) + + expect_equal(length(r_mz1@i), length(r_mz2@i)) + expect_equal(length(r_mz1@x), length(r_mz2@x)) + expect_equal(length(r_mz1@p), length(r_mz2@p)) + + expect_equal(length(r_mz1@i), length(r_mz3@i)) + expect_equal(length(r_mz2@i), length(r_mz3@i)) + + expect_equal(length(r_mz1@x), length(r_mz3@x)) + expect_equal(length(r_mz2@x), length(r_mz3@x)) + + expect_equal(length(r_mz1@p), length(r_mz3@p)) + expect_equal(length(r_mz2@p), length(r_mz3@p)) + } +}) + +test_that("MZ twins coded at relatedness 1 via twinID column (minimal data.frame)", { + # Simple pedigree: two parents and two MZ twin children + ped <- data.frame( + ID = c(1, 2, 3, 4), + momID = c(NA, NA, 2, 2), + dadID = c(NA, NA, 1, 1), + sex = c("M", "F", "M", "M"), + twinID = c(NA, NA, 4, 3), + zygosity = c(NA, NA, "MZ", "MZ") + ) + + # Without mz_twins: siblings get 0.5 + r_no_mz <- ped2add(ped, mz_twins = FALSE, sparse = FALSE) + expect_equal(r_no_mz["3", "4"], 0.5) + expect_equal(r_no_mz["4", "3"], 0.5) + + # With mz_twins: MZ twins get 1.0 + r_mz <- ped2add(ped, mz_twins = TRUE, sparse = FALSE) + expect_equal(r_mz["3", "4"], 1.0) + expect_equal(r_mz["4", "3"], 1.0) + + # Self-relatedness should still be 1 + expect_equal(r_mz["3", "3"], 1.0) + expect_equal(r_mz["4", "4"], 1.0) + + # Parent-child relatedness unchanged + expect_equal(r_mz["3", "1"], 0.5) + expect_equal(r_mz["4", "1"], 0.5) + expect_equal(r_mz["3", "2"], 0.5) + expect_equal(r_mz["4", "2"], 0.5) +}) + +test_that("MZ twins without zygosity column assumes all twinID pairs are MZ", { + ped <- data.frame( + ID = c(1, 2, 3, 4), + momID = c(NA, NA, 2, 2), + dadID = c(NA, NA, 1, 1), + sex = c("M", "F", "M", "M"), + twinID = c(NA, NA, 4, 3) + ) + + r_mz <- ped2add(ped, mz_twins = TRUE, sparse = FALSE) + expect_equal(r_mz["3", "4"], 1.0) + expect_equal(r_mz["4", "3"], 1.0) +}) + +test_that("DZ twins with zygosity column are NOT modified", { + ped <- data.frame( + ID = c(1, 2, 3, 4), + momID = c(NA, NA, 2, 2), + dadID = c(NA, NA, 1, 1), + sex = c("M", "F", "M", "F"), + twinID = c(NA, NA, 4, 3), + zygosity = c(NA, NA, "DZ", "DZ") + ) + + r_mz <- ped2add(ped, mz_twins = TRUE, sparse = FALSE) + # DZ twins remain at sibling relatedness = 0.5 + expect_equal(r_mz["3", "4"], 0.5) + expect_equal(r_mz["4", "3"], 0.5) +}) + +test_that("MZ twins: downstream child relatedness is correct", { + # 3-generation pedigree: parents -> MZ twins -> twin2 has a child + ped <- data.frame( + ID = c(1, 2, 3, 4, 5, 6), + momID = c(NA, NA, 2, 2, NA, 4), + dadID = c(NA, NA, 1, 1, NA, 5), + sex = c("M", "F", "M", "M", "F", "M"), + twinID = c(NA, NA, 4, 3, NA, NA), + zygosity = c(NA, NA, "MZ", "MZ", NA, NA) + ) + + r_mz <- ped2add(ped, mz_twins = TRUE, sparse = FALSE) + + # MZ twins at 1.0 + expect_equal(r_mz["3", "4"], 1.0) + + # Child of twin2 (ID=4) should be 0.5 to twin2 (parent) + expect_equal(r_mz["6", "4"], 0.5) + + # Child of twin2 should ALSO be 0.5 to twin1 (genetically identical to parent) + expect_equal(r_mz["6", "3"], 0.5) + + # Diagonal for both twins should be clean (no inflation) + expect_equal(r_mz["3", "3"], 1.0) + expect_equal(r_mz["4", "4"], 1.0) +}) + test_that(".assignParentValue works", { expect_equal(.assignParentValue("generation"), .5) expect_equal(.assignParentValue("additive"), .5) diff --git a/tests/testthat/test-helpTwins.R b/tests/testthat/test-helpTwins.R new file mode 100644 index 00000000..e19643bb --- /dev/null +++ b/tests/testthat/test-helpTwins.R @@ -0,0 +1,100 @@ +test_that("fuse twins behaves", { + # Simple pedigree: two parents and two MZ twin children + ped1 <- potter + ped1$ID <- ped1$personID + ped2 <- ped1 + ped2$famID <- 2 + + ped2$ID <- ped2$personID + 100 + ped2$momID <- ped2$momID + 100 + ped2$dadID <- ped2$dadID + 100 + ped2$twinID <- ped2$twinID + 100 + + ped <- rbind(ped1, ped2) + remove(ped2) + remove(ped1) + # returnRows = TRUE, + # returnIDs = FALSE, + # returnAsList = TRUE + + returnedRowsList <- findMZtwins(ped, returnRows = TRUE, returnIDs = FALSE, returnAsList = TRUE) + returnIDsList <- findMZtwins(ped, returnRows = FALSE, returnIDs = TRUE, returnAsList = TRUE) + returnedBothList <- findMZtwins(ped, returnRows = T, returnIDs = T, returnAsList = TRUE) + + # no error should be thrown when running fuseTwins with any of the above outputs as arguments + expect_no_error( + fuseTwins(ped, + test_df_twins = TRUE, + mz_id_pairs = NULL, + mz_row_pairs = NULL + ) + ) + + expect_no_error( + fuseTwins(ped, + test_df_twins = TRUE, + mz_id_pairs = NULL, + mz_row_pairs = returnedRowsList + ) + ) + expect_no_error( + fuseTwins(ped, + test_df_twins = TRUE, + mz_id_pairs = returnIDsList, + mz_row_pairs = NULL + ) + ) + expect_no_error( + fuseTwins(ped, + test_df_twins = TRUE, + mz_id_pairs = returnedBothList$pair_ids, + mz_row_pairs = returnedBothList$pair_rows + ) + ) + + + df_null <- tryCatch( + fuseTwins(ped, + test_df_twins = TRUE, + mz_id_pairs = NULL, + mz_row_pairs = NULL + ), + error = function(e) e + ) + + + df_returnedRows <- tryCatch( + fuseTwins(ped, + test_df_twins = TRUE, + mz_id_pairs = NULL, + mz_row_pairs = returnedRowsList + ), + error = function(e) e + ) + + + df_returnedIDs <- tryCatch( + fuseTwins(ped, + test_df_twins = TRUE, + mz_id_pairs = returnIDsList, + mz_row_pairs = NULL + ), + error = function(e) e + ) + + + df_returnedBoth <- tryCatch( + fuseTwins(ped, + test_df_twins = TRUE, + mz_id_pairs = returnedBothList$pair_ids, + mz_row_pairs = returnedBothList$pair_rows + ), + error = function(e) e + ) + + + expect_equal(df_returnedRows, df_returnedIDs) + expect_equal(df_returnedRows, df_returnedBoth) + expect_equal(df_returnedRows, df_null) + expect_equal(nrow(df_returnedRows), 2) # One pair of twins should returned +}) diff --git a/tests/testthat/test-segmentPedigree.R b/tests/testthat/test-segmentPedigree.R index a0f64e08..cf1cce02 100644 --- a/tests/testthat/test-segmentPedigree.R +++ b/tests/testthat/test-segmentPedigree.R @@ -1,9 +1,11 @@ test_that("ped2fam is smart about string ids", { data(hazard) + # ped2fam should work with numeric IDs and produce numeric IDs ds_num <- ped2fam(hazard, famID = "newFamID") expect_true(is.numeric(ds_num$ID)) expect_true(is.numeric(ds_num$newFamID)) hazard$ID_og <- hazard$ID + # ped2fam should work with string IDs and produce string IDs hazard$ID <- paste0("ID", hazard$ID) hazard$dadID <- paste0("ID", hazard$dadID) hazard$dadID[hazard$dadID == "IDNA"] <- NA diff --git a/tests/testthat/test-simulatePedigree.R b/tests/testthat/test-simulatePedigree.R index ea4880af..a13ff5d1 100644 --- a/tests/testthat/test-simulatePedigree.R +++ b/tests/testthat/test-simulatePedigree.R @@ -20,7 +20,7 @@ test_that("simulated pedigree generates expected data structure", { if (isFALSE(beta) || (isTRUE(beta) && beta_match_base)) { expect_equal(length(results$ID), base_length, tolerance = strict_tolerance) } else { - expect_true(length(results$ID) >= base_length - base_length_tol * base_length && length(results$ID) <= base_length_tol * base_length + base_length, + expect_true(length(results$ID) >= base_length - base_length_tol && length(results$ID) <= base_length + base_length_tol, info = paste0("Beta=TRUE: Expected 45-70 individuals, got ", length(results$ID)) ) } diff --git a/tests/testthat/test-tweakPedigree.R b/tests/testthat/test-tweakPedigree.R index b0e63987..47625e7f 100644 --- a/tests/testthat/test-tweakPedigree.R +++ b/tests/testthat/test-tweakPedigree.R @@ -163,17 +163,24 @@ test_that("makeInbreeding - Inbred mates specified by generation and sibling", { marR <- .7 gen_inbred <- 2 type_inbred <- "sibling" + prefer_unmated <- c(TRUE, FALSE) + ped <- simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR) # - result <- makeInbreeding(ped, gen_inbred = gen_inbred, type_inbred = type_inbred) - expect_equal(names(result), c("famID", "ID", "gen", "dadID", "momID", "spID", "sex")) + for (prefer in prefer_unmated) { - # do we have the same people? - expect_equal(result$ID, ped$ID) + result <- makeInbreeding(ped, gen_inbred = gen_inbred, type_inbred = type_inbred, + prefer_unmated = prefer, + verbose = TRUE) + expect_equal(names(result), c("famID", "ID", "gen", "dadID", "momID", "spID", "sex")) - # did we get more spID values than we started with? - expect_gt(sum(!is.na(result$spID)), sum(!is.na(ped$spID))) + # do we have the same people? + expect_equal(result$ID, ped$ID) + + # did we get more spID values than we started with? + expect_gt(sum(!is.na(result$spID)), sum(!is.na(ped$spID))) + } }) test_that("makeInbreeding - Inbred mates specified by generation and cousin", { @@ -221,6 +228,9 @@ test_that("dropLink - Drop specified by ID", { # are the dataframes the same in both the undropped and dropepd relationships for all but the dropped ID? + expect_equal(colnames(result), c("famID", "ID", "gen", "dadID", "momID", "spID", "sex")) + names(ped) <- c("famID", "ID", "gen", "dadID", "momID", "spID", "sex") + expect_equal(result[result$ID != ID_drop, ], ped[ped$ID != ID_drop, ]) # are the families of the dropped ID in the original? @@ -242,6 +252,9 @@ test_that("dropLink - Drop specified by ID", { ped <- simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR) result <- dropLink(ped, ID_drop = ID_drop) + expect_equal(colnames(result), c("famID", "ID", "gen", "dadID", "momID", "spID", "sex")) + names(ped) <- c("famID", "ID", "gen", "dadID", "momID", "spID", "sex") + # are the dataframes the same in both the undropped and dropped relationships for all but the dropped ID? expect_equal(result[result$ID != ID_drop, ], ped[ped$ID != ID_drop, ]) @@ -262,7 +275,8 @@ test_that("dropLink - Drop specified by generation", { ped <- simulatePedigree(kpc = kpc, Ngen = Ngen, sexR = sexR, marR = marR) result <- dropLink(ped, gen_drop = gen_drop) - + expect_equal(colnames(result), c("famID", "ID", "gen", "dadID", "momID", "spID", "sex")) + names(ped) <- c("famID", "ID", "gen", "dadID", "momID", "spID", "sex") # are the dataframes the same in both the undropped and dropped relationships for all but the dropped gen? expect_equal(result[result$gen != gen_drop, ], ped[ped$gen != gen_drop, ]) @@ -399,3 +413,55 @@ test_that("addPersonToPed works as expected with zygosity", { expect_equal(updated3$personID[4], max(ped$personID, na.rm = TRUE) + 2) }) + +# Tests for single-ID specification (auto-find the other) + +test_that("makeTwins - specify only ID_twin1, auto-find twin2", { + ped <- data.frame( + famID = c(1, 1, 1, 1), + ID = c(1, 2, 3, 4), + gen = c(1, 1, 2, 2), + dadID = c(NA, NA, 1, 1), + momID = c(NA, NA, 2, 2), + spID = c(NA, NA, NA, NA), + sex = c("M", "F", "M", "F") + ) + # Person 3 (M) and 4 (F) are siblings. With DZ zygosity, either could be auto-selected. + result <- makeTwins(ped, ID_twin1 = 3, zygosity = "DZ") + expect_equal(sum(!is.na(result$twinID)), 2) + # Twin1 should be person 3 + expect_equal(result$twinID[result$ID == 3], 4) + expect_equal(result$twinID[result$ID == 4], 3) +}) + +test_that("makeTwins - specify only ID_twin2, auto-find twin1", { + ped <- data.frame( + famID = c(1, 1, 1, 1), + ID = c(1, 2, 3, 4), + gen = c(1, 1, 2, 2), + dadID = c(NA, NA, 1, 1), + momID = c(NA, NA, 2, 2), + spID = c(NA, NA, NA, NA), + sex = c("M", "F", "M", "F") + ) + result <- makeTwins(ped, ID_twin2 = 4, zygosity = "DZ") + expect_equal(sum(!is.na(result$twinID)), 2) + expect_equal(result$twinID[result$ID == 4], 3) + expect_equal(result$twinID[result$ID == 3], 4) +}) + +test_that("makeInbreeding - specify only ID_mate1, auto-find mate2", { + ped <- data.frame( + famID = c(1, 1, 1, 1), + ID = c(1, 2, 3, 4), + gen = c(1, 1, 2, 2), + dadID = c(NA, NA, 1, 1), + momID = c(NA, NA, 2, 2), + spID = c(NA, NA, NA, NA), + sex = c("M", "F", "M", "F") + ) + # Person 3 (M) should auto-find person 4 (F) as opposite-sex sibling + result <- makeInbreeding(ped, ID_mate1 = 3) + expect_equal(result$spID[result$ID == 3], 4) + expect_equal(result$spID[result$ID == 4], 3) +}) diff --git a/vignettes/v5_ASOIAF.Rmd b/vignettes/v5_ASOIAF.Rmd index d6c5a1ee..58c051f6 100644 --- a/vignettes/v5_ASOIAF.Rmd +++ b/vignettes/v5_ASOIAF.Rmd @@ -69,7 +69,8 @@ add <- ped2com(df_got, isChild_method = "partialparent", component = "additive", adjacency_method = "direct", - sparse = TRUE + sparse = TRUE, + mz_twins = TRUE ) mt <- ped2com(df_got, diff --git a/vignettes/v5_ASOIAF.html b/vignettes/v5_ASOIAF.html index 73674c37..933c0ac9 100644 --- a/vignettes/v5_ASOIAF.html +++ b/vignettes/v5_ASOIAF.html @@ -48,8 +48,9 @@ }