GEint: Misspecified Models for Gene-Environment Interaction

Documented in GE_bias_normal_squaredmis

#' GE_bias_normal_squaredmis.R
#'
#' A function to calculate the bias in testing for GxE interaction, making many more
#' assumptions than GE_bias().  The additional assumptions are added to simplify the process
#' of calculating/estimating many higher order moments which the user may not be familiar with. \cr
#' The following assumptions are made: \cr
#' (1) All fitted covariates besides G (that is, E, all Z, and all W) have a marginal standard 
#' normal distribution with mean 0 and variance 1.  This corresponds to the case of the researcher
#' standardizing all of their fitted covariates. \cr
#' (2) All G are generated by means of thresholding two independent normal RVs and are centered to have mean 0. \cr
#' (3) The joint distributions of E, Z, W, and the thresholded variables underlying G can be described
#' by a multivariate normal distribution. \cr
#' (4) The misspecification is of the form f(E)=h(E)=E^2, and M_j=W_j^2 for all j. In particular,
#' W always has the same length as M here. \cr
#' 
#' @param beta_list A list of the effect sizes in the true model.
#' Use the order beta_0, beta_G, beta_E, beta_I, beta_Z, beta_M.
#' If G or Z or M is a vector, then beta_G/beta_Z/beta_M should be vectors.
#' If Z and/or M/W do not exist in your model, then set beta_Z and/or beta_M = 0.
#' @param rho_list A list of expectations (which happen to be covariances if all covariates
#' are centered at 0) in the order specified by GE_enumerate_inputs().
#' If Z and/or M/W do not exist in your model, then treat them as constants 0. For example,
#' if Z doesn't exist and W includes 2 covariates, then set cov(EZ) = 0 and cov(ZW) = (0,0).
#' If describing expectations relating two vectors, i.e. Z includes two covariates and W
#' includes three covariates, sort by the first term and then the second. Thus in the 
#' example, the first three terms of cov(ZW) are cov(Z_1,W_1),cov(Z_1,W_2), cov(Z_1,W_3), 
#' and the last three terms are cov(Z_3,W_1), cov(Z_3,W_2), cov(Z_3,W_3).
#' @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. Should be a d*1 vector.
#' @param cov_Z Should be a matrix equal to cov(Z) or NULL if no Z.
#' @param cov_W Should be a matrix equal to cov(W) or NULL if no W.  
#' @param corr_G Should be a matrix giving the *pairwise correlations* between each SNP
#' in the set, or NULL. Must be specified if G is a vector.  For example, the [2,3] element
#' of the matrix would be the pairwise correlation between SNP2 and SNP3.

#' 
#' @return A list with the elements:
#' \item{alpha_list}{The asymptotic values of the fitted coefficients alpha.}
#' \item{beta_list}{The same beta_list that was given as input.}
#' \item{cov_list}{The list of all covariances (both input and calculated) for use with GE_nleqslv() 
#'	and GE_bias().}
#' \item{mu_list}{List of calculated means for f(E), h(E), Z, M, and W for use with GE_nleqslv() 
#' and GE_bias().}
#' \item{HOM_list}{List of calculated Higher Order Moments for use with GE_nleqslv() and GE_bias().}
#'
#' @export
#' @examples 
#' GE_bias_normal_squaredmis( beta_list=as.list(runif(n=6, min=0, max=1)),
#' rho_list=as.list(rep(0.3,6)), cov_Z=1, cov_W=1, prob_G=0.3)

