[](https://joss.theoj.org/papers/10.21105/joss.06203)
[](https://www.repostatus.org/#active)
[](https://cran.r-project.org/package=BGmisc)
diff --git a/README.md b/README.md
index 6b8357fe..eede5556 100644
--- a/README.md
+++ b/README.md
@@ -5,16 +5,20 @@
+
[](https://joss.theoj.org/papers/10.21105/joss.06203)
[](https://www.repostatus.org/#active)
-[](https://cran.r-project.org/package=BGmisc)
+[](https://cran.r-project.org/package=BGmisc)
[](https://cran.r-project.org/package=BGmisc)
[](https://github.com/R-Computing-Lab/BGmisc/actions/workflows/R-CMD-check.yaml)
-[](https://github.com/R-Computing-Lab/BGmisc/actions/workflows/R-CMD-dev_maincheck.yaml)
-[](https://app.codecov.io/gh/R-Computing-Lab/BGmisc)
+[](https://github.com/R-Computing-Lab/BGmisc/actions/workflows/R-CMD-dev_maincheck.yaml)
+[](https://app.codecov.io/gh/R-Computing-Lab/BGmisc)
diff --git a/data-raw/HouseofHabsburg.ged b/data-raw/HouseofHabsburg.ged
index 98b7100c..fd358bf8 100644
--- a/data-raw/HouseofHabsburg.ged
+++ b/data-raw/HouseofHabsburg.ged
@@ -1,4 +1,4 @@
-0 HEAD
+0 HEAD
1 SOUR PAF
2 NAME Personal Ancestral File
2 VERS 5.2.18.0
@@ -260,7 +260,7 @@
0 @I17@ INDI
1 NAME /Frederick I/
2 SURN Frederick I
-2 _AKA Frederick the Handsome or the Fair
+2 _AKA Frederick the Handsome or the Fair
1 SEX M
1 BIRT
2 DATE 1289
diff --git a/data-raw/Targaryen tree Dance.txt b/data-raw/Targaryen tree Dance.txt
new file mode 100644
index 00000000..082c263c
--- /dev/null
+++ b/data-raw/Targaryen tree Dance.txt
@@ -0,0 +1,17 @@
+{{familytree/start|style=font: 90% sans-serif;|state={{{state|}}}|title={{{title|}}}}}
+{{familytree | | | | | | | | | | | | | JAE |v| ALY | |JAE=[[Jaehaerys I Targaryen|Jaehaerys I]]|ALY=[[Alysanne Targaryen|Alysanne]]|}}
+{{familytree | | | | | | | |,|-|-|-|v|-|-|-|^|-|-|-|v|-|-|-|.| | | | | | | | | | | | | | | | | | | | | }}
+{{familytree | | |JOC|v|AEM| |DAE|v|ROD| |BAE|v|ALY| |JOC=[[Jocelyn Baratheon|JocelynJust how related are Jon Snow and Daenerys Targaryen? This vignette
-walks through how to quantify their genetic relatedness using functions
-from the BGmisc package. While the Game of Thrones canon gives us some
-clues, we can use a formal pedigree-based approach to quantify their
-genetic relatedness. This vignette demonstrates how to compute
-coefficients of relatedness using the BGmisc package, along
-with basic data manipulation from tidyverse. We will also handle
-incomplete parental information programmatically and generate a plot of
-the reconstructed pedigree.
Just how closely related are Jon Snow and Daenerys Targaryen? +According to the lore of A Song of Ice and Fire, Daenerys is +Jon’s paternal aunt. This would suggest a theoretical genetic +relatedness of 0.25, assuming a simple pedigree and no inbreeding. But +with tangled ancestries and potentially missing information, how +confident can we be in that estimate?
+In this vignette, we use the BGmisc package to
+reconstruct the ASOIAF pedigree, handle incomplete parentage
+data, and compute additive genetic and common nuclear relatedness. We’ll
+focus on Jon and Daenerys as a case study, but the methods generalize to
+any characters in the provided dataset.
We begin by loading the necessary packages and accessing the built-in
-ASOIAF pedigree dataset included with
-BGmisc.
We begin by loading the required libraries and examining the
+structure of the built-in ASOIAF pedigree.
The ASOIAF data contains character IDs, family identifiers, and -parent-child links extracted from A Song of Ice and Fire lore.
+The ASOIAF dataset includes character IDs, names, family identifiers, +and parent identifiers for a subset of characters drawn from the A +Song of Ice and Fire canon.
## id famID momID dadID name sex
## 1 1 1 NA NA Walder Frey M
@@ -370,9 +375,10 @@ Load Packages and Data
We use checkSex() to ensure that all individuals have
-valid sex codes, repairing as needed. This is important for correct
-pedigree plotting and downstream calculations.
Many pedigree-based algorithms rely on biological sex for downstream
+calculationss and visualization. We use checkSex() to
+inspect the sex variable, repairing inconsistencies
+programmatically.
We convert the component matrices into a long-format table of
-pairwise relationships using com2links(). This gives us a
-long dataframe where each row represents a pair of individuals and their
-relatedness. The function can return the entire matrix or just the lower
-triangular part, which is often sufficient for our purposes. We set
-writetodisk = FALSE to keep the data in memory.
For interpretability, we convert these square matrices into
+long-format tables using com2links(). This function returns
+a dataframe where each row represents a unique pair of individuals,
+including their additive and common nuclear coefficients.
df_links <- com2links(
writetodisk = FALSE,
- ad_ped_matrix = add, cn_ped_matrix = cn,
+ ad_ped_matrix = add, cn_ped_matrix = cn, mit_ped_matrix = mt,
drop_upper_triangular = TRUE
-)# %>%
+) # %>%
# filter(ID1 != ID2)The function can return the entire matrix or just the lower
+triangular part, which is often sufficient for our purposes. Setting
+drop_upper_triangular = TRUE ensures we only retain one
+entry per dyad, since the matrices are symmetric. We also keep the data
+in memory by setting writetodisk = FALSE.
Next, we extract the IDs corresponding to Jon Snow and Daenerys
-Targaryen. We use the filter() function to find the rows in
-the df_links dataframe where either ID1 or ID2 corresponds
-to Jon Snow, and then filter again to find the row where the other ID
-corresponds to Daenerys Targaryen.
We next identify the rows in the pairwise relatedness table that +correspond to Jon Snow and Daenerys Targaryen. First, we retrieve their +individual IDs:
# Find the IDs of Jon Snow and Daenerys Targaryen
jon_id <- df_got %>%
@@ -435,35 +461,48 @@ Locate Jon and Daenerys
dany_id <- df_got %>%
filter(name == "Daenerys Targaryen") %>%
pull(ID)We then filter the pairwise table to retrieve the row containing -their relationship.
+Then we isolate their dyad:
jon_dany_row <- df_links %>%
filter(ID1 == jon_id | ID2 == jon_id) %>%
- filter(ID1 %in% dany_id| ID2 %in% dany_id)
+ filter(ID1 %in% dany_id | ID2 %in% dany_id)
-jon_dany_row ## ID1 ID2 addRel cnuRel
-## 1 206 211 0.31274414 0
-## 2 211 304 0.01953125 0This row contains the additive relatedness coefficient between Jon -and Daenerys, which allows us to assess how closely related they are -genetically. We’d expect to see a value of 0.25 for an Aunt-Nephew -relationship, which is what Jon and Daenerys are in the show. However, -the value is 0.3127441, indicating a more complex relationship.
+jon_dany_row## ID1 ID2 addRel mitRel cnuRel
+## 1 206 211 0.31274414 0 0
+## 2 211 304 0.01953125 0 0This table contains the additive and nuclear relatedness estimates +for Jon and Daenerys. If the pedigree reflects their canonical +aunt-nephew relationship and is free from inbreeding, we’d expect to see +an additive coefficient close to 0.25. However, the value is 0.3127441, +indicating a more complex relationship.
Many real-world and fictional pedigrees contain individuals with
+unknown or partially known parentage. In such cases, plotting tools
+typically fail unless these gaps are handled. We use
+checkParentIDs() to:
Identify individuals with one known parent and one +missing
Create “phantom” placeholders for the missing parent
-Optionally repair and harmonize parent fields
To facilitate plotting, we check for individuals with one known parent but a missing other. For those cases, we assign a placeholder ID to the missing parent.
-df_repaired <- checkParentIDs(df_got,addphantoms=TRUE,
- repair=TRUE,
- parentswithoutrow=FALSE,
- repairsex=FALSE
- ) %>% mutate(fam=1,
- affected = case_when(ID %in% c(jon_id,dany_id) ~ 1,
- TRUE ~ 0)
- )df_repaired <- checkParentIDs(df_got,
+ addphantoms = TRUE,
+ repair = TRUE,
+ parentswithoutrow = FALSE,
+ repairsex = FALSE
+) %>% mutate(
+ famID = 1,
+ affected = case_when(
+ ID %in% c(jon_id, dany_id, "365") ~ 1,
+ TRUE ~ 0
+ )
+)## REPAIR IN EARLY ALPHAThis code creates new IDs for individuals with one known parent and a
missing other. It checks if either momID or
@@ -473,9 +512,13 @@
#fixParents(id=df_got$ID, dadid=df_got$dadID, momid=df_got$momID, sex=df_got$sex, missid = NA)
-
-plotPedigree(df_repaired,affected=df_repaired$affected,verbose=FALSE)We can now visualize the repaired pedigree using the
+plotPedigree() function. This function generates a plot of
+the pedigree, with individuals colored based on their affected status.
+In this case, we highlight Jon and Daenerys as “affected” individuals.
+Otherwise they would be difficult to distinguish from the rest of the
+pedigree.
## Did not plot the following people: 85 88 125 142 228 229 258 259 274 275 305 336 357 381 388 405 409 418 420 424 428 451 487## named list()This vignette provides a detailed guide to specific functions within
+the BGmisc package that aid in the identification and
+fitting of variance component models common in behavior genetics. We
+will explore key functions such as identifyComponentModel,
+providing practical examples and theoretical background. Identification
+ensures a unique set of parameters that define the model-implied
+covariance matrix, preventing free parameters from trading off one
+another.
Ensure that the BGmisc package is installed and
+loaded.
Ensure that the following dependencies are installed before +proceeding as they provide us with behavior genetic data and models:
+EasyMx
OpenMx
Note: If any of the libraries are not installed, you can install them
+using install.packages(“package_name”).
In this section, we will demonstrate core functions related to the +identification and fitting of variance component models.
+comp2vech FunctionThe comp2vech function is used to vectorize a components
+model. The function is often used in conjunction with the identification
+process. In this example, we apply it to a list of matrices:
comp2vech(list(
+ matrix(c(1, .5, .5, 1), 2, 2),
+ matrix(1, 2, 2)
+))
+#> [1] 1.0 0.5 1.0 1.0 1.0 1.0The result showcases how the matrices have been transformed, +reflecting their role in subsequent variance component analysis.
+identifyComponentModel FunctionThe identifyComponentModel function helps determine if a
+variance components model is identified. It accepts relatedness
+component matrices and returns information about identified and
+non-identified parameters.
Here’s an example using the classical twin model with only MZ +twins:
+identifyComponentModel(
+ A = list(matrix(1, 2, 2)),
+ C = list(matrix(1, 2, 2)),
+ E = diag(1, 2)
+)
+#> Component model is not identified.
+#> Non-identified parameters are A, C
+#> $identified
+#> [1] FALSE
+#>
+#> $nidp
+#> [1] "A" "C"As you can see, the model is not identified. We need to add an +additional group so that we have sufficient information. Let us add the +rest of the classical twin model, in this case DZ twins.
+identifyComponentModel(
+ A = list(matrix(c(1, .5, .5, 1), 2, 2), matrix(1, 2, 2)),
+ C = list(matrix(1, 2, 2), matrix(1, 2, 2)),
+ E = diag(1, 4)
+)
+#> Component model is identified.
+#> $identified
+#> [1] TRUE
+#>
+#> $nidp
+#> character(0)As you can see the model is identified, now that we’ve added another +group. Let us confirm by fitting a model. First we prepare the data.
+require(dplyr)
+#> Loading required package: dplyr
+#>
+#> Attaching package: 'dplyr'
+#> The following objects are masked from 'package:BGmisc':
+#>
+#> between, first, last
+#> The following objects are masked from 'package:stats':
+#>
+#> filter, lag
+#> The following objects are masked from 'package:base':
+#>
+#> intersect, setdiff, setequal, union
+# require(purrr)
+
+data(twinData, package = "OpenMx")
+selVars <- c("ht1", "ht2")
+
+mzdzData <- subset(
+ twinData, zyg %in% c(1, 3),
+ c(selVars, "zyg")
+)
+
+mzdzData$RCoef <- c(1, NA, .5)[mzdzData$zyg]
+
+
+mzData <- mzdzData %>% filter(zyg == 1)Let us fit the data with MZ twins by themselves.
+run1 <- emxTwinModel(
+ model = "Cholesky",
+ relatedness = "RCoef",
+ data = mzData,
+ use = selVars,
+ run = TRUE, name = "TwCh"
+)
+#> Running TwCh with 4 parameters
+#> Warning: In model 'TwCh' Optimizer returned a non-zero status code 5. The
+#> Hessian at the solution does not appear to be convex. See
+#> ?mxCheckIdentification for possible diagnosis (Mx status RED).
+
+summary(run1)
+#> Summary of TwCh
+#>
+#> The Hessian at the solution does not appear to be convex. See ?mxCheckIdentification for possible diagnosis (Mx status RED).
+#>
+#> free parameters:
+#> name matrix row col Estimate Std.Error A lbound ubound
+#> 1 sqrtA11 sqrtA 1 1 0.05122646 NA 1e-06
+#> 2 sqrtC11 sqrtC 1 1 0.03518629 NA 0!
+#> 3 sqrtE11 sqrtE 1 1 0.02325722 0.0007017955 ! 0!
+#> 4 Mht1 Means ht1 1 1.62974908 0.0027023907
+#>
+#> Model Statistics:
+#> | Parameters | Degrees of Freedom | Fit (-2lnL units)
+#> Model: 4 1112 -3693.148
+#> Saturated: 5 1111 NA
+#> Independence: 4 1112 NA
+#> Number of observations/statistics: 569/1116
+#>
+#>
+#> ** Information matrix is not positive definite (not at a candidate optimum).
+#> Be suspicious of these results. At minimum, do not trust the standard errors.
+#>
+#> Information Criteria:
+#> | df Penalty | Parameters Penalty | Sample-Size Adjusted
+#> AIC: -5917.148 -3685.148 -3685.078
+#> BIC: -10747.543 -3667.773 -3680.471
+#> To get additional fit indices, see help(mxRefModels)
+#> timestamp: 2025-04-21 16:43:57
+#> Wall clock time: 0.04138112 secs
+#> optimizer: SLSQP
+#> OpenMx version number: 2.21.13
+#> Need help? See help(mxSummary)As you can see the model was unsuccessful because it was not +identified. But when we add another group, so that the model is +identified, the model now fits.
+run2 <- emxTwinModel(
+ model = "Cholesky",
+ relatedness = "RCoef",
+ data = mzdzData,
+ use = selVars,
+ run = TRUE, name = "TwCh"
+)
+#> Running TwCh with 4 parameters
+
+summary(run2)
+#> Summary of TwCh
+#>
+#> free parameters:
+#> name matrix row col Estimate Std.Error A lbound ubound
+#> 1 sqrtA11 sqrtA 1 1 0.06339271 0.0014377690 1e-06
+#> 2 sqrtC11 sqrtC 1 1 0.00000100 0.0250260004 ! 0!
+#> 3 sqrtE11 sqrtE 1 1 0.02330040 0.0007015267 0!
+#> 4 Mht1 Means ht1 1 1.63295540 0.0020511844
+#>
+#> Model Statistics:
+#> | Parameters | Degrees of Freedom | Fit (-2lnL units)
+#> Model: 4 1803 -5507.092
+#> Saturated: 5 1802 NA
+#> Independence: 4 1803 NA
+#> Number of observations/statistics: 920/1807
+#>
+#> Information Criteria:
+#> | df Penalty | Parameters Penalty | Sample-Size Adjusted
+#> AIC: -9113.092 -5499.092 -5499.048
+#> BIC: -17811.437 -5479.794 -5492.498
+#> To get additional fit indices, see help(mxRefModels)
+#> timestamp: 2025-04-21 16:43:57
+#> Wall clock time: 0.03366184 secs
+#> optimizer: SLSQP
+#> OpenMx version number: 2.21.13
+#> Need help? See help(mxSummary)ped_add_partial_complete <- ped2com(df,
isChild_method = "partialparent",
component = "additive",
- adjacency_method = "direct"
-)
-ped_add_classic_complete <- ped2com(df,
- isChild_method = "classic",
- component = "additive", adjacency_method = "direct"
-)The following plots display the full additive matrices. These matrices should be identical.
This can be confirmed visually and numerically.
@@ -392,25 +394,31 @@
-corrplot(as.matrix(ped_add_partial_complete),
- method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE, title = "Additive component - Partial parent method"
-)
To verify this, we subtract one matrix from the other and calculate RMSE. The difference should be numerically zero. Indeed, it is 0.
-corrplot(as.matrix(ped_add_classic_complete - ped_add_partial_complete),
+corrplot((as.matrix(ped_add_classic_complete) - as.matrix(ped_add_partial_complete)),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE
-)
-#> Warning in corrplot(as.matrix(ped_add_classic_complete -
-#> ped_add_partial_complete), : col.lim interval too wide, please set a suitable
-#> value
-
+ is.corr = FALSE, order = "hclust",
+ tl.pos = "l", tl.col = "black", tl.srt = 5, tl.cex = 0.2,
+ col = COL1("Reds", 100), mar = c(0, 0, 2, 0)
+)
ped_add_partial_mom <- ped_add_partial <- ped2com(df,
isChild_method = "partialparent",
component = "additive",
- adjacency_method = "direct"
-)
-
-ped_add_classic_mom <- ped_add_classic <- ped2com(df,
- isChild_method = "classic",
- component = "additive", adjacency_method = "direct"
-)The two methods now treat individual 4 differently in the parent adjacency matrix. The classic method applies a fixed contribution because one parent remains. The partial parent method inflates the @@ -436,26 +446,34 @@
corrplot(as.matrix(ped_add_classic),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE, title = "Classic (mother removed)"
-)
-corrplot(as.matrix(ped_add_partial),
- method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE, title = "Partial (mother removed)"
-)
corrplot(as.matrix(ped_add_partial),
+ method = "color", type = "lower", col.lim = c(0, 1),
+ is.corr = FALSE, title = "Partial (mother removed)",
+ order = "hclust",
+ tl.pos = "l", tl.col = "black", tl.srt = 5, tl.cex = 0.2,
+ col = COL1("Reds", 100), mar = c(0, 0, 2, 0)
+)
We quantify the overall matrix difference:
-sqrt(mean((ped_add_classic - ped_add_partial)^2))
+
Next, we compare each method to the matrix from the complete
pedigree. This evaluates how much each method deviates from the correct
additive structure.
-corrplot(as.matrix(ped_add_classic_complete - ped_add_classic),
+corrplot(as.matrix(ped_add_classic_complete) - as.matrix(ped_add_classic),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE
-)
-
+ is.corr = FALSE,
+ order = "hclust",
+ tl.pos = "l", tl.col = "black", tl.srt = 5, tl.cex = 0.2,
+ col = COL1("Reds", 100), mar = c(0, 0, 2, 0)
+)
+
@@ -463,9 +481,12 @@ Introducing Missingness: Remove a Parent
is 0.0299137.
corrplot(as.matrix(ped_add_classic_complete - ped_add_partial),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE
-)
-
+ is.corr = FALSE,
+ order = "hclust",
+ tl.pos = "l", tl.col = "black", tl.srt = 5, tl.cex = 0.2,
+ col = COL1("Reds", 100), mar = c(0, 0, 2, 0)
+)
corrplot(as.matrix(ped_add_classic_dad),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE, title = "Classic (father removed)"
-)
corrplot(as.matrix(ped_add_partial_dad),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE, title = "Partial (father removed)"
-)
Again, we compare to the true matrix from the complete pedigree:
corrplot(as.matrix(ped_add_classic_complete - ped_add_classic),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE
-)
corrplot(as.matrix(ped_add_classic_complete - ped_add_partial),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE
-)
The partial parent method again yields a matrix closer to the full-data version.
@@ -538,7 +575,7 @@For each one, we construct the additive relationship matrix under complete information and then simulate two missingness scenarios:
inbreeding_list <- list()
results <- data.frame(
- FamIDs = FamIDs,
- RMSE_partial_dad = rep(NA, length(FamIDs)),
- RMSE_partial_mom = rep(NA, length(FamIDs)),
- RMSE_classic_dad = rep(NA, length(FamIDs)),
- RMSE_classic_mom = rep(NA, length(FamIDs)),
- max_R_classic_dad = rep(NA, length(FamIDs)),
- max_R_partial_dad = rep(NA, length(FamIDs)),
- max_R_classic_mom = rep(NA, length(FamIDs)),
- max_R_partial_mom = rep(NA, length(FamIDs)),
- max_R_classic = rep(NA, length(FamIDs))
+ famIDs = famIDs,
+ RMSE_partial_dad = rep(NA, length(famIDs)),
+ RMSE_partial_mom = rep(NA, length(famIDs)),
+ RMSE_classic_dad = rep(NA, length(famIDs)),
+ RMSE_classic_mom = rep(NA, length(famIDs)),
+ max_R_classic_dad = rep(NA, length(famIDs)),
+ max_R_partial_dad = rep(NA, length(famIDs)),
+ max_R_classic_mom = rep(NA, length(famIDs)),
+ max_R_partial_mom = rep(NA, length(famIDs)),
+ max_R_classic = rep(NA, length(famIDs))
)The loop below performs this procedure for all families in the dataset and stores both the RMSEs and the maximum relatedness values.
-for (i in 1:length(FamIDs)) {
+for (i in 1:length(famIDs)) {
# make three versions to filter down
- df_fam_dad <- df_fam_mom <- df_fam <- df[df$FamID == FamIDs[i], ]
+ df_fam_dad <- df_fam_mom <- df_fam <- df[df$famID == famIDs[i], ]
results$RMSE_partial_mom[i] <- sqrt(mean((ped_add_classic_complete - ped_add_partial_mom)^2))
@@ -638,7 +675,85 @@ Inbreeding Dataset: Family-Level Evaluation
ped_add_partial_mom = ped_add_partial_mom,
ped_add_classic_mom = ped_add_classic_mom
)
-}
+}
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_dad)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_partial_mom)^2):
+#> argument is not numeric or logical: returning NA
+#> Warning in mean.default((ped_add_classic_complete - ped_add_classic_mom)^2):
+#> argument is not numeric or logical: returning NA1To understand what these matrices look like, we visualize them for @@ -651,57 +766,84 @@
1
corrplot(as.matrix(fam1$ped_add_classic_mom),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE, title = "Classic - Mom Missing"
-)
corrplot(as.matrix(fam1$ped_add_partial_mom),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE, title = "Partial - Mom Missing"
-)
corrplot(as.matrix(fam1$ped_add_classic_dad),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE, title = "Classic - Dad Missing"
-)
corrplot(as.matrix(fam1$ped_add_partial_dad),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE, title = "Partial - Dad Missing"
-)
To visualize the differences from the true matrix:
corrplot(as.matrix(fam1$ped_add_classic_complete - fam1$ped_add_classic_mom),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE, title = "Classic Mom Diff from Complete"
-)
corrplot(as.matrix(fam1$ped_add_classic_complete - fam1$ped_add_partial_mom),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE, title = "Partial Mom Diff from Complete"
-)
corrplot(as.matrix(fam1$ped_add_classic_complete - fam1$ped_add_classic_dad),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE, title = "Classic Dad Diff from Complete"
-)
corrplot(as.matrix(fam1$ped_add_classic_complete - fam1$ped_add_partial_dad),
method = "color", type = "lower", col.lim = c(0, 1),
- is.corr = FALSE, title = "Partial Dad Diff from Complete"
-)
These plots show how each method responds to missing data, and whether it maintains consistency with the complete pedigree. We observe that the partial parent method typically introduces smaller deviations. @@ -728,13 +870,14 @@
results %>%
select(RMSE_diff_mom, RMSE_diff_dad) %>%
summary()
-#> RMSE_diff_mom RMSE_diff_dad
-#> Min. :0.001222 Min. :0.001222
-#> 1st Qu.:0.001869 1st Qu.:0.002036
-#> Median :0.002538 Median :0.002520
-#> Mean :0.005763 Mean :0.005786
-#> 3rd Qu.:0.005625 3rd Qu.:0.005625
-#> Max. :0.024221 Max. :0.024221In all families, both RMSE_diff_mom and
RMSE_diff_dad are positive—indicating that the classic
method produces larger the errors relative to the partial method. This
@@ -742,9 +885,9 @@
To verify this directly:
mean(results$RMSE_diff_mom > 0, na.rm = TRUE)
-#> [1] 1
+#> [1] NaN
mean(results$RMSE_diff_dad > 0, na.rm = TRUE)
-#> [1] 1These proportions show how often the partial method produces a lower RMSE across the dataset. This confirms the earlier findings: when pedigree data are incomplete, the partial parent method more faithfully @@ -752,17 +895,18 @@
results %>%
as.data.frame() %>%
select(
- -FamIDs, -RMSE_diff_mom, -RMSE_diff_dad, -max_R_classic_dad,
+ -famIDs, -RMSE_diff_mom, -RMSE_diff_dad, -max_R_classic_dad,
-max_R_partial_dad, -max_R_classic_mom, -max_R_partial_mom, -max_R_classic
) %>%
summary()
-#> RMSE_partial_dad RMSE_partial_mom RMSE_classic_dad RMSE_classic_mom
-#> Min. :0.04773 Min. :0.04773 Min. :0.04895 Min. :0.04895
-#> 1st Qu.:0.05570 1st Qu.:0.05349 1st Qu.:0.05774 1st Qu.:0.05555
-#> Median :0.06206 Median :0.06899 Median :0.06457 Median :0.07158
-#> Mean :0.07545 Mean :0.07686 Mean :0.08124 Mean :0.08262
-#> 3rd Qu.:0.08237 3rd Qu.:0.08323 3rd Qu.:0.08866 3rd Qu.:0.08866
-#> Max. :0.15547 Max. :0.15547 Max. :0.17969 Max. :0.17969