Nothing
#' GE_translate_inputs_old.R
#'
#' Mostly for internal use, function called by GE_bias_normal_old() and GE_scoreeq_sim_old()
#' to translate the rho_list inputs and return a total covariance matrix for simulation/
#' checking validity of covariance structure. If invalid covariance structure, will stop
#' and return an error message.
#'
#' @param rho_list A list of the 6 pairwise covariances between the
#' covariates. These should be in the order (1) cov_GE (2) cov_GZ (3) cov_EZ
#' (4) cov_GW (5) cov_EW (6) cov_ZW. If Z or M are vectors then terms like cov_GZ should be vectors
#' (in the appropriate order).
#' If Z or M are vectors, then cov_ZW should be a vector in the order (cov(Z_1,W_1),...,cov(Z_1,W_q),
#' cov(Z_2,W_1),........,cov(Z_p,W_q) where Z is a vector of length p and W is a vector of length q.
#' @param cov_Z Only used if Z is a vector, gives the covariance matrix of Z (remember by assumption
#' Z has mean 0 and variance 1). The (i,j) element of the matrix should be the (i-1)(i-2)/2+j element
#' of the vector.
#' @param cov_W Only used if W is a vector, gives the covariance matrix of W (remember by assumption
#' W has mean 0 and variance 1). The (i,j) element of the matrix should be the (i-1)(i-2)/2+j element
#' of the vector.
#' @param prob_G Probability that each allele is equal to 1. Since each SNP has
#' two alleles, the expectation of G is 2*prob_G.
#'
#' @return A list with the elements:
#' \item{sig_mat_total}{The sigma parameter for rmvnorm call to generate our data.}
#' \item{sig_mat_ZZ}{The covariance matrix of Z, i.e. E[ZZ^T]}
#' \item{sig_mat_WW}{The covariance matrix of W, i.e. E[WW^T]}
#'
#' @keywords internal
#' @export
#' @examples
#' GE_translate_inputs_old( beta_list=as.list(runif(n=6, min=0, max=1)),
#' rho_list=as.list(rep(0.3,6)), prob_G=0.3)
GE_translate_inputs_old <- function(beta_list, rho_list, prob_G, cov_Z=NULL, cov_W=NULL)
{
# First, make sure we got good inputs
if (length(beta_list) != 6 | length(rho_list) != 6 | !inherits(beta_list, "list") | !inherits(rho_list, "list"))
{
stop('Input vectors not the right size!')
}
# How long are these vectors? Remember that W is the same length as M by assumption.
# Check for Z and M/W nonexistence LATER
num_Z <- length(beta_list[[5]])
num_W <- length(beta_list[[6]])
# Make sure we have compatible lengths for rho_list
if (length(rho_list[[2]]) != num_Z | length(rho_list[[3]]) != num_Z | length(rho_list[[4]]) != num_W
| length(rho_list[[5]]) != num_W | length(rho_list[[6]]) != num_Z*num_W) {
stop('Incompatible number of elements in beta/rho_list')
}
# Fill in our covariances.
rho_GE <- rho_list[[1]]; rho_GZ <- rho_list[[2]]; rho_EZ <- rho_list[[3]]
rho_GW <- rho_list[[4]]; rho_EW <- rho_list[[5]]; rho_ZW <- rho_list[[6]]
# Now check for Z and M/W nonexistence, after confirming rho_list is correct length
# If no Z or M/W, then we need to set num_Z and num_W=0 in this function so that
# the indicies will be correct later for building the covariance matrix.
if (length(beta_list[[5]]) == 1 & beta_list[[5]][1] == 0) { num_Z <- 0 }
if (length(beta_list[[6]]) == 1 & beta_list[[6]][1] == 0) { num_W <- 0 }
################################################################
# Build our covariance matrix in steps.
# The 3x3 in the top left is always the same to build, vectors or not.
w <- qnorm(1-prob_G) # Threshold for generating G
r_GE <- rho_GE / (2*dnorm(w))
sig_mat_GE <- matrix(data=c(1, 0, r_GE, 0, 1, r_GE, r_GE, r_GE, 1), nrow=3)
# Build the p*3 matrix that describes Z with G1,G2,E
if (num_Z != 0) {
sig_mat_Z_column <- matrix(data=NA, nrow=num_Z, ncol=3)
r_GZ <- rho_GZ / (2*dnorm(w))
for (i in 1:num_Z) {
sig_mat_Z_column[i,] <- c(r_GZ[i], r_GZ[i], rho_EZ[i])
}
}
# Build the q*3 matrix that describes W with G1,G2,E
if (num_W != 0)
{
sig_mat_W_column <- matrix(data=NA, nrow=num_W, ncol=3)
r_GW <- rho_GW / (2*dnorm(w))
for (i in 1:num_W) {
sig_mat_W_column[i,] <- c(r_GW[i], r_GW[i], rho_EW[i])
}
}
# Build the p*q matrix that describes Z with W.
if (num_Z != 0 & num_W != 0) {
sig_mat_Z_W <- matrix(data=NA, nrow=num_Z, ncol=num_W)
for (i in 1:num_Z) {
start_ind <- (i-1)*num_W+1
end_ind <- i*num_W
sig_mat_Z_W[i,] <- rho_ZW[start_ind:end_ind]
}
} # No else here, just don't build it later.
# If Z or W vectorized, build the ZZ and WW covariance matrices too
if (num_Z > 1) {
sig_mat_ZZ <- matrix(data=0, nrow=num_Z, ncol=num_Z)
sig_mat_ZZ[upper.tri(sig_mat_ZZ)] <- cov_Z
sig_mat_ZZ <- sig_mat_ZZ + t(sig_mat_ZZ)
diag(sig_mat_ZZ) <- 1
} else if (num_Z == 0) { # If no Z, then return NULL for sig_mat_ZZ.
sig_mat_ZZ <- NULL
} else if (num_Z == 1) {
sig_mat_ZZ <- matrix(data=1, nrow=1, ncol=1)
}
if (num_W > 1) {
sig_mat_WW <- matrix(data=0, nrow=num_W, ncol=num_W)
sig_mat_WW[upper.tri(sig_mat_WW)] <- cov_W
sig_mat_WW <- sig_mat_WW + t(sig_mat_WW)
diag(sig_mat_WW) <- 1
} else if (num_W == 0) { # If no W, then return NULL for sig_mat_WW.
sig_mat_WW <- NULL
} else if (num_W == 1) {
sig_mat_WW <- matrix(data=1, nrow=1, ncol=1)
}
# Now put it all together.
# Top left 3x3 always there.
sig_mat_total <- matrix(data=NA, nrow=(3+num_Z+num_W), ncol=(3+num_Z+num_W))
sig_mat_total[1:3, 1:3] <- sig_mat_GE
# If we have a Z
if (num_Z > 0) {
sig_mat_total[4:(3+num_Z), 1:3] <- sig_mat_Z_column
sig_mat_total[1:3, 4:(3+num_Z)] <- t(sig_mat_Z_column)
sig_mat_total[4:(3+num_Z), 4:(3+num_Z)] <- sig_mat_ZZ
}
# If we have a W
if (num_W > 0) {
sig_mat_total[(4+num_Z):(3+num_Z+num_W), 1:3] <- sig_mat_W_column
sig_mat_total[1:3, (4+num_Z):(3+num_Z+num_W)] <- t(sig_mat_W_column)
sig_mat_total[(4+num_Z):(3+num_Z+num_W), (4+num_Z):(3+num_Z+num_W)] <- sig_mat_WW
}
# If we have both Z and W
if (num_Z > 0 & num_W > 0)
{
sig_mat_total[4:(3+num_Z), (4+num_Z):(3+num_Z+num_W)] <- sig_mat_Z_W
sig_mat_total[(4+num_Z):(3+num_Z+num_W), 4:(3+num_Z)] <- t(sig_mat_Z_W)
}
# Final sanity check that the building went correctly.
if (!isSymmetric(sig_mat_total)) {stop("Problem building covariance matrix!")}
# Now make sure we can actually generate data with this structure
test_data <- tryCatch(mvtnorm::rmvnorm(n=1, sigma=sig_mat_total),
warning=function(w) w, error=function(e) e)
if (class(test_data)[1] != 'matrix') {stop('You specified an impossible covariance matrix!')}
return(list(sig_mat_total=sig_mat_total, sig_mat_ZZ=sig_mat_ZZ, sig_mat_WW=sig_mat_WW))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.