Permalink
Cannot retrieve contributors at this time
Join GitHub today
GitHub is home to over 50 million developers working together to host and review code, manage projects, and build software together.
Sign up| #' Scale phenotype component. | |
| #' | |
| #' The function scales the specified component such that the average column | |
| #' variance is equal to the user-specified proportion of variance. | |
| #' | |
| #' @param component [N x P] Phenotype matrix [double] where [N] are the number | |
| #' of samples and P the number of phenotypes | |
| #' @param propvar Number [double] specifying the proportion of variance that | |
| #' should be explained by this phenotype component | |
| #' @return If propvar != 0, a named list with the [N x P] matrix of the scaled | |
| #' component (component) and its scale factor [double] (scale_factor) else | |
| #' returns NULL | |
| #' @export | |
| #' @examples | |
| #' x <- matrix(rnorm(100), nc=10) | |
| #' x_scaled <- rescaleVariance(x, propvar=0.4) | |
| rescaleVariance <- function(component, propvar) { | |
| if (!is.numeric(propvar)) { | |
| stop("propvar needs to be numeric") | |
| } | |
| if (propvar < 0 || propvar > 1) { | |
| stop("propvar cannot be less than 0 and or greater than 1") | |
| } | |
| if (!is.null(component) && !is.matrix(component)){ | |
| stop("component needs to be a matrix") | |
| } | |
| if (!is.null(component) && propvar != 0) { | |
| var_component <- var(component) | |
| mean_var <- mean(diag(var_component)) | |
| scale_factor <- sqrt(propvar/mean_var) | |
| component_scaled <- component * scale_factor | |
| colnames(component_scaled) <- colnames(component) | |
| rownames(component_scaled) <- rownames(component) | |
| return(list(component=component_scaled, | |
| scale_factor=scale_factor)) | |
| } else { | |
| return(NULL) | |
| } | |
| } | |
| #' Phenotype transformation. | |
| #' | |
| #' Transformation of phenotype component by applying a user-specified | |
| #' non-linear transformation to the phenotype component. | |
| #' | |
| #' @param component [N x P] Phenotype matrix [double] where [N] are the number | |
| #' of samples and P the number of phenotypes | |
| #' @param method [string] one of exp (exponential), log (logarithm), poly | |
| #' (polynomial), sqrt (squareroot) or custom (user-supplied function) | |
| #' @param transformNeg [string] one of abs (absolute value) or set0 (set all | |
| #' negative values to zero). If method==log and transformNeg==set0, negative | |
| #' values set to 1e-5 | |
| #' @param alpha [double] weighting scalar for non-linearity: alpha==0 fully | |
| #' linear phenotype, alpha==1 fully non-linear phenotype. See @details. | |
| #' @param logbase [int] when method==log, sets the log base for transformation | |
| #' @param expbase [double] when method==exp, sets the exp base for | |
| #' transformation. | |
| #' @param power [double] when method==poly, sets the power to raise to. | |
| #' @param f [function] function accepting component as a single argument. | |
| #' @param verbose [boolean]; If TRUE, progress info is printed to standard out. | |
| #' @details transformNonlinear takes a phenotype component as input and | |
| #' transforms it according to the specified transformation method. The user can | |
| #' choose how strongly non-linear the resulting phenotype component should be, | |
| #' by specifying the weighting parameter alpha: | |
| #' component_transformed = (1 - alpha) \* component + | |
| #' alpha \* transformfunction(component) | |
| #' @return [N x P] transformed phenotype matrix [double] | |
| #' @export | |
| #' @examples | |
| #' # Simulate non-genetic covariate effects | |
| #' cov_effects <- noiseFixedEffects(N=100, P=5) | |
| #' # Transform logarithmically | |
| #' covs_log <- transformNonlinear(cov_effects$shared, alpha=0.5, method="log", | |
| #' transformNeg="abs") | |
| #' # Transform custom | |
| #' f_custom <- function(x) {x^2 + 3*x} | |
| #' covs_custom <- transformNonlinear(cov_effects$shared, alpha=0.5, | |
| #' method="custom", f=f_custom) | |
| transformNonlinear <- function(component, alpha, method, logbase=10, power=2, | |
| expbase=NULL, transformNeg="abs", f=NULL, | |
| verbose=TRUE) { | |
| testNumerics(numbers=c(alpha, expbase, logbase), proportions=alpha, | |
| positives=logbase) | |
| nonlinear <- function(x, method, expbase, logbase, power, f) { | |
| if (!is.null(transformNeg)) { | |
| if (transformNeg == "abs") { | |
| x <- abs(x) | |
| } else if (transformNeg == "set0") { | |
| if (method == "log") { | |
| x[x < 0] <- 1e-5 | |
| } else { | |
| x[x < 0] <- 0 | |
| } | |
| } else { | |
| stop("Negative transformation method not known") | |
| } | |
| } | |
| vmessage(c("Use", method, "as transformation method"), verbose=verbose) | |
| if (method == "exp") { | |
| if (is.null(expbase)) y <- exp(x) | |
| if (!is.null(expbase)) y <- expbase^x | |
| } else if (method == "log") { | |
| if (any(x <= 0)) { | |
| stop(paste("For transformation log, all values have to", | |
| " be greater than zero")) | |
| } | |
| y <- log(x, base=logbase) | |
| } else if (method == "poly") { | |
| y <- x^power | |
| } else if (method == "sqrt") { | |
| if (any(x < 0)) { | |
| stop(paste("For transformation square root, all values have to", | |
| "be greater or equal to zero")) | |
| } | |
| y <- sqrt(x) | |
| } else if (method == "custom") { | |
| y <- f(x) | |
| } else { | |
| stop("Transformation method not known") | |
| } | |
| return(y) | |
| } | |
| component_trans <- (1 - alpha) * component + | |
| alpha * nonlinear(component, method=method, logbase=logbase, | |
| expbase=expbase, power=power, f=f) | |
| return(component_trans) | |
| } | |
| #' Set simulation model. | |
| #' | |
| #' Based on parameters provided, this function sets the name for the phenotype | |
| #' simulation. It carries out compatibiltiy checks of the specifie parameters | |
| #' and checks for any missing information. | |
| #' | |
| #' @param genVar Total genetic variance [double]. | |
| #' @param h2s Proportion [double] of variance of genetic variant effects. | |
| #' @param h2bg Proportion [double] of variance of infinitesimal genetic effects | |
| #' i.e. correlation introduced by sample kinship). | |
| #' @param theta Proportion [double] of variance of shared genetic variant | |
| #' effects. | |
| #' @param eta Proportion [double] of variance of shared infinitesimal genetic | |
| #' effects. | |
| #' @param noiseVar Total noise variance [double]. | |
| #' @param rho Proportion [double] of variance of correlated noise effects. | |
| #' @param pcorr Correlation [double] between phenotypes. | |
| #' @param delta Proportion [double] of variance of non-genetic covariate effect. | |
| #' @param gamma Proportion [double] of variance of shared non-genetic covariate | |
| #' effects. | |
| #' @param phi Proportion [double] of variance of observational noise effects. | |
| #' @param alpha Proportion [double] of variance of shared observational noise | |
| #' effect. | |
| #' @param pIndependentConfounders Proportion [double] of non-genetic covariate | |
| #' to have a trait-independent effect. | |
| #' @param pTraitIndependentConfounders Proportion [double] of traits influenced | |
| #' by independent non-genetic covariate effects. | |
| #' @param pIndependentGenetic Proportion [double] of genetic variant effects to | |
| #' have a trait-independent fixed effect. | |
| #' @param pTraitIndependentGenetic Proportion [double] of traits influenced by | |
| #' independent genetic variant effects. | |
| #' @param proportionNonlinear [double] proportion of the phenotype to be non- | |
| #' linear | |
| #' @param cNrSNP Number [integer] of causal SNPs; used as genetic variant | |
| #' effects. | |
| #' @param NrConfounders Number [integer] of non-genetic covariates; used as | |
| #' non-genetic covariate effects. | |
| #' @param verbose [boolean]; If TRUE, progress info is printed to standard out. | |
| #' @return Named list containing the genetic model (modelGenetic), the noise | |
| #' model (modelNoise) and the input parameters (h2s, h2bg, noiseVar, rho, delta, | |
| #' phi, gamma, theta, eta, alpha, pcorr, proportionNonlinear). Model options | |
| #' are: modelNoise: "noNoise", "noiseFixedOnly", "noiseBgOnly", | |
| #' "noiseCorrelatedOnly", "noiseFixedAndBg","noiseCorrelatedAndBg", | |
| #' "noiseFixedAndCorrelated", "noiseFixedAndBgAndCorrelated" | |
| #' modelGenetic: "noGenetic","geneticBgOnly", "geneticFixedOnly", | |
| #' "geneticFixedAndBg" | |
| #' @export | |
| #' @examples | |
| #' #genetic fixed effects only | |
| #' model <- setModel(genVar=1, h2s=1) | |
| #' | |
| #' #genetic fixed and bg effects | |
| #' model <- setModel(genVar=1, h2s=0.01) | |
| #' | |
| #' #genetic and noise fixed effects only | |
| #' model <- setModel(genVar=0.4, h2s=1, delta=1) | |
| setModel <- function(genVar=NULL, h2s=NULL, theta=0.8, h2bg=NULL, eta=0.8, | |
| noiseVar=NULL, delta=NULL, gamma=0.8, rho=NULL, phi=NULL, | |
| alpha=0.8, pcorr=0.6, pIndependentConfounders=0.4, | |
| pTraitIndependentConfounders=0.2, pIndependentGenetic=0.4, | |
| pTraitIndependentGenetic=0.2, proportionNonlinear=0, | |
| cNrSNP=NULL, NrConfounders=10, | |
| verbose=TRUE) { | |
| if (is.null(c(genVar, noiseVar, h2bg, h2s, delta, rho, phi))) { | |
| stop("No variance components specified") | |
| } | |
| if (!is.null(NrConfounders) && all(NrConfounders != 0)) { | |
| # for testing purposes set to anything other integer than 0 | |
| NrConfounders=10 | |
| } | |
| numbers <- list(genVar=genVar, h2s=h2s, h2bg=h2bg, theta=theta, | |
| eta=eta, noiseVar=noiseVar, delta=delta, gamma=gamma, | |
| rho=rho, phi=phi, alpha=alpha, pcorr=pcorr, | |
| pIndependentConfounders=pIndependentConfounders, | |
| pTraitIndependentConfounders= | |
| pTraitIndependentConfounders, | |
| pIndependentGenetic=pIndependentGenetic, | |
| pTraitIndependentGenetic=pTraitIndependentGenetic, | |
| proportionNonlinear=proportionNonlinear, | |
| cNrSNP=cNrSNP, NrConfounders=NrConfounders) | |
| proportions <- list(genVar=genVar, h2s=h2s, h2bg=h2bg, theta=theta, | |
| eta=eta, noiseVar=noiseVar, delta=delta, gamma=gamma, | |
| rho=rho, phi=phi, alpha=alpha, pcorr=pcorr, | |
| pIndependentConfounders=pIndependentConfounders, | |
| pTraitIndependentConfounders= | |
| pTraitIndependentConfounders, | |
| pIndependentGenetic=pIndependentGenetic, | |
| pTraitIndependentGenetic=pTraitIndependentGenetic, | |
| proportionNonlinear=proportionNonlinear) | |
| testNumerics(numbers=numbers, proportions=proportions) | |
| if (is.null(genVar)) { | |
| if (is.null(noiseVar)) { | |
| stop(paste("Neither genVar nor noiseVar are provided, thus", | |
| "proportion of variance from genetics cannot be deduced") | |
| ) | |
| } else { | |
| genVar <- 1 - noiseVar | |
| } | |
| } | |
| if (is.null(noiseVar)) { | |
| noiseVar <- 1 - genVar | |
| } | |
| if (noiseVar + genVar > 1) { | |
| stop("Sum of genetic variance and noise variance is greater than 1") | |
| } | |
| if (noiseVar == 0 && any(!is.null(phi), !is.null(rho), !is.null(delta))) { | |
| stop(paste("The noise variance is set to 0 (or genetic variance set to", | |
| " 1) but noise variance components are supplied") | |
| ) | |
| } | |
| if (genVar == 0 && any(!is.null(h2s), !is.null(h2bg))) { | |
| stop(paste("The genetic variance is set to 0 (or noise variance set to", | |
| " 1) but genetic variance components are supplied") | |
| ) | |
| } | |
| vmessage(c("The total noise variance (noiseVar) is:", noiseVar), | |
| verbose=verbose) | |
| if ((noiseVar) == 0 ) { | |
| modelNoise="noNoise" | |
| delta <- 0 | |
| rho <- 0 | |
| phi <- 0 | |
| vmessage(c("The noise model is:", modelNoise), verbose=verbose) | |
| } else { | |
| if (all(c(is.null(delta), is.null(rho), is.null(phi)))) { | |
| stop(paste("Neither delta (non-genetic covariate effect variance)", | |
| "nor rho (correlated noise effect variance) or phi", | |
| "(observational noise effect variance) are provided, at", | |
| "least one is required") | |
| ) | |
| } | |
| if (length(c(delta, rho, phi)) >=2) { | |
| if (sum(c(delta, rho, phi)) > 1 ) { | |
| stop(paste("Sum of the proportion of the variance of noise", | |
| "effects is greater than 1; change noiseVar, delta", | |
| "(non-genetic covariate effect variance),", | |
| "rho (correlated noise effect variance) or phi", | |
| "(observational noise effect variance) such that", | |
| "delta + rho + phi = 1") | |
| ) | |
| } | |
| if (length(c(delta, rho, phi)) == 3 && | |
| sum(c(delta, rho, phi)) < 1) { | |
| stop(paste("Sum of the proportion of the variance of noise", | |
| "effects is less than 1; change noiseVar, delta", | |
| "(non-genetic covariate effect variance),", | |
| "rho (correlated noise effect variance) or phi", | |
| "(observational noise effect variance) such that", | |
| "delta + rho + phi = 1") | |
| ) | |
| } | |
| } | |
| if (is.null(phi) && all(c(!is.null(delta), !is.null(rho)))) { | |
| phi <- 1 - delta - rho | |
| } | |
| if (is.null(rho) && all(c(!is.null(delta), !is.null(phi)))) { | |
| rho <- 1 - delta - phi | |
| } | |
| if (is.null(delta) && all(c(!is.null(rho), !is.null(phi)))) { | |
| delta <- 1 - rho - phi | |
| } | |
| if (is.null(delta)) { | |
| delta <- 0 | |
| } | |
| if (is.null(phi)) { | |
| phi <- 0 | |
| } | |
| if (is.null(rho)) { | |
| rho <- 0 | |
| } | |
| if(delta + rho + phi != 1) { | |
| stop(paste("Not enough components provided to set proportions of", | |
| "noise variance correctly; if noise variance is only", | |
| "explained by one component, its proportion of variance", | |
| "needs to be set to 1; otherwise, the proportions of", | |
| "variance of at least 2 components need to be specified") | |
| ) | |
| } | |
| if (delta != 0 && NrConfounders == 0) { | |
| stop(paste("Proportion of of non-genetic covariate variance ", | |
| "(delta) is", delta, "but number of", | |
| "NrConfounder is set to zero")) | |
| } | |
| if (gamma == 1) { | |
| pIndependentConfounders=0 | |
| pTraitIndependentConfounders=0 | |
| } | |
| if (gamma == 0) { | |
| pIndependentConfounders=1 | |
| pTraitIndependentConfounders=1 | |
| } | |
| if (phi == 1) { | |
| modelNoise="noiseBgOnly" | |
| vmessage(c("The noise model is:", modelNoise), verbose=verbose) | |
| vmessage(c("Proportion of observational noise variance (phi):", | |
| phi), verbose=verbose) | |
| vmessage(c("Variance of shared observational noise effect (alpha):", | |
| alpha), verbose=verbose) | |
| } else if (delta == 1) { | |
| modelNoise="noiseFixedOnly" | |
| vmessage(c("The noise model is:", modelNoise), verbose=verbose) | |
| vmessage(c("Proportion of non-genetic covariate variance (delta):", | |
| delta), verbose=verbose) | |
| vmessage(c("Proportion of variance of shared non-genetic covariate", | |
| "effects (gamma):", gamma), verbose=verbose) | |
| vmessage(c("Proportion of non-genetic covariates to have", | |
| "a trait-independent effect (pIndependentConfounders" | |
| , "):", pIndependentConfounders), | |
| verbose=verbose) | |
| vmessage(c("Proportion of traits influenced by independent", | |
| "non-genetic covariate effects", | |
| "(pTraitIndependentConfounders):", | |
| pTraitIndependentConfounders), verbose=verbose) | |
| } else if (rho == 1) { | |
| modelNoise="noiseCorrelatedOnly" | |
| vmessage(c("The noise model is:", modelNoise), verbose=verbose) | |
| vmessage(c("Proportion of variance of correlated noise effects", | |
| "(rho):", rho), verbose=verbose) | |
| } else if (rho == 0) { | |
| modelNoise="noiseFixedAndBg" | |
| vmessage(c("The noise model is:", modelNoise), verbose=verbose) | |
| vmessage(c("Proportion of non-genetic covariate variance (delta):", | |
| delta), verbose=verbose) | |
| vmessage(c("Proportion of variance of shared non-genetic covariate", | |
| "effects (gamma):", gamma), verbose=verbose) | |
| vmessage(c("Proportion of non-genetic covariates to have", | |
| "a trait-independent effect (pIndependentConfounders" | |
| , "):", pIndependentConfounders), | |
| verbose=verbose) | |
| vmessage(c("Proportion of traits influenced by independent", | |
| "non-genetic covariate effects", | |
| "(pTraitIndependentConfounders):", | |
| pTraitIndependentConfounders), verbose=verbose) | |
| vmessage(c("Proportion of observational noise variance (phi):", | |
| phi), verbose=verbose) | |
| vmessage(c("Variance of shared observational noise effect (alpha):", | |
| alpha), verbose=verbose) | |
| } else if (phi == 0) { | |
| modelNoise="noiseFixedAndCorrelated" | |
| vmessage(c("The noise model is:", modelNoise), verbose=verbose) | |
| vmessage(c("Proportion of non-genetic covariate variance (delta):", | |
| delta), verbose=verbose) | |
| vmessage(c("Proportion of variance of shared non-genetic covariate", | |
| "effects (gamma):", gamma), verbose=verbose) | |
| vmessage(c("Proportion of non-genetic covariates to have", | |
| "a trait-independent effect (pIndependentConfounders" | |
| , "):", pIndependentConfounders), | |
| verbose=verbose) | |
| vmessage(c("Proportion of traits influenced by independent", | |
| "non-genetic covariate effects", | |
| "(pTraitIndependentConfounders):", | |
| pTraitIndependentConfounders), verbose=verbose) | |
| vmessage(c("Proportion of variance of correlated noise effects", | |
| "(rho):", rho), verbose=verbose) | |
| vmessage(c("Correlation between phenotypes (pcorr):", pcorr), | |
| verbose=verbose) | |
| } else if (delta == 0 ) { | |
| modelNoise="noiseCorrelatedAndBg" | |
| vmessage(c("The noise model is:", modelNoise), verbose=verbose) | |
| vmessage(paste("Proportion of variance of correlated noise effects", | |
| "(rho):", rho), verbose=verbose) | |
| vmessage(c("Proportion of observational noise variance (phi):", | |
| phi), verbose=verbose) | |
| vmessage(c("Variance of shared observational noise effect (alpha):", | |
| alpha), verbose=verbose) | |
| } else { | |
| modelNoise="noiseFixedAndBgAndCorrelated" | |
| vmessage(c("The noise model is:", modelNoise), verbose=verbose) | |
| vmessage(c("Proportion of non-genetic covariate variance (delta):", | |
| delta), verbose=verbose) | |
| vmessage(c("Proportion of variance of shared non-genetic covariate", | |
| "effects (gamma):", gamma), verbose=verbose) | |
| vmessage(c("Proportion of non-genetic covariates to have", | |
| "a trait-independent effect (pIndependentConfounders" | |
| , "):", pIndependentConfounders), | |
| verbose=verbose) | |
| vmessage(c("Proportion of traits influenced by independent", | |
| "non-genetic covariate effects", | |
| "(pTraitIndependentConfounders):", | |
| pTraitIndependentConfounders), verbose=verbose) | |
| vmessage(c("Proportion of variance of correlated noise effects", | |
| "(rho):", rho), verbose=verbose) | |
| vmessage(c("Proportion of observational noise variance (phi):", | |
| phi), verbose=verbose) | |
| vmessage(c("Variance of shared observational noise effect (alpha):", | |
| alpha), verbose=verbose) | |
| } | |
| vmessage("\n", verbose=verbose) | |
| } | |
| vmessage(c("The total genetic variance (genVar) is:", genVar), | |
| verbose=verbose) | |
| if ( genVar == 0 ) { | |
| modelGenetic="noGenetic" | |
| h2s <- 0 | |
| h2bg <- 0 | |
| vmessage(c("The genetic model is:", modelGenetic), verbose=verbose) | |
| } else { | |
| if (all(c(is.null(h2bg), is.null(h2s)))) { | |
| stop(paste("Neither genetic variant effect variance (h2s) nor", | |
| "infinitesimal genetic effect variance (h2bg) provided", | |
| ", at least one is required") | |
| ) | |
| } | |
| if (length(c(h2bg, h2s)) == 2 && sum(c(h2bg, h2s)) != 1) { | |
| stop(paste("Sum of the proportion of the variance of genetic", | |
| "effects is not equal to 1; change h2s (genetic variant", | |
| "effect variance) or h2bg (infinitesimal genetic effect", | |
| "variance) such that h2s + h2bg = 1") | |
| ) | |
| } | |
| if (is.null(h2s)) { | |
| h2s <- 1 - h2bg | |
| } else { | |
| h2bg <- 1 - h2s | |
| } | |
| if (!is.null(cNrSNP)) { | |
| if (h2s != 0 && cNrSNP == 0) { | |
| stop(paste("Proportion of variants of genetic variant effects ", | |
| "(h2s) is", h2s, "but number of", | |
| "cNrSNP is set to zero")) | |
| } | |
| } | |
| if (theta == 1) { | |
| pIndependentGenetic=0 | |
| pTraitIndependentGenetic=0 | |
| } | |
| if (theta == 0) { | |
| pIndependentGenetic=1 | |
| pTraitIndependentGenetic=1 | |
| } | |
| if ( h2s == 0) { | |
| modelGenetic="geneticBgOnly" | |
| vmessage(c("The genetic model is:", modelGenetic), verbose=verbose) | |
| vmessage(c("Proportion of variance of infinitesimal genetic", | |
| "effects (h2bg):", h2bg), verbose=verbose) | |
| vmessage(c("Proportion of variance of shared infinitesimal genetic", | |
| "effects (eta):", eta), verbose=verbose) | |
| } else if ( h2s == 1) { | |
| modelGenetic="geneticFixedOnly" | |
| vmessage(c("The genetic model is:", modelGenetic), verbose=verbose) | |
| vmessage(c("Proportion of variance of genetic variant effects", | |
| "(h2s):", h2s), verbose=verbose) | |
| vmessage(c("Proportion of variance of shared genetic variant", | |
| "effects (theta):", theta), verbose=verbose) | |
| vmessage(c("Proportion of genetic variant effects to have a", | |
| "trait-independent fixed effect", | |
| "(pIndependentGenetic):", pIndependentGenetic), | |
| verbose=verbose) | |
| vmessage(c("Proportion of traits influenced by independent", | |
| "genetic variant effects (pTraitIndependentGenetic):", | |
| pTraitIndependentGenetic), | |
| verbose=verbose) | |
| } else { | |
| modelGenetic="geneticFixedAndBg" | |
| vmessage(c("The genetic model is:", modelGenetic), verbose=verbose) | |
| vmessage(c("Proportion of variance of genetic variant effects", | |
| "(h2s):", h2s), verbose=verbose) | |
| vmessage(c("Proportion of variance of shared genetic variant", | |
| "effects (theta):", theta), verbose=verbose) | |
| vmessage(c("Proportion of genetic variant effects to have a", | |
| "trait-independent fixed effect", | |
| "(pIndependentGenetic):", pIndependentGenetic), | |
| verbose=verbose) | |
| vmessage(c("Proportion of traits influenced by independent", | |
| "genetic variant effects (pTraitIndependentGenetic):", | |
| pTraitIndependentGenetic), | |
| verbose=verbose) | |
| vmessage(c("Proportion of variance of infinitesimal genetic", | |
| "effects (h2bg):", h2bg), verbose=verbose) | |
| vmessage(c("Proportion of variance of shared infinitesimal genetic", | |
| "effects (eta):", eta), verbose=verbose) | |
| } | |
| vmessage(c("Proportion of non-linear phenotype transformation", | |
| "(proportionNonlinear):", proportionNonlinear), | |
| verbose=verbose) | |
| vmessage("\n", verbose=verbose) | |
| } | |
| return(list(modelGenetic=modelGenetic, modelNoise=modelNoise, | |
| genVar=genVar, h2s=h2s, h2bg=h2bg, noiseVar=noiseVar, rho=rho, | |
| delta=delta, phi=phi, gamma=gamma, theta=theta, eta=eta, | |
| alpha=alpha, pcorr=pcorr, | |
| pTraitIndependentGenetic=pTraitIndependentGenetic, | |
| pIndependentGenetic=pIndependentGenetic, | |
| pTraitIndependentConfounders=pTraitIndependentConfounders, | |
| pIndependentConfounders=pIndependentConfounders, | |
| proportionNonlinear=proportionNonlinear | |
| )) | |
| } | |
| #' Run phenotype simulation. | |
| #' | |
| #' runSimulation wraps around setModel, the phenotype component functions | |
| #' (genFixedEffects, genBgEffects, noiseBgEffects, noiseFixedEffects and | |
| #' correlatedBgEffects), rescales each component and combines them into the | |
| #' final phenotype. For details to all parameters, see the respective functions. | |
| #' | |
| #' @param P Number [integer] of phenotypes to simulate. | |
| #' @param N Number [integer] of samples to simulate. | |
| #' @param genVar Proportion [double] of total genetic variance. | |
| #' @param h2s Proportion [double] of genetic variance of genetic variant effects. | |
| #' @param h2bg Proportion [double] of genetic variance of infinitesimal genetic | |
| #' effects; either h2s or h2bg have to be specified and h2s + h2bg = 1. | |
| #' @param theta Proportion [double] of variance of shared genetic variant | |
| #' effects. | |
| #' @param eta Proportion [double] of variance of shared infinitesimal genetic | |
| #' effects. | |
| #' @param noiseVar Proportion [double] of total noise variance. | |
| #' @param rho Proportion [double] of noise variance of correlated effects; sum | |
| #' of rho, delta and phi has to be equal 1. | |
| #' @param delta Proportion [double] of noise variance of non-genetic covariate | |
| #' effects; sum of rho, delta and phi has to be equal 1. | |
| #' @param gamma Proportion [double] of variance of shared non-genetic covariate | |
| #' effects. | |
| #' @param phi Proportion [double] of noise variance of observational noise | |
| #' effects; sum of rho, delta and phi has to be equal 1. | |
| #' @param alpha Variance [double] of shared observational noise effect. | |
| #' @param tNrSNP Total number [integer] of SNPs to simulate; these SNPs are used | |
| #' for kinship estimation. | |
| #' @param cNrSNP Number [integer] of causal SNPs; used as genetic variant | |
| #' effects. | |
| #' @param SNPfrequencies Vector of allele frequencies [double] from which to | |
| #' sample. | |
| #' @param genotypefile Needed when reading external genotypes (into memory), | |
| #' path/to/genotype file [string] in format specified by \link{format}. | |
| #' @param format Needed when reading external genotypes, specifies | |
| #' the format of the genotype data; has to be one of plink, oxgen, genome, | |
| #' bimbam and delim when reading files into memory, or one of oxgen, bimbam or | |
| #' delim if sampling genetic variants from file; for details see | |
| #' \link{readStandardGenotypes} and \link{getCausalSNPs}. | |
| #' @param genoFilePrefix Needed when sampling cuasal SNPs from file, full | |
| #' path/to/chromosome-wise-genotype-file-ending-before-"chrChromosomeNumber" | |
| #' (no '~' expansion!) [string] | |
| #' @param genoFileSuffix Needed when sampling causal SNPs from file, | |
| #' following chromosome number including fileformat (e.g. ".csv") [string] | |
| #' @param genoDelimiter Field separator [string] of genotypefile or genoFile if | |
| #' format == delim. | |
| #' @param skipFields Number [integer] of fields (columns) in to skip in | |
| #' genoFilePrefix-genoFileSuffix-file. See details in \link{getCausalSNPs} if | |
| #' format == delim. | |
| #' @param header [logical] Can be set to indicate if | |
| #' genoFilePrefix-genoFileSuffix file has a header for format == 'delim'. | |
| #' See details in \link{getCausalSNPs}. | |
| #' @param probabilities [bool]. If set to TRUE, the genotypes in the files | |
| #' described by genoFilePrefix and genoFileSuffix are provided as triplets of | |
| #' probablities (p(AA), p(Aa), p(aa)) and are converted into their expected | |
| #' genotype frequencies by 0*p(AA) + p(Aa) + 2p(aa) via \link{probGen2expGen}. | |
| #' @param chr Numeric vector of chromosomes [integer] to chose NrCausalSNPs | |
| #' from; only used when external genotype data is sampled i.e. | |
| #' !is.null(genoFilePrefix) | |
| #' @param NrSNPsOnChromosome Specifies the number of SNPs [integer] per entry in | |
| #' chr (see above); has to be the same length as chr. If not provided, lines in | |
| #' genoFilePrefix-genoFileSuffix file will be counted (which can be slow for | |
| #' large files). | |
| #' @param NrChrCausal Number [integer] of causal chromosomes to chose | |
| #' NrCausalSNPs from (as opposed to the actual chromosomes to chose from via chr | |
| #' ); only used when external genotype data is sampled i.e. | |
| #' !is.null(genoFilePrefix). | |
| #' @param kinshipfile path/to/kinshipfile [string]; if provided, | |
| #' kinship for simulation of genetic backgound effect will be read from file. | |
| #' @param kinshipHeader [boolean] If TRUE kinship file has header information. | |
| #' @param kinshipDelimiter Field separator [string] of kinship file. | |
| #' @param standardise [boolean] If TRUE genotypes will be standardised for | |
| #' kinship estimation (recommended). | |
| #' @param distBetaGenetic Name [string] of distribution to use to simulate | |
| #' effect sizes of genetic variants; one of "unif" or "norm". | |
| #' @param mBetaGenetic Mean/midpoint [double] of normal/uniform distribution | |
| #' for effect sizes of genetic variants. | |
| #' @param sdBetaGenetic Standard deviation/extension from midpoint [double] | |
| #' of normal/uniform distribution for effect sizes of genetic variants. | |
| #' @param pIndependentGenetic Proportion [double] of genetic variant effects to | |
| #' have a trait-independent fixed effect. | |
| #' @param pTraitIndependentGenetic Proportion [double] of traits influenced by | |
| #' independent genetic variant effects. | |
| #' @param keepSameIndependentSNPs [boolean] If set to TRUE, the | |
| #' independent SNPs effects always influence the same subset of traits. | |
| #' @param pTraitsAffectedGenetics Proportion [double] of traits affected by the | |
| #' genetic variant effect. For non-integer results of pTraitsAffected*P, the | |
| #' ceiling of the result is used. Allows to simulate for instance different | |
| #' levels of pleiotropy. | |
| #' @param NrFixedEffects Number [integer] of different non-genetic covariate | |
| #' effects to simulate; allows to simulate non-genetic covariate effects from | |
| #' different distributions or with different parameters. | |
| #' @param NrConfounders Number [integer] of non-genetic covariates; used as | |
| #' non-genetic covariate effects. | |
| #' @param distConfounders Vector of name(s) [string] of distributions to use to | |
| #' simulate confounders; one of "unif", "norm", "bin", "cat_norm", "cat_unif". | |
| #' @param mConfounders Vector of mean(s)/midpoint(s) [double] of | |
| #' normal/uniform distribution for confounders. | |
| #' @param sdConfounders Vector of standard deviation(s)/extension from | |
| #' midpoint(s) [double] of normal/uniform distribution for confounders. | |
| #' @param catConfounders Vector of confounder categories [factor]; required if | |
| #' distConfounders "cat_norm" or "cat_unif". | |
| #' @param probConfounders Vector of probability(ies) [double] of binomial | |
| #' confounders (0/1); required if distConfounders "bin". | |
| #' @param distBetaConfounders Vector of name(s) [string] of distribution to use | |
| #' to simulate effect sizes of confounders; one of "unif" or "norm". | |
| #' @param mBetaConfounders Vector of mean(s)/midpoint(s) [double] of | |
| #' normal/uniform distribution for effect sizes of confounders. | |
| #' @param sdBetaConfounders Vector of standard deviation(s)/extension from | |
| #' midpoint(s) [double] of normal/uniform distribution for effect sizes of | |
| #' confounders. | |
| #' @param pIndependentConfounders Vector of proportion(s) [double] of | |
| #' non-genetic covariate effects to have a trait-independent effect. | |
| #' @param pTraitIndependentConfounders Vector of proportion(s) [double] of | |
| #' traits influenced by independent non-genetic covariate effects. | |
| #' @param keepSameIndependentConfounders [boolean] If set to TRUE, the | |
| #' independent confounder effects always influence the same subset of traits. | |
| #' @param pTraitsAffectedConfounders Proportion(s) [double] of traits | |
| #' affected by the non-genetic covariates. For non-integer results of | |
| #' pTraitsAffected*P, the ceiling of the result is used. | |
| #' @param meanNoiseBg Mean [double] of the normal distributions for the | |
| #' simulation observational noise effects. | |
| #' @param sdNoiseBg Standard deviation [double] of the normal distributions for | |
| #' the simulations of the observational noise effects. | |
| #' @param pcorr Correlation [double] between phenotypes. | |
| #' @param corrmatfile path/to/corrmatfile.csv [string] with comma-separated | |
| #' [P x P] numeric [double] correlation matrix; if provided, correlation matrix | |
| #' for simulation of correlated backgound effect will be read from file; | |
| #' file should NOT contain an index or header column. | |
| #' @param nonlinear nonlinear transformation method [string]; one exp | |
| #' (exponential), log (logarithm), poly (polynomial), sqrt (squareroot) or | |
| #' custom (user-supplied function); if log or exp, base can be specified; if | |
| #' poly, power can be specified; if custom, a custom function (see for details). | |
| #' Non-linear transformation is optional, default is NULL ie no transformation | |
| #' (see details). | |
| #' @param logbase [int] base of logarithm for non-linear phenotype | |
| #' transformation (see details). | |
| #' @param expbase [int] base of exponential function for non-linear phenotype | |
| #' transformation (see details). | |
| #' @param power [double] power of polynomial function for non-linear phenotype | |
| #' transformation. | |
| #' @param transformNeg [string] transformation method for negative values in non | |
| #' linear phenotype transformation. One of abs (absolute value) or set0 (set all | |
| #' negative values to zero). If nonlinear==log and transformNeg==set0, negative | |
| #' values set to 1e-5 | |
| #' @param customTransform [function] custom transformation function accepting | |
| #' a single argument. | |
| #' @param proportionNonlinear [double] proportion of the phenotype to be non- | |
| #' linear (see details) | |
| #' @param sampleID Prefix [string] for naming samples (will be followed by | |
| #' sample number from 1 to N when constructing sample IDs); only used if | |
| #' genotypes/kinship are simulated/do not have sample IDs. | |
| #' @param phenoID Prefix [string] for naming traits (will be followed by | |
| #' phenotypes number from 1 to P when constructing phenotype IDs). | |
| #' @param snpID Prefix [string] for naming SNPs (will be followed by | |
| #' SNP number from 1 to NrSNP when constructing SNP IDs). | |
| #' @param seed Seed [integer] to initiate random number generation. | |
| #' @param verbose [boolean]; If TRUE, progress info is printed to standard out | |
| #' @return Named list of i) dataframe of proportion of variance | |
| #' explained for each component (varComponents), | |
| #' ii) a named list with the final simulated phenotype components | |
| #' (phenoComponentsFinal), iii) a named list with the intermediate simulated | |
| #' phenotype components (phenoComponentsIntermediate), iv) a named list of | |
| #' parameters describing the model setup (setup) and v) a named list of raw | |
| #' components (rawComponents) used for genetic effect simulation (genotypes | |
| #' and/or kinship, eigenvalues and eigenvectors of kinship) | |
| #' @details Phenotypes are modeled under a linear additive model where | |
| #' Y = WA + BX + G + C + Phi, with WA the non-genetic covariates, BX the genetic | |
| #' variant effects, G the infinitesimal genetic effects, C the correlated | |
| #' background effects and the Phi the observational noise. For more information | |
| #' on these components look at the respective function descriptions (see also) | |
| #' Optionally the phenotypes can be non-linearly transformed via: | |
| #' Y_trans = (1-alpha) x Y + alpha x f(Y). Alpha is the proportion of non- | |
| #' linearity of the phenotype and f is a non-linear transformation, and one of | |
| #' exp, log or sqrt. | |
| #' @export | |
| #' @seealso \link{setModel}, \link{geneticFixedEffects}, | |
| #' \link{geneticBgEffects}, \link{noiseBgEffects}, \link{noiseFixedEffects}, | |
| #' \link{correlatedBgEffects} and \link{rescaleVariance}. | |
| #' @examples | |
| #' # simulate phenotype of 100 samples, 10 traits from genetic and noise | |
| #' # background effects, with variance explained of 0.2 and 0.8 respectively | |
| #' genVar = 0.2 | |
| #' simulatedPhenotype <- runSimulation(N=100, P=5, cNrSNP=10, | |
| #' genVar=genVar, h2s=1, phi=1) | |
| runSimulation <- function(N, P, | |
| genVar=NULL, h2s=NULL, theta=0.8, h2bg=NULL, eta=0.8, | |
| noiseVar=NULL, rho=NULL, delta=NULL, gamma=0.8, | |
| phi=NULL, alpha=0.8, | |
| tNrSNP=5000, cNrSNP=20, | |
| SNPfrequencies=c(0.1, 0.2, 0.4), | |
| genotypefile=NULL, format='delim', | |
| genoFilePrefix=NULL, genoFileSuffix=NULL, | |
| genoDelimiter=",", skipFields=NULL, | |
| header=FALSE, | |
| probabilities=FALSE, | |
| chr=NULL, NrSNPsOnChromosome=NULL, | |
| NrChrCausal=NULL, | |
| kinshipfile=NULL, | |
| kinshipHeader=FALSE, kinshipDelimiter=",", | |
| standardise=TRUE, | |
| distBetaGenetic="norm", mBetaGenetic=0, | |
| sdBetaGenetic=1,pTraitsAffectedGenetics=1, | |
| pIndependentGenetic=0.4, pTraitIndependentGenetic=0.2, | |
| keepSameIndependentSNPs=FALSE, | |
| NrFixedEffects=1, NrConfounders=10, | |
| distConfounders="norm", mConfounders=0, | |
| sdConfounders=1,catConfounders=NULL, | |
| probConfounders=NULL, distBetaConfounders="norm", | |
| mBetaConfounders=0, sdBetaConfounders=1, | |
| pTraitsAffectedConfounders=1, | |
| pIndependentConfounders=0.4, | |
| pTraitIndependentConfounders=0.2, | |
| keepSameIndependentConfounders=FALSE, | |
| pcorr=0.8, corrmatfile=NULL, | |
| meanNoiseBg=0, sdNoiseBg=1, | |
| nonlinear=NULL, logbase=10, expbase=NULL, power=NULL, | |
| customTransform=NULL, transformNeg="abs", | |
| proportionNonlinear=0, | |
| sampleID="ID_", phenoID="Trait_", snpID="SNP_", | |
| seed=219453, verbose=FALSE) { | |
| vmessage(c("Set seed:", seed), verbose=verbose) | |
| set.seed(seed) | |
| model <- setModel(genVar=genVar, h2s=h2s, h2bg=h2bg, theta=theta, eta=eta, | |
| noiseVar=noiseVar, rho=rho, delta=delta, gamma=gamma, | |
| phi=phi, alpha=alpha, pcorr=pcorr, | |
| pIndependentConfounders=pIndependentConfounders, | |
| pTraitIndependentConfounders=pTraitIndependentConfounders, | |
| pIndependentGenetic=pIndependentGenetic, | |
| pTraitIndependentGenetic=pTraitIndependentGenetic, | |
| proportionNonlinear=proportionNonlinear, | |
| cNrSNP=cNrSNP, NrConfounders=NrConfounders, | |
| verbose=verbose) | |
| id_snps <- NULL | |
| id_samples <- NULL | |
| id_phenos <- paste(phenoID, 1:P, sep="") | |
| ### create simulated phenotypes | |
| # 1. Simulate genetic terms | |
| vmessage(c("Simulate genetic effects (genetic model: ", | |
| model$modelGenetic,")"), sep="", verbose=verbose) | |
| if (grepl('Fixed', model$modelGenetic)) { | |
| if (is.null(genoFilePrefix) && is.null(genotypefile)) { | |
| if (!grepl('Bg', model$modelGenetic) && tNrSNP != cNrSNP) { | |
| warning(paste("The genetic model does not contain infinitesimal", | |
| "genetic effects but the total number of SNPs to", | |
| "simulate (tNrSNP:", | |
| tNrSNP, ") is larger than the number of genetic | |
| variant effects SNPs (cNrSNP:", cNrSNP, | |
| "). If genotypes are not needed, consider setting | |
| tNrSNPs=cNrSNPs to speed up computation")) | |
| } | |
| genotypes <- simulateGenotypes(N=N, NrSNP=tNrSNP, | |
| frequencies=SNPfrequencies, | |
| sampleID=sampleID, | |
| snpID=snpID, | |
| verbose=verbose) | |
| id_samples <- genotypes$id_samples | |
| id_snps <- genotypes$id_snps | |
| } else if (! is.null(genotypefile)) { | |
| genotypes <- readStandardGenotypes(N=N, filename=genotypefile, | |
| format=format, | |
| verbose=verbose, | |
| sampleID=sampleID, | |
| snpID=snpID, | |
| delimiter=genoDelimiter) | |
| id_samples <- genotypes$id_samples | |
| id_snps <- genotypes$id_snps | |
| } else { | |
| genotypes <- NULL | |
| } | |
| causalSNPs <- getCausalSNPs(N=N, NrCausalSNPs=cNrSNP, chr=chr, | |
| NrChrCausal=NrChrCausal, | |
| NrSNPsOnChromosome=NrSNPsOnChromosome, | |
| genotypes=genotypes$genotypes, | |
| genoFilePrefix=genoFilePrefix, | |
| genoFileSuffix=genoFileSuffix, | |
| format=format, | |
| probabilities=probabilities, | |
| skipFields=skipFields, | |
| header=header, | |
| delimiter=genoDelimiter, | |
| sampleID=sampleID, | |
| verbose=verbose) | |
| if (is.null(id_snps)) { | |
| id_snps <- colnames(causalSNPs) | |
| } | |
| if (is.null(id_samples)) { | |
| id_samples <- rownames(causalSNPs) | |
| } | |
| vmessage("Simulate genetic variant effects", verbose=verbose) | |
| genFixed <- geneticFixedEffects(X_causal=causalSNPs, N=N, P=P, | |
| pTraitsAffected= | |
| pTraitsAffectedGenetics, | |
| pIndependentGenetic= | |
| model$pIndependentGenetic, | |
| pTraitIndependentGenetic= | |
| model$pTraitIndependentGenetic, | |
| keepSameIndependent= | |
| keepSameIndependentSNPs, | |
| distBeta=distBetaGenetic, | |
| mBeta=mBetaGenetic, | |
| sdBeta=sdBetaGenetic, | |
| id_phenos=id_phenos, | |
| id_samples=id_samples, | |
| phenoID=phenoID, | |
| verbose=verbose) | |
| var_genFixed_shared <- model$theta * model$h2s * model$genVar | |
| var_genFixed_independent <- (1 - model$theta) * model$h2s * model$genVar | |
| genFixed_shared_rescaled <- rescaleVariance(genFixed$shared, | |
| var_genFixed_shared) | |
| genFixed_independent_rescaled <- | |
| rescaleVariance(genFixed$independent, var_genFixed_independent) | |
| Y_genFixed <- addNonNulls(list(genFixed_shared_rescaled$component, | |
| genFixed_independent_rescaled$component)) | |
| } else { | |
| genFixed <- NULL | |
| genFixed_shared_rescaled <- NULL | |
| genFixed_independent_rescaled <- NULL | |
| genotypes <- NULL | |
| var_genFixed_shared <- 0 | |
| var_genFixed_independent <- 0 | |
| Y_genFixed <- NULL | |
| cNrSNP <- 0 | |
| } | |
| if (grepl('Bg', model$modelGenetic)) { | |
| if (is.null(kinshipfile)) { | |
| if (is.null(genotypes) && !is.null(genotypefile)){ | |
| genotypes <- readStandardGenotypes(genotypefile, format=format, | |
| verbose=verbose, | |
| sampleID = sampleID, | |
| snpID = snpID, | |
| delimiter = genoDelimiter) | |
| } | |
| if (is.null(genotypes)) { | |
| genotypes <- simulateGenotypes(N=N, NrSNP=tNrSNP, | |
| frequencies=SNPfrequencies, | |
| sampleID=sampleID, | |
| verbose=verbose) | |
| } | |
| id_samples <- genotypes$id_samples | |
| id_snps <- genotypes$id_snps | |
| kinship <- getKinship(X=genotypes$genotypes, | |
| N=N, standardise=standardise, | |
| id_samples=id_samples, | |
| sampleID=sampleID, | |
| verbose=verbose) | |
| } else { | |
| vmessage("Read kinship from file", verbose=verbose) | |
| kinship <- getKinship(kinshipfile=kinshipfile, | |
| sep=kinshipDelimiter, | |
| N=N, header=kinshipHeader, verbose=verbose, | |
| sampleID=sampleID, id_samples=id_samples) | |
| if(!is.null(genotypes)) { | |
| if(colnames(kinship) != genotypes$id_samples) { | |
| stop(paste("Sample names in kinship file do not match", | |
| "sample names of genotypes", sep="")) | |
| } | |
| } else { | |
| id_samples <- colnames(kinship) | |
| } | |
| } | |
| if (eta == 1) { | |
| genBgShared <- TRUE | |
| genBgIndependent <- FALSE | |
| } | |
| else if (eta == 0) { | |
| genBgShared <- FALSE | |
| genBgIndependent <- TRUE | |
| } else { | |
| genBgShared <- TRUE | |
| genBgIndependent <- TRUE | |
| } | |
| vmessage("Simulate infinitesimal genetic effects", verbose=verbose) | |
| genBg <- geneticBgEffects(N=N, P=P, kinship=kinship, | |
| shared=genBgShared, | |
| independent=genBgIndependent, | |
| id_phenos=id_phenos) | |
| eval_kinship <- genBg$eval_kinship | |
| evec_kinship <- genBg$evec_kinship | |
| var_genBg_shared <- model$eta * model$h2bg * model$genVar | |
| var_genBg_independent <- (1 - model$eta) * model$h2bg * model$genVar | |
| genBg_shared_rescaled <- rescaleVariance(genBg$shared, var_genBg_shared) | |
| genBg_independent_rescaled <- rescaleVariance(genBg$independent, | |
| var_genBg_independent) | |
| Y_genBg <- addNonNulls(list(genBg_shared_rescaled$component, | |
| genBg_independent_rescaled$component)) | |
| cov_genBg_shared <- genBg$cov_shared | |
| cov_genBg_shared_rescaled <- cov_genBg_shared * | |
| genBg_shared_rescaled$scale_factor^2 | |
| cov_genBg_independent <- genBg$cov_independent | |
| cov_genBg_independent_rescaled <- cov_genBg_independent * | |
| genBg_independent_rescaled$scale_factor^2 | |
| cov_genBg <- cov_genBg_shared_rescaled + cov_genBg_independent_rescaled | |
| } else { | |
| genBg <- NULL | |
| genBg_shared_rescaled <- NULL | |
| genBg_independent_rescaled <- NULL | |
| kinship <- NULL | |
| var_genBg_shared <- 0 | |
| var_genBg_independent <- 0 | |
| Y_genBg <- NULL | |
| cov_genBg <- NULL | |
| cov_genBg_shared <- NULL | |
| cov_genBg_independent <- NULL | |
| eval_kinship <- NULL | |
| evec_kinship <- NULL | |
| } | |
| # 1. Simulate noise terms | |
| vmessage(c("Simulate noise terms (noise model: ", model$modelNoise, ")"), | |
| sep="", verbose=verbose) | |
| if (grepl('Correlated', model$modelNoise)) { | |
| corr_mat <- NULL | |
| vmessage("Simulate correlated background effects", verbose=verbose) | |
| if (!is.null(corrmatfile)){ | |
| vmessage("Read file with correlation matrix for correlated", | |
| "background effect", verbose=verbose) | |
| corr_mat <- as.matrix(data.table::fread(corrmatfile, sep=",", | |
| data.table=FALSE, | |
| stringsAsFactors=FALSE)) | |
| } | |
| correlatedBg <- correlatedBgEffects(N=N, P=P, corr_mat=corr_mat, | |
| pcorr=model$pcorr, | |
| id_phenos=id_phenos, | |
| id_samples=id_samples, | |
| sampleID=sampleID, | |
| phenoID=phenoID) | |
| var_noiseCorrelated <- model$rho * model$noiseVar | |
| correlatedBg_rescaled <- rescaleVariance(correlatedBg$correlatedBg, | |
| var_noiseCorrelated) | |
| Y_correlatedBg <- correlatedBg_rescaled$component | |
| cov_correlatedBg <- correlatedBg$cov_correlated * | |
| correlatedBg_rescaled$scale_factor^2 | |
| } else { | |
| correlatedBg <- NULL | |
| var_noiseCorrelated <- 0 | |
| Y_correlatedBg <- NULL | |
| cov_correlatedBg <- NULL | |
| } | |
| if (grepl('Bg', model$modelNoise)) { | |
| if (alpha == 1) { | |
| noiseBgShared <- TRUE | |
| noiseBgIndependent <- FALSE | |
| } | |
| else if (alpha == 0) { | |
| noiseBgShared <- FALSE | |
| noiseBgIndependent <- TRUE | |
| } else { | |
| noiseBgShared <- TRUE | |
| noiseBgIndependent <- TRUE | |
| } | |
| vmessage("Simulate observational noise effects", verbose=verbose) | |
| noiseBg <- noiseBgEffects(N=N, P=P, mean=meanNoiseBg, sd=sdNoiseBg, | |
| shared=noiseBgShared, | |
| independent=noiseBgIndependent, | |
| id_phenos=id_phenos, id_samples=id_samples, | |
| sampleID=sampleID, phenoID=phenoID) | |
| var_noiseBg_shared <- model$alpha * model$phi * model$noiseVar | |
| var_noiseBg_independent <- (1 - model$alpha) * model$phi * model$noiseVar | |
| noiseBg_shared_rescaled <- rescaleVariance(noiseBg$shared, | |
| var_noiseBg_shared) | |
| noiseBg_independent_rescaled <- rescaleVariance(noiseBg$independent, | |
| var_noiseBg_independent) | |
| Y_noiseBg <- addNonNulls(list(noiseBg_shared_rescaled$component, | |
| noiseBg_independent_rescaled$component)) | |
| cov_noiseBg_shared <- noiseBg$cov_shared | |
| cov_noiseBg_shared_rescaled <- cov_noiseBg_shared * | |
| noiseBg_shared_rescaled$scale_factor^2 | |
| cov_noiseBg_independent <- noiseBg$cov_independent | |
| cov_noiseBg_independent_rescaled <- cov_noiseBg_independent * | |
| noiseBg_independent_rescaled$scale_factor^2 | |
| cov_noiseBg <- cov_noiseBg_shared_rescaled + | |
| cov_noiseBg_independent_rescaled | |
| } else { | |
| noiseBg <- NULL | |
| noiseBg_shared_rescaled <- NULL | |
| noiseBg_independent_rescaled <- NULL | |
| var_noiseBg_shared <- 0 | |
| var_noiseBg_independent <- 0 | |
| Y_noiseBg <- NULL | |
| cov_noiseBg <- NULL | |
| cov_noiseBg_shared <- NULL | |
| cov_noiseBg_independent <- NULL | |
| } | |
| if (grepl('Fixed', model$modelNoise)) { | |
| vmessage("Simulate confounder effects", verbose=verbose) | |
| noiseFixed <- noiseFixedEffects(P=P, N=N, | |
| NrFixedEffects = NrFixedEffects, | |
| NrConfounders=NrConfounders, | |
| pTraitsAffected= | |
| pTraitsAffectedConfounders, | |
| pIndependentConfounders= | |
| model$pIndependentConfounders, | |
| pTraitIndependentConfounders= | |
| model$pTraitIndependentConfounders, | |
| keepSameIndependent= | |
| keepSameIndependentConfounders, | |
| distConfounders=distConfounders, | |
| mConfounders=mConfounders, | |
| sdConfounders=sdConfounders, | |
| catConfounders=catConfounders, | |
| probConfounders = probConfounders, | |
| distBeta=distBetaConfounders, | |
| mBeta=mBetaConfounders, | |
| sdBeta=sdBetaConfounders, | |
| id_phenos=id_phenos, | |
| id_samples=id_samples, | |
| sampleID=sampleID, | |
| phenoID=phenoID, | |
| verbose=verbose) | |
| var_noiseFixed_shared <- model$gamma * model$delta * model$noiseVar | |
| var_noiseFixed_independent <- (1 - model$gamma) * model$delta * | |
| model$noiseVar | |
| noiseFixed_shared_rescaled <- rescaleVariance(noiseFixed$shared, | |
| var_noiseFixed_shared) | |
| noiseFixed_independent_rescaled <- rescaleVariance( | |
| noiseFixed$independent, | |
| var_noiseFixed_independent) | |
| Y_noiseFixed <- | |
| addNonNulls(list(noiseFixed_shared_rescaled$component, | |
| noiseFixed_independent_rescaled$component)) | |
| } else { | |
| noiseFixed <- NULL | |
| noiseFixed_shared_rescaled <- NULL | |
| noiseFixed_independent_rescaled <- NULL | |
| var_noiseFixed_shared <- 0 | |
| var_noiseFixed_independent <- 0 | |
| Y_noiseFixed <- NULL | |
| } | |
| # 3. Construct final simulated phenotype | |
| vmessage("Construct final simulated phenotype" | |
| , verbose=verbose) | |
| components <- list(Y_genFixed, Y_genBg, Y_noiseFixed, Y_noiseBg, | |
| Y_correlatedBg) | |
| Y <- addNonNulls(components) | |
| # 4. Transformation | |
| if (!is.null(nonlinear)) { | |
| vmessage("Transform phenotypes", verbose=verbose) | |
| Y_transformed <- transformNonlinear(Y, method=nonlinear, | |
| alpha=proportionNonlinear, | |
| logbase=logbase, expbase=expbase, | |
| power=power, f=customTransform, | |
| transformNeg=transformNeg) | |
| Y_transformed <- scale(Y_transformed) | |
| } else { | |
| Y_transformed <- NULL | |
| } | |
| Y <- scale(Y) | |
| varComponents <- data.frame(genVar=model$genVar, h2s=model$h2s, | |
| h2bg=model$h2bg, | |
| proportionNonlinear=model$proportionNonlinear, | |
| var_genFixed_shared=var_genFixed_shared, | |
| var_genFixed_independent= | |
| var_genFixed_independent, | |
| var_genBg_shared=var_genBg_shared, | |
| var_genBg_independent=var_genBg_independent, | |
| noiseVar=model$noiseVar, | |
| var_noiseFixed_shared=var_noiseFixed_shared, | |
| var_noiseFixed_independent= | |
| var_noiseFixed_independent, | |
| var_noiseBg_shared=var_noiseBg_shared, | |
| var_noiseBg_independent=var_noiseBg_independent, | |
| var_noiseCorrelated=var_noiseCorrelated) | |
| phenoComponentsFinal <- list(Y=Y, Y_genFixed=Y_genFixed, Y_genBg=Y_genBg, | |
| Y_noiseFixed=Y_noiseFixed, Y_noiseBg=Y_noiseBg, | |
| Y_correlatedBg=Y_correlatedBg, | |
| Y_transformed=Y_transformed, | |
| cov_genBg=cov_genBg, | |
| cov_noiseBg=cov_noiseBg, | |
| cov_correlatedBg = cov_correlatedBg) | |
| phenoComponentsIntermediate <- list( | |
| Y_genFixed_shared= | |
| genFixed_shared_rescaled$component, | |
| Y_genFixed_independent= | |
| genFixed_independent_rescaled$component, | |
| Y_noiseFixed_shared= | |
| noiseFixed_shared_rescaled$component, | |
| Y_noiseFixed_independent= | |
| noiseFixed_independent_rescaled$component, | |
| Y_genBg_shared= | |
| genBg_shared_rescaled$component, | |
| Y_genBg_independent= | |
| genBg_independent_rescaled$component, | |
| Y_noiseBg_shared= | |
| noiseBg_shared_rescaled$component, | |
| Y_noiseBg_independent= | |
| noiseBg_independent_rescaled$component, | |
| cov_genBg_shared=cov_genBg_shared, | |
| cov_genBg_independent=cov_genBg_independent, | |
| cov_noiseBg_shared=cov_noiseBg_shared, | |
| cov_noiseBg_independent=cov_noiseBg_independent, | |
| genFixed=genFixed, | |
| noiseFixed=noiseFixed) | |
| setup <- list(P=P, N=N, NrCausalSNPs=cNrSNP, | |
| modelGenetic=model$modelGenetic, | |
| modelNoise=model$modelNoise, | |
| id_samples=id_samples, id_phenos=id_phenos, id_snps=id_snps) | |
| rawComponents <- list(kinship=kinship, genotypes=genotypes, | |
| eval_kinship=eval_kinship, evec_kinship=evec_kinship) | |
| return(list(varComponents=varComponents, | |
| phenoComponentsFinal=phenoComponentsFinal, | |
| phenoComponentsIntermediate=phenoComponentsIntermediate, | |
| setup=setup, rawComponents=rawComponents)) | |
| } | |