R/helpers_simulation.R

In smasongarrison/discord: Functions for Discordant Kinship Modeling

#' rmvn
#' @keywords internal
#' @param n Sample Size
#' @param sigma Covariance matrix
#' @return Generates multivariate normal data from a covariance matrix (\code{sigma}) of length \code{n}


rmvn <- function(n, sigma) {
  Sh <- with(
    svd(sigma),
    v %*% diag(sqrt(d)) %*% t(u)
  )
  matrix(stats::rnorm(ncol(sigma) * n),
    ncol = ncol(sigma)
  ) %*% Sh
}

#' Simulate Biometrically informed Univariate Data
#'
#' @description Generate paired univariate data, given ACE parameters.
#' @param r Levels of relatedness; default is MZ and DZ twins c(1,.5)
#' @param npg Sample size per group; default is 100.
#' @param npergroup List of sample sizes by group; default repeats \code{npg} for all groups.
#' @param mu Mean for generated variable; default is 0.
#' @param ace Vector of variance components, ordered by c(a, c, e); default is c(1,1,1).
#' @param r_vector Alternative, give vector of relatedness coefficients for entire sample.
#' @param ... Optional pass on additional inputs.

#' @return Returns \code{data.frame} with the following:
#' \item{id}{id}
#' \item{A1}{genetic component for kin1}
#' \item{A2}{genetic component for kin2}
#' \item{C1}{shared-environmental component for kin1}
#' \item{C2}{shared-environmental component for kin2}
#' \item{E1}{non-shared-environmental component for kin1}
#' \item{E2}{non-shared-environmental component for kin2}
#' \item{y1}{generated variable for kin1 with mean of \code{mu}}
#' \item{y2}{generated variable for kin2 with mean of \code{mu}}
#' \item{r}{level of relatedness for the kin pair}


kinsim_internal <- function(
    r = c(1, .5),
    npg = 100,
    npergroup = rep(npg, length(r)),
    mu = 0,
    ace = c(1, 1, 1),
    r_vector = NULL,
    ...) {
  sA <- ace[1]^0.5
  sC <- ace[2]^0.5
  sE <- ace[3]^0.5

  S2 <- matrix(c(
    0, 1,
    1, 0
  ), 2)
  datalist <- list()

  if (is.null(r_vector)) {
    id <- 1:sum(npergroup)

    for (i in 1:length(r)) {
      n <- npergroup[i]

      A.r <- sA * rmvn(n,
        sigma = diag(2) + S2 * r[i]
      )
      C.r <- stats::rnorm(n,
        sd = sC
      )
      C.r <- cbind(
        C.r,
        C.r
      )
      E.r <- cbind(
        stats::rnorm(n,
          sd = sE
        ),
        stats::rnorm(n,
          sd = sE
        )
      )

      y.r <- mu + A.r + C.r + E.r


      r_ <- rep(r[i], n)

      data.r <- data.frame(A.r, C.r, E.r, y.r, r_)
      names(data.r) <- c("A1", "A2", "C1", "C2", "E1", "E2", "y1", "y2", "r")
      datalist[[i]] <- data.r
      names(datalist)[i] <- paste0("datar", r[i])
    }
    merged.data.frame <- Reduce(function(...) merge(..., all = T), datalist)
    merged.data.frame$id <- id
  } else {
    id <- 1:length(r_vector)
    data_vector <- data.frame(id, r_vector)
    data_vector$A.r1 <- as.numeric(NA)
    data_vector$A.r2 <- as.numeric(NA)
    unique_r <- matrix(unique(r_vector))
    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]
      )
      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]
    }
    n <- length(r_vector)
    A.r <- matrix(c(
      data_vector$A.r1,
      data_vector$A.r2
    ), ncol = 2)
    C.r <- stats::rnorm(n, sd = sC)
    C.r <- cbind(C.r, C.r)
    E.r <- cbind(
      stats::rnorm(n,
        sd = sE
      ),
      stats::rnorm(n,
        sd = sE
      )
    )

    y.r <- mu + A.r + C.r + E.r

    data.r <- data.frame(id, A.r, C.r, E.r, y.r, r_vector)
    names(data.r) <- c("id", "A1", "A2", "C1", "C2", "E1", "E2", "y1", "y2", "r")
    datalist[[i]] <- data.r
    names(datalist)[i] <- paste0("datar", r[i])

    merged.data.frame <- data.r
  }

  return(merged.data.frame)
}