add CIs

Author: smasongarrison

R/buildmxPedigrees.R

@@ -287,10 +287,12 @@ buildFamilyGroups <- function(
 #' @param model_name Name of the overall pedigree model.
 #' @param vars A named list or vector of initial variance component values.
 #' @param group_models A list of OpenMx models for each family group.
+#' @param ci Logical. If TRUE, include confidence interval computations for the variance components. Default is FALSE
 #' @return An OpenMx pedigree model combining variance components and family groups.
 #' @export
 
-buildPedigreeMx <- function(model_name, vars, group_models) {
+buildPedigreeMx <- function(model_name, vars, group_models,
+                            ci = FALSE) {
   if (!requireNamespace("OpenMx", quietly = TRUE)) {
     stop("OpenMx package is required for buildPedigreeMx function. Please install it.")
   }
@@ -335,7 +337,12 @@ buildPedigreeMx <- function(model_name, vars, group_models) {
       Ver = isTRUE(flags$Ver)
     ),
     group_models,
-    OpenMx::mxFitFunctionMultigroup(group_names)
+    OpenMx::mxFitFunctionMultigroup(group_names),
+    ci = if (ci & any(flags$Vad, flags$Vdd, flags$Vcn, flags$Vce, flags$Vmt, flags$Vam, flags$Ver)) {
+      OpenMx::mxCI(c("vad", "vdd", "vcn", "vce", "vmt", "vam", "ver")[c(flags$Vad, flags$Vdd, flags$Vcn, flags$Vce, flags$Vmt, flags$Vam, flags$Ver)])
+    } else {
+      NULL
+    }
   )
 }
 

vignettes/articles/v61_pedigree_model_fitting.Rmd

@@ -156,21 +156,24 @@ group_models_mace <- lapply(seq_len(n_families), function(i) {
 multi_model_ace <- buildPedigreeMx(
   model_name   = "MultiPedigreeModel",
   vars         = start_vars,
-  group_models = group_models_ace
+  group_models = group_models_ace,
+  ci = TRUE
 )
 
 multi_model_mace <- buildPedigreeMx(
   model_name   = "MultiPedigreeModel",
   vars         = start_vars,
-  group_models = group_models_mace
+  group_models = group_models_mace,
+  ci = TRUE
 )
 ```
 ```{r, include = FALSE}
 fitted_multi_ace <- NULL # ensure variable exists for later use
 fitted_multi_mace <- NULL # ensure variable exists for later use
 ```
 ```{r, eval = run_models}
-fitted_multi_ace <- mxRun(multi_model_ace, intervals= TRUE)
+fitted_multi_ace <- mxRun(multi_model_ace,
+                          intervals= TRUE)
 saveRDS(fitted_multi_ace, "inst/extdata/fitted_multi_ace.rds")
 
 fitted_multi_mace <- mxRun(multi_model_mace, intervals= TRUE)
@@ -220,3 +223,12 @@ cat("Common nuclear  (Vcn):", fitted_multi_mace$ModelOne$Vcn$values/total_var_ma
 cat("Mitochondrial (Vmt):", fitted_multi_mace$ModelOne$Vmt$values/total_var_mace, "\n")
 cat("Unique environ. (Ver):", fitted_multi_mace$ModelOne$Ver$values/total_var_mace, "\n")
 ```
+Now we can compare the ACE and MACE models to see if including the mitochondrial component changes our estimates of heritability and environmental contributions. This can provide insights into the role of mitochondrial inheritance in the trait of interest (LRS) in these squirrels.
+
+```{r, eval = has_openmx && !is.null(fitted_multi_ace) && !is.null(fitted_multi_mace)}
+
+mxCompare(fitted_multi_mace, fitted_multi_ace)
+
+fitted_multi_ace_ci <- omxRunCI(fitted_multi_ace, verbose = 0, optimizer = "SLSQP")
+
+```