##############################################
# Auxiliary functions for simulation studies #
##############################################
# Methodology of Rebonato and Jackel (2000) to create a
# positive definite matrix out of a non-positive definite matrix.
# Reference: Rebonato and Jackel, "The most general methodology for creating a
# valid correlation matrix for risk management and option pricing purposes",
# Journal of Risk, Vol 2, No 2, 2000
fixNonPositiveDefiniteMatrixIter <- function(origMat) {
cholStatus <- try(u <- chol(origMat), silent = TRUE)
cholError <- ifelse(class(cholStatus) == "try-error", TRUE, FALSE)
newMat <- origMat
iter <- 0
while (cholError) {
iter <- iter + 1
cat("iteration ", iter, "\n")
# replace -ve eigen values with small +ve number
newEig <- eigen(newMat)
newEig2 <- ifelse(newEig$values < 0, 1e-10, newEig$values)
# create modified matrix eqn 5 from Brissette et al 2007,
# inv = transp for eig vectors
newMat <- newEig$vectors %*% diag(newEig2) %*% t(newEig$vectors)
# normalize modified matrix eqn 6 from Brissette et al 2007
newMat <- newMat #/sqrt(diag(newMat) %*% t(diag(newMat)))
# try chol again
cholStatus <- try(u <- chol(newMat), silent = TRUE)
cholError <- ifelse(class(cholStatus) == "try-error", TRUE, FALSE)
}
return(list(newMat, iter))
}
fixNonPositiveDefiniteMatrix <- function(origMat) {
return((fixNonPositiveDefiniteMatrixIter(origMat))[[1]])
}
##########################
# Generating the dataset #
##########################
computeSinglePhenotype <- function(Phi_list, fam.nf, phen.mu, phen.var_g, phen.var_e) {
phen_df <- c()
for (f in 1:length(fam.nf)) {
if (f == 1) {
firstIdx = 1
} else {
firstIdx = sum(fam.nf[1:(f-1)]) + 1
}
lastIdx = firstIdx + fam.nf[f] - 1
CovYf <- kronecker(2*Phi_list[[f]], phen.var_g) + kronecker(diag(1,fam.nf[f]), phen.var_e)
Sigma_f = fixNonPositiveDefiniteMatrix(CovYf)
phen_f <- mvrnorm(1, mu=phen.mu[firstIdx:lastIdx], Sigma=Sigma_f)
phen_f <- matrix(phen_f, byrow=T, ncol=1)
phen_df <- rbind(phen_df, phen_f)
}
return(phen_df)
}
#' @importFrom MASS mvrnorm
generateFamilyDependentSEM <- function(betayx, betazx.y, betazy.x,
varx_g, varx_e, vary.x_g, vary.x_e, varz.xy_g, varz.xy_e,
fam.nf, pedigrees) {
Phi_list <- getKinshipList(pedigrees)
N <- sum(fam.nf)
muX <- rep(0, N) # as.numeric(kronecker(rep(1, nf), phen.mu0))
X <- computeSinglePhenotype(Phi_list, fam.nf, muX, varx_g, varx_e)
muY.X <- betayx * X
Y <- computeSinglePhenotype(Phi_list, fam.nf, muY.X, vary.x_g, vary.x_e)
muZ.XY <- betazx.y * X + betazy.x * Y
Z <- computeSinglePhenotype(Phi_list, fam.nf, muZ.XY, varz.xy_g, varz.xy_e)
return(data.frame(X=X, Y=Y, Z=Z))
}
generateFamilyDependentFCI_SEM <- function(betazx.u, betazu.x, betawy.u, betawu.y,
varx_g, varx_e, varu_g, varu_e, vary_g, vary_e,
varz.xu_g, varz.xu_e, varw.yu_g, varw.yu_e,
fam.nf, pedigrees) {
Phi_list <- getKinshipList(pedigrees)
N <- sum(fam.nf)
muX <- rep(0, N) # as.numeric(kronecker(rep(1, nf), phen.mu0))
X <- computeSinglePhenotype(Phi_list, fam.nf, muX, varx_g, varx_e)
muU <- rep(0, N) # as.numeric(kronecker(rep(1, nf), phen.mu0))
U <- computeSinglePhenotype(Phi_list, fam.nf, muU, varu_g, varu_e)
muY <- rep(0, N) # as.numeric(kronecker(rep(1, nf), phen.mu0))
Y <- computeSinglePhenotype(Phi_list, fam.nf, muY, vary_g, vary_e)
muZ.XU <- betazx.u * X + betazu.x * U
Z <- computeSinglePhenotype(Phi_list, fam.nf, muZ.XU, varz.xu_g, varz.xu_e)
muW.YU <- betawy.u * Y + betawu.y * U
W <- computeSinglePhenotype(Phi_list, fam.nf, muW.YU, varw.yu_g, varw.yu_e)
return(data.frame(X=X, Y=Y, Z=Z, W=W, U=U))
}
getScenarioParametersFromList <- function(parametersTable, rowId) {
betayx <- parametersTable[rowId, 1]
betazx.y <- parametersTable[rowId, 2]
betazy.x <- parametersTable[rowId, 3]
varx_g <- parametersTable[rowId, 4]
varx_e <- parametersTable[rowId, 5]
vary.x_g <- parametersTable[rowId=6]
vary.x_e <- parametersTable[rowId, 7]
varz.xy_g <- parametersTable[rowId, 8]
varz.xy_e <- parametersTable[rowId, 9]
}
getFCIScenarioParameters <- function(scenario) {
if (scenario == 1) {
betazx.u = 0.75
betazu.x = 0.80
betawy.u = 0.75
betawu.y = 0.80
varx_g = 0.875
varx_e = 0.250
varu_g = 0.425
varu_e = 0.625
vary_g = 0.875
vary_e = 0.250
varz.xu_g = 0.250
varz.xu_e = 0.275
varw.yu_g = 0.250
varw.yu_e = 0.275
}
varx = varx_g + varx_e
varu = varu_g + varu_e
vary = vary_g + vary_e
varz.xu <- varz.xu_g + varz.xu_e
varz_g = betazx.u^2 * varx_g + betazu.x^2 * varu_g + varz.xu_g
varz_e = betazx.u^2 * varx_e + betazu.x^2 * varu_e + varz.xu_e
varz = varz_g + varz_e
varw.yu <- varw.yu_g + varw.yu_e
varw_g = betawy.u^2 * vary_g + betawu.y^2 * varu_g + varw.yu_g
varw_e = betawy.u^2 * vary_e + betawu.y^2 * varu_e + varw.yu_e
varw = varw_g + varw_e
covxz = betazx.u * (varx_g + varx_e)
covxz_g = betazx.u * (varx_g)
covxz_e = betazx.u * (varx_e)
rhoxz = covxz/(sqrt(varx)*sqrt(varz))
rhoxz_g = covxz_g/(sqrt(varx_g)*sqrt(varz_g))
rhoxz_e = covxz_e/(sqrt(varx_e)*sqrt(varz_e))
covzw = betazu.x * betawu.y * (varu_g + varu_e)
covzw_g = betazu.x * betawu.y * (varu_g)
covzw_e = betazu.x * betawu.y * (varu_e)
rhozw = covzw/(sqrt(varz)*sqrt(varw))
rhozw_g = covzw_g/(sqrt(varz_g)*sqrt(varw_g))
rhozw_e = covzw_e/(sqrt(varz_e)*sqrt(varw_e))
covzu = betazu.x * (varu_g + varu_e)
covzu_g = betazu.x * (varu_g)
covzu_e = betazu.x * (varu_e)
rhozu = covzu/(sqrt(varz)*sqrt(varu))
rhozu_g = covzu_g/(sqrt(varz_g)*sqrt(varu_g))
rhozu_e = covzu_e/(sqrt(varz_e)*sqrt(varu_e))
covuw = betawu.y * (varu_g + varu_e)
covuw_g = betawu.y * (varu_g)
covuw_e = betawu.y * (varu_e)
rhouw = covuw/(sqrt(varu)*sqrt(varw))
rhouw_g = covuw_g/(sqrt(varu_g)*sqrt(varw_g))
rhouw_e = covuw_e/(sqrt(varu_e)*sqrt(varw_e))
covwy = betawy.u * (vary_g + vary_e)
covwy_g = betawy.u * (vary_g)
covwy_e = betawy.u * (vary_e)
rhowy = covwy/(sqrt(varw)*sqrt(vary))
rhowy_g = covwy_g/(sqrt(varw_g)*sqrt(vary_g))
rhowy_e = covwy_e/(sqrt(varw_e)*sqrt(vary_e))
betaxz = covxz/varz
betauz = covzu/varz
betauw = covuw/varw
betayw = covwy/varw
betawz = covzw/varz
betazw = covzw/varw
covwx.z = - betaxz * covzw - betawz * covxz + betawz * betaxz * varz
covwx.z_g = - betaxz * covzw_g - betawz * covxz_g + betawz * betaxz * varz_g
covwx.z_e = - betaxz * covzw_e - betawz * covxz_e + betawz * betaxz * varz_e
covux.z = - betaxz * covzu - betauz * covxz + betauz * betaxz * varz
covux.z_g = - betaxz * covzu_g - betauz * covxz_g + betauz * betaxz * varz_g
covux.z_e = - betaxz * covzu_e - betauz * covxz_e + betauz * betaxz * varz_e
varx.z = varx - 2 * betaxz * covxz + betaxz^2 * varz
varx.z_g = varx_g - 2 * betaxz * covxz_g + betaxz^2 * varz_g
varx.z_e = varx_e - 2 * betaxz * covxz_e + betaxz^2 * varz_e
varu.z = varu - 2* betauz * covzu + betauz^2 * varz
varu.z_g = varu_g - 2* betauz * covzu_g + betauz^2 * varz_g
varu.z_e = varu_e - 2* betauz * covzu_e + betauz^2 * varz_e
varw.z = varw - 2* betawz * covzw + betawz^2 * varz
varw.z_g = varw_g - 2* betawz * covzw_g + betawz^2 * varz_g
varw.z_e = varw_e - 2* betawz * covzw_e + betawz^2 * varz_e
rhoux.z = covux.z/(sqrt(varx.z) * sqrt(varu.z))
rhoux.z_g = covux.z_g/(sqrt(varx.z_g) * sqrt(varu.z_g))
rhoux.z_e = covux.z_e/(sqrt(varx.z_e) * sqrt(varu.z_e))
covuy.w = - betayw * covuw - betauw * covwy + betauw * betayw * varw
covuy.w_g = - betayw * covuw_g - betauw * covwy_g + betauw * betayw * varw_g
covuy.w_e = - betayw * covuw_e - betauw * covwy_e + betauw * betayw * varw_e
covzy.w = - betayw * covzw - betazw * covwy + betazw * betayw * varw
covzy.w_g = - betayw * covzw_g - betazw * covwy_g + betazw * betayw * varw_g
covzy.w_e = - betayw * covzw_e - betazw * covwy_e + betazw * betayw * varw_e
varu.w = varu - 2* betauw * covuw + betauw^2 * varw
varu.w_g = varu_g - 2* betauw * covuw_g + betauw^2 * varw_g
varu.w_e = varu_e - 2* betauw * covuw_e + betauw^2 * varw_e
varz.w = varz - 2* betazw * covzw + betazw^2 * varw
varz.w_g = varz_g - 2* betazw * covzw_g + betazw^2 * varw_g
varz.w_e = varz_e - 2* betazw * covzw_e + betazw^2 * varw_e
vary.w = vary - 2* betayw * covwy + betayw^2 * varw
vary.w_g = vary_g - 2* betayw * covwy_g + betayw^2 * varw_g
vary.w_e = vary_e - 2* betayw * covwy_e + betayw^2 * varw_e
rhouy.w = covuy.w/(sqrt(vary.w) * sqrt(varu.w))
rhouy.w_g = covuy.w_g/(sqrt(vary.w_g) * sqrt(varu.w_g))
rhouy.w_e = covuy.w_e/(sqrt(vary.w_e) * sqrt(varu.w_e))
rhozy.w = covzy.w/(sqrt(vary.w) * sqrt(varz.w))
rhozy.w_g = covzy.w_g/(sqrt(vary.w_g) * sqrt(varz.w_g))
rhozy.w_e = covzy.w_e/(sqrt(vary.w_e) * sqrt(varz.w_e))
rhowx.z = covwx.z/(sqrt(varx.z) * sqrt(varw.z))
rhowx.z_g = covwx.z_g/(sqrt(varx.z_g) * sqrt(varw.z_g))
rhowx.z_e = covwx.z_e/(sqrt(varx.z_e) * sqrt(varw.z_e))
return(list(betazx.u=betazx.u, betazu.x=betazu.x, betawy.u=betawy.u, betawu.y=betawu.y,
varx_g=varx_g, varx_e=varx_e, varu_g=varu_g, varu_e=varu_e, vary_g=vary_g, vary_e=vary_e,
varz.xu_g=varz.xu_g, varz.xu_e=varz.xu_e, varw.yu_g=varw.yu_g, varw.yu_e=varw.yu_e,
rhoxz=rhoxz, rhoxz_g=rhoxz_g, rhoxz_e=rhoxz_e, #rho(X,Z)
rhozu=rhozu, rhozu_g=rhozu_g, rhozu_e=rhozu_e,
rhouw=rhouw, rhouw_g=rhouw_g, rhouw_e=rhouw_e,
rhowy=rhowy, rhowy_g=rhowy_g, rhowy_e=rhowy_e, #rho(W,Y)
rhowx.z=rhowx.z, rhowx.z_g=rhowx.z_g, rhowx.z_e=rhowx.z_e,
rhoux.z=rhoux.z, rhoux.z_g=rhoux.z_g, rhoux.z_e=rhoux.z_e,
rhouy.w=rhouy.w, rhouy.w_g=rhouy.w_g, rhouy.w_e=rhouy.w_e,
rhozy.w=rhozy.w, rhozy.w_g=rhozy.w_g, rhozy.w_e=rhozy.w_e,
rhozw=rhozw, rhozw_g=rhozw_g, rhozw_e=rhozw_e) #rho(Z,W)
)
}
getScenarioParameters <- function(scenario) {
if (scenario == 1) {
parameters <- c(0.0, 0.75, 0.80, 0.875, 0.250, 0.425, 0.625, 0.250, 0.275)
betayx <- parameters[1]
betazx.y <- parameters[2]
betazy.x <- parameters[3]
varx_g <- parameters[4]
varx_e <- parameters[5]
vary.x_g <- parameters[6]
vary.x_e <- parameters[7]
varz.xy_g <- parameters[8]
varz.xy_e <- parameters[9]
pCor_t <- diag(3)
pCor_t[1,2] = pCor_t[2,1] = -0.5539314
pCor_t[1,3] = pCor_t[3,1] = 0.7392961
pCor_t[2,3] = pCor_t[3,2] = 0.7492686
pCor_g <- diag(3)
pCor_g[1,2] = pCor_g[2,1] = -0.6112120
pCor_g[1,3] = pCor_g[3,1] = 0.8143451
pCor_g[2,3] = pCor_g[3,2] = 0.7218537
pCor_e <- diag(3)
pCor_e[1,2] = pCor_e[2,1] = -0.4949355
pCor_e[1,3] = pCor_e[3,1] = 0.5816751
pCor_e[2,3] = pCor_e[3,2] = 0.7698004
#r12t, r13t, r23t, r12g, r13g, r23g, r12e, r13e, r23e
marg_corr <- c(0.0000000, 0.5880770, 0.6060122,
0.0000000, 0.6966364, 0.5178755,
0.0000000, 0.4152274, 0.7003010)
#r12.3t, r13.2t, r23.1t, r12.3g, r13.2g, r23.1g, r12.3e, r13.2e, r23.1e
part_corr <- c( -0.5539314, 0.7392961, 0.7492686,
-0.6112120, 0.8143451, 0.7218537,
-0.4949355, 0.5816751, 0.7698004)
trueAdjDAG_t <- matrix(0, 3, 3)
trueAdjDAG_t[3,1] <- 1
trueAdjDAG_t[3,2] <- 1
trueAdjUDG_t <- matrix(1, 3, 3)
diag(trueAdjUDG_t) <- 0
trueAdjDAG_g <- trueAdjDAG_t
trueAdjDAG_e <- trueAdjDAG_t
trueAdjUDG_g <- trueAdjUDG_t
trueAdjUDG_e <- trueAdjUDG_t
# > print(c(varx, vary, varz))
# parameters parameters parameters
# 1.125000 1.050000 1.829813
# > print(c(varx.y, varx.z, vary.z, vary.x, varz.x, varz.y))
# parameters parameters parameters parameters parameters parameters
# 1.1250000 0.7359361 0.6643867 1.0500000 1.1970000 1.1578125
# >
# > print(c(varx_g, vary_g, varz_g))
# parameters parameters parameters
# 0.875000 0.425000 1.014188
# > print(c(varx.y_g, varx.z_g, vary.z_g, vary.x_g, varz.x_g, varz.y_g))
# parameters parameters parameters parameters parameters parameters
# 0.8750000 0.4854311 0.3265660 0.4250000 0.5220000 0.7421875
# >
# > print(c(varx_e, vary_e, varz_e))
# parameters parameters parameters
# 0.250000 0.625000 0.815625
# > print(c(varx.y_e, varx.z_e, vary.z_e, vary.x_e, varz.x_e, varz.y_e))
# parameters parameters parameters parameters parameters parameters
# 0.2500000 0.2505050 0.3378208 0.6250000 0.6750000 0.4156250
# >
# > print(c(rhoyx, rhoxz, rhoyz))
# parameters parameters parameters
# 0.0000000 0.5880770 0.6060122
# > print(c(rhoyx_g, rhoxz_g, rhoyz_g))
# parameters parameters parameters
# 0.0000000 0.6966364 0.5178755
# > print(c(rhoyx_e, rhoxz_e, rhoyz_e))
# parameters parameters parameters
# 0.0000000 0.4152274 0.7003010
# >
# > print(c(varx_g/varx, vary_g/vary, varz_g/varz))
# parameters parameters parameters
# 0.7777778 0.4047619 0.5542576
# > print(c(rhoyx.z, rhozx.y, rhozy.x))
# parameters parameters parameters
# -0.5539314 0.7392961 0.7492686
# > print(c(rhoyx.z_g, rhozx.y_g, rhozy.x_g))
# parameters parameters parameters
# -0.6112120 0.8143451 0.7218537
# > print(c(rhoyx.z_e, rhozx.y_e, rhozy.x_e))
# parameters parameters parameters
# -0.4949355 0.5816751 0.7698004
} else if (scenario == 2) {
parameters <- c(-0.8, 0.80, 0.00, 0.425, 0.600, 0.250, 0.725, 0.475, 0.275)
betayx <- parameters[1]
betazx.y <- parameters[2]
betazy.x <- parameters[3]
varx_g <- parameters[4]
varx_e <- parameters[5]
vary.x_g <- parameters[6]
vary.x_e <- parameters[7]
varz.xy_g <- parameters[8]
varz.xy_e <- parameters[9]
pCor_t <- diag(3)
pCor_t[1,2] = pCor_t[2,1] = -5.139177e-01
pCor_t[1,3] = pCor_t[3,1] = 5.859564e-01
pCor_t[2,3] = pCor_t[3,2] = -2.596613e-16
pCor_g <- diag(3)
pCor_g[1,2] = pCor_g[2,1] = -0.6478292
pCor_g[1,3] = pCor_g[3,1] = 0.4739953
pCor_g[2,3] = pCor_g[3,2] = 0.0000000
pCor_e <- diag(3)
pCor_e[1,2] = pCor_e[2,1] = -0.4349207
pCor_e[1,3] = pCor_e[3,1] = 0.6932896
pCor_e[2,3] = pCor_e[3,2] = 0.0000000
#r12t, r13t, r23t, r12g, r13g, r23g, r12e, r13e, r23e
marg_corr <- c( -0.6341981, 0.6830606, -0.4331958,
-0.7218537, 0.6034262, -0.4355854,
-0.5884368, 0.7633486, -0.4491824)
#r12.3t, r13.2t, r23.1t, r12.3g, r13.2g, r23.1g, r12.3e, r13.2e, r23.1e
part_corr <- c( -5.139177e-01, 5.859564e-01, -2.596613e-16,
-0.6478292, 0.4739953, 0.0000000,
-0.4349207, 0.6932896, 0.0000000)
trueAdjDAG_t <- matrix(0, 3, 3)
trueAdjDAG_t[2,1] <- 1
trueAdjDAG_t[3,1] <- 1
trueAdjUDG_t <- matrix(0, 3, 3)
trueAdjDAG_t[1,2] <- 1
trueAdjDAG_t[1,3] <- 1
trueAdjDAG_t[2,1] <- 1
trueAdjDAG_t[3,1] <- 1
trueAdjDAG_g <- trueAdjDAG_t
trueAdjDAG_e <- trueAdjDAG_t
trueAdjUDG_g <- trueAdjUDG_t
trueAdjUDG_e <- trueAdjUDG_t
# print(c(varx, vary, varz))
# parameters parameters parameters
# 1.025 1.631 1.406
# > print(c(varx.y, varx.z, vary.z, vary.x, varz.x, varz.y))
# parameters parameters parameters parameters parameters parameters
# 0.6127376 0.5467639 1.3249289 0.9750000 0.7500000 1.1421521
# >
# > print(c(varx_g, vary_g, varz_g))
# parameters parameters parameters
# 0.425 0.522 0.747
# > print(c(varx.y_g, varx.z_g, vary.z_g, vary.x_g, varz.x_g, varz.y_g))
# parameters parameters parameters parameters parameters parameters
# 0.2150680 0.2824981 0.4307988 0.2500000 0.4750000 0.6126435
# >
# > print(c(varx_e, vary_e, varz_e))
# parameters parameters parameters
# 0.600 1.109 0.659
# > print(c(varx.y_e, varx.z_e, vary.z_e, vary.x_e, varz.x_e, varz.y_e))
# parameters parameters parameters parameters parameters parameters
# 0.3976695 0.2642657 0.8941301 0.7250000 0.2750000 0.5295085
# >
# > print(c(rhoyx, rhoxz, rhoyz))
# parameters parameters parameters
# -0.6341981 0.6830606 -0.4331958
# > print(c(rhoyx_g, rhoxz_g, rhoyz_g))
# parameters parameters parameters
# -0.7218537 0.6034262 -0.4355854
# > print(c(rhoyx_e, rhoxz_e, rhoyz_e))
# parameters parameters parameters
# -0.5884368 0.7633486 -0.4491824
# >
# > print(c(varx_g/varx, vary_g/vary, varz_g/varz))
# parameters parameters parameters
# 0.4146341 0.3200490 0.5312945
# > print(c(rhoyx.z, rhozx.y, rhozy.x))
# parameters parameters parameters
# -5.139177e-01 5.859564e-01 -2.596613e-16
# > print(c(rhoyx.z_g, rhozx.y_g, rhozy.x_g))
# parameters parameters parameters
# -0.6478292 0.4739953 0.0000000
# > print(c(rhoyx.z_e, rhozx.y_e, rhozy.x_e))
# parameters parameters parameters
# -0.4349207 0.6932896 0.0000000
} else if (scenario == 3) {
parameters <- c(-0.8, -0.9, 0.9, 0.6, 0.55, 0.925, 0.3, 0.25, 0.975)
betayx <- parameters[1]
betazx.y <- parameters[2]
betazy.x <- parameters[3]
varx_g <- parameters[4]
varx_e <- parameters[5]
vary.x_g <- parameters[6]
vary.x_e <- parameters[7]
varz.xy_g <- parameters[8]
varz.xy_e <- parameters[9]
pCor_t <- diag(3)
pCor_t[1,2] = pCor_t[2,1] = 0.004686763
pCor_t[1,3] = pCor_t[3,1] = -0.567485469
pCor_t[2,3] = pCor_t[3,2] = 0.668964732
pCor_g <- diag(3)
pCor_g[1,2] = pCor_g[2,1] = 0.4244711
pCor_g[1,3] = pCor_g[3,1] = -0.7658174
pCor_g[2,3] = pCor_g[3,2] = 0.8659171
pCor_e <- diag(3)
pCor_e[1,2] = pCor_e[2,1] = -0.4077895
pCor_e[1,3] = pCor_e[3,1] = -0.4348273
pCor_e[2,3] = pCor_e[3,2] = 0.4466625
#r12t, r13t, r23t, r12g, r13g, r23g, r12e, r13e, r23e
marg_corr <- c(-0.6126326, -0.7592639, 0.8092371,
-0.5416214, -0.7821601, 0.8771786,
-0.7347634, -0.7364439, 0.7460787)
#r12.3t, r13.2t, r23.1t, r12.3g, r13.2g, r23.1g, r12.3e, r13.2e, r23.1e
part_corr <- c( 0.004686763, -0.567485469, 0.668964732,
0.4244711, -0.7658174, 0.8659171,
-0.4077895, -0.4348273, 0.4466625)
# notation col -> row
trueAdjDAG_t <- matrix(0, 3, 3)
trueAdjDAG_t[3,1] <- 1
trueAdjDAG_t[2,1] <- 1
trueAdjDAG_t[3,2] <- 1
trueAdjUDG_t <- matrix(1, 3, 3)
diag(trueAdjUDG_t) <- 0
trueAdjDAG_g <- trueAdjDAG_t
trueAdjDAG_e <- trueAdjDAG_t
trueAdjUDG_g <- trueAdjUDG_t
trueAdjUDG_e <- trueAdjUDG_t
# > print(c(varx, vary, varz))
# [1] 1.15000 1.96100 5.23531
# > print(c(varx.y, varx.z, vary.z, vary.x, varz.x, varz.y))
# [1] 0.7183835 0.4870461 0.6768104 1.2250000 2.2172500 1.8068906
# >
# > print(c(varx_g, vary_g, varz_g))
# [1] 0.60000 1.30900 2.57389
# > print(c(varx.y_g, varx.z_g, vary.z_g, vary.x_g, varz.x_g, varz.y_g))
# [1] 0.4377287 0.2341570 0.3454861 0.9250000 0.9992500 0.6045603
# >
# > print(c(varx_e, vary_e, varz_e))
# [1] 0.55000 0.65200 2.66142
# > print(c(varx.y_e, varx.z_e, vary.z_e, vary.x_e, varz.x_e, varz.y_e))
# [1] 0.2806548 0.2528892 0.3313243 0.3000000 1.2180000 1.2023304
# >
# > print(c(rhoyx, rhoxz, rhoyz))
# [1] -0.6126326 -0.7592639 0.8092371
# > print(c(rhoyx_g, rhoxz_g, rhoyz_g))
# [1] -0.5416214 -0.7821601 0.8771786
# > print(c(rhoyx_e, rhoxz_e, rhoyz_e))
# [1] -0.7347634 -0.7364439 0.7460787
# >
# > print(c(varx_g/varx, vary_g/vary, varz_g/varz))
# [1] 0.5217391 0.6675166 0.4916404
# > print(c(rhoyx.z, rhozx.y, rhozy.x))
# [1] 0.004686763 -0.567485469 0.668964732
# > print(c(rhoyx.z_g, rhozx.y_g, rhozy.x_g))
# [1] 0.4244711 -0.7658174 0.8659171
# > print(c(rhoyx.z_e, rhozx.y_e, rhozy.x_e))
# [1] -0.4077895 -0.4348273 0.4466625
} else if (scenario == 4) {
###################################################
# Defining zero genetic variability for X. #
# There is an unshielded collider in the #
# environmental network and only the edge #
# Y->Z in the genetic network. #
# The total network han an unshielded collider. #
# Note that the effect for the rhoyx.z_e is #
# negative even though there is no negative coef. #
# Zero heritability for X and high heritability #
# for Y and Z. #
###################################################
betayx = 0
betazx.y = 0.25
betazy.x = 0.55
varx_e = 2.5
varx_g = 0
vary.x_e <- 0.75
vary.x_g <- 1.75
varz.xy_e <- 0.25
varz.xy_g <- 0.75
pCor_t <- diag(3)
pCor_t[1,2] = pCor_t[2,1] = -0.2412257
pCor_t[1,3] = pCor_t[3,1] = 0.3676073
pCor_t[2,3] = pCor_t[3,2] = 0.6562050
pCor_g <- diag(3)
pCor_g[1,2] = pCor_g[2,1] = -0.03723522
pCor_g[1,3] = pCor_g[3,1] = 0.00000000
pCor_g[2,3] = pCor_g[3,2] = 0.64325443
pCor_e <- diag(3)
pCor_e[1,2] = pCor_e[2,1] = -0.4258620
pCor_e[1,3] = pCor_e[3,1] = 0.6201737
pCor_e[2,3] = pCor_e[3,2] = 0.6897489
#r12t, r13t, r23t, r12g, r13g, r23g, r12e, r13e, r23e
marg_corr <- c( 0.0000000, 0.2858310, 0.6288281,
0.0000000, 0.0000000, 0.6432544,
0.0000000, 0.4967813, 0.5986164)
#r12.3t, r13.2t, r23.1t, r12.3g, r13.2g, r23.1g, r12.3e, r13.2e, r23.1e
part_corr <- c( -0.2412257, 0.3676073, 0.6562050,
-0.03723522, 0.00000000, 0.64325443,
-0.4258620, 0.6201737, 0.6897489)
trueAdjDAG_t <- matrix(0, 3, 3)
trueAdjDAG_t[3,1] <- 1
trueAdjDAG_t[3,2] <- 1
trueAdjDAG_g <- trueAdjDAG_t
trueAdjDAG_g[3,1] <- 0
trueAdjDAG_e <- trueAdjDAG_t
trueAdjUDG_t <- matrix(1, 3, 3)
diag(trueAdjUDG_t) <- 0
trueAdjUDG_g <- trueAdjDAG_g
trueAdjDAG_t[2,3] <- 1
trueAdjUDG_e <- trueAdjUDG_t
# print(c(varx, vary, varz))
# 2.5000 2.5000 1.9125
# print(c(varx.y, varx.z, vary.z, vary.x, varz.x, varz.y))
# 2.500000 2.295752 1.511438 2.500000 1.756250 1.156250
# print(c(varx_g, vary_g, varz_g))
# 0.000000 1.750000 1.279375
# print(c(varx.y_g, varx.z_g, vary.z_g, vary.x_g, varz.x_g, varz.y_g))
# 0.0000000 0.1366328 1.0273159 1.7500000 1.2793750 0.7500000
# print(c(varx_e, vary_e, varz_e))
# 2.500000 0.750000 0.633125
# print(c(varx.y_e, varx.z_e, vary.z_e, vary.x_e, varz.x_e, varz.y_e))
# 2.500000 2.159119 0.484122 0.750000 0.476875 0.406250
# print(c(rhoyx, rhoxz, rhoyz))
# 0.0000000 0.2858310 0.6288281
# print(c(rhoyx_g, rhoxz_g, rhoyz_g))
# 0.0000000 0.0000000 0.6432544
# print(c(rhoyx_e, rhoxz_e, rhoyz_e))
# 0.0000000 0.4967813 0.5986164
# print(c(varx_g/varx, vary_g/vary, varz_g/varz))
# 0.0000000 0.7000000 0.6689542
# print(c(rhoyx.z, rhozx.y, rhozy.x))
# -0.2412257 0.3676073 0.6562050
# print(c(rhoyx.z_g, rhozx.y_g, rhozy.x_g))
# -0.03723522 0.00000000 0.64325443
# print(c(rhoyx.z_e, rhozx.y_e, rhozy.x_e))
# -0.4258620 0.6201737 0.6897489
} else {
return(NULL)
}
return(list(betayx=betayx, betazx.y=betazx.y, betazy.x=betazy.x,
varx_g=varx_g, varx_e=varx_e,
vary.x_g=vary.x_g, vary.x_e=vary.x_e,
varz.xy_g=varz.xy_g, varz.xy_e=varz.xy_e,
pCor_t=pCor_t, pCor_g=pCor_g, pCor_e=pCor_e,
marg_corr=marg_corr, part_corr=part_corr,
DAG_t=trueAdjDAG_t, DAG_g=trueAdjDAG_g, DAG_e=trueAdjDAG_e,
UDG_t=trueAdjUDG_t, UDG_g=trueAdjUDG_g, UDG_e=trueAdjUDG_e))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.