diff --git a/.Rbuildignore b/.Rbuildignore index 31a9a74..73b20e6 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -3,6 +3,8 @@ ^LICENSE\.md$ ^R/\.Rhistory$ ^R/paper$ +^checkParents\.X +^test-clean\.X CITATION.cff$ ^doc$ ^data-raw$ diff --git a/.github/workflows/R-CMD-dev_maincheck.yaml b/.github/workflows/R-CMD-dev_maincheck.yaml index 1db28e4..fe93bc4 100644 --- a/.github/workflows/R-CMD-dev_maincheck.yaml +++ b/.github/workflows/R-CMD-dev_maincheck.yaml @@ -33,22 +33,22 @@ jobs: steps: - uses: actions/checkout@v3 - - uses: r-lib/actions/setup-pandoc@v2 + - uses: r-lib/actions/setup-pandoc@v3 - - uses: r-lib/actions/setup-r@v2 + - uses: r-lib/actions/setup-r@v3 with: r-version: ${{ matrix.config.r }} rtools-version: ${{ matrix.config.rtools }} http-user-agent: ${{ matrix.config.http-user-agent }} use-public-rspm: true - - uses: r-lib/actions/setup-r-dependencies@v2 + - uses: r-lib/actions/setup-r-dependencies@v3 with: cache-version: 3 extra-packages: | any::rcmdcheck upgrade: 'TRUE' needs: check - - uses: r-lib/actions/check-r-package@v2 + - uses: r-lib/actions/check-r-package@v3 with: upload-snapshots: true diff --git a/.gitignore b/.gitignore index 9d9975f..3f8093b 100644 --- a/.gitignore +++ b/.gitignore @@ -4,13 +4,9 @@ R/plane.txt R/.Rhistory .Rhistory paper/paper.html -inst/doc -/doc/ /Meta/ *.knit.md vignettes/articles/paper.html BGmisc.code-workspace -.lintr - tests/testthat/Rplots.pdf -docs + diff --git a/.lintr b/.lintr new file mode 100644 index 0000000..10e49fe --- /dev/null +++ b/.lintr @@ -0,0 +1,2 @@ +linters: linters_with_defaults(line_length_linter(120),commented_code_linter = NULL,object_name_linter = object_name_linter(styles = c("snake_case", "symbols"))) # see vignette("lintr") +encoding: "UTF-8" diff --git a/CODE_OF_CONDUCT.md b/CODE_OF_CONDUCT.md index f6002f7..3e8d15a 100644 --- a/CODE_OF_CONDUCT.md +++ b/CODE_OF_CONDUCT.md @@ -60,7 +60,7 @@ representative at an online or offline event. Instances of abusive, harassing, or otherwise unacceptable behavior may be reported to the community leaders responsible for enforcement at -garrissm@wfu.edu. +[garrissm@wfu.edu](mailto:garrissm@wfu.edu). All complaints will be reviewed and investigated promptly and fairly. All community leaders are obligated to respect the privacy and security of the @@ -116,7 +116,7 @@ the community. This Code of Conduct is adapted from the [Contributor Covenant][homepage], version 2.0, available at -https://www.contributor-covenant.org/version/2/0/code_of_conduct.html. +. Community Impact Guidelines were inspired by [Mozilla's code of conduct enforcement ladder](https://github.com/mozilla/diversity). @@ -124,5 +124,5 @@ enforcement ladder](https://github.com/mozilla/diversity). [homepage]: https://www.contributor-covenant.org For answers to common questions about this code of conduct, see the FAQ at -https://www.contributor-covenant.org/faq. Translations are available at -https://www.contributor-covenant.org/translations. +. Translations are available at +. diff --git a/DESCRIPTION b/DESCRIPTION index a32c0b1..fb92558 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,17 +1,17 @@ Package: BGmisc Title: An R Package for Extended Behavior Genetics Analysis -Version: 1.2.0 +Version: 1.2.1 Authors@R: c( - person("S. Mason", "Garrison", , "garrissm@wfu.edu", role = c("aut", "cre"), + person("S. Mason", "Garrison", , email= "garrissm@wfu.edu", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-4804-6003")), person(c("Michael", "D."), "Hunter", role = "aut", comment = c(ORCID = "0000-0002-3651-6709")), person("Xuanyu", "Lyu", role = "aut", comment = c(ORCID = "0000-0002-2841-5529")), - person(c("Rachel", "N."), "Good", role = "ctb", - comment = c(ORCID = "0000-0000-0000-0000")), - person(c("Jonathan", "D."), "Trattner", role = "aut", - comment = c(ORCID = "0000-0002-1097-7603")), + person(c("Rachel", "N."), "Good", role = "ctb"), + person(c("Jonathan", "D."), "Trattner", role = "aut", email = "code@jdtrat.com", + comment = c(ORCID = "0000-0002-1097-7603", + url = "https://www.jdtrat.com/")), person(c("S.", "Alexandra"), "Burt", role = "aut", comment = c(ORCID = "0000-0001-5538-7431")) ) diff --git a/NEWS.md b/NEWS.md index 87892be..20f848d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,8 @@ +# BGmisc 1.2.1 + +* Added alternative transpose options for the matrix +* Added generalization of Falconer's formula + # BGmisc 1.2.0 * Added numerous code checks, increased code coverage to 85% diff --git a/R/checkParents.X b/R/checkParents.X new file mode 100644 index 0000000..d6cea52 --- /dev/null +++ b/R/checkParents.X @@ -0,0 +1,154 @@ +# Challenge: Missing parents: If one parent is missing and the other one isn't, this needs to be handled somehow. Firstly, I think it can cause certain ways of estimating relatedness to give wrong numbers. And secondly, it requires us to make guesses in cases where e.g. two people have the same mother and missing fathers: They could then be either half-sibs, if the missing fathers are different people, or full sibs if not. This then also affects relatedness for their descendants. + +#' Validates and Optionally Repairs Parent IDs in a Pedigree Dataframe +#' +#' This function takes a pedigree object and performs two main tasks: +#' 1. Checks for the validity of parent IDs, specifically looking for instances where only one parent ID is missing. +#' 2. Optionally repairs the missing parent IDs based on a specified logic. +#' +#' @param ped A dataframe representing the pedigree data with columns 'ID', 'dadID', and 'momID'. +#' @param verbose A logical flag indicating whether to print progress and validation messages to the console. +#' @param repair A logical flag indicating whether to attempt repairs on missing parent IDs. +#' +#' @return Depending on the value of `repair`, either a list containing validation results or a repaired dataframe is returned. +#' @examples +#' \dontrun{ +#' ped <- data.frame(ID = 1:4, dadID = c(NA, 1, 1, 2), momID = c(NA, NA, 2, 2)) +#' checkParentIDs(ped, verbose = TRUE, repair = FALSE) +#' } +#' @export +checkParentIDs <- function(ped, verbose = FALSE, repair = FALSE) { + # Standardize column names in the input dataframe + ped <- standardize_colnames(ped) + + # Initialize a list to store validation results + validation_results <- list() + + if (verbose) { + cat("Step 1: Checking for missing parents...\n") + } + + # Identify missing fathers and mothers + missing_fathers <- ped$ID[which(is.na(ped$dadID) & !is.na(ped$momID))] + missing_mothers <- ped$ID[which(!is.na(ped$dadID) & is.na(ped$momID))] + + # Update the validation_results list + if (length(missing_fathers) > 0) { + validation_results$missing_fathers <- missing_fathers + } + if (length(missing_mothers) > 0) { + validation_results$missing_mothers <- missing_mothers + } + + # If no missing parents are found + if (length(validation_results) == 0) { + if (verbose) { + cat("No missing single parents found.\n") + } + validation_results$missing_parents <- FALSE + } + # If missing parents are found + else { + if (verbose) { + cat("Missing single parents found.\n") + } + validation_results$missing_parents <- TRUE + } + + if (verbose) { + cat("Step 2: Determining the if moms are the same sex and dads are same sex\n") + } + + # Determine the most frequent sex for moms and dads + most_frequent_sex_mom <- names(sort(table(ped$sex[ped$ID %in% ped$momID]), decreasing = TRUE))[1] + most_frequent_sex_dad <- names(sort(table(ped$sex[ped$ID %in% ped$dadID]), decreasing = TRUE))[1] + + # are all moms/dads the same sex? + validation_results$mom_sex <- unique(ped$sex[ped$ID %in% ped$momID]) + validation_results$dad_sex <- unique(ped$sex[ped$ID %in% ped$dadID]) + + # Store the most frequent sex for moms and dads + if (is.numeric(ped$sex)) { + validation_results$female_var <- as.numeric(most_frequent_sex_mom) + validation_results$male_var <- as.numeric(most_frequent_sex_dad) + } else if (is.character(ped$sex) | is.factor(ped$sex)) { + validation_results$female_var <- most_frequent_sex_mom + validation_results$male_var <- most_frequent_sex_dad + } else { + print("You should never see this. If you do, then you have a problem with the data type of the sex variable") + } + + # verbose + if (length(validation_results$mom_sex) == 1) { + if (verbose) { + cat(paste0( + "All moms are '", + validation_results$female_var, + "'.\n" + )) + } + validation_results$female_moms <- TRUE + } else { + validation_results$female_moms <- FALSE + } + + if (length(validation_results$dad_sex) == 1) { + if (verbose) { + cat(paste0( + "All dads are '", + validation_results$male_var, + "'.\n" + )) + } + validation_results$male_dads <- TRUE + } else { + validation_results$male_dads <- FALSE + } + # Check for inconsistent gender roles + wrong_sex_moms <- ped$ID[which(ped$sex[ped$ID %in% ped$momID] != validation_results$female_var)] + wrong_sex_dads <- ped$ID[which(ped$sex[ped$ID %in% ped$dadID] != validation_results$male_var)] + + + + # Are any parents in both momID and dadID? + momdad <- intersect(ped$dadID, ped$momID) + if (length(momdad) > 0) { + validation_results$parents_in_both <- momdad + } + + + male_moms <- ped$ID[which(ped$dadID & !is.na(ped$momID))] + + if (repair) { + if (verbose) { + cat("Validation Results:\n") + print(validation_results) + cat("Step 2: Attempting to repair missing parents...\n") + } + cat("REPAIR IN EARLY ALPHA\n") + # Initialize a list to track changes made during repair + changes <- list() + + # [Insert logic to repair parent IDs here] + + # Update the pedigree dataframe after repair + repaired_ped <- ped + + if (verbose) { + cat("Changes Made:\n") + print(changes) + } + return(repaired_ped) + } else { + return(validation_results) + } +} +#' Repair Parent IDs +#' +#' This function repairs parent IDs in a pedigree. +#' @param ped A pedigree object +#' @param verbose A logical indicating whether to print progress messages +#' @return A corrected pedigree +repairParentIDs <- function(ped, verbose = FALSE) { + checkParentIDs(ped = ped, verbose = verbose, repair = TRUE) +} diff --git a/R/convertPedigree.R b/R/convertPedigree.R index 428b282..4123c4c 100644 --- a/R/convertPedigree.R +++ b/R/convertPedigree.R @@ -9,6 +9,8 @@ #' @param gc logical. If TRUE, do frequent garbage collection via \code{\link{gc}} to save memory #' @param flatten.diag logical. If TRUE, overwrite the diagonal of the final relatedness matrix with ones #' @param standardize.colnames logical. If TRUE, standardize the column names of the pedigree dataset +#' @param tcross.alt.crossprod logical. If TRUE, use alternative method of using Crossprod function for computing the transpose +#' @param tcross.alt.star logical. If TRUE, use alternative method of using \%\*\% for computing the transpose #' @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". #' @export @@ -20,6 +22,8 @@ ped2com <- function(ped, component, gc = FALSE, flatten.diag = FALSE, standardize.colnames = TRUE, + tcross.alt.crossprod = FALSE, + tcross.alt.star = FALSE, ...) { # Validate the 'component' argument and match it against predefined choices component <- match.arg(tolower(component), @@ -169,7 +173,15 @@ ped2com <- function(ped, component, if (verbose) { cat("Doing tcrossprod\n") } - r <- Matrix::tcrossprod(r2) + if(tcross.alt.crossprod){ + cat("Doing alt tcrossprod crossprod t \n") + r <- crossprod(t(as.matrix(r2))) + }else if(tcross.alt.star){ + cat("Doing alt tcrossprod %*% t \n") + r <- r2 %*% t(as.matrix(r2)) + }else{ + r <- Matrix::tcrossprod(r2) + } if (component == "generation") { return(gen) } else { @@ -194,7 +206,8 @@ ped2com <- function(ped, component, #' @export #' -ped2add <- function(ped, max.gen = Inf, sparse = FALSE, verbose = FALSE, gc = FALSE, flatten.diag = FALSE) { +ped2add <- function(ped, max.gen = Inf, sparse = FALSE, verbose = FALSE, gc = FALSE, flatten.diag = FALSE, standardize.colnames = TRUE, + tcross.alt.crossprod = FALSE, tcross.alt.star = FALSE) { ped2com( ped = ped, max.gen = max.gen, @@ -202,7 +215,10 @@ ped2add <- function(ped, max.gen = Inf, sparse = FALSE, verbose = FALSE, gc = FA verbose = verbose, gc = gc, component = "additive", - flatten.diag = flatten.diag + flatten.diag = flatten.diag, + standardize.colnames = standardize.colnames, + tcross.alt.crossprod = tcross.alt.crossprod, + tcross.alt.star = tcross.alt.star ) } @@ -212,7 +228,7 @@ ped2add <- function(ped, max.gen = Inf, sparse = FALSE, verbose = FALSE, gc = FA #' @export #' @aliases ped2mt #' -ped2mit <- ped2mt <- function(ped, max.gen = Inf, sparse = FALSE, verbose = FALSE, gc = FALSE, flatten.diag = FALSE) { +ped2mit <- ped2mt <- function(ped, max.gen = Inf, sparse = FALSE, verbose = FALSE, gc = FALSE, flatten.diag = FALSE, standardize.colnames = TRUE, tcross.alt.crossprod = FALSE, tcross.alt.star = FALSE) { ped2com( ped = ped, max.gen = max.gen, @@ -220,19 +236,19 @@ ped2mit <- ped2mt <- function(ped, max.gen = Inf, sparse = FALSE, verbose = FALS verbose = verbose, gc = gc, component = "mitochondrial", - flatten.diag = flatten.diag + flatten.diag = flatten.diag, + standardize.colnames = standardize.colnames, + tcross.alt.crossprod = tcross.alt.crossprod, + tcross.alt.star = tcross.alt.star ) } - - - #' Take a pedigree and turn it into a common nuclear environmental relatedness matrix #' @inheritParams ped2com #' @details The algorithms and methodologies used in this function are further discussed and exemplified in the vignette titled "examplePedigreeFunctions". #' @export #' -ped2cn <- function(ped, max.gen = Inf, sparse = FALSE, verbose = FALSE, gc = FALSE, flatten.diag = FALSE) { +ped2cn <- function(ped, max.gen = Inf, sparse = FALSE, verbose = FALSE, gc = FALSE, flatten.diag = FALSE, standardize.colnames = TRUE, tcross.alt.crossprod = FALSE, tcross.alt.star = FALSE) { ped2com( ped = ped, max.gen = max.gen, @@ -240,7 +256,10 @@ ped2cn <- function(ped, max.gen = Inf, sparse = FALSE, verbose = FALSE, gc = FAL verbose = verbose, gc = gc, component = "common nuclear", - flatten.diag = flatten.diag + flatten.diag = flatten.diag, + standardize.colnames = standardize.colnames, + tcross.alt.crossprod = tcross.alt.crossprod, + tcross.alt.star = tcross.alt.star ) } diff --git a/R/formula.R b/R/formula.R index 03c897f..cae0f88 100644 --- a/R/formula.R +++ b/R/formula.R @@ -87,3 +87,45 @@ inferRelatedness <- function(obsR, aceA = .9, aceC = 0, sharedC = 0) { calc_r <- (obsR - sharedC * aceC) / aceA return(calc_r) } + +#' Falconer's Formula +#' Use Falconer's formula to solve for H using the observed correlations for two groups of any two relatednesses. +#' @param r1 Relatedness coefficient of the first group. +#' @param r2 Relatedness coefficient of the second group. +#' @param obsR1 Observed correlation between members of the first group. +#' @param obsR2 Observed correlation between members of the second group. +#' +#' @return Heritability estimates (`heritability_estimates`). + +calculateH <- function(r1, r2, obsR1, obsR2) { + # Check for equal relatedness coefficients to avoid division by zero + if (any(r1 - r2 == 0)) { + stop("Relatedness coefficients r1 and r2 must not be equal for any pair.") + } + if (any(abs(obsR1) > 1) || any(abs(obsR2) > 1)) { + warning("The observed correlations should be between -1 and 1.") + } + + if (any(obsR1 * obsR2 < 0)) { + warning("The correlations should not have opposite signs.") + } + + if(any(obsR1 < 0 & obsR2 < 0)){ + message("Your scale might be reverse coded because you have negative correlations. Please check your data. ") + } + + + # Calculate heritability estimates (H^2) for all pairs + heritability_estimates <- (obsR1 - obsR2) / (r1 - r2) + # Check for unrealistic heritability estimates and warn the user + if (any(heritability_estimates < 0)) { + warning("Some calculated heritability values are negative, which may indicate assumption violations or questions about directionality.") + } + + if (any(heritability_estimates > 1)) { + warning("Some calculated heritability values are greater than 1, which may suggest overestimation or errors in the observed correlations or relatedness coefficients.") + } + + + return(heritability_estimates) +} diff --git a/README.md b/README.md index 11fd019..95ac055 100644 --- a/README.md +++ b/README.md @@ -13,6 +13,7 @@ developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.re version](https://www.r-pkg.org/badges/version/BGmisc)](https://cran.r-project.org/package=BGmisc) [![Package downloads](https://cranlogs.r-pkg.org/badges/grand-total/BGmisc)](https://cran.r-project.org/package=BGmisc)
+ [![R-CMD-check](https://github.com/R-Computing-Lab/BGmisc/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/R-Computing-Lab/BGmisc/actions/workflows/R-CMD-check.yaml) [![Dev Main branch](https://github.com/R-Computing-Lab/BGmisc/actions/workflows/R-CMD-dev_maincheck.yaml/badge.svg)](https://github.com/R-Computing-Lab/BGmisc/actions/workflows/R-CMD-dev_maincheck.yaml) diff --git a/cran-comments.md b/cran-comments.md index 8196af5..0314d3b 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,12 +1,9 @@ -Description ------------------------------------------------ +# Description This update reflects a substantial improvement in the codebase as part of the peer review process for JOSS, including the addition of numerous function checks, increased code coverage to 85%, and the replacement of sapply usage. We also added a new function to simulate twins and the ability to trace paternal and maternal lines. We also added a Harry Potter pedigree. - -Test Environments ------------------------------------------------ +# Test Environments 1. Local OS: Windows 11 x64 (build 22635), R version 4.3.2 (2023-10-31 ucrt) 2. **GitHub Actions**: diff --git a/man/BGmisc-package.Rd b/man/BGmisc-package.Rd index 963688f..55e1b84 100644 --- a/man/BGmisc-package.Rd +++ b/man/BGmisc-package.Rd @@ -24,10 +24,14 @@ Authors: \itemize{ \item Michael D. Hunter (\href{https://orcid.org/0000-0002-3651-6709}{ORCID}) \item Xuanyu Lyu (\href{https://orcid.org/0000-0002-2841-5529}{ORCID}) - \item Rachel N. Good (\href{https://orcid.org/0000-0000-0000-0000}{ORCID}) - \item Jonathan D. Trattner (\href{https://orcid.org/0000-0002-1097-7603}{ORCID}) + \item Jonathan D. Trattner \email{code@jdtrat.com} (\href{https://orcid.org/0000-0002-1097-7603}{ORCID}) (https://www.jdtrat.com/) \item S. Alexandra Burt (\href{https://orcid.org/0000-0001-5538-7431}{ORCID}) } +Other contributors: +\itemize{ + \item Rachel N. Good [contributor] +} + } \keyword{internal} diff --git a/man/calculateH.Rd b/man/calculateH.Rd new file mode 100644 index 0000000..9fe2415 --- /dev/null +++ b/man/calculateH.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/formula.R +\name{calculateH} +\alias{calculateH} +\title{Falconer's Formula +Use Falconer's formula to solve for H using the observed correlations for two groups of any two relatednesses.} +\usage{ +calculateH(r1, r2, obsR1, obsR2) +} +\arguments{ +\item{r1}{Relatedness coefficient of the first group.} + +\item{r2}{Relatedness coefficient of the second group.} + +\item{obsR1}{Observed correlation between members of the first group.} + +\item{obsR2}{Observed correlation between members of the second group.} +} +\value{ +Heritability estimates (`heritability_estimates`). +} +\description{ +Falconer's Formula +Use Falconer's formula to solve for H using the observed correlations for two groups of any two relatednesses. +} diff --git a/man/ped2add.Rd b/man/ped2add.Rd index ee6e054..ad87012 100644 --- a/man/ped2add.Rd +++ b/man/ped2add.Rd @@ -10,7 +10,10 @@ ped2add( sparse = FALSE, verbose = FALSE, gc = FALSE, - flatten.diag = FALSE + flatten.diag = FALSE, + standardize.colnames = TRUE, + tcross.alt.crossprod = FALSE, + tcross.alt.star = FALSE ) } \arguments{ @@ -27,6 +30,12 @@ generations as there are in the data.} \item{gc}{logical. If TRUE, do frequent garbage collection via \code{\link{gc}} to save memory} \item{flatten.diag}{logical. If TRUE, overwrite the diagonal of the final relatedness matrix with ones} + +\item{standardize.colnames}{logical. If TRUE, standardize the column names of the pedigree dataset} + +\item{tcross.alt.crossprod}{logical. If TRUE, use alternative method of using Crossprod function for computing the transpose} + +\item{tcross.alt.star}{logical. If TRUE, use alternative method of using \%\*\% for computing the transpose} } \description{ Take a pedigree and turn it into an additive genetics relatedness matrix diff --git a/man/ped2cn.Rd b/man/ped2cn.Rd index af35226..0d7e3eb 100644 --- a/man/ped2cn.Rd +++ b/man/ped2cn.Rd @@ -10,7 +10,10 @@ ped2cn( sparse = FALSE, verbose = FALSE, gc = FALSE, - flatten.diag = FALSE + flatten.diag = FALSE, + standardize.colnames = TRUE, + tcross.alt.crossprod = FALSE, + tcross.alt.star = FALSE ) } \arguments{ @@ -27,6 +30,12 @@ generations as there are in the data.} \item{gc}{logical. If TRUE, do frequent garbage collection via \code{\link{gc}} to save memory} \item{flatten.diag}{logical. If TRUE, overwrite the diagonal of the final relatedness matrix with ones} + +\item{standardize.colnames}{logical. If TRUE, standardize the column names of the pedigree dataset} + +\item{tcross.alt.crossprod}{logical. If TRUE, use alternative method of using Crossprod function for computing the transpose} + +\item{tcross.alt.star}{logical. If TRUE, use alternative method of using \%\*\% for computing the transpose} } \description{ Take a pedigree and turn it into a common nuclear environmental relatedness matrix diff --git a/man/ped2com.Rd b/man/ped2com.Rd index 7c4c5f0..67daf5f 100644 --- a/man/ped2com.Rd +++ b/man/ped2com.Rd @@ -13,6 +13,8 @@ ped2com( gc = FALSE, flatten.diag = FALSE, standardize.colnames = TRUE, + tcross.alt.crossprod = FALSE, + tcross.alt.star = FALSE, ... ) } @@ -35,6 +37,10 @@ generations as there are in the data.} \item{standardize.colnames}{logical. If TRUE, standardize the column names of the pedigree dataset} +\item{tcross.alt.crossprod}{logical. If TRUE, use alternative method of using Crossprod function for computing the transpose} + +\item{tcross.alt.star}{logical. If TRUE, use alternative method of using \%\*\% for computing the transpose} + \item{...}{additional arguments to be passed to \code{\link{ped2com}}} } \description{ diff --git a/man/ped2mit.Rd b/man/ped2mit.Rd index af6f739..530ca4e 100644 --- a/man/ped2mit.Rd +++ b/man/ped2mit.Rd @@ -11,7 +11,10 @@ ped2mit( sparse = FALSE, verbose = FALSE, gc = FALSE, - flatten.diag = FALSE + flatten.diag = FALSE, + standardize.colnames = TRUE, + tcross.alt.crossprod = FALSE, + tcross.alt.star = FALSE ) } \arguments{ @@ -28,6 +31,12 @@ generations as there are in the data.} \item{gc}{logical. If TRUE, do frequent garbage collection via \code{\link{gc}} to save memory} \item{flatten.diag}{logical. If TRUE, overwrite the diagonal of the final relatedness matrix with ones} + +\item{standardize.colnames}{logical. If TRUE, standardize the column names of the pedigree dataset} + +\item{tcross.alt.crossprod}{logical. If TRUE, use alternative method of using Crossprod function for computing the transpose} + +\item{tcross.alt.star}{logical. If TRUE, use alternative method of using \%\*\% for computing the transpose} } \description{ Take a pedigree and turn it into a mitochondrial relatedness matrix diff --git a/tests/testthat/test-clean.X b/tests/testthat/test-clean.X new file mode 100644 index 0000000..ef8aebe --- /dev/null +++ b/tests/testthat/test-clean.X @@ -0,0 +1,15 @@ +# Test that cleaning functions work +test_that("Check Sex", { + data(potter) + df_fam <- potter + df_fam$sex[df_fam$name=="Vernon Dursley"] <- 0 + + result <- checkSex( + df_fam, + code_male = 1, + verbose = FALSE, + repair = TRUE + ) + expect_equal(result[[2]]$sex[result[[2]][["name"]]=="Vernon Dursley"], "M", tolerance = 1e-8) +}) + diff --git a/tests/testthat/test-network.R b/tests/testthat/test-network.R index d5ecbfe..656dfd4 100755 --- a/tests/testthat/test-network.R +++ b/tests/testthat/test-network.R @@ -52,6 +52,25 @@ test_that("ped2add produces correct matrix dims, values, and dimnames for hazard expect_equal(dn[[1]], as.character(hazard$ID)) }) +test_that("ped2add produces correct matrix dims, values, and dimnames for alternative transpose", { + data(hazard) + add <- ped2add(hazard, tcross.alt.crossprod = TRUE) + # Check dimension + expect_equal(dim(add), c(nrow(hazard), nrow(hazard))) + # Check several values + expect_true(all(diag(add) == 1)) + expect_equal(add, t(add)) + expect_equal(add[2, 1], 0) + expect_equal(add[10, 1], .25) + expect_equal(add[9, 1], 0) + expect_equal(add["5", "6"], .5) + # Check that dimnames are correct + dn <- dimnames(add) + expect_equal(dn[[1]], dn[[2]]) + expect_equal(dn[[1]], as.character(hazard$ID)) +}) +# to do, combine the sets that are equalivant. shouldn't need to run 1000 expect equals + test_that("ped2add produces correct matrix dims, values, and dimnames for inbreeding data", { data(inbreeding) add <- ped2add(inbreeding) @@ -69,6 +88,25 @@ test_that("ped2add produces correct matrix dims, values, and dimnames for inbree expect_equal(dn[[1]], dn[[2]]) expect_equal(dn[[1]], as.character(inbreeding$ID)) }) + + +test_that("ped2add produces correct matrix dims, values, and dimnames for inbreeding data with alternative transpose", { + data(inbreeding) + add <- ped2add(inbreeding, tcross.alt.star = TRUE) + # Check dimension + expect_equal(dim(add), c(nrow(inbreeding), nrow(inbreeding))) + # Check several values + expect_true(all(diag(add) >= 1)) + expect_equal(add, t(add)) + expect_equal(add[2, 1], 0) + expect_equal(add[6, 1], .5) + expect_equal(add[113, 113], 1.1250) + expect_equal(add["113", "112"], 0.62500) + # Check that dimnames are correct + dn <- dimnames(add) + expect_equal(dn[[1]], dn[[2]]) + expect_equal(dn[[1]], as.character(inbreeding$ID)) +}) test_that("ped2add flattens diagonal for inbreeding data", { data(inbreeding) add <- ped2add(inbreeding, flatten.diag = TRUE) diff --git a/tests/testthat/test-relatedness-functions.R b/tests/testthat/test-relatedness-functions.R index a2f4108..387c5bc 100644 --- a/tests/testthat/test-relatedness-functions.R +++ b/tests/testthat/test-relatedness-functions.R @@ -33,3 +33,60 @@ test_that("inferRelatedness performs as expected", { "aceA and aceC must be proportions between 0 and 1" ) }) + + +# Test 1: Basic Functionality Test +test_that("calculateH returns correct heritability estimates", { + expect_equal(calculateH(0.5, 0.25, 0.4, 0.2), 0.8) + expect_equal(calculateH(0.5, 0.125, 0.5, 0.25), 2/3, tolerance = 1e-8) +}) + +# Test 2: unusual warning +test_that("provides warning for unusual H value", { + r1 <- 0.25 + r2 <- 0.125 + obsR1 <- 0.3 + obsR2 <- 0.1 + expected <- 1.6 + expect_warning(expect_equal(calculateH(r1, r2, obsR1, obsR2), expected)) +}) + + +# Test 3: Vectorized Input Test +test_that("calculateH handles vectorized inputs correctly", { + r1 <- c(0.5, 0.5) + r2 <- c(0.25, 0.125) + obsR1 <- c(0.4, 0.5) + obsR2 <- c(0.2, 0.25) + expected <- c(0.8, 2/3) + expect_equal(calculateH(r1, r2, obsR1, obsR2), expected, tolerance = 1e-8) +}) + + +# Test 4: Equal Relatedness Coefficients Test +test_that("calculateH stops for equal relatedness coefficients", { + expect_error(calculateH(0.5, 0.5, 0.4, 0.2), + "Relatedness coefficients r1 and r2 must not be equal for any pair.") +}) + +# Test 5: Negative and Positive Correlation Test +test_that("calculateH handles both negative and positive correlations", { + # Test for negative correlations with expected warnings about negative heritability values + expect_warning( + expect_equal(calculateH(0.5, 0.25, -0.4, -0.2), -0.8), + regexp = "Some calculated heritability values are negative" + ) + # Test for a scenario leading to a positive heritability estimate + + expect_warning( + expect_warning( + expect_equal(calculateH(0.5, 0.25, 0.2, -0.1), 1.2), + regexp = "The correlations should not have opposite signs."), + regexp = "Some calculated heritability values are greater than 1") +}) + +# Test 6: illegal correlation values +test_that("calculateH stops for illegal coefficients", { + expect_warning(calculateH(0.5, 0.25, 1.4, 1.4), + "The observed correlations should be between -1 and 1") +}) diff --git a/vignettes/network.Rmd b/vignettes/network.Rmd index 3b5d0c0..375b727 100644 --- a/vignettes/network.Rmd +++ b/vignettes/network.Rmd @@ -181,6 +181,6 @@ However, this subset does not plot the relationship between spouses (such as the subset_rows <- c(1:5, 31:36) subset_potter <- potter[subset_rows, ] ``` -```{r, echo=FALSE, results='hide', out.width='50%', fig.cap="Potter Subset Pedigree"} +```{r, echo=FALSE, results='hide', out.width='50%'} plotPedigree(subset_potter, code_male = 1, verbose = TRUE) ``` diff --git a/vignettes/network.html b/vignettes/network.html index 5c7ed08..07fd015 100644 --- a/vignettes/network.html +++ b/vignettes/network.html @@ -529,12 +529,7 @@

Subsetting Pedigrees

are not children to connect the two individuals together yet.

subset_rows <- c(1:5, 31:36)
 subset_potter <- potter[subset_rows, ]
-
-Potter Subset Pedigree -

-Potter Subset Pedigree -

-
+

diff --git a/vignettes/pedigree.Rmd b/vignettes/pedigree.Rmd index aa029e1..3548c17 100644 --- a/vignettes/pedigree.Rmd +++ b/vignettes/pedigree.Rmd @@ -49,7 +49,6 @@ The columns represents the individual's family ID, the individual's personal ID, ## Plotting Pedigree - Pedigrees are visual diagrams that represent family relationships across generations. They are commonly used in genetics to trace the inheritance of specific traits or conditions. This vignette will guide you through visualizing simulated pedigrees using the `plotPedigree` function. This function is a wrapper function for `Kinship2`'s base R plotting. ### Single Pedigree Visualization