add a test and rename internal function

Author: smasongarrison

R/func_kinsim.R

@@ -137,7 +137,7 @@ kinsim <- function(
       sigma_a[4, 1] <- cov_a * r_all[i]
       sigma_a[3, 2] <- cov_a * r_all[i]
       sigma_a[2, 3] <- cov_a * r_all[i]
-      A.r <- rmvn(n,
+      A.r <- .rmvn(n,
         sigma = sigma_a
       )
 
@@ -154,7 +154,7 @@ kinsim <- function(
       sigma_c[4, 1] <- cov_c * 1
       sigma_c[3, 2] <- cov_c * 1
       sigma_c[2, 3] <- cov_c * 1
-      C.r <- rmvn(n,
+      C.r <- .rmvn(n,
         sigma = sigma_c
       )
       C.r[, 1:2] <- C.r[, 1:2] * sC[1]
@@ -166,7 +166,7 @@ kinsim <- function(
       sigma_e[3, 1] <- cov_e
       sigma_e[2, 4] <- cov_e
       sigma_e[4, 2] <- cov_e
-      E.r <- rmvn(n,
+      E.r <- .rmvn(n,
         sigma = sigma_e
       )
       E.r[, 1:2] <- E.r[, 1:2] * sE[1]
@@ -225,7 +225,7 @@ kinsim <- function(
       sigma_a[3, 2] <- cov_a * r_val
       sigma_a[2, 3] <- cov_a * r_val
 
-      A_tmp <- rmvn(n_sub, sigma = sigma_a)
+      A_tmp <- .rmvn(n_sub, sigma = sigma_a)
 
       A.r[idx, 1:2] <- A_tmp[, 1:2] * sA[1]
       A.r[idx, 3:4] <- A_tmp[, 3:4] * sA[2]
@@ -241,7 +241,7 @@ kinsim <- function(
       sigma_c[3, 2] <- cov_c * 1
       sigma_c[2, 3] <- cov_c * 1
 
-      C_tmp <- rmvn(n_sub, sigma = sigma_c)
+      C_tmp <- .rmvn(n_sub, sigma = sigma_c)
       C.r[idx, 1:2] <- C_tmp[, 1:2] * sC[1]
       C.r[idx, 3:4] <- C_tmp[, 3:4] * sC[2]
 
@@ -252,7 +252,7 @@ kinsim <- function(
       sigma_e[2, 4] <- cov_e
       sigma_e[4, 2] <- cov_e
 
-      E_tmp <- rmvn(n_sub, sigma = sigma_e)
+      E_tmp <- .rmvn(n_sub, sigma = sigma_e)
       E.r[idx, 1:2] <- E_tmp[, 1:2] * sE[1]
       E.r[idx, 3:4] <- E_tmp[, 3:4] * sE[2]
     }

R/helpers_simulation.R

@@ -8,7 +8,7 @@
 #' @return Matrix of dimension \code{n × ncol(sigma)} containing random samples
 #'   from the multivariate normal distribution.
 #' @keywords internal
-rmvn <- function(n, sigma) {
+.rmvn <- function(n, sigma) {
   Sh <- with(
     svd(sigma),
     v %*% diag(sqrt(d)) %*% t(u)
@@ -89,16 +89,17 @@ kinsim_internal <- function(
 
 
     # Generate data for each relatedness group
-    for (i in 1:length(r)) {
+#    for (i in 1:length(r)) {
+    for (i in seq_along(r)) {
       n <- npergroup[i]
 
       # Generate correlated genetic components based on relatedness
-      A.r <- sA * rmvn(n, sigma = diag(2) + S2 * r[i])
+      A.r <- sA * .rmvn(n, sigma = diag(2) + S2 * r[i])
 
       # Generate shared environmental components (same for both members)
       #     C.r <- stats::rnorm(n,sd = sC)
       #   C.r <- cbind(C.r,C.r )
-      C.r <- sC * rmvn(n, sigma = diag(2) + S2 * c_rel)
+      C.r <- sC * .rmvn(n, sigma = diag(2) + S2 * c_rel)
 
       # Generate non-shared environmental components (different for each member)
       E.r <- cbind(
@@ -140,7 +141,7 @@ kinsim_internal <- function(
     # Generate genetic components for each unique relatedness value
     for (i in 1:length(unique_r)) {
       n <- length(r_vector[r_vector == unique_r[i]])
-      A.rz <- sA * rmvn(n, sigma = diag(2) + S2 * unique_r[i])
+      A.rz <- sA * .rmvn(n, sigma = diag(2) + S2 * unique_r[i])
       data_vector$A.r1[data_vector$r_vector == unique_r[i]] <- A.rz[, 1]
       data_vector$A.r2[data_vector$r_vector == unique_r[i]] <- A.rz[, 2]
     }
@@ -149,7 +150,7 @@ kinsim_internal <- function(
       data_vector$A.r1,
       data_vector$A.r2
     ), ncol = 2)
-    C.r <- sC * rmvn(n, sigma = diag(2) + S2 * c_rel)
+    C.r <- sC * .rmvn(n, sigma = diag(2) + S2 * c_rel)
 
     E.r <- cbind(
       stats::rnorm(n, sd = sE),

man/rmvn.Rd man/dot-rmvn.Rdman/rmvn.Rd renamed to man/dot-rmvn.Rd

@@ -161,11 +161,13 @@ test_that("genetic correlation between variables is present when cov_a is non-ze
   expect_gt(cor_val, 0.1) # minimal expected correlation
 })
 
-test_that("kinsim handles r_vector input correctly", {
+test_that("kinsim handles r_vector and c_vecttor input correctly", {
   r_vec <- rep(c(1, 0.5), each = 100)
-  df <- kinsim(r_vector = r_vec)
+  c_vector <- rep(1, 200)
+  df <- kinsim(r_vector = r_vec, c_vector = c_vector)
   expect_equal(nrow(df), 200)
   expect_equal(df$r, r_vec)
+  expect_equal(df$C1_1, df$C1_2)
 })
 
 test_that("output has correct ID range", {