| Title: | Causal Inference Test |
|---|---|
| Description: | A likelihood-based hypothesis testing approach is implemented for assessing causal mediation. Described in Millstein, Chen, and Breton (2016), <doi:10.1093/bioinformatics/btw135>, it could be used to test for mediation of a known causal association between a DNA variant, the 'instrumental variable', and a clinical outcome or phenotype by gene expression or DNA methylation, the potential mediator. Another example would be testing mediation of the effect of a drug on a clinical outcome by the molecular target. The hypothesis test generates a p-value or permutation-based FDR value with confidence intervals to quantify uncertainty in the causal inference. The outcome can be continuous, binary, or a censored survival time; the potential mediator is continuous and may be multivariate; and the instrumental variable can be continuous or binary and is not limited to a single variable but may be a design matrix representing multiple variables. |
| Authors: | Joshua Millstein [aut, cre] (ORCID: <https://orcid.org/0000-0001-7961-8943>), Mingzhi Ye [aut] |
| Maintainer: | Joshua Millstein <[email protected]> |
| License: | Artistic-2.0 |
| Version: | 3.0.0 |
| Built: | 2026-11-01 18:52:58 UTC |
| Source: | https://github.com/uscbiostats/cit |
This package implements a formal statistical hypothesis test for causal mediation. For example, it could be used to test for mediation of a known causal association between a DNA variant, the 'instrumental variable', and a clinical outcome or phenotype by gene expression or DNA methylation, the potential mediator. Another example would be testing mediation of the effect of a drug on a clinical outcome by the molecular target. The hypothesis test generates a p-value or permutation-based false discovery rate (FDR) value with confidence intervals to quantify uncertainty in the causal inference. The outcome can be continuous (cit.cp), binary (cit.bp) or a censored survival time (cit.cox); the potential mediator is continuous and may be multivariate; and the instrumental variable can be continuous or binary and is not limited to a single variable but may be a design matrix representing multiple variables.
This package implements a statistical framework in which existing notions of causal mediation are formalized into a hypothesis test, thus providing a standard quantitative measure of uncertainty in the form of a p-value and optionally, permutation-based FDR. We treat the causal inference as a 'chain' of mathematical conditions that must be satisfied to conclude that the potential mediator is causal for the trait, where the inference is only as good as the weakest link in the chain. P-values are computed for the component conditions. The intersection-union test, in which a series of statistical tests are combined to form an omnibus test, is then employed to generate the overall test result as the maximum of the component p-values. If we let L denote a locus of interest or other instrumental variable, T denote an outcome trait, and G denote a potential mediator, then the four component conditions are 1) L and T are associated, 2) G and T are associated given L, 3) L and G are associated given T, and 4) the association between L and T given G is weaker than expected if G were not a mediator. Condition 4 is tested against null data generated without mediation (Millstein et al. 2009), so the test addresses mediation and can detect it when L also has a direct effect on T.
Version 3.0.0 adds survival outcomes, multiple mediators, separate covariates for the outcome and mediator models, CIT v2 (robust = TRUE), which prevents a pair of variables from being called causal in both directions and is recommended when reverse causation is a concern, and a choice between empirical and parametric versions of test 4.
Joshua Millstein, Mingzhi Ye
Maintainer: Joshua Millstein <[email protected]>
Millstein J, Chen GK, Breton CV. 2016. cit: hypothesis testing software for mediation analysis in genomic applications. Bioinformatics. btw135. PMID: 27153715.
Millstein J, Zhang B, Zhu J, Schadt EE. 2009. Disentangling molecular relationships with a causal inference test. BMC Genetics, 10:23.
This function implements a formal statistical hypothesis test, resulting in a p-value, to quantify uncertainty in a causal inference pertaining to a measured factor, e.g. a molecular species, which potentially mediates a known causal association between a locus or other instrumental variable and a binary trait. If the number of permutations is greater than zero, then the results can be used with fdr.cit to generate permutation-based FDR values (q-values) that are returned with confidence intervals to quantify uncertainty in the estimate. The outcome is binary, the potential mediator is continuous and may be multivariate, and the instrumental variable can be continuous, discrete (such as coding a SNP 0, 1, 2), or binary and is not limited to a single variable but may be a design matrix representing multiple variables.
cit.bp(L, G, T, CT = NULL, maxit = 10000, n.perm = 0, perm.index = NULL, rseed = NULL, CG = NULL, robust = FALSE, test4 = c("empirical", "parametric"), n.resampl = 50, sign.pattern = TRUE, C = NULL)cit.bp(L, G, T, CT = NULL, maxit = 10000, n.perm = 0, perm.index = NULL, rseed = NULL, CG = NULL, robust = FALSE, test4 = c("empirical", "parametric"), n.resampl = 50, sign.pattern = TRUE, C = NULL)
L |
Vector or n x p design matrix representing the instrumental variable(s). |
G |
Continuous vector, or n x k matrix of k continuous potential mediators tested jointly. |
T |
Binary vector, coded 0/1, representing the clinical trait or outcome of interest. |
CT |
Vector or n x p matrix of adjustment covariates for the models with |
n.resampl |
Number of null instances of the test 4 statistic used to estimate the null distribution when |
n.perm |
If greater than 0, |
perm.index |
Optional n x n.perm matrix whose columns are permutations of |
rseed |
Optional integer seed. When supplied, results are reproducible and the caller's random number stream is left unchanged. When |
CG |
Vector or n x p matrix of adjustment covariates for the models with |
robust |
If |
test4 |
How the test 4 p-value is obtained from null data: |
maxit |
Maximum number of null statistics for the empirical test 4. |
sign.pattern |
If |
C |
Deprecated name for |
Let L denote the instrumental variable, G the potential mediator and T the outcome. The CIT is an intersection-union test of four component conditions, each with its own p-value:
T is associated with L (p_TassocL; likelihood-ratio test of L in the logistic regression T ~ CT + L).
T is associated with G given L (p_TassocGgvnL; likelihood-ratio test of G in T ~ CT + L + G).
G is associated with L given T (p_GassocLgvnT; F test of L in G ~ T + CG + L, or Wilks' lambda when G has several columns).
The association between L and T given G is weaker than expected if G were not a mediator (p_LindTgvnG).
The omnibus p-value, p_cit, is the maximum of the component p-values, representing the probability of the data if at least one of the component null hypotheses is true.
Test 4 compares the observed likelihood-ratio statistic for L in the logistic regression T ~ CT + G + L with statistics computed from null data in which G is replaced by G* = fitted(G ~ L + CG) + permuted residuals. G* keeps its association with L but has no link to T except through L, so these null data represent no mediation. The test is therefore a test for mediation, not for the absence of a direct effect of L on T: it can reject when part of the effect of L on T is mediated by G and part is direct. With test4 = "empirical", the p-value is (r + 1) / (B + 1), where r of B null statistics are at or below the observed statistic; resampling begins with 1000 null statistics and continues, up to maxit, only while the test 4 p-value could determine p_cit. With test4 = "parametric", n.resampl null statistics give a method-of-moments estimate of the non-centrality, lambda = mean(X*) - df, for a non-central chi-square distribution; each statistic is mapped to a normal score through the fitted distribution, and the p-value is the normal probability of the observed score given the mean and standard deviation of the null scores.
With robust = TRUE (CIT v2), the test 3 statistic is referred to a non-central F distribution whose non-centrality is the method-of-moments estimate, lambda = X - df, from the likelihood-ratio statistic X for L in the logistic regression T ~ CT + G + L. For a continuous outcome (cit.cp) the corresponding change guarantees that the CIT cannot conclude that G mediates the effect of L on T and also, with the roles of G and T exchanged (and CT and CG exchanged), that T mediates the effect of L on G; for a binary outcome the same construction is used, but the guarantee is approximate.
With sign.pattern = TRUE (the default; single instrument and single mediator), mediation is also required to be consistent in sign: the estimates in tests 1 to 3 must satisfy sign(T ~ L) = sign(T ~ G given L) x sign(G ~ L given T), as they do when L affects G and G affects T. The alternative hypothesis is the union of the four sign patterns that satisfy this rule. For each pattern the three tests are one-sided in the pattern's direction, and the reported components are those of the pattern whose largest component is smallest. When the estimates agree in sign, tests 1 to 3 therefore report half their two-sided p-values; when they do not, at least one component is 0.5 or more and mediation is not inferred. Association with the wrong sign, as produced for example by a hidden confounder of G and T that reverses the association of G with L given T, no longer counts as evidence for mediation. With robust = TRUE, the components keep their two-sided p-values and a component whose sign disagrees is set to 1. When L or G has several columns the signs are not defined and sign.pattern has no effect. Use sign.pattern = FALSE for the two-sided components of cit 2.3.2.
Each component test uses the observations complete for its own variables: test 1 uses L, T and CT; test 2 L, G, T and CT; test 3 L, G, T and CG; test 4 all of them. With robust = TRUE, test 3 uses the test 4 observations, so that the test 3 statistic and its non-centrality come from the same data.
If permutations are conducted, the results are returned as a data frame in which each row represents an analysis using a unique permutation, except the first row (perm = 0), which has results from the observed data. Each permutation is a single reordering of the n observations, applied to L in tests 1 and 3 (and, with robust = TRUE, in the statistic that gives the test 3 non-centrality), to G in test 2, and to the mediator-model residuals in test 4, so that dependence among the component tests, and among CITs that share a perm.index, is preserved. In tests 1 to 3 the variable is permuted before observations with missing values are removed, so missing values move with their rows; in test 4 the permutation is restricted to the test 4 observations, keeping each pairing whose source observation is complete. These results can be aggregated across tests and passed to fdr.cit.
A data frame. If n.perm = 0 it has one row; otherwise it has a first column perm (0 for the observed data, then 1, ..., n.perm). Columns:
p_cit |
CIT (omnibus) p-value. |
p_TassocL |
component p-value for the test of association between T and L. |
p_TassocGgvnL |
component p-value for the test of association between T and G given L. |
p_GassocLgvnT |
component p-value for the test of association between G and L given T. |
With sign.pattern = TRUE, the components for tests 1 to 3 are the one-sided p-values of the selected sign pattern (see Details), and p_cit is the maximum of the reported components.
p_LindTgvnG |
component p-value for test 4, that the association between L and T given G is weaker than expected without mediation. |
Joshua Millstein, Mingzhi Ye
Millstein J, Chen GK, Breton CV. 2016. cit: hypothesis testing software for mediation analysis in genomic applications. Bioinformatics. btw135. PMID: 27153715.
Millstein J, Zhang B, Zhu J, Schadt EE. 2009. Disentangling molecular relationships with a causal inference test. BMC Genetics, 10:23.
set.seed(1) ss <- 200 L <- matrix(rbinom(ss * 3, 2, .5), ncol = 3) G <- matrix(apply(.3 * L, 1, sum) + rnorm(ss), ncol = 1) T <- rbinom(ss, 1, plogis(.8 * G)) CT <- matrix(rnorm(ss * 2), ncol = 2) cit.bp(L, G, T) cit.bp(L, G, T, CT, robust = TRUE) cit.bp(L, G, T, CT, n.perm = 5)set.seed(1) ss <- 200 L <- matrix(rbinom(ss * 3, 2, .5), ncol = 3) G <- matrix(apply(.3 * L, 1, sum) + rnorm(ss), ncol = 1) T <- rbinom(ss, 1, plogis(.8 * G)) CT <- matrix(rnorm(ss * 2), ncol = 2) cit.bp(L, G, T) cit.bp(L, G, T, CT, robust = TRUE) cit.bp(L, G, T, CT, n.perm = 5)
This function implements a formal statistical hypothesis test, resulting in a p-value, to quantify uncertainty in a causal inference pertaining to a measured factor, e.g. a molecular species, which potentially mediates a known causal association between a locus or other instrumental variable and a censored survival time. If the number of permutations is greater than zero, then the results can be used with fdr.cit to generate permutation-based FDR values (q-values) that are returned with confidence intervals to quantify uncertainty in the estimate. The outcome is a possibly censored survival time, the potential mediator is continuous and may be multivariate, and the instrumental variable can be continuous, discrete (such as coding a SNP 0, 1, 2), or binary and is not limited to a single variable but may be a design matrix representing multiple variables.
cit.cox(L, G, T, E, CT = NULL, n.resampl = 50, n.perm = 0, perm.index = NULL, rseed = NULL, CG = NULL, test4 = c("empirical", "parametric"), maxit = 10000, test3 = TRUE)cit.cox(L, G, T, E, CT = NULL, n.resampl = 50, n.perm = 0, perm.index = NULL, rseed = NULL, CG = NULL, test4 = c("empirical", "parametric"), maxit = 10000, test3 = TRUE)
L |
Vector or n x p design matrix representing the instrumental variable(s). |
G |
Continuous vector, or n x k matrix of k continuous potential mediators tested jointly. |
T |
Numeric vector of follow-up times. |
E |
Event indicator, coded 1 for an event and 0 for censoring. |
CT |
Vector or n x p matrix of adjustment covariates for the models with |
n.resampl |
Number of null instances of the statistic for the last test used to estimate the null distribution when |
n.perm |
If greater than 0, |
perm.index |
Optional n x n.perm matrix whose columns are permutations of |
rseed |
Optional integer seed. When supplied, results are reproducible and the caller's random number stream is left unchanged. When |
CG |
Vector or n x p matrix of adjustment covariates for the models with |
test4 |
How the p-value for the last test is obtained from null data: |
maxit |
Maximum number of null statistics for the empirical version of the last test. |
test3 |
If |
Let L denote the instrumental variable, G the potential mediator and T the outcome. For a survival outcome the CIT is an intersection-union test of four component conditions, each with its own p-value:
T is associated with L (p_TassocL; likelihood-ratio test for L in the Cox proportional hazards model T ~ CT + L, Efron ties).
T is associated with G given L (p_TassocGgvnL; likelihood-ratio test for G in the Cox model T ~ CT + L + G).
G is associated with L given the event indicator E (p_GassocLgvnT; F test for L in the linear model G ~ E + CG + L).
The association between L and T given G is weaker than expected if G were not a mediator (p_LindTgvnG).
The third test guards against reverse causation, here the event changing G, as when tissue is collected at or after the event. If G depends on L only through the event, L and G are independent given E and the test does not reject. If G is measured before any event can occur, reverse causation is ruled out by design, and the test can be omitted with test3 = FALSE; p_cit is then the maximum of the other three components.
The omnibus p-value, p_cit, is the maximum of the component p-values, representing the probability of the data if at least one of the component null hypotheses is true.
The last test compares the observed likelihood-ratio statistic for L in the Cox model T ~ CT + G + L with statistics computed from null data in which G is replaced by G* = fitted(G ~ L + CG) + permuted residuals. G* keeps its association with L but has no link to T except through L, so these null data represent no mediation. The test is therefore a test for mediation, not for the absence of a direct effect of L on T: it can reject when part of the effect of L on T is mediated by G and part is direct. With test4 = "empirical", the p-value is (r + 1) / (B + 1), where r of B null statistics are at or below the observed statistic; resampling begins with 1000 null statistics and continues, up to maxit, only while this p-value could determine p_cit. With test4 = "parametric", n.resampl null statistics give a method-of-moments estimate of the non-centrality, lambda = mean(X*) - df, for a non-central chi-square distribution; each statistic is mapped to a normal score through the fitted distribution, and the p-value is the normal probability of the observed score given the mean and standard deviation of the null scores.
Each component test uses the observations complete for its own variables: the first test uses L, T, E and CT; the second L, G, T, E and CT; the third L, G, T, E and CG; the last all of L, G, T, E, CT and CG.
If permutations are conducted, the results are returned as a data frame in which each row represents an analysis using a unique permutation, except the first row (perm = 0), which has results from the observed data. Each permutation is a single reordering of the n observations, applied to L in the first and third tests, to G in the second, and to the mediator-model residuals in the last, so that dependence among the component tests, and among CITs that share a perm.index, is preserved. In the first three tests the variable is permuted before observations with missing values are removed, so missing values move with their rows; in the last test the permutation is restricted to its observations, keeping each pairing whose source observation is complete. These results can be aggregated across tests and passed to fdr.cit.
A data frame. If n.perm = 0 it has one row; otherwise it has a first column perm (0 for the observed data, then 1, ..., n.perm). Columns:
p_cit |
CIT (omnibus) p-value. |
p_TassocL |
component p-value for the test of association between T and L. |
p_TassocGgvnL |
component p-value for the test of association between T and G given L. |
p_GassocLgvnT |
component p-value for the test of association between G and L given the event indicator (absent when |
p_LindTgvnG |
component p-value for the last test, that the association between L and T given G is weaker than expected without mediation. |
Joshua Millstein, Mingzhi Ye
Millstein J, Chen GK, Breton CV. 2016. cit: hypothesis testing software for mediation analysis in genomic applications. Bioinformatics. btw135. PMID: 27153715.
Millstein J, Zhang B, Zhu J, Schadt EE. 2009. Disentangling molecular relationships with a causal inference test. BMC Genetics, 10:23.
set.seed(1) ss <- 200 L <- matrix(rbinom(ss * 2, 2, .5), ncol = 2) G <- drop(L %*% c(.3, .3)) + rnorm(ss) T <- rexp(ss, exp(.5 * G)) E <- rbinom(ss, 1, .8) cit.cox(L, G, T, E, maxit = 2000) cit.cox(L, G, T, E, test4 = "parametric") cit.cox(L, G, T, E, test4 = "parametric", test3 = FALSE) # G measured before any eventset.seed(1) ss <- 200 L <- matrix(rbinom(ss * 2, 2, .5), ncol = 2) G <- drop(L %*% c(.3, .3)) + rnorm(ss) T <- rexp(ss, exp(.5 * G)) E <- rbinom(ss, 1, .8) cit.cox(L, G, T, E, maxit = 2000) cit.cox(L, G, T, E, test4 = "parametric") cit.cox(L, G, T, E, test4 = "parametric", test3 = FALSE) # G measured before any event
This function implements a formal statistical hypothesis test, resulting in a p-value, to quantify uncertainty in a causal inference pertaining to a measured factor, e.g. a molecular species, which potentially mediates a known causal association between a locus or other instrumental variable and a quantitative trait. If the number of permutations is greater than zero, then the results can be used with fdr.cit to generate permutation-based FDR values (q-values) that are returned with confidence intervals to quantify uncertainty in the estimate. The outcome is continuous, the potential mediator is continuous and may be multivariate, and the instrumental variable can be continuous, discrete (such as coding a SNP 0, 1, 2), or binary and is not limited to a single variable but may be a design matrix representing multiple variables.
cit.cp(L, G, T, CT = NULL, n.resampl = 50, n.perm = 0, perm.index = NULL, rseed = NULL, CG = NULL, robust = FALSE, test4 = c("empirical", "parametric"), maxit = 10000, sign.pattern = TRUE, C = NULL)cit.cp(L, G, T, CT = NULL, n.resampl = 50, n.perm = 0, perm.index = NULL, rseed = NULL, CG = NULL, robust = FALSE, test4 = c("empirical", "parametric"), maxit = 10000, sign.pattern = TRUE, C = NULL)
L |
Vector or n x p design matrix representing the instrumental variable(s). |
G |
Continuous vector, or n x k matrix of k continuous potential mediators tested jointly. |
T |
Continuous vector representing the clinical trait or outcome of interest. |
CT |
Vector or n x p matrix of adjustment covariates for the models with |
n.resampl |
Number of null instances of the test 4 statistic used to estimate the null distribution when |
n.perm |
If greater than 0, |
perm.index |
Optional n x n.perm matrix whose columns are permutations of |
rseed |
Optional integer seed. When supplied, results are reproducible and the caller's random number stream is left unchanged. When |
CG |
Vector or n x p matrix of adjustment covariates for the models with |
robust |
If |
test4 |
How the test 4 p-value is obtained from null data: |
maxit |
Maximum number of null statistics for the empirical test 4. |
sign.pattern |
If |
C |
Deprecated name for |
Let L denote the instrumental variable, G the potential mediator and T the outcome. The CIT is an intersection-union test of four component conditions, each with its own p-value:
T is associated with L (p_TassocL; F test of L in T ~ CT + L).
T is associated with G given L (p_TassocGgvnL; F test of G in T ~ CT + L + G).
G is associated with L given T (p_GassocLgvnT; F test of L in G ~ T + CG + L, or Wilks' lambda when G has several columns).
The association between L and T given G is weaker than expected if G were not a mediator (p_LindTgvnG).
The omnibus p-value, p_cit, is the maximum of the component p-values, representing the probability of the data if at least one of the component null hypotheses is true.
Test 4 compares the observed F statistic for L in T ~ CT + G + L with statistics computed from null data in which G is replaced by G* = fitted(G ~ L + CG) + permuted residuals. G* keeps its association with L but has no link to T except through L, so these null data represent no mediation. The test is therefore a test for mediation, not for the absence of a direct effect of L on T: it can reject when part of the effect of L on T is mediated by G and part is direct. With test4 = "empirical", the p-value is (r + 1) / (B + 1), where r of B null statistics are at or below the observed statistic; resampling begins with 1000 null statistics and continues, up to maxit, only while the test 4 p-value could determine p_cit. With test4 = "parametric", n.resampl null statistics give a method-of-moments estimate of the non-centrality, lambda = mean(F*) df1 (df2 - 2) / df2 - df1; each statistic is mapped to a normal score through the fitted non-central F distribution, and the p-value is the normal probability of the observed score given the mean and standard deviation of the null scores.
With robust = TRUE (CIT v2), the test 3 statistic is referred to a non-central F distribution whose non-centrality is the method-of-moments estimate from the F statistic for L in T ~ CT + G + L. With this change the CIT cannot conclude that G mediates the effect of L on T and also, with the roles of G and T exchanged (and CT and CG exchanged), that T mediates the effect of L on G.
With sign.pattern = TRUE (the default; single instrument and single mediator), mediation is also required to be consistent in sign: the estimates in tests 1 to 3 must satisfy sign(T ~ L) = sign(T ~ G given L) x sign(G ~ L given T), as they do when L affects G and G affects T. The alternative hypothesis is the union of the four sign patterns that satisfy this rule. For each pattern the three tests are one-sided in the pattern's direction, and the reported components are those of the pattern whose largest component is smallest. When the estimates agree in sign, tests 1 to 3 therefore report half their two-sided p-values; when they do not, at least one component is 0.5 or more and mediation is not inferred. Association with the wrong sign, as produced for example by a hidden confounder of G and T that reverses the association of G with L given T, no longer counts as evidence for mediation. With robust = TRUE, the components keep their two-sided p-values and a component whose sign disagrees is set to 1. When L or G has several columns the signs are not defined and sign.pattern has no effect. Use sign.pattern = FALSE for the two-sided components of cit 2.3.2.
Each component test uses the observations complete for its own variables: test 1 uses L, T and CT; test 2 L, G, T and CT; test 3 L, G, T and CG; test 4 all of them. With robust = TRUE, test 3 uses the test 4 observations, so that the test 3 statistic and its non-centrality come from the same data.
If permutations are conducted, the results are returned as a data frame in which each row represents an analysis using a unique permutation, except the first row (perm = 0), which has results from the observed data. Each permutation is a single reordering of the n observations, applied to L in tests 1 and 3 (and, with robust = TRUE, in the statistic that gives the test 3 non-centrality), to G in test 2, and to the mediator-model residuals in test 4, so that dependence among the component tests, and among CITs that share a perm.index, is preserved. In tests 1 to 3 the variable is permuted before observations with missing values are removed, so missing values move with their rows; in test 4 the permutation is restricted to the test 4 observations, keeping each pairing whose source observation is complete. These results can be aggregated across tests and passed to fdr.cit.
A data frame. If n.perm = 0 it has one row; otherwise it has a first column perm (0 for the observed data, then 1, ..., n.perm). Columns:
p_cit |
CIT (omnibus) p-value. |
p_TassocL |
component p-value for the test of association between T and L. |
p_TassocGgvnL |
component p-value for the test of association between T and G given L. |
p_GassocLgvnT |
component p-value for the test of association between G and L given T. |
With sign.pattern = TRUE, the components for tests 1 to 3 are the one-sided p-values of the selected sign pattern (see Details), and p_cit is the maximum of the reported components.
p_LindTgvnG |
component p-value for test 4, that the association between L and T given G is weaker than expected without mediation. |
Joshua Millstein, Mingzhi Ye
Millstein J, Chen GK, Breton CV. 2016. cit: hypothesis testing software for mediation analysis in genomic applications. Bioinformatics. btw135. PMID: 27153715.
Millstein J, Zhang B, Zhu J, Schadt EE. 2009. Disentangling molecular relationships with a causal inference test. BMC Genetics, 10:23.
set.seed(1) ss <- 100 L <- matrix(rbinom(ss * 3, 2, .5), ncol = 3) G <- matrix(apply(.3 * L, 1, sum) + rnorm(ss), ncol = 1) T <- matrix(.3 * G + rnorm(ss), ncol = 1) CT <- matrix(rnorm(ss * 2), ncol = 2) cit.cp(L, G, T) cit.cp(L, G, T, CT, robust = TRUE) cit.cp(L, G, T, CT, test4 = "parametric") perm.index <- replicate(5, sample(ss)) cit.cp(L, G, T, CT, perm.index = perm.index)set.seed(1) ss <- 100 L <- matrix(rbinom(ss * 3, 2, .5), ncol = 3) G <- matrix(apply(.3 * L, 1, sum) + rnorm(ss), ncol = 1) T <- matrix(.3 * G + rnorm(ss), ncol = 1) CT <- matrix(rnorm(ss * 2), ncol = 2) cit.cp(L, G, T) cit.cp(L, G, T, CT, robust = TRUE) cit.cp(L, G, T, CT, test4 = "parametric") perm.index <- replicate(5, sample(ss)) cit.cp(L, G, T, CT, perm.index = perm.index)
Estimates the total, direct and indirect effects of the instrumental variable L on the outcome T through the mediator G, with delta-method confidence intervals, using the regression estimators of VanderWeele and Vansteelandt in the form given by Valeri and VanderWeele (2013). It is intended to accompany cit.cp, cit.bp and cit.cox: the CIT tests whether the data support mediation, and these estimates describe its size.
cit.effects(L, G, T, E = NULL, CT = NULL, CG = NULL, outcome = c("continuous", "binary", "survival"), conf.level = 0.95)cit.effects(L, G, T, E = NULL, CT = NULL, CG = NULL, outcome = c("continuous", "binary", "survival"), conf.level = 0.95)
L |
Vector or n x p design matrix representing the instrumental variable(s). Effects are reported for each column. |
G |
Continuous vector: a single potential mediator. |
T |
Outcome: a continuous vector, a binary vector coded 0/1, or a vector of follow-up times. |
E |
Event indicator (1 = event, 0 = censored), for |
CT |
Vector or n x p matrix of adjustment covariates for the outcome model. |
CG |
Vector or n x p matrix of adjustment covariates for the mediator model. |
outcome |
Type of outcome: |
conf.level |
Confidence level of the intervals. |
Two models are fitted to the observations complete for every variable used: the mediator model G ~ L + CG (linear) and the outcome model T ~ L + G + CT (linear, logistic or Cox). As in the CIT, the models have no L x G interaction. Let b1 be the coefficient of a column of L in the mediator model, and t1 and t2 the coefficients of that column and of G in the outcome model. For a one-unit increase in that column, with the other columns held fixed:
direct effect = t1;
indirect effect = t2 b1;
total effect = t1 + t2 b1.
For a continuous outcome these are mean differences. For a binary outcome they are log odds ratios and for a survival outcome log hazard ratios, reported as odds ratios or hazard ratios. The ratio-scale decomposition assumes a rare outcome (VanderWeele and Vansteelandt 2010; VanderWeele 2011), and a warning is given when more than 10% of observations have the event. The proportion mediated is indirect / total for a continuous outcome and OR_D (OR_I - 1) / (OR_D OR_I - 1) on the ratio scales (with HR for OR).
Standard errors are by the delta method, treating the two models as independent (Valeri and VanderWeele 2013). The p-values are Wald tests of a zero effect (on the log scale for ratios). The p-value for the indirect effect is a Sobel-type test; it is not the CIT, which also tests the direction of the relationship.
These estimates rely on stronger assumptions than the CIT. Interpreting them as natural direct and indirect effects requires no unmeasured confounding of the L-T, L-G and G-T relationships given the covariates, and no confounder of G and T that is itself affected by L. When L is a genetic variant, its relationships with G and T are protected by random assortment, but the G-T relationship is not. Measurement error in G biases the indirect effect towards zero and the direct effect away from zero.
A data frame with one row per column of L and effect: L (column name), effect ("total", "direct", "indirect", "prop.mediated"), scale ("difference", "OR", "HR" or "proportion"), estimate, lower, upper, se (on the log scale for ratios), p and n (observations used).
VanderWeele TJ, Vansteelandt S. 2009. Conceptual issues concerning mediation, interventions and composition. Statistics and Its Interface 2:457-468.
VanderWeele TJ, Vansteelandt S. 2010. Odds ratios for mediation analysis for a dichotomous outcome. American Journal of Epidemiology 172:1339-1348.
VanderWeele TJ. 2011. Causal mediation analysis with survival data. Epidemiology 22:582-585.
Valeri L, VanderWeele TJ. 2013. Mediation analysis allowing for exposure-mediator interactions and causal interpretation: theoretical assumptions and implementation with SAS and SPSS macros. Psychological Methods 18:137-150.
set.seed(1) n <- 500 L <- rbinom(n, 2, 0.3) G <- 0.5 * L + rnorm(n) T <- 0.3 * G + 0.1 * L + rnorm(n) cit.effects(L, G, T)set.seed(1) n <- 500 L <- rbinom(n, 2, 0.3) G <- 0.5 * L + rnorm(n) T <- 0.3 * G + 0.1 * L + rnorm(n) cit.effects(L, G, T)
Results including permutation results from cit.cp, cit.bp or cit.cox from multiple tests are used in a permutation-based approach to compute false discovery rates (q-values) for each test. q-values are computed for the omnibus test as well as all component tests. q-values are returned with confidence intervals to quantify uncertainty in the estimate.
fdr.cit( cit.perm.list, cl=.95, c1=NA )fdr.cit( cit.perm.list, cl=.95, c1=NA )
cit.perm.list |
List, where each item is a dataframe output from cit.cp, cit.bp or cit.cox, which must be run with n.perm > 0 and the same perm.index. |
cl |
Confidence level for the q-value confidence interval (default is .95). |
c1 |
Over-dispersion parameter. Setting the over-dispersion parameter to one corresponds to the assumption that all tests (referring to multiple omnibus tests conducted using different sets of variables) are independent. The default (c1=NA) is to use an empirically estimated over-dispersion parameter based on the permutation results. |
For the initial cit.cp and cit.bp test computations, we suggest setting n.perm to 100 or greater for moderate numbers of tests, say 10s or 100s. Larger numbers of permutations will generate more precise FDR estimates but require substantial computation time. If a large number of tests, say 1000s, have been conducted, then n.perm can be as small as 10 and still generate good FDR estimates and confidence intervals. FDR confidence intervals are especially useful to quantify uncertainty when FDR value > 0.05. The over-dispersion parameter should be set to one if it is known that all tests are independent, however, this is rarely the case in typical multiple testing settings. The fdr.cit function can also be used with simulated data in order to generate power calculations. In this case the over-dispersion parameter may be estimated from pilot data or other datasets and fixed to the prior estimate using the c1 argument. The omnibus CIT FDR value combines the component test FDR values (four component tests, or three for cit.cox with test3 = FALSE) in an intersection-union type of approach analogous to the method used to generate the omnibus CIT p-value; see below.
The omnibus q-value. The component p-values can have different null distributions, for example when a parametric null is only approximate, so they are first put on a common scale. Each component p-value is mapped to its null distribution function , estimated from the permuted p-values of that component pooled over all tests, giving the calibrated value ; the calibrated omnibus value is p.cit.cal . For a test with p.cit.cal , each component FDR is estimated by permutation (Millstein and Volfson 2013) for the tests with , , and the estimates are combined as ; q.cit is the minimum of over the values of p.cit.cal . A common threshold lies at the same null quantile for every component. The rejection regions {p.cit.cal } are nested, so selecting the tests with q.cit selects the tests with p.cit.cal below a single cutoff, as with ordinary q-values (Storey 2002, 2003). Tests whose calibrated values tie, for instance because their component p-values are all smaller than every permuted p-value, share a q-value; more permutations separate them. When the component p-values are uniform under the null, p.cit.cal is close to p.cit. The component q-values (q.TaL, q.TaGgvL, q.GaLgvT, q.LiTgvG) are each evaluated at the test's own component p-value.
The combination is justified under the two-groups model (Efron et al. 2001; Efron 2008), in which each hypothesis is null or non-null at random and its p-value is drawn given that status; under this model the positive FDR (pFDR) of a rejection region equals the probability that a hypothesis in the region is null (Storey 2003). Let be the null statuses of the K component hypotheses of a test. A CIT call is a false discovery when at least one component null holds (the composite null). Suppose (1) over the tests analysed together, the K component null statuses are mutually independent, and (2) given , the component p-values are independent, with each depending only on its own status . Then equals the pFDR of the call set {p.cit.cal }, and q.cit estimates the corresponding q-value consistently when the component estimates are consistent. These conditions concern how component null statuses co-occur among the tests analysed together, which depends on how those tests were selected, and they can fail. If many tests are null for one component and strongly non-null for the others, the omnibus q-value can understate the FDR; if tests under the composite null tend to be null for every component, it can overstate it.
It is important to be aware that there are certain conditions under which this permutation-based FDR is not stable. This includes the case where all observed tests achieve the significance level and the case where there are zero positive tests among the permuted results. In the latter case, we take the conservative approach of setting the number of positive permutation tests to one. Both these cases produce an NA for the upper confidence limit. The latter problem may be alleviated by increasing the number of permutations conducted.
A dataframe which includes the following columns:
p.cit |
CIT (omnibus intersection/union test, IUT) p-value, max of component p-values |
q.cit |
CIT (omnibus FDR intersection/union test, IUT) q-value: permutation-based q-value for p.cit.cal, combining the component FDRs evaluated at a common null quantile (see Details) |
q.cit.ll |
Lower 95 percent confidence limit (if cl=.95) for q.cit estimate |
q.cit.ul |
Upper 95 percent confidence limit (if cl=.95) for q.cit estimate |
q.TaL |
component q-value (FDR) for the test of association between T and L |
q.ll.TaL |
Lower 95 percent confidence limit (if cl=.95) for q.TaL estimate |
q.ul.TaL |
Upper 95 percent confidence limit (if cl=.95) for q.TaL estimate |
q.TaGgvL |
component q-value (FDR) for the test of association between T and G|L |
q.ll.TaGgvL |
Lower 95 percent confidence limit (if cl=.95) for q.TaGgvL estimate |
q.ul.TaGgvL |
Upper 95 percent confidence limit (if cl=.95) for q.TaGgvL estimate |
q.GaLgvT |
component q-value (FDR) for the test of association between G and L|T (for cit.cox results, G and L given the event indicator; absent with test3 = FALSE) |
q.ll.GaLgvT |
Lower 95 percent confidence limit (if cl=.95) for q.GaLgvT estimate |
q.ul.GaLgvT |
Upper 95 percent confidence limit (if cl=.95) for q.GaLgvT estimate |
q.LiTgvG |
component q-value (FDR) for test 4, that the association between L and T given G is weaker than expected without mediation |
q.ll.LiTgvG |
Lower 95 percent confidence limit (if cl=.95) for q.LiTgvG estimate |
q.ul.LiTgvG |
Upper 95 percent confidence limit (if cl=.95) for q.LiTgvG estimate |
The component p-values are also returned, in columns named as in the input, followed by
p.cit.cal |
calibrated omnibus value, the largest of the component p-values each mapped to its permutation null distribution function; tests are ranked by p.cit.cal for q.cit (see Details) |
Joshua Millstein
Millstein J, Chen GK, Breton CV. 2016. cit: hypothesis testing software for mediation analysis in genomic applications. Bioinformatics. btw135. PMID: 27153715. Millstein J, Volfson D. 2013. Computationally efficient permutation-based confidence interval estimation for tail-area FDR. Frontiers in Genetics | Statistical Genetics and Methodology 4(179):1-11.
Efron B, Tibshirani R, Storey JD, Tusher V. 2001. Empirical Bayes analysis of a microarray experiment. Journal of the American Statistical Association 96(456):1151-1160.
Efron B. 2008. Microarrays, empirical Bayes and the two-groups model. Statistical Science 23(1):1-22.
Storey JD. 2002. A direct approach to false discovery rates. Journal of the Royal Statistical Society, Series B 64(3):479-498.
Storey JD. 2003. The positive false discovery rate: a Bayesian interpretation and the q-value. Annals of Statistics 31(6):2013-2035.
# Sample Size ss = 100 n.perm = 20 perm.index = matrix(NA, nrow=ss, ncol=n.perm ) for( j in 1:ncol(perm.index) ) perm.index[, j] = sample( 1:ss ) n.tests = 20 myresults = vector('list', n.tests) for( tst in 1:n.tests ){ # Errors e1 = matrix(rnorm(ss),ncol=1) e2 = matrix(rnorm(ss),ncol=1) # Simulate genotypes, gene expression, and clinical traits L = matrix(rbinom(ss,2,.5),ncol=1) G = matrix(.5*L + e1,ncol=1) T = matrix(.3*G + e2,ncol=1) T = ifelse( T > median(T), 1, 0 ) myresults[[ tst ]] = cit.bp(L, G, T, perm.index=perm.index, n.perm=n.perm) } fdr.cit( myresults )# Sample Size ss = 100 n.perm = 20 perm.index = matrix(NA, nrow=ss, ncol=n.perm ) for( j in 1:ncol(perm.index) ) perm.index[, j] = sample( 1:ss ) n.tests = 20 myresults = vector('list', n.tests) for( tst in 1:n.tests ){ # Errors e1 = matrix(rnorm(ss),ncol=1) e2 = matrix(rnorm(ss),ncol=1) # Simulate genotypes, gene expression, and clinical traits L = matrix(rbinom(ss,2,.5),ncol=1) G = matrix(.5*L + e1,ncol=1) T = matrix(.3*G + e2,ncol=1) T = ifelse( T > median(T), 1, 0 ) myresults[[ tst ]] = cit.bp(L, G, T, perm.index=perm.index, n.perm=n.perm) } fdr.cit( myresults )
This function can be used to estimate FDR, corresponding confidence interval, and pi0, the proportion of true null hypotheses, given a selected significance threshold, and results from permuted data.
fdr.od(obsp, permp, pnm, ntests, thres, cl=.95, od=NA)fdr.od(obsp, permp, pnm, ntests, thres, cl=.95, od=NA)
obsp |
observed vector of p-values. |
permp |
list of dataframes that include a column of permutation p-values (or statistics) in each. The length of the list permp = number of permutations. |
pnm |
name of column in each list component dataframe that includes p-values (or statistics). |
ntests |
total number of observed tests, which is usually the same as the length of obsp and the number of rows in each permp dataframe. However, this may not be the case if results were filtered by a p-value threshold or statistic threshold. If filtering was conducted then thres must be smaller (more extreme) than the filtering criterion. |
thres |
significance threshold. |
cl |
confidence level for FDR confidence interval (default is .95). |
od |
Over-dispersion parameter. Setting the over-dispersion parameter to one corresponds to the assumption that all tests (referring to multiple omnibus tests conducted using different sets of variables) are independent. The default (c1=NA) is to use an empirically estimated over-dispersion parameter based on the permutation results. |
If a very large number of tests are conducted, it may be useful to filter results, that is, save only results of those tests that meet some relaxed nominal significance threshold. This alleviates the need to record results for tests that are clearly non-significant. Results from fdr_od() are valid as long as thres < the relaxed nomimal significance threshold for both observed and permuted results. It is not necessary for the input to fdr_od() to be p-values, however, fdr_od() is designed for statistics in which smaller values are more extreme than larger values as is the case for p-values. Therefore, if raw statistics are used, then a transformation may be necessary to insure that smaller values are more likely associated with false null hypotheses than larger values.
If there are zero positive tests among any of the permuted results (s.perm = 0) at the specified significance threshold (thres), then a conservative approximation is conducted by setting s.perm = 1 for the FDR calculation. This approximation is necessary for the permutation-based FDR to be estimable. This situation can be identified by noting that s.perm = 0 in the output. Increasing the number of permutations conducted may result in positive permutation tests, and thereby overcome the problem. In certain situations, for instance when a large proportion of tests meet the significance threshold, pi0 is estimated to be very small, and thus has a large influence on the FDR estimate. To limit this influence, pi0 is constrained to be .5 or greater, resulting in a more conservative estimate under these conditions.
A vector which includes:
FDR |
FDR point estimate |
ll |
lower confidence limit |
ul |
upper confidence limit |
pi0 |
proportion of true null hypotheses |
od |
overdispersion parameter |
s.obs |
observed number of positive tests |
s.perm |
total number of positive tests summed across all permuted result sets |
Joshua Millstein
Millstein J, Chen GK, Breton CV. 2016. cit: hypothesis testing software for mediation analysis in genomic applications. Bioinformatics. btw135. PMID: 27153715. Millstein J, Volfson D. 2013. Computationally efficient permutation-based confidence interval estimation for tail-area FDR. Frontiers in Genetics | Statistical Genetics and Methodology 4(179):1-11.
nrow_=100 ncol_=100 X = as.data.frame(matrix(rnorm(nrow_*ncol_),nrow=nrow_,ncol=ncol_)) Y = as.data.frame(matrix(rnorm(nrow_*ncol_),nrow=nrow_,ncol=ncol_)) nperm = 10 myanalysis = function(X,Y){ ntests = ncol(X) rslts = as.data.frame(matrix(NA,nrow=ntests,ncol=2)) names(rslts) = c("ID","pvalue") rslts[,"ID"] = 1:ntests for(i in 1:ntests){ fit = cor.test(X[,i],Y[,i],na.action="na.exclude", alternative="two.sided",method="pearson") rslts[i,"pvalue"] = fit$p.value } return(rslts) } # End myanalysis # Generate observed results obs = myanalysis(X,Y) ## Generate permuted results perml = vector('list',nperm) for(p_ in 1:nperm){ X1 = X[order(runif(ncol_)),] perml[[p_]] = myanalysis(X1,Y) } ## FDR results fdr.od(obs$pvalue,perml,"pvalue",ncol_,.05)nrow_=100 ncol_=100 X = as.data.frame(matrix(rnorm(nrow_*ncol_),nrow=nrow_,ncol=ncol_)) Y = as.data.frame(matrix(rnorm(nrow_*ncol_),nrow=nrow_,ncol=ncol_)) nperm = 10 myanalysis = function(X,Y){ ntests = ncol(X) rslts = as.data.frame(matrix(NA,nrow=ntests,ncol=2)) names(rslts) = c("ID","pvalue") rslts[,"ID"] = 1:ntests for(i in 1:ntests){ fit = cor.test(X[,i],Y[,i],na.action="na.exclude", alternative="two.sided",method="pearson") rslts[i,"pvalue"] = fit$p.value } return(rslts) } # End myanalysis # Generate observed results obs = myanalysis(X,Y) ## Generate permuted results perml = vector('list',nperm) for(p_ in 1:nperm){ X1 = X[order(runif(ncol_)),] perml[[p_]] = myanalysis(X1,Y) } ## FDR results fdr.od(obs$pvalue,perml,"pvalue",ncol_,.05)
Given a vector of p-values, this function will use a parametric expression for tail-area FDR to assign q-values, minimum tail-area FDR, for each p-value. This function is used by the function fdr.cit to generate q-values for the forth test, LindTgvnG.
fdr.q.para( pvals )fdr.q.para( pvals )
pvals |
Vector of p-values. |
The parametric expression is constructed by setting the average number of positive tests among permuted results to p*m, where p is the p-value and m is the number of tests, the length of vector pvals. Following Storey and Tibshirani (2003), the q-value is set to the minimum FDR for p-values greater than or equal to the observed p-value.
A vector of q-values.
Joshua Millstein
Millstein J, Chen GK, Breton CV. 2016. cit: hypothesis testing software for mediation analysis in genomic applications. Bioinformatics. btw135. Millstein J, Volfson D. 2013. Computationally efficient permutation-based confidence interval estimation for tail-area FDR. Frontiers in Genetics | Statistical Genetics and Methodology 4(179):1-11. Storey DJ, Tibshirani R. 2003. Statistical significance for genomewide studies. PNAS 100(16):9440-9445.
# Sample Size ss = 100 n.tests = 20 mypvals = rep(NA, n.tests) for( tst in 1:n.tests ){ # Errors e1 = matrix(rnorm(ss),ncol=1) e2 = matrix(rnorm(ss),ncol=1) # Simulate genotypes, gene expression, and clinical traits L = matrix(rbinom(ss,2,.5),ncol=1) G = matrix(.5*L + e1,ncol=1) T = matrix(.3*G + e2,ncol=1) T = ifelse( T > median(T), 1, 0 ) mypvals[ tst ] = cit.bp(L, G, T)$p_cit } fdr.q.para( mypvals )# Sample Size ss = 100 n.tests = 20 mypvals = rep(NA, n.tests) for( tst in 1:n.tests ){ # Errors e1 = matrix(rnorm(ss),ncol=1) e2 = matrix(rnorm(ss),ncol=1) # Simulate genotypes, gene expression, and clinical traits L = matrix(rbinom(ss,2,.5),ncol=1) G = matrix(.5*L + e1,ncol=1) T = matrix(.3*G + e2,ncol=1) T = ifelse( T > median(T), 1, 0 ) mypvals[ tst ] = cit.bp(L, G, T)$p_cit } fdr.q.para( mypvals )
Given a vector of p-values, this function uses fdr.od() to compute tail-area FDR to assign q-values, minimum tail-area FDR, for each p-value. This function is used by the function fdr.cit to generate q-values for the three component tests, TassocL, TassocGgvnL, and GassocLgvnT.
fdr.q.perm( obs.p, perml, pname, ntests, cl=.95, od=NA )fdr.q.perm( obs.p, perml, pname, ntests, cl=.95, od=NA )
obs.p |
observed vector of p-values. |
perml |
list of dataframes that include a column of permutation p-values (or statistics) in each. The length of the list permp = number of permutations. |
pname |
name of column in each list component dataframe that includes p-values (or statistics). |
ntests |
total number of observed tests, which is usually the same as the length of obsp and the number of rows in each permp dataframe. However, this may not be the case if results were filtered by a p-value threshold or statistic threshold. If filtering was conducted then thres must be smaller (more extreme) than the filtering criterion. |
cl |
confidence level for FDR confidence interval (default is .95). |
od |
Over-dispersion parameter. Setting the over-dispersion parameter to one corresponds to the assumption that all tests (referring to multiple omnibus tests conducted using different sets of variables) are independent. The default (c1=NA) is to use an empirically estimated over-dispersion parameter based on the permutation results. |
Following Storey and Tibshirani (2003), the q-value is set to the minimum FDR for p-values greater than or equal to the observed p-value.
A vector of q-values.
Joshua Millstein
Millstein J, Chen GK, Breton CV. 2016. cit: hypothesis testing software for mediation analysis in genomic applications. Bioinformatics. btw135. Millstein J, Volfson D. 2013. Computationally efficient permutation-based confidence interval estimation for tail-area FDR. Frontiers in Genetics | Statistical Genetics and Methodology 4(179):1-11. Storey DJ, Tibshirani R. 2003. Statistical significance for genomewide studies. PNAS 100(16):9440-9445.
rowno=100 colno=100 X = as.data.frame(matrix(rnorm(rowno*colno),nrow=rowno,ncol=colno)) Y = as.data.frame(matrix(rnorm(rowno*colno),nrow=rowno,ncol=colno)) nperm = 10 myanalysis = function(X,Y){ ntests = ncol(X) rslts = as.data.frame(matrix(NA,nrow=ntests,ncol=2)) names(rslts) = c("ID","pvalue") rslts[,"ID"] = 1:ntests for(i in 1:ntests){ fit = cor.test(X[,i],Y[,i],na.action="na.exclude", alternative="two.sided",method="pearson") rslts[i,"pvalue"] = fit$p.value } return(rslts) } # End myanalysis # Generate observed results obs = myanalysis(X,Y) ## Generate permuted results perml = vector('list',nperm) for(perm in 1:nperm){ X1 = X[order(runif(colno)),] perml[[perm]] = myanalysis(X1,Y) } ## FDR results fdr.q.perm(obs$pvalue,perml,"pvalue",colno)rowno=100 colno=100 X = as.data.frame(matrix(rnorm(rowno*colno),nrow=rowno,ncol=colno)) Y = as.data.frame(matrix(rnorm(rowno*colno),nrow=rowno,ncol=colno)) nperm = 10 myanalysis = function(X,Y){ ntests = ncol(X) rslts = as.data.frame(matrix(NA,nrow=ntests,ncol=2)) names(rslts) = c("ID","pvalue") rslts[,"ID"] = 1:ntests for(i in 1:ntests){ fit = cor.test(X[,i],Y[,i],na.action="na.exclude", alternative="two.sided",method="pearson") rslts[i,"pvalue"] = fit$p.value } return(rslts) } # End myanalysis # Generate observed results obs = myanalysis(X,Y) ## Generate permuted results perml = vector('list',nperm) for(perm in 1:nperm){ X1 = X[order(runif(colno)),] perml[[perm]] = myanalysis(X1,Y) } ## FDR results fdr.q.perm(obs$pvalue,perml,"pvalue",colno)
Tail-area false discovery rates (FDRs) or 'q-values' are combined across multiple component tests to estimate a single q-value that represents the intersection of alternative hypotheses or the union of null hypotheses. In other words, this is an estimate of the rate at which tests called 'significant' at the specified omnibus q-value include at least one component condition that is called 'significant' due to chance alone.
iuq( qvec )iuq( qvec )
qvec |
vector of q-values, each corresponding to a component test. |
The omnibus q-value is greater than or equal to the maximum of the component q-values.
A single value that represents the rate at which tests called 'significant' at the specified omnibus q-value include at least one component condition that is called 'significant' due to chance alone.
Joshua Millstein
Millstein J, Chen GK, Breton CV. 2016. cit: hypothesis testing software for mediation analysis in genomic applications. Bioinformatics. btw135. PMID: 27153715.
# component q-values qvec = c( .1, .04, .02, .11 ) # omnibus q-value iuq( qvec )# component q-values qvec = c( .1, .04, .02, .11 ) # omnibus q-value iuq( qvec )
This function is used by cit.cp to compute F test given a continuous outcome and full vs reduced sets of covariates
linreg( nms.full, nms.redu=NULL, nm.y, mydat )linreg( nms.full, nms.redu=NULL, nm.y, mydat )
nms.full |
vector of variable names for all covariates included in the full model. |
nms.redu |
vector of variable names for all covariates included in the reduced model. If nms.redu is NULL then the reduced model is fitted with the intercept only. |
nm.y |
character string, which is the name of the outcome variable. |
mydat |
the dataframe that includes all variables with each variable in a column. |
An F test is conducted using the glm function by comparing the full and reduced models. This function is called by cit.cp.
A single p-value is returned.
Joshua Millstein
Millstein J, Zhang B, Zhu J, Schadt EE. 2009. Disentangling molecular relationships with a causal inference test. BMC Genetics, 10:23.
ss = 500 cols = 6 nm.y = "y" nms.full = paste( "x", 1:(cols-1), sep="" ) nms.redu = paste( "x", 1:2, sep="" ) mydat = as.data.frame( matrix( rnorm( ss*cols ), ncol=cols ) ) names( mydat ) = c( nm.y, nms.full ) linreg(nms.full, nms.redu, nm.y, mydat)ss = 500 cols = 6 nm.y = "y" nms.full = paste( "x", 1:(cols-1), sep="" ) nms.redu = paste( "x", 1:2, sep="" ) mydat = as.data.frame( matrix( rnorm( ss*cols ), ncol=cols ) ) names( mydat ) = c( nm.y, nms.full ) linreg(nms.full, nms.redu, nm.y, mydat)