updated MphEM function and its documentation and tests.

Author: fboehm

DESCRIPTION

@@ -12,7 +12,7 @@ Suggests: covr,
     testthat,
     knitr,
     rmarkdown
-RoxygenNote: 6.0.1
+RoxygenNote: 6.1.0
 VignetteBuilder: knitr
 Imports: readr,
     Matrix

R/MphEM.R

@@ -7,10 +7,12 @@
 #' @param Y matrix of phenotype values
 #' @param V_g genetic covariance matrix
 #' @param V_e error covariance matrix
+#' @param verbose_output logical indicating whether to output entire collection of intermediate values for all iterations. Default is FALSE.
 #' @return a list of lists. Length of list corresponds to number of EM iterations
 #' @export
 MphEM <- function(max_iter = 10000, max_prec = 1 / 1000000,
-                  eval, X, Y, V_g, V_e){
+                  eval, X, Y, V_g, V_e,
+                  verbose_output = FALSE){
   n_size <- length(eval)
   # be careful with defining c_size here
 
@@ -74,13 +76,23 @@ MphEM <- function(max_iter = 10000, max_prec = 1 / 1000000,
     # update V_g and V_e
     uv_out[[1]] -> V_g
     uv_out[[2]] -> V_e
-    out[[t]] <- list(logl_new = logl_new, Vg = V_g, Ve = V_e, Sigma_uu = Sigma_uu,
+    if (verbose_output) {
+      out[[t]] <- list(logl_new = logl_new, Vg = V_g, Ve = V_e, Sigma_uu = Sigma_uu,
                      Sigma_ee = Sigma_ee, B = B,
                      U_hat = U_hat, E_hat = E_hat,
                      OmegaU = OmegaU, OmegaE = OmegaE, logdet_Ve = logdet_Ve,
                      UltVeh = UltVeh, UltVehi = UltVehi,
                      Dl = D_l, xHiy = xHiy, logl_const = logl_const, UltVehiU = UltVehiU
                      )
+    } else {
+      out[[1]] <- list(logl_new = logl_new, Vg = V_g, Ve = V_e, Sigma_uu = Sigma_uu,
+                       Sigma_ee = Sigma_ee, B = B,
+                       U_hat = U_hat, E_hat = E_hat,
+                       OmegaU = OmegaU, OmegaE = OmegaE, logdet_Ve = logdet_Ve,
+                       UltVeh = UltVeh, UltVehi = UltVehi,
+                       Dl = D_l, xHiy = xHiy, logl_const = logl_const, UltVehiU = UltVehiU
+      )
+    }
   }
   if (length(out) == max_iter){warning("EM algorithm didn't converge.")}
   return(out)

man/MphEM.Rd

@@ -8,9 +8,16 @@ as.matrix(pheno[, c(1, 6)]) -> phe16
 MphEM(max_iter = 10, eval = e_out$values, X = matrix(c(-10, rep(0, 99)), nrow = 1),
       Y = t(phe16) %*% e_out$vectors,
       V_g = matrix(c(1.91352, 0, 0, 0.530827), nrow = 2),
-      V_e = matrix(c( 0.320028, 0, 0, 0.561589), nrow = 2)
+      V_e = matrix(c( 0.320028, 0, 0, 0.561589), nrow = 2), verbose_output = TRUE
 )-> foo
 
+MphEM(max_iter = 10, eval = e_out$values, X = matrix(c(-10, rep(0, 99)), nrow = 1),
+      Y = t(phe16) %*% e_out$vectors,
+      V_g = matrix(c(1.91352, 0, 0, 0.530827), nrow = 2),
+      V_e = matrix(c( 0.320028, 0, 0, 0.561589), nrow = 2), verbose_output = FALSE
+)-> bar
+
+
 context("Testing MphEM() with two phenotypes and intercept only (no genetic data)")
 
 
@@ -21,3 +28,11 @@ test_that("First calculations of for Vg and for Ve are accurate", {
   expect_equal(foo[[5]][[2]], matrix(c(1.89405, 0.112242, 0.112242, 0.531075), nrow = 2), tolerance = 0.00001)
   expect_equal(foo[[5]][[3]], matrix(c(0.324097, 0.153797, 0.153797, 0.561429), nrow = 2), tolerance = 0.00001)
 })
+
+
+test_that("MphEM gives same output, regardless of value of input `verbose_output`", {
+  expect_equal(foo[[10]][[2]], bar[[1]][[2]])
+  expect_equal(foo[[10]][[3]], bar[[1]][[3]])
+  expect_length(foo, 10)
+  expect_length(bar, 1)
+})