GE_bias_normal_squaredmis <- function(beta_list, rho_list, prob_G, cov_Z=NULL, cov_W=NULL, corr_G=NULL)
{
  # Need survival function.
  surv <- function(x) {1-pnorm(x)}
  
  # Record some initial quantities
  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]]
  	
  beta_0 <- beta_list[[1]]; beta_G <- beta_list[[2]]; beta_E <- beta_list[[3]]
  beta_I <- beta_list[[4]]; beta_Z <- beta_list[[5]]; beta_M <- beta_list[[6]]
  
  # How long are they
  num_G <- length(beta_list[[2]])
  num_W <- length(beta_list[[6]]); if (num_W == 1 & beta_list[[6]][1] == 0) {num_W <- 0}
  num_Z <- length(beta_list[[5]]); if (num_Z == 1 & beta_list[[5]][1] == 0) {num_Z <- 0}
  
  # Some error checking, make sure the covariance matrix is ok
  translated_inputs <- GE_translate_inputs(beta_list=beta_list, rho_list=rho_list, 
                                               prob_G=prob_G, cov_Z=cov_Z, cov_W=cov_W, corr_G=corr_G)
  sig_mat <- translated_inputs$sig_mat_total
  
  # First, the covariance matrix of G
  if (num_G > 1) 
  {
    mu_GG <- matrix(data=NA, nrow=num_G, ncol=num_G)
    for (i in 1:num_G) {
      for (j in 1:num_G) {
       sd1 <- sqrt(2*prob_G[i]*(1-prob_G[i]))
       sd2 <- sqrt(2*prob_G[j]*(1-prob_G[j]))
       mu_GG[i, j] <- corr_G[i,j] * sd1 * sd2
     }
    }
  } else {
    mu_GG <- 2*prob_G*(1-prob_G)
  }
  
  # Start with various GE moments
  w_vec <- qnorm(1-prob_G)
  r_GE <- rho_GE / (2*dnorm(w_vec))
  mu_f <- 1     # By assumption
  mu_h <- 1
  mu_GE <- rho_GE
  mu_Gf <- 2*r_GE^2*w_vec*dnorm(w_vec) + 2*surv(w_vec) - 2*prob_G
  mu_Gh <- mu_Gf
  mu_EE <- 1
  mu_Ef <- 0
  mu_EM <- 	rep(0, max(1, num_W))				# Because third moment of E is 0
  mu_fW <- 	rep(0, max(1, num_W))       # Same as mu_EM
  mu_EW <- rho_EW			
  mu_EZ <- rho_EZ	
  mu_fZ <- 	rep(0, max(1, num_Z))	
  mu_GEE <- mu_Gf
  mu_GEf <- 2*(r_GE^3*w_vec^2*dnorm(w_vec) - r_GE^3*dnorm(w_vec) + 3*r_GE*dnorm(w_vec))
  mu_GEh <- mu_GEf
  
  # Now do covariances with Z
  if (num_Z == 0) {
    mu_Z <- as.matrix(0)
    mu_GZ <- matrix(data=0, nrow=num_G, ncol=1)
    mu_ZG <- matrix(data=0, nrow=1, ncol=num_G)
    mu_ZW <- matrix(data=0, nrow=1, ncol=max(1, num_W))
    mu_WZ <- matrix(data=0, nrow=max(1, num_W), ncol=1)
    mu_ZM <- matrix(data=0, nrow=1, ncol=max(1, num_W))
    mu_MZ <- matrix(data=0, nrow=max(1, num_W), ncol=1)
    r_GZ <- matrix(data=0, nrow=num_G, ncol=1)
  } else {
    mu_Z <- rep(0, num_Z)     # By assumption
    r_GZ <- matrix(data=NA, nrow=num_G, ncol=num_Z)
    mu_GZ <- matrix(data=NA, nrow=num_G, ncol=num_Z)
    mu_ZG <- matrix(data=NA, nrow=num_Z, ncol=num_G)
    for (i in 1:num_G) {
      r_GZ[i, ] <- rho_GZ[((i-1)*num_Z+1):(i*num_Z)] / (2*dnorm(w_vec[i]))
      mu_GZ[i, ] <- rho_GZ[((i-1)*num_Z+1):(i*num_Z)]
      mu_ZG[, i] <- rho_GZ[((i-1)*num_Z+1):(i*num_Z)]
    }
    
    # Do ZW covariances
    # MZ is 0 for the same reason fW and fZ are 0.
    mu_MZ <- matrix(data=0, nrow=max(1, num_W), ncol=num_Z)
    mu_ZM <- matrix(data=0, nrow=num_Z, ncol=max(1, num_W))
    if (num_W == 0) {
      mu_ZW <- matrix(data=0, nrow=num_Z, ncol=1)
      mu_WZ <- matrix(data=0, nrow=1, ncol=num_Z)
    } else {
      mu_ZW <- matrix(data=NA, nrow=num_Z, ncol=num_W)
      mu_WZ <- matrix(data=NA, nrow=num_W, ncol=num_Z)
      for (i in 1:num_Z) {
        mu_ZW[i, ] <- rho_ZW[((i-1)*num_W+1):(i*num_W)]
        mu_WZ[, i] <- rho_ZW[((i-1)*num_W+1):(i*num_W)]
      }
    }
  }
  
  # Now do covariances with W
  if (num_W == 0) {
    mu_W <- as.matrix(0)
    mu_M <- as.matrix(0)
    mu_WM <- as.matrix(0)
    mu_GW <- matrix(data=0, nrow=num_G, ncol=1)
    mu_WG <- matrix(data=0, nrow=1, ncol=num_G)
    r_GW <- matrix(data=0, nrow=num_G, ncol=1)
    mu_GM <- matrix(data=0, nrow=num_G, ncol=1)
  } else {
    mu_M <- rep(1, num_W)		
    mu_W <- rep(0, num_W)	
    mu_WM <- matrix(data=0, nrow=num_W, ncol=num_W)
    r_GW <- matrix(data=NA, nrow=num_G, ncol=num_W)
    mu_GW <- matrix(data=NA, nrow=num_G, ncol=num_W)
    mu_WG <- matrix(data=NA, nrow=num_W, ncol=num_G)
    mu_GM <- matrix(data=NA, nrow=num_G, ncol=num_W)
    for (i in 1:num_G) {
      r_GW[i, ] <- rho_GW[((i-1)*num_W+1):(i*num_W)] / (2*dnorm(w_vec[i]))
      mu_GW[i, ] <- rho_GW[((i-1)*num_W+1):(i*num_W)]
      mu_WG[, i] <- rho_GW[((i-1)*num_W+1):(i*num_W)]
      mu_GM[i, ]  <- 	2*r_GW[i, ]^2*w_vec[i]*dnorm(w_vec[i]) + 2*surv(w_vec[i]) - 2*prob_G[i]	
    }
  }
  
  ########################
  # Higher order moments with G+E
  mu_GEG <- matrix(data=NA, nrow=num_G, ncol=num_G)
  mu_GhG <- matrix(data=NA, nrow=num_G, ncol=num_G)
  mu_GEEG <- matrix(data=NA, nrow=num_G, ncol=num_G)
  mu_GEfG <- matrix(data=NA, nrow=num_G, ncol=num_G)
  mu_GEhG <- matrix(data=NA, nrow=num_G, ncol=num_G)
  
  # Fill the diagonals first
  # Here G = G1 + G2 where G1 and G2 are independent binaries
  for (i in 1:num_G) {
    temp_w <- w_vec[i]
    temp_r_GE <- r_GE[i]
    temp_prob_G <- prob_G[i]
    
    # We need as intermediate quantities E[G_1E], E[G_1E^2], E[G_1E^3], E[G_1G_2E], E[G_1G_2E^2], E[G1EZ], 
    # E[G1EW], E[G1EW^2], E[G1WE^2], E[G1ZE^2], E[G1G2E^2]
    mu_G1_E <- temp_r_GE*dnorm(temp_w)
    mu_G1_EE <- temp_r_GE^2*temp_w*dnorm(temp_w) + surv(temp_w)
    mu_G1_EEE <- temp_r_GE^3*temp_w^2*dnorm(temp_w) - temp_r_GE^3*dnorm(temp_w) + 3*temp_r_GE*dnorm(temp_w)
    
    # E[G1G2E] requires numerical integration
    temp_sig <- matrix(data=c(1-temp_r_GE^2, -temp_r_GE^2, -temp_r_GE^2, 1-temp_r_GE^2), nrow=2)
    f_G1_G2_E <- function(x,w,r_GE) {
      x*dnorm(x)*mvtnorm::pmvnorm(lower=c(w,w), upper=c(Inf,Inf), mean=c(r_GE*x, r_GE*x), sigma=temp_sig)
    }
    mu_G1_G2_E <- pracma::quadinf(f=f_G1_G2_E, xa=-Inf, xb=Inf, w=temp_w, r_GE=temp_r_GE)$Q[1]
    
    # E[G1G2E^2] requires numerical integration
    temp_sig <- matrix(data=c(1-temp_r_GE^2, -temp_r_GE^2, -temp_r_GE^2, 1-temp_r_GE^2), nrow=2)
    f_G1_G2_EE <- function(x,w,r_GE) {
      x^2*dnorm(x)*mvtnorm::pmvnorm(lower=c(w,w), upper=c(Inf,Inf), mean=c(r_GE*x, r_GE*x), sigma=temp_sig)
    }
    mu_G1_G2_EE <- pracma::quadinf(f=f_G1_G2_EE, xa=-Inf, xb=Inf, w=temp_w, r_GE=temp_r_GE)$Q[1]
    
    # E[G1G2E^3] requires numerical integration
    temp_sig <- matrix(data=c(1-temp_r_GE^2, -temp_r_GE^2, -temp_r_GE^2, 1-temp_r_GE^2), nrow=2)
    f_G1_G2_EEE <- function(x,w,r_GE) {
      x^3*dnorm(x)*mvtnorm::pmvnorm(lower=c(w,w), upper=c(Inf,Inf), mean=c(r_GE*x, r_GE*x), sigma=temp_sig)
    }
    mu_G1_G2_EEE <- pracma::quadinf(f=f_G1_G2_EEE, xa=-Inf, xb=Inf, w=temp_w, r_GE=temp_r_GE)$Q[1]
    
    # See gen_cor_bin_normal to see how to do these
    mu_GEG[i, i] <- 2*mu_G1_E + 2*mu_G1_G2_E - 8*temp_prob_G*mu_G1_E
    mu_GhG[i, i] <- 2*mu_G1_EE + 2*mu_G1_G2_EE + 4*temp_prob_G^2*1 - 8*temp_prob_G*mu_G1_EE
    mu_GEEG[i, i] <- 2*mu_G1_EE + 2*mu_G1_G2_EE + 4*temp_prob_G^2*1 - 8*temp_prob_G*mu_G1_EE
    mu_GEfG[i, i] <- 2*mu_G1_EEE + 2*mu_G1_G2_EEE + 4*temp_prob_G^2*0 - 8*temp_prob_G*mu_G1_EEE
    mu_GEhG[i, i] <- mu_GEfG[i, i]
  }

  
  # Now fill the off diagonals
  # Here G = G1 + G2 where G1 and G2 are independent binaries
  f_G11_G21_X <- function(x, w1, w2, r_G1X, r_G2X, sig_mat_int)
  {
    x*dnorm(x)*mvtnorm::pmvnorm(lower=c(w1,w2), upper=c(Inf,Inf), mean=c(r_G1X*x,r_G2X*x), sigma=sig_mat_int)
  }
  f_G11_G21_XX <- function(x, w1, w2, r_G1X, r_G2X, sig_mat_int)
  {
    x^2*dnorm(x)*mvtnorm::pmvnorm(lower=c(w1,w2), upper=c(Inf,Inf), mean=c(r_G1X*x,r_G2X*x), sigma=sig_mat_int)
  }
  f_G11_G21_XXX <- function(x, w1, w2, r_G1X, r_G2X, sig_mat_int)
  {
    x^3*dnorm(x)*mvtnorm::pmvnorm(lower=c(w1,w2), upper=c(Inf,Inf), mean=c(r_G1X*x,r_G2X*x), sigma=sig_mat_int)
  }
  
  # Now do the off-diagonals
  if (num_G > 1) {
    # Use bindata to get the r_G1G2 values 
    G_cor_struct <- corr_G
    cprob <- bindata::bincorr2commonprob(margprob = c(prob_G), bincorr=G_cor_struct)
    sigma_struct <- bindata::commonprob2sigma(commonprob=cprob)
    
    for (i in 2:num_G) {
      for (j in 1:(i-1)) {
        
        r_G1X <- r_GE[i]
        r_G2X <- r_GE[j]
        w1 <- w_vec[i]
        w2 <- w_vec[j]
        r_G1G2 <- sigma_struct[i, j]
        prob_G1 <- prob_G[i]
        prob_G2 <- prob_G[j]
        
        sig_mat_int <- matrix(data=c(1-r_G1X^2, r_G1G2-r_G1X*r_G2X, r_G1G2-r_G1X*r_G2X, 1-r_G2X^2), nrow=2)
        A <- pracma::quadinf(f=f_G11_G21_X, xa=-Inf, xb=Inf, w1=w1, w2=w2, r_G1X=r_G1X, r_G2X=r_G2X, sig_mat_int=sig_mat_int)$Q[1] 
        sig_mat_int <- matrix(data=c(1-r_G1X^2, -r_G1X*r_G2X, -r_G1X*r_G2X, 1-r_G2X^2), nrow=2)
        B <- pracma::quadinf(f=f_G11_G21_X, xa=-Inf, xb=Inf, w1=w1, w2=w2, r_G1X=r_G1X, r_G2X=r_G2X, sig_mat_int=sig_mat_int)$Q[1] 
        2*A + 2*B - 4*prob_G2*r_G1X*dnorm(w1) - 4*prob_G1*r_G2X*dnorm(w2)
        mu_GEG[i, j] <- 2*A + 2*B - 4*prob_G[j]*r_G1X*dnorm(w1) - 4*prob_G[i]*r_G2X*dnorm(w2)
        mu_GEG[j, i] <- mu_GEG[i, j] 
        
        sig_mat_int <- matrix(data=c(1-r_G1X^2, r_G1G2-r_G1X*r_G2X, r_G1G2-r_G1X*r_G2X, 1-r_G2X^2), nrow=2)
        A <- pracma::quadinf(f=f_G11_G21_XX, xa=-Inf, xb=Inf, w1=w1, w2=w2, r_G1X=r_G1X, r_G2X=r_G2X, sig_mat_int=sig_mat_int)$Q[1] 
        sig_mat_int <- matrix(data=c(1-r_G1X^2, -r_G1X*r_G2X, -r_G1X*r_G2X, 1-r_G2X^2), nrow=2)
        B <- pracma::quadinf(f=f_G11_G21_XX, xa=-Inf, xb=Inf, w1=w1, w2=w2, r_G1X=r_G1X, r_G2X=r_G2X, sig_mat_int=sig_mat_int)$Q[1] 
        mu_GhG[i, j] <- 2*A + 2*B - 4*prob_G2*(r_G1X^2*w1*dnorm(w1)+(1-pnorm(w1))) - 
          4*prob_G1*(r_G2X^2*w2*dnorm(w2)+(1-pnorm(w2))) + 4*prob_G1*prob_G2
        mu_GhG[j, i] <- mu_GhG[i, j] 
        
        mu_GEEG[i, j]  <- mu_GhG[i, j] 
        mu_GEEG[j, i] <- mu_GhG[i, j] 
        
        sig_mat_int <- matrix(data=c(1-r_G1X^2, r_G1G2-r_G1X*r_G2X, r_G1G2-r_G1X*r_G2X, 1-r_G2X^2), nrow=2)
        A <- pracma::quadinf(f=f_G11_G21_XXX, xa=-Inf, xb=Inf, w1=w1, w2=w2, r_G1X=r_G1X, r_G2X=r_G2X, sig_mat_int=sig_mat_int)$Q[1] 
        sig_mat_int <- matrix(data=c(1-r_G1X^2, -r_G1X*r_G2X, -r_G1X*r_G2X, 1-r_G2X^2), nrow=2)
        B <- pracma::quadinf(f=f_G11_G21_XXX, xa=-Inf, xb=Inf, w1=w1, w2=w2, r_G1X=r_G1X, r_G2X=r_G2X, sig_mat_int=sig_mat_int)$Q[1] 
        mu_GEfG[i, j] <- 2*A + 2*B - 4*prob_G2*(r_G1X^3*w1^2*dnorm(w1) - r_G1X^3*dnorm(w1) + 3*r_G1X*dnorm(w1)) - 
          4*prob_G1*(r_G2X^3*w2^2*dnorm(w2) - r_G2X^3*dnorm(w2) + 3*r_G2X*dnorm(w2))
        mu_GEfG[j, i] <- mu_GEfG[i, j]
        mu_GEhG[i, j] <- mu_GEfG[i, j]
        mu_GEhG[j, i] <- mu_GEhG[i, j]
      }
    }
  }

  
  
  ################################################################
  # Now do the higher order matrix terms
  f_G1_E_Z <- function(x, w, r_EZ, r_GE, r_GZ) {
    ( r_EZ * x * surv( (w-x*r_GE) / sqrt(1-r_GE^2) ) + dnorm( (w-r_GE*x) / sqrt(1-r_GE^2) ) * 
        (r_GZ-r_GE*r_EZ) / sqrt(1-r_GE^2) ) * x* dnorm(x)
  }
  
  f_G1_E_W <- function(x, w, r_EW, r_GE, r_GW) {
    ( r_EW * x * surv( (w-x*r_GE) / sqrt(1-r_GE^2) ) + dnorm( (w-r_GE*x) / sqrt(1-r_GE^2) ) * 
        (r_GW-r_GE*r_EW) / sqrt(1-r_GE^2) ) * x* dnorm(x)
  }
  
  f_G1_E_WW <- function(x, w, r_GE, r_GW, r_EW) {
    ( r_EW * x* surv( (w-x*r_GW) / sqrt(1-r_GW^2) ) + dnorm( (w-r_GW*x) / 
                                                               sqrt(1-r_GW^2) ) * (r_GE-r_GW*r_EW) / sqrt(1-r_GW^2) ) * x^2 * dnorm(x)
  }
  
  # E[G1WE^2] requires numerical integration
  f_G1_W_EE <- function(x, w, r_GE, r_GW, r_EW) {
    ( r_EW * x* surv( (w-x*r_GE) / sqrt(1-r_GE^2) ) + dnorm( (w-r_GE*x) / 
                                                               sqrt(1-r_GE^2) ) * (r_GW-r_GE*r_EW) / sqrt(1-r_GE^2) ) * x^2 * dnorm(x)
  }
  
  # E[G1ZE^2] requires numerical integration
  f_G1_Z_EE <- function(x, w, r_GE, r_GZ, r_EZ) {
    ( r_EZ * x* surv( (w-x*r_GE) / sqrt(1-r_GE^2) ) + dnorm( (w-r_GE*x) / 
                                                               sqrt(1-r_GE^2) ) * (r_GZ-r_GE*r_EZ) / sqrt(1-r_GE^2) ) * x^2 * dnorm(x)
  }
  
  # First G and Z terms
  if (num_Z > 0) {
    mu_GEZ <- matrix(data=NA, nrow=num_G, ncol=num_Z)
    mu_ZEG <- matrix(data=NA, nrow=num_Z, ncol=num_G)
    mu_GhZ <- matrix(data=NA, nrow=num_G, ncol=num_Z)
    mu_ZhG <- matrix(data=NA, nrow=num_Z, ncol=num_G)
    
    for (i in 1:num_G) {
     for (j in 1:num_Z) {
        temp_r_EZ <- rho_EZ[j]
        temp_r_GE <- r_GE[i]
        temp_r_GZ <- r_GZ[i, j]
        temp_w <- w_vec[i]
      
        temp_G1_E_Z <- pracma::quadinf(f= f_G1_E_Z, xa=-Inf, xb=Inf, w=temp_w, r_EZ=temp_r_EZ, r_GE=temp_r_GE, r_GZ=temp_r_GZ)$Q
        mu_GEZ[i, j] <- 2*temp_G1_E_Z - 2*prob_G[i]*rho_EZ[j]
        mu_ZEG[j, i] <- 2*temp_G1_E_Z - 2*prob_G[i]*rho_EZ[j]
      
        temp_G1_Z_EE <- pracma::quadinf(f=f_G1_Z_EE, xa=-Inf, xb=Inf, w=temp_w, r_GE=temp_r_GE, r_GZ=temp_r_GZ, r_EZ=temp_r_EZ)$Q
        mu_GhZ[i, j] <- 2*temp_G1_Z_EE
        mu_ZhG[j, i] <- 2*temp_G1_Z_EE
      }
    }
  } else {
    mu_GEZ <- matrix(data=0, nrow=num_G, ncol=1)
    mu_ZEG <- matrix(data=0, nrow=1, ncol=num_G)
    mu_GhZ <- matrix(data=0, nrow=num_G, ncol=1)
    mu_ZhG <- matrix(data=0, nrow=1, ncol=num_G)
  }
  
  # Now G and W terms
  if (num_W > 0) {
    mu_GEW <- matrix(data=NA, nrow=num_G, ncol=num_W)
    mu_WEG <- matrix(data=NA, nrow=num_W, ncol=num_G)
    mu_GhW <- matrix(data=NA, nrow=num_G, ncol=num_W)
    mu_WhG <- matrix(data=NA, nrow=num_W, ncol=num_G)
    mu_GEM <- matrix(data=NA, nrow=num_G, ncol=num_W)
    
    for (i in 1:num_G) {
      for (j in 1:num_W) {
        temp_r_EW <- rho_EW[j]
        temp_r_GE <- r_GE[i]
        temp_r_GW <- r_GW[i, j]
        temp_w <- w_vec[i]
        
        temp_mu_G1_E_W <- pracma::quadinf(f= f_G1_E_W, xa=-Inf, xb=Inf, w=temp_w, r_EW=temp_r_EW, r_GE=temp_r_GE, r_GW=temp_r_GW)$Q
        mu_GEW[i, j] <- 2*temp_mu_G1_E_W	- 2*prob_G[i]*rho_EW[j]	
        mu_WEG[j, i] <- 2*temp_mu_G1_E_W	- 2*prob_G[i]*rho_EW[j]	
        
        temp_mu_G1_E_WW <- pracma::quadinf(f=f_G1_E_WW, xa=-Inf, xb=Inf, w=temp_w, r_GE=temp_r_GE, r_GW=temp_r_GW, r_EW=temp_r_EW)$Q
        mu_GEM[i, j] <- 2*temp_mu_G1_E_WW
        
        temp_mu_G1_W_EE <- pracma::quadinf(f=f_G1_W_EE, xa=-Inf, xb=Inf, w=temp_w, r_GE=temp_r_GE, r_GW=temp_r_GW, r_EW=temp_r_EW)$Q
        mu_GhW[i, j] <- 2*temp_mu_G1_W_EE
        mu_WhG[j, i] <- 2*temp_mu_G1_W_EE
      }
    }
  } else {
    mu_GEW <- matrix(data=0, nrow=num_G, ncol=1)
    mu_WEG <- matrix(data=0, nrow=1, ncol=num_G)
    mu_GhW <- matrix(data=0, nrow=num_G, ncol=1)
    mu_WhG <- matrix(data=0, nrow=1, ncol=num_G)
    mu_GEM <- matrix(data=0, nrow=num_G, ncol=1)
  }
  
  
  # Now put it all into GE_bias_set
  GE_bias_results <- GE_bias(beta_list=beta_list, 
                                 cov_list=list(mu_GE, mu_Gf, mu_Gh,
                                                  mu_EE, mu_Ef, mu_EZ,
                                                  mu_EM, mu_EW, mu_fZ, mu_fW),
                                 cov_mat_list=list(mu_GG, mu_GZ, mu_GM, mu_GW,
                                                   cov_Z, mu_ZW, mu_ZM, cov_W, mu_WM),
                                 mu_list=list(mu_f, mu_h, rep(0, max(1, num_Z)), mu_M, mu_W),
                                 HOM_list=list(mu_GEG, mu_GhG, mu_GEE, mu_GEf, mu_GEh,
                                               mu_GEZ, mu_GEM, mu_GEW, mu_GhW, mu_GhZ, 
                                               mu_GEEG, mu_GEfG, mu_GEhG))
  
  # Do not return \alpha_Z or \alpha_W if no \beta_Z or \beta_M
  if (num_W == 0 & num_Z == 0) {
    GE_bias_results <- GE_bias_results[-c(5,6)]
  } else if (num_Z == 0) {
    GE_bias_results <- GE_bias_results[-5]
  } else if (num_W == 0) {
    GE_bias_results <- GE_bias_results[-6]
  }
  
  # Return 
  return(list(alpha_list=GE_bias_results,
              beta_list = beta_list,
              cov_list=list(mu_GE, mu_Gf, mu_Gh,
                            mu_EE, mu_Ef, mu_EZ,
                            mu_EM, mu_EW, mu_fZ, mu_fW),
              cov_mat_list=list(mu_GG, mu_GZ, mu_GM, mu_GW,
                                cov_Z, mu_ZW, mu_ZM, cov_W, mu_WM),
              mu_list=list(mu_f, mu_h, rep(0, max(1, num_Z)), mu_M, mu_W),
              HOM_list=list(mu_GEG, mu_GhG, mu_GEE, mu_GEf, mu_GEh,
                            mu_GEZ, mu_GEM, mu_GEW, mu_GhW, mu_GhZ, 
                            mu_GEEG, mu_GEfG, mu_GEhG),
              sig_mat=sig_mat))
}