#' 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))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.